Out of Refugia: Population Genetic Structure and Evolutionary History of the Alpine Medicinal Plant Gentiana lawrencei var. farreri (Gentianaceae)

Understanding the genetic structure and evolutionary history of plants contributes to their conservation and utilization and helps to predict their response to environmental changes. The wildflower and traditional Chinese and Tibetan medicinal plant Gentiana lawrencei var. farreri is endemic to the Qinghai-Tibetan Plateau (QTP). To explore its genetic structure and evolutionary history, the genetic diversity, divergence, and demographics were analyzed in individuals from 31 locations across the QTP using 1 chloroplast marker and 10 nuclear microsatellite loci. High genetic diversity was detected in G. lawrencei var. farreri, and most of the genetic variance was found within populations. Values of FST in G. lawrencei var. farreri from nuclear microsatellite and chloroplast data were 0.1757 and 0.739, respectively. The data indicated the presence of isolation by distance. The southeast edge of the QTP was the main refugium for G. lawrencei var. farreri, and one microrefugium was also detected in the plateau platform of the QTP. Both nuclear microsatellite and chloroplast data indicated that the populations were divided into two geographically structured groups, a southeast group and a northwest group. The current genetic pattern was mainly formed through recolonization from the two independent refugia. Significant melt was detected at the adjacent area of the two geographically structured groups. Approximate Bayesian computation showed that the northwest group had diverged from the southeast group, which then underwent population expansion. Our results suggest that the two-refugia pattern had a significant impact on the genetic structure and evolutionary history of G. lawrencei var. farreri.


INTRODUCTION
As the highest and largest plateau in the world, the Qinghai-Tibetan Plateau (QTP) has experienced drastic environmental and climate changes during the past several million years (Zhang et al., 2000;Zheng et al., 2002). The fragile ecosystem of the plateau is undergoing drastic changes at present due to the effects of global climate change (e.g., Che et al., 2014;Ma et al., 2017). Exploring the ways in which organisms responded to environmental and climate changes in the past may help us understand and predict how organisms will respond to such changes in the future.
The QTP is one of the world's biodiversity hot spots, with a high density of endemic plant species (Myers et al., 2000). The QTP has thus recently received a great deal of attention from botanists. The evolutionary dynamics of a number of plant species in the QTP have been explored; these include Saxifraga (Ebersbach et al., 2017;Gao et al., 2017), Rhodiola (Li et al., 2018), and Spiraea (Khan et al., 2018). Additional examples have been cited in a number of recent reviews (Qiu et al., 2011;Liu et al., 2014;Wen et al., 2014;Favre et al., 2015;Sun et al., 2017;Xing and Ree, 2017;Yu et al., 2017). These studies indicated that evolutionary patterns varied among different taxa, and no common evolutionary models have emerged. However, it is undoubtedly true that geological and climatic changes in the QTP had triggered plant speciation and diversification in this area (Wen et al., 2014;Favre et al., 2016). Although more and more taxa have been studied in recent years, the QTP is still poorly studied in comparison to mountainous areas in Europe and the Americas (Favre et al., 2015;Sun et al., 2017). Therefore, more research is required in order to explore the dynamics of speciation and the evolutionary history of plants in the QTP.
As a typical alpine taxon, Gentiana L. occurs in numerous mountain systems of the world and encompasses 362 species (Ho and Liu, 2001). Hosting c. 250 species (Ho and Pringle, 1995), the mountain ranges surrounding the QTP are the main center of diversity for Gentiana and acted as the primary source area for dispersals to other parts of the world . The diversification within Gentiana has been influenced by geological and climatic changes in the QTP , similar to many other alpine taxa (Wen et al., 2014;Matuszak et al., 2016;Yu et al., 2017). The previous studies focused mainly on the phylogenetics and evolution in or above the section level (e.g., Yuan et al., 1996;Zhang et al., 2009;Favre et al., 2016;Sun et al., 2018b). Few studies have explored intraspecific population genetics. Population genetics was considered in just one species in the section Cruciata Gaudin (Lu et al., 2015), which is endemic to the QTP, and seven species in the section Ciminalis (Adanson) Dumorti (Diadema et al., 2005;Kropf et al., 2006;Christe et al., 2014), which is endemic to Europe. These population genetics studies examined how species in the genus diverged with the changing environments and offered insights into their mode of speciation. As a diversity center and source area, more taxa from the QTP should be studied to understand the speciation pattern and evolutionary dynamics of Gentiana.
In this study, we focused on Gentiana lawrencei var. farreri T. N. Ho (Gentianaceae), which belongs to section Kudoa (Masamune) Satake and Toyokuni ex Toyokuni and is endemic to the QTP. It lives in alpine meadows and wet meadows with altitude from 2,400 to 4,600 m. Its flower and fruit period is from August to October and pollinators are insects, mostly Bombus, Thripidae, and Formica (Hou et al., 2009), and seeds are not dispersed very far from the mother plant. The species is a perennial wildflower that has been domesticated for horticultural gardening (Rybczyński et al., 2015). It is also used in traditional Chinese and Tibetan medicine (Ho and Liu, 2001;Yang et al., 2012). In the QTP, G. lawrencei var. farreri is a very common species in its section Kudoa, in which most taxa have geographically limited distributions (Ho and Pringle, 1995). Exploring the intraspecies genetic diversification and evolutionary history of the common species can contribute to its utilization and will be helpful to understand the divergence and speciation of the section as well. Using both chloroplast DNA (cpDNA) and nuclear microsatellites, we investigated the population genetics structure and evolutionary history of G. lawrencei var. farreri from the point of matrilineal inheritance as well as nuclear inheritance.

