We sampled leaf material from 996 individuals in 60 sites of Quercus aquifolioides and 661 individuals in 46 sites of Quercus spinosa throughout the species distribution range. The study sites were at least 30 km apart, and individuals were at least 100 m apart from each other to avoid sampling clone individuals. All leaf materials were rapidly dried in silica gel and stored for DNA isolation. The detailed information on sampling sites is depicted in Figure 1 and Supplementary Table 1.

Total genomic DNA was extracted from leaf samples for each individual using an improved cetyltrimethylammonium bromide (CTAB) method (Richards et al., 1994). We randomly selected one individual from each of six distant sites for pre-amplification experiments with 25 nuclear microsatellite (nSSR) loci developed for other oak species (Dow et al., 1995; Steinkellner et al., 1997; Kampfer et al., 1998; Ueno et al., 2008; Durand et al., 2010; Supplementary Table 2). We excluded loci harboring null alleles using MICRO-CHECKER 2.2.3 (Van Oosterhout et al., 2004). Departure from Hardy–Weinberg equilibrium (HWE) and linkage disequilibrium (LD) was evaluated using GenALEx 6 (Peakall and Smouse, 2006) and FDIST2 (Beaumont and Nichols, 1996). Finally, fifteen successfully amplified SSR loci were retained for subsequent analyses. The reaction procedures are modified from Du et al. (2017). The allele sizes were subsequently scored using GeneMarker v. 2.2 (Softgenetics, United States), and the genotypes were checked visually two times. A subset of the data, 959 individuals from 58 study sites of Q. aquifolioides at 15 nSSRs, were from Du et al. (2017) and Li Y. et al. (2021), and the additional data were first reported in this study.

We estimated genetic diversity indices including mean observed heterozygosity (H_O), mean expected heterozygosity (H_E), mean unbiased expected heterozygosity (uH_E), mean effective population size (N_E), and mean Shannon index (I) by GenAlEx 6 (Peakall and Smouse, 2006). The significance of genetic diversity was evaluated by t-test in SPSS 22 (SPSS Inc., Chicago, IL, United States) with a significance level of 0.05. In order to formulate optimal conservation strategies by revealing priority areas for conservation, we applied spatial analysis to improve the understanding of the geographic distribution of genetic diversity across the oak distribution range. We calculated and mapped the AR and LCA according to the circular neighborhood methodology described by Van Zonneveld et al. (2012). AR, also referred to mean number of alleles per locus, is a straightforward measure of genetic diversity based on molecular markers that aim at selecting populations for conservation (Frankel et al., 1995; Petit et al., 1998). LCA are alleles that occur in 25% or less of all grid cells and with a frequency of at least 5% in a grid cell per locus. Population with high LCA indicate the presence of genotypes adapted to specific environments; therefore, priority for conservation should be given to those populations (Frankel et al., 1995). After applying circular neighborhood to all samples, we calculated the AR and LCA for all 10-minute grid cells by GenAlEx 6 (Peakall and Smouse, 2006). AR was corrected by rarefaction to a minimum sample size of 10 re-sampled trees per cell with the HP-RARE software (Kalinowski, 2005).

We examined the genetic differentiation using hierarchical analysis of molecular variance (AMOVA, Excoffier et al., 1992) in Arlequin 3.5 (Excoffier and Lischer, 2010). The significance of fixation indices was tested using 10,000 permutations in Arlequin 3.5. We used a model-based clustering program implemented in STRUCTURE 2.3 (Pritchard et al., 2000) to infer the genetic clustering without consideration of sampling information. The program was run with the number of clusters (K) varied from 1 to 10 with 20 independent replicates conducted for each K-value, and the length of the burn-in period was set to 100,000 steps followed by the number of Markov chain Monte Carlo (MCMC) after burn-in of 100,000. We selected the optimal K-value by ΔK statistics performed in the web-based program STRUCTURE HARVESTER (Earl and Vonholdt, 2012). Graphic visualization of the STRUCTURE results was produced using DISTRUCT 1.1 (Rosenberg, 2004). We also conducted a principal component analysis (PCA) to visualize the genetic relatedness among individuals by calculating principal components (Novembre and Stephens, 2008) using “adegenet” R package (Jombart and Ahmed, 2011). The first two eigenvectors were plotted, and the discrete points reflect the real structure of populations. In addition, we conducted a principal coordinate analysis (PCoA, Gower, 1966) based on genetic covariance among populations in GenAlEx 6 (Peakall and Smouse, 2006) and plotted the first two eigenvectors to visualize genetic relatedness.

