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 F_ST 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 (Favre et al., 2016). The diversification within Gentiana has been influenced by geological and climatic changes in the QTP (Favre et al., 2016), 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.

Materials and Methods

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.

TABLE 1

TABLE 1. Summary genetic statistics for Gentianan lawrencei var. farreri.

FIGURE 1

FIGURE 1. Geographical distribution of 16 chloroplast haplotypes (H1–H16) identified in Gentiana lawrencei var. farreri. (A) The distributions of G. lawrencei var. farreri on the Qinghai-Tibetan Plateau. (B) Geographical distribution of haplotypes across sampled populations. Pie charts display haplotype frequencies in each locality. (C) Phylogeny of the 16 chloroplast haplotypes detected in G. lawrencei var. farreri. The detailed tree with outgroups is presented in Supplementary Figure S1.

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 ycf1 (Forward: TGCTTGTATTATTACTTTGTGG; Reverse: TTGGATTGGTATTAGTCTGGAT) and was verified in the sampled populations. The ycf1 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₂, 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).

All of the samples were genotyped at 10 microsatellite loci (Law4, Law5, Law24, Law32, Law37, Law41, Law45, Law57, Law71, and Law77), which were amplified and genotyped easily, as described by Sun et al. (2018a). The forward primer 5′ ends were labeled with a fluorescent dye (FAM, HEX, or ROX) [Sangon Biotech (Shanghai) Co., Ltd.]. PCR reactions and profiles followed Sun et al. (2018a). The PCR products were subsequently detected using a 3730XL Genetic Analyzer Sequencer (Applied Biosystems). Allele sizing was performed using GENEMAPPER version 4.0 (Applied Biosystems) by comparing alleles with a GeneScan-500LIZ size standard (Applied Biosystems).

Data Analysis

Chloroplast Sequence Data

Sequences were aligned and edited with GENEIOUS PRO 3.5.6 (Kearse et al., 2012). Haplotypes were identified in DnaSP 5.1 (Librado and Rozas, 2009) and were deposited in GenBank (accession nos. MH481820–MH481835). Gene diversity (h), nucleotide diversity (π) indices, F_ST and a hierarchical AMOVA were performed in ARLEQUIN 3.5 (Excoffier and Lischer, 2010). The coefficients of differentiation G_ST were calculated using the software PERMUT (Pons and Petit, 1996) to estimate differentiation among populations. Spatial genetic structure of cpDNA haplotypes was analyzed with SAMOVA 1.0 (Dupanloup et al., 2002).

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. (Wang et al., 2017) 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¹ 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 ycf1 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.

Results

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 diverged approximately 4.89 Ma (95% highest posterior density: 1.61–10.77 Ma).

FIGURE 2

FIGURE 2. TCS statistical parsimony network of G. lawrencei var. farreri 16 chloroplast haplotypes.

FIGURE 3

FIGURE 3. Majority rule consensus phylogenetic tree of the 16 chloroplast haplotypes detected in G. lawrencei var. farreri based on Bayesian inference. Numbers on the branches indicate Bayesian posterior probabilities. Node age estimates are marked with black arrows. Grey bars represent 95% highest posterior densities.

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.

TABLE 2

TABLE 2. Analyses of molecular variance in G. lawrencei var. farreri based on 10 SSR loci and cpDNA. d.f., degrees of freedom.

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² = 0.485, P < 0.01).

FIGURE 4

FIGURE 4. Genetic structure across the distribution range of Gentiana lawrencei var. farreri. (A) Geographical distribution of genetic clusters estimated by Bayesian clustering. Pie charts show the frequencies of clusters in each population; colors correspond to each cluster. (B) The bar plot shows the probabilities of ancestral clusters of each sample estimated by Bayesian clustering. The name of each population is shown below the bar plot. The white dash line marks the southeast group and the northwest group mentioned in the text.

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.

FIGURE 5

