Long‐distance dispersal or postglacial contraction? Insights into disjunction between Himalaya–Hengduan Mountains and Taiwan in a cold‐adapted herbaceous genus, Triplostegia

Abstract Current disjunct patterns can result from long‐distance dispersal or postglacial contraction. We herein investigate the evolutionary history of Triplostegia to elucidate the disjunction between the Himalaya–Hengduan Mountain region (HHM) and Taiwan (TW). Genetic structure of Triplostegia was investigated for 48 populations using sequences from five chloroplast loci and the ribosomal nuclear internal transcribed spacer. Divergence time estimation, ancestral area reconstruction, and species distribution modeling (SDM) were employed to examine the biogeographic history of Triplostegia. Substantial genetic differentiation among populations from southwestern China (SW), Central China (CC), and TW was detected. Triplostegia was inferred to have originated in SW, and diversification began during the late Miocene; CC was colonized in the mid‐Pliocene, and TW was finally colonized in the early Pleistocene. SDM suggested an expansion of climatically suitable areas during the Last Glacial Maximum and range contraction during the Last interglacial in Triplostegia. Disjunction between HHM and TW in Triplostegia is most likely the consequence of topographic isolation and postglacial contraction. The potential climatic suitability areas for Triplostegia by 2070s (2061–2080) are predicted to slightly shrink and move northward. With continued global warming and human‐induced deforestation, extinction risk may increase for the cold‐adapted species, and appropriate strategies should be employed for ecosystem conservation.