Sampling
We sampled 31 populations of G. lawrencei var. farreri throughout the QTP, amounting to 423 individuals (Table 1 and Figure 1). For large populations (>100 individuals; N pop = 28), 10-30 mature plants were randomly sampled along a single transect. Sampled individuals were at least 20 m apart. For small populations (<100 individuals; N pop = 3), 25-50% of the plants (3-9 individuals) were sampled. Voucher specimens were deposited in the herbarium of the College of Life Science, Luoyang Normal University.

Molecular Data Acquisition
Total genomic DNA was extracted from dry leaves with a Dzup plant genomic DNA extraction kit (Sangon, Shanghai, China). Through comparison analysis of chloroplast genomes (Fu et al., 2016;Sun et al., 2018b), one divergence hot spot was identified in ycf 1 (Forward: TGCTTGTATTATTACTTTGTGG; Reverse: TTGGATTGGTATTAGTCTGGAT) and was verified in the sampled populations. The ycf 1 was then amplified in all of the individuals of G. lawrencei var. farreri. PCRs were performed in 25-µL volumes containing approximately 20 ng of template DNA, 1× PCR Buffer, 1.5 mM MgCl 2 , 0.25 mM of each dNTP, 0.20 mM of each primer, and 1 unit of Taq DNA polymerase (Takara, Dalian, China). The PCR cycling profile included an initial step of 5 min at 95 • C followed by 35 cycles of denaturation at 94 • C for 50 s, 45 s of annealing at 53 • C, 20 s at 72 • C, and a final extension at 72 • C for 6 min. PCR products were sequenced in the reverse direction on an ABI 3730 xl automated capillary sequencer (Applied Biosystems, Foster City, CA, United States) with BigDye v3.1 (Applied Biosystems).
To obtain the phylogenetic relationship among haplotypes, maximum-likelihood (ML) analyses were conducted in PhyML 3.0 (Guindon and Gascuel, 2003) using the subtree pruning and regrafting (SPR) method with the HKY model which was estimated in jModelTest 2.1.7 (Darriba et al., 2012). The sequences of outgroups were extracted from the complete chloroplast genomes of Gentiana straminea Maxim. (Ni et al., 2016), G. macrophy Pall.  and G. stipitata Edgeworth (Sun et al., 2018b). The robustness of the ML trees was tested with 1000 bootstrap replicates. A statistical parsimony network was constructed and visualized utilizing PopArt 1.7 (Leigh and Bryant, 2015) to elucidate matrilineal evolutionary relationships.
The divergence times of cpDNA haplotypes were estimated using the Bayesian method implemented in BEAST 1.7.5 (Drummond et al., 2012) under the HKY substitution model, the Yule model, and an uncorrelated lognormal clock model (Drummond et al., 2006). The outgroups used in the ML analysis were used here as well. We constrained one of the nodes with a seed fossil of G. straminea by lognormal priors with an offset at 5.0 Ma, a mean of 0.7, and a standard deviation of 1.0, as applied by Pirie et al. (2015) and Favre et al. (2016). We ran three independent MCMC chains with 50,000,000 generations, sampling every 5,000th generation and discarding the initial 20% as a burn-in. Convergence was confirmed in TRACER 1.5 1 and judged by ESS values (>200). Trees were summarized using TreeAnnotator 1.7.5 (Drummond et al., 2012).
To assess historical population demography, Extended Bayesian Skyline Plots (EBSP) were generated in BEAST. The Bayesian MCMC was run for 100,000,000 generations, discarding the initial 10% as a burn-in. Tajima's D (Tajima, 1989) and Fu's Fs (Fu, 1997) were also estimated using 10,000 coalescent simulations in ARLEQUIN.

