Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Open Access

Peer-reviewed

Research Article

Phylogeography of Rhodiola kirilowii (Crassulaceae): A Story of Miocene Divergence and Quaternary Expansion

Abstract

The evolution and current distribution of the Sino-Tibetan flora have been greatly affected by historical geological events, such as the uplift of the Qinghai-Tibetan Plateau (QTP), and Quaternary climatic oscillations. Rhodiola kirilowii, a perennial herb with its distribution ranging from the southeastern QTP and the Hengduan Mountains (HM) to adjacent northern China and central Asia, provides an excellent model to examine and disentangle the effect of both geological orogeny and climatic oscillation on the evolutionary history of species with such distribution patterns. We here conducted a phylogeographic study using sequences of two chloroplast fragments (trnL-F and trnS-G) and internal transcribed spacers in 29 populations of R. kirilowii. A total of 25 plastid haplotypes and 12 ITS ribotypes were found. Molecular clock estimation revealed deep divergence between the central Asian populations and other populations from the HM and northern China; this split occurred ca. 2.84 million year ago. The majority of populations from the mountains of northern China were dominated by a single haplotype or ribotype, while populations of the HM harbored both high genetic diversity and high haplotype diversity. This distribution pattern indicates that HM was either a diversification center or a refugium for R. kirilowii during the Quaternary climatic oscillations. The present distribution of this species on mountains in northern China may have resulted from a rapid glacial population expansion from the HM. This expansion was confirmed by the mismatch distribution analysis and negative Tajima's D and Fu's FS values, and was dated to ca. 168 thousand years ago. High genetic diversity and population differentiation in both plastid and ITS sequences were revealed; these imply restricted gene flow between populations. A distinct isolation-by-distance pattern was suggested by the Mantel test. Our results show that in old lineages, populations may harbour divergent genetic forms that are sufficient to maintain or even increase overall genetic diversity despite fragmentation and low within-population variation.

Citation: Zhang J-Q, Meng S-Y, Rao G-Y (2014) Phylogeography of Rhodiola kirilowii (Crassulaceae): A Story of Miocene Divergence and Quaternary Expansion. PLoS ONE 9(11): e112923. https://doi.org/10.1371/journal.pone.0112923

Editor: Zhong-Jian Liu, The National Orchid Conservation Center of China; The Orchid Conservation & Research Center of Shenzhen, China

Received: July 23, 2014; Accepted: October 16, 2014; Published: November 12, 2014

Copyright: © 2014 Zhang et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: The authors confirm that all data underlying the findings are fully available without restriction. All sequences data are available from the GenBank database with accession numbers KM459552-KM459594.

Funding: This research was supported by National Natural Science Foundation of China (Grant no. 31470313). JQZ has been supported by a scholarship for exchanging PhD candidates from the China Scholarship Council (No. 201206010113) for his studies at the Smithsonian Institution. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Introduction

The evolution and genetic structure of extant species have been greatly affected by both geological history and climate oscillations [1], [2]. The Quaternary ice age played a critical role in shaping the current distributions and geographic structure of the genetic diversity of organisms in the North Hemisphere [1], [3]. Many phylogeographic studies of various species have shown that the spatial structure of genetic diversity is strongly associated with refugia during glaciations [4]. Although there was no major continental ice sheet in eastern Asia during the Quaternary [5], the continuing uplift of the Qinghai-Tibet Plateau (QTP) from the middle Tertiary has greatly changed the geology and topography of Eastern Asia. The QTP is the highest and largest plateau in the world [6], which harbors more than 12,000 species from 1500 genera [7]–[9]. Its uplifting has affected the distribution and genetic diversity of organisms found in this area [10]–[12]. Although phylogeographic studies on various plant groups from diverse geographic regions in China have increased greatly in recent decades, few of them have tried to disentangle the relative influences of both geographic history (e.g., the uplift of the QTP) and Pleistocene climatic oscillations in shaping the current distribution and genetic pattern of species in the QTP and adjacent area [13]–[15].

The Hengduan Mountains (HM), which are located at the southern and eastern edge of the QTP [16]–[18], have been considered the core area of the Himalayan biodiversity hotspot, which is of global importance [19]–[22]. It has long been thought that the severe alteration in the topology and climate of this area since the uplift of the QTP between the Miocene and the Quaternary triggered the origin of many plant taxa, both genera and species [23]–[25]. The Hengduan Mountains have also been suggested to have acted as a refugium for plants on the QTP and in the north temperate zone of eastern Asia during the Quaternary ice age [7]. More than ten plant species in the QTP and adjacent areas have been studied, focusing on their geographic patterns of genetic diversity. The results have confirmed that the Hengduan Mountains were an important refugium during the Quaternary ice age [26]–[31], although it is suspected that some species survived in ice-free areas on the plateau itself during the glaciation [32]–[34]. Many plant species distribute from the HM to northern China (e.g., Rhodiola kirilowii (Regel) Maxim., Incarvillea sinensis Lam., Caragana jubata (Pall.) Poir., Hypecoum leptocarpum Hook. f. et Thoms.), some of which could have originated and expanded northeast from the Hengduan Mountains [35], [36]. It is believed that northern China was not covered by an ice sheet during the Quaternary [1], [37], so determining whether species found in both southwestern China (QTP and adjacent area) and northern China experienced northerly expansion and southerly retreat during the Quaternary climate oscillation is of great interest. Although intensive phylogenetic and phylogeographic investigations of plants of the QTP and adjacent area have been conducted in recent decades [23], [29], [38]–[43], few studies have discussed the floristic relationship between the HM and the vast adjacent northern plain in northern China, which is dominated by homogeneous grassland [44]. Furthermore, most studies have focused on tree species, while herbaceous species have received much less attention [31], [42]. Thus, a phylogeographic study of an herbaceous species with a distribution that covers both the QTP and its adjacent areas and northern China in order to understand the evolutionary history of plants that exhibit similar distribution patterns is particularly valuable.

Habitat fragmentation could dramatically influence the demographic history of plant species [45], [46]. An increasing number of studies indicate that fragmentation could reduce genetic diversity and increase genetic divergence between populations because of restricted gene flow, inbreeding and genetic drift in isolated small populations [47]. On the other hand, global genetic variation could be maintained or even increased by a fragmented population, and allopatric speciation can occur due to genetic isolation [15], [48]. The genetic and ecological consequences of habitat fragmentation depend on other factors that affect genetic diversity: population size, gene flow and the timescale of fragmentation [49]. Thus, analysis of genetic divergence and gene flow between fragmented populations would help us to understand processes of speciation and extinction resulting from environmental changes. Besides, genetic diversity studies of naturally fragmented populations may reveal the consequences of population fragmentation over long periods of time and provide references for predicting consequences of fragmentation caused by human activity [50].

