Geological and Climatic Factors Affect the Population Genetic Connectivity in Mirabilis himalaica (Nyctaginaceae): Insight From Phylogeography and Dispersal Corridors in the Himalaya-Hengduan Biodiversity Hotspot

Abstract

The genetic architecture within a species in the Himalaya-Hengduan Mountains (HHM) region was considered as the consolidated consequence of historical orogenesis and climatic oscillations. The visualization of dispersal corridors as the function of population genetic connectivity became crucial to elucidate the spatiotemporal dynamics of organisms. However, geodiversity and physical barriers created by paleo geo-climatic events acted vigorously to impact notable alterations in the phylogeographic pattern and dispersal corridors. Therefore, to achieve detailed phylogeography, locate dispersal corridors and estimate genetic connectivity, we integrated phylogeography with species distribution modelling and least cost path of Mirabilis himalaica (Edgew.) Heimerl in the HHM. We amplified four cpDNA regions (petL-psbE, rps16-trnK, rps16 intron, trnS-trnG), and a low copy nuclear gene (G3pdh) from 241 individuals of 29 populations. SAMOVA, genealogical relationships, and phylogenetic analysis revealed four spatially structured phylogroups for M. himalaica with the onset of diversification in late Pliocene (c. 3.64 Ma). No recent demographic growth was supported by results of neutrality tests, mismatch distribution analysis and Bayesian skyline plot. Paleo-distribution modelling revealed the range dynamics of M. himalaica to be highly sensitive to geo-climatic change with limited long-distance dispersal ability and potential evolutionary adaptation. Furthermore, river drainage systems, valleys and mountain gorges were identified as the corridors for population genetic connectivity among the populations. It is concluded that recent intense mountain uplift and subsequent climatic alterations including monsoonal changes since Pliocene or early Pleistocene formulated fragmented habitats and diverse ecology that governed the habitat connectivity, evolutionary and demographic history of M. himalaica. The integrative genetic and geospatial method would bring new implications for the evolutionary process and conservation priority of HHM endemic species.

Introduction

The Himalaya and the Hengduan Mountains (HHM) are one of the main biodiversity hotspots of the Northern Hemisphere formed as a result of the collision between the Indian and Eurasian plate, started c. 55–50 million years ago (Ma) (Yin and Harrison, 2000; Dupont-Nivet et al., 2010). Mulch and Chamberlain (2006) suggests that the orogeny of the Hengduan Mountains occurred as a final propagation of the uplift after late Miocene (c. 10 Ma). In addition, recent violent Qinghai-Tibetan Plateau (QTP) uplifts (c. 3.6–1.5 Ma), and associated intensification of the Asian monsoon (c. 2.6–3.6 Ma) caused dramatic changes in the biodiversity of HHM and QTP regions alongside alterations in regional landforms and environmental heterogeneity (Li et al., 1979; Li and Fang, 1999; An et al., 2001).

The incredible biodiversity and genetic architecture in the HHM region were considered as the aftereffect of historical orogenesis independently or in combination with climatic oscillations (Avise, 2009; Favre et al., 2015). Moreover, the geological and climatic factors enormously influence the spatiotemporal pattern of temperate plant species, including historical dispersal and gene flow among the populations (Yu et al., 2015). Such orogenic uplift of mountains results in the formation of the geographical barriers, leading to fragmented habitat, loss of dispersal corridors, and restricted gene flow among populations; such factors furnish circumstances for divergence owing to genetic drift and natural selection (Liu et al., 2006; Meng et al., 2017). Historical dispersal process is likely to be a key factor to determine the current spatial population structure of species affected by the geological and climatic alterations (Hewitt, 2000). Therefore, it is crucial to locate historical dispersal for a better understanding of the demographic and evolutionary history of species, as an outcome of paleo geo-climatic events (Brown and Yoder, 2015).

Phylogeographic approach integrated with ensemble species distribution modelling (SDM) analysis provides comprehension of how paleo-environmental alterations in landscape and climate have governed species distributions and demographical history (Avise, 2000; Wang and Yan, 2014). In addition, landscape genetics has been widely used for modelling the dispersal corridors of species (Wang et al., 2009; Yu et al., 2017). Least cost path (LCP) uses landscape genetic approach coupled with species distribution models and population genetic data to recognize population genetic connectivity in a spatially explicit framework (Chan et al., 2011; Yu et al., 2015). LCP contributes to the most suitable habitat and fewest movement barriers providing the best theoretical route for dispersal in organisms (Larkin et al., 2004).

Over the last decade, a number of phylogeographic studies have investigated the link between organismic evolution of plant species and the geologic uplift associated with climate changes in the HHM and QTP (Qiu et al., 2011; Jian et al., 2015; Luo et al., 2016). For instance, phylogeographic studies on Primula tibetica related to the effects of Quaternary climatic oscillations in QTP (Ren et al., 2017a), Metrocoris sichuanensis concerned with geological effects in Sichuan Basin (Ye et al., 2018), and Allium section Sikkimensia linked with Hengduan Mountains massif uplift (Xie et al., 2019). They have demonstrated that the dramatic uplift independently or in combination with Pleistocene glaciations influenced their patterns of genetic variation. However, few studies have compared the relative significance of geo-climatic mechanisms influencing the historical dispersal and the population genetic connectivity in this region (Rana et al., 2019).

The present study focuses on the HHM endemic Mirabilis himalaica (Nyctaginaceae) that distributes from the Western Himalaya to the Hengduan Mountains that stretch from N India, WC Nepal, Bhutan, and S Xizang, SW Gansu, N Sichuan, NW Yunnan in the HHM region (Lu et al., 2003; Wang et al., 2019). The genus Mirabilis L. constitutes ~60 spp., predominantly in temperate and tropical North America and South America and single species in Asia (M. himalaica; Lu et al., 2003; Spellenberg, 2003; Wang et al., 2019). Mirabilis himalaica grows in thickets, grasslands, dry and warm river valleys between 1,400–3,800 m (Lu et al., 2003; Wang et al., 2019). For the present work, we hypothesize that the geological events to be the most important driver, if the divergence occurs before mid-Pleistocene with relatively old divergence time. We further presumed that the Pleistocene climate along with geological events greatly influenced the lineage colonization and spatiotemporal distribution of the M. himalaica. Therefore, we set the following primary goals for our research: 1) to identify genetic diversity and phylogeographic structure; 2) to date the divergence time between lineages and locate population genetic connectivity during the late Quaternary in the HHM region; 3) and to elucidate how paleo geological change and climatic oscillations have affected comprehensive evolutionary and demographic history of M. himalaica. This study represents the integrative approach of maternally inherited (cpDNA) and biparentally inherited (G3pdh) population genetic data in combination with the SDM and LCP approach throughout the distribution range of species to outreach the role played by the historical processes on the present-day spatial population structure of organisms in the HHM.