Microsatellite Data
The presence of scoring errors was checked using MICROCHECKER 2.2.3 (van Oosterhout et al., 2004). The raw microsatellite data was uploaded to Figshare (doi: 10.6084/m9.figshare.7195511). The null allele frequency estimates were analyzed following the Expectation Maximization algorithm (Dempster et al., 1977) in FreeNa (Chapuis and Estoup, 2007). Departure from Hardy-Weinberg equilibrium (HWE) and tests of genotypic linkage disequilibrium (LD) were carried out using Genepop 4.0 (Rousset, 2008). We calculated gene diversity according to Nei (1987) and observed heterozygosity (H O ), expected heterozygosity (H E ), and pairwise F ST using ARLEQUIN 3.5 (Excoffier and Lischer, 2010). Allele richness (R S ) was calculated using FSTAT 2.9 (Goudet, 2001). A hierarchical analysis of molecular variance (AMOVA; Excoffier et al., 1992) was used to further quantify genetic differentiation in ARLEQUIN with 1000 permutations. We looked for significant correlations between genetic and geographical distances using ARLEQUIN 3.5. For checking isolation by distance, F ST was regressed against straight geographical distance and tested for significance using the Mantel test with 10,000 permutations. Recent changes of effective population sizes were tested with the Wilcoxon signed-rank test implemented in BOTTLENECK 1.2.02 (Piry et al., 1999) using the two-phase mutational model (TPM) with 10, 000 iterations.
We investigated the genetic structure among populations and individuals using the Bayesian clustering algorithm implemented in STRUCTURE 2.3.1 (Pritchard et al., 2000). We ran 10 independent replicates testing for 1-12 clusters (K). Each run started with a burn-in of 100,000 followed by 1 million iterations. The most likely true value of K was determined using the second rate of change in the likelihood distribution ( K; Evanno et al., 2005) in STRUCTURE HARVESTER 0.6.94 (Earl and vonHoldt, 2012). Graphical representation of individual cluster assignments was performed using DISTRUCT 1.1 (Rosenberg, 2004).

Microsatellite and Chloroplast Sequence Data
Demographic history was modeled through coalescent-based approximate Bayesian computation (ABC; for reviews see Bertorelle et al., 2010;Csillery et al., 2010). We divided populations of G. lawrencei var. farreri into two lineages identified by STRUCTURE. The first lineage (named SE) contained populations from the southeast; the second (named NW) contained populations from the northwest. In this study, we tested whether population divergence was followed by changes in effective population size. The ABC analyses were conducted using DIYABC 2.0 (Cornuet et al., 2014). First, three population divergence scenarios were tested. Scenario 1 assumed one ancestral population with an effective population size of Na that split into two populations of effective size N1 (SE) and N2 (NW) at time t1. These daughter populations were then simulated to evolve independently. Scenario 2 assumed that SE was the ancestral population and that NW diverged from SE at time t1. Scenario 3 was the same as scenario 2 but swapped SE and NW. Based on the scenario with the highest support from the population divergence data (scenario 2; see section "Results"), three scenarios with effective population size changes were then tested. Compared with scenario 2, scenario 4 assumed that NW had undergone population expansion at time t2 after diverging from SE at time t1. Scenario 5 assumed that a small number of founders (N1f) originated from SE at time t1 and that the population had undergone expansion at time t2. Scenario 6 assumed that both SE and NW had undergone population expansion. Scenarios 4, 5, 6, and 2 were compared together. In the ABC analysis, the prior distributions of all parameters were uniform, and the same range was used for common parameters between models. For all six scenarios, we assumed that microsatellites evolved under a stepwise-mutation model and that the ycf 1 evolved under a default model, and that all individual loci had the same value. A total of 500,000 coalescent simulations were run for each scenario. The posterior probabilities of the scenarios were subsequently estimated based on (i) a direct estimate, in which the 500 data sets with summary statistics closest to the target values were extracted, and (ii) a polychotomous logistic approach that recovers the 1% closest to the simulated data sets. For the scenario with the highest support, posterior distributions for parameters were estimated using a local linear regression on the closest 1% of the simulated data sets.