Rhodiola L. (Crassulaceae) consists of about 70 species mainly found at high- elevation in cold regions of the Northern Hemisphere, with the QTP and adjacent area representing its distribution center [51], [52]. A recent study [53] explored the origin and diversification history of Rhodiola and examined the biogeographic relationships between the QTP and other regions of the Northern Hemisphere. The results showed that the origin and diversification times of this group are largely correlated with the extensive uplifts of the Qinghai-Tibetan Plateau, and the ancestral area of Rhodiola is on the QTP, from where it dispersed to other regions of the Northern Hemisphere. Rhodiola kirilowii, one of the most widespread species of this genus, is distributed through the HM, northwest to the QTP platform and the Tianshan Mountains, and northeast to the plain of northern China [51]. Its distribution is naturally fragmented, especially in northern China, as the species often grows on alpine screes at the forest margin, usually 2000–5000 m in elevation. A range-wide phylogeographic analysis of R. kirilowii would be very helpful for exploring the floristic relationship between the QTP plus its neighboring area and the vast adjacent northeastern plain in northern China, and to reveal the consequences of natural habitat fragmentation. Our objectives included: (1) inferring the historical demography that could explain the current distribution pattern of R. kirilowii by examining the spatial pattern of plastid DNA (pDNA) and nrDNA ITS variation; (2) assessing the floristic relationship between the Hengduan Mountains and northern China; (3) revealing the consequences of natural habitat fragmentation on genetic diversity of R. kirilowii.

Materials and Methods

Ethics statement

No specific permits were required for the described locations in China because all researchers collecting the samples had introduction letters from College of Life Sciences, Peking University, Beijing. The field studies did not involve endangered or protected species. All locations of R. kirilowii populations sampled were shown in Table S1.

Population sampling

Population sampling was conducted throughout the whole distribution range of R. kirilowii during the summers of 2009 to 2012, with the exact GPS coordinates and altitude shown in Table S1. Fresh leaves were collected from 29 populations and, with few exceptions, 10–12 individuals were sampled from each population. To avoid sampling closely related plants, the individuals sampled were at least 30 m apart. In total, 306 individuals of R. kirilowii were sampled and leaves were dried with silica gel. In addition, four individuals of R. rosea from different populations were sampled as outgroups.

DNA extraction, PCR amplification, cloning and sequencing

DNA was extracted from silica-gel dried leaves with a Plant Genomic DNA Kit (TianGen Biotech, Beijing, China). The amplification primers used for ITS were ITS-1 and ITS-4 [52], for trnL-F were c and f [54] and for trnS-G were trnS and trnG [55]. Polymerase Chain Reaction (PCR) was performed in 20 µl reaction mixture volumes containing 2 µl 10× buffer, 0.5 µl of each primer, 0.4 µl of dNTP mixture, 1 U of Taq polymerase (TianGen Biotech, Beijing, China) and 1 µl template genomic DNA. The PCR cycling programs followed Liu et al. [56]. Direct sequencing was conducted for the pDNA fragments using the same primers as for the amplification. Most of the ITS samples were also directly sequenced, although a few of them were ligated onto pGEM-T Easy Vector using a Promega Kit (Promega Corporation, Madison, WI, USA). Plasmids containing the right insertion were chosen for sequencing reactions. All sequencing was conducted with a 3730 automatic DNA sequencer by Biomed Biotech, Beijing, China. Contigs were edited and assembled using ContigExpress (a component of Vector NTI Suite 6.0, InforMax). Sequence alignment was mostly done using ClustalW version 1.7 [57] and checked by eye in BioEdit version 7.0.1.

Phylogeny and divergence time

Four accessions from R. rosea were chosen as the outgroup for the phylogenetic analysis based on previous phylogeny studies [58] and morphological classifications. Sequences from two pDNA regions (trnL-F and trnS-G) were concatenated into a matrix. Chloroplast haplotypes and ITS ribotypes were determined from both indels and nucleotide substitutions using DnaSP v5 [59], respectively. GenBank accession numbers of unique sequences of ITS and the two plastid fragments were KM459552–KM459594. Maximum parsimony and Bayesian inference analysis were implemented to reconstruct the phylogenetic relationships among haplotypes. Parsimony analysis was conducted using PAUP* version 4.0b10 [60]. In all parsimony analysis all characters were weighted equally and indels were treated as missing data. Heuristic searches with MULTREES and TBR branch swapping were conducted. Starting trees were constructed using 1000 replicates of random addition sequences. Support for individual nodes was assessed by bootstrapping [61]. For the bootstrap analysis PAUP* was set to run 1000 replicates with ten replicates of random addition sequences and NNI branch swapping, saving a maximum of 1000 trees of the random addition replicates. Nucleotide substitution model (TPM2uf+G for plastid data set, and SYM+G for the ITS data set) parameters were determined for Bayesian inference (BI) using the Akaike information criterion (AIC) in Modeltest version 3.7 [62], [63]. Bayesian inference was conducted using MrBayes version 3.2.1 (nst  = 6, rates  =  gamma for the plastid data set, nst  = 6, rates  =  gamma for the ITS data set) [64]. The search was started from a random tree using four Metropolis-coupled Markov chain Monte Carlo chains with 5,000,000 generations. The sampling rate of trees was every 1000 generations. To assess the convergence of two runs, the average standard deviation of split frequencies was used. The first ca. 10% of generations was discarded as the burn-in. The remaining data were used to construct a 50% majority rule consensus tree and the proportion of bifurcations found in this consensus tree was given as posterior clade probabilities (PP) to estimate the robustness of the BI trees.

In order to detect genealogical relationships among sequences with shallow genetic divergences, we constructed pDNA and ITS haplotype networks with NETWORK ver. 4.2.0.1 [65]. We hypothesized that both site mutations and indels evolved with equal likelihood and each indel was assumed to have originated independently of all other indels.

To infer the divergence time between lineages, we used the ITS data set to conduct a dating analysis with the BEAST software [66]. To examine the constancy of molecular evolution rate among lineages of the phylogenies, a likelihood ratio test was performed with PAUP version 4.0b10, in which likelihood scores of the trees with and without an enforced molecular clock were compared. Significance was calculated by comparing twice the difference in log likelihoods to a χ2 distribution with n-2 degrees of freedom, where n is the number of taxa (haplotypes plus outgroups). A molecular clock could not be rejected for the ITS data (2logeLR  = 20.84, df  = 14, P = 0.11). We then used BEAST v.1.7.5 [66] to estimate divergence times and confidence intervals. We used GTR + G substitution models, a fixed molecular clock for the ITS data set, a constant population size coalescent tree prior and a UPGMA starting tree for both data sets. We sampled all parameters once every 2000 generations from 20,000,000 MCMC generations with a burn-in of 5,000,000 generations. We then used Tracer [66] to examine convergence of chains to the stationary distribution and the analysis was repeated before combining the two independent runs. As no fossils of R. kirilowii and its relatives were available, we used substitution rates to estimate divergence times. The substitution rate of nrITS in shrubs and herbs varies from 3.46×10−9 to 8.69×10−9 s s−1 yr−1 [67]. Considering the uncertainty of the rates, we used a normal distribution prior to cover these ranges with a 95% confidence interval, which is a mean of 6.075×10−9, and SD of 1.590×10−9.

Population genetics

