Deep Intraspecific Divergence in the Endemic Herb Lancea tibetica (Mazaceae) Distributed Over the Qinghai-Tibetan Plateau

Qinghai-Tibetan Plateau (QTP) is an important biodiversity hub, which is very sensitive to climate change. Here in this study, we investigated genetic diversity and past population dynamics of Lancea tibetica (Mazaceae), an endemic herb to QTP and adjacent highlands. We sequenced chloroplast and nuclear ribosomal DNA fragments for 429 individuals, collected from 29 localities, covering their major distribution range at the QTP. A total of 19 chloroplast haplotypes and 13 nuclear genotypes in two well-differentiated lineages, corresponding to populations into two groups isolated by Tanggula and Bayangela Mountains. Meanwhile, significant phylogeographical structure was detected among sampling range of L. tibetica, and 61.50% of genetic variations was partitioned between groups. Gene flow across the whole region appears to be restricted by high mountains, suggesting a significant role of geography in the genetic differences between the two groups. Divergence time between the two lineages dated to 8.63 million years ago, which corresponded to the uplifting of QTP during the late Miocene and Pliocene. Ecological differences were found between both the lineages represent species-specific characteristics, sufficient to keep the lineages separated to a high degree. The simulated distribution from the last interglacial period to the current period showed that the distribution of L. tibetica experienced shrinkage and expansion. Climate changes during the Pleistocene glacial-interglacial cycles had a dramatic effect on L. tibetica distribution ranges. Multiple refugia of L. tibetica might have remained during the species history, to south of the Tanggula and north of Bayangela Mountains, both appeared as topological barrier and contributed to restricting gene flow between the two lineages. Together, geographic isolation and climatic factors have played a fundamental role in promoting diversification and evolution of L. tibetica.