Sequence Characters and Genetic Diversity
All of the samples were amplified at the ycf1 region. The aligned sequence of the ycf1 region in this study was 408 bp in length. The G. lawrencei var. farreri dataset consisted of 17 base substitutions (Supplementary Table S1) that identified 16 haplotypes (H1-H16) among the 423 individuals. Within the total population of G. lawrencei var. farreri, 10 populations included a single haplotype, and 9 haplotypes were exclusive to single populations (Figure 1). The values of h and π ranged from 0 to 0.714 and 0 to 0.0083, respectively ( Table 1). The averages of h and π of G. lawrencei var. farreri were 0.263 and 0.0018, respectively. Populations from the southeast generally had higher levels of genetic diversity (h = 0.414 and π = 0.0026) than the overall average values.
All 423 of the samples were amplified at 10 nuclear microsatellites. Among all populations, the frequencies of null allele in the 10 microsatellite loci were low (ranged from 0.002 to 0.066) thus all loci were included in following analysis. Overall, population-level indices of genetic diversity, including Na, R S , H O , and H E , were variable across populations ( Table 1). The mean Na ranged from 3.4 to 10.4 with an average of 6.48. The mean R S ranged from 1.94 to 3.89 with an average of 3.17. The mean H O and H E ranged from 0.460 to 0.820 and 0.450 to 0.790 and averaged 0.692 and 0.674, respectively. Populations from the southeast generally had higher levels of genetic diversity (Na = 7.3, R S = 3.54, H O = 0.731, and H E = 0.720) compared to the overall average values. Among the 10 loci, the number of loci/location that deviated from HW at the population level ranged from 0 to 3 with an average of 1.10. Pairwise comparisons among the 10 microsatellite loci revealed no evidence of linkage disequilibrium (P > 0.05). The data from the 10 microsatellite loci used in G. lawrencei var. farreri are summarized in Supplementary Table S2.

Phylogenetic Analysis
The ML trees of the 16 cpDNA haplotypes (Figure 1 and Supplementary Figure S1) supported two clades; this was approximately consistent with the parsimony network (Figure 2). Clade I contained 8 haplotypes (H9-H16) that were mainly distributed in the southeast. Clade II included 7 haplotypes (H1-H8); within these, H2 was a distinct high-frequency haplotype located in the center of the network. H12 was also a high-frequency haplotype located in the center of Clade I. The molecular dating analysis of the cpDNA dataset (Figure 3) estimated that the two lineages in G. lawrencei var. farreri

Genetic Structure
In the 16 cpDNA haplotypes detected from G. lawrencei var. farreri, H2 was dominant and appeared in 27 (87.10%) populations. All of the other haplotypes appeared in less than 10 populations. The southeast of the distribution area (populations DF, KD, XC, and XGLL) contained 8 (50%) haplotypes and YuShu-ChenDuo region (populations QML, CD and YSb) contained 4 haplotypes. SAMOVA results indicated that F CT reached a plateau when K = 2 (Supplementary Table S3). Populations GZa, DF, XC, and XGLL from the south were clustered into one group, while the remaining populations formed a separate group. According to the clustering, AMOVA based on the cpDNA dataset revealed that most of the genetic variation occurred among groups (67.05%), with little further subdivision among populations within groups (6.88%) or within populations (26.07%) ( Table 2). Based on cpDNA data, F ST in G. lawrencei var. farreri was 0.739, while the value of G ST was 0.474. The G ST of two possible hot spots was also calculated. The values of G ST in the southeast populations (including DF, KD, XC, and XGLL) and YuShu-ChenDuo populations (including QML, CD, and YSb) were 0.374 and 0.111, respectively.
In the Bayesian individual-based clustering analysis based on the microsatellite dataset, the log-likelihood values for the number of clusters (K) increased with each K-value and did not plateau (Supplementary Figure S2). The method proposed by Evanno et al. (2005) suggested that the number of clusters was 2 ( K = 2, Supplementary Figure S2). Assuming two genetic clusters, G. lawrencei var. farreri individuals sampled from the southeast were predominantly assigned to Cluster 1, while individuals from the northwest were predominantly assigned to Cluster 2 (Figure 4). AMOVA revealed that 8.36% of the genetic variance was found between groups (Va = 0.330), 12.65% among populations within species (Vb = 0.500), and 78.99% within populations (Vc = 3.121) ( Table 2). The value F ST was 0.2101, and the pairwise F ST between the southeast edge group and the plateau platform group was 0.1533 (P = 0.000). The population pairwise values of F ST are listed in Supplementary  Table S4. "Isolation By Distance" showed a significant correlation between genetic and geographical distance among G. lawrencei var. farreri populations (r 2 = 0.485, P < 0.01).