| INTRODUCTION
Organismal genetic structure and distributions have been greatly impacted by climatic oscillations, especially the Quaternary glacial/interglacial cycles (Comes & Kadereit, 1998;Hewitt, 2000;Lomolino, Riddle, Whittaker, & Brown, 2010). Compared with North America and northern Europe, East Asia was inferred to be less affected by glaciation except in the higher mountains (Pinot et al., 1999;Shi & Yao, 2002). Numerous species show disjunct distributions in the high mountains of East Asia Liu et al., 2014;Meng et al., 2017;Wang, 1989;Wang et al., 2013). In particular, disjunctions between the Himalaya-Hengduan Mountain region (HHM) and Taiwan (TW) have received extensive attention because of the geographically distant separation of the two regions by the Taiwan Strait and a terrestrial distance of over 2,000 km (Chen, Ying, & Lu, 2012;Wu & Wang, 2001;Ye, Chen, Liu, Qin, & Yang, 2012).
Long-distance dispersal in plants usually involves rare events driven by dispersal vectors such as animals, currents, winds, or water (Nathan, 2006). Phylogeographic study of Cunninghamia revealed that populations in Taiwan were derived from continental Asia via longdistance seed dispersal (Lu et al., 2001). The postglacial contraction hypothesis states that many species had a continuous distribution from Mainland China to TW during the Quaternary glacial period, and disjunction between HHM and TW occurred when taxa migrated to higher altitudes as the Earth became warmer Wang, 1989;Wang et al., 2013;Ye et al., 2014). Most disjunct genera and species between HHM and TW are herbaceous (Chen, Ying, & Lu, 2012;Ye et al., 2012). However, few studies on the HHM-TW disjunction focused on herbaceous plants, so it remains unclear whether similar mechanisms and timing shaped the disjunct patterns of woody and herbaceous plants between HHM and TW.
Recently, species distribution modeling (SDM) has been used to infer range shifts of plants and animals in response to the Quaternary climate oscillations (Jakob, Ihlow, & Blattner, 2007;Luo et al., 2016;Peterson, Ammann, & Diniz-Filho, 2013). In general, most temperate species retracted into smaller refugial areas during glacial periods and would experience range expansion during interglacial periods (Bueno et al., 2016;Hewitt, 1996Hewitt, , 1999Hewitt, , 2000. However, several recent studies suggested a different biogeographic history for cold-adapted coniferous species that inhabit high mountains, with range expansion during glacial periods and contraction during interglacial periods (Liu et al., 2013;Shao & Xiang, 2015;Tian, López-Pujol, Wang, Ge, & Zhang, 2010;Wang et al., 2017). Compared with lowland species, mountainous species might be more likely to survive during climatic oscillations, because short-distance vertical migrations allow them to track their optimal temperature niche (Hoorn, Mosbrugger, Mulch, & Antonelli, 2013;Sandel et al., 2011). Nevertheless, continued global warming in the near future may be a major threat to cold-adapted species, especially those that are already confined to mountaintops or islands, because migration to higher elevations may be not possible or range shifts are not fast enough to track suitable climate (Wiens, 2016). How organisms responded to the Quaternary climatic oscillations and predictions of response to future climate can likely be inferred with short-lived herbaceous plants, which have far more life cycles within a given time period and thus may respond more quickly to rapidly changing environments (Comes & Kadereit, 1998;Faye et al., 2016;Zanne et al., 2014).
As an herbaceous genus of Caprifoliaceae s.l. (The Angiosperm Phylogeny Group, 2016), Triplostegia Wall. ex DC. has a disjunct distribution between HHM and TW (see Figure S1), with only a few sporadic records in Indonesia (Hong, Ma, & Barrie, 2011;Zhang, Chen, Li, Chen, & Tang, 2003). Triplostegia is sensitive to habitat and climatic changes and requires a specific forest habitat that has fertile soil and an altitude ranging from 1,500 to 4,000 m; these features make this genus an ideal model for studying formation of HHM-TW disjunctions and how herbs respond to Quaternary climate oscillations and future global warming.
By sequencing five chloroplast DNA (cpDNA) regions and the ribosomal nuclear internal transcribed spacer (ITS) for 397 individuals from 48 populations, estimating divergence time, reconstructing ancestral areas, and conducting SDM, we (1) investigated the genetic structure and identified possible refugia for Triplostegia; (2) reconstructed the biogeographic history of Triplostegia and explored how the HHM-TW disjunction formed; and (3) predicted the future climatically suitable areas of the cold-adapted Triplostegia under a global warming scenario.

| Population sampling
Triplostegia includes two species: T. glandulifera Wall. ex DC. and T. grandiflora Gagnep. (Peng, Tobe, & Takahashi, 1995;Pyck & Smets, 2004;Zhang, Chen, Li, Chen, & Tang, 2003). Triplostegia grandiflora is endemic to China (confined to Yunnan and Sichuan Provinces), and its geographical range overlaps with that of T. glandulifera (Hong, Ma, & Barrie, 2011). Both species are perennial herbs only with minor differences in flower size. We herein treated Triplostegia as a monotypic genus based on the evidence from pollen morphology and phylogenomic analysis of nine complete chloroplast genomes (Niu et al., unpublished).
From September 2011 to August 2014, we collected 397 individuals from 48 natural populations, which covered the entire geographical range of Triplostegia in China. Generally, 10-12 individuals were arbitrarily selected from a population, whereas three to six individuals were collected from populations with few individuals (Table S1). For each individual, leaves were collected and dried in silica gel for DNA extraction. Voucher specimens were deposited at the Herbarium of Institute of Botany, Chinese Academy of Sciences (PE).
All haplotype sequences identified in this study were deposited in GenBank (https://www.ncbi.nlm.nih.gov/; see Table S3 for accession numbers). Genealogical relationships among haplotypes were investigated using the median-joining network method as implemented in NETWORK 5.0.0 (Bandelt, Forster, & Röhl, 1999).
The average gene diversity within populations (H S ), total gene diversity (H T ), and the coefficients of differentiation (G ST , a population differentiation estimate based solely on haplotype frequencies, and N ST, a parameter that considers both haplotype frequencies and their genetic divergence) were estimated using PERMUT 1.0.0 (Pons & Petit, 1996) with 1,000 random permutations. When N ST is larger than G ST , closely related haplotypes occur more often in the same area than less closely related haplotypes, which indicates strong phylogeographic structure.
Spatial genetic structure of cpDNA and ITS haplotypes was analyzed by SAMOVA 1.0 (Dupanloup, Schneider, & Excoffier, 2002). This program uses a simulated annealing approach to identify groups of populations (K) that are geographically homogeneous and maximally differentiated from each other. To maximize the proportion of total genetic variance observed between groups (Φ CT ), 1,000 iterations were run for K = (2, …, 8) from each of 100 random initial conditions.
To further quantify genetic differentiation among populations and groups, the hierarchical analysis of molecular variance (AMOVA) was estimated in ARLEQUIN 3.5.1 (Excoffier & Lischer, 2010) with significance tests based on 1,000 permutations.

| Demographic reconstruction
Mismatch distribution analysis with a spatial expansion model (Rogers & Harpending, 1992) was conducted in ARLEQUIN to infer the historical demography of Triplostegia based on both cpDNA and ITS datasets. The goodness of fit based on the sum of squared deviation (SSD) and Harpending's raggedness index (H Rag ) was tested using 1,000 parametric bootstrap replicates. If the expansion model was not rejected, the formula t = τ/2u was used to estimate the expansion time (t), where u = μkg [μ, substitution rate, substitutions per site per year, subs/site/year; k, average sequence length of cpDNA; g, generation time in years (i.e., age of first reproduction; approximated as 1 year, J. F. Ye, personal observation)]. In addition, statistical significance of population expansion was further examined with two neutrality tests, Tajima's D (Tajima, 1989) and Fu's F S (Fu, 1997), as implemented in DNASP, in which significantly negative values indicate population expansion.

| Phylogenetic relationships among haplotypes
Phylogenetic analyses for cpDNA and ITS haplotypes were performed using the maximum likelihood (ML), Bayesian inference (BI), and neighbor-joining (NJ) methods, with gaps treated as missing data. The ML analysis was conducted using RAxML 8.2.8 (Stamatakis, 2006(Stamatakis, , 2014  generations. The average standard deviation of the split frequencies approached 0.01, which indicated that two runs converged to a stationary distribution. After discarding the first 25% trees as burn-in, a 50% majority-rule consensus tree was constructed from the remaining trees to estimate posterior probabilities (PP). NJ analysis for ITS was carried out in MEGA 6.0.0 (Tamura, Stecher, Peterson, Filipski, & Kumar, 2013) with 1,000 BS replicates.

| Divergence time estimation
To estimate the crown age of Triplostegia, we downloaded sequences of three cpDNA regions (trnH-psbA, trnL-F, and trnS-trnG) for 41 Dipsacales species (including five Adoxaceae species as outgroups) from GenBank (Table S4) that represent all major lineages of Dipsacales (Bell & Donoghue, 2005;Carlson, Linder, & Donoghue, 2012;Wang, Landrein et al., 2015). Sequences were assembled, with 16 Triplostegia haplotypes identified from the same three cpDNA regions, to represent all 20 cpDNA haplotypes identified in the phylogeographic analyses. Bayesian searches for tree topologies and node ages of this cpDNA dataset were conducted in BEAST 1.7.3 (Drummond & Rambaut, 2007), which was run on the CIPRES Science Gateway Portal. An uncorrelated lognormal relaxed clock (Drummond, Nicholls, Rodrigo, & Solomon, 2002) was applied with the GTR + G substitution model and a Yule process tree prior. A secondary calibration point (52.8-89.0 Ma; Magallón, Gómez-Acevedo, Sánchez-Reyes, & Hernández-Hernández, 2015) was used to constrain the crown group of Dipsacales that assumed a normal distribution with a mean (± SD) of 70.94 ± 11 Ma (node 1 in Figure S2).
Three fossil dates were used to assign minimum age constraints on three internal nodes of Dipsacales with a lognormal distribution (nodes 2-4 in Figure S2). The earliest fossil record of Viburnum from the late Paleocene to early Eocene (c. 56 Ma; Baskin et al., 2006;Wing et al., 1995) was used to constrain the crown group of Adoxaceae (node 2 in Figure S2; mean = 0, SD = 1.0, offset = 56 Ma) following Moore and Donoghue (2007). Fossil seeds of Weigela have been reported from the Miocene and Pliocene in Poland (Łańcucka-Środoniowa, 1967), the Miocene in Mammoth Mts., Eastern Russia, the Oligocene and Miocene of Western Siberia (Dorofeev, 1963), and the Miocene in Denmark (Friis, 1985). We thus constrained divergence between Weigela and its sister Diervilla to 23 Ma (node 3 in Figure   S2; mean = 0, SD = 1.0, offset = 23 Ma) following Wang, Landrein et al. (2015). Finally, we used the fossil fruits of Diplodipelta from the late Eocene Florissant flora of Colorado (36-35 Ma;Manchester & Donoghue, 1995) to constrain the stem age of Dipelta (node 4 in Figure S2; mean = 0, SD = 1.0, offset = 36 Ma) following Bell and Donoghue (2005) and Wang, Landrein et al. (2015). We assessed estimate robustness by sequentially removing each calibration while retaining the other three calibrations (Wahlberg, Wheat, & Peña, 2013).
To determine the intraspecific node ages within Triplostegia, we applied a similar BEAST analysis to a dataset that included all 20 haplotypes of five cpDNA regions, with Dipsacus asper as the out-group.
We used the median stem age and crown age of Triplostegia obtained from the Dipsacales chronogram to constrain the root and crown age of Triplostegia, respectively (nodes A and B in Figure S2). The normal distribution and a constant-size coalescent tree prior were applied.
MCMC runs were performed, each of 10 8 generations, with sampling every 10 4 generations, following a burn-in of the initial 10% of cycles.

| Ancestral area reconstruction
We used the Bayesian Binary MCMC (BBM) method implemented in Plants in Taiwan (SDNPT, http://www.hast.biodiv.tw/Announce/ newsC.aspx). The out-group (Dipsacus asper) is found in both SW and CC, so we set its distribution as SW/CC. The number of maximum areas at each node was set to three. To account for phylogenetic uncertainty, 1,000 trees from BEAST were randomly chosen for BBM analysis. We set the root distribution to null, applied 10 MCMC chains with the JC + G model, and sampled the posterior distribution every 100 generations for 10 6 generations.

| Species distribution modeling
We used MaxEnt species distribution modeling (henceafter MaxEnt) to predict past, current, and future climatically suitable areas for Triplostegia (Phillips, Anderson, Dudík, Schapire, & Blair, 2017). These nine climatic variables (i.e., the maximum, minimum, mean and standard deviation of temperature and precipitation) can influence the distribution and physiological performance of plant species. Altitude is an important factor to predict suitable habitats for Triplostegia, but we did not take it into consideration due to its close relationship with temperature and precipitation.
Current distribution of Triplostegia was developed using the same nine bioclimatic variables for current climate  We corrected all the occurrence data to reduce the negative effect of sampling bias on the performance of the SDMs. For records from GBIF, CVH, and SDNPT, we excluded the duplicated entries within a 2.5 arc-minute pixel (4.3 km at the equator) and omitted entries with implausible geographical coordinates using ArcGIS 10.2 (ESRI Inc., Redlands, CA, USA; http://www.esri.com/software/arcgis/arcgis-for-desktop) and Google Earth (http://earth.google.com/; Gueta & Carmel, 2016;Maldonado et al., 2015). We also checked voucher photos for records from CVH and SDNPT to make sure that the species are correctly identified. Finally, a total of 250 cleaned occurrence records were obtained from GBIF, CVH, SDNPT, and our own collections (Table S1 and Figure S3), which well covered the distribution range of Triplostegia (Hong, Ma, & Barrie, 2011).
As suggested by previous studies (Halvorsen et al., 2016;Martinez-Freiria, Velo-Antón, & Brito, 2015;Merow, Smith, & Silander, 2013;Phillips et al., 2017;Radosavljevic, Anderson, & Araújo, 2014), our modeling was set as follows: (1) the regularization multiplier level was at two to produce smooth and general response curves that represent a biologically realistic model; (2) cloglog output format was used to evaluate distribution probability of Triplostegia; (3) ten independent replicates of cross-validation procedures for model testing; (4) the maximum number of background points was set to 10 4 ; and (5) jackknife analyses of the regularized gain using training data were conducted to elucidate the importance of each predictor (Fawcett, 2006). The other settings followed the suggestions of Merow, Smith, and Silander (2013). No obvious differences were detected among the predictions based on the three GCMs ( Figure S4). To decrease the prediction uncertainty, mean consensus for the pixels of distribution probabilities modeled by three GCMs was calculated to assess the distributions of Triplostegia during LGM and the future (RCPs 4.5 and 8.5 for both 2050s and 2070s). We evaluated the model performance by calculating the area under the receiver operating characteristic (ROC) curve (AUC). AUC was used to assess the accuracy of predictions, and models with AUC values larger than 0.7 were considered useful for our study (Phillips & Dudík, 2008). However, using the traditional AUC may be not sufficient to assess the performance of SDMs. We herein used a binomial test (based on the training omission rate) for model validation based on the occurrence data only. The training omission rate is the proportion of training occurrence records among pixels of predicted absences (Phillips & Dudík, 2008;Phillips, Anderson, & Schapire, 2006). One-sided test was used for the null hypothesis, and the test points are predicted no better than those by a random prediction (Phillips & Dudík, 2008;Phillips, Anderson, & Schapire, 2006). We used the 11 common thresholds defaulted by MaxEnt to define binomial probabilities (Phillips, Anderson, & Schapire, 2006). The models with <17% training omission rate were useful for our study (Phillips, Anderson, & Schapire, 2006).
H d and π for the cpDNA dataset were estimated to be 0.87 and 0.24 × 10 −2 , respectively (Table S1). H T and H S were 0.88 and 0.18, respectively. The coefficients of differentiation measured over the 48 populations were G ST = 0.80 and N ST = 0.96. Permutation test revealed that N ST was significantly higher than G ST (p < .05), which indicates strong phylogeographic structure for the Triplostegia cpDNA haplotypes. The highest levels of cpDNA sequence diversity (π) were recorded within the SW group and the lowest levels in the CC group (Table S1). Differences in cpDNA sequences were highly significant among the three geographical regions (F CT = 0.85, p < .001, AMOVA;  Table 1). Differences in cpDNA sequences were significant between every two of the three geographical regions (Table S5). The highest differentiation index was found between SW and TW (F CT = 0.88 *** ), and the lowest between CC and TW (F CT = 0.68 *** ).
The matrix for aligned ITS had 573 bp with 36 polymorphic sites and 14 ITS haplotypes (see Table S3 for accession numbers).
H d and π for ITS were estimated to be 0.67 and 0.011, respectively (Table S1). Network and geographical distribution analysis based on ITS also recognized three groups in SW, CC, and TW ( Figure 2).
The highest level of ITS π was recorded within the SW group and the lowest in the CC group (Table S1). Differences in ITS sequences were highly significant among the three geographical regions (F CT = 0.89, p < .001, AMOVA; Table 1), between populations within regions (F SC = 0.98, p < .001), and among all populations (F ST = 0.99, p < .001). Differences in ITS sequences were significant between each pair of geographical regions (Table S6). The highest differentiation index was observed between SW and TW (F CT = 0.91 *** ) and the lowest between CC and TW (F CT = 0.84 *** ). Analysis of the ITS variation across all populations revealed that the total gene diversity (H T = 0.67) was much higher than the within-population diversity (H S = 0.025). The coefficients of differentiation measured over the 48 populations were G ST = 0.96 and N ST = 0.99. Permutation test revealed that N ST was significantly higher than G ST (p < .05), which indicated a strong phylogeographic structure for ITS haplotypes.
SAMOVA indicated that values of Φ CT reached a plateau at 0.88 for three groups (K = 3; Figure 2c), which corresponded to the three clades in SW, CC, and TW.

| Phylogenetic relationships and divergence times
Both ML and BI analyses based on the 20 cpDNA haplotypes supported three lineages within Triplostegia that corresponded to their distributions in SW, CC, and TW ( Figure S5). However, incongruence existed regarding the relationships among three lineages between the BI (SW(CC, TW)) and ML (TW(SW,CC)) analyses, although the support values were not high. The ML, BI, and NJ analyses based on ITS retrieved the same topology (SW(CC, TW)) as the BI tree based on cpDNA ( Figure S5).
The chronogram for Dipsacales derived from cpDNA supported the monophyly of the 16 Triplostegia haplotypes (PP = 1.00; Figure   S2). The stem age of Triplostegia was estimated to be c. 48 T A B L E 1 Molecular variances of cpDNA and ITS for 48 Triplostegia populations F I G U R E 1 (a) Geographical distribution of the cpDNA haplotypes detected in Triplostegia (see Table S1 for population codes) and (b) network of genealogical relationships among 20 haplotypes. Each circle represents a single haplotype with sizes in proportion to frequency. Black dots in network represent missing haplotypes. Potential interglacial refugia recognized by this study are encircled by dashed lines showed that the coalescent time of all 20 haplotypes was during the late Miocene (6.79 Ma, 95% HPD: 2.50-13.6 Ma; node I in Figure   S5)  Figure S5). For this chronogram, BEAST provided an average substitution rate of 0.49 × 10 −9 subs/site/year, which is similar to that of Tetrastigma hemsleyanum derived from three F I G U R E 2 Analysis based on ITS dataset. (a) Geographical distribution of haplotypes in Triplostegia (see Table S1 for population codes); (b) network of genealogical relationships between 14 haplotypes. Each circle represents a single haplotype sized in proportion to its frequency. Black dots in the network represent missing haplotypes; (c) results of SAMOVA (the red point indicated the optimized number of groups K); and (d) neighbor-joining tree without out-group. Potential interglacial refugia recognized by this study are encircled by dashed lines cpDNA markers (petL-psbE, trnK-matK, and rbcL; 0.50 × 10 −9 subs/ site/year; Wang, Jiang et al., 2015).

| Demographical analyses
The overall observed mismatch distribution of pairwise nucleotide differences was not unimodal and rejected the hypothesis that there was recent demographic population growth in Triplostegia (Figures 3a and S6) Figure S6). Based on the value of τ and assuming a substitution rate of 0.49 × 10 −9 subs/site/year for cpDNA, the demographic expansion in SW populations was inferred to have occurred c. 0.33 Ma (Table 2).

| Ancestral area reconstruction
Based on the topology of the intraspecific chronogram ( Figure S5), BBM analysis supported an ancestral distribution in SW (node I in

| Species distribution modeling
The training and test AUC values of the models were larger than 0.9, and all the training omission rates were <17%, indicating a good performance of the predictive model. The jackknife results showed that temperature was the most influential factor.
During the LIG, high climatic suitability areas for Triplostegia were more restricted than its current range in southwestern China (e.g., the

| Genetic diversity and potential refugia
Results from PERMUT (N ST > G ST , p < .05), AMOVA (Table 1), and SAMOVA ( Figure 2c) revealed strong genetic structure for the Triplostegia populations from SW, CC, and TW (Figures 1 and 2). The three distinct haplotype groups could be distinguished by a strong underlying phylogenetic hierarchy that supports three major clades in Triplostegia that correspond to their distributions in SW, CC, and TW (Figures 2d and S5). The strong genetic structure among Triplostegia  Strong geography-haplotype-correlated genetic structures were also detected in birds (Qu et al., 2015;Wang et al., 2013), insects (Ye et al., 2014), and plants (Gao et al., 2007;Meng et al., 2017), which highlight the importance of topographic complexity in promoting species differentiation both by increasing habitat diversity and limiting gene flow between elevation-restricted populations (Hoorn, Mosbrugger, Mulch, & Antonelli, 2013;Verboom, Bergh, Haiden, Hoffmann, & Britton, 2015).
Potential refugia are generally predicted to have a widespread ancestral haplotype, and other mutationally derived and unique haplotypes (Crandall & Templeton, 1993). According to coalescent theory, most ancient haplotypes should be located at the interior nodes of a haplotype network, whereas the most recent haplotypes should be at the tips (Posada & Crandall, 2001). In the cpDNA network of  (Peng, Wan, & Luo, 2000), which might have provided suitable habitats for Triplostegia during the Pleistocene climate oscillations.

| Biogeographic history
Our analyses revealed that coalescence of all Triplostegia cpDNA haplotypes is most likely to have occurred in SW, and diversification  (current, 1960-1990); and (d) 2070s (2061-2080) with the representative concentration pathways (RCP) of 4.5. Climatic suitability increases with color from blue to red. Resolution for the potential distribution map is 2.5 arc-minutes initiated during the late Miocene (c. 6.79 Ma; 95% HPD: 2.50-13.6 Ma; node I in Figures 4 and S5). The late Miocene has been recognized as an important period of diversification for both woody and herbaceous plants (e.g., Heterobalanus, Meng et al., 2017;Cercidiphyllum, Qi et al., 2012;Tetracentron, Sun et al., 2014;Tetrastigma, Wang, Jiang et al., 2015). Differentiation within Triplostegia might have been initiated and driven by both uplift of the Hengduan Mts. and global cooling after the late Miocene (Nagalingum et al., 2011;Zachos, Dickens, & Zeebe, 2008). The Hengduan Mts., at the southeastern margin of the QTP, are proposed to have experienced rapid uplift during the late Miocene (Clark et al., 2005;Kirby et al., 2002). By investigating the evolutionary histories of multiple plant groups, Xing and Ree (2017) detected an increase in the rate of in situ diversification in the Hengduan Mts. during the late Miocene (c. 8.00 Ma). The SW origin of Triplostegia is further supported by our phylogenetic reconstruction based on nine complete chloroplast genomes (Niu et al., unpub-lished) that represented populations from SW, CC, and TW, which strongly supported (PP = 1.00; BS = 100%) initial divergence of the SW clade, and CC and TW populations grouped together ( Figure   S8).  (Favre et al., 2015;Ge et al., 2013;Mulch & Chamberlain, 2006). The absence of shared haplotypes between SW and CC indicates a lack of gene flow between the two regions that was probably due to habitat isolation caused by elevation gradients (average elevation for populations in SW is c. 900 m higher than those from CC; Table S1). Morphologically, CC populations also differed from SW populations based on the presence of a serrate leaf blade (vs. pinnatifid leaf blade in SW; Figure 4) Gao et al., 2007) was congruently inferred to be mid-Pliocene, which was shortly after Taiwan attained its modern form. The lowest differentiation index between CC and TW (Tables S5 and S6) and the more similar leaf blades (Figure 4) further support the idea that TW populations are descendants of CC.
Our SDM analysis provides a detailed picture of the last glacial cycle (Figure 5a,b) and glimpses of the preceding cycles. Based on the paleoclimate record, temperature in each glacial period was similar (Bloom, 2010

| Conservation implications
Identifying areas with stable niches and understanding their importance in determining the current distribution of species represent a pivotal task for biodiversity conservation (Marta, Lacasella, Gratton, Cesaroni, & Sbordoni, 2016). The predicted high climatic suitability areas of Triplostegia in 2070s are more restricted than that of the present, especially in southern HHM (Figure 5c,d). Additionally, areas with high climatic suitability in Mainland China are predicted to show a slightly northward migration, whereas climatic suitability in Taiwan may suffer from declining ( Figure 5d). The potential reason may be that Mainland China has very northern mountain habitats that are suitable for Triplostegia; however, Taiwan has neither northern cold areas nor higher mountains for Triplostegia dispersal with continued global warming ( Figure 5d). Moreover, global mean annual temperatures are likely to rise by an additional 1-4°C by the year 2100 (Stocker et al., 2013). It is undoubtedly a big challenge for the survival of many organisms, especially cold-adapted taxa that inhabit high mountains, such as Triplostegia, which may experience local extinction at the warm edge, and their suitable habitat distribution in the near future will be as limited as or more restricted than during the LIG (Beever, Ray, Wilkening, Brussard, & Mote, 2011;Wiens, 2016).
Species distribution modeling only uses species occurrence data and environmental variables to predict the potential geographical distribution of species (Guisan & Thuiller, 2005). However, climate change is not the sole factor affecting species' habitat suitability. Human activities (e.g., overharvesting, urbanization, fire suppression, and nitrogen pollution) also have profound effects on the distribution of organisms (MacDougall, McCann, Gellner, & Turkington, 2013). Although we predicted a northward range shift under continuing global warming, the potential habitats may be no longer suitable for survival because of human activities. Indeed, we failed to discover Triplostegia in numerous places with detailed specimen records. Triplostegia represents series of organisms that are already confined to mountaintops (QTP) and/or islands (TW). Therefore, conservation actions are needed to protect these taxa, which are sensitive to climate change and habitat disturbance.

| CONCLUSION
In this study, we explored how geology and climate together influenced the evolutionary history of Triplostegia, which is an ecologically sensi-