The historical gene flow of two oak species was assessed by Migrate-n 3.6 (Beerli and Felsenstein, 2001; Beerli, 2006) based on the Bayes factor value. First, we generated initial θ (4Neμ, four times effective population size multiplied by mutation rate per site per generation) and M (immigration rate divided by the mutation rate) to estimate the amount and direction of gene flow. A continuous Brownian motion model and the default genetic differentiation were used to generate initial theta and migration values. Then, we started three independent MCMC chains with 500,000 iterations, respectively. We sampled every 100 steps under a constant mutation model and discarded the first 10,000 records as burn-in. After checking the model convergence, we calculated the mode value and 95% posterior probability.

We obtained climatic variables of the current conditions (∼1970–2000) from WorldClim¹, a database of high spatial resolution global weather and climate data (Fick and Hijmans, 2017). A total of 31 climatic variables, including the full suite of 19 mean annual bioclimatic variables and 12 average monthly climate data for precipitation, were downloaded. We excluded climatic variables that were highly correlated with the threshold values of 0.7 using a variance inflation factor (VIF) test in “usdm” R package (Naimi et al., 2014). After avoiding the high multicollinearity bias, four climatic variables, namely precipitation seasonality (bio15, coefficient of variation), mean temperature of the driest quarter (bio09), temperature annual range (bio07, between the minimum temperature of the coldest month and the maximum temperature of the warmest month), and precipitation during June (prec06), were finally remained for downstream analyses (Supplementary Table 1).

The linear relationships analysis can integrate environmental variables and spatial genetic structure into the analytical framework to assess the contributions of geography and environment in driving genetic differentiation (Feng and Du, 2022). The loading results of this analysis can be interpreted as the response proportion of environmental factors to genetic variation. In this study, we performed isolation-by-resistance (IBR) to illustrate the effects of the heterogeneous landscapes on the population genetic connectivity of two oak species in “ade4” R package (Dray and Dufour, 2007). We first predicted the potential distribution of two oak species based on the current ecological niche model (ENM) in MAXENT (Phillips and Dudik, 2008) and then transformed the environmental rasters into resistance surfaces. We generated the resistance distance based on circuit theory in CIRCUITSCAPE 4.0.5 (McRae, 2006; McRae et al., 2008) and “ResistanceGA” R package (Peterman, 2018). We performed Mantel tests of isolation by distance (IBD; Van Strien et al., 2015) and isolation by environment (IBE; Manthey and Moyle, 2015) to test the linear relationships between geographic or environmental distance and genetic distance using “ecodist” R package (Goslee and Urban, 2007). To distinguish the impact of IBD and IBE, a partial Mantel test was used to evaluate IBE/IBD by controlling the linear influence of geographic/environmental distance (Smouse et al., 1986). In addition, we performed multiple regression on distance matrices (MRM, Lichstein, 2007) to test the multivariate correlation between genetic distance matrix and climate distance using “ecodist” R package (Goslee and Urban, 2007). The significance for Mantel tests and MRM was evaluated by 10,000 permutation tests with the significance level set to 0.05.

We performed redundancy analyses (RDAs) to detect the multivariate relationship between genetic variation and climate variation (Van den Wollenberg, 1977; Legendre and Legendre, 1998) using “vegan” R package (Oksanen et al., 2017). A partial redundancy analysis (pRDAs, Legendre and Legendre, 1998, 2012) was performed to avoid the linear influence of geographic/climate variables when analyzing the climate/geographic variables. Statistical significance was evaluated from 999 permutations.