Haplotype diversity (h) and nucleotide diversity (π) for each population were calculated using DnaSP v5 [59]. Average gene diversity with a population (HS), total gene diversity (HT) and between population divergence (GST, NST) were estimated by the program PERMUT (available at www.pierroton.inra.fr/genetics/labo/Software/) with 1000 permutations tests. A higher NST than GST indicates the presence of phylogeographic structure, which means that closely related haplotypes will be found more often in the same area than more divergent haplotypes [68]. Measures of DNA divergence between the three geographically defined groups (i.e. the Hengduan Mountains, northern China, and central Asia) (FST) were calculated with Arlequin version 3.5 [69], and the significance was tested using 10,000 permutations. To test the significance of isolation by distance between populations, the Mantel test with 1,000 random permutations on matrices of pair wise population FST values and the geographical distances was performed using Arlequin ver. 3.5 [69].

The spatial genetic pattern was examined by spatial analysis of molecular variance (SAMOVA) using SAMOVA v. 1.0 [70]. The program seeks the composition of a user-defined number (K) of groups of geographically adjacent populations that maximize the FCT value (differences among groups of populations). We determined the value of K by repeatedly running the program with K = 2–10, and choosing the one which gave the maximum FCT value. The amount of variation among populations within a region and within a population was calculated by the hierarchical analysis of molecular variance (AMOVA) framework carried out using Arlequin version 3.5 [69]. Significant differences were identified by a nonparametric permutation procedure with 1000 permutations.

We calculated Tajima's D and Fu's Fs to identify signatures of demographic expansion of each clade recovered based on the plastid data set [71], [72]. If an expansion hypothesis was true, the D and Fs statistic should have large negative values as a result of an excess of rare new mutations. We also conducted a mismatch distribution analysis [73], [74] to detect the population expansion of R. kirilowii. We pooled all haplotypes of each clade together without considering their frequencies because population structure has little effect on mismatch distribution [75]. We used 1000 parametric bootstrap replicates to test the fitness of observed mismatch distributions to the theoretical distribution under a sudden expansion model [74], [76], examining the sum of squared deviations (SSD) between observed and expected mismatch distributions, and the raggedness index (HRag) of Harpending [77]. When we identified that a group had experienced expansion, we used the parameter-value for the mode of the mismatch distribution (τ) to estimate time (in generations) since expansion with the equation t = τ/2u [74], [75]. Here, u was calculated as u = µ×k×g, where µ is the number of substitutions per site per year (s s−1 yr−1), k is the average sequence length in the present study and g is the generation time in years (age when first reproduction occurs). In the present study, our value for k was 1549 bp. The pDNA substitution rates for most angiosperm species have been estimated to be in the range 1−3×10−9 substitutions per site per year (s s−1 yr−1) [78]. As the two fragments we used are non-coding regions of the pDNA genome, we assumed an evolutionary rate of 1.52×10−9 s s−1 yr−1 for the plastid data set [79]. Ten years was used as an approximation for g (J. Q. Zhang, pers. obs.).

Results

Plastid sequence variation and haplotype distribution

The length of aligned sequences was 941 bp for trnL-F, and 609 bp for trnS-G. After combining the two data sets, we detected 25 haplotypes (H1–H25) from the 306 sampled individuals of R. kirilowii, determined by 36 single nucleotide polymorphic sites and three indels. The length of these haplotypes varied between 1424 and 1491 bp. Eight haplotypes were unique to single populations, while the rest were shared by different populations (Table S2). Among these, H2 was the most common haplotype – occurring in 13 of the 29 populations – followed by H1, H3, H5, and H12 – occurring in five populations (Table S2). At the species level, the estimated haplotype diversity was h = 0.913. Haplotype diversity varied across populations, ranging from 0.000 to 0.788, with the LHX and QES populations having the highest h values (Table S1). Nucleotide diversity was π = 0.0036 at the species level, and it ranged from 0.000 to 0.0039, with population BS having the highest π value (Table S1). Within-population gene diversity (HS) was much lower than total gene diversity (HT) (0.367 and 0.914, respectively, Table 1). Of the 29 sampled populations, 11 harbored only a single haplotype, while the other 18 were polymorphic, with two or more haplotypes (Fig. 1; Table S1). Population SK harbored six haplotypes which was the most in any of the populations.

thumbnail
Download:
Figure 1. A map of the sampling sites and the geographic distribution of Rhodiola kirilowii haplotypes.

Pie charts show the proportion of chlorotypes within each population. The numbers beside the circles represent population numbers listed in Table S1. Dash line on the map indicate the distribution area of R. kirilowii.

https://doi.org/10.1371/journal.pone.0112923.g001

thumbnail
Download:
Table 1. Genetic diversity and genetic differentiation of 29 populations of Rhodiola kirilowii at the species level and group levels.

https://doi.org/10.1371/journal.pone.0112923.t001

NrDNA ITS variation and ribotype distribution

The length of aligned ITS sequences was 618 bp. Twelve different ITS sequences (ribotypes) were recovered (R1–R12), three of which occurred in not more than one population. These ribotypes were determined by a total of 30 substitutions, and their length was in the range 610–618 bp. Among the ribotypes determined, only three occurred in a single population, while the other nine were shared by at least two populations (Table S3). The most widespread ribotypes were R7 and R8, both of which were found in eight populations (Table S3). Two populations in Xinjiang (WLMQ and NS) have the same ribotype, R11 (Table S1). Ribotype diversity was estimated to be h = 0.866 at the species level, ranging from 0.000 to 0.667 in different populations. Population DLS3 has the highest h value (Table S1). At the species level, nucleotide diversity was 0.010 (π = 0.010), but it varied across populations, ranging from 0.000 to 0.0016, with the LHS population having the highest π value (Table S1). Like the pDNA data, within-population gene diversity (HS) was much lower than total gene diversity (HT) (0.156 and 0.879, respectively; Table 1). Notably, 20 of the 29 populations only harbored a single ribotype, and only nine populations harbored more than one ribotype. Population LHS had the most ribotypes: three of them.

Phylogenetic relationships and lineages divergence

For the plastid data, our phylogenetic analyses based on MP and Bayesian methods yielded trees with largely congruent topologies (Fig. 2a). The monophyly of R. kirilowii was well supported. The 25 haplotypes clustered into three multiple-haplotype clades (A, B, and C) (Fig. 2a) and four monotypic ones. The haplotype network (Fig. 2b) depicted relationships between haplotypes in more detail. The three haplotypes (H22, 23, 24) in the central Asian populations (WLMQ and NS) resided at the base of the phylogeny. Three haplotypes found in northern China were nested in clade A, which comprised both the HM and northern China haplotypes.

thumbnail
Download:
Figure 2. Phylogenetic relationships obtained by analysis of pDNA haplotypes.

(a) Bayesian tree with numbers above the branches indicating bootstrap values greater than 50% for MP analysis and Bayesian posterior probabilities. (b) NETWORK-derived genealogical relationship. The sizes of the circles in the network are proportional to the observed frequencies of the haplotypes. The small black bars represent mutation steps and the red dots represent missing chlorotypes. For both subfigures, different colors represent different geographic origins of haplotypes: red, the Qinghai-Tibetan Plateau and Hengduan Mountains; blue, northern China; green, central Asia.

