Jian-Qiang Zhang
Da-Lv Zhong- College of Life Sciences, Shaanxi Normal University, Xi’an, China
How geological events and climate oscillations in the Pleistocene glaciation shaped the geographic distribution of genetic variation of species on the Qinghai-Tibetan Plateau (QTP) and its adjacent areas has been extensively studied. However, little studies have investigated whether closely related species in the same genus with similar physiological and life history traits responded similarly to the glacial climatic oscillations. If this is not the case, we would expect that the population histories of studied species were not driven by extrinsic environmental changes alone. Here we conducted a phylogeographic study of a succulent alpine plant Rhodiola fastigiata, using sequences from chloroplast genome and nrITS region, as well as ecological niche modeling. The results of R. fastigiata were compared to other congeneric species that have been studied, especially to R. alsia and R. crenulata. We found that for both markers, two geographic groups could be revealed, corresponding to the QTP plateau and the Hengduan Mountains, respectively, indicating isolated refugia in those two areas. The two groups diverged 1.23 Mya during the Pleistocene. We detected no significant population expansion by mismatch distribution analysis and Bayesian Skyline Plot. We found that even these similar species with similar physiological and life history traits have had different demographic histories in the Quaternary glacial periods. Our comparative phylogeographic study sheds new lights into phylogeographic research that extrinsic environmental changes are not the only factor that can drive population demography, and other factors, such as coevolved interactions between plants and their specialized pathogens, that probably played a role need to be examined with more case studies.
Introduction
How geological events and climate oscillations in the Pleistocene glaciation affected the distribution of genetic variation of species over time and space is a central issue in phylogeographic studies (Avise, 2009). As more and more species are surveyed phytogeographically, scientists are interested in whether organisms that inhabit the same community or have similar distributions experienced the same phylogeographic history, which is the focus of comparative phylogeography (Bermingham and Moritz, 1998; Avise et al., 2016). These studies utilized statistical phylogeography and ecological modeling to provide important insights into how congruence and incongruence happen (e.g., Barrow et al., 2017). Recent development includes statistically comparing different hypothetical models about the structure of refugia and post-glacial re-colonization routes, as well as spatially explicit phylogeographic analyses (Lemmon and Lemmon, 2008; Barrow et al., 2015). However, less studies have considered how closely related species in the same genus or lower taxonomic group with similar physiology and life histories respond to climatic oscillations. Other studies demonstrated that species with different characters (e.g., stress tolerance, life cycle, and dispersal ability) might have experienced different glacial histories (Shafer et al., 2010; Stewart et al., 2010). Nevertheless, if species with similar physiological or life-history traits showed different demographic histories, climate and intrinsic traits are not enough to explain demographic changes of species. Thus it is of great interest to explore how species in the same genus responded to the Quaternary climatic oscillations and if other factors affected species demographic changes.
Phylogeographic studies in the Qinghai-Tibetan Plateau (QTP; including the plateau platform, Himalaya and Hengduan Mountains) have increased in recent years (Liu J.Q. et al., 2012), as this area is an important alpine biodiversity hotspot (Myers et al., 2000; Mittermeier et al., 2004) with a very complex geological history and ecological heterogeneity. Recent studies have showed complex evolutionary scenarios for species in this area and different patterns in response to the geological events and climatic oscillations (Liu J.Q. et al., 2012). It is common for plants to survive in the refugium in Hengduan Mountain area (Yang et al., 2008; Zhang et al., 2010), while there might be some other small regional refugia on the plateau platform (Wang et al., 2009a,b, 2010; Opgenoorth et al., 2010). Two recent comparative studies showed that distribution of focal species remained more or less stable during the Quaternary (Luo et al., 2016, 2017), raising the question of what attributes of plants decided their reaction to glacial climatic oscillations.
Rhodiola, a genus mainly distributed in the QTP area, provides us an ideal model to study how much plants from the same genus differs in glacial history. There are several species in Rhodiola that have been studied phylogeographically: Rhodiola alsia (Gao et al., 2009), R. dumulosa (Hou and Lou, 2014), R. kirilowii (Zhang et al., 2014), R. crenulata (Zhang et al., Unpublished) and sect. Trifida (Li et al., 2018). These species represent different morphological space of Rhodiola: dioecious or hermaphrodite, having elongated rhizomes or not, bearing persistent old flowering stems or not. Here we added an example of R. fastigiata, which is morphologically very similar to R. alsia, to see if closely related species with similar distribution and morphology responded the same to the Quaternary climatic oscillations.
Rhodiola fastigiata (Hook. f. & Thomson) S.H. Fu is a perennial herb with a conspicuously elongated rhizome and dense persistent flower stems. It is distributed on rocky slopes at ca. 3500–5400 m in elevation of the QTP area (Fu and Ohba, 2001). The distribution pattern of the species is island-like, as the populations are isolated by forests at lower elevation. Here we used ITS and two plastid sequences, as well as ecological niche models (ENMs) to infer demographic history of this species. We also compared the results with other congeneric species to test the hypothesis that closely related species should have similar reactions to glacial climatic oscillations as they bear the same intrinsic traits and experienced the same extrinsic environmental changes.
Materials and Methods
Population Sampling
Through 2016–2017, we collected 22 populations of R. fastigiata in Yunnan, Sichuan, and Xizang provinces, covering all its distribution range. For each population, we sampled 7–20 individuals at least 20 m apart. The details of the collected samples are shown in Supplementary Table S1. Fresh leaves were directly put into silica gel for quick dehydration. We included a total of 287 individuals of R. fastigiata. R. coccinea were also collected for outgroup in the analyses.
DNA Extraction, PCR Amplification, Cloning and Sequencing
We used Plant Genomic DNA Kit (TianGen Biotech, Xi’an, China) to extract DNA from silica-gel dried leaves. ITS-1 and ITS-4 (Mayuzumi and Ohba, 2004) primers were used for amplifying ITS, c and f (Taberlet et al., 1991) for trnL-F, and trnS and trnG for trnS-G (Hamilton, 1999). Polymerase Chain Reaction (PCR) mixture were 20 μl with 2 μl 10× buffer, 0.5 μl of each primer, 0.4 μl of dNTP mixture, 1 U of Taq polymerase (TianGen Biotech, Xi’an, China) and 1 μl template genomic DNA. The PCR cycling programs followed Liu P.L.et al. (2012). The same amplification primers were used for sequencing of the chloroplast fragments and most ITS sequences. For individuals that have multiple peaks, we ligated the PCR products into pGEM-T Easy Vector using a Promega Kit (Promega Corporation, Madison, WI, United States). Then we chose plasmid that containing the right insertion to sequence. Sequencing was conducted on a 3730 automatic DNA sequencer in Tsingke Biotech, Xi’an, China. We used ContigExpress (a component of Vector NTI Suite 6.0, InforMax) to check and assemble contigs, and ClustalW v 1.7 (Thompson et al., 1994) to align the sequences. Finally, we manually checked the aligned matrix in BioEdit v 7.0.1.
Phylogeny, Divergence and Biogeography Analysis
As chloroplast genome evolves as a whole in plant, we concatenated the two plastid sequences. We extracted unique sequences (haplotypes and ribotypes) using DnaSP v5 (Librado and Rozas, 2009). We ignored indels as the two indels in trnS-G only occurred in one individual, respectively. To infer phylogenetic relationships among haplotypes and ribotypes, we implemented Maximum parsimony and Bayesian inference to reconstruct phylogeny with PAUP∗ 4.0b10 (Swofford, 2003) and MrBayes 3.2.1 (Ronquist and Huelsenbeck, 2003), respectively. One individual of R. coccinea was used as outgroup according to previous study (Zhang et al., 2014). For the MP tree, we set heuristic searches with MULTREES and TBR branch swapping. One thousand replicates of random addition sequences were used to construct staring trees. Bootstrap values (Felsenstein, 1985) were applied to evaluate support for nodes. One thousand replicates with ten replicates of random addition sequence and NNI branching swapping were applied. As Bayesian method needs an appropriate substitution model, we first determined it (TPM2uf + I + G for plastid dataset, and SYM + G for the ITS dataset) with Akaike information criterion (AIC) in Modeltest 3.7 (Posada and Buckley, 2004; Posada, 2006). We set four Metropolis-coupled Markov chain Monte Carlo chains with 10,000,000 generations to get the best convergence of the chain with a sampling sequence of 1000 generations. The average standard deviation of split frequencies was used to assess the convergence of the two runs. We discarded the first 20% as burn-in, and constructed a 50%-majority rule consensus tree with the remaining trees. We also constructed networks of haplotypes and ribotypes with NETWORK 4.2.0.1 (Bandelt et al., 1999) to detect the network connection of haplotypes.
Divergence time between lineages were inferred by BEAST software (Drummond and Rambaut, 2007). Before the BEAST run, we tested if the strict clock model was suitable using a likelihood ratio test in PAUP 4.0b10. As the results showed that molecular clock could not be rejected (2logeLR = 2.76, df = 18, P > 0.05) for the ITS data set, we set the inference parameter as GTR + G model and a strict molecular clock, a constant population size coalescent tree prior and a UPGMA starting tree. We set 20,000,000 MCMC generations with sampling frequency as every 2,000 generations, and the first 5,000,000 generations as burn-in. Tracer (Drummond and Rambaut, 2007) was used to examine convergence of chains. We did not find reliable fossils that could be used for dating analysis of R. fastigiata, so we conducted dating based on substitution rate for ITS in shrubs and herbs reported by previous study (6.075 ± 1.590 × 10-9 s/s/y; Richardson et al., 2001).
Population Genetics
Molecular diversity index including haplotype diversity (h) and nucleotide diversity (π) at both population and species level were calculated with DnaSP v5 (Librado and Rozas, 2009). We used PERMUT1 to estimate average gene diversity with a population (HS), total gene diversity (HT) and two population divergence index (GST, NST). We also performed 1,000 permutations replication to test the significant level. If NST value is significantly higher than GST, more close related haplotypes will cluster (Zhang et al., 2005), indicating the existence of phylogeographic structure. To infer if geographic distance is correlated with genetic distance (isolation by distance), we conducted a Mental test on matrices of pair-wise geographic distance and FST using ARLEQUIN v 3.5 (Excoffier and Lischer, 2010) with 1,000 random permutations.
We performed a spatial analysis of molecular variance (SAMOVA) with SAMOVA v 1.0 (Dupanloup et al., 2002). This program finds the best K-value (number of geographic groups) by maximizing FCT value between K groups of geographically adjacent populations. In our study, we set K from 2 to 10. After defining geographic groups, the amount of variation among populations within a group and a population were calculated by the hierarchical analysis of molecular variance (AMOVA) in ARLEQUIN v 3.5 (Excoffier and Lischer, 2010). A non-parametric permutation procedure with 1,000 permutations was conducted to test significant difference. Average FST between revealed geographic groups were also calculated with the same software.
To identify signatures of demographic expansion of populations and clades on the plastid tree, we estimated Tajima’s D and Fu’s Fs values (Tajima, 1989; Fu, 1997). Significantly negative D and FS values are expected if the focal population experienced expansion, because excessive rare and new mutation will appear in an expansion scenario. We also conducted a mismatch distribution analysis (Rogers and Harpending, 1992; Schneider and Excoffier, 1999) to detect the population expansion scenario of R. fastigiata. We pooled the whole haplotypes of each clade because evidence showed that population structure had little effect on mismatching distribution (Rogers, 1995). The fitting degree of observed mismatch distributions to the expected distribution under a recent expansion model (Rogers and Harpending, 1992; Excoffier et al., 2005) were tested using 1,000 parametric bootstrap replicates with sum of squared deviations (SSD) and the raggedness index (HRag) of Harpending (1994). If one group has the signal of expansion, we used the parameter-value (τ) to estimate when this expansion happened with equation t = τ/2u (Rogers and Harpending, 1992; Rogers, 1995). Here u = μ × k × g, where μ is the substitution rate (s/s/y), k is the average length of sequence data used, and g is the generation time (y). In this study, k was 1,593 bp, and the substitution rate was set to 2 × 10-9 s/s/y (Yamane et al., 2006). We assumed 10 years for generation time according to previous study (Zhang et al., 2014).
We also performed an Extended Bayesian Skyline Plot (EBSP) analysis using BEAST2 (Bouckaert et al., 2014) to estimate the demographic change of R. fastigiata. Effective population size through time were estimated based on coalescent process. As our test did not reject the molecular clock hypothesis, we used a strict clock model. We run the four MCMC chains for 50,000,000 generations with a sampling frequency of every 5,000 generations. After visualizing in Tracer (Rambaut et al., 2018), we drew the plot using a custom script in R software (R Core Team, 2017).
Ecological Niche Modeling
To better understand the potential range shift of R. fastigiata, we employed ENMs to model distribution change in response to glacial climatic oscillations. By combining our own data and online herbarium records (e.g., Chinese Virtual Herbarium, and Global Biodiversity Information Facility), we obtained a total of 149 spatially unique localities for the modeling. Data from online databases were checked by the author to exclude misidentification. MAXENT 3.3.3e (Phillips and Dudík, 2008) was used to model potential distribution area of R. fastigiata of two time points: the LGM and the present. We downloaded environmental layers of 19 bioclimatic variables (Supplementary Table S4) for the Last Glacial Maximum (LGM) and the current time from the WorldClim website at a spatial resolution of 2.5 arc-minutes (Hijmans et al., 2005). After excluding highly correlated climate variables by examining pairwise correlations, seven variables (Supplementary Table S4) with pairwise Pearson correlation coefficients below 0.7 were used. As other study stated, areas under the “receiver operating characteristic (ROC) curve” (AUC) (Peterson et al., 2008; Elith and Leathwick, 2009) values were used to evaluate the accuracy of each model prediction. The threshold for good performance was set to 0.7 (Fielding and Bell, 1997). Finally, we implemented DIVA-GIS 7.5 (Hijmans et al., 2001) to draw the suitable distributions ranges.
Results
Haplotype Variation and Distribution
After alignment, the total length of the combined cpDNA sequence was 1,593 bp. Thirty-six haplotypes based on 36 nucleotide substitutions were detected (including one outgroup individual; Supplementary Table S2). Unique sequences from each fragment were deposited in GenBank with accession number MH023238-MH023277. Haplotype diversity at the species level was h = 0.852, varying in different populations from 0.000 to 0.803 (SJL-4) (Supplementary Table S1). At the species level, nucleotide diversity was π = 0.0029, with a range from 0.000 to 0.0036 in different populations (Supplementary Table S1). Within-population gene diversity (HS) was significantly lower than total gene diversity (HT) (0.494 and 0.877, respectively, P < 0.05; Table 1). Twenty-three of the 36 haplotypes were found in only one population (Supplementary Table S2), and others occurred in at least two or more populations with H2 found in 13 of the 22 sampled populations (Figure 1). In the 22 populations, 5 only harbored a single haplotype, while SJL-1 harbored 6 haplotypes (Figure 1).
TABLE 1. Genetic diversity and genetic differentiation of 22 populations of Rhodiola fastigiata at the species level and group levels.
FIGURE 1. Map showing the sampling localities and the geographic distribution of haplotypes of Rhodiola fastigiata based on the cpDNA dataset. Pie charts show the proportion of haplotypes within each population. The numbers besides the circles represent population number as in Supplementary Table S1. Dashed line on the map indicates the distribution area of R. fastigiata. The two geographic groups defined by the SAMOVA analysis are also indicated by colored dashed lines.
Ribotype Variation and Distribution
The aligned ITS sequences was 640 bp in length. We detected 19 ribotypes, defined by 18 nucleotide substitutions (including one outgroup individual; Supplementary Table S3). Sequences were uploaded to GenBank with accession number MH 023219-MH023237. At the species level, ribotype diversity was 0.703, ranging from 0.000 to 0.833 (Supplementary Table S1). Nucleotide diversity π = 0.0018 at the species level, varying from 0.000 to 0.0027 in different populations. Population BDL had the highest ribotype diversity, while BM-2 had the highest π value. Similar to the plastid dataset, within-population gene diversity (HS) was significantly lower than total gene diversity (HT) (0.216 and 0.746, respectively; Table 1). Of the 19 ribotypes, 13 only occurred in one population (Supplementary Table S3). 14 out of the 22 populations only harbored one ribotypes, with population BDL had 6 ribotypes, the most in all populations (Figure 2).
FIGURE 2. Map showing the sites of sampled populations and the geographic distribution of ribotypes of R. fastigiata based on the ITS dataset. Pie charts show the proportion of ribotypes within each population. The numbers besides the circles represent population number as in Supplementary Table S1. Dashed line on the map indicates the distribution area of R. fastigiata. The two geographic groups defined by the SAMOVA analysis are also indicated by colored dashed lines.
Population Structure
The permutation tests of plastid data demonstrated that NST was significantly higher than GST (Table 1), indicating presence of phylogeographic structure. However, we failed to detect a significant phylogeographic structure from the ITS data. In the SAMOVA analysis for both cpDNA and ITS datasets FCT value reached a plateau when K = 3 (Supplementary Figure S1). In the plastid dataset, the first group consisted of 11 populations on the QTP plateau, and group two comprised 10 populations in the Hengduan Mountains. Group three only had one population on the plateau and was thus merged to group 1 in the discussion hereafter (Figure 1). For the ITS data, we also detected three geographic groups, two of them corresponding to the QTP plateau and the Hengduan Mountains, respectively, and the third group with only two ribotypes was merged to group 1 hereafter as they were also on the QTP plateau. Average FST between revealed geographic groups are shown in Table 2. AMOVA analysis showed that for the cpDNA dataset, 31.89% of the total variations occurred among groups decided by SAMOVA, and 44.77% of the total variations were with populations. For the ITS dataset, variations among groups also accounted for 71.26% of the total variation (Table 3). The Mantel test based on both datasets showed a significant correlation between genetic distance and geographic distance: a pattern of isolation-by-distance (cpDNA: r = 0.387, P < 0.05; ITS: r = 0.423, P < 0.01).
TABLE 2. Pairwise comparisons of FST among regions estimated from internal transcribed spacer (ITS) sequences (upper part) and cpDNA sequences (lower part) of Rhodiola fastigiata.
TABLE 3. Analysis of molecular variance (AMOVA) of cpDNA haplotypes and ITS ribotypes for Rhodiola fastigiata populations.
Phylogenetic Relationships and Lineages Divergence
We only showed the Bayesian tree as MP and Bayesian methods yield largely congruent topology for both datasets (Figures 3A, 4A). In the plastid data tree, all haplotypes of R. fastigiata formed a well-supported clade (Figure 3A). Two clades (Clades A and B) could be found. Clade A comprised 13 haplotypes and clade B only 3. The haplotype network (Figure 3B) depicted relationships between haplotypes more clearly: all haplotypes can be divided into two groups, corresponding to the two geographic groups determined by SAMOVA analysis (Although our SAMOVA analysis revealed a K-value of 3, group 3 only had one haplotype).
FIGURE 3. Phylogenetic relationships based on cpDNA haplotypes. (A) The Bayesian tree topology of the 36 cpDNA haplotypes detected in R. fastigiata. Numbers above the branches are MP bootstrap support values (left) and Bayesian posterior probability (right). Different branch color represents geographic groups defined by the SAMOVA analysis. H9 represents outgroup sequence of R. coccinea. (B) NETWORK-derived genealogic relationships. The sizes of the circles in the network were proportional to the observed frequencies of the haplotypes in each population. The white diamonds represent missing chlorotypes. Different color of circles represents geographic groups defined by the SAMOVA analysis.
FIGURE 4. Phylogenetic relationships based on ITS ribotypes. (A) The Bayesian tree of the 19 ITS ribotypes detected in R. fastigiata. Numbers above the branches are MP bootstrap support values (left) and Bayesian posterior probability (right). R1 represents outgroup sequence of R. coccinea. (B) NETWORK-derived genealogic relationships. The size of the circles in the network is proportional to the observed frequencies of the haplotypes in each population. The white diamonds represent missing ribotypes. Different color of circles represents geographic groups defined by the SAMOVA analysis.
For the ITS data, ribotypes were also divided into two geographic groups (indicated by different color in Figure 4). Similar to the plastid data, although our SAMOVA analysis revealed three geographic groups, group 3 only comprised two ribotypes that had the same distribution area with group 1. Dating analysis based on the ITS data showed that R. fastigiata diverged from its closest relatives at 3.47 Mya (Node 0; 95% HPD: 1.82–5.51 Mya) in the Pliocene (Figure 5). Further divergence occurred after 2.41 May (Node 1; 95% HPD: 1.20–3.72 Mya), in the Pleistocene.
FIGURE 5. Divergence time of R. fastigiata based on ITS ribotypes estimated with BEAST. Blue bars indicate 95% highest posterior density intervals. Numbers on the branch represent posterior probability of each node. The vertical red line on the tree tips indicates the Holocene.
Demographic Analyses
We did not detect any signal of population expansion of R. fastigiata at the species level or for any clades of the tree: under a population expansion model, all haplotypes as a whole showed a multi modal mismatch distribution in the mismatch distribution analysis (Supplementary Figure S2). Furthermore, SSD and the raggedness index significantly differed from expected values if an expansion hypothesis is true (Table 4). No significantly negative value of Tajima’s D and Fu’s Fs were also revealed (Table 4). Our Extended Bayesian Skyline Plot (EBSP) analysis revealed a similar pattern: the effective population size of R. fastigiata remained stable in last 0.82 million years, with a slight increasing toward the present time (Supplementary Figure S3).
TABLE 4. Results of the mismatch distribution analysis and neutrality tests of the three multiple-haplotype pDNA clades.
Species Distribution Models for R. fastigiata
The AUC value (AUC > 0.98; standard deviation [SD] < 0.01) showed that our modeling was effectively predictive according to 20 MAXENT replicating runs. The predicted distribution of R. fastigiata of present day is generally the same as the actual distribution today (Figure 6B). Our result also showed that the distribution of R. fastigiata during the LGM was slightly smaller than present distribution, especially in the eastern QTP and northern part of the Hengduan Mountains (Figure 6A). This pattern is in accordance with the EBSP results shown in the last section.
FIGURE 6. Distribution dynamics of R. fastigiata during the present day and the LGM based on species distribution modeling using MAXENT. Predicted distributions are shown for (A) the LGM model, (B) present. White squares represent presence locations (field survey data). LGM, last glacial maximum.
Discussion
Phylogeographic History of R. fastigiata
Rhodiola fastigiata inhabits the subnival belt a little below the permanent snowline, which represents the highest terrestrial habitat that an organism can occupy (Körner, 2003). As species growing in these habitats are cold adapted, we hypothesized that these species would expand their distribution area as temperature went down in the glacial times, and vice versa. Repeated expansion and contraction of distribution area will create an island-like distribution pattern. Our data was consistent with this hypothesis. In the AMOVA analysis, both markers showed a high FST value (ITS, 0.78; cpDNA, 0.55), indicating great difference among populations, and most variation was at group and population level (Table 3). In addition, NST was significantly greater than GST for the cpDNA data (Table 1), indicating significant phylogeographic structure. However, such phylogeographic structure was not detected in the ITS data. This results indicated that gene flow via pollen is more efficient than via seeds in R. fastigiata, which is congruent with its congeneric species R. kirilowii (Zhang et al., 2014). Another reason for lack of significant phylogeographic structure in the ITS data is that 63.6% of the sampled population harbored only one haplotype (Figure 2).
In the spatial analysis of molecular variance (SAMOVA), we detected two groups in both markers. The two groups largely corresponded to the Hengduan Mountains area and the QTP plateau platform. A high genetic differentiation among groups were detected by the AMOVA analysis (Table 3), and confirmed by our pairwise FST calculation among three groups (Table 2). Molecular dating analysis showed that these two groups diverged at 1.23 Mya (95% HPD: 0.56–2.01 Mya) based on the ITS dataset. Thus we can infer that the Pleistocene climatic oscillations might have triggered the divergence of the two groups. Populations of R. fastigiata might have retreated to isolated glacial refugia in the HM and QTP area during the glacial period. We detected no range expansion both at the species level and in each of the clades of the cpDNA tree (Figure 3 and Table 4). Ecological niche modeling results also support the scenario that later glacial time has little to do with the distribution of the studies species (Figure 6). This result is in congruence with two recent studies focused on several subnival species in the same area (Luo et al., 2016, 2017), in which ENM demonstrated that those species’ distribution area remained stable in the Quaternary glaciations.
Comparison to Other Species of Rhodiola
Gao et al. (2012) studied the Quaternary history of R. alsia, which is the first phylogeographic study in Rhodiola. R. fastigiata and R. alsia are very similar in morphology except tiny difference in leaf shape and length of carpel. Their work with two cpDNA markers (rpl20-rps12 and trnS-G) and ITS markers showed that there were probably three potential refugia, one in the Hengduan Mountain area and the other two on the QTP platform (Gao et al., 2009, 2012). The dated divergence time of three clades was 0.35–0.87 Mya ago, during the Pleistocene. The phylogeographic pattern of R. alsia is very similar to R. fastigiata with regard to refugia’s location. We also revealed a divergence between the QTP plateau group and the Hengduan Mountains group, indicating that refugia existed in both regions. However, the divergence of the two groups was 1.23 Mya (95% HPD: 0.56–2.01 Mya). Although there is overlap between the 95% HPDs, the intra-specific divergence time of R. fastigiata is slightly earlier than that of R. alsia. Besides, Gao et al. (2009) used mismatch distribution analysis and showed that population expansion had occurred in the history of R. alsia, whereas such pattern was not revealed for R. fastigiata in current study. A relatively stable population demography was also demonstrated in our ENM data (Figure 6) and EBSP analysis (Supplementary Figure S3). A recent study also showed that the Quaternary glacial had had little impact on four subnival plants on the QTP (Luo et al., 2016). Study of another morphological similar species, R. crenulata, which also inhabits subnival areas of the QTP, showed that there had been two glacial refugia, i.e., the HM and QTP platform (Zhang et al., Unpublished). However, mismatch distribution analysis showed that a recent population expansion happened in ca. 0.31 Mya. This expansion was also detected in the EBSP analysis.
Even though R. fastigiata and R. alsia have similar physiological and life history traits, they had different demographic histories in the Quaternary glacial periods, at least based on our analyses. Although previous literature has demonstrated that glacial histories of species with different features (e.g., cold tolerance, drought tolerance, life cycle, and dispersal ability) were also distinct (Shafer et al., 2010; Stewart et al., 2010), our results showed that even species with similar physiology traits in the same genus could have different response to glacial climatic oscillations. This indicates that factors other than climate would affect population histories of plants. One hypothesis is that the changing balance of coevolved interactions between hosts and their specialized pathogens could drive population dynamics (Ricklefs, 2015). Further studies are needed to test this hypothesis. However, this conclusion needs to be treated with caveats, as it should be tested using statistical phylogeography or spatially explicit phylogeographic analyses with balanced sampling and molecular markers (Lemmon and Lemmon, 2008; Barrow et al., 2015).
Refugia on the QTP was also identified for R. dumulosa (Hou and Lou, 2014), a cold adapted species with elongated rhizomes and persistent flowering stems. These traits may confer it the ability to survive in the refugia of the QTP platform. On the contrary, R. kirilowii, a species with no persistent flower stems, only occupied the Hengduan Mountain areas as major refugia. Refugia were also detected in other area of species’ distribution in Tianshan area in Xinjiang. Interestingly, a recent study of R. chrysanthemifolia, a woodland- and shrubbery-inhabited herbaceous species, used two cpDNA markers and ITS data and uncovered that there had been multiple microrefugia during the LGM, even earlier glaciations, a pattern also found in two sympatric tree and shrub species (Gao et al., 2016). These wood species might have provided it suitable habitats for surviving in the glacial periods. This pattern was further supported by a recent study including more species and populations (Li et al., 2018). Although the distribution area of R. chrysanthemifolia is similar to that of R. crenulata, R. fastigiata and R. alsia, their mirco-habitats are distinct: R. chrysanthemifolia is lower in elevation and mainly grows on the forest floor and in alpine shrubberies. This example showed that micro-inhabits play an important role in how plants respond to glacial climatic oscillations.
It is interesting to compare genetic diversity of hermaphrodite species and dioecious species, because Diocy is thought to evolve to promote outcrossing and increase heterozygosity (Barrett, 2002). We compared cpDNA and ITS diversity between hermaphrodite species and dioecious species, and no significant difference for both datasets was detected (ITS: p = 0.23; cpDNA: p = 0.36). Although this result needs to be treated with caution because of different sampling strategy in different studies, it provides evidence that dioecy in Rhodiola was not evolved to promote cross pollination. Therefore, the alternative hypothesis proposed by Darwin (1877) that plants evolve dioecy to allocate reproductive energy might be right. This makes sense because species like Rhodiola often grow on very high altitude and harsh environment where the growing season is very short.
Author Contributions
J-QZ conceived the ideas and wrote the manuscript. D-LZ, R-WZ, and W-YS conducted the experiments and collected the data. W-JS and D-LZ analyzed the data.
Funding
This study was supported by National Natural Science Foundation of China (No. 31500177), Shaanxi Science and Technology Program (No. 2016JQ3025), and the Fundamental Research Funds for the Central Universities (No. GK201603065 to J-QZ).
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 Dr. Yu-Qu Zhang for his support in sample collection and Dr. Gang Liu in the ENM analysis.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2018.00462/full#supplementary-material
FIGURE S1 | Correlation between the F statistics and grouping number (K = 2–10) from the SAMOVA results. (a) cpDNA; (b) ITS.
FIGURE S2 | Historical demography for overall populations and in each regional group based on the plastid DNA dataset. Clades A and B corresponded to the Bayesian phylogenetic tree in Figure 3. Mismatch distribution showing histogram of observed mismatch frequencies and best-fit curve of the sudden expansion model.
FIGURE S3 | Extended Bayesian Skyline Plot (EBSP) results calculated by BEAST2. Grey shade indicates 95% CPD of the effective population size.
TABLE S1 | Locations of populations of Rhodiola sect. Trifida 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.
TABLE S2 | Haplotype composition of 22 sampled based on cpDNA dataset for populations.
TABLE S3 | Ribotype composition of 21 sampled based on ITS dataset for populations.
TABLE S4 | Bioclimatic variables (named BIO1 to BIO19) from WorldClim (Hijmans et al., 2005). Variables marked with an asterisk (∗) were used for the climatic niche models of Rhodiola fastigiata.
Footnote
References
Avise, J. C. (2009). Phylogeography: retrospect and prospect. J. Biogeogr. 36, 3–15. doi: 10.1111/j.1365-2699.2008.02032.x
Avise, J. C., Bowen, B. W., and Ayala, F. J. (2016). In the light of evolution X: comparative phylogeography. Proc. Natl. Acad. Sci. U.S.A. 113, 7957–7961. doi: 10.1073/pnas.1604338113
Bandelt, H. J., Forster, P., and Rohl, A. (1999). Median-joining networks for inferring intraspecific phylogenies. Mol. Biol. Evol. 16, 37–48. doi: 10.1093/oxfordjournals.molbev.a026036
Barrett, S. C. (2002). The evolution of plant sexual diversity. Nat. Rev. Genet. 3, 274–284. doi: 10.1038/nrg776
Barrow, L. N., Bigelow, A. T., Phillips, C. A., and Lemmon, E. M. (2015). Phylogeographic inference using Bayesian model comparison across a fragmented chorus frog species complex. Mol. Ecol. 24, 4739–4758. doi: 10.1111/mec.13343
Barrow, L. N., Soto-Centeno, J. A., Warwick, A. R., Lemmon, A. R., and Moriarty Lemmon, E. (2017). Evaluating hypotheses of expansion from refugia through comparative phylogeography of south-eastern Coastal Plain amphibians. J. Biogeogr. 44, 2692–2705. doi: 10.1111/jbi.13069
Bermingham, E., and Moritz, C. (1998). Comparative phylogeography: concepts and applications. Mol. Ecol. 7, 367–369. doi: 10.1046/j.1365-294x.1998.00424.x
Bouckaert, R., Heled, J., Kuhnert, D., Vaughan, T., Wu, C. H., Xie, D., et al. (2014). BEAST 2: a software platform for Bayesian evolutionary analysis. PLoS Comput. Biol. 10:e1003537. doi: 10.1371/journal.pcbi.1003537
Darwin, C. R. (1877). The Different Forms of Flowers on Plants of the Same Species. London: Murray. doi: 10.5962/bhl.title.110054
Drummond, A. J., and Rambaut, A. (2007). BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol. Biol. 7:214. doi: 10.1186/1471-2148-7-214
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
Elith, J., and Leathwick, J. R. (2009). Species distribution models: ecological explanation and prediction across space and time. Annu. Rev. Ecol. Evol. Syst. 40, 677–697. doi: 10.1146/annurev.ecolsys.110308.120159
Excoffier, L., Laval, G., and Schneider, S. (2005). Arlequin (version 3.0): an integrated software package for population genetics data analysis. Evol. Bioinform. 1, 47–50. doi: 10.1177/117693430500100003
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
Felsenstein, J. (1985). Confidence-limits on phylogenies - an approach using the bootstrap. Evolution 39, 783–791. doi: 10.1111/j.1558-5646.1985.tb00420.x
Fielding, A. H., and Bell, J. F. (1997). A review of methods for the assessment of prediction errors in conservation presence/absence models. Environ. Conserv. 24, 38–49. doi: 10.1017/S0376892997000088
Fu, K. T., and Ohba, H. (2001). “Crassulaceae,” in Flora of China, Vol. 8, eds C. Y. Wu and P. H. Raven (Beijing: Science Press), 202–268.
Fu, Y. X. (1997). Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics 147, 915–925.
Gao, Q. B., Zhang, D. J., Chen, S. Y., Duan, Y. Z., Zhang, F. Q., Li, Y. H., et al. (2009). Chloroplast DNA phylogeography of Rhodiola alsia (Crassulaceae) in the Qinghai-Tibet Plateau. Botany 87, 1077–1088. doi: 10.1139/B09-059
Gao, Q. B., Zhang, D. J., Duan, Y. Z., Zhang, F. Q., Li, Y. H., Fu, P. C., et al. (2012). Intraspecific divergences of Rhodiola alsia (Crassulaceae) based on plastid DNA and internal transcribed spacer fragments. Bot. J. Linn. Soc. 168, 204–215. doi: 10.1111/j.1095-8339.2011.01193.x
Gao, Q. B., Zhang, F. Q., Xing, R., Gornall, R. J., Fu, P. C., Li, Y., et al. (2016). Phylogeographic study revealed microrefugia for an endemic species on the Qinghai-Tibetan Plateau: Rhodiola chrysanthemifolia (Crassulaceae). Plant Syst. Evol. 302, 1179–1193. doi: 10.1007/s00606-016-1324-4
Hamilton, M. B. (1999). Four primer pairs for the amplification of chloroplast intergenic regions with intraspecific variation. Mol. Ecol. 8, 521–523.
Harpending, H. C. (1994). Signature of ancient population growth in a low resolution mitochondrial DNA mismatch distribution. Hum. Biol. 66, 591–600.
Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G., and Jarvis, A. (2005). Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965–1978. doi: 10.1002/joc.1276
Hijmans, R. J., Guarino, L., Cruz, M., and Rojas, E. (2001). Computer tools for spatial analysis of plant genetic resources data: 1, DIVA-GIS. Plant Genet. Resour. Newslett. 127, 15–19.
Hou, Y., and Lou, A. (2014). Phylogeographical patterns of an alpine plant, Rhodiola dumulosa (Crassulaceae), inferred from chloroplast DNA sequences. J. Hered. 105, 101–110. doi: 10.1093/jhered/est072
Körner, C. (2003). Alpine Plant Life: Functional Plant Ecology of High Mountain Ecosystems. Berlin: Springer. doi: 10.1007/978-3-642-18970-8
Lemmon, A. R., and Lemmon, E. M. (2008). A likelihood framework for esti- mating phylogeographic history on a continuous landscape. Syst. Biol. 57, 544–561. doi: 10.1080/10635150802304761
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. Phylogen. Evol. 121, 110–120. doi: 10.1016/j.ympev.2018.01.001
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
Liu, J. Q., Sun, Y. S., Ge, X. J., Gao, L. M., and Qiu, Y. X. (2012). Phylogeographic studies of plants in China: advances in the past and directions in the future. J. Syst. Evol. 50, 267–275. doi: 10.1111/j.1759-6831.2012.00214.x
Liu, P. L., Wan, Q., Guo, Y. P., Yang, J., and Rao, G. Y. (2012). Phylogeny of the genus Chrysanthemum L.: evidence from single-copy nuclear gene and chloroplast DNA sequences. PLoS One 7:e48970. doi: 10.1371/journal.pone.0048970
Luo, D., Xu, B., Li, Z. M., and Sun, H. (2017). The ‘Ward Line–Mekong–Salween Divide’is an important floristic boundary between the eastern Himalaya and Hengduan Mountains: evidence from the phylogeographical structure of subnival herbs Marmoritis complanatum (Lamiaceae). Bot. J. Linn. Soc. 185, 482–496. doi: 10.1093/botlinnean/box067
Luo, D., Yue, J. P., Sun, W. G., Xu, B., Li, Z. M., Comes, H. P., et al. (2016). Evolutionary history of the subnival flora of the Himalaya-Hengduan Mountains: first insights from comparative phylogeography of four perennial herbs. J. Biogeogr. 43, 31–43. doi: 10.1111/jbi.12610
Mayuzumi, S., and Ohba, H. (2004). The phylogenetic position of eastern Asian Sedoideae (Crassulaceae) inferred from chloroplast and nuclear DNA sequences. Syst. Bot. 29, 587–598. doi: 10.1600/0363644041744329
Mittermeier, R. Z., Gil, P. R., Hoffman, M., Pilgrim, J., Brooks, T., Mittermeier, C. G., et al. (2004). Hotspots revisited: Earth’s Biologically Richest and Most Endangered Terrestrial Ecoregions. Mexico City: CEMEX.
Myers, N., Mittermeier, R. A., Mittermeier, C. G., Da Fonseca, G. A., and Kent, J. (2000). Biodiversity hotspots for conservation priorities. Nature 403, 853–858. doi: 10.1038/35002501
Opgenoorth, L., Vendramin, G. G., Mao, K., Miehe, G., Miehe, S., Liu, J. L., et al. (2010). Tree endurance on the Tibetan Plateau marks the world’s highest known tree line of the Last Glacial Maximum. New Phytol. 186, 780–780. doi: 10.1111/j.1469-8137.2009.03007.x
Peterson, A. T., Papeʂ, M., and Soberón, J. (2008). Rethinking receiver operating characteristic analysis applications in ecological niche modeling. Ecol. Model. 213, 63–72. doi: 10.1016/j.ecolmodel.2007.11.008
Phillips, S. J., and Dudík, M. (2008). Modeling of species distributions with Maxent: new extensions and a comprehensive evaluation. Ecography 31, 161–175. doi: 10.1111/j.0906-7590.2008.5203.x
Posada, D. (2006). ModelTest Server: a web-based tool for the statistical selection of models of nucleotide substitution online. Nucleic Acids Res. 34, W700–W703. doi: 10.1093/nar/gkl042
Posada, D., and Buckley, T. R. (2004). Model selection and model averaging in phylogenetics: advantages of akaike information criterion and Bayesian approaches over likelihood ratio tests. Syst. Biol. 53, 793–808. doi: 10.1080/10635150490522304
R Core Team (2017). R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing.
Rambaut, A., Drummond, A. J., Xie, D., Baele, G., and Suchard, M. A. (2018). Tracer v1.7. Available at: http://tree.bio.ed.ac.uk/software/tracer
Richardson, J. E., Pennington, R. T., Pennington, T. D., and Hollingsworth, P. M. (2001). Rapid diversification of a species-rich genus of neotropical rain forest trees. Science 293, 2242–2245. doi: 10.1126/science.1061421
Ricklefs, R. E. (2015). Intrinsic dynamics of the regional community. Ecol. Lett. 18, 497–503. doi: 10.1111/ele.12431
Rogers, A. R. (1995). Genetic evidence for a Pleistocene population explosion. Evolution 49, 608–615. doi: 10.1111/j.1558-5646.1995.tb02297.x
Rogers, A. R., and Harpending, H. (1992). Population growth makes waves in the distribution of pairwise genetic differences. Mol. Biol. Evol. 9, 552–569.
Ronquist, F., and Huelsenbeck, J. P. (2003). MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19, 1572–1574. doi: 10.1093/bioinformatics/btg180
Schneider, S., and 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.
Shafer, A., Cullingham, C. I., Cote, S. D., and Coltman, D. W. (2010). Of glaciers and refugia: a decade of study sheds new light on the phylogeography of northwestern North America. Mol. Ecol. 19, 4589–4621. doi: 10.1111/j.1365-294X.2010.04828.x
Stewart, J. R., Lister, A. M., Barnes, I., and Dalén, L. (2010). Refugia revisited: individualistic responses of species in space and time. Proc. Biol. Sci. 277, 661–671. doi: 10.1098/rspb.2009.1272
Swofford, D. L. (2003). PAUP∗. Phylogenetic Analysis Using Parsimony (∗ and Other Methods). Version 4. Sunderland, MA: Sinauer Associates,
Taberlet, P., Gielly, L., Pautou, G., and Bouvet, J. (1991). Universal primers for amplification of 3 noncoding regions of chloroplast DNA. Plant Mol. Biol. 17, 1105–1109. doi: 10.1007/BF00037152
Tajima, F. (1989). Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 123, 585–595.
Thompson, J. D., Higgins, D. G., and Gibson, T. J. (1994). Clustal-W - improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. Nucleic Acids Res. 22, 4673–4680. doi: 10.1093/nar/22.22.4673
Wang, H., Qiong, L., Sun, K., Lu, F., Wang, Y., Song, Z., et al. (2010). Phylogeographic structure of Hippophae tibetana (Elaeagnaceae) highlights the highest microrefugia and the rapid uplift of the Qinghai-Tibetan Plateau. Mol. Ecol. 19, 2964–2979. doi: 10.1111/j.1365-294X.2010.04729.x
Wang, L. Y., Abbott, R. J., Zheng, W., Chen, P., Wang, Y. J., and Liu, J. Q. (2009a). History and evolution of alpine plants endemic to the Qinghai-Tibetan Plateau: Aconitum gymnandrum (Ranunculaceae). Mol. Ecol. 18, 709–721. doi: 10.1111/j.1365-294X.2008.04055.x
Wang, L. Y., Ikeda, H., Liu, T. L., Wang, Y. J., and Liu, J. Q. (2009b). Repeated range expansion and glacial endurance of Potentilla glabra (Rosaceae) in the Qinghai-Tibetan Plateau. J. Integr. Plant Biol. 51, 698–706. doi: 10.1111/j.1744-7909.2009.00818.x
Yamane, K., Yano, K., and Kawahara, T. (2006). Pattern and rate of indel evolution inferred from whole chloroplast intergenic regions in sugarcane, maize and rice. DNA Res. 13, 197–204. doi: 10.1093/dnares/dsl012
Yang, F. S., Li, Y. F., Ding, X., and Wang, X. Q. (2008). Extensive population expansion of Pedicularis longiflora (Orobanchaceae) on the Qinghai-Tibetan Plateau and its correlation with the Quaternary climate change. Mol. Ecol. 17, 5135–5145. doi: 10.1111/j.1365-294X.2008.03976.x
Zhang, J. Q., Meng, S. Y., and Rao, G. Y. (2014). Phylogeography of Rhodiola kirilowii (Crassulaceae): a story of Miocene divergence and Quaternary expansion. PLoS One 9:e112923. doi: 10.1371/journal.pone.0112923
Zhang, Q., Yang, R., Wang, Q., and Liu, J. Q. (2005). Phylogeography of Juniperus przewalskii (Cupressaceae) inferred from the chloroplast DNA trnT-trnF sequence variation. Acta Phytotaxon. Sin. 43, 503–512. doi: 10.1360/aps040148
Keywords: Crassulaceae, Hengduan Mountains, Qinghai-Tibetan Plateau, Quaternary climatic oscillations, Rhodiola fastigiata
Citation: Zhang J-Q, Zhong D-L, Song W-J, Zhu R-W and Sun W-Y (2018) Climate Is Not All: Evidence From Phylogeography of Rhodiola fastigiata (Crassulaceae) and Comparison to Its Closest Relatives. Front. Plant Sci. 9:462. doi: 10.3389/fpls.2018.00462
Received: 25 January 2018; Accepted: 23 March 2018;
Published: 10 April 2018.
Edited by:
Genlou Sun, Saint Mary’s University, CanadaReviewed by:
Daniel Pinero, Universidad Nacional Autónoma de México, MexicoYingjuan Su, Sun Yat-sen University, China
Copyright © 2018 Zhang, Zhong, Song, Zhu and Sun. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jian-Qiang Zhang, jqzhang@snnu.edu.cn; zjqsbzh@gmail.com