Genetic Differentiation and Demographic History of Three Cerris Oak Species in China Based on Nuclear Microsatellite Makers

Knowledge of interspecific divergence and population expansions/contractions of dominant forest trees in response to geological events and climatic oscillations is of major importance to understand their evolution and demography. However, the interspecific patterns of genetic differentiation and spatiotemporal population dynamics of three deciduous Cerris oak species (Q. acutissima, Q. variabilis and Q. chenii) that are widely distributed in China remain poorly understood. In this study, we genotyped 16 nuclear loci in 759 individuals sampled from 44 natural populations of these three sibling species to evaluate the plausible demographical scenarios of the closely related species. We also tested the hypothesis that macroand microevolutionary processes of the three species had been triggered and molded by Miocene–Pliocene geological events and Quaternary climatic change. The Bayesian cluster analysis showed that Q. acutissima and Q. chenii were clustered in the same group, whereas Q. variabilis formed a different genetic cluster. An approximate Bayesian computation (ABC) analyses suggested that Q. variabilis and Q. acutissima diverged from their most common ancestor around 19.84 Ma, and subsequently Q. chenii diverged from Q. acutissima at about 9.6 Ma, which was significantly associated with the episodes of the Qinghai–Tibetan Plateau (QTP). In addition, ecological niche modeling and population history analysis showed that these three Cerris oak species repeatedly underwent considerable ‘expansion–contraction’ during the interglacial and glacial periods of the Pleistocene, although they have varying degrees of tolerance for the climatic change. Overall, these findings indicated geological and climatic changes during the Miocene–Pliocene and Pleistocene as causes of species divergence and range shifts of dominant tree species in the subtropical and warm temperature areas in China.