https://doi.org/10.1371/journal.pone.0112923.g002

For the ITS data, phylogenetic relationships reconstructed by Bayesian and MP methods were consistent (Fig. 3a). All 12 ribotypes formed a well-supported clade. Ribotype 11 from central Asia was highly diverged from the other ribotypes, and the other 11 ribotypes formed a well-support clade (BP  = 100, PP  = 1.00); several major clades with high bootstrap values and posterior probabilities were recognized (I, II, and III). The ribotype network (Fig. 3b) contained the same relationships as the phylogenetic trees (Fig. 3a), although the ribotype relationships were shown with more details. The dating analysis showed that Ribotype 11 diverged from the other ribotypes at 2.84 Mya (95% HPD: 1.75–4.10 Mya), i.e. in the late Pliocene (Fig. 4). Further diversification of the multiple-ribotype clade took place at 0.96 Mya (95% HPD: 0.48–1.56 Mya), during the Pleistocene.

thumbnail
Download:
Figure 3. Phylogenetic relationships obtained by analysis of ITS ribotypes.

(a) Bayesian tree with numbers above the branches indicating bootstrap values greater than 50% for MP analysis and Bayesian posterior probabilities. (b) NETWORK-derived genealogical relationship. The sizes of the circles in the network are proportional to the observed frequencies of the ribotypes. The small black bars represent mutation steps and the red dots represent missing ribotypes. For both subfigures, different colors represent different geographic origins of ribotypes: red, the Qinghai-Tibetan Plateau and Hengduan Mountains; blue, northern China; green, central Asia.

https://doi.org/10.1371/journal.pone.0112923.g003

thumbnail
Download:
Figure 4. Maximum clade credibility tree of the ribotypes generated from BEAST.

Gray bars indicate 95% highest posterior density intervals.

https://doi.org/10.1371/journal.pone.0112923.g004

Population structure

Population differentiation of R. kirilowii was particularly high, with GST values of 0.599 and 0.822 for pDNA and ITS data, respectively (Table 1). A significantly higher value of NST than GST was detected by the permutation tests examining the ITS data (P<0.01), indicating a strong phylogeographic structure in R. kirilowii. In contrast, the permutation tests of pDNA showed that NST was not significantly higher than GST (Table 1). In the SAMOVA analysis for both pDNA and ITS data, the FCT value increased progressively as the value of K increased from two to ten (Fig. S1), making identifying the number K of groups ambiguous. However, the phylogeny based on both pDNA and ITS data sets showed consistency between the three geographic regions (i.e. the Hengduan Mountains, northern China, and central Asia) and lineage divergence. Thus, we divided the 29 populations collected into three geographic groups as mentioned above. Pairwise FST between the three geographic regions varied from 0.5387 to 0.9163, and from 0.1937 to 0.5768 for pDNA and ITS data, respectively (Table 2), suggesting high genetic differentiation between the regions. The Mantel test based on both pDNA and ITS data revealed a significant pattern of isolation-by-distance (pDNA: r = 0.452, P<0.01; ITS: r = 0.489, P<0.01).

thumbnail
Download:
Table 2. Pairwise comparisons of FST among regions estimated from internal transcribed spacer (ITS) sequences (upper part) and pDNA sequences (lower part) of Rhodiola kirilowii.

https://doi.org/10.1371/journal.pone.0112923.t002

Demographic analyses

Mismatch distribution analysis suggested that, under a model of population expansion, only Clade A of the pDNA tree, which included most haplotypes of the QTP and all haplotypes of northern China, showed a unimodal mismatch distribution (Fig. S2). The observed variance (SSD) and the raggedness index for Clade A were not significantly different from the expectation under the population expansion model (Table 4). In addition, Tajima's D and Fu's Fs were also significantly negative (Table 4). Thus, there is a strong indication that Clade A underwent a rapid expansion. The estimated time of this expansion was 0.168 Mya (Table 4).

Discussion

Pliocene intraspecific divergence, Pleistocene range expansions, and glacial refugia of R. kirilowii

Both our phylogenetic and NETWORK analysis showed that the central Asia populations (Population WLMQ and NS) are clearly divergent from the other populations of R. kirilowii. The results of BEAST analysis based on ITS data showed that this divergence took place in the middle to late Pliocene (Fig. 4), during which the most recent uplift of the Qinghai-Tibetan Plateau occurred [80]. Although the uplift of the QTP began around 40 Mya [81]–[84], the recent intensive uplift took place 3.4 to 1.6 Mya and resulted in the QTP reaching its present elevation, i.e. more than 4,000 m [80]. Meanwhile, the recent uplift of the QTP created the Asian monsoon climate, which resulted in the rapid expansion of dry habitats on the platform of the QTP and in the interior of Asia [85]. Rhodiola kirilowii often grows at forest margins, on grassy slopes, often in partial shade near mountain peaks [51]. The vegetation of the QTP before its last intensive uplift was mostly composed of deciduous broadleaved forest and coniferous mixed forest [86], [87], habitats where R. kirilowii usually occurs. Thus, R. kirilowii could have been widely spread from central Asia (i.e. Tianshan Mountain) to the QTP area. However, after the recent intensive uplifts, the flora of the northern and western part of the QTP was replaced by alpine scrub grasslands and desert [88] due to its drier and colder environment. Thus, the distribution of R. kirilowii was fragmented by the intensive uplift of the QTP in the late Pliocene. Pliocene uplift of the QTP has been suggested as a driving force for intraspecific divergence of several alpine species in this region [33], [34], [41], [42], and most recently for two conifer species, Taxus wallichiana Zucc. [89] and Picea likiangensis (Franch.) E. Pritz. [13]. However, most of the previous studies revealed divergence between eastern Himalaya and Hengduan Mountains [13], [89]; in contrast, our results demonstrate a rare divergence between central Asia and the Hengduan Mountains area. Numerous studies on other shrubs, herbs and animal groups with dated molecular phylogenies have also indicated extensive species diversification in the QTP and adjacent area within the Pliocene [15], [23], [32]. Therefore, uplifts of the QTP caused species' distribution fragmentation, which promoted both intraspecific and interspecific divergence on a large scale in this region.

The HM area has long been recognized not only as an important species center for Tertiary elements, but also a principal glacial refugium for many plants in East Asia [26], [29], [31], [41]. A glacial refugium is characterized by high genetic diversity and haplotype uniqueness [90]. These features were clearly revealed in the populations of the HM area, especially populations BM, DF, and LHX (Fig. 1; Table S1). All three of these populations are characterized by high genetic diversity and haplotype uniqueness. Thus, our results confirm the HM area as an important Quaternary glacial refugium, as documented in previous studies [26], [29], [31], [34], [41], [91], [92]. The populations of central Asia diverged in the Pliocene as discussed above, thus they may have survived the Last Glacial Maximum (LGM) in the East Tianshan Mountains, which other studies have also identified as a potential major refugium for plants in northwest China [93].