Materials and Methods

Sampling, DNA Extraction, PCR Amplification, and Sequencing

In total, 241 individuals, from 29 populations covering the possible geographical range i.e. from the Western Himalaya to the Hengduan Mountains were sampled for the phylogeographic study of Mirabilis himalaica (Figure 1 and Supplementary Table S1). The fresh sampled leaves of 6–10 individuals from each population (Table 1) were at least 30 m distance apart and dried in silica gel. Voucher specimens were deposited at National Herbarium and Plant Laboratories (KATH), Nepal, and Kunming Institute of Botany (KUN), China. Total genomic DNA was extracted from c. 30 mg silica dried leaves using a DNAsecure Plant Kit (Tiangen Biotechnology Co. Ltd., Beijing, China) following the manufacturer's protocol. After preliminary screening, we amplified four cpDNA regions i.e. petL-psbE, rps16-trnK, rps16 intron, and trnS-trnG, and a low copy nuclear gene i.e. G3pdh for each individual. PCRs were conducted in a 30 µl containing 2 µl DNA template (10–50 ng/µl), 15 µl 2x Taq Plus Master Mix with dye (Tiangen Biotech.), 1 µl 10 µM of each primer (see Supplementary Methods S1 for detail).

Figure 1

The sequences of four cpDNA regions were concatenated, revised manually and aligned in Geneious 7.0.2 (Kearse et al., 2012). Chromatograms of the G3pdh with “double peaks” at polymorphic sites were further analyzed by inferring the identity of haplotypes within a heterogeneous individual through haplotype subtraction (Clark, 1990; Christe et al., 2014) using PHASE algorithm in DnaSP 5.0 (Librado and Rozas, 2009). The unique cpDNA/G3pdh haplotypes within populations were identified using DnaSP 5.0 by considering the polymorphic sites only. Newly generated haplotype sequence data have been deposited in GenBank (petL-psbE: MK792906–MK792916; rps16-trnK: MK792917–MK792925; trnS-trnG: MK792926–MK792935; rps16 intron: MK803108–MK803115; G3pdh: MK803090–MK803107).

Genetic Polymorphism and Structure Analyses

Unique chlorotypes of cpDNA and haplotypes of G3pdh were determined in DnaSP 5.0 (Librado and Rozas, 2009). Geographical distributions of chlorotypes/haplotypes were plotted using ArcGIS 10.2.1 (ESRI, Inc., Redlands, CA, USA). To quantify the level of genetic variation, total haplotype diversity (H_T) and average within-population diversity (H_S) (Nei, 1975) were calculated. G_ST and N_ST were used to estimate differentiation between populations (Nei, 1975). N_ST was compared to G_ST using U-statistics; an observed value of N_ST > G_ST generally indicates the presence of phylogeographical structure (Pons and Petit, 1996). We computed all these parameters employing Haplonst (Pons and Petit, 1996). DnaSP was used to analyze the genetic diversity parameters, including the haplotype diversity (H_d; Nei and Tajima, 1981) and nucleotide diversity (π; Nei and Li, 1979).

Spatial analysis of molecular variance (Samova v2.0; Dupanloup et al., 2002), was implemented to define the number of groups of populations (K). We also performed an analysis of molecular variance (AMOVA) to examine the genetic variation using Arlequin v3.5 (Excoffier and Lischer, 2010). The F-statistic (F_ST/F_SC/F_CT) was calculated, and significance was tested for overall as well as regional populations. “Isolation by distance” (IBD) was evaluated by regressing the net nucleotide divergence between populations (D_A) in contrary to their geographical distance, applying Mantel's (1967) test with 999 permutations in GenAlEx 6.5 (Peakall and Smouse, 2012). Network v5.0.0 (Bandelt et al., 1999; http://www.fluxus-engineering.com) was employed to construct median-joining networks of the cpDNA and G3pdh sequences to assess their genealogical relationships. All variable sites were included and weighted equally.

Phylogenetic Analyses and Divergence Dating

Phylogenetic relationships among chlorotypes of cpDNA along with sequences of three closest species from Mirabilis L. (M. jalapa, EF079612; M. multiflora, EF079603; M. albida, EF079602/KR014118) as outgroups were constructed by Bayesian inference (BI) in MrBayes v3.2 (Ronquist et al., 2012) and Beast v1.8.4 (Drummond and Rambaut, 2007; Drummond et al., 2012) under the GTR+I nucleotide substitution model selected by Akaike information criterion (AIC) in jModelTest v2.1.6 (Guindon and Gascuel, 2003; Posada and Buckley, 2004) (see Supplementary Methods S1 for detail parameters). Beast v1.8.4 was employed to estimate the temporal intraspecific divergence times (crown ages) of chlorotypes/haplotypes. We assumed a strict clock (P = 0.85; i.e. P > 0.05), based on a likelihood-ratio test (Felsenstein, 1988) in Paup* 4.0b10 (Swofford, 2002) and constant population size for the coalescent tree prior to the distribution of divergence times. We used two secondary calibration points (Figure 2) to estimate the lineage divergence times based on Wang et al., 2019 (A: 13.13 ± 4.2 Ma; 95% highest posterior density [HPD] intervals: 6.91–20.62 Ma; i.e. a crown age for Mirabilis, and B: 5.22 ± 1.7 Ma; 95% 2.53–8.18 Ma i.e. divergence time of M. himalaica from its North American counterparts; Figure 2) (see Supplementary Methods S1 for detail parameters). Three independent runs of 20 million generations were carried out, with number of chains = 4, and sampling every 1,000 generations, where first 20% of the trees were discarded as burn-in. TreeAnnotator v1.8.4 (Drummond and Rambaut, 2007; Drummond et al., 2012) was used to obtain a maximum-credibility tree and FigTree v1.4.0 (Rambaut, 2012) was employed to view the resulting tree.

Figure 2

Historical Demographic Analyses

Tajima's (1989)D and Fu's (1997)F_S of neutrality tests were performed using Arlequin v3.5 (Excoffier and Lischer, 2010), with 1,000 simulated samples. Pairwise mismatch distribution analysis (MDA; Rogers and Harpending, 1992) was performed in Arlequin v3.5 to detect demographic expansion events of M. himalaica. The MDA compared the observed frequency distribution of pairwise differences among haplotypes with their theoretical distribution expected under the ‘pure population growth’ and ‘spatial expansion’ models, respectively. With the sum of squared deviations (SSD) and Harpending's (1994) raggedness index (H_Rag) ‘goodness of fit’ was tested, using 1,000 parametric bootstrap replicates.

In addition, Bayesian skyline plot (BSP: Drummond et al., 2005) generated in Beast v1.8.4 was used to reconstruct the effective population size fluctuations since the time of the most recent common ancestor for each marker and the combined markers (cpDNA and G3pdh). The MCMC chains (nchains = 4) were run for 20 million generations sampling every 100 generations, with effective sample size (ESS) greater than 200. We used the same settings for three independent runs to ensure the consistency of the results. The demographic history through time was reconstructed using TRACER v1.7.1 (Rambaut et al., 2018).