Introduction
Historical and ecological factors such as geological and/or climatic processes have a profound impact on biodiversity and sustainability at every level of the biota, from genes to spatial distribution [1][2][3][4]. The interactions among these processes would create multiple ecological niches, favoring high rates of intra-and interspecific diversification [5,6]. In parallel, climate buffering associated with complex topography seems to have contributed to a particularly low extinction rate of species/lineages adapted to contrasting environmental conditions since the Tertiary era [7,8]. Mountains uplift and climatic oscillations have exerted a significant influence on current geographic distribution of biodiversity, especially for the interspecific divergence and population demography of plants [9,10]. However, a knowledge of interspecific divergence and range dynamics for sedentary and long-lived tree species in response to the past geological and climate change is informative not only to know how forests were affected by historically environmental factors, but also to elucidate the components of such forests might respond to future environmental change [11].
Quercus L. (oak) is the most widespread and species-rich tree genus in the Northern Hemisphere [12]. This genus characterized by a strong potential for ecological adaptation, and species in the genus are usually have a high potential for introgression and reticulate evolution (e.g., [8,13,14]). Genus Quercus has recently been formalized as two subgenera and eight sections based on the restriction-site associated DNA sequencing (RAD-sequencing) data [15]. The predominantly subgenus Quercus, namely the New World oaks, includes the sections such as Quercus, Lobatae, Protobalanus, Virentes and Ponticae. The subgenus Cerris, namely the Old World oaks, includes Cerris, Ilex and Cyclobalanopsis sections. Among these sections, the Eurasian sect. Cerris comprises about 15 species, most of which are confined to western Eurasia [16], and three deciduous species (Q. acutissima, Q. chenii and Q. variabilis) that are widely distributed in Eastern Asia [17][18][19][20]. In recent years, detailed phylogeographic inferences of the three Cerris oaks in Eastern Asia had been generally examined with plastid DNA and nuclear markers; the results suggested climatic factors and environmental heterogeneity might play key roles in driving the intraspecific genetic divergence [21][22][23][24][25]. However, limited information was obtained to resolve interspecific divergence and demographic history of the three closely related Cerris oak species (Q. acutissima, Q. chenii and Q. variabilis) that widely distributed in China.
According to the State Forestry Administration of China, oaks comprise 13% of the natural forests in China. Q. acutissima Carruthers is a dominant species occupying a wide range of environmental conditions in subtropical and warm temperature zones in China [25]. This species is often found in deciduous broad-leaved forests on hillsides, which is regarded as the most common oak in deciduous secondary forests in the lowlands and a component of deciduous evergreen forests at low latitudes [26]. The hardwood of Q. acutissima could be used as the excellent building material, and the fresh leaves are a favorite food for Chinese oak silkworms. Another noteworthy usage of Q. acutissima is edible fungi cultivation [27]. Due to the high ecological and economic values, Q. acutissima has been listed as a precious timber species in China and also listed as the least concern according to the IUCN Red List of Threatened Species v3.1. (https://www.iucnredlist. org/species/194048/2295121) (accessed on 24 May 2021). Q. variabilis Blume is also the dominant species for the forests spread in the northern region of the subtropical areas of China, it usually has the sympatric distribution with Q. acutissima [28]. However, the geography of Q. chenii Nakai (Q. acutissima ssp. chenii Camus) is mainly restricted to central and eastern China [28]. This species can form pure forests and is distributed in the deciduous broad-leaves forests in hilly landscape [29]. Due to the widespread distribution and dominant position in warm temperate deciduous forests in China, these three Cerris oaks could provide a useful model to evaluate the impact of the geographical change and climatic oscillations on the interspecific divergence and demographic history of sibling tree species at a larger geographic scale. Understanding these abiotic factors in driving interspecific differentiation and species dynamics is a major task in ecological and evolutionary research [8,30,31].
In this study, we collected 759 individuals of the three Cerris oaks (Q. acutissima, Q. chenii and Q. variabilis) covering their main distribution ranges across China. Sixteen expressed sequence tag-single sequence repeat (EST-SSR) loci were surveyed for population genetic analyses, and ecological niche modelling on these three Cerris oaks were also conducted. We aim to obtain detailed scenarios of phylogeography and historical demography of the three important oak species in China. In particular, (1) we analyzed patterns of interspecific genetic differentiation among the three East Asia Cerris oaks and employed an approximate Bayesian computation (ABC) framework to compare different plausible scenarios of species divergence, which allowed us to infer the evolutionary and demographic history of the three relatives in section Cerris. (2) To test the hypothesis that interspecific divergence occurred progressively in the related oak species in response to the mountain uplifting, such as QTP uplifts. Meanwhile, we want to infer whether pairwise species genetic differentiation is associated with niche overlap. (3) To test the hypothesis that these related species underwent historical 'expansion-contraction' following the recent uplifts and/or the Pleistocene glacial/interglacial periods in this region.

Population Sampling
Leaf materials were collected from 44 wild populations (N = 759) of the three Cerris oaks, which distributed in subtropical China, including 19 populations (N = 336) of Q. variabilis, 20 populations (N = 338) of Q. acutissima, and 5 populations (N = 85) of Q. chenii (Table S1, Figure 1). For each population, mature foliage was sampled from 6-20 individuals spaced at least >100 m apart and dried rapidly with silica gel and stored at −10 • C until DNA isolation. Voucher specimens of each population were temporarily stored in the Evolutionary Botany Lab of Northwest University.

DNA Extraction and Microsatellite Genotyping
Total genomic DNA for each individual was extracted from the dried leaf tissue using the Plant Genomic DNA Extraction Kit (Tiangen, Beijing, China). A total of 16 nSSR loci were designed using Batch Primer 3 [32] based on Q. variabilis microsatellite sequences deposited in GenBank [33]. All 16 nSSR loci have high polymorphism and stably amplified among the three species (Table S2). Polymerase chain reactions (PCR) were carried out in a SimpliAmp TM Thermal Cycler (Applied Biosystems, Carlsbad, CA, USA) with a volume of 20 µL containing 10 µL 2 × Taq PCR Master Mix, 1 µL of primer template (10-30 ng), and 0.3 µL of each primer, 8.4 µL double-distilled H 2 O. The thermo-cycling conditions were set as follows: 5 min at 94 • C, followed by 30 cycles of 30 s at 94 • C, 40 s annealing at 55 • C, 40 s extension at 72 • C, with a final extension cycle of 5 min at 72 • C. PCR products were examined on an ABI 3730 XL analyzer (Applied Biosystems, Foster City, CA, USA), and the data were scored and compiled using GeneMarker ® v2.2 (SoftGenetics, State College, PA, USA).

Genetic Diversity, Genetic Structure and Gene Flow Analysis
We calculated the genetic diversity indexes, including observed heterozygosity (H O ), expected heterozygosity (H E ), and polymorphism information content (PIC) for each locus using software CERVUS v3.0.7 and FSTAT v2.9.3 [34,35]. For each population, the genetic diversity parameters were estimated using GENEALEX v6.5 [36] and FSTAT v2.9.3, including number of different alleles (Na), number of effective alleles (Ne), inbreeding coefficient (F IS ), H O , and H E . We also using FSTAT v2.9.3 to tested Hardy-Weinberg equilibrium (HWE) for all locus in each population.
Genetic structure of populations for the three oak species were surveyed using Bayesian clustering approach implemented in STRUCTURE v2.2.4 [37]. The analyses were run using the admixture model. A total of 10 independent runs were performed for the K value from one to eight with 200,000 burn-in steps followed by 500,000 Markov Chain Monte Carlo (MCMC) steps. The STRUCTURE output was analyzed by Structure Harvester (http://taylor0.biology.ucla.edu/structureHarvester/ accessed on 13 August 2020) and the most likely numbers of clusters (K) was estimated according to Pritchard et al. (2000) and Evanno et al. (2005) [37,38]. Graphics were displayed in the DISTRUCT program [39].

DNA Extraction and Microsatellite Genotyping
Total genomic DNA for each individual was extracted from the dried leaf tissue using the Plant Genomic DNA Extraction Kit (Tiangen, Beijing, China). A total of 16 nSSR loci were designed using Batch Primer 3 [32] based on Q. variabilis microsatellite sequences deposited in GenBank [33]. All 16 nSSR loci have high polymorphism and stably amplified among the three species (Table S2)  We also performed AMOVA analysis by pooling all the three oaks as the same species samples, and three specific analyses considering separately species that were identified with the analysis in ARLEQUIN v3.11 [40]. Pairwise F ST was used to assess population differentiation within and among species and populations [41].
In order to infer historical gene flow (Nm) patterns, we used a maximum likelihood procedure and a Brownian motion microsatellite model in MIGRATE-N v3.6 [45] to estimate the effective population sizes (θ) and mutation scaled immigration (M) among the groups with 10 short chains following 5000 steps and three long chains of 50,000 iterations, sampling every 100 steps under a constant mutation model and discarding the first 10,000 records as a burn-in. In addition, to assess patterns of recent migration, we also estimated inter-species migration rates using a Bayesian approach in BAYESASS v3.0 [46].

Species Divergence and Demographic History
To estimate plausible of divergence and populations dynamics of the three closely oak species, Approximate Bayesian Computation (ABC) analysis was performed in software DIYABC v2.0.4 [47]. Four possible historical divergence scenarios for these three oak species were built by the DIYABC (Figure 2A). The model 1 consisted of a simultaneous divergence of these three oak species (QA: Q. acutissima, QV: Q. variabilis and QC: Q. chenii). The differences between model 2 and model 3 was that QC derived from QA or QV, respectively. The model 4 was considered to be admixture model that QC was hypothesized to be the admixture of QA and QV. We ran ABC analysis with 1,000,000 simulations for each model. To identify the most likely divergence scenario, the posterior probability of model was assessed using a logistic regression on the 1% of the stimulated data sets closest to the observed data and estimated the relative posterior probability (PP) with 95% confidence intervals (95% CI) for each model. We also assess confidence in power model by estimating type I and type II error rates with 500 pseudo-observed data sets (PODs) under each scenario. The performance of parameter estimation was determined by calculating the median of the relative median of absolute errors (RAME) using 500 PODs. We also applied a principal component analysis (PCA) for checking model and visualize the fit between simulated and observed data sets. In addition, DIYABC was used to simulate and examine the effective population sizes of each species. We separately tested the four scenarios of demographic changes for the three oak species: continuous expansion, continuous shrinkage, shrinkage-expansion, and expansion-shrinkage ( Figure 2B). The prior distributions of demographic parameters were listed in Tables S3 and S4. under each scenario. The performance of parameter estimation was determined by calculating the median of the relative median of absolute errors (RAME) using 500 PODs. We also applied a principal component analysis (PCA) for checking model and visualize the fit between simulated and observed data sets. In addition, DIYABC was used to simulate and examine the effective population sizes of each species. We separately tested the four scenarios of demographic changes for the three oak species: continuous expansion, continuous shrinkage, shrinkage-expansion, and expansion-shrinkage ( Figure 2B). The prior distributions of demographic parameters were listed in Tables S3 and S4.
To test the recent dynamics of effective population size (bottle-necks), we tested the genetic bottleneck among all populations in software BOTTLENECK v1.2.02 [48]. We used the Wilcoxon signed-rank necks under the infinite allele model (IAM) and the stepwise mutation model (SMM) [49].

Niche Comparison Analyses on Geographic Space (G-Space) Environment Space (E-Space)
In order to capture ecological differences of three oak species, we conducted comparative analysis in both geographic (G) and environment (E) space for the three oaks, these two types of niche analyses would complement each other in the niche comparison analyses [50]. Historical shifts in geographic distribution across different times, including the Last Interglacial (LIG, ~130 Ka), Last Glacial Maximum (LGM, ~22 Ka) and present, which were performed in the MAXENT v3.4.1 [51]. In these analyses, we selected 14 bioclimatic factors according to the previous studies (Table S5) [21,22,24]. Climatic variables at 2.5 arc- To test the recent dynamics of effective population size (bottle-necks), we tested the genetic bottleneck among all populations in software BOTTLENECK v1.2.02 [48]. We used the Wilcoxon signed-rank necks under the infinite allele model (IAM) and the stepwise mutation model (SMM) [49].

Niche Comparison Analyses on Geographic Space (G-Space) Environment Space (E-Space)
In order to capture ecological differences of three oak species, we conducted comparative analysis in both geographic (G) and environment (E) space for the three oaks, these two types of niche analyses would complement each other in the niche comparison analyses [50]. Historical shifts in geographic distribution across different times, including the Last Interglacial (LIG,~130 Ka), Last Glacial Maximum (LGM,~22 Ka) and present, which were performed in the MAXENT v3.4.1 [51]. In these analyses, we selected 14 bioclimatic factors according to the previous studies (Table S5) [21,22,24]. Climatic variables at 2.5 arcmin resolution were obtained from the WorldClim database (http://www.worldclim.org/) (accessed on 1 July 2019). In order to estimate the suitable area gains and losses (or unchanged areas) for past scenario with present, the overlapped areas were predicted using Intersect tool of ArcGIS v. 10.2 (ESRI, Redlands, CA, USA).
In addition, to measure ecological differences between each pair of species, we used the ENMTOOLS v1.3 [52] to calculate two alternative statistics: I and D [52,53], which range from 0 (no overlap) to 1 (complete overlap). To test the null hypothesis that these three oak species have identical environmental space, we used the niche identity test initially introduced by Warren et al. (2008) [52]. In this test, we compared the overlap of a species pair's actual niches to the distribution of niche overlaps obtained from pairs of pseudoniches (n = 100 pseudoreplicates) (see [52]).
To assess the degree of the niche overlap, we first assessed ENMs in spatial environment using the principal components analysis (PCA; hereafter standard PCA) to examine the climatic variability of the realized niches in the genetically studied populations across the total climatic space. A simple agglomerative hierarchical clustering method (UPGMA) was conducted based on the means of the principal component (PC). After that, following the approach developed by Broennimann et al. [54], and the R script from Herrando-Moraira et al. [55], the PCA was used to translate occurrences and climate data into environment axes (PCA-env), and densities of points in multidimensional E space were used to test the niche overlap among these three oak species, and calculated the D statistics.

Genetic Diversity
The genetic diversity indexes for each locus were summarized in Tables S6 and S7. The values of Ho were much lower than H E (Table S1), and these three oak species have different degrees of deviation from HWE in most populations (Table S8), which were associated with excess of homozygous caused by inbreeding or mating among relatives (Table S1). In addition, these populations SPU, SPg, SPI and SPWS of Q. variabilis, E, J. HF, CY, P, YS, YK and PS of Q. acutissima (Table S1) had higher genetic diversity, and most of them were located in Yunnan-Guizhou Plateau (e.g., SPU/E), Qinling mountains (e.g., SPg/J, SPWS) and Nanling mountains (e.g., YS, PS).

Genetic Differentiation and Structure Analysis
The AMOVA analysis detected significant genetic differentiation among the three species groups (F CT = 0.195, p < 0.001), but only 19.52% of the genetic variation was partitioned. At the species level, genetic variation mainly occurred within populations (91.03%, 78.63% and 94.55%, F ST = 0.090, 0.214 and 0.055, respectively) with all F ST values having highly significant estimations (p < 0.001) ( Table 1). According to the STRUCTURE analysis, the results of the ∆K method suggested that the genetic structure of the three Cerris oak species can be explained by two genetic clusters ( Figure S1). The populations of Q. acutissima and Q. chenii formed one group, whereas the other cluster comprised all populations from Q. variabilis, and both groups seemed to be mixed with a little gene flow (Figure 1a). When K = 3, the former group was further subdivided, as Q. acutissima and Q. chenii (Figure 1b). Meanwhile, we detected a shared signal of genetic components between Q. acutissima and Q. chenii populations (Figure 1b). Subsequent hierarchical analyses within Q. acutissima and Q. variabilis (except for Q. chenii, due to the limited populations) showed a further genetic subdivision ( Figure S2).
The UPGMA dendrograms based on Nei's (1972) genetic distance (Ds) suggested that the three closely related oak species were completely separated (Figure 1c). A similar result was obtained in the PCoA analysis (Figure 1d), where three obviously differentiated groups (Q. variabilis, Q. acutissima and Q. chenii) were found with few individuals mixed among the three species.

Genetic Migration among Species
The BAYESASS analysis detected recent but very low levels of genetic migration from Q. acutissima to Q. chenii (m = 0.024) and from Q. chenii to Q. acutissima (m = 0.004), Q. acutissima to Q. variabilis (m = 0.001) and from variabilis to Q. acutissima (m = 0.002) (Table S9). In addition, the Migrate-n analysis identified historical asymmetric gene flow (Nm) among the three oak species, with greater Nm between Q. acutissima and Q. variabilis (Table 2).
Posterior parameter estimates (median and 95% confidence intervals) for the best supported demographical model 2.   Table S4). Under the TPM model and the SMM model in the BOTTLENECK, the Wilcoxon's signed rank test revealed population SPF, belonged to the Q. variabilis and CY, YS belonged to the Q. acutissima had undergone possible bottleneck effect (Table S11).

Ecological Niche Modeling
All ENM simulations had high AUC values (>0.85) and predicted well for current species distributions for the three oak species (Figure S8a-i). The projected distributions into past environments suggested that the distribution ranges of the three related species contracted during the LGM periods, whereas with considerable expansion at present, and these three oak species' distributions were more restricted during the LGM than at LIG (Figures S8a-i and 3a-f). It is also remarkable the isolated distribution of Q. chenii during the LIG and LGM, showed very little range overlap with Q. variabilis and Q. acutissima (Figure 3). In addition, all ENMs indicated that the predicated distributions of Q. variabilis and Q. acutissima were considerably overlapped during three periods (Figures S8 and 3). According to a standard PCA analysis, the first three principal components accounted for 67.87% of the total climatic variation (PC1 = 40.50%, PC2 = 27.37%) ( Figure  S9A). The first principal component (PC1) was mainly explained by BIO 17 (Precipitation of the driest quarter), and the second (PC2) was explained by BIO 4 (Temperature seasonality) ( Figure S9A). According to the PCA-env analysis, the first two axes explained 75.23% of total variation in the climatic conditions across the population distributions of the three oak species (PC1 = 49.46% and PC2 = 25.77%), where the first component (PC1) was mainly explained by BIO 1 (Annual precipitation) and the second was explained by BIO 10 (Mean temperature of the warmest quarter) (Figures S9B,C and S10). Multiple niche plots for the 20% occurrence density showed that Q. variabilis and Q. acutissima shared their climate niche to a great extent ( Figure S8B). The occurrence density of 100% plotted points in the PCA-env space showed similar results with 20% of the occurrence density ( Figure S9B,C). Besides, ecological niche revealed differentiation between Q. chenii and Q. variabilis/Q. acutissima ( Figure S9B,C). The overlap D values between Q. variabilis and Q. acutissima was relatively high (D = 0.687); however, the generally low overlap D values were found between Q. variabilis/Q. acutissima and Q. chenii (D = 0.025; D = 0.038). According to a standard PCA analysis, the first three principal components accounted for 67.87% of the total climatic variation (PC1 = 40.50%, PC2 = 27.37%) ( Figure S9A). The first principal component (PC1) was mainly explained by BIO 17 (Precipitation of the driest quarter), and the second (PC2) was explained by BIO 4 (Temperature seasonality) ( Figure S9A). According to the PCA-env analysis, the first two axes explained 75.23% of total variation in the climatic conditions across the population distributions of the three oak species (PC1 = 49.46% and PC2 = 25.77%), where the first component (PC1) was mainly explained by BIO 1 (Annual precipitation) and the second was explained by BIO 10 (Mean temperature of the warmest quarter) (Figures S9B,C and S10). Multiple niche plots for the 20% occurrence density showed that Q. variabilis and Q. acutissima shared their climate niche to a great extent ( Figure S8B). The occurrence density of 100% plotted points in the PCA-env space showed similar results with 20% of the occurrence density ( Figure S9B,C). Besides, ecological niche revealed differentiation between Q. chenii and Q. variabilis/Q. acutissima ( Figure S9B,C). The overlap D values between Q. variabilis and Q. acutissima was relatively high (D = 0.687); however, the generally low overlap D values were found between Q. variabilis/Q. acutissima and Q. chenii (D = 0.025; D = 0.038).

Species Divergence
Accumulating evidence has shown that the tectonic-climatic interactions since the Miocene have significantly contributed to the evolutionary process of plant species in East Asia [24,57,58]. By using population-genetic data obtained from nuclear loci, our Bayesian clustering analysis results indicated that Q. chenii and Q. acutissima were closely related (K = 2), Q. acutissima and Q. chenii were clustered in the same group, whereas the Q. variabilis formed its own genetic cluster; when K = 3, Q. acutissima and Q. chenii were split into different clusters, with some individuals of these two species behave as mosaics (Figure 1a,b). The PCoA and UPGMA results also indicated that there were three genetic groups, which highly consistent with the Bayesian clustering analysis (Figure 1c,d). Furthermore, the scenario of ancient divergence was strongly supported by the model 2 in the approximate Bayesian computation analysis (Table S10, Figure S3). Assuming an average generation time of 80 years for Quercus species, Q. variabilis and Q. acutissima was inferred to have diverged from a most common ancestor around 19.84 Ma (95% CI: 16.24-23.76), whereas Q. chenii diverged from Q. acutissima at about 9.6 Ma (95% CI: 5.07-15.20 Ma) (Figure 2A, Table 3). A recent molecular dating based on chloroplast-nuclear data matrix showed that these three Cerris oaks had initially diverged around 24.56 Ma [59], to some extent, it was earlier than our molecular dating (~19.84 Ma) using the nuclear microsatellite data, whereas the 95% HPD were relatively identical (18.47-30.44 Ma for chloroplast-nuclear data matrix and 16.24-23.76 Ma for present study). In addition, the divergence of Q. chenii (9.6 Ma, 95% CI: 5.07-15.20 Ma) is comparable with previous estimates for chloroplast DNA haplotypes (11.63-6.21Ma) [24]. Although we must be cautious when regarding the results of molecular dating, these divergence times consistent with the first fruit specimen of Quercus belonging to the Cerris group found in the middle Miocene deposits in Shanwang, China [28]. Our analyses indicated that the divergence of the East Asia Cerris oaks took place in the Miocene (including divergence time among Q. variabilis and Q. acutissima at 19.84 Ma, divergence of Q. chenii around 9.6 Ma), which coincided with a period of intense uplift of the Hengduan Mountain massif. It has been recognized that the QTP uplifting began around 40 Ma, and more extensive uplifts further occurred at 20-22 Ma, followed by several additional uplifts during the Miocene [60][61][62][63][64]. Such geological changes may have produced many small fragmented habitats with different microclimates, thereby impacting the direction of natural selection [65]. Indeed, the rise of the QTP has dramatically modified the global and East Asian climate by triggering and intensifying the Asian monsoon. These new climatic conditions and heterogeneity environment would promote the species differentiation. An association between Miocene uplift of QTP and interspecific/intraspecific divergence time has previously been noted for other studies of tree, shrub and herb, such as Myricaria (~20 Ma) [66], Cyclocarya paliurus (~16.69 Ma) [57] and Notopterygium (~10.90 Ma) [67]. In addition, the niche overlap analysis showed a remarkable niche differentiation among Q. chenii and Q. variabilis/Q. acutissima by the identity test (D = 0.406, I = 0.698; D = 0.389, I = 0.678) and E-space analyses (D = 0.025, D = 0.038) (Figures S8k,l and S9). Ecological niche partitioning would have reinforced the divergence of the lineages following initial spatial isolation and might bring differentiated adaptability to their respective environmental conditions [68]. Furthermore, geographical isolation and environmental heterogeneity also have the impact on the genetic differentiation of other Quercus species, such as three alpine sclerophyllous oaks (Q. aquifolioides, Q. spinosa and Q. rehderiana) [69] and Californian scrub white oak species (Q. berberidifolia and Q. durata) [8]. In fact, our G-space analysis indicated that the potential distributions of Q. chenii were significantly differ from the other two oak species during the LIG and LGM (Figures S8d-f and 3). However, the E-space and Gspace analyses (Figures S8 and S9) could not meticulously tease apart geographic isolation signals of potential distribution and particular niche possession for the two related species Q. variabilis and Q. acutissima (e.g., Figure S8a,b,d,e,j), which probably reflects multiple factors that jointly promote genetic differentiation among the three closely related taxa during the early stage of species formation. Therefore, we still need to use genetic makers derived high throughout sequencing and available bioclimatic models to address in more detailed environmental selection and divergent adaptation of these related Cerris oak species in the future.
BAYESASS and Migrate-N analysis provided unambiguous evidences of the interspecific bidirectional asymmetric gene flow among the three oak species (Table 2 and S9). The historical gene flow from Q. acutissima and Q. variabilis to Q. chenii was lower than that from Q. chenii to Q. acutissima and Q. variabilis (Table 2). In contrast, the contemporary gene flow from Q. acutissima to Q. chenii was higher than that in the opposite direction (Table S9). This pattern might be explained by the fact that Q. acutissima and Q. variabilis were the earliest differentiated species, and Q. chenii was the most recently diverged species. In addition, oak species were paraded as typical examples of the intricate concept of biological species due to their high potential for interspecies introgression. Simultaneously, oak species usually brought difficulties to the phylogenetic tree inference (e.g., [16,70]). In our studies, the Bayesian clustering analysis also have shown considerable rates of hybridization among the three related species (Figure 1a,b). The appearance of introgression was commonly detected in regions of sympatry/parapatry among related species, as is well exemplified in the codistributed populations of Q. berberidifolia and Q. Cornelius-mulleri in southern California, and Q. berberidifolia and Q. durata in central-northern California [8,71]. Albeit some populations showed genetic hybridization (e.g., SPW, SPY, SPU, SPS, CS and CY), all the putative species basically maintain their taxonomic and keep lineages separated to a high degree (Figure 1a,b). Asynchronous flowering phenology or selection against hybrids with lower survivability or tending to parent environment might be contribution of maintaining the distinctness of three studied species in the presence of gene flow [8,72,73]. Furthermore, the climate seasonality (e.g., temperature and precipitation seasonality) as the most important ecological factors that might affect the phenology (including flowering time and growing season) [74,75]. In fact, these environmental variables detected in the E-space analysis, e.g., BIO 4 (temperature seasonality) BIO 10 (mean temperature of the warmest quarter) and BIO 17 (precipitation of the driest quarter), were considered as potential abiotic factors in triggering the genetic differentiation among the three oak species ( Figure S9). In addition, our analyses also indicated a delicate intraspecific genetic structure of Q. acutissima and Q. variabilis across their main distribution areas ( Figure S2), similar to previous studies that suggested genetic subdivision in these two species and other oak species [21,23,25,76]. Therefore, these results suggested that widespread gene flow in the wind-pollinated oak trees is a remarkable homogenizing factor within species, while selection against hybrids or assortative mating might plausibly prevent the three closely related species from converging into mosaic/hybrid regions [8,72].

Demographic History
It is generally assumed that plant species have experienced important distributional shifts driven by glacial-interglacial cycles during the Pleistocene [77,78]. The current ABC simulations of population-genetic data indicated that Q. variabilis has experienced an earlier population contraction (~0.95 Ma), followed by a recent population expansion (~0.19 Ma) ( Figure 2B, Table S4). Such demographic events could be linked to the glacial/interglacial recorded in the QTP, and the most severe glaciation on the QTP reached its maximum between 0.8 and 0.6 Ma when the ice sheet is thought to have covered an area five to seven larger areas than it does today [79]. In fact, the uplift of the QTP might have caused habitat disruption, and the unfavorable environment might lead to extinction of some populations while producing population contraction for others [67]. In terms of the recent population expansion models, previous studies found that the antepenultimate interglacial period and the last interglacial period took place between 0.5 and 0.3 Ma, and between 0.14 and 0.12 Ma, respectively [80]. Besides, the ABC analysis suggested that Q. chenii experienced a population expansion around 0.32 to 0.12 Ma (Table S4). Indeed, the results of our ENM analysis indicated that Q. chenii and Q. variabilis expanded their ranges in the relatively warm interglacial period (~0.14-0.12 Ma) and distribution shrank in the cold LGM (~22 Ka) (Figure 3d,f,g,i). In addition, the expansion of Q. acutissima began at 0.66 Ma and the further expansion in the range of this species occurred around 0.42 Ma (Table S4). This period comprises several interglacial intervals in the middle Pleistocene, which is an extralong interglacial period may have persisted throughout the Marine isotope stages (MISs) 15-13 (0.62-0.48 Ma) [81]. It was reported that southern Loess Plateau was likely under a subtropical semi-humid climate in MIS 13, with an annual mean temperature at least 4-6 • C higher and annual rainfall 200-300 mm higher than at present [23,82]. Moreover, the temperature and precipitation are especially important for plant species like other Quercus species that dwells in the northern region of the subtropical China [24,76]. Recently, Zhang et al. (2018) reported that an intraspecific admixture event of Q. acutissima is the causes for expansion from refugia into the mountainous area of central China occurred about 0.69 Ma, implying that the warm and humid climate triggered the genetic admixture along with the duration of recolonization [23]. Moreover, the ENM results suggested that this species was also predicated to have undergone a range contraction during the last glacial periods (Figure 3a,c,e). Overall, our findings highlight that the importance of climatic oscillations at the Pleistocene, had a remarkable effect on the distribution of the three Cerris oaks and supported the 'contraction-expansion' scenario usually reported in other deciduous oak species (e.g., [75,83]) in China.
In addition, our present study suggested that Q. acutissima and Q. variabilis were likely to expand from multiple glacial refugia after LGM (~22 Ka), in Qinling mountains and Yungui Plateau (Figure 3a,c), while Q. chenii survived in the mountainous areas in East China (Figure 3e). In fact, mountains have often served as refugia for many species distributed in the temperate and tropical regions of the world, because the diverse topography of mountains may offer a wide array of climatic and edaphic conditions, even in the harshest periods (e.g., glacial periods of the Pleistocene) [84]. Nuclear microsatellite data in our study revealed several populations (SPU/E, SPg/J, SPWS) located in these regions possessed high genetic diversity (Table S1). Besides, Palaeovegetation reconstructions of the LGM period supported a hypothesis that trees might have similar refugia in the cold periods [21,22]. Hence, our inference strengthened this view and confirmed that tree species might have persisted in multiple refugium located at subtropical and warm temperature zones in China. Strikingly, such tree species have been isolated for a long time among multiple habitats within the mountains area; these aspects might further promote the genetic differentiation among lineages. Meanwhile, our results also indicated that a probable scenario of interspecific introgression due to the shared refugium and/or secondary contact when expanded from the refugium in the warm periods (Figure 3a,c,e).

Conclusions
The results of our study on clearly indicated that evolution and demography of three East Asia Cerris oak species (Q. acutissima, Q. variabilis and Q. chenii) was strongly influenced by the topography change and climate oscillation during the Miocene-Pliocene and Pleistocene. Furthermore, like other deciduous oak species in Eastern Asia (e.g., Q. fabri and Q. liaotungensis [75,83]), these three Cerris oak species in the Eastern Asia share a similar performance of 'expansion-contraction' in response to past climatic oscillations.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/f12091164/s1, Supporting Tables: Table S1. Genetic diversity of 44 populations of three Cerris oaks (Q. variabilis, Q. acutissima and Q. chenii) estimated based on 16 microsatellite loci; Table S2. Primer sequences and characterization of the 16 microsatellite loci used in this study; Table S3. Prior distributions for model parameters used in divergence model comparisons; Table S4. Posterior probability under logistic and descriptions of four scenarios, and prior setting for all parameters used in approximate Bayesian computation (ABC) for demographic analysis of three species (Figure 3b); Table S5. The 14 bioclimatic variables that were used in Ecological niche modeling; Table S6. Diversity indices of 16 EST-SSR loci for 759 individuals; Table S7. Diversity indices of 16 EST-SSR loci for three related species; Table S8 Hardy-Weinberg equilibrium (HWE) test in this study; Table S9. Estimation of rates of contemporary gene flow per generation among three species (QV, QA and QC) using the programs BAYESASS v3.0; Table S10. Posterior probability of each scenario and 95% confidence intervals (CI) based on the logistic regression approach for approximate Bayesian computation analyses (ABC) for all three related species (Q. acutissima, Q. variabilis and Q. chenii). Type I and type II errors for the best supported scenario (in bold) are indicated; Table S11. Indices of bottleneck test of the three related species. Supporting figures; Figure S1. Bayesian inference analysis of 16 EST-SSR data determining the most likely number of clusters (K) for three Cerris oaks (Q. variabili, Q. acutissima and Q. chenii). The distribution of the likelihood L(K) values (a) and DeltaK values (b) are presented for K = 1-8 (10 replicates); Figure S2. Bayesian inference analysis of nuclear data for Q. variabili and Q. acutissima separately, STRUCTURE plots are presented for K = 2.; Figure S3. Comparison of posterior probabilities for four simulated scenarios of species differentiation obtained by logistic regression from 1% of the closest data set in DIYABC; Figure S4 Principal component analysis results obtained along the first three axes for the pre-evaluate models and prior distributions of three oak species divergence in DIYABC; Figure S5. Model checking in DIYABC for the best-supported (species divergence: model 2); Figure S6. Principal component analysis results obtained along the first three axes for the pre-evaluate models and prior distributions of three oak species demographic models in DIYABC. (a-c) representing the species of Q. variabilis, Q. acutissima, and Q. chenii, respectively; Figure S7. Model checking in DIYABC for the best-supported demographic models; Figure Figure S10. Contribution of each environmental variable to spatial distribution of the PCA-env.