Few phylogeographic studies have discussed the floristic relationship of HM with the vast adjacent area in northern China, which is dominated by homogeneous grassland [44]. With R. kirilowii's distribution extending from the HM to the mountains of northern China, we have the chance to discuss this issue. Both our pDNA and ITS results show that most populations from northern China harbor only one haplotype or ribotype (Fig. 1). Only three closely related haplotypes (H12, H13, and H 25) and ribotypes (R6, R7, and R12) were found in the northern China groups. In contrast, most populations from the HM area harbored more than one haplotype or ribotype. On the regional level, the total gene diversity of the HM was also higher than that of northern China (Table 1). Furthermore, our phylogenetic analysis showed that the three haplotypes and ribotypes of the northern China group were nested within a large clade which includes most of the HM haplotypes (Figs. 2 & 3). All the above evidence indicates that the northern China populations originated from ones in the HM as a result of population expansion. Our mismatch distribution analysis revealed a unimodal pattern for Clade A of the pDNA phylogenetic tree (Fig. S2), indicating a northeastward range expansion. Furthermore, negative values of Fu's FS and Tajima's D were also observed, which again reinforce the expansion scenario (Table 4). No such range expansions were indicated for the other two clades (Table 4). Using a mean pDNA mutation rate of 1.52×10−9 substitutions per site per year [79], the range expansion of Clade A was estimated to have occurred 168 kya, which fell within the Guxiang (the Penultimate) glacial period of the QTP and the HM area. During the glacial period, the cold-adapted R. kirilowii could respond to these oscillations by moving down the mountains and migrating along the mountain ridges of northern China, e.g., the Qingling Mountains, the Taihang Mountains and the Lvliang Mountains. When the temperature increased after the glacial period, this species retreated into the habitats on mountain tops. The distribution pattern exhibited by R. kirilowii is similar to that of Stellera chamaejasme L., a perennial herb mainly occurring in the HM and northern China [42]. The present distribution of this species in northern China resulted from a rapid post-glaciation expansion from the HM [42]. However, another species with a similar distribution pattern displayed a different pattern: the two highly divergent lineages corresponding to the eastern QTP and away from the QTP (northern China) of Incarvillea sinensis diverged 4.4 Mya, which was much earlier than R. kirilowii and S. chamaejasme [38].

Population structure and genetic diversity of naturally isolated populations

For both pDNA and ITS data sets, we detected a significant high total genetic diversity (pDNA, HT0.914, ITS, HT = 0.879; Table 1) based on 29 populations and 306 individuals of R. kirilowii. Meanwhile, the estimate of genetic differentiation among the 29 populations was also high based on both pDNA and ITS data (FST = 0.598 and 0.826, separately). A high total genetic diversity has also been documented in two previous phylogeographic studies involving two other congeneric species (i.e. R. alsia (Fröd.) S. H. Fu and R. dumulosa (Franch.) S. H. Fu; HT = 0.950 and 0.981, respectively) [91], [94]. However, the total genetic diversities of these three species are all higher than other alpine species found in the QTP and adjacent area, for example, Pedicularis longiflora Rudolph (HT = 0.770) [31], Aconitum gymnandrum Maxim. (HT = 0.739) [34], Incarvillea sinensis Lam. (HT = 0.677) [38], and Taxus wallichiana Zucc. (HT = 0.809) [89]. Thus, in the case of R. kirilowii, the hypothesis that genetic variability decreased as a result of population fragmentation was not supported. Although numerous studies have revealed that habitat fragmentation can reduce genetic diversity due to restricted gene flow, inbreeding and genetic drift in isolated small populations [47], some studies have told another story, namely that overall genetic diversity can be maintained or even increased by fragmentation, as a result of allopatric speciation occurring due to isolation [15], [48], [50], [95]. Rhodiola kirilowii grows exclusively at forest margins, on grassy slopes near mountain peaks at an elevation of 2,000 to 5,600 m, effectively occupying ‘ecological islands’ which are naturally fragmented. These conditions are found especially in the HM area, where severe uplift created high mountains and deeps valleys, providing great opportunities for species like R. krilowii to undergo allopatric divergence. A high genetic differentiation among populations could be the result of spatially isolated populations with restricted gene flow [96], [97]. The scenario of high genetic diversity caused by isolation was also supported by the Mantel test between the pairwise FST and geographic distance (for pDNA, r = 0.452, P = 0.0000; for ITS, r = 0.489, P = 0.0000), indicating a significant and strong isolation-by-distance pattern.

Compared to the high total genetic diversity as mentioned above, the mean within-population genetic diversity at both species and regional level was low (Table 1). The AMOVA analysis results also suggested that most genetic variations should be attributed to among-region or among-population levels for the ITS data (Table 3). Nevertheless, the AMOVA analysis of the pDNA data revealed that within-population variation is greater than among-population variation (Table 3). The pattern was confirmed by the results of the permutation test, which revealed a significantly higher NST value than GST value in the ITS data but not in the pDNA data (Table 1). This phenomenon could be ascribed to repeated population expansion and contraction responding to climate oscillations during the Quaternary, especially in the HM area. Rhodiola kirilowii is a cold-adapted species, so when temperatures rose during the inter-glacial period, populations of this species were able to retreat to mountain peaks, which served as refuges. Such refuges may harbor high genetic diversity, as indicated in this study. However, this high genetic diversity was not revealed in the ITS data set, which suggests frequent gene flow via pollen within populations and between close populations. More strikingly, all individuals from central Asia (populations WLMQ and NS; Table S1) all bear the same ribotype. On the other hand, these populations harbor three pDNA haplotypes. Our results indicated that in R. kirilowii, gene flow via pollen can be very efficient in closely distributed populations, while the dispersal ability of seeds is limited. Once the capsules of R. kirilowii split, the seeds are ejected and dispersed by wind. This kind of dispersal mechanism could be easily hampered by deep valleys and high mountains between populations. At a larger scale, the AMOVA analysis showed that the variation at the regional level accounted for 79.26% and 45.28% of the ITS and pDNA data, respectively (Table 3). Pairwise FST among the three geographic regions also suggested a high level of divergence, which may be the result of the long isolation of these three regions.

thumbnail
Download:
Table 3. Analysis of molecular variance (amova) of chlorotypes and ITS ribotypes for Rhodiola kirilowii populations, partitioned by subspecies and regions, respectively.

https://doi.org/10.1371/journal.pone.0112923.t003

thumbnail
Download:
Table 4. Results of the mismatch distribution analysis and neutrality tests of the three multiple-haplotype pDNA clades.

https://doi.org/10.1371/journal.pone.0112923.t004

Supporting Information

Figure S1.

Correlation between the F statistics and grouping number (K = 2–10) from the SAMOVA results. (a) results based on pDNA haplotypes; (b) results based on ITS ribotypes.

https://doi.org/10.1371/journal.pone.0112923.s001

(TIF)

Figure S2.

Mismatch distribution analyses of the three multiple-haplotype pDNA clades. The histogram of observed mismatch frequencies and the best-fit curve of the sudden expansion model was shown.