A limitation of the linear relationships analysis is the inability to fit the variation in the rate of compositional turnover along environmental gradients and the curvilinear relationship between genetic distance and environmental and geographic distance. Therefore, a non-linear relationship is essential for applying the associated turnover function to each mapped environmental variable (Fitzpatrick and Keller, 2015). In this study, we performed generalized dissimilarity modeling (GDM) to identify non-linear relationships between genetic distance matrix (response variable) and geographic/environmental distances (predictors) using “gdm” package (Ferrier, 2002; Ferrier et al., 2007). We also evaluated the variation in the rate of allelic compositional change along environmental gradients by fitting splines (Fitzpatrick and Keller, 2015). Genetic distances among individuals were calculated based on allele frequency, and geographic distance was based on Euclidean distance among coordinates. We assessed the variable significance by randomization tests and assessed uncertainty due to sampling error by simulating 1,000 bootstrap iterations (Ferrier et al., 2007; Fitzpatrick et al., 2013).

We found that the genetic diversity was higher in Q. aquifolioides than in Q. spinosa (H_O: 0.59 vs. 0.41; H_E: 0.58 vs. 0.49; uH_E: 0.61 vs. 0.52; P < 0.01) (Supplementary Table 3). We also identified significantly higher genetic diversity in Q. aqu-East than Q. aqu-West lineage (H_O: 0.60 vs. 0.53; H_E: 0.61 vs. 0.52; uH_E: 0.64 vs. 0.54; P < 0.01) and a slightly higher genetic diversity in Q. spi-East lineage than Q. spi-West lineage (H_O: 0.42 vs. 0.41; H_E: 0.51 vs. 0.47; uH_E: 0.53 vs. 0.50; P < 0.01) (Supplementary Table 3).

We applied a circular neighborhood re-sampling technique to ensure sufficiently and more evenly distributed data points for spatial diversity analysis. A total dataset of 31,872 trees for Q. aquifolioides and 21,152 trees for Q. spinosa was used for further AR and LCA analyses (Supplementary Figure 1). Our results showed that the enriched regions of AR and LCA of two oak species were different (Figure 2). For Q. aquifolioides, the populations located at HDM (Q. aqu-East lineage) in southwest Sichuan province and northwest Yunnan province contained higher AR and LCA than QTP (Q. aqu-West lineage) (Figures 2a,c). The marginal population PW of Q. aquifolioides located at the easternmost end of the HDM with a lower AR than in other areas of the HDM (Figure 2a). For Q. spinosa, populations from QM (Q. spi-East lineage) revealed higher AR and LCA than HDM (Figures 2b,d).

Bayesian clustering identified K = 2 as the optimal number of evolutionary clusters (Supplementary Figure 2), subdivided all individuals into two clusters, one corresponded to Q. aquifolioides and the other to Q. spinosa. When K = 3, Q. aquifolioides was maintained unchanged while Q. spinosa is further subdivided into two geographically related lineages: Q. spi-West (22 sites from HDM) and Q. spi-East lineage (24 sites from QM and lowlands in East China). When K = 4, Q. aquifolioides is divided into Q. aqu-West (17 sites from QTP) and Q. aqu-East lineage (43 sites from HDM) (Figure 3). The results of PCA and PCoA were largely consistent with the STRUCTURE analysis with clear separation in interspecific and intraspecific levels (Figure 4). AMOVA showed a high level of genetic differentiation between Q. aquifolioides and Q. spinosa, and most of the variation occurred within populations (F_ST = 0.35, 65.4%) (Table 2).

The Migrate-n analysis generated θ and M values greater than zero, which revealed an asymmetric historical gene flow between two species, mainly from Q. aquifolioides to Q. spinosa (56.1 vs. 47.5) (Table 3). Moreover, we found gene movements occurred predominantly from Q. aqu-West into Q. aqu-East lineage (46.6 vs. 34.5) and from Q. spi-West into Q. spi-East lineage (65.6 vs. 43.0) (Table 3).