Ensemble Species Distribution Modelling and Visualizing Dispersal Corridors

An ensemble of 'species distribution modelling (SDM; Guisan and Zimmermann, 2000) was carried out using “Biomod 2” (Thuiller et al., 2014) package implemented in R-programming language (R v3.4.1; R Development Core Team, 2016) for past (Last Interglacial, LIG; 120–140 Ka and Last Glacial Maximum, LGM; c. 22 Ka) (Otto-Bliesner et al., 2006), current (1990–2000) and future (2070; RCP 4.5). The model used 72 spatially filtered occurrence points to prevent from spatial autocorrelation that were checked in the 4-min grid (Rana et al., 2019). The eight predictive bioclimatic variables were selected based on iterative calculations of Variance Inflation Factors (VIF < 10; Fox and Weisberg, 2011) (Table 2), Pearson correlation (r < 0.8; Supplementary Table S2), followed by confirmatory test Principal component analysis (PCA) analysis implemented in R-Programming (Supplementary Table S3). The predictive performances of the 10 simulated models were calibrated and evaluated using 25% of the data that uses AUC (Area under ROC curve) > 0.8, TSS (True Skill Statistics); Cohen's Kappa > 0.7 (Supplementary Table S4). These models were then projected onto plaeo and future (2070, RCP4.5) climatic scenarios based on different General Circulation Models (one LIG, two each for LGM and future) to determine the distribution range shifts and suitable habitats of M. himalaica. The ensemble consensus model was converted into binary presence (1)/absence (0) applying thresholds that allow a maximum of 50% probability of the suitable habitat (Forester et al., 2013). We also reclassified changes in LIG, LGM, and future conditions compared to current suitability into retracted, stable, and expanded areas. (see Supplementary Methods S1 for details).

The dispersal corridors were identified following the least cost path (LCP) methods using SDMtoolbox v2.0 (Brown et al., 2017) under the four climatic scenarios LIG, LGM, current, and future. For this analysis, we adopted the haplotype information of 29 localities, and populations that share haplotypes from two molecular markers (cpDNA and G3pdh) (Figures 1B, D). Firstly, we obtained a dispersal cost layer (resistance layer) by inverting the SDMs, and subsequently, we created a cost distance raster for each sampling locality using the resistance layer. Corridor layers were established based on the cost distance raster between two localities that shared haplotypes. To avoid oversimplifying landscape processes, we classified the value of each corridor layer into four intervals (three cutoff values: 1%, 2%, 5%) and reclassified as new values (5, 2, 1, 0, respectively). Finally, we summed up and standardized all of the pairwise reclassified corridor layers and identified dispersal maps of M. himalaica in an explicit landscape.

Results

Haplotype Diversity and Distributions

The total aligned sequences of the four combined cpDNA regions are 3290 bp (petL-psbE/838 bp, rps16-trnK/602 bp, rps16 intron/883 bp, trnS-trnG/967 bp) with 28 polymorphic sites and 71 indels (Supplementary Table S5). In total, 29 chlorotypes (C1–C29) were identified among 241 individuals (Table 1 and Figure 1B). The most common chlorotype was C1 (shared by 11 populations) followed by C12 (shared by six populations) (Table 1 and Figure 1B). Out of total chlorotypes, WHI and QTP groups (grouping by Samova) harboured two chlorotypes each (6.9%), WHN group contained four chlorotypes (13.8%), and the remaining 21 were occupied by HM group (72.4%); indicating a very high level of molecular diversity in HM region. Nineteen chlorotypes (65.5% of the total) were relatively rare, each of them restricted to a single population, whereas 15 of them were distributed in the HM region.

Out of 241 individuals, only about 5% (11 individuals) are heterozygous for partial G3pdh sequence under haplotypes H1, H4, H8, H11, H16, and H18. Population YY is the most diverse with six heterozygous individuals. The length of the aligned sequence of G3pdh was 941 bp, containing 18 nuclear haplotypes with 20 polymorphic sites (Supplementary Table S6). Thirteen haplotypes (72.2%) were unique, each endemic to a single population; 11 of them (84.6%) were restricted to the HM region. H1 was the most common haplotype within 16 populations of HM regions (Figure 1C). The chlorotype/haplotype frequencies and distributions in populations are listed and shown in Table 1 and Figure 1.

Genetic Polymorphism and Population Structure

Total haplotype diversity (H_d)/nucleotide diversity (π) for the cpDNA and G3pdh sequences were 0.901/0.0013 and 0.777/0.00173, respectively. For cpDNA, the highest value of H_d and π falls in population HS (H_d = 0.786; π = 0.0006); whereas for G3pdh, population JS exhibited particularly a high value (H_d = 0.607; π = 0.0007) (Table 1). Total genetic diversity (H_T) in the overall population of M. himalaica were 0.917/0.803 (cpDNA/G3pdh) (Table 3). The average within-population diversity for both cpDNA/G3pdh (H_S = 0.303/0.340) was also relatively high. Both, H_T and H_S are highest in the HM group (for WHI group not calculated, because of only one population). The inter-population differentiation (G_ST) of the two datasets i.e. cpDNA and G3pdh was 0.669 and 0.577, respectively. For both datasets, N_ST was significantly greater than G_ST (U > 1.96, p < 0.05), indicating the existing phylogeographical structure in M. himalaica (Table 3).

Samova of cpDNA reached a platform when K = 5 (F_CT = 0.872) and revealed a substantial spatial population genetic structure with five groups, i.e. HM (Hengduan Mountain) group with 18 populations, QTP (Qinghai-Tibetan Plateau) group with three populations, WHN (Western Himalaya Nepal: WC Nepal) I group with three populations (CJ, LJ, KM), WHN II group with four populations (MM, JM, LM, TM) and WHI (Western Himalaya India: Himachal Pradesh, India) group with one population (Supplementary Figure S1). Samova of G3pdh data reached a platform when K = 6 (F_CT= 0.6996, Supplementary Figure S1) and divided two same groups (QTP and WHI) like that of the cpDNA dataset. While the WHN was separated into two different subgroups (WHN1: CJ, LJ, MM, KM; WHN2: JM, LM, TM), and the HM was also divided into two subgroups (GS vs. other 17 populations). Considering the geographical unit of the HHM and the unstable population MM, we merged the divided subgroups and regarded four structured phylogroups (HM, QTP, WHN, and WHI; Figure 1) for M. himalaica. This grouping is also consistent with the Samova grouping (K=4) of cpDNA. Thus, based on four phylogroups, there is greater genetic partition among groups (78.58%/52.92%), than among populations within groups (17.30%/25.05%) or within populations (4.11%/22.03%) (Table 4). Examining the genealogical relationships/MJ network for chlorotypes and haplotypes, four clusters (HM, QTP, WHN, and WHI) were formed corresponding to the defined phylogroups (Figures 1A, C). The significant F-statistics from AMOVA suggested structured population within groups (Table 4). Additionally, the Mantel test revealed a statistically significant pattern of (IBD) for both cpDNA and G3pdh (r_M= 0.62/0.43; P = 0.001) (Table 4).