https://doi.org/10.1371/journal.pone.0112923.s002

(TIF)

Table S1.

Locations of populations of R. kirilowii sampled, sample sizes (n), frequencies of cpDNA haplotypes and ITS sequences per population, and estimates of haplotype diversity and nucleotide diversity for chlorotypes and ribotypes within populations.

https://doi.org/10.1371/journal.pone.0112923.s003

(DOCX)

Table S2.

Haplotype composition of 29 sampled populations of Rhodiola kirilowii.

https://doi.org/10.1371/journal.pone.0112923.s004

(DOCX)

Table S3.

Ribotypes composition of 29 sampled populations of Rhodiola kirilowii.

https://doi.org/10.1371/journal.pone.0112923.s005

(DOCX)

Acknowledgments

The authors thank two anonymous reviewers for their insightful suggestions, which made the paper much better. We also thank Yiming An, Kai Wu, and Dongqing Zhang for field assistance; and Jing Li and Yongshuai Sun for their assistance in data analysis. The study represents part of Jianqiang Zhang's dissertation research.

Author Contributions

Conceived and designed the experiments: JQZ GYR. Performed the experiments: JQZ. Analyzed the data: JQZ. Contributed reagents/materials/analysis tools: JQZ SYM. Wrote the paper: JQZ GYR.