INTRODUCTION
The Qinghai-Tibetan Plateau (QTP) is one of the largest and youngest plateaus in the world, formed by several uplift events after the collision of the Indian tectonic plate with the Asian plate, about 40 Ma (million years ago) (Guo et al., 2002;Spicer et al., 2003). Further significant uplift of the QTP occurred during the periods of the East Asian summer and winter monsoons about 15 Ma (Wan et al., 2007), 10-8 Ma (Molnar et al., 1993;Harrison et al., 1995) and 3.6-2.6 Ma (Zhisheng et al., 2001). The monsoon system interacted with the glacialinterglacial cycle and produced a more variable monsoon climate during the Pleistocene (Zhisheng et al., 2001). In recent decades, considerable disagreement has arisen on the relation between the uplifting of the plateau and the East Asian monsoons. Some scholars suggest that the uplift of the QTP modified the global and East Asian climate dramatically (Bloemendal, 1989;Zhisheng et al., 2001) and triggered and intensified the Asian monsoon, which in turn strongly influenced biological processes in the region (Li and Fang, 1999). In contrast, other scholars, such as Renner (2016) and Spicer (2017), suggest that there was no obvious impact on the East Asian monsoon from the uplifting of the QTP, even they hold the opinion that uplift having reached average heights of 4-5 km since the mid-Eocene. However, these climatic oscillations did affect the demography of some species, leading to their range shifting or extinction. Furthermore, the harsh climate of this region may have improved the adaptability of some local organisms (Davis and Shaw, 2001;Hewitt, 2004;Wan et al., 2016).
Numerous endemic species occur in the QTP and adjacent highlands, which represent centers of the preservation of ancient species and the differentiation of young species (Wu, 1988;Myers et al., 2000;Liu et al., 2012). A popular but rarely proved hypothesis is that the uplift of mountains creates environmental conditions (such as dispersed barriers or new habitats) that increase the rate of speciation (Xing and Ree, 2017). However, several phylogeographic studies have shown that certain species may have retreated during the ice age to refugia located at edge of the QTP, and recolonized the QTP and the adjacent highlands after the ice age, e.g., Juniperus przewalskii (Zhang et al., 2005), Metagentiana striata (Chen et al., 2008), and Pedicularis longiflora (Yang et al., 2008). Recent studies on Aconitum gymnandrum (Wang et al., 2009), Hippophae rhamnoides (Jia et al., 2012), and Spiraea alpina  suggest that some species also survived in the QTP at high altitude during the glacial period. For those species, multiple refugia may have remained during the glacial period, some on the QTP and others on its edge . For every species that has been researched, there is a species-specific feature in their evolutionary histories, even in some closely related species, such as S. alpina and S. mongolica . Therefore, further phylogeographic studies of a wider range of species are necessary to improve and refine the model for differentiation and formation of species in the region.
According to the present taxonomical treatment, Lancea Hook. f. et Thoms. is a small genus of the Mazaceae with only two species, L. tibetica and L. hirsuta (Hong et al., 1998). As a traditional Tibetan medicinal plant, L. tibetica has been used in the treatment of leukemia, intestinal angina, heart disease, and cough (Song et al., 2011b). Phytochemical studies on L. tibetica suggest that it contains more than 71 compounds that have pharmacodynamic effects, including anti-tumor, antioxidant, and hypoglycemic-inhibiting activities (Song et al., 2011a;Liu et al., 2014Liu et al., , 2015. L. tibetica is a perennial species endemic to the QTP, widely distributed in alpine meadows at altitudes of 2,000-4,500 m (Hong et al., 1998). The generation time for L. tibetica is 2 years according to our preliminary field observations. Under the inferior living condition, local human harvest the wild populations puts extra pressure on L. tibetica threatened with extinction (Tian et al., 2016). In this study, by combining ecological niche modeling and molecular data, we investigated the historical, genealogical and promoted diversity of L. tibetica, to gain insights into its intraspecific divergence and spatiotemporal population dynamics. Our study provides an important advance in knowledge of the population dynamics of endemic species on the QTP.

Population Sampling and Experimental Protocols
According to the Flora of China (Hong et al., 1998) and herbarium records from the Chinese Virtual Herbarium (CVH 1 ), L. tibetica mainly occurs in Gansu, Qinghai, Sichuan, Xizang, and Yunnan in China. It should be noted that few CVH herbarium records of L. tibetica are mainly before the 1960s. Some locations with just one or two records are difficult to sample again in our field investigation. In this study, a total of 429 individuals were collected from 29 populations covering the major distribution range of L. tibetica (Table 1 and Figure 1). Fresh leaves were sampled, dried in silica gel and kept at −20 • C until DNA extraction. Pedicularis rhinanthoides and P. chinensis were also sampled as outgroups (sample information may be found under the GenBank accession numbers MH628332-MH628339). All the sampling locations were geo-referenced, and voucher specimens deposited into the Herbarium of Northwest Plateau Institute of Biology (HNWP), Chinese Academy of Sciences.

Genetic Variation and Population Genetic Structure
During all analyses, insertion-deletion polymorphisms (indels) were coded as presence/absence characters. After alignment, cpDNA haplotypes and ITS genotypes were identified and distinguished using DnaSP v5.0 (Librado and Rozas, 2009). The level of genetic variation, total haplotype diversity (Hd) and nucleotide diversity (Pi) were also calculated in DnaSP. The program PERMUT (Pons and Petit, 1996) was used to estimate within-population diversity (H S ), total gene diversity (H T ), genetic differentiation (G ST ) and population subdivision of phylogenetically ordered alleles (N ST ) (Nei, 1987;Grivet and Petit, 2002). The G ST value represents the degree of genetic differentiation among the population and was calculated as G ST = (H T − H S )/H T (Raymond and Rousset, 1995). The U-statistical method was used to compare G ST and N ST (using 1,000 repeat replacement tests) and to check the geographical distribution pattern.
Population subdivision analysis was performed in the program SAMOVA v1.0 (Dupanloup et al., 2010), to define groups that are geographically homogeneous and genetically differentiated. The analysis used the data from cpDNA, run for K = 2-10, starting from 1,000 random initial conditions for each run, to obtain the maximal value of F CT for the most appropriate K-value and grouping method. Genetic differentiation based on cpDNA was estimated through analysis of molecular variance (AMOVA) as implemented in the program ARLEQUIN v3.5 (Excoffier and Lischer, 2010). To calculate the average effective gene flow, we used the formula N m = ([1/F ST ]−1)/2. Genetic distances    were estimated in ARLEQUIN with 1,000 permutation tests.

Phylogeny and Demographic History Based on cpDNA Sequences
Relationships among cpDNA haplotypes were constructed via a maximum-parsimony median-joining network using NETWORK 4.6 with default parameters (Bandelt et al., 1999). Phylogeny of cpDNA haplotypes was estimated using MrBayes 3.1.2 (Drummond et al., 2012). P. rhinanthoides and P. chinensis were selected as outgroups, as Pedicularis and Lancea were formerly in the same family Scrophulariaceae (Hong et al., 1998). The best-fitting GTR + G + I model was selected by MrModeltest 2.3 (Nylander, 2004) using the Akaike Information Frontiers in Genetics | www.frontiersin.org FIGURE 2 | Geographic distribution of haplotypes detected from the combined cpDNA sequences of L. tibetica (population codes as detailed in Table 1).

Criterion. For MrBayes, two independent Markov-chain Monte
Carlo analyses for 100,000,000 generations were performed with a random starting tree. One cold and three heated chains were run simultaneously, with trees sampled every 1,000 generations, and discarding the first 25% as burn-in. FIGTREE 1.3.1 (Rambaut, 2009) was used to display the tree. Tajima's D and Fu's F S statistics were calculated to test for evidence of range expansion (Tajima, 1989;Fu, 1997;Jaeger et al., 2005). A significant value for D or a significantly large negative value for F S may be the result of population expansion (ArisBrosou and Excoffier, 1996). To analyze the dynamic size of the populations, we performed mismatch distribution as implemented in ARLEQUIN. The observed and expected mismatch distribution of the sum of squared deviation (SSD) and Harpending's raggedness index (HRI) were used as test statistics. A unimodal shape of the mismatch distribution provides evidence of sudden population expansion during the history of a species. All the tests were implemented in ARLEQUIN v3.01 (Rogers and Harpending, 1992) with 1,000 significant permutations. When the sudden expansion model was accepted, the formula τ = 2ut was used to estimate the age of expansion (t), where τ is the total number of mutations and u is the mutation rate per generation for the whole analyzed sequence. The value of u is calculated as u = µkg, where µ is the substitution rate per nucleotide site in 1 year, k is mean sequence length of the analyzed DNA region and g is the generation time of the plant. We used the substitution rates (2 × 10 −9 s s −1 year −1 ) of cpDNA (Wolfe et al., 1987) to estimate the expansion time of both clades.

Divergence Time Analysis Based on ITS Sequences
There are few reports of fossil data of Lamiales, and only a few fossil records [Fraxinus L. (Call and Dilcher, 1992;Magallon, 2000); Catalpa Bur. (Meyer and Manchester, 1997)] are reliable (Manchester, 1999). Based on these fossil records, Nie et al. (2006) analyzed the divergence time of Lamiales, showed the divergence between Mazus reptans and L. tibetica was around 25 Ma. The sequence of M. reptans (LC027734) Frontiers in Genetics | www.frontiersin.org FIGURE 3 | Geographic distribution of genotypes detected from the ITS sequences of L. tibetica (population codes as detailed in Table 1).
was used as the outgroup in analyzing ITS data (Nie et al., 2006). The GTR + I base substitution model was selected with a loose molecular clock model of the uncorrelated index in BEAST 1.5.0 (Drummond and Rambaut, 2007). The two independent models were analyzed and combined using LogCombiner v1.5.3. Convergence was traced using TRACER v1.7 (Rambaut et al., 2018). The program TreeAnnotator v1.5.3 (Rambaut et al., 2018) was used to summarize the maximum credible tree. Finally, a tree showing ages for each branch was displayed in FigTree v1.3.1 (Rambaut, 2009).
To evaluate the effectiveness of these algorithms we used TSS and ROC values >0.7 to assemble the raster layers using median values. A total of 13 bioclimatic variables were chosen (annual mean temperature, mean diurnal range, isothermality, temperature seasonality, maximum temperature of warmest month, minimum temperature of coldest month, temperature annual range, mean temperature of wettest quarter, mean temperature of driest quarter, mean temperature of warmest quarter, annual precipitation, precipitation of wettest month, and precipitation of driest month) with low correlation and high informativeness after a jackknife procedure on the 19 bioclimatic variables downloaded from the WorldClim database (Robert et al., 2005). We selected the maximum entropy model and machine learning algorithm as implemented in MAXENT v3.3.3k (Phillips et al., 2006;Phillips and Dudík, 2008) to predict suitable climate models for both lineages. MAXENT can produce a useful model with a small sample size (Hernandez et al., 2006;Pearson et al., 2007;Wisz et al., 2008;Anderson and Gonzalez, 2011). We used all the 19 bioclimatic variables from 1950 to 2000, downloaded from the WorldClim database (Hijmans et al., 2005). In addition, we selected ten environmental variables (annual mean temperature, mean diurnal range, isothermality, maximum temperature of warmest month, minimum temperature of coldest month, mean temperature of driest quarter, mean temperature of warmest quarter, mean temperature of the coldest quarter, precipitation seasonality, precipitation of coldest quarter) to perform the tests. The restricted dataset was used to avoid including highly correlated variables and prevent potential overfitting (Peterson and Nakazawa, 2008). Model performance was evaluated by the area under the receiver operating characteristic curve (AUC) using the program MAXENT. We used a jackknife (or 'leave-one-out') procedure to train and test the model. Values between 0.7 and 0.9 indicated good discrimination (Swets, 1988).
To measure the niche similarity between lineages, we used ENMTools 1.3 (Warren et al., 2008(Warren et al., , 2010 to calculate Schoener's D and Warren's I indices (Warren et al., 2008) and quantify niche overlap: a value of 0 means ecological niches do not overlap at all, and 1 means the habitats are estimated to be equally suitable for both lineages. The overlap test was performed in layers using the program MAXENT. A niche identity test was obtained based on 200 pseudo-replicates to generate a distribution of the expected values of each index. The significance of observed and expected indices were estimated using SPSS v20.0 (IBM Corp, 2013).

Sequencing, Genetic Variation, and Population Genetic Structure
The total alignment length of four chloroplast gene regions (trnH-psbA, matK, trnL-F, and rbcL) in all individuals was 2,179 bp, included 14 substitutions and four indels (also coded as substitutions during analysis, Table 2). Based on those polymorphisms, we identified a total of 19 haplotypes (A-S), which were asymmetrically distributed across the 29 populations (Figure 2). The total estimated haplotype diversity (Hd) was 0.6247 and nucleotide diversity (Pi) was 0.00108 (Table 1). At the population level, the populations 9-15 showed a higher Hd and Pi. Haplotypes A and D were widely distributed in the south and north ranges, respectively (Figure 2). Populations  Estimates were acquired under a model of spatial expansion using ARLEQUIN. τ, time in several generations elapsed since the sudden expansion episode; HRag, Harpending's raggedness index; SSD, the sum of squared deviations; NC, not calculated; Ma, million years ago.  9, 10, and 12-15 showed higher haplotype and nucleotide diversities.
The total alignment length of ITS in all individuals was 693 bp, which included eight substitutions that enabled us to identify thirteen genotypes (G1-G13; Table 3).
In combination with the geographical distribution of L. tibetica, our results indicated that the G1-G4 genotypes were mainly distributed to the south of the Tanggula Mountains, while the other genotypes were found to the north (Figure 3).
The program SAMOVA divided all the populations into two groups based on the chloroplast sequences, corresponding to the south and north lineages. The south lineage included populations 1-7, 9, 10, 12, and 13 while the north lineage included populations 8, 11, and 14-29, although the F CT was not the highest. The F CT values changed very little with increasing number of groups (K) and were highest at K = 5 when populations 6, 14, and 15 formed three groups (these populations had a high proportion of private haplotypes). The average genetic diversity (H S ) was 0.311, while the total genetic diversity (H T ) FIGURE 7 | Distribution models for L. tibetica, simulated based on current, mid-Holocene (6 ka), last glacial maximum (20 ka), and last interglacial maximum (135 ka) periods.
was 0.644. The N ST (0.662) was significantly higher than G ST (0.507), as shown by a U-test (P < 0.01), indicating significant phylogeographical structure in L. tibetica. AMOVA revealed that 61.50% of genetic variation was partitioned among groups, 15.50% among populations within the group, and 22.94% within populations (Table 4). Moreover, the average gene flow among the populations and between the two groups of L. tibetica was 0.249 and 0.149, respectively.

Phylogeny and Demographic History Based on cpDNA Sequences
The Bayesian inference tree topology of the 19 cpDNA haplotypes strongly supported the hypothesis of two lineages (south and north; Figure 4A). Haplotypes in the south lineage occurred in populations from the south of the QTP, and haplotypes in the north lineage occurred in populations from the north. The maximum-parsimony median-joining network also grouped all the cpDNA haplotypes into two major groups (south and north) separated by two mutational steps ( Figure 4B).
The results of Tajima's D and Fu's FS was not significant. However, the observed mismatch distributions of haplotypes for each lineage failed to reject the spatial expansion model (SSD, H Rag values P > 0.05; Table 5 and Figure 5). The observed mismatch distributions of the whole population rejected the spatial expansion model (SSD, H Rag values P < 0.05; Table 5). Based on the range of the plastid DNA substitution rate, a haplotype sequence length of 2,179 bp and 2-year generation time, the expansion of the south lineage was estimated to have occurred at 0.727 Ma (with a confidence interval 0.024-1.412 Ma), and that of the north at 0.172 Ma (with a confidence interval 0.021-0.201 Ma).

Divergence Time Analysis Based on ITS Sequences
According to preliminary calculations using ITS sequence data, L. tibetica diverged from M. reptans around 25 Ma, and the divergence of L. tibetica between the major north and south lineages was dated at around 8.63 Ma (Figure 6). These estimates of dates of origin need to be treated with caution but the estimated divergence times correspond well with the geological evidence of the QTP uplift during the late Miocene and Pliocene (Li and Fang, 1999;Zheng et al., 2000;Mulch and Chamberlain, 2006).

Ecological Niche Modeling
The predicted distribution of L. tibetica underwent significant changes during the glacial-interglacial period (Figure 7). From the last interglacial maximum to the last glacial maximum to the mid-Holocene, the range of the predicted distribution of L. tibetica experienced successive reduction and expansion. There was no significant change from the mid-Holocene to the current period (Figure 7). The AUC values for ecological niche modeling of the north and south lineages were 0.986 and 0.960, respectively, indicating far better than a random prediction. A test of identity between the two lineages showed that there was distinct niche differentiation (P < 0.05). A background test of both lineages also showed that the ecological niches of the two lineages are well differentiated (Figures 8a,c,d). Values of Schoener's D and Warren's I indices suggested significant niche divergence between the south and north lineages (P < 0.05, Figure 8b).

Genetic Structure and Intraspecific Divergence
The average effective gene flow within the distribution range of L. tibetica is low, compared with that of previous studies on other species in the area, e.g., Spiraea mongolica (0.41) (Wang et al., 2014) and Camellia flavida (0.35) (Wei et al., 2017). Higher gene flow in those other species might have resulted in higher genetic differentiation among their populations. We found that the average effective gene flow among the two lineages of L. tibetica (0.149) was lower than that among the different populations (0.249). The seeds of L. tibetica are small and wingless and disperse near the parent plants, a feature that is likely to have enhanced the degree of genetic differentiation by restricting gene flow (Hong et al., 1998). However, we found high genetic differentiation among the populations, and most genetic variation was distributed among the populations and groups, based on SAMOVA. The geographic isolation of populations within species and variation in ecological factors are major driving forces to cryptic speciation (Hoskin et al., 2005;Liu et al., 2013).
The results of SAMOVA, the Bayesian inference tree and parsimony network analysis showed that L. tibetica comprises two major cpDNA groups. One group has its main geographic distribution to the north of the Tanggula and Bayangela Mountains, while the other group lies mainly to the south of the QTP. Gene flow across the whole region appears to be restricted by high mountains, suggesting a significant role of geography in the genetic differences between the two groups. Similarly, the ITS sequence variation showed clearly that the divergence of L. tibetica between the major north and south lineages was around 8.63 Ma. Although the estimates of dates of origin need to be treated with caution, they correspond well with geological evidence that the QTP experienced uplift during the late Miocene and Pliocene periods (Li and Fang, 1999;Zheng et al., 2000;Mulch and Chamberlain, 2006). This evidence suggests that the Tanggula and Bayangela Mountains appear to act as a geographical barrier for L. tibetica, probably imposed significant barriers on gene flow and divided the species into north and south lineages.
As indicated by the cpDNA, the ecological differences between the two lineages seem to represent species-specific characteristics that would be sufficient to keep the lineages separated to a high degree. Such distinct ecological niches would have reinforced the divergence of the two lineages following their initial spatial isolation. Thus, each of the two lineages may have given rise to some degree of differential adaptation to its respective environmental conditions. It is likely that the extensive QTP uplifts created fragmentation and isolation of habitats and niche differentiation, and provided the preconditions for the adaptive divergence of fragmented populations and subsequent speciation (Hewitt, 1996;Abbott and Brennan, 2014). In addition, about 9-8 Ma, enhanced aridity in the Asian interior and the onset of Indian and East Asian monsoons (Zhisheng et al., 2001) might have provided different ecological niches for the different lineages of L. tibetica. Some studies have reported that if related species live in significantly different niches, ecological divergence would likely be important in facilitating speciation, even in the presence of gene exchange (Nosil, 2008;Nosil et al., 2009a;Anacker and Strauss, 2014;Wan et al., 2016). Although ecological divergences in this case have not resulted in the emergence of new species, the initial divergence demonstrates the potential for ecological speciation. If geographic isolation and restricted gene flow are maintained long enough, they may eventually lead to reproductive isolation (Rieseberg and Burke, 2001;Nosil et al., 2009b;Thorpe et al., 2010), resulting in the formation of new species. Our results support the conclusion from previous studies that the uplift of the QTP and its associated climatic changes were most likely the main cause of plant diversification (Cun and Wang, 2010;Xu et al., 2010;Yang et al., 2012).

Quaternary Demographic History and Glacial Refugia of L. tibetica
Climate changes during the Pleistocene glacial-interglacial cycles had a dramatic effect on species distribution ranges (Comes and Kadereit, 1998;Hewitt, 2004), causing migration and/or extinction of populations, followed by periods of isolation, divergence and subsequent expansion (Taberlet et al., 1998;Cun and Wang, 2010). During the Pleistocene period, continuing climatic oscillations caused repeated shifts in the abundance of alpine species (Tang and Shen, 1995;Herzschuh et al., 2010). There were some glacial refugia on the QTP platform, giving some plant species a chance to survive the changing climate (Yang et al., 2008;Wang et al., 2009Wang et al., , 2015Li et al., 2011;Gao et al., 2016;Liu et al., 2018). In the present study, the simulated distribution from the last interglacial period to the current period showed that the distribution of L. tibetica experienced shrinkage and expansion (Figure 7). In the last glacial maximum period, extreme cold and dry weather substantially reduced its distribution from a continuous geographical distribution to a more fragmented pattern. Our distribution simulations for the mid-Holocene and the current period showed that the distribution ranges of L. tibetica were the same, and apparently more extensive than those in the last glacial maximum periods. Taken together, our results reveal that L. tibetica did experience population expansion.
Molecular data also provided further support for the above hypothesis. Based on cpDNA sequence variation, the north and south lineages of L. tibetica experienced a rapid range expansion at 0.172 Ma and 0.727 Ma, respectively (Table 5), consistent with the Pleistocene (Hewitt, 2000;Zheng et al., 2002). The largely open alpine regions that became available following the end of the major glaciation would have provided extensive opportunities for L. tibetica to expand its range. Indeed, such expansions of geographical range into alpine regions of the QTP have been reported previously for several plant species and are likely to have been common during the largest Pleistocene glaciation on the QTP (Liu et al., , 2013Gulzar et al., 2018;Lin et al., 2018). The QTP has been shown to be sensitive to climatic shifts, when plants were profoundly affected by alpine glaciation (Meng et al., 2007;Chen et al., 2008).
According to Taberlet and Cheddadi (2002), localities with high levels of genetic variation and unique haplotypes have often been recognized as possible refugia or as centers of diversification for species, whereas localities with low levels of genetic variation represent recent colonization. Some reports suggest that the mountainous areas of subtropical China may have provided refugia for warm-temperate evergreen species through periods of adverse climatic conditions (Bennett and Provan, 2008;Qiu et al., 2011;Wang et al., 2015). Our results showed that since the last glacial maximum, the south and north lineages experienced population expansion, while the population as a whole showed no expansion. In addition, even though L. tibetica has only two major lineages, some of its populations (9, 10, 12, 13, 14, and 15) contained higher haplotype and nucleotide diversities. It is likely that populations located in the Tanggula Mountains could survive in alternative habitats within a short distance, allowing biodiversity to persist during climate modifications (Hoorn et al., 2013). These patterns collectively suggest that areas south of the Tanggula and north of the Bayangela Mountains harbored refugia during the early Pleistocene, and then the southern and northern distribution ranges expanded rapidly at 0.727 and 0.172 Ma, respectively, during the interglacial periods. It might be a possible explanation for our finding that these populations have lower haplotype and nucleotide diversities than do some other species. Previous reports have suggested that rapid range expansion should decrease intra-population genetic diversity in the direction of spread (Hewitt, 1996;Soltis et al., 1997;Nason et al., 2002).

CONCLUSION
In conclusion, analyses of L. tibetica from our sample range, bringing together molecular phylogeography and species distribution modeling, indicate that a combination of geographic isolation and climatic factors have played a fundamental role in promoting diversification and evolution of this species. This study provides valuable evidence that advances research on genetic differentiation on the QTP.

AUTHOR CONTRIBUTIONS
FZ and SC conceived and designed the experiments. MX and ZT performed the experiments. MX, ZT, FZ, and GK analyzed the data. QG, RX, YZ, and JY contributed reagents, materials, and analysis tools. MZ wrote the paper. FZ, SC, and GK reviewed and edited the paper. All authors approved the final version of the manuscript.