Our analyses revealed highly significant correlations between pairwise genetic distances and resistance distance in Q. aquifolioides and Q. spinosa, but not in Q. aqu-West lineage (P = 0.745). The pattern of IBR in Q. spi-East was stronger than in Q. spi-West lineage (R² = 0.65 vs. R² = 0.42; Figure 5). In addition, we found significant patterns of IBD and IBE by Mantel and partial Mantel tests in Q. aquifolioides and Q. spinosa (Table 4 and Supplementary Figure 3). However, significant IBD was only detected in Q. aqu-West and Q. spi-East lineage; significant IBE was detected in Q. aqu-East and Q. spi-West lineage (Table 4 and Supplementary Figure 3). More specifically, genetic distance was significantly associated with annual range temperature (bio07) and seasonal precipitation (bio15) in Q. aquifolioides and Q. aqu-East lineage; bio15 and precipitation during June (prec06) in Q. spinosa and Q. spi-West lineage (Supplementary Table 4). These findings were consistent with optimal MRM models that yielded similar results (Supplementary Table 5).

The percentages of variance explained by RDA and pRDA were similar, and we thus report results for pRDA. Geography (4.0 and 5.2%) explained more genetic variance than climate variables (1.5 and 3.0%) in Q. aqu-West and Q. spi-East lineage, whereas climate variables (2.7 and 5.6%) explained more genetic variance than geography (1.4 and 5.1%) in Q. aqu-East and Q. spi-West lineage (Table 5). Partitioning of the total genetic variance revealed that bio07 and prec06 explained most genetic variance in Q. aqu-East lineage (44.5 and 23.5%), while in Q. spi-West lineage, prec06 and bio15 were the two most explanatory environmental variables (46.0% and 20.3%) (Table 5 and Supplementary Figure 4).

Generalized dissimilarity modeling analyses suggested that geography was the most important predictor among all variables considered in Q. aquifolioides, Q. aqu-West, Q. spinosa, and Q. spi-East lineage (59.3, 62.9, 53.4, and 8.8%) (Supplementary Table 6). However, there was almost no contribution of geography in Q. aqu-East and Q. spi-West lineage (2.8 and 7.2%), while bio15 was the most important environmental factor (36.2 and 30.8%) (Figure 6 and Supplementary Table 6). These results were consistent with I-spline analysis (Supplementary Figure 5).

Quercus aquifolioides and Quercus spinosa were assigned to two distinct genetic clusters (Figures 3, 4), and the direction of the interspecific and intraspecific gene flow was from west to east (Table 3). These results showed a strong West-to-East migration pattern and likely reflect long-term geographic isolation due to the orogenic history of QTP and adjacent regions (Qiu et al., 2011; Wen et al., 2014). In addition, the amount of historical gene flow was asymmetric, mainly from Q. aquifolioides to Q. spinosa at the interspecific level (56.1 vs. 47.5) and West lineages to East lineages at the intraspecific level (46.6 vs. 34.5 and 65.6 vs. 43.0) indicating the asymmetric introgression of the species (Currat et al., 2008; Du et al., 2011).

There are already several migration routes for the out-of-QTP hypothesis (Table 1), and all of the studies suggested that migration, orographic, and climate oscillations catalyzed intraspecific differentiation, diversification, and adaptation of species in this region (see the summary of Wen et al., 2014). It is suggested that orographic and climatic oscillations might result in lots of small fragmented habitats with different microclimates, which could influence the direction of natural selection, and might promote intraspecific high differentiation of species (Sobel et al., 2010). Our results represent a very typical case of a West-to-East migration pattern, which might be triggered by extensive uplifts of the QTP (see the summary of Favre et al., 2015). The QTP uplift events provided opportunities for the ancestral population in this region continually expanded to its eastward ranges and gradually triggered and facilitated speciation and diversifications of oak species (Zhou, 1992). Meanwhile, the West-to-East migration pattern indicated that migrating to new favorable locations might be a survival strategy of species to rapid climate changes as in Sophora davidii (Fan et al., 2013) and Gentiana (Favre et al., 2016).