References

  1. 1. Hewitt G (2000) The genetic legacy of the Quaternary ice ages. Nature 405: 907–913.
  2. 2. Mayr E (1942) Systematics and the Origin of Species. New York: Columbia University Press.
    • 3. Hewitt GM (2004) Genetic consequences of climatic oscillations in the Quaternary. Philosophical Transactions of the Royal Society of London Series B-Biological Sciences 359: 183–195.
    • 4. Avise JC (2009) Phylogeography: retrospect and prospect. Journal of Biogeography 36: 3–15.
    • 5. Zhang DF, Fengquan L, Jianmin B (2000) Eco-environmental effects of the Qinghai-Tibet Plateau uplift during the Quaternary in China. Environmental Geology 39: 1352–1358.
    • 6. Zheng D (1996) The system of physico-geographical regions of the Qinghai-Tibetan (Xizang) Plateau. Science China, Series D Earth Science 39: 410–417.
    • 7. Wu ZY (1987) Origin and evolution of the flora of Tibet. In: Sinica A, editor. Flora Xizangica. Beijing: Science Press. pp.874–902.
      • 8. Wang WC (1993) Vascular Plants of the Hengduan Mountains, Vol. 1. Beijing: Science Press.
        • 9. Wang WC (1994) Vascular Plants of the Hengduan Mountains, Vol. 2. Beijing: Science Press.
          • 10. Tapponnier P, Xu ZQ, Roger F, Meyer B, Arnaud N, et al. (2001) Geology-Oblique stepwise rise and growth of the Tibet plateau. Science 294: 1671–1677.
          • 11. Molnar P, England P, Martinod J (1993) Mantle dynamics, uplift of the Tibetan Plateau, and the Indian monsoon. Reviews of Geophysics 31: 357–396.
          • 12. Wang CS, Zhao XX, Liu ZF, Lippert PC, Graham SA, et al. (2008) Constraints on the early uplift history of the Tibetan Plateau. Proceedings of the National Academy of Sciences of the United States of America 105: 4987–4992.
          • 13. Li L, Abbott RJ, Liu B, Sun Y, Li L, et al. (2013) Pliocene intraspecific divergence and Plio-Pleistocene range expansions within Picea likiangensis (Lijiang spruce), a dominant forest tree of the Qinghai-Tibet Plateau. Molecular Ecology 22: 5237–5255.
          • 14. Liu L, Hao ZZ, Liu YY, Wei XX, Cun YZ, et al. (2014) Phylogeography of Pinus armandii and its relatives: heterogeneous contributions of geography and climate changes to the genetic differentiation and diversification of Chinese white pines. PloS ONE 9: e85920.
          • 15. Xu TT, Abbott RJ, Milne RI, Mao K, Du FK, et al. (2010) Phylogeography and allopatric divergence of cypress species (Cupressus L.) in the Qinghai-Tibetan Plateau and adjacent regions. BMC Evolutionary Biology 10
          • 16. Wu ZY (1988) Hengduan Mountain flora and her significance. Journal of Japanese Botany 63: 297–311.
          • 17. Li BY (1987) On the boundaries of the Hengduan Mountains. Journal of Mountain Research 2: 001.
          • 18. Li BY (1989) Geomorphologic regionalization of the Hengduan mountainous region. Journal of Mountain Research 7: 13–20.
          • 19. Boufford DE, Dijk PPV (2000) South-Central China. In: Mittermeier RA, Myers N, Mittermeier CG, editors. Hotspots: Earth's Biologically Richest and Most Endangered Terrestrial Ecoregions. Mexico: Cemex. pp.338–351.
            • 20. Boufford DE, Dijk PPV, Zhi L (2004) Mountains of southwest China. In: Mittermeier RA, Robles-Gil P, Hoffman M, Pilgrim JD, Brooks TM, et al., editors. Hotspots revisited: Earth's Biologically Richest and Most Endangered Ecoregions. Mexico: Cemex. pp.159–164.
              • 21. Myers N, Mittermeier RA, Mittermeier CG, da Fonseca GAB, Kent J (2000) Biodiversity hotspots for conservation priorities. Nature 403: 853–858.
              • 22. Wilson EO (1992) The Diversity of Life. Cambridge, MA: Harvard University Press.
                • 23. Liu JQ, Wang YJ, Wang AL, Hideaki O, Abbott RJ (2006) Radiation and diversification within the Ligularia-Cremanthodium-Parasenecio complex (Asteraceae) triggered by uplift of the Qinghai-Tibetan Plateau. Molecular Phylogenetics and Evolution 38: 31–49.
                • 24. Sun YS, Wang AL, Wan DS, Wang Q, Liu JQ (2012) Rapid radiation of Rheum (Polygonaceae) and parallel evolution of morphological traits. Molecular Phylogenetics and Evolution 63: 150–158.
                • 25. Wang YJ, Liu JQ, Miehe G (2007) Phylogenetic origins of the Himalayan endemic Dolomiaea, Diplazoptilon and Xanthopappus (Asteraceae: Cardueae) based on three DNA regions. Annals of Botany 99: 311–322.
                • 26. Cun YZ, Wang XQ (2010) Plant recolonization in the Himalaya from the southeastern Qinghai-Tibetan Plateau: Geographical isolation contributed to high population differentiation. Molecular Phylogenetics and Evolution 56: 972–982.
                • 27. Qiu YX, Fu CX, Comes HP (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. Molecular Phylogenetics and Evolution 59: 225–244.
                • 28. Zhang TC, Comes HP, Sun H (2011) Chloroplast phylogeography of Terminalia franchetii (Combretaceae) from the eastern Sino-Himalayan region and its correlation with historical river capture events. Molecular Phylogenetics and Evolution 60: 1–12.
                • 29. Chen SY, Wu GL, Zhang DJ, Gao QB, Duan YZ, et al. (2008) Potential refugium on the Qinghai-Tibet Plateau revealed by the chloroplast DNA phylogeography of the alpine species Metagentiana striata (Gentianaceae). Botanical Journal of the Linnean Society 157: 125–140.
                • 30. Meng LH, Yang R, Abbott RJ, Miehe G, Hu TH, et al. (2007) Mitochondrial and chloroplast phylogeography of Picea crassifolia Kom. (Pinaceae) in the Qinghai-Tibetan Plateau and adjacent highlands. Molecular Ecology 16: 4128–4137.
                • 31. Yang FS, Li YF, Ding X, Wang XQ (2008) Extensive population expansion of Pedicularis longiflora (Orobanchaceae) on the Qinghai-Tibetan Plateau and its correlation with the Quaternary climate change. Molecular Ecology 17: 5135–5145.
                • 32. Jia DR, Abbott RJ, Liu TL, Mao KS, Bartish IV, et al. (2012) Out of the Qinghai-Tibet Plateau: evidence for the origin and dispersal of Eurasian temperate plants from a phylogeographic study of Hippophae rhamnoides (Elaeagnaceae). New Phytologist 194: 1123–1133.
                • 33. Jia DR, Liu TL, Wang LY, Zhou DW, Liu JQ (2011) Evolutionary history of an alpine shrub Hippophae tibetana (Elaeagnaceae): allopatric divergence and regional expansion. Biological Journal of the Linnean Society 102: 37–50.
                • 34. Wang LY, Abbott RJ, Zheng W, Chen P, Wang YJ, et al. (2009) History and evolution of alpine plants endemic to the Qinghai-Tibetan Plateau: Aconitum gymnandrum (Ranunculaceae). Molecular Ecology 18: 709–721.
                • 35. Wang WC (1992) On some distribution patterns and some migration routes found in the eastern Asiatic region. Acta Phytotaxonomica Sinica 30: 1–24.
                • 36. Wang WC (1992) On some distribution patterns and some migration routes found in the eastern Asiatic region. Acta Phytotaxonomica Sinica 30: 97–117.
                • 37. Liu KB (1988) Quaternary history of the temperate forests of China. Quaternary Science Reviews 7: 1–20.
                • 38. Chen ST, Xing YW, Su T, Zhou ZK, Dilcher EDL, et al. (2012) Phylogeographic analysis reveals significant spatial genetic structure of Incarvillea sinensis as a product of mountain building. BMC Plant Biology 12: 58.
                • 39. Meng LH, Yang HL, Wu GL, Wang YJ (2008) Phylogeography of Hippophae neurocarpa (Elaeagnaceae) inferred from the chloroplast DNA trnL-F sequence variation. Journal of Systematics and Evolution 46: 32–40.
                • 40. Wang LY, Ikeda H, Liu TL, Wang YJ, Liu JQ (2009) Repeated range expansion and glacial endurance of Potentilla glabra (Rosaceae) in the Qinghai-Tibetan Plateau. Journal of Integrative Plant Biology 51: 698–706.
                • 41. Yang FS, Qin AL, Li YF, Wang XQ (2012) Great genetic differentiation among populations of Meconopsis integrifolia and its implication for plant speciation in the Qinghai-Tibetan Plateau. PLoS ONE 7
                • 42. Zhang YH, Volis S, Sun H (2010) Chloroplast phylogeny and phylogeography of Stellera chamaejasme on the Qinghai-Tibet Plateau and in adjacent regions. Molecular Phylogenetics and Evolution 57: 1162–1172.
                • 43. Gao Y, Harris A, Zhou S, He X (2013) Evolutionary events in Lilium (including Nomocharis, Liliaceae) are temporally correlated with orogenies of the QT plateau and the Hengduan Mountains. Molecular Phylogenetics and Evolution 68: 443–460.
                • 44. Wu ZY (1980) The Vegetation of China. Beijing: Science Press.
                  • 45. Fahrig L (2003) Effects of habitat fragmentation on biodiversity. Annual Review of Ecology Evolution and Systematics 34: 487–515.
                  • 46. Young A, Boyle T, Brown T (1996) The population genetic consequences of habitat fragmentation for plants. Trends in Ecology and Evolution 11: 413–418.
                  • 47. Aguilar R, Quesada M, Ashworth L, Herrerias-Diego Y, Lobo J (2008) Genetic consequences of habitat fragmentation in plant populations: susceptible signals in plant traits and methodological approaches. Molecular Ecology 17: 5177–5188.
                  • 48. Carson HL (1990) Increased genetic variance after a population bottleneck. Trends in Ecology and Evolution 5: 228–230.
                  • 49. Young AG, Merriam HG, Warwick SI (1993) The effects of forest fragmentation on genetic variation in Acer Saccharum Marsh (Sugar Maple) Populations. Heredity 71: 277–289.
                  • 50. Hou Y, Lou A (2011) Population genetic diversity and structure of a naturally isolated plant species, Rhodiola dumulosa (Crassulaceae). PLoS ONE 6
                  • 51. Fu KT, Ohba H (2001) Crassulaceae. In: Wu ZY, Raven PH, editors. Beijing: Science Press. pp.202–268.
                    • 52. Mayuzumi S, Ohba H (2004) The phylogenetic position of eastern Asian Sedoideae (Crassulaceae) inferred from chloroplast and nuclear DNA sequences. Systematic Botany 29: 587–598.
                    • 53. Zhang JQ, Meng SY, Allen GA, Wen J, Rao GY (2014) Rapid radiation and dispersal out of the Qinghai-Tibetan Plateau of an alpine plant lineage Rhodiola (Crassulaceae). Molecular phylogenetics and evolution 77: 147–158.
                    • 54. Taberlet P, Gielly L, Pautou G, Bouvet J (1991) Universal primers for amplification of 3 noncoding regions of chloroplast DNA. Plant Molecular Biology 17: 1105–1109.
                    • 55. Hamilton MB (1999) Four primer pairs for the amplification of chloroplast intergenic regions with intraspecific variation. Molecular Ecology 8: 521–523.
                    • 56. Liu PL, Wan Q, Guo YP, Yang J, Rao GY (2012) Phylogeny of the genus Chrysanthemum L.: evidence from single-copy nuclear gene and chloroplast DNA sequences. PLoS ONE 7
                    • 57. Thompson JD, Higgins DG, Gibson TJ (1994) ClustalW - improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. Nucleic Acids Research 22: 4673–4680.
                    • 58. Zhang JQ, Meng SY, Wen J, Rao GY (2014) Phylogenetic relationships and character evolution of Rhodiola (Crassulaceae) based on nuclear ribosomal ITS and plastid trnL-F and psbA-trnH sequences. Systematic Botany 39: 441–451.
                    • 59. Librado P, Rozas J (2009) DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25: 1451–1452.
                    • 60. Swofford DL (2003) PAUP*. Phylogenetic analysis using parsimony (* and other methods). Version 4
                    • 61. Felsenstein J (1985) Confidence limits on phylogenies-an approach using the bootstrap. Evolution 39: 783–791.
                    • 62. Posada D (2006) ModelTest Server: a web-based tool for the statistical selection of models of nucleotide substitution online. Nucleic Acids Research 34: W700–W703.
                    • 63. Posada D, Buckley TR (2004) Model selection and model averaging in phylogenetics: advantages of akaike information criterion and Bayesian approaches over likelihood ratio tests. Systematic Biology 53: 793–808.
                    • 64. Ronquist F, Huelsenbeck JP (2003) MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19: 1572–1574.
                    • 65. Bandelt HJ, Forster P, Rohl A (1999) Median-joining networks for inferring intraspecific phylogenies. Molecular Biology and Evolution 16: 37–48.
                    • 66. Drummond AJ, Rambaut A (2007) BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evolutionary Biology 7
                    • 67. Richardson JE, Pennington RT, Pennington TD, Hollingsworth PM (2001) Rapid diversification of a species-rich genus of neotropical rain forest trees. Science 293: 2242–2245.
                    • 68. Zhang Q, Yang R, Wang Q, Liu JQ (2005) Phylogeography of Juniperus przewalskii (Cupressaceae) inferred from the chloroplast DNA trnT-trnF sequence variation. Acta Phytotaxonomica Sinica 43: 503–512.
                    • 69. Excoffier L, Lischer HEL (2010) Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Molecular Ecology Resources 10: 564–567.
                    • 70. Dupanloup I, Schneider S, Excoffier L (2002) A simulated annealing approach to define the genetic structure of populations. Molecular Ecology 11: 2571–2581.
                    • 71. Tajima F (1989) Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 123: 585–595.
                    • 72. Fu YX (1997) Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics 147: 915–925.
                    • 73. Schneider S, Excoffier L (1999) Estimation of past demographic parameters from the distribution of pairwise differences when the mutation rates very among sites: application to human mitochondrial DNA. Genetics 152: 1079–1089.
                    • 74. Rogers AR, Harpending H (1992) Population growth makes waves in the distribution of pairwise genetic differences. Molecular Biology and Evolution 9: 552–569.
                    • 75. Rogers AR (1995) Genetic evidence for a Pleistocene population explosion. Evolution 49: 608–615.
                    • 76. Excoffier L, Laval G, Schneider S (2005) Arlequin (version 3.0): an integrated software package for population genetics data analysis. Evolutionary Bioinformatics 1: 47–50.
                    • 77. Harpending HC (1994) Signature of ancient population growth in a low-resolution mitochondrial DNA mismatch distribution. Human Biology 66: 591–600.
                    • 78. Wolfe KH, Li WH, Sharp PM (1987) Rates of nucleotide substitution vary greatly among plant mitochondrial, chloroplast, and nuclear DNAs. Proceedings of the National Academy of Sciences of the United States of America 84: 9054–9058.
                    • 79. Yamane K, Yano K, Kawahara T (2006) Pattern and rate of indel evolution inferred from whole chloroplast intergenic regions in sugarcane, maize and rice. DNA research 13: 197–204.
                    • 80. Li JJ, Shi YF, Li BY (1995) Uplift of the Qinghai-Xizang (Tibet) Plateau and Global Change. Lanzhou.: Lanzhou University Press.
                      • 81. Guo ZT, Ruddiman WF, Hao QZ, Wu HB, Qiao YS, et al. (2002) Onset of Asian desertification by 22 Myr ago inferred from loess deposits in China. Nature 416: 159–163.
                      • 82. Harrison TM, Copeland P, Kidd WSF, Yin A (1992) Raising Tibet. Science 255: 1663–1670.
                      • 83. Shi YF, LI JJ, Li BY (1998) Uplift and Environmental Changes of Qinghai-Titetan Plateau in the Late Cenozoic. Guangzhou: Guangdong Science and Technology Press.
                        • 84. Spicer RA, Harris NBW, Widdowson M, Herman AB, Guo SX, et al. (2003) Constant elevation of southern Tibet over the past 15 million years. Nature 421: 622–624.
                        • 85. An ZS, Kutzbach JE, Prell WL, Porter SC (2001) Evolution of Asian monsoons and phased uplift of the Himalayan Tibetan plateau since Late Miocene times. Nature 411: 62–66.
                        • 86. Singh J, Singh S (1987) Forest vegetation of the Himalaya. The Botanical Review 53: 80–192.
                        • 87. Sun H (2002) Tethys retreat and Himalayan-Hengduanshan Mountains uplift and their significance on the origin and development of the Sino-Himalayan elements and alpine flora. Acta Botanica Yunnanica 24: 273–288.
                        • 88. Tang LY, Shen CM (1996) Late Cenozoic vegetational history and climatic characteristics of the Qinghai-Tibetan Plateau. Acta Micropalaeontologica Sinica 44: 872–877.
                        • 89. Liu J, Möller M, Provan J, Gao LM, Poudel RC, et al. (2013) Geological and ecological factors drive cryptic speciation of yews in a biodiversity hotspot. New Phytologist 199: 1093–1108.
                        • 90. Petit RJ, Aguinagalde I, de Beaulieu JL, Bittkau C, Brewer S, et al. (2003) Glacial refugia: hotspots but not melting pots of genetic diversity. Science 300: 1563–1565.
                        • 91. Gao Q, Zhang D, Duan Y, Zhang F, Li Y, et al. (2012) Intraspecific divergences of Rhodiola alsia (Crassulaceae) based on plastid DNA and internal transcribed spacer fragments. Botanical journal of the Linnean Society 168: 204–215.
                        • 92. Li Y, Zhai SN, Qiu YX, Guo YP, Ge XJ, et al. (2011) Glacial survival east and west of the 'Mekong-Salween Divide' in the Himalaya-Hengduan Mountains region as revealed by AFLPs and cpDNA sequence variation in Sinopodophyllum hexandrum (Berberidaceae). Molecular Phylogenetics and Evolution 59: 412–424.
                        • 93. Meng HH, Gao XY, Huang JF, Zhang ML (2014) Plant phylogeography in arid Northwest China: Retrospectives and perspectives. Journal of Systematics and Evolution DOI:https://doi.org/10.1111/jse.12088.
                        • 94. Hou Y, Lou A (2013) Phylogeographical patterns of an alpine plant, Rhodiola dumulosa (Crassulaceae), inferred from chloroplast DNA sequences. Journal of Heredity 105: 101–110.
                        • 95. Gao LM, Moeller M, Zhang XM, Hollingsworth ML, Liu J, et al. (2007) High variation and strong phylogeographic pattern among cpDNA haplotypes in Taxus wallichiana (Taxaceae) in China and North Vietnam. Molecular Ecology 16: 4684–4698.
                        • 96. Frei ES, Scheepens J, Stöcklin J (2012) High genetic differentiation in populations of the rare alpine plant species Campanula thyrsoides on a small mountain. Alpine Botany 122: 23–34.
                        • 97. Shimono Y, Watanabe M, Hirao AS, Wada N, Kudo G (2009) Morphological and genetic variations of Potentilla matsumurae (Rosaceae) between fellfield and snowbed populations. American Journal of Botany 96: 728–737.