Phylogenetic Relationships and Divergence Time

The cpDNA topologies generated in MrBayes and Beast indicated M. himalaica was monophyletic with high posterior probabilities values which clustered the haplotypes into four lineages (Figure 2), corresponding to the four spatial population genetic groupings and genealogical relationships. Based on molecular dating, the onset of diversification for M. himalaica appeared in late Pliocene (crown age: 3.64 Ma; 95% HPD: 2.42–5.62 Ma; Figure 2) i.e. the first splitting event between the WHI lineage and the ancestor of the WHN, QTP, and HM lineages. Subsequently, the second splitting event between WHN lineage and the ancestor of QTP and HM lineage is 2.33 Ma (95% HPD: 1.37–3.87 Ma). Ultimately, the third splitting event between QTP and HM lineages occur during 1.66 Ma (95% HPD: 0.96–2.86 Ma).

Historical Demography

Tajima’s D and Fu’s Fs values were not significant, suggesting that M. himalaica conforms to the neutral hypothesis and these populations did not experience recent demographic expansion events (Table 5). The mismatch distribution for the overall chlorotypes/haplotypes were multimodal or bimodal (Supplementary Figure S2). While, the SSD and H_Rag statistics indicated a good statistical fit to the pure population growth and/or the spatial expansion model, excepting pure population growth in cpDNA data (p > 0.05; Table 5). This result suggested M. himalaica has experienced demographic events with structured populations that are shrinking in size or demographic equilibrium (Rogers and Harpending, 1992; Harpending, 1994).

BSP reconstructions based on combined markers (cpDNA and G3pdh) showed that the population of M. himalaica was quite stable after an ever-increasing phase during 0.6–0.075 Ma (Supplementary Figure S3), indicating no recent demographic expansion events. This case is similar for cpDNA marker when performed BSP reconstruction, while BSP of G3pdh showed a prolonged phase of demographic stability or no distinct population expansion (Supplementary Figure S3).

Ensemble Species Distribution Modelling and Visualization of Dispersal Corridors

The projected current distributions were generally the representations of the actual distributions and suitable habitat, which is consistent with present occurrence records (Figure 3A) except certain parts of the Eastern Himalaya. The paleodistribution reconstruction showed that during LIG, M. himalaica presented fragmented distribution patterns in the HM, and parts of the Western Himalaya (Figure 3B). Later during LGM, the species became more prominent towards the South of the HM (North Yunnan), WHI, and few in the southeastern QTP region (Figures 3C, D). The LGM distribution range of M. himalaica predicted using two models (CCSM4, MIROC-ESM) varies in some parts of the Western Himalaya and the HM (Figures 3C, D). In terms of the stability, we found that no matter which model was used, the distribution of M. himalaica during the LGM was farther south than the present one, suggesting northward expansion of populations after the LGM (Figures 4C, D). In addition, the predicted future distributions based on the MIROC-ESM model showed more migration to the north-west in the HM and west of the Himalaya compared to the present (Figures 3F and 4F). The result varies from the CCSM4 model, which was showing relatively more expansion towards the west in all regions (Figures 3E and 4E). It was predicted that species range expanded to southwards from LIG to LGM, but the range expansion is towards north-westwards after LGM up to current and the future (Figures 3 and 4).

Figure 3

Figure 4

Based on the least cost path analysis, putative dispersal corridors for the four periods (i.e. LIG, LGM, current, and future) were visualized using cpDNA and G3pdh markers (Figure 5). When comparing the dispersal routes across different time periods and genetic markers, dispersal generally followed isolated corridors between regional populations. There were no traces of connectivity between the populations of HM with the QTP and the Western Himalayas, and vice versa. This might be due to significant ridges and peaks of the local landscape in the Himalaya which formed a spatial barrier that intensified and fixed genetic isolation. The result exhibited continuous patterns of landscape connectivity among the populations present in the HM, indicating the existence of a dispersal corridor along mountains in the Hengduan regions. Besides the HM region, partial dispersal routes were also identified in the QTP and the WHN with average to low connectivity level (Figure 5). The cpDNA data did not show connectivity between two population groups from Jumla and Mustang–Manang of WHN; but G3pdh data revealed an additional area of dispersal in WHN beyond the dispersal corridor evidenced from the cpDNA data for the region (Figure 5). However, in this geologically dynamic region, the corridor shifts its route in different periods based on the distribution pattern, indicating that the dispersal corridor is not static for M. himalaica in the HHM. The distinct dispersal corridor in the middle of the HM during LIG through gorges of mountains to Yalongjiang and Daduhe rivers, shift to the south in the LGM showing strong connectivity through the middle Jinshajiang following Yalonjiang and paleo-route of Daduhe river. Similarly, in the current condition the high-value corridor passes through upper Jinshajiang, Yalongjiang, Daduhe rivers and rivulets or streams; whereas, in the future, the species is likely to lose its corridor in the south HM and shift to northwards following the spatial distribution range. In addition, a relatively larger area of dispersal was identified during the LGM period for both the makers (Figures 5B, C, H, I) compared to others. The strongest population genetic connectivity for M. himalaica occurred along with the river systems in the HM region.

Figure 5

Discussion

Evolutionary and Demographic History

It was speculated that the ancestor of Mirabilis himalaica might have migrated from North America to the Himalayas either by Beringia or through long-distance dispersal, evolving allopatrically into the extant species (Wang et al., 2019). Mirabilis himalaica diverged from the closest taxa during Early Pliocene c. 5.99 Ma (stem age; 95% HPD: 4.79–8.13 Ma), and the separation of the lineages spanning late Pliocene to early Pleistocene between 3.64–2.33 Ma (Figure 2). Late Miocene QTP uplift might have triggered speciation of M. himalaica in combination with the plentiful environmental alterations in the QTP and adjacent regions (Li et al., 1979). During the recent violent QTP uplifts (c. 3.6–1.5 Ma), intensification of the Asian monsoon (c. 3.6–2.6 Ma) (Li et al., 1979; An et al., 2001; Sun and Wang, 2005) were in progress, involving several mountain ranges in the HHM biodiversity hotspot (Harrison et al., 1992; Royden et al., 2008; Xie et al., 2019). Furthermore, tectonic events and climate change have set off the lineage divergence and rapidly occupying the newly available terrain during the Pleistocene. The uplift process might have coupled with the Pleistocene climatic oscillations leading to habitat diversification, restricted gene flow, and finally differentiation of the species (Xie et al., 2014; Luo et al., 2017; Shahzad et al., 2017; Ren et al., 2017a).

