Li Feng1
Zhong-Hu Li
Gui-Fang Zhao- 1Key Laboratory of Resource Biology and Biotechnology in Western China (Ministry of Education), College of Life Sciences, Northwest University, Xi'an, China
- 2College of Life Sciences, Shaanxi Normal University, Xi'an, China
The East Himalaya-Hengduan Mountains (EH-HM) region has a high biodiversity and harbors numerous endemic alpine plants. This is probably the result of combined orographic and climate oscillations occurring since late Tertiary. Here, we determined the genetic structure and evolutionary history of alpine oak species (including Quercus spinosa, Quercus aquifolioides, and Quercus rehderiana) using both cytoplasmic-nuclear markers and ecological niche models (ENMs), and elucidated the impacts of climate oscillations and environmental heterogeneity on their population demography. Our results indicate there were mixed genetic structure and asymmetric contemporary gene flow within them. The ENMs revealed a similar demographic history for the three species expanded their ranges from the last interglacial (LIG) to the last glacial maximum (LGM), which was consistent with effective population sizes changes. Effects of genetic drift and fragmentation of habitats were responsible for the high differentiation and the lack of phylogeographic structure. Our results support that geological and climatic factors since Miocene triggered the differentiation, evolutionary origin and range shifts of the three oak species in the studied area and also emphasize that a multidisciplinary approach combining molecular markers, ENMs and population genetics can yield deep insights into diversification and evolutionary dynamics of species.
Introduction
Historical processes such as geographic and climate changes have profoundly shaped the population genetic structure and demographic history of extant species (Hewitt, 2000). Global cyclical cooling-warming events during the Quaternary Period have resulted in periodic expansions and contractions of the ranges of species. For instance, in western Eurasia, species retreated into southern refuges during glacial episodes and recolonized during warmer interglacials; this periodical process created a frequently observed south-north gradient in genetic diversity (southern richness, northern purity) during postglacial re-expansion (Hewitt, 1996). Other researches have suggested that environmental heterogeneity may also contribute to the genetic differentiation among species or populations (e.g., Ortego et al., 2012, 2015; Guichoux et al., 2013). In China, phylogeographic studies have focus mainly on the roles of historical orogenesis, climatic oscillations and environmental heterogeneity in evolutionary history of biotas in Qinghai-Tibet Plateau (QTP) and adjacent regions (Qiu et al., 2011; Liu et al., 2012; Wen et al., 2014). However, recent study has suggested that the demographic history of flora occurred in the QTP must be revised because the QTP has been 4–5 km high since the mid-Eocene while numerous studies linked young inferred divergence times to recent QTP uplift phases (also see review by Renner (2016) and references therein).
The QTP and adjacent regions have long been considered as biodiversity hotspots with numerous endemic species (Myers et al., 2000). The uplift of QTP as a consequence of the collision of the Indian subcontinent with the Eurasian plate 40–50 million years ago (Yin and Harrison, 2000) and subsequent Quaternary climatic oscillations altered the distribution patterns of organisms as well as affected biodiversity there (Qiu et al., 2011). It is hypothesized that uplifts of the QTP, which started during middle Eocene and continued into middle Pliocene, created conditions for the origin and divergence of multitudinous species in the QTP and adjacent regions (Qiu et al., 2011). During the Quaternary glacial periods, these locations likely became refuges for plants belonging to the North Temperate Zone of East Asia (Wu, 1987; Qiu et al., 2011; Liu et al., 2012). During interglacials, species previously retreated into a refuge might recolonize high-altitude regions and also areas outside the QTP (Qiu et al., 2011, and references therein). Nonetheless, some species in these areas responded to Quaternary climate changes in other ways. For instance, Taxus wallichiana and Picea likiangensis supposed to originate in the HHM regions experienced ranges expansion during the Last Glacial Maximum (LGM) rather than during the Last Interglacial (LIG) (Li et al., 2013; Liu et al., 2013). It has yet unknown if this unexpected pattern of range shift is valid also for other cold-tolerant trees in this area. In addition, previous studies suggested that the accumulated genetic differences between species due to genetic admixture and introgression might be eroded by a shift of distribution ranges introduced by geological and climatic changes (Du et al., 2011; Ma et al., 2014), however, it remains poorly understood.
The genus Quercus comprises 531 species of trees and shrubs, and includes keystone taxa in the temperate and (sub-)tropical areas of the Northern Hemisphere (Nixon, 1993). Hybridization between related oak species is common and frequently genetic admixture has been reported (Aldrich and Cavender-Bares, 2011, and references therein). Weak reproductive isolation, frequent interspecies gene exchange and high phenotypic variation result in diffuse species boundaries within the genus causing much debate about species concepts (Burger, 1975; Rushton, 1993). Up to date, numerous studies in Europe and America have focused on the impact of climate change on their evolutionary history (e.g., Ortego et al., 2012, 2015; Guichoux et al., 2013; Cavender-Bares et al., 2015). Within our scope, few reports concerned oak species in China (Zeng et al., 2010, 2011, 2015), and even less alpine sclerophyllous oaks.
In the present study, we focused on species of the Quercus group ilex (synonyms Quercus subgenus Heterobalanus), i.e., Quercus spinosa, Quercus aquifolioides, and Quercus rehderiana (see details in Note S1). Previous studies based on molecular and morphological evidence of pollen suggested this group corresponded to the subgenus Heterobalanus (Manos et al., 2001; Denk and Grimm, 2010; Denk and Tekleva, 2014). However, no information is available on the evolutionary history and demography of them. In this study we use an integrative approach to determine the evolutionary history of the three oak species to infer the most plausible scenario(s) of speciation and their responses to climatic changes. The major objectives of this study were to investigate (i) whether the uplift of the QTP as well as Quaternary climate shifts gave rise to the current differences of the three oak species; (ii) if so, whether they have responded similarly to climatic changes and the uplift of the QTP; and (iii) whether and how the gene flow effected their population genetic structure. Obviously, knowledge of the population structure and evolutionary history of the evergreen oak species in temperate China is needed to understand the complicated evolutionary history of species in the East Himalaya-Hengduan Mountains (EH-HM) and adjacent regions.
Materials and Methods
Sampling and Genotyping
Samples were collected strictly according to morphological description from 608 individual trees in 33 natural populations and one individual of Q. spinosa in Taiwan (Table S1), covering the majority of their distribution ranges. Generally, 4–20 samples were collected from each population depending on the population size. All sampled individuals from each population were spaced at least 100 m apart. After DNA extraction, three chloroplast DNA (cpDNA), one nuclear ribosomal internal transcribed spacer (ITS) and 12 nuclear simple sequence repeat (nSSR) microsatellite loci were amplified (see details in Note S2, and Table S10).
DNA Sequence Analysis
Population genetic parameters, including the number of segregating sites (S), the number of haplotypes (NH), the haplotype diversity (Hd), the nucleotide diversity (ϕ), and neutrality test statistics of Tajima's D (Tajima, 1989) and Fu's (Fs) (Fu and Li, 1993) to assess possible expansions with associated significance values for species were estimated by DnaSP v5.00.04 (Librado and Rozas, 2009). The program PERMUT v1.2.1 (Pons and Petit, 1996) was employed to estimate the average gene diversity within populations (HS), total gene diversity (HT) and the differentiation for unordered alleles (GST) and for ordered alleles (NST) based on 1000 random permutations.
The phylogenetic relationships of cpDNA haplotypes were reconstructed via Bayesian and maximum likelihood (ML) methods with Castanea mollissima as the outgroup, and we utilized the program BEAST v1.7.5 (Drummond et al., 2012) to estimate their divergence times. In tree-related analysis, we employed jMODELTEST v1.0 (Posada, 2008) to evaluate the best-fit model of nucleotide substitution for maximum likelyhood (ML) method to infer haplotype relationships(GTR+I+G in the present case). The program MrBayes v3.1 (Ronquist and Huelsenbeck, 2003) was used for the Bayesian analysis with a burn-in of 3 × 106 Markov chain Monte Carlo (MCMC) repetitions for cpDNA sequences. We also used maximum parsimony (MP) and ML in PAUP* v4.0 beta 10 (Swofford, 2002) to determine the relationships of cpDNA haplotypes and visualized results using FIGTREE v1.3.1 (http://tree.bio.ed.ac.uk/software/figtree/).
The program BEAST v1.7.5 was used to estimate the divergence times of the three related species based on an uncorrelated lognormal relaxed molecular clock model when considering the difference of substitution rates among three cpDNA fragments and the different evolutionary rate along tree branch during their evolutionary history. Three independent runs of 5 × 108 MCMC steps were carried out, with sampling at every 5000 generations, following a burn-in of the initial 10% cycles. MCMC samples were inspected in TRACER v1.5 to confirm sampling adequacy and convergence of the chains to a stationary distribution. The results were visualized using FIGTREE v.1.3.1. Generally, previous studies utilized the cpDNA substitutions rates [1.1 × 10−9 to 2.9 × 10−9 substitutions per site per year (s/s/y)] (Wolfe et al., 1987) when estimating divergence time, while, Petit and Hampe (2006) suggested that woody plants shared a slow evolutionary clock due to a long generation time. Hence, we estimated the genetic divergence times between the clades by using cpDNA substitution rates of 7.905 × 10−10 s/s/y in Q. variabilis (Chen et al., 2012).
We also used the locus-by-locus analysis of molecular variance (AMOVA) approach as implemented in ARLEQUIN v3.5 (Excoffier and Lischer, 2010) to examine variance within and between population groups with significance differences on 1000 permutations. The program NETWORK v4.2.0.1 (Bandelt et al., 1999) was used to infer the relationships of haplotypes of cpDNA and rDNA with median-joining networks.
Microsatellite Data Analysis
We used MICROCHECKER v2.2.3 (Van Oosterhout et al., 2004) to test the presence of null alleles for all loci. For each microsatellite locus, genetic diversity indices (total number of alleles, AO, observed heterozygosity, HO, expected heterozygosity over all populations, HE, gene diversity within populations, HS, and total genetic diversity, HT) were estimated by POPGENE v1.31 (Yeh et al., 1997). Linkage disequilibrium (LD) and departure from Hardy-Weinberg equilibrium (HWE) were also evaluated using FSTAT v2.9.3 (Goudet, 2001). Significance levels were corrected by the sequential Bonferroni method (Rice, 1989).
Genetic differentiation among populations was evaluated using θ (FST) (Weir and Cockerham, 1984) and the standardized genetic differentiation G'ST(Hedrick, 2005) across loci with the program SMOGD (Crawford, 2010). Compared to traditional measures, G'ST is a more suitable measure for highly polymorphic markers such as microsatellites. In addition, we performed AMOVA analysis in ARLEQUIN v3.5 for these three related species, respectively. The significance of fixation indices was tested using 104 permutations.
The program STRUCTURE v2.3.3 (Pritchard et al., 2000) was used to determine whether these three named species are genetically distinct based on the Bayesian clustering analysis using the microsatellite data. The analyses were run using the admixture model with independent allele frequencies. A total of 10 independent simulations were run for K from 1 to 30 with 1.0 × 105 burn-in steps followed by 2.0 × 105 MCMC steps. Two alternative methods were utilized to estimate the most likely number (K) of genetic clusters with the program STRUCTURE HARVSTER (Earl and vonHoldt, 2012), i.e., by tracing the change in the average of log-likelihood L(K) (Pritchard et al., 2000) and by calculating delta K (ΔK) (Evanno et al., 2005). As with the whole dataset, two clusters were identified (Figure 3) and we repeated the structure analysis independently except for running K from 1 to 15 for populations occurred in the EH-HM region (green cluster in Figure 3C).
In order to detect the relationship between genetic and geographic distances among populations within species, we conducted a Mantel test (Mantel, 1967) with 10,000 permutations to determine the possible role of isolation-by-distance (IBD) in the formation of current population structure using IBDW v3.23 (Jensen et al., 2005).
In order to determine whether there has been subsequent population expansion within the species, we used BOTTLENECK v1.2.02 (Piry et al., 1999) to detect genetic bottlenecks in all populations. First of all, we evaluated any excess in heterozygosity in a population at mutation-drift equilibrium through Wilcoxon sign-rank test with 104 replications. Then, we performed the analysis under stepwise mutation model (SMM) and two-phase model of microsatellite evolution with the proportion of single step stepwise mutation set to 95% and the variance set to 5% with 104 simulations as recommended by Piry et al. (1999). The program 2MOD v0.2 (Ciofi et al., 1999) was used to evaluate the relative likelihood of migration-drift equilibrium, i.e., the relative importance of gene flow vs. genetic drift in forming current population structure. We ran the program six times with each model for all three species; 105 iterations were performed, and the first 10% of the iterations were discarded as burn-in. We used Bayesian factors as calculated by TRACER v1.5 to describe the probability of the most likely model.
To explore the historical and contemporary gene flow within the three oak species, we employed MIGRATE-N v3.6 (Beerli, 2006) and BAYESASS v3.0 (Wilson and Rannala, 2003) to estimate the effective number of migrants (2Nem, where Ne = effective population size, m = migration rates per generation) and contemporary migration rates (over the last few generations, mc), respectively. We defined four clusters as follows: one for Q. aquifolioides, one for Q. rehderiana, and the other two for Q. spinosa as suggested by the results of STRUCTURE (see results section). We performed maximum-likelihood analyses in MIGRATE-N v3.6 using 10 short chains (104 trees) and three long chains (105 trees); 104 trees were discarded as burn-in. To increase the efficiency of the MCMC, we applied the following parameters: replicated = YES, long chains, randomtree = YES, heating = ADAPTIVE: 1 (1, 1.2, 1.5, 3.0). In order to ensure the consistency of estimates, we repeated this procedure five times and reported average maximum-likelihood estimates with 95% confidence intervals. When estimated the contemporary gene flow using BAYESASS v3.0, We examined the parameters including migration rates (ΔM), allele frequencies (ΔA) and inbreeding coefficients (ΔF) to ensure the optimal acceptance rates of the three parameters fell within the range of 20–60%. The delta values for the three parameters were 0.26, 0.37, and 0.49, respectively. Following that, we conducted the analyses with 107 iterations after a burn-in of 106 iterations, setting 1000 as the sampling frequency. Ten independent runs were executed to minimize the convergence problem. The result with the lowest deviance was adopted according to the method of Meirmans (Meirmans, 2014). We estimated the 95% credible interval as mc ± 1.96 • standard deviation.
In order to evaluate the plausible scenario of evolutionary history (i.e., divergent scenario) and population dynamics (i.e., changes in effective population sizes) of these oak species, We used the approximate Bayesian computation (ABC) procedure in DIYABC v2.04 (Cornuet et al., 2014). Firstly, we used a multispecies coalescent model to estimate the species tree based on data of three cpDNA fragments under the Star BEAST (*BEAST) option implemented in BEAST v1.7.5. The result of species tree suggested the Q. spinosa was more ancient than other two species. Furthermore, Ma (2006) suggested that the Q. rehderiana was the youngest species within the three oak species. Additionally, Zhou (1992, 1993) suggested that Chinese oak species originated from the Indo-Chinese region, which implies the migration routs of oak species in this study are west to east. In other words, the oak species or populations occurred in the EH-HM region are more ancient than that in other regions. Hence, based on the above results and the STRUCTURE results (Figure 3C and Figure S3), nine alternative scenarios of population history for the four lineages were summarized (Figure 5A, Tables S6, S7). We assumed uniform priors on all parameters and used a goodness-of-fit test to check the priors of all parameters before implementing the simulations (Table S7). Each simulation was summarized by the following summary statistics: the mean number of alleles and the mean genic diversity for each lineage, the FST, the mean classification index, and shared allele distance between pairs of lineages. To select the model that best explains the evolutionary history observed in these species, 1 million simulated data sets for each scenario were run and total 9 million simulations were run for all scenarios, of which the 1% closest to the observed data was used to estimate the relative posterior probabilities of each scenario via a logistic regression approach and posterior parameter distributions according to the most likely scenario (Cornuet et al., 2008, 2014). In addition, DIYABC was used to simulate and examine demographic changes in the recent past. We tested the following six scenarios of demographic changes of the oak species in this study based on the discordant results of neutrality test of cpDNA and ecological niche modeling (see Results Section): continuous shrinkage or expansion; shrinkage-expansion or expansion-shrinkage; and shrinkage-expansion-shrinkage or expansion-shrinkage-expansion (Figure 4B and Table S7). DIYABC allows selection of the demographic scenario that best fits the data and parameters of interest. Summary statistics included the mean number of alleles across loci, the mean gene diversity across loci, and the mean allele size range variation across loci.
Ecological Niche Modeling
Ecological niche modeling was employed to predict suitable palaeo- and current distribution ranges of the three species via MAXENT v.3.3.3k (Phillips et al., 2006) and DESKTOP GARP v1.1.6 (available at http://www.nhm.ku.edu/desktopgarp/). All the bioclimatic variables (19 environmental variables) were downloaded from WorldClim website with a resolution of 2.5 arc-min (http://www.worldclim.org/) (Hijmans et al., 2005). The LGM data used in this study were under Community Climate System Model (CCSM; Collins et al., 2006). Because of the environmental layers were of 2.5 arc-min spatial resolution for the “present-day” and LGM periods, the available 30 arc-sec spatial resolution layers of LIG were resampled using ARC GIS 10.0 to obtain same level of resolution. In addition, to reduce the effects of overfitting, ecological niche models (ENMs) based on current climatic data were trained using the methods of Sheppard (2013). Subsequently, only seven variables (BIO4, BIO7, BIO9, BIO11, BIO15, BIO18, BIO19) were used to develop the current distributions for the three oak species (Table S8), while the remaining 12 were discarded because of high autocorrelation. Information on the geographic distribution of related species was based on the 33 populations included in this study, and sites were removed that were separated from each other by <0.1°, so as to reduce the effect of spatial autocorrelation, after that we totally obtained 170 records of the three oak species (50 for Q. aquifolioides, 40 for Q. rehderiana and 80 for Q. spinosa) from Chinese Virtual Herbarium (http://www.cvh.org.cn/cms/). We used the default parameters of MAXENT and included 80% of species records for training and 20% for testing the model and the default convergence threshold (10−5) and set the number of maximum iterations to 5000 and the number of replicates to 10. The output format was set to logistic. Models no better than having random predictive ability have an area under the receiver operating characteristic (ROC) curve (AUC) value of 0.5, whereas models with perfect predictive ability have an AUC approximating 1. An AUC score above 0.7 was considered as good model performance (Fielding and Bell, 1997). For GARP, we changed the total run to 100, maximum iterations to 1000 and training proportion to 80. Training data were divided into true training data (for model rule development) and intrinsic testing data (for rule testing intrinsic to GARP processing) (Figure S4). The “best subset” of the total data was generated according to the softomission threshold of 20% and commission threshold of 50%. We predicted the occurrence of each species based on present data and a uniform probability distribution. Finally we used DIVA-GIS v7.5 (Hijmans et al., 2001) and ARC GIS v10.0 to draw the range of suitable distributions.
Furthermore, to assess levels of niche divergence within the three oak species, we used ENMTOOLS v1.3 (Warren et al., 2008) to calculated the niche overlap statistic Schoener's D (Schoener, 1968) and standardized Hellinger distance (calculated as I). We then used the background similarity test as implemented in ENMTOOLS v1.3 to determine whether the ENMs of different pairs of species are more or less similar than expected from the differences in the environmental backgrounds of the regions where they occur (Warren et al., 2008). The background area should include accessible areas for the three oak species, not just the observed niche or an area tightly delimited by species occurrence (McCormack et al., 2010). The test was run 100 times, and the observed niche overlap was compared to the distribution of overlap values from the runs in each direction to determine whether species' niches were significantly more or less divergent than expected (two-tailed test). Rejection of the null hypothesis demonstrates that observed niche differences are a function of habitat selection and/or availability and may be interpreted as evidence for niche conservatism or niche divergence. Failure to reject the null hypothesis indicates that niche differentiation is explained by variation in the environmental conditions available to each within their respective ranges (Warren et al., 2008, 2010; McCormack et al., 2010).
Results
Cp DNA and Its Sequence Polymorphism
The concatenated sequences of the three cpDNA fragments were 1785 bp in length, a total of 25 chlorotypes were identified. Q. aquifolioides (C1, C3, and C5) and Q. rehderiana (C23–C25) had three private haplotypes, respectively, while two haplotypes (C2 and C4) were shared between them. Q. spinosa possessed 15 private haplotypes (C7–C19, C21–C22); in addition, it shared C6 with Q. aquifolioides and C20 with Q. rehderiana, respectively (Figure 1 and Table S2). Estimates of average within-species haplotype diversity was highest in Q. spinosa (Hd = 0.927). Values of NST and GST were both very high in all three species, and in no instance was NST significantly higher than GST (Table 1).
Figure 1. (A) Network of the chloroplast (cp) DNA haplotypes detected in three oak species. Different species are denoted by different colors of the circle, each sector of a circle is in proportional to the frequency of each chlorotype. (B) Geographic distribution of the chloroplast (cp) DNA haplotypes detected in three oak species. Haplotype frequencies of each population are denoted by the pie charts with population IDs in the circle. Green, blue and red dots indicate the sample locations of Quercus spinosa, Q. aquifolioides, and Q. rehderiana, respectively.
Table 1. Genetic diversity parameters based on cpDNA and ITS sequences in the population groups of the three species.
We obtained a 501-bp fragment of ITS regions, and a total of 14 ribotypes were identified. Among these genotypes, N3 and N5 occurred across the three species, N1 and N4 were unique to Q. aquifolioides, while Q. rehderiana and Q. spinosa possessed two (N13 and N14) and five private ribotypes (N8-N12), respectively (Figure S1A and Table S3). Values of Hd and pi were similar for the three related species. Estimates for average within-population and total ITS diversity were higher in Q. aquifolioides (HS = 0.306, HT = 0.883) than that in the other two species. In Q. aquifolioides and Q. spinosa, NST (0.848 and 0.958, respectively) was significantly larger than GST (0.654 and 0.872, respectively) (p < 0.05, Table 1), implying the presence of phylogeographic structure.
The AMOVA analyses for cpDNA and ITS data showed that genetic variation mainly occurred among populations within species (77.42 and 76.56%, respectively, Table 2). For nuclear and chloroplastic loci, the Tajima's D values were positive in all cases except for the ITS data in Q. spinosa; all these values were not statistically significant (Table 1).
Table 2. Analysis of molecular variance (AMOVA) for groups of related species.
The structure obtained from NETWORK for cpDNA was basically consistent with the topology of the phylogenetic tree (Figures 1A, 2). For cpDNA data, C20, C12, and C25 took central positions, and might be more ancient. Furthermore, almost all populations harbored one private haplotype. When referred to the ITS subsets (Figure S1A), N5 and N7 took central positions in the network with a higher proportion of individuals (16.54 and 9.02%, respectively). N3 was found to be the predominant ribotypes, and was mainly distributed from the central China outward to the adjacent EH-HM regions with the highest proportion of individuals (32.33%).
Figure 2. BEAST-derived chronograms of three oak species based on cpDNA (psbA-trnH, psbB-psbF, and matK) sequences. Blue bars indicate the 95% highest posterior density (HPD) credibility intervals for node ages (in Myr ago, Ma). Bootstrap values (>50%) based on ML, MP, and BI analysis and posterior probabilities are labeled above and below nodes, respectively. Mean divergence dates and 95% HPDs for major nodes (a–g) are summarized in the upper left figure.
Molecular Dating Based on cpDNA Data
The BEAST-derived chronograms of oak species based on three cpDNA sequences suggested that the divergence could date back to 11.45 million years ago (Ma) (95% highest posterior density, HPD: 5.72–23.85 Ma, PP = 1.00), indicating a Mid-Miocene split between the species and the outgroup. And the coalescent time of all chlorotypes within species (c. 8.54 Ma, 95% HPD: 4.46–14.93 Ma, PP = 1.00) indicated that the oak species diverged in Late Miocene (Figure 2). In addition, strong differentiation occurred within taxon during the periods of Quaternary climatic oscillations.
Nuclear Microsatellite Diversity and Population Structure
The null allele test indicated a lower frequency of null allele at each of the 12 loci than the threshold frequency (□ = 0.15) across the 33 populations, and there was no evidence for LD. After Bonferroni correction, significant deviation from HWE induced by homozygote excess was detected in all 12 loci when all samples were treated as a single population within species. However, there were no HWE deviations within each population after Bonferroni correction.
Genetic diversity indices were summarized for each loci and species (Table 3 and Table S4 in). Amplification of all 12 loci in 609 individuals revealed a total of 211 alleles with a range of 7–33 alleles per locus, and the mean HT and HS were 0.699 and 0.673 with a range of 0.279–0.905 and 0.278–0.891, respectively. Population differentiation was significant for all 12 loci (p < 0.05, Table S4) with the mean FST value for multilocus estimate being 0.371 (range: 0.120–0.543). As for the standardized genetic differentiation, due to the high HT and variable nSSR makers used in this study, G'ST was higher than FST across all loci (Table S4), which was in line with previous studies [Meirmans and Hedrick, 2010; also see review by Edelaar and Björklund (2011)]. Within species, mean AS and HE were similar, ranging from 12.5 to 13.8 and from 0.642 to 0.694, respectively. On average, Q. rehderiana had the highest FIS value whereas Q. aquifolioides had the highest HO value.
Table 3. Diversity measures in the related species.
The nSSR-derived AMOVA demonstrated significant genetic differentiation (RST = 0.3628, p < 0.01), and the majority (63.72%) of the variation partitioned within populations (Table 2). There was no obvious IBD for any of the three species (Figure S2). For the Bayesian analysis of population structure, the most likely number of genetic clusters was estimated at 2 (Figure 3). Populations of Q. spinosa, including one population (DL) from Yunnan Province, those from the east of Mts. Hengduan (except GY and QL) and one population of Q. rehderiana (WN), formed one group; the other group comprised all populations of Q. aquifolioides, six populations of Q. rehderiana (except WN), and the remaining populations of Q. spinosa from Mts. Hengduan and Yungui Platea. (Figure S3A). Populations occurred in the EH-HM region were further subdivided into three sub-clusters, which was in line with the delimitation of three oak species occurred there (Figure S3), however, the low ΔK value implied, indeed, the genetic structure of these three sub-clusters had some similarities.
Figure 3. Bayesian inference analysis of microsatellite data for determining the most likely number of clusters (K) for the three species. The distribution of the likelihood L(K) values (A) and ΔK values (B) are presented for K = 1–30 (10 replicates). STRUCTURE plots (C) are presented for best K = 2 and K = 3, respectively (1 = SHS; 2 = LS; 3 = SQS).
The Wilcoxon test under the SMM and TPM models failed to detect any recent genetic bottleneck in any of the populations, and all populations showed an L-sharp allele distribution when ruling out two populations (Table S5, SHS and LS with smaller population sizes). According to the 2MOD analysis, the most likely model of population history which led to the current population structure was gene flow-drift model rather than pure drift model (P = 1.0, Bayesian factor = 100,000).
Gene Flow within the Three Oak Species
Our study revealed interesting patterns of historical and contemporary gene flow among the four clusters. Almost all the historical gene flow of the related pairs were symmetrical with slight differences, however, there were significant asymmetrical contemporary gene flow in related pairs except for that within Q. spinosa (Table 4). It should be pointed out that we used the program IMa2 (Hey, 2010) to estimate the gene flow between related pairs in the beginning, however, we couldn't obtain ideal results due to its disequilibrium results of parameter settings in our study and it's slow computational speed.
Table 4. Rates of historical and contemporary gene flows per generation among the four clusters as estimated using microsatellite data with the programs MIGRATE and BAYESASS, respectively.
Evolutionary Dynamics and Changes in Effective Population Size
In the DIYABC analysis, the posterior probabilities for scenarios 2 was 0.5539 [95% confidence interval (CI) 0.5139–0.5940], much higher than other eight scenarios. The median values of the effective population sizes of QS1, QS2, QA, QR and NA were 3.75 × 105, 1.19 × 104, 2.57 × 105, 5.22 × 105, and 1.68 × 105, respectively. The estimated median time of divergence between QS1, QS2, and QR (t1), and the divergence time between QA and the ancestor of QS1 and QS2 (t2) were 3.24 × 104 and 7.58 × 104 generations ago, respectively. Assuming the generation time to be 150 years, the divergence times of t1 and t2 corresponded to 4.86 Ma and 11.37 Ma, respectively. The estimated median mutation rates and the proportion of multiple step mutations in the generalized stepwise model of microsatellites were estimated to be 1.75 × 10−5 and 0.564, respectively (Table 5).
Table 5. Posterior median estimate and 95% highest posterior density interval (HPDI) for demographic parameters in scenarios 2 based on the nuclear multilocus microsatellite data for three oak species.
Our nuclear microsatellite data for the three closely related oak species were best fitted with the contraction-expansion model [posterior probability = 0.8114 (95% CI: 0.7933–0.8295), 0.7247 (95% CI: 0.7093–0.8401) and 0.6583 (95% CI: 0.6440–0.6725) for the Q. spinosa, Q. aquifolioides and Q. rehderiana, respectively]. As indicated by the contraction-expansion model, the point estimates of the current population sizes for the Q. spinosa, Q. aquifolioides and Q. rehderiana were 5.72 × 105, 5.43 × 105, and 5.72 × 105, respectively, the population sizes between current and ancient were 7.41 × 102, 2.00 × 103, and 1.56 × 103, respectively, and the ancient population sizes were 5.31 × 105, 4.84 × 105, and 4.60 × 105, respectively. The contraction time was 3.17 × 104, 3.22 × 104, and 3.21 × 104 generations ago for the Q. spinosa, Q. aquifolioides and Q. rehderiana, respectively (c. 4.80, 4.83, and 4.82 Ma if a generation length of 150 years is assumed), and the expansion time was 3.03 × 103, 3.67 × 103, and 3.53 × 103 generations ago for the Q. spinosa, Q. aquifolioides, and Q. rehderiana, respectively (c. 0.45, 0.55, and 0.53 Ma if a generation length of 150 years is assumed) (Table S9).
Ecological Niche Modeling
There were similar change tendencies of suitable distributions of the three oak species obtained by MAXENT and GARP (Figure 5A and Figure S4). All models had high predictive ability (AUC > 0.9). In addition, the projected present-day distribution is consistent with collection records (Figure 4A). The ENMs suggested the three oak species continued their expansions during the LGM. Since the LGM, however, except for the Q. spinosa has distinct expansion, there have been little change in the predicted spatial distributions of the other two species despite climate warming (Figure 4A).
Figure 4. (A) The nine scenarios of the population history of three oak species in DIYABC V2.04. QS1 and QS2 represent two groups of Quercus spinosa; QA and QR represent Q. aquifolioides and Q. rehderiana, respectively. The current population sizes of QS1, QS2, QA, and QR were denoted as NQS1, NQS2, NQA, and NQR, respectively, while NA represents the population size of the ancestral lineage at time t3 (or t2). t1–t3, divergence times for the depicted event. (B) Schematic representation of six demographic models of changes in population size tested within three oak species (Quercus spinosa, Q. aquifolioides, and Q. rehderiana) using DIY-ABC. NA, ancestral population size; N1, current population size; N1b and N2, populations sizes between NA and N1; db, duration of bottleneck.
For background tests, the observed niche overlap values for I and D were significantly higher (p < 0.001) than the predicted scores under the null hypothesis both for Q. aquifolioides vs. Q. spinosa, Q. spinosa vs. Q. rehderiana, and Q. rehderiana vs. Q. aquifolioides (Figure 4B), indicating that the niches of all three species were very different and more conserved than expected based on the environments available for occupation by each. Additionally, the values of I and D greater than 0.5 in all three pairs of species, which indicated that the niches of all the species pairs were slightly overlapped (Figure 5B).
Figure 5. (A) Predicted distributions of Quercus aquifolioides, Q. spinosa, and Q. rehderiana based on ecological niche modeling using MAXENT, black dots indicate extant occurrence points. (B) The background tests between different pairs of species. Null distributions are shown by dotted blue bars for D and solid red bars for I. The x-axis indicates values of I and D, and y-axis indicates number of randomizations. Red and blue arrows indicate values in actual MAXENT runs, respectively.
Discussion
High Genetic Diversity and Strong Genetic Differentiation
The three species investigated in our study are characterized by long generation times, wind pollination, and seed dispersal through animals and gravity; all these characteristics are associated with high genetic diversity and differentiation (Hamrick et al., 1979, 1992). Our study did observe high genetic diversities for all three species based on both cpDNA and ITS datasets (HT = 0.971–1 and 0.807–0.883, respectively), and the estimated differentiation values (GST) were larger than the mean genetic differentiation calculated for maternally inherited plastid markers in white oaks of Europe (0.975–1 vs. 0.848) (Petit et al., 2002), implying the existence of strong barriers to gene flow. The high mountains, gorges, and complex paleo-drainage basins in EH-HM region have served as effective barriers to seed and/or pollen dispersal and impacted genetic diversification as well as differentiation of organisms, which has also been reported in other phylogeographic studies (e.g., Li et al., 2013; Liu et al., 2013). In addition, the low ability of seeds is also an important factor contributing to high population differentiation. Most of the Fagus seeds simply drop to the ground near their parent trees and a few may roll down on steep terrain, sometimes animals (e.g., jays or squirrels) may also carry some seeds short distances (Gómez, 2003; Xiao et al., 2009). This may also be the case for the three oak species, although their seed dispersal modes are still to be studied in detail.
Our results of AMOVA (Table 2) demonstrated that both cpDNA and ITS marker systems are strongly differentiated among populations but almost homogenous within populations. When combined with the general lack of IBD, The island-like population genetic structure of three oak species in accordance with their highly fragmented habitats located in the EH-HM region, and largely driven by historical effects of genetic drift rather than currently limited gene flow alone (Hutchison and Templeton, 1999). In fact, similar genetic structure have also been detected in other species in this region, such as Garrulax elliotii from the eastern Himalayas (Qu et al., 2011) and four herbs from HHM region (Luo et al., 2015).
Demographic History and Ranges Expansion of the Three Oak Species
The uplift of the QTP since the Cenozoic had a significant impact on species differentiation and genetic structure (Qiu et al., 2011). Some studies have shown that the main uplifting of QTP occurred in 10 Ma (Mulch and Chamberlain, 2006), reaching peak elevation already shortly before the Late Pliocene (Sun et al., 2011). Our molecular dating by cpDNA (Figure 2) and SSR data both fell into the late Miocene. The split between the related oak species and C. mollissima (Figure 2) might due to the enhanced monsoon climate have triggered differentiation between these species after the uplifting of the QTP at 15–13 Ma (An et al., 2001). Initial differentiation within taxa at 8.54 Ma (95% HPD: 4.46–14.93 Ma) may be explained by expanding arid areas in Asian inland, the effect of East Asian monsoon before 9–8 Ma, and the tempestuous uplift of the QTP during this period (An et al., 2001). This may have resulted in lots of small fragmented habitats with different microclimate, thereby impacting the direction of natural selection (Sobel et al., 2010). Furthermore, our ABC simulation indicated the divergence time of Q. rehderiana was 4.86 Ma (95%HPD: 3.36–5.91 Ma), implying a Pliocene split. Likewise, molecular dating of cpDNA data revealed dramatic divergence since 5.08 Ma (95%HPD: 2.29–9.44 Ma, Figure 2). Our divergence time estimation coincides with a period of intense uplift of the Hengduan Mountain massif (Sun et al., 2011; Favre et al., 2015), suggesting that the split could similarly have occurred during the colonization of the newly available terrain by these oak species. In fact, an association between Pliocene uplift of the QTP and time of intraspecific and/or interspecific divergence has previously been noted in this region (e.g., Li et al., 2013; Sun et al., 2014). It is plausible, that species had fragmented habitats caused by the QTP uplifting during Pliocene may have fostered both intraspecific and interspecific divergence on a large scale in the region.
However, the divergence time presented here should be treated with caution because we utilized the substitution rates of Q. variabilis, which is a deciduous tree and has shorter generation time than the three species, might be underestimate the divergence time. Another nuanced caveat here is some inherent problems of SSR data, for instance, homoplasy, trending to underestimate divergence time on large time scales (Selkoe and Toonen, 2006). Moreover, the assumption of no gene flow in DIYABC also leads to the underestimation of the divergence time between species (Leaché et al., 2013). However, many Miocene fossils belonging to group Ilex were discovered in the QTP could not be classified as any extant species (Zhou, 1993), it seems likely that the divergence time of the three oak species fell into the Miocene.
Our ABC demographic analysis detected a contraction-expansion model for all the three oak species, their effective population sizes decreased within the past c. 4.80 Ma and increased around c. 0.5 Ma. Although the EH-HM region became colder during the Pliocene due to the extensive uplift of the QTP (Shi et al., 1999), which facilitated the growing of three oak species. However, other alpine species (e.g., spruce, Sun et al., 2014) in this region could also adapt this climate, occupying some of the niches originally belonged to oaks in the present study and enlarged their ranges via successful competition. Based on the process above, the three closely oak species might decrease their population sizes. Although the temperature on the QTP increased at the end of the largest glaciation (ca. 1.2–0.6 Ma), the relatively cold climate may have continued until the late Ionian stage (0.3–0.126 Ma) (Shi et al., 1998). Thus, because these sclerophyllous oaks could adapt cold habitats (Zhou et al., 1994), it is feasible that they expanded their ranges and increased their population sizes during this period. The fact that a moderately cold climate has prevailed on the EH-HM region as the LIG will have provided opportunities for these three species to have continued their ranges expansion. Indeed, our ENM analyses also suggested that all the three species had experienced ranges expansion from LIG (0.14–0.12 Ma) to LGM, which was similar to other alpine plants reported in this region, e.g., P. likiangensis (Li et al., 2013) and T. wallichiana (Liu et al., 2013). It seems therefore, that trees such as mentioned above, which occur in cold environments nowadays, may have exhibited different range dynamics during past climatic oscillations relative to species associated with warmer environments.
Genetic Admixture and Asymmetric Contemporary Gene Flow
Although these three oak species were morphologically differentiated, they were genetically admixed in the present study due to incomplete lineage sorting and interspecific hybridization. Our results indicated that the haplotypes N2, N3, and N5 shared among species are ancient and have disjunctive distribution (Figure S1B), and shared haplotypes C2, C4, C6, C20, and N6 between sympatric species (Figure 1 and Figure S1B) as expected for incomplete lineage sorting (Maddison and Knowles, 2006) and introgressive hybridization, respectively.
Generally, numerous studies have revealed shared DNA polymorphisms between closely related species or species complex (e.g., Du et al., 2009; Li et al., 2012; Simeone et al., 2016). This situation can be divided into two categories: firstly, retention of ancestral polymorphisms which caused by incomplete lineage sorting during and following speciation (Heckman et al., 2007; Wilyard et al., 2009); secondly, introgression or introgressive hybridization which caused by genetic exchange after secondary contact between previously geographically separated species (Liston et al., 1999; Gay et al., 2007). Incomplete lineage sorting occurs when lineages fail to coalesce in the ancestral population of species (Pamilo and Nei, 1988; Maddison, 1997; Degnan and Salter, 2005). Therefore, the probability of incomplete lineage sorting depends on both the effective population size (Ne) in the ancestral population of species, which determines the rate of coalescence of lineages, and the time between successive speciation events (Hudson, 1990; Degnan and Salter, 2005). While the introgression is mostly detected in zones of sympatry/parapatry between two or more species, in other words, it occurs in co-distributed populations belong to different species [Ortego et al., 2015; also see review by Abbott et al. (2013) and references therein].
Our results of 2MOD strongly supported the gene flow-drift model, indicating that interspecific hybridization significantly contributed to the observed genetic admixture between the three species. Despite these three oaks have weak boundaries and close phylogenetic relationships (Figures 1A, 2), the possibility that the shared ancestral polymorphism can't be completely rejected. In this regard, the real cause of genetic admixture within these oak species deserves further study.
Our Bayesian cluster analysis (Figure 3C) indicated that populations from EH-HM region could be considered as a single cluster disregarding the species boundaries. One plausible explanation is that, despite existing genetic barriers in the EH-HM region, due to changes of their distribution ranges that give rise to sympatric distribution as well as existing incomplete reproductive isolation within oak species, gene flow (Table 4) among them will blur species boundaries. Similar situations were reported in other studies (Du et al., 2011; Ma et al., 2014). Populations in the east of Mts. Qinling (Figure S3B) could be another cluster due to analogous environments resulting in phenotypic convergence or fixed similar alleles during population expansions (Excoffier and Ray, 2008).
Our results indicated there were asymmetric contemporary gene flow for five related pairs whereas the historical gene flow for all pairs seemed symmetric (Table 4). A possible explanation for the asymmetric contemporary gene flow is the species richness in the EH-HM region. Generally, at the present time, Q. spinosa occurs in low-altitude regions, while Q. aquifolioides grows in high-altitude regions; both are common keystone species in local forests, however, Q. rehderiana occupies moderate-altitude areas with scattered distribution, such distribution patterns usually lead to stronger gene flow from core populations to peripheral ones as suggested by Alleaume-Benharira et al. (2006). In addition, human activities in mid-/low-altitude regions during the past decades would more significantly affect the population sizes of Q. rehderiana, giving rise to the current model of gene flow.
In general, gene flow is a vital factor for population structure over time, it may reduce local adaptation by homogenizing populations found in differing environments or by spreading detrimental alleles across populations, and it also serves to introduce potentially adaptive alleles to populations, and increases genetic diversity (Sexton et al., 2011; Epps and Keyghobadi, 2015; Welt et al., 2015). Previous study have suggested that contemporary and historical gene flow were generally low and similar among populations occurred in highly fragmented habitats (Chiucchi and Gibbs, 2010), and if genetic structure was weak or admixture, it might due to strong historical gene flow (Epps et al., 2013). Indeed, our results (Table 4) indicated the historical and contemporary gene flow among these oak species were generally low, however, whether the genetic admixture is due to historical gene flow rather than contemporary ones, we need further investigation because oak species themselves have weak boundaries.
Conclusions
Our analyses of the three related species reveal similar demographic history according to ENMs and ABC analyses. The neutral genetic markers did not depict the profile of differentiation entirely due to introgression or incomplete lineage sorting. Geological and climatic changes since Miocene have markedly affected the observed patterns of genetic variation and divergence among and within species. Our study indicates that a combination of orographic and bioclimatic analyses can yield deep insights into the diversification and evolutionary history of species in the EH-HM and adjacent regions.
Author Contributions
GZ and ZL designed the research; LF and JY collected samples; LF, QZ, and YZ performed experiments; LF and ZQ analyzed data; LF and ZQ prepared figures and tables; LF, ZQ, and ZL wrote the manuscript; and LF, QZ, ZQ, JY, YZ, ZL, and GZ revised the manuscript.
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 acknowledge Hans-Peter Comes and Peterman B Peter for their constructive comments on earlier versions of the manuscript. This research was financially supported by the Program for Changjiang Scholars and Innovative Research Team in University (No. IRT1174) and the Key Laboratory Program Founded by the Education Department of Shaanxi Province (12JS083). The authors are grateful to P. F. Dai, Z. H. Huang, X. Y. Wang, and T. Zhou for their assistance in sample collections, and to Z. L. Liu, P. Zhao, and D. Duan for their valuable suggestions on data analyses.
Supplementary Material
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2016.01688/full#supplementary-material
Figure S1. (A) Network of the ITS haplotypes detected in three oak species. Different species are denoted by different colors of the circle, each sector of a circle is in proportional to the frequency of each haplotype. (B) Geographic distribution of the ITS haplotypes detected in three oak species. Haplotype frequencies of each population are denoted by the pie charts with population IDs in the circle. Green, blue and red dots indicate the sample locations of Quercus spinosa, Q. aquifolioides, and Q. rehderiana, respectively.
Figure S2. Mantel tests between microsatellite-based genetic distances and geographic distances. (A) Quercus spinosa; (B) Q. aquifolioides; (C) Q. rehderiana.
Figure S3. Bayesian inference analysis of microsatellite data for determining the most likely number of clusters (K) for the three lineages of oak species occurred in the EH-HM region. The distribution of the likelihood L(K) values (A) and ΔK values (B) are presented for K = 1–15 (10 replicates). STRUCTURE plots (C) are presented for best K = 3 (QS1: populations of the Quercus spinosa with green in Figure 3C; QA: Q. aquifolioides; QR: Q. rehderiana).
Figure S4. Ecological niche models predicted distributions of the three species using GARP during three periods. Different colors corresponded to different fitting indices with low in blue and high in red for the current, LGM and LIG distribution.
Table S1. Sampling details and haplotypes detected in each population of the three related species.
Table S2. The variable sites in aligned cpDNA sequences of psbA-trnH, psbB-psbF, and matK that yielded 25 chlorotypes (C1–C25) recorded across the three species from 34 populations.
Table S3. The variable sites in aligned nrDNA sequences detected in the ITS4-ITS5 region of the three related species from 33 populations, which identified 14 haplotypes (N1–N14).
Table S4. Genetic diversity parameters estimated at 12 microsatellite loci in 33 populations.
Table S5. Bottleneck analysis for 33 populations of the three related species.
Table S6. The nine scenarios for four lineages.
Table S7. Prior distributions for model parameters used in model comparisons.
Table S8. The 19 bioclimatic variables used in the full set of variables, with √ indicating the variables used in the customized sets for Q. spinosa (QS), Q. rehderiana (QR), and Q. aquifolioides (QA); the variables used in the reduced set.
Table S9. Posterior median estimate and 95% highest posterior density interval (HPDI) for demographic parameters of an contraction-expansion model based on the nuclear multilocus microsatellite data of three oak species, respectively.
Table S10. The SSR dataset of these three oak species in this study.
Note S1. Description of the oak species and their habitats in our study.
Note S2. Methodological details for sequences and microsatellite genotyping.
References
Abbott, R., Albach, D., Ansell, S., Arntzen, J. W., Baird, S. J. E., Bierne, N., et al. (2013). Hybridization and speciation. J. Evol. Biol. 26, 229–246. doi: 10.1111/j.1420-9101.2012.02599.x
Aldrich, P. R., and Cavender-Bares, J. (2011). “Quercus,” in Wild Crop Relatives: Genomic and Breeding Resources, ed C. Kole (Berlin, Springer), 89–129.
Alleaume-Benharira, M., Pen, I., and Ronce, O. (2006). Geographical patterns of adaptation within a species' range: interactions between drift and gene flow. J. Evol. Biol. 19, 203–215. doi: 10.1111/j.1420-9101.2005.00976.x
An, Z. S., Kutzbach, J. E., Prell, W. L., and Porter, S. C. (2001). Evolution of Asian monsoons and phased uplift of the Himalaya-Tibetan plateau since Late Miocene times. Nature 411, 62–66. doi: 10.1038/35075035
Bandelt, H. J., Forster, P., and Röhl, A. (1999). Median-joining networks for inferring intraspecific phylogenies. Mol. Biol. Evol. 16, 37–48.
Beerli, P. (2006). Comparison of Bayesian and maximum-likelihood inference of population genetic parameters. Bioinformatics 22, 341–345. doi: 10.1093/bioinformatics/bti803
Cavender-Bares, J., González-Rodríguez, A., Eaton, D. A. R., Hipp, A. A. L., Beulke, A., and Manos, P. S. (2015). Phylogeny and biogeography of the American live oaks (Quercus subsection Virentes): a genomic and population genetics approach. Mol. Ecol. 24, 3668–3687. doi: 10.1111/mec.13269
Chen, D., Zhang, X., Kang, H., Sun, X., Yin, S., Du, H., et al. (2012). Phylogeography of Quercus variabilis based on Chloroplast dna sequence in East Asia: multiple glacial refugia and mainland-migrated island populations. PLoS ONE 7:e47268. doi: 10.1371/journal.pone.0047268
Chiucchi, J. E., and Gibbs, H. L. (2010). Similarity of contemporary and historical gene flow among highly fragmented populations of an endangered rattlesnake. Mol. Ecol. 19, 5345–5358. doi: 10.1111/j.1365-294X.2010.04860.x
Ciofi, C., Beaumontf, M. A., Swingland, I. R., and Bruford, M. W. (1999). Genetic divergence and units for conservation in the Komodo dragon Varanus komodoensis. Proc. R. Soc. Lond. B Biol. Sci. 266, 2269–2274. doi: 10.1098/rspb.1999.0918
Collins, W. D., Bitz, C. M., Blackmon, M. L., Bonan, G. B., Bretherton, C. S., Carton, J. A., et al. (2006). The community climate system model version 3 (CCSM3). J. Climate 19, 2122–2143. doi: 10.1175/JCLI3761.1
Cornuet, J. M., Pudlo, P., Veyssier, J., Dehne-Garcia, A., Gautier, M., Leblois, R., et al. (2014). DIYABC v2.0: a software to make approximate Bayesian computation inferences about population history using single nucleotide polymorphism, DNA sequence and microsatellite data. Bioinformatics 30, 1187–1189. doi: 10.1093/bioinformatics/btt763
Cornuet, J. M., Santos, F., Beaumont, M. A., Robert, C. P., Marin, J. M., Balding, D. J., et al. (2008). Inferring population history with DIY ABC: a user-friendly approach to approximate Bayesian computation. Bioinformatics 24, 2713–2719. doi: 10.1093/bioinformatics/btn514
Crawford, N. G. (2010). SMOGD: software for the measurement of genetic diversity. Mol. Ecol. Resour. 10, 556–557. doi: 10.1111/j.1755-0998.2009.02801.x
Degnan, J. H., and Salter, L. A. (2005). Gene tree distributions under the coalescent process. Evolution 59, 24–37. doi: 10.1111/j.0014-3820.2005.tb00891.x
Denk, T., and Grimm, G. W. (2010). The oaks of western Eurasia: traditional classifications and evidence from two nuclear markers. Taxon 59, 351–366. doi: 10.2307/25677595
Denk, T., and Tekleva, M. V. (2014). Pollen morphology and ultrastructure of Quercus with focus on Group Ilex (= Quercus Subgenus Heterobalanus (Oerst.) Menitsky): implications for oak systematics and evolution. Grana 53, 255–282. doi: 10.1080/00173134.2014.918647
Drummond, A. J., Suchard, M. A., Xie, D., and Rambaut, A. (2012). Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol. Biol. Evol. 29, 1969–1973. doi: 10.1093/molbev/mss075
Du, F. K., Peng, X. L., Liu, J. Q., Lascoux, M., Hu, F. S., and Petit, R. J. (2011). Direction and extent of organelle DNA introgression between two spruce species in the Qinghai-Tibetan Plateau. New Phytol. 192, 1024–1033. doi: 10.1111/j.1469-8137.2011.03853.x
Du, F. K., Petit, R. J., and Liu, J. Q. (2009). More introgression with less gene flow: chloroplast vs. mitochondrial DNA in the Picea asperata complex in China, and comparison with other conifers. Mol. Ecol. 18, 1396–1407. doi: 10.1111/j.1365-294X.2009.04107.x
Earl, D., and vonHoldt, B. (2012). STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv. Genet. Resour. 4, 359–361. doi: 10.1007/s12686-011-9548-7
Edelaar, P. I. M., and Björklund, M. (2011). If FST does not measure neutral genetic differentiation, then comparing it with QST is misleading. Or is it?. Mol. Ecol. 20, 1805–1812. doi: 10.1111/j.1365-294X.2011.05051.x
Epps, C. W., Castillo, J. A., Schmidt-Küntzel, A., du Preez, P., Stuart-Hill, G., Jago, M., et al. (2013). Contrasting historical and recent gene flow among African buffalo herds in the Caprivi Strip of Namibia. J. Hered. 104, 172–181. doi: 10.1093/jhered/ess142
Epps, C. W., and Keyghobadi, N. (2015). Landscape genetics in a changing world: disentangling historical and contemporary influences and inferring change. Mol. Ecol. 24, 6021–6040. doi: 10.1111/mec.13454
Evanno, G., Regnaut, S., and Goudet, J. (2005). Detecting the number of clusters of individuals using the software structure: a simulation study. Mol. Ecol. 14, 2611–2620. doi: 10.1111/j.1365-294X.2005.02553.x
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
Excoffier, L., and Ray, N. (2008). Surfing during population expansions promotes genetic revolutions and structuration. Trends Ecol. Evol. 23, 347–351. doi: 10.1016/j.tree.2008.04.004
Favre, A., Päckert, M., Pauls, S. U., Jähnig, S. C., Uhl, D., Michalak, I., et al. (2015). The role of the uplift of the Qinghai-Tibetan Plateau for the evolution of Tibetan biotas. Biol. Rev. Camb. Philos. Soc. 90, 236–253. doi: 10.1111/brv.12107
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.
Fu, Y. X., and Li, W. H. (1993). Statistical tests of neutrality of mutations. Genetics 133, 693–709.
Gay, L., Neubauer, G., Zagalska-Neubauer, M., Debain, C., Pons, J. M., David, P., et al. (2007). Molecular and morphological patterns of introgression between two large white-headed gull species in a zone of recent secondary contact. Mol. Ecol. 16, 3215–3227. doi: 10.1111/j.1365-294X.2007.03363.x
Gómez, J. M. (2003). Spatial patterns in long-distance dispersal of Quercus ilex acorns by jays in a heterogeneous landscape. Ecography 26, 573–584. doi: 10.1034/j.1600-0587.2003.03586.x
Goudet, J. (2001). FStat, A Program to Estimate and Test Gene Diversities and Fixation Indices (Version 2.9.3.2). Available online at: http://www2.unil.ch/popgen/softwares/fstat.htm
Guichoux, E., Garnier-Géré, P., Lagache, L., Lang, T., Boury, C., and Petit, R. J. (2013). Outlier loci highlight the direction of introgression in oaks. Mol. Ecol. 22, 450–462. doi: 10.1111/mec.12125
Hamrick, J. L., Godt, M. J. W., and Sherman-Broyles, S. L. (1992). “Factors influencing levels of genetic diversity in woody plant species,” in Population Genetics of Forest Trees, eds W. T. Adams, S. H. Strauss, D. L. Copes, and A. R. Griffin (Berlin, Springer), 95–124.
Hamrick, J., Linhart, Y., and Mitton, J. (1979). Relationships between life history characteristics and electrophoretically detectable genetic variation in plants. Annu. Rev. Ecol. Evol. Syst. 10, 173–200. doi: 10.1146/annurev.es.10.110179.001133
Heckman, K. L., Mariani, C. L., Rasoloarison, R., and Yoder, A. D. (2007). Multiple nuclear loci reveal patterns of incomplete lineage sorting and complex species history within western mouse lemurs (Microcebus). Mol. Phylogenet. Evol. 43, 353–367. doi: 10.1016/j.ympev.2007.03.005
Hedrick, P. W. (2005). A standardized genetic differentiation measure. Evolution 59, 1633–1638. doi: 10.1111/j.0014-3820.2005.tb01814.x
Hewitt, G. (2000). The genetic legacy of the Quaternary ice ages. Nature 405, 907–913. doi: 10.1038/35016000
Hewitt, G. M. (1996). Some genetic consequences of ice ages, and their role in divergence and speciation. Biol. J. Linn. Soc. 58, 247–276. doi: 10.1111/j.1095-8312.1996.tb01434.x
Hey, J. (2010). Isolation with migration models for more than two populations. Mol. Biol. Evol. 27, 905–920. doi: 10.1093/molbev/msp296
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 toolsfor spatial analysis of plant genetic resources data: 1, DIVA-GIS. Plant Genet. Res. Newslett. 127, 15–19.
Hudson, R. R. (1990). “Gene genealogies and the coalescent process,” in Oxford Surveys in Evolutionary Biology, eds D. J. Futuyma, J. D. Antonovics (New York, NY: Oxford University Press), 1–44.
Hutchison, D. W., and Templeton, A. R. (1999). Correlation of pairwise genetic and geographic distance measures: the relative influences of gene flow and drift on the distribution of genetic variability. Evolution 53, 1898–1914. doi: 10.2307/2640449
Jensen, J. L., Bohonak, A. J., and Kelley, S. T. (2005). Isolation by distance, web service. BMC genet. 6:13. doi: 10.1186/1471-2156-6-13
Leaché, A. D., Harris, R. B., Rannala, B., and Yang, Z. (2013). The influence of gene flow on species tree estimation: a simulation study. Syst. Biol. 63, 17–30. doi: 10.1093/sysbio/syt049
Li, L., Abbott, R. J., Liu, B. B., Sun, Y. S., Li, L. L., Zou, J. B., 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. Mol. Ecol. 22, 5237–5255. doi: 10.1111/mec.12466
Li, Z. H., Zou, J. B., Mao, K. S., Lin, K., Li, H. P., Liu, J. Q., et al. (2012). Population genetic evidence for complex evolutionary histories of four high altitude juniper species in the Qinghai-Tibetan Plateau. Evolution 66, 831–845. doi: 10.1111/j.1558-5646.2011.01466.x
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
Liston, A., Robinson, W. A., Piñero, D., and Alvarez-Buylla, E. R. (1999). Phylogenetics of Pinus (Pinaceae) based on nuclear ribosomal DNA internal transcribed spacer region sequences. Mol. Phylogenet. Evol. 11, 95–109. doi: 10.1006/mpev.1998.0550
Liu, J., Möller, M., Provan, J., Gao, L. M., Poudel, R. C., and Li, D. Z. (2013). Geological and ecological factors drive cryptic speciation of yews in a biodiversity hotspot. New Phytol. 199, 1093–1108. doi: 10.1111/nph.12336
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
Luo, D., Yue, J. P., Sun, W. G., Xu, B., Li, Z. M., Comes, H. P., et al. (2015). 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
Ma, C. L. (2006). Phylogeny and Biogeography of Quercus Sect. Heterobalanus. dissertation/PhD's thesis. [Kunming]: Kunming Institute of Botany, Chinese Academy of Sciences.
Ma, Y. Z., Li, Z. H., Wang, X., Shang, B. L., Wu, G. L., and Wang, Y. J. (2014). Phylogeography of the genus Dasiphora (Rosaceae) in the Qinghai-Tibetan Plateau: divergence blurred by expansion. Biol. J. Linn. Soc. 111, 777–788. doi: 10.1111/bij.12246
Maddison, W. P. (1997). Gene trees in species trees. Syst. Biol. 46, 523–536. doi: 10.1093/sysbio/46.3.523
Maddison, W. P., and Knowles, L. L. (2006). Inferring phylogeny pespite incomplete lineage sorting. Syst. Biol. 55, 21–30. doi: 10.1080/10635150500354928
Manos, P. S., Zhou, Z. K., and Cannon, C. H. (2001). Systematics of Fagaceae: phylogenetic tests of reproductive trait evolution. Int. J. Pl. Sci. 162, 1361–1379. doi: 10.1086/322949
Mantel, N. (1967). The detection of disease clustering and a generalized regression approach. Cancer Res. 27, 209–220.
McCormack, J. E., Zellmer, A. J., and Knowles, L. L. (2010). Does niche divergence accompany allopatric divergence in Aphelocoma jays as predicted under ecological speciation? Insights from tests with niche models. Evolution 64, 1231–1244. doi: 10.1111/j.1558-5646.2009.00900.x
Meirmans, P. G. (2014). Nonconvergence in Bayesian estimation of migration rates. Mol. Ecol. Resour. 14, 726–733. doi: 10.1111/1755-0998.12216
Meirmans, P. G., and Hedrick, P. W. (2010). Assessing population structure: FST and related measures. Mol. Ecol. Resour. 11, 5–18. doi: 10.1111/j.1755-0998.2010.02927.x
Mulch, A., and Chamberlain, C. P. (2006). Earth science: the rise and growth of Tibet. Nature 439, 670–671. doi: 10.1038/439670a
Myers, N., Mittermeier, R. A., Mittermeier, C. G., da Fonseca, G. A. B., and Kent, J. (2000). Biodiversity hotspots for conservation priorities. Nature 403, 853–858. doi: 10.1038/35002501
Nixon, K. C. (1993). Infrageneric classification of Quercus (Fagaceae) and typification of sectional names. Ann. For. Sci. 50, 25s–34s. doi: 10.1051/forest:19930701
Ortego, J., Noguerales, V., Gugger, P. F., and Sork, V. L. (2015). Evolutionary and demographic history of the Californian scrub white oak species complex: an integrative approach. Mol. Ecol. 24, 6188–6208. doi: 10.1111/mec.13457
Ortego, J., Riordan, E. C., Gugger, P. F., and Sork, V. L. (2012). Influence of environmental heterogeneity on genetic diversity and structure in an endemic southern Californian oak. Mol. Ecol. 21, 3210–3223. doi: 10.1111/j.1365-294X.2012.05591.x
Pamilo, P., and Nei, M. (1988). Relationships between gene trees and species trees. Mol. Biol. Evol. 5, 568–583.
Petit, R. J., Csaikl, U. M., Bordács, S., Burg, K., Coart, E., Cottrell, J., et al. (2002). Chloroplast DNA variation in European white oaks: phylogeography and patterns of diversity based on data from over 2600 populations. For. Ecol. Manag. 156, 5–26. doi: 10.1016/S0378-1127(01)00645-4
Petit, R. J., and Hampe, A. (2006). Some evolutionary consequences of being a tree. Annu. Rev. Ecol. Evol. Syst. 37, 187–214. doi: 10.1146/annurev.ecolsys.37.091305.110215
Phillips, S. J., Anderson, R. P., and Schapire, R. E. (2006). Maximum entropy modeling of species geographic distributions. Ecol. Model. 190, 231–259. doi: 10.1016/j.ecolmodel.2005.03.026
Piry, S., Luikart, G., and Cornuet, J.-M. (1999). BOTTLENECK: a program for detecting recent effective population size reductions from allele data frequencies. J. Hered. 86, 502–503.
Pons, O., and Petit, R. J. (1996). Measwring and testing genetic differentiation with ordered versus unordered alleles. Genetics 144, 1237–1245.
Posada, D. (2008). jModelTest: Phylogenetic model averaging. Mol. Biol. Evol. 25, 1253–1256. doi: 10.1093/molbev/msn083
Pritchard, J. K., Stephens, M., and Donnelly, P. (2000). Inference of population structure using multilocus genotype data. Genetics 155, 945–959.
Qiu, Y. X., Fu, C. X., and Comes, H. P. (2011). Plant molecular phylogeography in China and adjacent regions: tracing the genetic imprints of Quaternary climate and environmental change in the world's most diverse temperate flora. Mol. Phylogenet. Evol. 59, 225–244. doi: 10.1016/j.ympev.2011.01.012
Qu, Y. H., Luo, X., Zhang, R. Y., Song, G., Zou, F. S., and Lei, F. M. (2011). Lineage diversification and historical demography of a montane bird Garrulax elliotii-implications for the Pleistocene evolutionary history of the eastern Himalayas. BMC Evol Biol. 11–174. doi: 10.1186/1471-2148-11-174
Renner, S. S. (2016). Available data point to a 4-km-high Tibetan Plateau by 40 Ma, but 100 molecular-clock papers have linked supposed recent uplift to young node ages. J. Biogeogr. 43, 1479–1487. doi: 10.1111/jbi.12755
Rice, W. R. (1989). Analyzing tables of statistical tests. Evolution 43, 223–225. doi: 10.2307/2409177
Ronquist, F., and Huelsenbeck, J. P. (2003). MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19, 1572–1574. doi: 10.1093/bioinformatics/btg180
Rushton, B. S. (1993). Natural hybridization within the genus Quercus L. Ann. For. Sci. 50, 73s–90s. doi: 10.1051/forest:19930707
Schoener, T. W. (1968). The Anolis lizards of Bimini: resource partitioning in a complex fauna. Ecology 49, 704–726. doi: 10.2307/1935534
Selkoe, K. A., and Toonen, R. J. (2006). Microsatellites for ecologists: a practical guide to using and evaluating microsatellite markers. Ecol. Lett. 9, 615–629. doi: 10.1111/j.1461-0248.2006.00889.x
Sexton, J. P., Strauss, S. Y., and Rice, K. J. (2011). Gene flow increases fitness at the warm edge of a species' range. Proc. Natl. Acad. Sci. U.S.A. 108, 11704–11709. doi: 10.1073/pnas.1100404108
Sheppard, C. S. (2013). How does selection of climate variables affect predictions of species distributions? A case study of three new weeds in New Zealand. Weed Res. 53, 259–268. doi: 10.1111/wre.12021
Shi, Y. F., Li, J. J., and Li, B. Y. (1998). Uplift and Environmental Changes of Qinghai-Tibetan Plateau in the Late Cenozoic. Guangzhou: Guangdong Science and Technology Press.
Shi, Y. F., Li, J. J., Li, B. Y., Yao, T. D., Wang, S. M., Li, S. J., et al. (1999). Uplift of the Qinghai-Xizang (Tibetan) plateau and east Asia environmental change during late Cenozoic. J. Geog. Sci. 1, 10–20.
Simeone, M. C., Grimm, G. W., Papini, A., Vessella, F., Cardoni, S., Tordoni, E., et al. (2016). Plastome data reveal multiple geographic origins of Quercus Group Ilex. PeerJ 4:e1897. doi: 10.7717/peerj.1897
Sobel, J. M., Chen, G. F., Watt, L. R., and Schemske, D. W. (2010). The biology of speciation. Evolution 64, 295–315. doi: 10.1111/j.1558-5646.2009.00877.x
Sun, B. N., Wu, J. Y., Liu, Y. S. C., Ding, S. T., Li, X. C., Xie, S. P., et al. (2011). Reconstructing Neogene vegetation and climates to infer tectonic uplift in western Yunnan, China. Palaeogeogr. Palaeoclimatol. Palaeoecol. 304, 328–336. doi: 10.1016/j.palaeo.2010.09.023
Sun, Y. S., Li, L., Li, L. L., Zou, J. B., and Liu, J. Q. (2014). Distributional dynamics and interspecific gene flow in Picea likiangensis and P. wilsonii triggered by climate change on the Qinghai-Tibet Plateau. J. Biogeogr. 42, 475–484. doi: 10.1111/jbi.12434
Swofford, D. L. (2002). PAUP* version 4.0 b10. Phylogenetic Analysis Using Parsimony (* and Other Methods), Sunderland, MA: Sinauer.
Tajima, F. (1989). Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 123, 585–595.
Van Oosterhout, C., Hutchinson, W. F., Wills, D. P., and Shipley, P. (2004). MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data. Mol. Ecol. Notes 4, 535–538. doi: 10.1111/j.1471-8286.2004.00684.x
Warren, D. L., Glor, R. E., and Turelli, M. (2008). Environmental niche equivalency versus conservatism: quantitative approaches to niche evolution. Evolution 62, 2868–2883. doi: 10.1111/j.1558-5646.2008.00482.x
Warren, D. L., Glor, R. E., and Turelli, M. (2010). ENMTools: a toolbox for comparative studies of environmental niche models. Ecography 33, 607–611. doi: 10.1111/j.1600-0587.2009.06142.x
Weir, B. S., and Cockerham, C. C. (1984). Estimating F-statistics for the analysis of population structure. Evolution 38, 1358–1370. doi: 10.2307/2408641
Welt, R. S., Litt, A., and Franks, S. J. (2015). Analysis of population genetic structure and gene flow in an annual plant before and after a rapid evolutionary response to drought. AoB Plants 7:plv026. doi: 10.1093/aobpla/plv026
Wen, J., Zhang, J., Nie, Z. L., Zhong, Y., and Sun, H. (2014). Evolutionary diversifications of plants on the Qinghai-Tibetan Plateau. Front. Genet. 5:4. doi: 10.3389/fgene.2014.00004
Wilson, G. A., and Rannala, B. (2003). Bayesian inference of recent migration rates using multilocus genotypes. Genetics 163, 1177–1191.
Wilyard, A., Cronn, R., and Liston, A. (2009). Reticulate evolution and incomplete lineage sorting among the Ponderosa pines. Mol. Phylogenet. Evol. 52, 498–511. doi: 10.1016/j.ympev.2009.02.011
Wolfe, K. H., Li, W. H., and Sharp, P. M. (1987). Rates of nucleotide substitution vary greatly among plant mitochondrial, chloroplast, and nuclear DNAs. Proc. Natl. Acad. Sci. U.S.A. 84, 9054–9058.
Wu, C. Y. (1987). “Origin and evolution of flora of Xizang,” in Flora Xizangica ed C. Y. Wu (Beijing: Science Press), 879–902.
Xiao, Z. S., Gao, X., Jiang, M. M., and Zhang, Z. B. (2009). Behavioral adaptation of Pallas's squirrels to germination schedule and tannins in acorns. Behav. Ecol. 20, 1050–1055. doi: 10.1093/beheco/arp096
Yeh, F. C., Yang, R. C., Boyle, T. B. J., Ye, Z. H., and Mao, J. X. (1997). POPGENE, the User-Friendly Shareware for Population Genetic Analysis. Edmonton: Molecular Biology and Biotechnology Centre, University of Alberta.
Yin, A., and Harrison, T. M. (2000). Geologic evolution of the Himalayan-Tibetan orogen. Annu. Rev. Earth Planet. Sci. 28, 211–280. doi: 10.1146/annurev.earth.28.1.211
Zeng, Y. F., Liao, W. J., Petit, R. J., and Zhang, D. Y. (2010). Exploring species limits in two closely related Chinese oaks. PLoS ONE 5:e15529. doi: 10.1371/journal.pone.0015529
Zeng, Y. F., Liao, W. J., Petit, R. J., and Zhang, D. Y. (2011). Geographic variation in the structure of oak hybrid zones provides insights into the dynamics of speciation. Mol. Ecol. 20, 4995–5011. doi: 10.1111/j.1365-294X.2011.05354.x
Zeng, Y. F., Wang, W. T., Liao, W. J., Wang, H. F., and Zhang, D. Y. (2015). Multiple glacial refugia for cool-temperate deciduous trees in northern East Asia: the Mongolian oak as a case study. Mol. Ecol. 24, 5676–5691. doi: 10.1111/mec.13408
Zhou, Z. K. (1992). Origin, phylogeny and dispersal of Quercus from China. Acta Botanica Yunnanica 14, 227–236.
Zhou, Z. K. (1993). A taxonomical revision of fossil evergreen Sclerophyllous oaks from China. J. Integr. Plant Biol. 12, 954–961.
Keywords: East Himalaya-Hengduan Mountains, ecological niche models, environmental heterogeneity, evolutionary history, gene flow, Quercus
Citation: Feng L, Zheng Q-J, Qian Z-Q, Yang J, Zhang Y-P, Li Z-H and Zhao G-F (2016) Genetic Structure and Evolutionary History of Three Alpine Sclerophyllous Oaks in East Himalaya-Hengduan Mountains and Adjacent Regions. Front. Plant Sci. 7:1688. doi: 10.3389/fpls.2016.01688
Received: 25 July 2016; Accepted: 26 October 2016;
Published: 11 November 2016.
Edited by:
Tian Tang, Sun Yat-sen University, ChinaReviewed by:
Gong Xun, Chinese Academic of Sciences, ChinaPei-Chun Liao, National Taiwan Normal University, Taiwan
Copyright © 2016 Feng, Zheng, Qian, Yang, Zhang, Li and Zhao. 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) or licensor 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: Gui-Fang Zhao, Z2Z6aGFvQG53dS5lZHUuY24=