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

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 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.

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).

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 (p; 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.fluxusengineering.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) 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. 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  Table 1 for population codes). In network figure, the size of circles corresponds to the frequency of each chloro/haplotypes and the light green squared box represents chloro/haplotypes missing from the dataset and each branch represents one mutation. Pie charts indicate the sample of each population and divisions within correspond to chloro/haplotypes with a number of individuals. The dotted line delineates the phylogroups (HM, Hengduan Mountains group; QTP, Qinghai-Tibetan Plateau group; WHN, Western Himalaya Nepal group; and WHI, Western Himalaya India group). Map source: http://www.diva-gis.org/ and http://www. worldclim.org. 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.

Historical Demographic Analyses
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 .

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 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 S3). The predictive performances of the 10 simulated models were calibrated and evaluated using 25% of the data that uses AUC (Area under ROC   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.

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 (p) 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 p falls in population HS (H d = 0.786; p = 0.0006); whereas for G3pdh, population JS exhibited particularly a high value (H d = 0.607; p = 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 interpopulation 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).

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)

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). 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.

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 . During the recent violent QTP uplifts (c. 3.6-1.5 Ma), intensification of the Asian monsoon (c. 3.6-2.6 Ma) An et al., 2001;Sun and Wang, 2005) were in progress, involving several mountain ranges in the HHM biodiversity hotspot (Harrison et al., 1992; TABLE 5 | Results showing mismatch distribution analysis under models of (A) pure population growth and (B) spatial expansion, and neutrality test for the cpDNA and G3pdh sequence of M. himalaica; tests were conducted in ARLEQUIN with the sum of squared deviations (SSD) and Harpending's (1994) 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 northwestward 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.

AUTHOR CONTRIBUTIONS
HR and HS planned and designed the research. HR and SR carried out sampling and performed experiments. HR and DL analyzed data. HR and DL conceived the manuscript with the support of SR and HS.