FIGURE 5. Schematic representation of models tested for G. lawrencei var. farreri population history on the Qinghai-Tibetan Plateau using approximate Bayesian computation. (A) Tested population divergence scenarios and their posterior probabilities. (B) 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.”

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 post-glacial 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) and Liu 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 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), the Open Project of Qinghai Provincial Key Laboratory of Crop Molecular Breeding (Grant #2017-ZJ-Y14), Science and Technology Project of Henan Province, China (Grant No. 162102110097), and Educational Commission of Henan Province (Nos. 16A210033, 17A180008).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Ru-Chong Xu and Shuang Wang of College of Life Science, Luoyang Normal University for the work in sample collecting.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2018.00564/full#supplementary-material

Footnotes

1. ^ http://tree.bio.ed.ac.uk/software/tracer/

References

Bertorelle, G., Benazzo, A., and Mona, S. (2010). ABC as a flexible framework to estimate demography over space and time: some cons, many pros. Mol. Ecol. 19, 2609–2625. doi: 10.1111/j.1365-294X.2010.04690.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Chapuis, M. P., and Estoup, A. (2007). Microsatellite null alleles and estimation of population differentiation. Mol. Biol. Evol. 24, 621–631. doi: 10.1093/molbev/msl191

PubMed Abstract | CrossRef Full Text | Google Scholar

Che, M., Chen, B., Innes, J. L., Wang, G., Dou, X., Zhou, T., et al. (2014). Spatial and temporal variations in the end date of the vegetation growing season throughout the Qinghai–Tibetan Plateau from 1982 to 2011. Agric. Forest Meteorol. 189, 81–90. doi: 10.1016/j.agrformet.2014.01.004

CrossRef Full Text | Google Scholar

Christe, C., Caetano, S., Aeschimann, D., Kropf, M., Diadema, K., and Naciri, Y. (2014). The intraspecific genetic variability of siliceous and calcareous Gentiana species is shaped by contrasting demographic and re-colonization processes. Mol. Phylogenet. Evol. 70, 323–336. doi: 10.1016/j.ympev.2013.09.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Cornuet, J. M., Pudlo, P., Veyssier, J., Dehne-Garcia, A., Gautier, M., Leblois, R., et al. (2014). DIYABC v2.0: a software to make approximate Bayesian computation inferences about population history using single nucleotide polymorphism, DNA sequence and microsatellite data. Bioinformatics 30, 1187–1189. doi: 10.1093/bioinformatics/btt763

PubMed Abstract | CrossRef Full Text | Google Scholar

Csillery, K., Blum, M. G. B., Gaggiotti, O. E., and Francois, O. (2010). Approximate Bayesian Computation (ABC) in practice. Trends Ecol. Evol. 25, 410–418. doi: 10.1016/j.tree.2010.04.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Darriba, D., Taboada, G. L., Doallo, R., and Posada, D. (2012). jModelTest 2: more models, new heuristics and parallel computing. Nat. Methods 9:772. doi: 10.1038/nmeth.2109

PubMed Abstract | CrossRef Full Text | Google Scholar

Dempster, A. P., Laird, N. M., and Rubin, D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. J. R. Stat. Soc. B 39, 1–38.

Google Scholar

Diadema, K., Bretagnolle, F., Affre, L., Yuan, Y.-M., and Médail, F. (2005). Geographic structure of molecular variation of Gentiana ligustica (Gentianaceae) in the maritime and ligurian regional hotspot, inferred from ITS sequences. Taxon 54, 887–894. doi: 10.2307/25065475

CrossRef Full Text | Google Scholar

Drummond, A. J., Ho, S. Y., Phillips, M. J., and Rambaut, A. (2006). Relaxed phylogenetics and dating with confidence. PLoS Biol. 4:e88. doi: 10.1371/journal.pbio.0040088

PubMed Abstract | CrossRef Full Text | Google Scholar

Drummond, A. J., Suchard, M. A., Xie, D., and Rambaut, A. (2012). Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol. Biol. Evol. 29, 1969–1973. doi: 10.1093/molbev/mss075

PubMed Abstract | CrossRef Full Text | Google Scholar

Dupanloup, I., Schneider, S., and Excoffier, L. (2002). A simulated annealing approach to define the genetic structure of populations. Mol. Ecol. 11, 2571–2581. doi: 10.1046/j.1365-294X.2002.01650.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Earl, D. A., and vonHoldt, B. M. (2012). Structure harvester: a website and program for visualizing structure output and implementing the evanno method. Conserv. Genet. Resour. 4, 359–361. doi: 10.1007/s12686-011-9548-7

CrossRef Full Text | Google Scholar

Ebersbach, J., Muellner-Riehl, A. N., Michalak, I., Tkach, N., Hoffmann, M. H., Röser, M., et al. (2017). In and out of the qinghai-tibet plateau: divergence time estimation and historical biogeography of the large arctic-alpine genus Saxifraga L. J. Biogeogr. 44, 900–910. doi: 10.1111/jbi.12899

CrossRef Full Text | Google Scholar

Evanno, G., Regnaut, S., and Goudet, J. (2005). Detecting the number of clusters of individuals using the software structure: a simulation study. Mol. Ecol. 14, 2611–2620. doi: 10.1111/j.1365-294X.2005.02553.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Excoffier, L., and Lischer, H. E. L. (2010). Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under linux and windows. Mol. Ecol. Resour. 10, 564–567. doi: 10.1111/j.1755-0998.2010.02847.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Excoffier, L., Smouse, P. E., and Quattro, J. M. (1992). Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data. Genetics 131, 479–491.

PubMed Abstract | Google Scholar

Favre, A., Michalak, I., Chen, C. H., Wang, J. C., Pringle, J. S., Matuszak, S., et al. (2016). Out-of-Tibet: the spatio-temporal evolution of Gentiana (Gentianaceae). J. Biogeogr. 43, 1967–1978. doi: 10.1111/jbi.12840

CrossRef Full Text | Google Scholar

Favre, A., Päckert, M., Pauls, S. U., Jähnig, S. C., Uhl, D., Michalak, I., et al. (2015). The role of the uplift of the Qinghai-Tibetan Plateau for the evolution of Tibetan biotas. Biol. Rev. 90, 236–253. doi: 10.1111/brv.12107

PubMed Abstract | CrossRef Full Text | Google Scholar

Fu, P. C., Gao, Q. B., Zhang, F. Q., Xing, R., Khan, G., Wang, J. L., et al. (2016). Responses of plants to changes in Qinghai–Tibetan Plateau and glaciations: evidence from phylogeography of a Sibiraea (Rosaceae) complex. Biochem. Syst. Ecol. 65, 72–82. doi: 10.1016/j.bse.2016.01.006

CrossRef Full Text | Google Scholar

Fu, Y. X. (1997). Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics 147, 915–925.

PubMed Abstract | Google Scholar

Gao, Q. B., Li, Y., Gengji, Z. M., Gornall, R. J., Wang, J. L., Liu, H. R., et al. (2017). Population genetic differentiation and taxonomy of three closely related species of saxifraga (Saxifragaceae) from southern tibet and the hengduan mountains. Front. Plant Sci. 8:1325. doi: 10.3389/fpls.2017.01325

PubMed Abstract | CrossRef Full Text | Google Scholar

Goudet, J. (2001). FSTAT, A Program to Estimate and Test Gene Diversities and Fixation Indices (Version 2.9. 3). Available at: www.unil.ch/izea/softwares/fstat.html

Google Scholar

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

CrossRef Full Text | Google Scholar

Hewitt, G. M. (2000). The genetic legacy of the Quaternary ice ages. Nature 405, 907–913. doi: 10.1038/35016000

PubMed Abstract | CrossRef Full Text | Google Scholar

Ho, T. N., and Liu, S. W. (2001). A Worldwide Monograph of Gentiana. Beijing: Science Press.

Google Scholar

Ho, T. N., and Pringle, J. S. (1995). “Gentianaceae,” in Flora of China, Vol. 16, eds Z. Y. Wu and P. H. Raven (Beijing: Science Press), 1–140.

Google Scholar

Hou, Q. Z., Duan, Y. W., Si, Q. W., and Yang, H. L. (2009). Pollination ecology of Gentiana lawrencei var. farreri, a late-flowering Qinghai-Tibet plateau species. Chin. J. Plant Ecol. 33, 1156–1164. doi: 10.3773/j.issn.1005-264x

CrossRef Full Text | Google Scholar

Kearse, M., Moir, R., Wilson, A., Stones-Havas, S., Cheung, M., Sturrock, S., et al. (2012). Geneious basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics 28, 1647–1649. doi: 10.1093/bioinformatics/bts199

PubMed Abstract | CrossRef Full Text | Google Scholar

Khan, G., Zhang, F. Q., Gao, Q. B., Fu, P. C., Zhang, Y., and Chen, S. L. (2018). Spiroides shrubs on qinghai-tibetan plateau: multilocus phylogeography and palaeodistributional reconstruction of Spiraea alpina and S. mongolica (Rosaceae). Mol. Phylogenet. Evol. 123, 137–148. doi: 10.1016/j.ympev.2018.02.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Kropf, M., Comes, H. P., and Kadereit, J. W. (2006). Long-distance dispersal vs vicariance: the origin and genetic diversity of alpine plants in the Spanish Sierra Nevada. New Phytol. 172, 169–184. doi: 10.1111/j.1469-8137.2006.01795.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Leigh, J. W., and Bryant, D. (2015). POPART: full-feature software for haplotype network construction. Methods Ecol. Evol. 6, 1110–1116. doi: 10.1111/2041-210X.12410

CrossRef Full Text | Google Scholar

Li, Y. C., Zhong, D. L., Rao, G. Y., Wen, J., Ren, Y., and Zhang, J. Q. (2018). Gone with the trees: phylogeography of Rhodiola sect. Trifida (Crassulaceae) reveals multiple refugia on the Qinghai-Tibetan Plateau. Mol. Phylogenet. Evol. 121, 110–120. doi: 10.1016/j.ympev.2018.01.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Librado, P., and Rozas, J. (2009). DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25, 1451–1452. doi: 10.1093/bioinformatics/btp187

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J. Q., Duan, Y. W., Hao, G., Ge, X, J., and Sun, H. (2014). Evolutionary history and underlying adaptation of alpine plants on the Qinghai-Tibet Plateau. J. Syst. Evol. 52, 241–249. doi: 10.1111/jse.12094

CrossRef Full Text | Google Scholar

Lu, Z., Chen, P., Bai, X., Xu, J., He, X., Niu, Z., et al. (2015). Initial diversification, glacial survival, and continuous range expansion of Gentiana straminea (Gentianaceae) in the Qinghai–Tibet Plateau. Biochem. Syst. Ecol. 62, 219–228. doi: 10.1016/j.bse.2015.09.005

CrossRef Full Text | Google Scholar

Ma, Z., Liu, H., Mi, Z., Zhang, Z., Wang, Y., Xu, W., et al. (2017). Climate warming reduces the temporal stability of plant community biomass production. Nat. Commun. 8:15378. doi: 10.1038/ncomms15378

PubMed Abstract | CrossRef Full Text | Google Scholar

Machado, A. P., Clément, L., Uva, V., Goudet, J., and Roulin, A. (2018). The rocky mountains as a dispersal barrier between barn owl (Tyto alba) populations in north america. J. Biogeogr. 45, 1288–1300. doi: 10.1111/jbi.13219

CrossRef Full Text | Google Scholar

Matuszak, S., Favre, A., Schnitzler, J., and Muellner-Riehl, A. N. (2016). Key innovations and climatic niche divergence as drivers of diversification in subtropical Gentianinae in southeastern and eastern Asia. Am. J. Bot. 103, 899–911. doi: 10.3732/ajb.1500352

PubMed Abstract | CrossRef Full Text | Google Scholar

Myers, N., Mittermeier, R. A., Mittermeier, C. G., da Fonseca, G. A. B., and Kent, J. (2000). Biodiversity hotspots for conservation priorities. Nature 403, 853–858. doi: 10.1038/35002501

PubMed Abstract | CrossRef Full Text | Google Scholar

Nei, M. (1987). Molecular Evolutionary Genetics. New York, NY: Columbia university press.

Google Scholar

Ni, L., Zhao, Z., Xu, H., Chen, S., and Dorje, G. (2016). The complete chloroplast genome of Gentiana straminea (Gentianaceae), an endemic species to the Sino-Himalayan subregion. Gene 577, 281–288. doi: 10.1016/j.gene.2015.12.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Petit, R. J., Aguinagalde, I., De Beaulieu, J. L., Bittkau, C., Brewer, S., Cheddadi, R., et al. (2003). Glacial refugia: hotspots but not melting pots of genetic diversity. Science 300, 1563–1565. doi: 10.1126/science.1083264

PubMed Abstract | CrossRef Full Text | Google Scholar

Pirie, M. D., Litsios, G., Bellstedt, D. U., Salamin, N., and Kissling, J. (2015). Back to Gondwanaland: can ancient vicariance explain (some) Indian Ocean disjunct plant distributions? Biol. Lett. 11:20150086. doi: 10.1098/rsbl.2015.0086

PubMed Abstract | CrossRef Full Text | Google Scholar

Piry, S., Luikart, G., and Cornuet, J. M. (1999). BOTTLENECK: a computer program for detecting recent reductions in the effective size using allele frequency data. J. Hered. 90, 502–503. doi: 10.1093/jhered/90.4.502

CrossRef Full Text | Google Scholar

Pons, O., and Petit, R. J. (1996). Measuring and testing genetic differentiation with ordered versus unordered alleles. Genetics 144, 1237–1245.

PubMed Abstract | Google Scholar

Pritchard, J. K., Stephens, M., and Donnelly, P. (2000). Inference of population structure using multilocus genotype data. Genetics 155, 945–959.

Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosenberg, N. E. (2004). DISTRUCT: a program for the graphical display of population structure. Mol. Ecol. Notes 4, 137–138. doi: 10.1046/j.1471-8286.2003.00566.x

CrossRef Full Text | Google Scholar

Rousset, F. (2008). GENEPOP’007: a complete reimplementation of the GENEPOP software for Windows and Linux. Mol. Ecol. Notes 8, 103–106. doi: 10.1111/j.1471-8286.2007.01931.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Rybczyński, J. J., Davey, M. R., and Mikuła, A. (eds). (2015). The Gentianaceae: Biotechnology and Applications, Vol. 2. Berlin: Springer. doi: 10.1007/978-3-642-54102-5

CrossRef Full Text | Google Scholar

Sun, H., Zhang, J., Deng, T., and Boufford, D. E. (2017). Origins and evolution of plant diversity in the Hengduan Mountains, China. Plant Divers. 39, 161–166. doi: 10.1016/j.pld.2017.09.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, S. S., Fu, P. C., Cheng, Y. W., Zhou, X. J., and Han, J. M. (2018a). Characterization and transferability of microsatellites for Gentiana lawrencei var. farreri (Gentianaceae). Appl. Plant Sci. 6:e1015. doi: 10.1002/aps3.1015

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, S. S., Fu, P. C., Zhou, X. J., Cheng, Y. W., Zhang, F. Q., Chen, S. L., et al. (2018b). The complete plastome sequences of seven species in Gentiana sect. Kudoa (Gentianaceae): insights into plastid gene loss and molecular evolution. Front. Plant Sci. 9:493. doi: 10.3389/fpls.2018.00493

PubMed Abstract | CrossRef Full Text | Google Scholar

Tajima, F. (1989). Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 123, 585–595.

Google Scholar

van Oosterhout, C., Hutchinson, W. F., Wills, D. P. M., and Shipley, P. (2004). Micro-checker: software for identifying and correcting genotyping errors in microsatellite data. Mol. Ecol. Notes 4, 535–538. doi: 10.1111/j.1471-8286.2004.00684.x

CrossRef Full Text | Google Scholar

Wang, X., Yang, N., Su, J., Zhang, H., and Cao, X. (2017). The complete chloroplast genome of Gentiana macrophylla. Mitochondr. DNA B 2, 395–396. doi: 10.1080/23802359.2017.1347831

PubMed Abstract | CrossRef Full Text | Google Scholar

Wen, J., Zhang, J. Q., Nie, Z. L., Zhong, Y., and Sun, H. (2014). Evolutionary diversifications of plants on the Qinghai-Tibetan Plateau. Front. Genet. 5:4. doi: 10.3389/fgene.2014.00004

CrossRef Full Text | Google Scholar

Xing, Y., and Ree, R. H. (2017). Uplift-driven diversification in the Hengduan Mountains, a temperate biodiversity hotspot. Proc. Nat. Acad. Sci. U.S.A. 114, E3444–E3451. doi: 10.1073/pnas.1616063114

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, A. M., Sun, J., Han, H., Shi, X. L., Xu, G. Q., and Zhang, X. R. (2012). Chemical constituents from Gentiana farreri Balf. f. Chin. Tradit. Patent Med. 4, 506–508.

Google Scholar

Yu, H., Zhang, Y., Wang, Z., Liu, L., Chen, Z., and Qi, W. (2017). Diverse range dynamics and dispersal routes of plants on the Tibetan Plateau during the late Quaternary. PLoS One 12:e0177101. doi: 10.1371/journal.pone.0177101

PubMed Abstract | CrossRef Full Text | Google Scholar

Yuan, Y. M., Kupfer, P., and Doyle, J. J. (1996). Infrageneric phylogeny of the genus Gentiana (Gentianaceae) inferred from nucleotide sequences of the internal transcribed spacers (ITS) of nuclear ribosomal DNA. Am. J. Bot. 83, 641–652. doi: 10.2307/2445924

CrossRef Full Text | Google Scholar

Zhang, D. F., Fengquan, L., and Jianmin, B. (2000). Eco-environmental effects of the Qinghai-Tibet Plateau uplift during the Quaternary in China. Environ. Geol. 39, 1352–1358. doi: 10.1007/s002540000174

CrossRef Full Text | Google Scholar

Zhang, X. L., Wang, Y. J., Ge, X. J., Yuan, Y. M., Yang, H. L., and Liu, J. Q. (2009). Molecular phylogeny and biogeography of Gentiana sect, Cruciata (Gentianaceae) based on four chloroplast DNA datasets. Taxon 58, 862–870.

Google Scholar

Zheng, B. X., Xu, Q. Q., and Shen, Y. P. (2002). The relationship between climate change and Quaternary glacial cycles on the Qinghai-Tibet: review and speculation. Quat. Int. 97, 93–101. doi: 10.1016/S1040-6182(02)00054-X

CrossRef Full Text | Google Scholar