Demographic History
Based on cpDNA data, the EBSP generated in BEAST showed that the historical population demography of G. lawrencei var. farreri, west and east groups were all constant within the most recent 200 Kya (Supplementary Figure S3). The values of Tajima's D were generally negative in G. lawrencei var. farreri populations but only significant in population MQ (Supplementary Table  S5). Values of Fu's F S were generally positive in G. lawrencei var. farreri populations, and all of the values were not significant (Supplementary Table S5). Based on microsatellite data, recent population declines were found within two (XC and MDa) of the sampled locations using the program BOTTLENECK (P < 0.05, Supplementary Table S6). We evaluated population divergence and effective population size changes from the combined nuclear and cp data. We first compared three population divergence scenarios; the ABC simulations showed slightly higher support for scenario 2 (subpopulation NW diverged from subpopulations SE) with direct estimation (not shown); there was stronger support with logistic regression, and thus scenario 2 was favored as the most probable ( Figure 5A). Effective population size changes were compared based on scenario 2. Direct estimation showed similar possibilities among different scenarios (not shown); however, logistic regression yielded strong support for scenario 5, suggesting that SE had undergone population expansion (Figure 5B). Parameter estimation based on scenario 5 showed that the mean value of N1, N2, N4, t1, and t3 were 99,300, 64,800, 82,300, 7,650, and 4,330, respectively. Although some of the estimated parameters showed low robustness, ABC analyses conclusively supported the scenario in which the northwest clade diverged from the southeast clade, which then experienced population expansion.

DISCUSSION
Our genetic investigation revealed that G. lawrencei var. farreri is divided into two geographically structured groups (Figures 1, 2). This genetic structure plausibly reflects the postglacial colonization history.

Genetic Divergence and Glacial Refugia
Due to maternal inheritance, plastid datasets can provide information on past changes in species distribution when colonization of new habitats occurs through seeds (Petit et al., 2003). Based on the cpDNA dataset, genetic diversity in the southeast edge of the QTP is significantly higher than in other regions. Since high diversity would be achieved through the mixing of already-existing genetic information but not created anew, the refugia should be identified not only by high genetic diversity but also by their genetic uniqueness (Petit et al., 2003). The southeast edge of the QTP contained half of the haplotypes, including several exclusive haplotypes, such as H10, H11, and H14-H16, which were clustered in one clade in the ML tree (Figure 1 and Supplementary Figure S1). The value of G ST in this area was also higher than in the remaining areas. The two central haplotypes in the network (H2 and H12), which were the ancient haplotypes, occurred in this area as well. The populations situated close to the southeast edge of the QTP (JD, DG, GZb, SD, AB, and REG) were also highly divergent. Meanwhile, the diversity declined away from the southeast edge of the QTP. These characteristics are well fitted by the refugia hypothesis (Hewitt, 2000;Petit et al., 2003). Therefore, the southeast edge of the QTP appears to have been the refugium of G. lawrencei var. farreri during glaciations. The nuclear microsatellites dataset also indicated that populations from the southeast edge of the QTP had higher genetic diversity than other areas. This area is well studied and is the most important refugium for alpine plants in the QTP (reviewed in Qiu et al., 2011;Liu et al., 2014;Sun et al., 2017). We also noted that the genetic diversity of the YuShu-ChenDuo region (populations QML, CD, and YSb) was higher than in other areas of the plateau platform. The haplotypes occurring in the plateau platform were observed from single clades in the ML tree except for H12, which occurred in both the YuShu-ChenDuo area and the southeast edge of the QTP. Although the value of G ST in this area was much lower than in the southeast populations, it was higher than in the adjacent areas of the plateau platform. Therefore, the YuShu-ChenDuo area may be a microrefugium of G. lawrencei var. farreri. Some microrefugia had already been detected on the QTP platform in hardy plants such as Sibiraea (Fu et al., 2016) and Rhodiola sect. Trifida (Li et al., 2018); these have been reviewed in Qiu et al. (2011) andLiu et al. (2014). The multi-refugia pattern will have obvious impacts on the genetic structure of plant populations.