In addition, the glacial cycles of the Quaternary period in the HHM region have likely affected the demographic history of focal species. The results of the neutrality tests and MDA in conjunction with BSPs analysis suggested that M. himalaica did not experienced recent expansion events. Considering the poor ability to seed disperse, associated with the intense altitudinal gradient characterizing the Himalaya. We speculated that there is no large-scale distribution range expansion/contraction, while four lineages of M. himalaica, after diverging from each other survived in situ or experienced restricted regional migration (SDM and LCP) to shape the current phylogeographical pattern and the effective population size of all lineages remained relatively constant throughout the evolutionary period. From SDM analysis, during LGM the species was progressively confined to southwards in the HM, and traces of habitat expansion occurs in WHI, southeastern QTP region (Figures 4C, D) compared to predicted current and LIG distributions. The reason behind the southward movement of populations lies in the fact that the northernmost part of subtropical Central and Northern China was cooler at least 7–10°C and dryer by 200–300 mm/year (Sun and Chen, 1991; Zhou et al., 1991). At the time both steppe and desert vegetation expanded in a west to south directions. The localized cooling caused by glacial advance and later, subsequent warming due to glacial retreat might have exposed new habitats for colonization (Matthews, 1992) in the Northern and Western regions of the HHM. Therefore, our results suggested geological events along with climatic oscillations including monsoonal system and the Quaternary glaciation could have facilitated the formation and divergence of lineage and led to the accumulation of genetic diversity and differentiation of M. himalaica in the HHM region.

Pleistocene Unstable Habitat Connectivity

In this study, we identified the possible dispersal routes of M. himalaica since LIG to the future conditions in the HHM. Despite the adequacy of these mountains as a barrier to dispersal (Oheimb et al., 2013), there are several gorges and passes through the Main Divide, which might have permitted the dispersal of M. himalaica. Consequently, the major River system and the river valleys act as a dispersal corridor or recolonization route within the structured populations. Strikingly, the connectivity of dynamic habitats for M. himalaica was supported by the ensemble SDM prediction, which indicated that the population habitat connectivity also experienced potent change since LIG. Further, population genetic connectivity analysis based on the genetic data and SDM resistance layer similarly indicated that population genetic connectivity was prominent from LIG to current, and later in future conditions too (Figure 5). Population genetic connectivity and historical demographic changes of extant organisms are often associated with cyclical Pleistocene glaciations (Hewitt, 2000); and appears to be an important driver of inconsistency in the population genetic connectivity in M. himalaica. Besides, Pleistocene glaciations likely were restricted to high or middle elevations and did not affect deep slopes or valleys in the south-west mountainous region of China (Qu et al., 2014). As a result, the migration of populations to the deep slopes or valleys might have formulated strong dispersal corridors in the southern Hengduan Mountains during LGM compared to LIG. Our result showed that the corridors shift northward after LGM up to the future, indicating a purely consistent pattern with the spatial distribution. However, the rapid range change during the LIG to the future transition probably contributed to the patterns of genetic connectivity of local populations.

Response to Climate Change and Implications for Conservation

According to the forecasted future potential distribution, the focal species shift north-westward losing their potential suitability in hilly and lower mountainous regions by 2070 (Figures 4E, F). This extensive elevational shift might be due to global warming which combined with other biotic/abiotic factors fabricating unsuitable habitats at lower elevations. The annual mean warming rates during the period 1901–2020 was 0.19°C/decade (Ren et al., 2017b; Krishnan et al., 2019) over the Hindu Kush Himalaya, which is expected to increase further in future. In addition, the rate of warming for the Himalaya has been reported approximately three times (i.e. 1.5 °C/year) than the global average (i.e. 0.06 °C/year) (Shrestha et al., 2012). This increased warming rate at Himalaya region may prompt a significant issue for montane species. Likewise, the northward elevational shift of the species in the future does not ensure an increase in plant production itself. Eventually, high mountains serve as a geographical blockade and obstruct the species to relocate upwards due to the climatic condition atop the summit (Salick et al., 2009; Rana et al., 2017; Rana et al., 2019), because of summit trap phenomena (Pertoldi and Bach, 2007). This finally leads to the extinction of species due to no habitat for survival atop the mountain. Future distributions show the north-westward elevational shift of the species habitat and the dispersal corridor as well, which represents more credible conservation priority areas in the north-west higher elevational region of the HHM for M. himalaica. Geo-climatic variables are by all account not the only component to cause species habitat or population destruction because in present world anthropogenic factors are equally contributing to the cause. In this way, predicted expansion regions could be prioritized to protect the species from serious destruction of genetic diversity.

References

• 1

  AnZ. S.KutzbachJ. E.PrellW. L.PorterS. C. (2001). Evolution of Asian monsoons and phased uplift of the Himalaya-Tibetan Plateau since Late Miocene times. Nature411, 62–66. doi: 10.1038/35075035

  • CrossRef
  • Google Scholar

• 2

  AviseJ. C. (2000). Phylogeography: the history and formation of species (Cambridge, MA: Harvard University Press).

  • Google Scholar

• 3

  AviseJ. C. (2009). Phylogeography: retrospect and prospect. J. Biogeogr.36, 3–15. doi: 10.1111/j.1365-2699.2008.02032.x

  • CrossRef
  • Google Scholar

• 4

  BandeltH. J.ForsterP.RohlA. (1999). Median-joining networks for inferring intraspecific phylogenies. Mol. Biol. Evol.16, 37–48. doi: 10.1093/oxfordjournals.molbev.a026036

  • CrossRef
  • Google Scholar

• 5

  BrownJ. L.YoderA. D. (2015). Shifting ranges and conservation challenges for lemurs in the face of climate change. Ecol. Evol.5, 1131–1142. doi: 10.1002/ece3.1418

  • CrossRef
  • Google Scholar

• 6

  BrownJ. L.BennettJ. R.FrenchC. M. (2017). SDMtoolbox 2.0: the next generation Python-based GIS toolkit for landscape genetic, biogeographic and species distribution model analyses. PeerJ5, e4095. doi: 10.7717/peerj.4095

  • CrossRef
  • Google Scholar

• 7

  ChanL. M.BrownJ. L.YoderA. D. (2011). Integrating statistical genetic and geospatial methods brings new power to phylogeography. Mol. Phylogenet. Evol.59, 523–537. doi: 10.1016/j.ympev.2011.01.020

  • CrossRef
  • Google Scholar

• 8

  ChristeC.CaetanoS.AeschimannD.KropfM.DiademaK.NaciriY. (2014). The intraspecific genetic variability of silicious and calcareous Gentiana species is shaped by contrasting demographic and re-colonization process. Mol. Phylogenet. Evol.70, 323–336. doi: 10.1016/j.ympev.2013.09.022

  • CrossRef
  • Google Scholar