Genotype–environment association (GEA) analysis, including mantel tests, redundancy analyses, and generalized dissimilarity modeling, integrates environmental variables and spatial genetic structure into the analytical framework to detect the adaptive variation (Feng and Du, 2022). GEA analysis is essential for understanding the mechanisms underlying the evolutionary responses to environments and was used to quantify patterns of interaction between genetic variation and climate conditions (Hansen et al., 2012). In addition, resistance analysis is important to understand how the species respond to different resistance distances (McRae, 2006). Based on this theory, we identified the landscape resistance matrix that was most highly correlated with pairwise genetic distances in Q. aquifolioides and Q. spinosa, especially in Q. aqu-East and Q. spi-East lineage (Figure 5), and this result may be related to increased habitat isolation in this area resulted from the disjunct distribution of the two oak species. However, the pairwise genetic distance was not correlated with resistance distance in Q. aqu-West lineage (Figure 5), indicating that populations from this lineage are low resistant to dispersal and might endure high genetic connectivity among populations. This result is confirmed with the lowest genetic differentiation in the Q. aqu-West lineage (Table 2).

The IBD and IBE results indicated the intraspecific lineages of oaks with different response patterns of genetic variation. We detected significant IBD patterns in Q. aqu-West and Q. spi-East lineage (Table 4 and Supplementary Figure 3, Supplementary Table 4). These results were consistent with GDM and RDA (Figure 6, Table 5 and Supplementary Figure 4, Supplementary Table 6), indicating that the genetic variation of Q. aqu-West and Q. spi-East lineage was mainly driven by selectively neutral evolutionary processes, not by strong selection pressure from the environment. The complex geological structure of mountains might form a natural geographic barrier for seeds or pollen dispersal, which can provide potential conditions for the formation and independent evolution of plants’ intraspecific lineages (Liu et al., 2013; Li et al., 2014).

By contrast, we detected significant IBE patterns in Q. aqu-East and Q. spi-West lineage (Table 4 and Supplementary Figure 5, Supplementary Tables 5, 6), where the extreme environmental conditions on the plateau might be regarded as a significant climatic barrier, rather than a geographic barrier. It also suggests that geographic isolation may cause interspecific and intraspecific differentiation; adaptation to local climate and environmental factors reinforces this differentiation and gradually forms this significant IBE pattern (Gao et al., 2021). Accordingly, GDM and RDA both suggesting temperature annual range (bio07) and precipitation during June (prec06) were the most significant environmental factor driving genetic variation in Q. aqu-East and Q. spi-West lineage (Figure 6, Table 5 and Supplementary Figure 4, Supplementary Table 6), respectively. The temperature may be the main driver of genetic variation for Q. aqu-East lineage, and it may influence the growth of plants in microhabitats by affecting the metabolic processes (Wahid et al., 2007). Precipitation might have a great impact on phenological and growth of oak from Q. spi-West lineage and then affect the ability to adapt to climate change.

A better understanding of the spatial distribution of genetic diversity is necessary for the formulation of effective and efficient conservation strategies (Petit et al., 1998). Priority for conservation should be given to populations that retain the highest AR and LCA because the likelihood to find interesting breeding materials is higher in the highest genetic diversity populations, which can indicate the presence of genotypes adapted to specific environments (Frankel et al., 1995; Tanksley and McCouch, 1997). We found that the priority conservation areas were different for the two species based on a large number of samples (1,657 individuals) across their distribution range (Figure 2). For Q. aquifolioides, a priority area for conservation should be the populations located at HDM (Q. aqu-East lineage), which contained the highest AR and LCA. In addition, the marginal population PW of Q. aquifolioides located at the easternmost end of the HDM with a lower AR and LCA than other areas. Risk of non-adaptedness (RONA) revealed that this marginal population might be at higher risk of extinction under future climate (Du et al., 2020). Therefore, populations isolated at the edge of the distribution range also should be considered in conservation activities to prevent the extinction of species in this area. For Q. spinosa, a priority area for conservation should be the populations located at QM (Q. spi-East lineage) with the highest genetic diversity. The second area of a higher diversity of Q. spinosa is located on the border zone between Sichuan and Yunnan province, probably related to the high species richness of this region. The higher AR and LCA in the Q. aqu-East and Q. spi-East lineage than in Q. aqu-West and Q. spi-West lineage suggested those populations received more immigrants and played an important role in evolution and diversification (López-Pujol et al., 2011; Favre et al., 2015), while the lower genetic diversity in Q. aqu-West and Q. spi-West lineage likely reflects population contraction or extinction–recolonization dynamics in this area (Qiu et al., 2011; Du et al., 2017; Meng et al., 2017).