Genetic Structure
Two evolutionary units were identified in G. lawrencei var. farreri. Based on the cpDNA data, phylogenetic relationships of haplotypes showed that the southeast and northwest areas were two geographically structured groups (Figures 1, 2). Consistent with the plastid data, the Bayesian clustering from the nuclear SSR data also indicated that populations of G. lawrencei var. farreri were divided into significantly geographical southeast and northwest groups (Figure 4), although the specific areas from the two kinds of data were not completely consistent. The inconsistence revealed by matrilineal and nuclear inheritance may be caused by pollen flow or introgression.
Mountains usually act as dispersal barriers for plants. The mountains in the distribution area of G. lawrencei var. farreri, for example, the Bayan Har Mountains and Tanggula Mountains, run northwest-southeast. If the mountains are dispersal barriers for G. lawrencei var. farreri, as, for example, the Rocky Mountains are for the barn owl (Machado et al., 2018), the plant populations would be divided into geographically northeast and southwest groups rather than the observed southeast and the northwest groups. The actual pattern could be explained by several hypotheses. One hypothesis is the two-refugia pattern of G. lawrencei var. farreri. Independent refugia during glaciation events is conducive to allopatric divergence (Petit et al., 2003). Recolonization after glaciation from independent refugia at the southeast edge of the QTP and the northwest plateau platform may have formed two genetic clades that laid the foundation for the current genetic structure. The other hypothesis is that the mountains in the QTP platform usually had low relative elevations and did not act as geographical barriers to dispersal (Yu et al., 2017). Plants that were pollinated by insects or animals and dispersed by the wind, as in G. lawrencei var. farreri (Hou et al., 2009), could spread their genetic information across the mountains at short geographical distances.
We also detected high genetic divergence and significant correlations between genetic and geographical distances in G. lawrencei var. farreri. Considering that the isolation by mountains in the plateau platform is limited for plants, we could infer that the independent refugia pattern had a significant influence on the evolutionary history of G. lawrencei var. farreri.

Evolutionary History
Molecular dating based on cpDNA data indicated that the southeast and northwest groups diverged around 4.89 Ma when the QTP has uplift already (Favre et al., 2015). Therefore, the divergence was not caused by the uplift of the QTP, but maybe caused by environmental changes, glaciation or geographic isolation after expansion. Compared with plastid data, nuclear data reflect biparental inheritance and can yield information on evolutionary events occurring via pollen flow. The ABC analysis  Tested scenarios with effective population size changes and their posterior probabilities. SE represents the southeast group, and NW represents the northwest group. N (N1, N2, . . .) represents effective population size, and t (t1, t2, and t3) represents time. Posterior probabilities were estimated based on a polychotomous logistic approach. More details are given in Section "Materials and Methods." based on the cpDNA and nuclear data showed that the two groups were not of independent origin; the northwest group was indicated as having diverged from the southeast group. When generation time of G. lawrencei var. farreri was set to 3 years, the divergence occurred about 22,950 years ago, when the QTP was in the Last Glacial period (Zheng et al., 2002). The analysis hinted that the northwest group may be the in situ remnant of the last colonization from the southeast group, possibly the remnant of the last colonization after glaciation from the southeast edge of the QTP.
The ABC analysis rejected the expansion of effective population size in both groups (Scenario 6) and only favored the scenario in which southeast groups had undergone population expansion after the northwest group had diverged. However, EBSP analysis based on cpDNA indicated that the effective population sizes of the two groups and the whole species were nearly constant during the most recent 200 Kya. Recent population expansions or declines were also detected in very limited areas (MQ, and MDa and XC, respectively). The demographic histories of these populations should be explored further.
The genetic dataset indicated recolonization routes for G. lawrencei var. farreri. The main route of recolonization from the southeast edge of the QTP after the most recent glaciation was northward, along the edge of the QTP (e.g., the Hehuang Valley) to the northern end of the distribution area of G. lawrencei var. farreri. This northward route was described by Yu et al. (2017) based on species distribution models of 20 plant species in the QTP. However, recolonization from the plateau platform of the QTP tended to occur both northward and southward. The recolonization from two locations in different directions led to population mixing. Populations such as GD, MQ, HN, JD, and ChD, located at the adjacent region of the two geographically structured groups, had intermediate cluster membership and thus acted as "melting points." AUTHOR CONTRIBUTIONS P-CF conceived and designed the experiments, collected the samples, performed the experiments, contributed reagents, materials, analysis tools, and wrote the paper. H-YY contributed reagents, materials, analysis tools, prepared figures and/or tables. Q-WL and H-MC performed the experiments and analyzed the data. S-LC analyzed the data and reviewed drafts of the paper.

FUNDING
We are grateful to financial support provided by the National Natural Science Foundation of China (Grant No. 31600296