• 9

  ClarkA. G. (1990). Inference of haplotypes from PCR-amplified samples of diploid populations. Mol. Biol. Evol.7, 111–122. doi: 10.1093/oxfordjournals.molbev.a040591

  • CrossRef
  • Google Scholar

• 10

  DrummondA. J.RambautA. (2007). BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol. Biol.7, 214. doi: 10.1186/1471-2148-7-214

  • CrossRef
  • Google Scholar

• 11

  DrummondA. J.RambautA.ShapiroB.PybusO. G. (2005). Bayesian coalescent inference of past population dynamics from molecular sequences. Mol. Biol. Evol.22 (5), 1185–1192. doi: 10.1093/molbev/msi103

  • CrossRef
  • Google Scholar

• 12

  DrummondA. J.SuchardM. A.XieD.RambautA. (2012). Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol. Biol. Evol.29 (8), 1969–1973. doi: 10.1093/molbev/mss075

  • CrossRef
  • Google Scholar

• 13

  DupanloupI.SchneiderS.ExcoffierL. (2002). A simulated annealing approach to define the genetic structure of populations. Mol. Ecol.11, 2571– 2581. doi: 10.1046/j.1365-294X.2002.01650.x

  • CrossRef
  • Google Scholar

• 14

  Dupont-NivetG.LippertP. C.Van HinsbergenD. J.MeijersM. J.KappP. (2010). Palaeolatitude and age of the Indo–Asia collision: palaeomagnetic constraints. Geophys. J. Int.182, 1189–1198. doi: 10.1111/j.1365-246X.2010.04697.x

  • CrossRef
  • Google Scholar

• 15

  ExcoffierL.LischerH. E. (2010). Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol. Ecol. Resour.10, 564–567. doi: 10.1111/j.1755-0998.2010.02847.x

  • CrossRef
  • Google Scholar

• 16

  FavreA.PackertM.PaulsS. U.JahnigS. C.UhlD.MichalakI.et al. (2015). The role of the uplift of the Qinghai-Tibetan Plateau for the evolution of Tibetan biotas. Biol. Rev. Camb. Philos. Soc.90, 236–253. doi: 10.1111/brv.12107

  • CrossRef
  • Google Scholar

• 17

  FelsensteinJ. (1988). Phylogenies from molecular sequences: inference and reliability. Annu. Rev. Genet.22, 521–565. doi: 10.1146/annurev.ge.22.120188.002513

  • CrossRef
  • Google Scholar

• 18

  ForesterB. R.DeChaineE. G.BunnA. G. (2013). Integrating ensemble species distribution modelling and statistical phylogeography to inform projections of climate change impacts on species distributions. Divers. Distrib.19, 1480–1495. doi: 10.1111/ddi.12098

  • CrossRef
  • Google Scholar

• 19

  FoxJ.WeisbergS. (2011). An R Companion to Applied Regression, second ed (Thousand Oaks, CA: Sage).

  • Google Scholar

• 20

  FuY. X. (1997). Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics147 (2), 915–925.

  • Google Scholar

• 21

  GuindonS.GascuelO. (2003). A simple, fast and accurate algorithm to estimate large phylogenies by maximum likelihood. Syst. Biol.52 (5), 696–704. doi: 10.1080/10635150390235520

  • CrossRef
  • Google Scholar

• 22

  GuisanA.ZimmermannN. E. (2000). Predictive habitat distribution models in ecology. Ecol. Model.135, 147–186. doi: 10.1016/s0304-3800(00)00354-9

  • CrossRef
  • Google Scholar

• 23

  HarpendingH. (1994). Signature of ancient population growth in a low-resolution mitochondrial DNA mismatch distribution. Hum. Biol.66, 591–600.

  • Google Scholar

• 24

  HarrisonT. M.CopelandP.KiddW.YinA. (1992). Raising Tibet. Science255, 1663–1670. doi: 10.1126/science.255.5052.1663

  • CrossRef
  • Google Scholar

• 25

  HewittG. (2000). The genetic legacy of the Quaternary ice ages. Nature405, 907–913. doi: 10.1038/35016000

  • CrossRef
  • Google Scholar

• 26

  JianH. Y.TangK. X.SunH. (2015). Phylogeography of Rosa soulieana (Rosaceae) in the Hengduan Mountains: refugia and ‘melting' pots in the Quaternary climate oscillations. Plant Syst. Evol.301 (7), 1819–1830. doi: 10.1007/s00606-015-1195-0

  • CrossRef
  • Google Scholar

• 27

  KearseM.MoirR.WilsonA.Stones-HavasS.CheungM.SturrockS.et al. (2012). Geneious basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics28 (12), 1647–1649. doi: 10.1093/bioinformatics/bts199

  • CrossRef
  • Google Scholar

• 28

  KrishnanR.ShresthaA. B.RenG.RajbhandariR.SaeedS.SanjayJ.et al. (2019). “Unravelling climate change in the Hindu kush himalaya: rapid warming in the mountains and increasing extremes,” in The Hindu Kush Himalaya Assessment. Eds. WesterP.MishraA.MukherjiA.ShresthaA. (Cham: Springer). doi: 10.1007/978-3-319-92288-1_3

  • CrossRef
  • Google Scholar

• 29

  LarkinJ. L.MaehrD. S.HoctorT. S.OrlandoM. A.WhitneyK. (2004). Landscape linkages and conservation planning for the black bear in west-central Florida. Anim. Conserv.7, 23–34. doi: 10.1017/S1367943003001100

  • CrossRef
  • Google Scholar

• 30

  LiJ. J.FangX. M. (1999). Uplift of the Tibetan Plateau and environmental changes. Sci. Bull.44 (23), 2117–2124. doi: 10.1007/BF03182692

  • CrossRef
  • Google Scholar

• 31

  LiJ. J.WenS. X.ZhangQ. S.WangF. B.ZhengB. X.LiB. Y. (1979). A discussion on the period, amplitude and type of the uplift of the Qinghai‐Xizang (Tibetan) Plateau. Sci. Sin.22, 1314–1328.

  • Google Scholar

• 32

  LibradoP.RozasJ. (2009). DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics25, 1451–1452. doi: 10.1093/bioinformatics/btp187

  • CrossRef
  • Google Scholar

• 33

  LiuJ. Q.WangY. J.WangA. L.HideakiO.AbbottR. J. (2006). Radiation and diversification within the Ligularia-Cremanthodium-Parasenecio complex (Asteraceae) triggered by uplift of the Qinghai-Tibetan Plateau. Mol. Phylogenet. Evol.38, 31–49. doi: 10.1016/j.ympev.2005.09.010

  • CrossRef
  • Google Scholar

• 34

  LuD. Q.GilbertM. G.WuZ. Y.RavenP. H. (2003). “Nyctaginaceae,” in Flora of China, vol. 5. (Beijing & St. Louis: Science Press & Missouri Botanical Garden Press), 432–433.

  • Google Scholar

• 35

  LuoD.YueJ. P.SunW. G.XuB.LiZ. M.ComesH. P.et al. (2016). Evolutionary history of the subnival flora of the Himalaya-Hengduan Mountains: first insights from comparative phylogeography of four perennial herbs. J. Biogeogr.43, 31–43. doi: 10.1111/jbi.12610

  • CrossRef
  • Google Scholar

• 36

  LuoD.XuB.LiZ. M.SunH. (2017). The ‘Ward Line-Mekong-Salween Divide' is an important floristic boundary between the eastern Himalaya and Hengduan Mountains: evidence from the phylogeographical structure of subnival herbs Marmoritis complanatum (Lamiaceae). Bot. J. Linn. Soc185, 482–496. doi: 10.1093/botlinnean/box067

  • CrossRef
  • Google Scholar

• 37

  MantelN. (1967). The detection of disease clustering and a generalized regression approach. Cancer Res.27, 209–220.

  • Google Scholar

• 38

  MatthewsJ. A. (1992). The Ecology of Recently-Deglaciated Terrain (Cambridge University Press, Cambridge). doi: 10.1002/jqs.3390080213

  • CrossRef
  • Google Scholar

• 39

  MengH. H.SuT.GaoX. Y.LiJ. I.JiangX. L.SunH.et al. (2017). Warm-cold colonization: response of oaks to uplift of the Himalaya-Hengduan Mountains. Mol. Ecol.26, 3276–3294. doi: 10.1111/mec.14092

  • CrossRef
  • Google Scholar

• 40

  MulchA.ChamberlainC. P. (2006). Earth science-The rise and growth of Tibet. Nature439, 670–671. doi: 10.1038/439670a

  • CrossRef
  • Google Scholar

• 41

  NeiM.LiW. H. (1979). Mathematical model for studying genetic variation in terms of restriction endonucleases. Proc. Natl. Acad. Sci. U.S.A.76, 5269–5273. doi: 10.1073/pnas.76.10.5269

  • CrossRef
  • Google Scholar

• 42

  NeiM.TajimaF. (1981). DNA polymorphism detectable by Restriction Endonucleases. Genetics97 (1), 145–163.

  • Google Scholar

• 43

  NeiM. (1975). Molecular population genetics and evolution (North-Holland Pub. Co., Amsterdam, Netherlands).

  • Google Scholar

• 44

  OheimbP. V.AlbrechtC.RiedelF.BössneckU.ZhangH.WilkeT. (2013). Testing the role of the Himalayan Mountains as a dispersal barrier in freshwater gastropods (Gyraulus spp.). Biol. J. Linn. Soc.109, 526–534. doi: 10.1111/bij.12068

  • CrossRef
  • Google Scholar

• 45

  Otto-BliesnerB. L.MarshallS. J.OverpeckJ. T.MillerG. H.HuA.CAPE Last Interglacial Project Members (2006). Simulating Arctic climate warmth and icefield retreat in the Last Interglaciation. Science311 (5768), 1751–1753. doi: 10.1126/science.1120808

  • CrossRef
  • Google Scholar

• 46

  PeakallR.SmouseP. E. (2012). GenAlEx 6.5: Genetic analysis in Excel. Population genetic software for teaching and research-an update. Bioinformatics28, 2537–2539. doi: 10.1093/bioinformatics/bts460

  • CrossRef
  • Google Scholar

• 47

  PertoldiC.BachL. A. (2007). Evolutionary aspects of climate-induced changes and the need for multidisciplinarity. J. Therm. Biol.32, 118–124. doi: 10.1016/j.jtherbio.2007.01.011

  • CrossRef
  • Google Scholar

• 48

  PonsO.PetitR. J. (1996). Measuring and testing genetic differentiation with ordered versus unordered alleles. Genetics144, 1237–1245.

  • Google Scholar

• 49

  PosadaD.BuckleyT. (2004). Model selection and model averaging in phylogenetics: advantages of Akaike information criterion and Bayesian approaches over likelihood ratio tests. Syst. Biol.53, 793–808. doi: 10.1093/bioinformatics/14.9.817

  • CrossRef
  • Google Scholar

• 50

  QiuY. X.FuC. X.ComesH. P. (2011). Plant molecular phylogeography in China and adjacent regions: tracing the genetic imprints of Quaternary climate and environmental change in the world's most diverse temperate flora. Mol. Phylogenet. Evol.59, 225–244. doi: 10.1016/j.ympev.2011.01.012

  • CrossRef
  • Google Scholar

• 51

  QuY. H.EricsonP. G. P.QuanQ.SongG.ZhangR. Y.GaoB.et al. (2014). Long-term isolation and stability explain high genetic diversity in the Eastern Himalaya. Mol. Ecol.23, 705–720. doi: 10.1111/mec.12619

  • CrossRef
  • Google Scholar

• 52

  R Development Core Team. (2016). R v3.4.1: A language and environment for statistical computing (Vienna, Austria: R Foundation for Statistical Computing).

  • Google Scholar

• 53

  RambautA.DrummondA. J.XieD.BaeleG.SuchardM. A. (2018). Posterior summarisation in Bayesian phylogenetics using Tracer 1.7. Syst. Biol. 67 (5), 901–904. doi: 10.1093/sysbio/syy032

  • CrossRef
  • Google Scholar

• 54

  RambautA. (2012). FigTree, tree figure drawing tool v.1.4.0 (University of Edinburgh: Institute of Evolutionary Biology). [WWW document] URL http://tree.bio.ed.ac.uk/software/figtree/[accessed 20 June 2018].

  • Google Scholar

• 55

  RanaS. K.RanaH. K.GhimireS. K.ShresthaK. K.RanjitkarS. (2017). Predicting the impact of climate change on the distribution of two threatened Himalayan medicinal plants of Liliaceae in Nepal. J. Mt. Sci.14, 558–570. doi: 10.1007/s11629-015-3822-1

  • CrossRef
  • Google Scholar

• 56

  RanaS. K.LuoD.RanaH. K.O'NeillA. R.SunH. (2019). Geoclimatic factors influence the population genetic connectivity of Incarvillea arguta (Bignoniaceae) in the Himalaya–Hengduan Mountains biodiversity hotspot. J. Syst. Evol. (Published online) doi: 10.1111/jse.12521

  • CrossRef
  • Google Scholar

• 57

  RenG. P.MateoR. G.LiuJ. Q.SuchanT.AlvarezN.GuisanA.et al. (2017a). Genetic consequences of Quaternary climatic oscillations in the Himalayas: Primula tibetica as a case study based on restriction site-associated DNA sequencing. New Phytol.213, 1500–1512. doi: 10.1111/nph.14221

  • CrossRef
  • Google Scholar

• 58

  RenY. Y.RenG. Y.SunX. B.ShresthaA. B.YouQ. L.ZhanY. J.et al. (2017b). Observed changes in surface air temperature and precipitation in the Hindu Kush Himalayan region over the last 100-plus years. Adv. Clim. Change. Res.8 (3), 148–156. doi: 10.1016/j.accre.2017.08.001

  • CrossRef
  • Google Scholar

• 59

  RogersA. R.HarpendingH. (1992). Population-growth makes waves in the distribution of pairwise genetic-differences. Mol. Biol. Evol.9, 552–569. doi: 10.1093/oxfordjournals.molbev.a040727

  • CrossRef
  • Google Scholar

• 60

  RonquistF.TeslenkoM.MarkP. V. D.AyresD.DarlingA.HöhnaS.et al. (2012). MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. Syst. Biol.61 (3), 539–542. doi: 10.1093/sysbio/sys029

  • CrossRef
  • Google Scholar

• 61

  RoydenL. H.BurchfielB. C.van der HilsT. R. D. (2008). The geological evolution of the Tibetan plateau. Science321, 1054–1058. doi: 10.1126/science.1155371

  • CrossRef
  • Google Scholar

• 62

  SalickJ.ZhendongF.BygA. (2009). Eastern Himalayan alpine plant ecology, Tibetan ethnobotany, and climate change. Global Environ. Change19, 147–155. doi: 10.1016/j.gloenvcha.2009.01.008

  • CrossRef
  • Google Scholar

• 63

  ShahzadK.JiaY.ChenF. L.ZebU.LiZ. H. (2017). Effects of mountain uplift and climatic oscillations on phylogeography and species divergence in four endangered Notopterygium Herbs. Front. Plant Sci.8, 19–29. doi: 10.3389/fpls.2017.01929

  • CrossRef
  • Google Scholar

• 64

  ShresthaU. B.GautamS.BawaK. S. (2012). Widespread climate change in the Himalayas and associated changes in local ecosystems. PloS One7, e36741. doi: 10.1371/journal.pone.0036741

  • CrossRef
  • Google Scholar

• 65

  SpellenbergR.Flora of North America Editorial Committee (2003). “Nyctaginaceae,” in Flora of North America north of Mexico (New York: Oxford University Press), 14–74.

  • Google Scholar

• 66

  SunX. J.WangP. X. (2005). How old is the Asian monsoon system? Palaeobotanical records from China. Palaeogeogr. Palaeoclimatol. Palaeoecol.222, 181–222. doi: 10.1016/j.palaeo.2005.03.005

  • CrossRef
  • Google Scholar

• 67

  SunX. J.ChenY. S. (1991). Palynological records of the last 11,000 years in China. Quat. Sci. Rev.10, 537–544. doi: 10.1016/0277-3791(91)90047-X

  • CrossRef
  • Google Scholar

• 68

  SwoffordD. L. (2002). Paup*: Phylogenetic analysis using parsimony (and other methods) v. 4.0b10 (Sunderland, Massachusetts, USA: Sinauer Associates).

  • Google Scholar

• 69

  TajimaF. (1989). Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics123, 585–595.

  • Google Scholar

• 70

  ThuillerW.GeorgesD.EnglerR. (2014). Biomod2: ensemble platform for species distribution modelling. R package version 3. 1–64, [WWW document] URL http://CRAN.R-project.org/package = biomod2. [accessed 20 August 2018].

  • Google Scholar

• 71

  WangY. L.YanG. Q. (2014). Molecular Phylogeography and Population Genetic Structure of O. longilobus and O. taihangensis (Opisthopappus) on the Taihang Mountains. PloS One9 (8), e104773. doi: 10.1371/journal.pone.0104773

  • CrossRef
  • Google Scholar

• 72

  WangI. J.SavageW. K.BradleyH. S. (2009). Landscape genetics and least-cost path analysis reveal unexpected dispersal routes in the California tiger salamander (Ambystoma californiense). Mol. Ecol.18, 1365–1374. doi: 10.1111/j.1365-294X.2009.04122.x

  • CrossRef
  • Google Scholar

• 73

  WangS. L.LiL.CiX. Q.ConranJ. G.LiJ. (2019). Taxonomic status and distribution of Mirabilis himalaica (Nyctaginaceae). J. Syst. Evol.57 (5), 431–439. doi: 10.1111/jse.12466

  • CrossRef
  • Google Scholar

• 74

  XieH.AshJ. E.LindeC. C.CunninghamS.NicotraA. (2014). Himalayan-Tibetan plateau uplift drives divergence of polyploid poppies: Meconopsis Viguier (Papaveraceae). PloS One9 (6), e99177. doi: 10.1371/journal.pone.0099177

  • CrossRef
  • Google Scholar

• 75

  XieC.XieD. F.ZhongY.GuoX. L.LiuQ.ZhouS. D.et al. (2019). The effect of Hengduan Mountains Region (HMR) uplift to environmental changes in the HMR and its eastern adjacent area: tracing the evolutionary history of Allium section Sikkimensia (Amaryllidaceae). Mol. Phylogenet. Evol.130, 380–396. doi: 10.1016/j.ympev.2018.09.011

  • CrossRef
  • Google Scholar

• 76

  YeZ.YuanJ. J.LiM.DamgaardJ.ChenP. P.ZhengC. G.et al. (2018). Geological effects influence population genetic connectivity more than Pleistocene glaciations in the water strider Metrocoris sichuanensis (Insecta: Hemiptera: Gerridae). J. Biogeogr.45, 690–701. doi: 10.1111/jbi.13148

  • CrossRef
  • Google Scholar

• 77

  YinA.HarrisonT. M. (2000). Geologic evolution of the Himalayan-Tibetan orogen. Annu. Rev. Earth Planet. Sci.28, 211–280. doi: 10.1146/annurev.earth.28.1.211

  • CrossRef
  • Google Scholar

• 78

  YuH. B.ZhangY. L.LiuL. S.QiW.LiS. C.HuZ. G. (2015). Combining the least cost path method with population genetic data and species distribution models to identify landscape connectivity during the late Quaternary in Himalayan hemlock. Ecol. Evol.5, 5781–5791. doi: 10.1002/ece3.1840

  • CrossRef
  • Google Scholar

• 79

  YuH. B.ZhangY. L.WangZ. F.LiuL. S.ChenZ.QiW. (2017). Diverse range dynamics and dispersal routes of plants on the Tibetan Plateau during the late Quaternary. PloS One12 (5), e0177101. doi: 10.1371/journal.pone.0177101

  • CrossRef
  • Google Scholar

• 80

  ZhouY.QiuG.GuoD. (1991). “Changes of permafrost in China during Quaternary,” in Quaternary Geology and Environment in China. Eds. LiuT. S. (Beijing, China: Science Press), 86–94.

  • Google Scholar