Sampling of oaks were granted and supported by National Forestry Bureau of China, Local National Nature Reserves, and Ministry of Environmental Conservation and Forestry, Myanmar.

Thirty-three Q. franchetii complex populations were sampled from Yunnan and Sichuan, China and Mount Victoria in Chin State, Myanmar. In total, 303 individuals from 33 populations were sampled for this study, covering the major known distribution range of the Q. franchetii complex (Table 1). Samples from the same population came from individuals that were separated by at least 50 m. At least 10 trees were sampled from each population, except for populations with very few individuals (WD, EY, YJ, XP, XN), in which case we sampled all the accessible adult trees in those populations. Fresh and healthy mature leaves were collected and put into containers with silica gel to dry them quickly until DNA extractions could be performed. The voucher specimens of the DNA samples were deposited in the herbarium of Shanghai Chenshan Botanical Garden (CSH).

Total genomic DNA was extracted using a modified cetyltrimethyl ammonium bromide (CTAB) protocol (Doyle and Doyle, 1987). Three pairs of cpDNA primers, namely psbA-trnH (Shaw et al., 2005), trnT-trnL (Taberlet et al., 1991), and atpI-atpH (Grivet et al., 2001), and eight highly polymorphic nuclear microsatellite loci, specifically nuclear simple sequence repeats (nSSRs), were selected for genotyping on all samples. Primer sequences and PCR amplification conditions are summarized in Supplementary Appendix 1. The cpDNA and nSSR amplification conditions followed the methods described by Xu et al. (2015) and An et al. (2016), respectively. PCR products of nSSR markers were 10 times diluted, then mixed with fluorescence size standards at a ratio of 6-FAM: HEX: ROX = 1:1:2, then genotyped by Shanghai Majorbio Bio-pharm Technology Co., Ltd. (Shanghai, China). The PCR products of cpDNAs were purified, then biodirectionally sequenced by the same company. All sequences obtained in this study have been deposited in GenBank (see “Data Availability”).

The cpDNA sequences were assembled and checked using Sequencher 4.01 (Gene Codes Corp., Ann Arbor, MI, United States). The ClustalW implementation in MEGA X (Kumar et al., 2018) was used for sequence alignment with manual adjustment. All microsatellite loci were checked using GeneMarker^® (Hulce et al., 2011) to detect and analyze allele sizes. Null alleles and stutter bands were checked with MicroChecker (Van Oosterhout et al., 2004).

Haplotypes and polymorphism statistics for cpDNA loci were calculated with DnaSP 6.0 (Rozas et al., 2017). The haplotype geographic distribution was projected onto a map using ArcGIS 10.5.¹ Total haplotype diversity (H_T), within-population diversity (Hs), and coefficients of differentiation (G_ST and N_ST) for cpDNA loci were estimated using PERMUT 2.0 (Pons and Petit, 1996). Haplotype diversity (h) and nucleotide diversity (π) for cpDNA loci were obtained with ARLEQUIN 3.5 (Excoffier and Lischer, 2010). The cpDNA haplotype network was constructed using the Median-Joining model in NETWORK 5.0.0.1 (Bandelt et al., 1999). GenAlEx 6.5 (Peakall and Smouse, 2012) was used to calculate the genetic diversity index of microsatellite data, including expected heterozygosity (H_e), observed heterozygosity (H_o), number of alleles (N_A), and effective number of alleles (N_E). Allelic richness (A_r) of nSSRs was determined using HP-RARE (Kalinowski, 2005).

The genetic differentiation values (F_ST) based on cpDNA and nSSR data were estimated using ARLEQUIN 3.5, respectively (Excoffier and Lischer, 2010). The Genetic Landscape GIS toolbox (Vandergast et al., 2011) in ArcGIS 10.5 was used to generate a geographical landscape map based on both genetic diversity (A_r) for cpDNA loci and genetic divergence (F_ST) based on cpDNA and nSSR data according to inverse distance weighted interpolation. For cpDNA and nSSR data, analysis of molecular variance (AMOVA) was performed to assess the genetic variation among populations and within populations. The Mantel test was also performed based on the genetic structure and geographic distance matrix with 1,000 random permutations to evaluate their relationship and test isolation by distance (IBD). For nSSR data, Principal coordinate analysis (PcoA) based on genetic distance was performed using GenAlex 6.5 (Peakall and Smouse, 2012) to assess differences among individuals or groups. Barrier (Manni et al., 2004) was used to set up geographic barriers according to sample locations to detect the existence of genetic barriers among populations.

We used GenePop v 4.2 (Rousset, 2008) to test departures from Hardy-Weinberg equilibrium (HWE) for each of the eight nSSR loci. As the nSSR loci significantly deviated from HWE (see Results section “Genetic Diversity and Structure”), Bayesian assignment probability methods using the programs InStruct (Gao H. et al., 2007) and STRUCTURE 2.3.4 (Pritchard et al., 2000) were both used to infer the population genetic structure. The number of clusters (K) was varied from 1 to 10, with 10 replicates at each value. Each run consisted of a burn-in length of 25,000 iterations with a run length of 500,000 MCMC (Markov chain Monte Carlo) iterations. The optimal K (number of clusters) was determined using the △K method (Evanno et al., 2005). CLUMPP (Jakobsson and Rosenberg, 2007) was used to align 10 runs of InStruct with the optimum K using a greedy algorithm.

Among the three models in the Bottleneck program (Luikart et al., 1998), the “Two-phase mutation model” (TPM) is the most suitable for microsatellite loci, which was thus selected to test whether the populations of the Q. franchetii complex had experienced a bottleneck. We used “Wilcoxon sign-rank test” method for a significance test.

The divergence time dating for cpDNA haplotypes of the Q. franchetii complex was estimated using a Bayesian approach as implemented in BEAST V2.4 (Suchard et al., 2018). Quercus glauca (Quercus section Cyclobalanopsis) was chosen as the outgroup to root the tree. An uncorrelated lognormal relaxed clock was applied with the K81uf + I substitution model, which was selected based on the Akaike information criterion (AIC) in Modeltest 3.7 (Posada and Crandall, 1998). The earliest conclusive leaf fossils of Quercus section Ilex were discovered in Tibet, dated to the late Eocene (ca. 34 Ma) (Su et al., 2019) was set as the minimum age to constrain the stem of the haplotype tree of the Q. franchetii complex with a lognormal distribution and a median of 37.8 Ma (95% HPD: 34.02–56.6 Ma). The MCMC chains were run for 100 million generations with a sampling frequency of once every 10,000 generations. Convergence was assessed using Tracer v1.7² (Rambaut et al., 2018), and the effective sample sizes for all parameters were calculated. The resulting tree and log files from the two replicate runs were combined with LogCombiner v1.8. Then, we used TreeAnnotator v. 1.8³ to generate the maximum clade credibility (MCC) tree after discarding the first 20% of the trees as burn-in. The results were visualized using FigTree v1.4.3.⁴

Pairwise mismatch distribution analysis for cpDNA loci, Tajima’s (1989) D and Fu’s (1997) Fs of neutrality tests were performed to detect possible demographic expansions of the Quercus franchetii complex using ARLEQUIN 3.5 (Excoffier and Lischer, 2010).

Ancestral range reconstruction was analyzed using statistical dispersal-vicariance (S-DIVA) analysis as implemented in RASP (Yu et al., 2015) and DEC models (Ree et al., 2005; Ree and Smith, 2008). Four areas were delimited for ancestral area reconstruction based on the geological characteristics and biogeographical division of southwestern China, including (A) the Nanpan River region (NPR); (B) the southwestern Red River (including southern Himalayas to Red River) (RR); (C) the Hengduan Mountains (HDM), and (D) the Yunnan-Guizhou Plateau (YGP) (Figure 2A).

We used MaxEnt 3.4 (Phillips et al., 2006) to simulate the potential distribution range of the Q. franchetii complex (including Q. franchetii and Q. lanata) under the past, current, and projected future climate scenarios based on the maximum entropy model (Phillips et al., 2004). In total, 145 accurate occurrence points were collected from the Global Biodiversity Information Facility (GBIF),⁵ Chinese Digital Herbarium (CVH),⁶ and the field collection records of our research team. Each voucher specimen in the distribution record was inspected and checked carefully. Nineteen bioclimatic variables with a 2.5-arc-min resolution for the present (1950–2000), the last glacial maximum (LGM) (CCSM) period, and the future (2060–2080, RCP2.6; RCP4.5; RCP6.0; RCP8.5) were downloaded from WorldClim 2.⁷ The nine environmental factors (bio1, Annual Mean Temperature; bio4, Temperature Seasonality; bio6, Min Temperature of Coldest Month; bio7, Temperature Annual Range; bio10, Mean Temperature of Warmest Quarter; bio11, Mean Temperature of Coldest Quarter; bio12, Annual Precipitation; bio13, Precipitation of Wettest Month; bio17, Precipitation of Driest Quarter) were selected after eliminating climatic variables such that none that were include were highly correlated (i.e., with correlation coefficients greater than 0.8) using the Dismo package in R.⁸ We also chose CHELSA climate data,⁹ which had a higher prediction power in mountain regions (specifically the Himalayas), to simulate the present distribution (Bobrowski and Schickhoff, 2017; Karger et al., 2017).

In order to improve the prediction accuracy of the model, it was necessary to optimize the model in MaxEnt 3.4 (Phillips et al., 2006) by setting the β multiplier and environmental characteristic parameters. The MaxEnt model captures five features: linear (L), quadratic (Q), hinge (H), product (P), and threshold (T). In this research, we used seven feature combinations (auto, L, H, LQ, LPQ, LQH, LQHP) and set the regularization multiplier from 0.5 to 10 with increments of 0.5. Then, we used ENMTools (Warren et al., 2010) to calculate the lambdas file for the maxent result, and selected the model with the smallest AIC value as the optimal model parameter for further analyses. The occurrence points of Q. franchetii and Q. lanata were randomly divided into 75 and 25% of the data for training and testing, respectively. The area under the curve (AUC) was used to evaluate the model’s accuracy. AUC values ranged from 0.5 to 1, where the higher AUC value indicated a better prediction. The maximum training sensitivity plus specificity thresholds for the presence or absence of species was used to draw a species distribution map in ArcGIS 10.5. In addition, we classified three categories equally between this threshold and 1, corresponding to low, medium, and high fitness areas, respectively. To compare the changes in species distribution across different periods, three indicators were calculated: locality habitat stability (N_Stab), habitat distribution area ratio (N_a), and habitat expansion extent (N_e). These values were calculated using the following formulas: N_Stab = 1- | N_Pre - N_LGM |, where N_Pre and N_LGM are the habitat suitability of the present and LGM distribution area; N_a = (present distribution areas)/(LGM or future distribution areas): a value close to 1 indicates a stable distribution of the species, while a value much higher or lower than 1 indicates that the distribution area of the species has expanded or contracted from the LGM to the present; N_e = [1 - (Distribution area overlapping between the LGM and present or the future and present/present distributions area)] × 100% represents the percentage of the distribution that has expanded from the LGM to the present.

The potential dispersal routes of the Q. franchetii complex in the past and present periods were inferred based on the least-cost path analysis method using SDM toolbox 2.0 (Brown, 2014) in ArcGIS 10.5. The specific steps are listed below: Firstly, we generated a resistance layer by inverting the SDMs (1-SDM). The resistance layer was used to create a cost distance raster for each sample locality. The corridor layers were built between two locations that only share haplotypes based on the cost distance raster. We used the categorical the least cost path (LCP) approach to better describe the habitat heterogeneity and its role in the dispersal. The value of each corridor layer was divided into low, medium, and high, and then these three intervals were re-divided into new values 5, 2, and 1. Finally, we reclassified all the corridor layers, summarized and standardized them from 0 to 1, and determined the dispersal corridors of the Q. franchetii complex.

The linear model in R 3.5¹⁰ was used to estimate the correlation of genetic diversity indexes (A_r and H_e) and genetic structure (cluster A from the Bayesian clustering, C_A) with habitat and geographic factors. Five variables—population longitude and latitude, habitat suitability for the present (N_Pre) and the LGM (N_LGM), and habitat stability (N_Stab) were used as explanatory covariates.

The psbA-trnH, atpI-atpH, and trnT-trnL sequence alignments were 482–520, 899–1030, and 796–838 bp in length, respectively. The combined length of the three aligned chloroplast fragments was 2,319 bp, with 94 polymorphisms across a total of 33 haplotypes among the 303 individuals analyzed. The cpDNA haplotype diversity was h = 0.956 and nucleotide diversity was π = 0.00536. The nucleotide diversity and haplotype diversity within the populations were 0–0.04184 and 0–0.644, respectively (Table 1 and Figure 1). The total diversity (H_T = 0.982) was much higher than the average within-population diversity (H_S = 0.123). The genetic diversity map showed that the NPR and HDM populations had high genetic diversity (Figure 3C). N_ST (0.959) was significantly greater than G_ST (0.874) (P < 0.05), which indicated that a clear phylogeographic structure existed among the Q. franchetii complex populations. The network analysis resolved four distinct cpDNA haplotype clades in the Q. franchetii complex, which was consistent with the haplotype structure determined by BEAST. The haplotypes in Clade IV exhibited a star-like structure. The most common haplotype, H9, had the widest distribution, mainly within the central YGP (Figure 2B).

Among the 303 individuals genotyped, a total of 115 alleles were identified using eight pairs of microsatellite primers. All eight microsatellite loci were highly polymorphic, with allele numbers varying from 8 to 18 per locus. The genetic diversity indexes N_A, N_E, H_e, and H_o were 4.182 (SE = 0.103), 2.578 (SE = 0.066), 0.537 (SE = 0.013), and 0.536 (SE = 0.016), respectively, while allele richness (A_r) was 3.18–4.34.

The AMOVA on cpDNA sequence data revealed greater genetic variation among populations (87%; F_ST = 0.87) than within populations (13%). In contrast, the nSSRs showed substantial genetic differences within populations (63%; F_ST = 0.369, Table 2). Pairwise F_ST values calculated based on chloroplast and SSR data were shown in Supplementary Appendix 2. In order to more visually observe genetic differentiation, F_ST values were projected onto a map and a genetic divergence map was made. The genetic divergence map showed that the population in southwestern Yunnan had a relatively high genetic divergence, while the divergence of the populations from the northeast was comparatively lower (Figures 3A,B). The Mantel tests of cpDNA data (r = 0.052, P = 0.001) and microsatellite data (r = 0.413, P = 0.001) both revealed significant correlations between genetic and geographic distances, but the correlations among chloroplast markers were weaker (Supplementary Appendix 3). Barriers inferred from nuclear genes showed a barrier between the YGP and HDM regions (Figure 4). As well as, distinct geographical isolation had been detected between the populations in southwestern Yunnan and NPR (Figure 4). Based on the TPM model, three populations were determined to have experienced bottleneck (PZH, XC, YJ; Table 1).

A total of three loci of nSSRs conformed to HWE, and five loci deviated from HWE (P < 0.05, Supplementary Appendix 1). In the Bayesian clustering analysis, the optimal K-value for STRUCTURE was K = 4 (Supplementary Appendix 4). However, the optimal K-value selected by InStruct was 2, and the second highest peak occurred for K = 4 (Supplementary Appendix 4). Thus, we performed cluster analysis on the Q. franchetii complex using K = 2 and K = 4, respectively for both STRUCTURE and InStruct. Following the InStruct runs with K = 2, cluster a was mainly composed of the YGP group, with some populations in the HDM (DL, BC, MD), and southwestern Yunnan (LC, SJ) groups. Cluster b was composed of 12 populations peripheral to those of cluster a. When K = 4, the YGP group was further divided into subgroups, and the populations from RR (PE, VM) separated from the rest of cluster b to form a new subgroup (Figure 4). The STRUCTURE results were identical to the results obtained from STRUCTURE (Supplementary Appendix 5). The principal coordinate analysis (PcoA) of the clustering results found that 33 natural populations can be divided into four regions, which is mostly consistent with the results obtained using InStruct. The first and second principal components explained 16.24 and 6.91% of the genetic variation, respectively (Supplementary Appendix 6).

The crown age of cpDNA haplotypes in the Q. franchetii complex was dated to the late Oligocene (30.7 Ma, 95% highest posterior density, HPD = 16.7–43.9 Ma), and Clade I (NPR) was the earliest derived. The haplotypes of the southwestern lineage (Clade II) began to diversify around 25.7 Ma (95% HPD = 11.9–37.2 Ma). The divergence of Clades III (HDM region) and IV (YGP region) was dated to ca. 24 Ma (95% HPD = 9.95–34.13 Ma), with subsequent rapid divergence of the haplotypes during the late Miocene (Figure 2).

Ecological niche modeling of the Q. franchetii complex in Maxent using WorldClim climate data revealed a high performance score (AUC = 0.9679–0.9719, standard deviation = 0.0157). Annual Mean Temperature (bio1) was the greatest contributor (36.46%, standard deviation = 1.3), followed by Temperature Seasonality (bio4) (34.7%, standard deviation = 1.3) and Annual Precipitation (bio12) (7.72%, standard deviation = 1.52) in identifying the areas of occurrence for Q. franchetii complex populations. The maximum sensitivity plus specificity value 0.13 was used as the species absence/presence threshold. The current distribution of the Q. franchetii complex was similar to the predicted distribution, except for quite a few occurrence sites in southwestern and southeastern Yunnan and Southern Himalayas that were located in predicted unsuitable areas (Figure 5B). The potential distribution range retreated to the south during the LGM (Figure 5A). The total distribution area of the LGM was greater than that of the present (N_a = 0.42, N_e = 21%). In the future period, the higher the RCP value was, the more obvious of an expansion in the suitable area to the northwest was detected. The ratios of the present range to the future range under RCP2.6, RCP4.5, RCP6.0, and RCP8.5 were 1.07, 1.15, 1.04, and 1.18, respectively. Accordingly, the highly suitable habitat would be reduced by 7.2, 7.2, 16.5, and 40.1%, respectively (Figures 5C–F). The predicted present distribution of Q. franchetii complex using CHELSA climate data was similar to that from WorldClim climate data (Supplementary Appendix 9). The suitable distribution area in the current period from the two data sources was largely overlapped (0.94). In order to better compare the distribution dynamics with other oak, of which WorldClim database were generally used for distribution simulation, we selected the results of WorldClim for the subsequent analyses.

Putative dispersal corridors in the two periods were visualized based on cpDNA haplotype diversity (Figure 6). The dispersal corridors during the two periods are consistent, and both showed that the YGP area had a higher dispersal ratio, but the populations of the peripheral areas were rather isolated.

The correlation analysis indicated that only the N_LGM was significantly correlated with GenPCoA1 (P < 0.05). A_r was associated with both longitude and latitude (P < 0.05). Cluster A (C_A) and GenPCoA1 were both associated with latitude (P < 0.05). However, H_e was not associated with N_Pre, N_stabLGM, N_LGM, latitude, nor longitude (P > 0.05) (Table 3).

Our study revealed a high genetic diversity estimate for the Q. franchetii complex (H_T, 0.982; A_r, 3.18–4.34), which is similar to the genetic diversity of other evergreen oaks in southwestern China, e.g., Q. schottkyana (H_T, 0.828; A_r, 4.83–7.78; Jiang et al., 2016), Q. kerrii (H_T, 0.71; A_r, 2.27–3.20; Jiang et al., 2018), Q. delavayi (H_T, 0.907; A_r, 3.750–5.237; Xu et al., 2020), and the Q. cocciferoides complex (H_T, 0.904; Liu, 2019). Comparatively, the genetic diversity for deciduous oaks seems lower, e.g., in the eight European white oaks (H_T, 0.635–0.847) (Petit et al., 2002), Quercus variabilis (H_T, 0.888, in 50 populations in East Asia) (Chen et al., 2012), Q. mongolica var. crispula (H_T, 0.827, in Japan) (Okaura et al., 2007), Q. acutissima (H_T, 0.791, in Southeast China) (Zhang et al., 2015). The deciduous species are mainly Miocene-Pliocene derived young lineages, but the divergence time of the evergreen oaks are much older that dated to the Eocene-Miocene (Hipp et al., 2020). It is easy to understand that oak taxa with long evolutionary histories and wide distribution can accumulate higher levels of genetic diversity than those “young” species. Additionally, the high environmental heterogeneity of southwestern China can buffer the climate extremes of the Quaternary episode (Zhang et al., 2010; Kai and Jiang, 2014; Jiang et al., 2018). As a result, the habitats of southwestern China have had long-term stability without significant regional extinction events or distribution range shifts to allow the genetic diversity of the species can be maintained. All these factors contributed to the high genetic diversity level found in evergreen oaks in southwestern China.

Comparing the genetic structure of the sympatric close related species can better illustrate the factors determining genetic diversity pattern. Notably, there was significant genetic differentiation among populations in the Q. franchetii complex, with F_ST estimates for cpDNA and nSSRs of 0.87 and 0.369, respectively, which is very similar to that found in Q. delavayi, as the cpDNA of the both species show significant phylogeographic structure, IBD pattern, and high differentiation among the populations (Xu et al., 2020). However, in two other sympatric/parapatric oaks, Q. kerrii (Jiang et al., 2018) and Q. schottkyana (Jiang et al., 2016), only low differentiation among the populations without phylogenetic structure were found (Table 4). The biological traits restrict gene flow and its efficiency, no doubt, can greatly impact the population genetic structure (Petit et al., 2003; Cavender-Bares et al., 2015; Xu et al., 2020). Pollen and seed mediated gene flow among the populations is different, as they have different dispersal efficiencies when barriers exist (e.g., rivers and high mountains). Generally, the seed-mediated gene flow among populations is more restricted in the species with instant germination seeds comparing to those species that seed with a period of dormancy. Consequently, the typical recalcitrant seed species shows significant phylogeographic structure in cpDNA makers. Vice versa, as pollen of oaks can disperse long distance, population differentiation revealed by biparental markers was much lower than that in maternal makers (Xu et al., 2020). However, the genetic structure of Q. franchetii seems not only determine by seed/pollen mediated gene flow efficiency. Although Q. franchetii has temporary seed dormancy (2–4 months; Xia et al., 2012), its seed size and tannin content similar to those of Q. schottkyana, but its population genetic structure is similar to that of Q. delavayi—a typical recalcitrant small seeds species. Thus, factors beyond seed germination schedule, seed size, and dispersal abilities also played important roles in shaping the population structure of these oaks in southwestern China, e.g., their ancestor distribution range and evolutionary history, etc.

Contemporary and historical factors shaped the genetic structure of organisms (Van Oppen et al., 2011; Hernawan et al., 2017; Li J. J. et al., 2017). Geological and climatic factors have been shown to influence the evolutionary histories of taxa and shape their genetic structures (Feng et al., 2014; Xing et al., 2014; Lu et al., 2018; Chen et al., 2020). Thus, the genetic structures of ancient lineages with long evolutionary histories and wide distribution range can essentially record more ancient geological events than in those of young lineages. The evergreen oak lineages in YGP and southwestern China, e.g., Q. schottkayana at 6.37 Ma (Jiang et al., 2016), Q. kerrii at 6–7 Ma (Jiang et al., 2018), and Q. cocciferoides at ca. 5 Ma (Jiang et al., 2019; Liu, 2019) are later derived “young” (originated at the late Neogene) and they generally show no (or very limited) IBD pattern among the populations. In contrasts, the early derived species Q. delavayi with crown node age at 10.92 Ma, and the Q. franchetii complex with crown node age at 30.7 Ma (95% HPD 16.7–43.9 Ma) had similar genetic structures and distinct IBD pattern. The different genetic structures detected in the “young” and “old” species might reflect the outcomes of past geological events on the biota at the different epochs, as the ancient geological events may imprint a genetic structure pattern in an “old” and widespread lineages, but not in those of the “young” lineages. However, our hypothesis requires further investigation, as the molecular markers we used in previous phylogeographical studies on oaks of southwestern China are not all the same. Regardless, these markers are informative to reveal the spatial genetic structures of these oaks, but it contains potential bias when comparing across the species to reveal the underlying mechanisms that shaped the genetic structures of these oaks. Further studies using universal high-throughput markers to scan oak populations in southwestern China and incorporating seed functional traits analyses are necessary to reveal the underlying drivers lead to the contemporary genetic structure found in these oaks.

The Oligocene was a key period in the formation of the modern topography of China. During this period, the high eastern and low western topographies of China were totally reversed (Ming, 2007; Di et al., 2018). The tectonic plates associated with the spreading of the East Asian marginal sea and the uplift of the eastern margin of the Tibetan Plateau exerted important influence on the topography and drainage of East Asia (Clark et al., 2004). This was accompanied by dramatic tectonically induced topographic changes and landscape development, which resulted in repeated river captures and reversals (Zheng, 2015). Paleomagnetic and stratigraphic evidence suggests that there was a wide Paleo-Red River drainage basin between the southeastern Tibetan Plateau and the South China Sea, which included Hoh-Xil, Songpan-Ganzi, northern Qiangtang, Yidun, and western Yangtze Terranes (Figure 7). Since the late Eocene, the hard collision between the Indian and Eurasian plates began. These massifs/blocks were extruded eastward under the resulting pressure (Chen et al., 2017). Meanwhile, the Ailao Shan-Red River fault zone began to slip rapidly. The Lanpin-Simao block began to rotate clockwise substantially, which blocked the upper and lower Yangtze River from continually flowing to the south along the paleo-drainage (Figure 7). As a result, the northern drainage basin of the paleo Red River disappeared, and the modern Red River began to become established (Replumaz and Tapponnier, 2003; Royden et al., 2008; Chen et al., 2017; He et al., 2021). These the Oligocene to early Miocene events greatly changed the regional biota, e.g., giving rise to the diversification of Cautleya, Roscoea (Zhao et al., 2016), Badidae (Rüber et al., 2004), and spiny frogs (Paini) (Che et al., 2010).

Our ancestral range reconstruction based on cpDNA data showed that the Q. franchetii complex once had a wide distribution in southwestern China and the southern Himalaya regions during the early Oligocene, followed by three vicariance events. Among these events, the NPR (Clade I) was first diversified during the mid-Oligocene (ca. 30 Ma), and then, during the late Oligocene and the early Miocene, the RR lineage and HDM lineage were derived, respectively. Within the main lineages, the fast divergence of the cpDNA haplotypes occurred during the late Miocene (Figure 7A). The NPR was located at the core area during early river capture; the hard collision between the Indian and Eurasian plates might have squeezed the plates in this region leading to river re-alignment, which induced the divergence of the NPR lineage (Clark et al., 2004; Yang et al., 2012; Figure 7B). During the late Oligocene, the Ailao Shan-Red River fault had an early left strike-slip that may have raised the barrier that blocked gene flow among the populations and promoted allopatric divergence (Hall et al., 2008; Fyhn et al., 2009; Lin et al., 2009; Deng et al., 2014). Thus, the RR lineage diverged (Clade II). During the early Miocene, large fault basins were established corresponding to the fast uplift of the YGP and HDM regions (Zhang, 2012), which might dramatically change the regional topography and climate, eventually blocking gene flow in the two regions and leading to the divergence of the YGP lineage (Clade III) and HDM lineage (Clade IV). Followed by the late Neogene period fast HDM uplifts, the complex topography of southwestern China was eventually formed, which further restricted regional seed-mediated gene flow and promoted the divergence of cpDNA haplotypes. A similar scenario was also detected in other oaks, e.g., Q. delavayi (Xu et al., 2020), and Q. aquifolioides (Du et al., 2017), as well as plant lineages with wide distribution in semi-humid evergreen broadleaved forests on YGP region, e.g., Primula secundiflora (Wang et al., 2008), Terminalia franchetii (Zhang et al., 2011), and Cycas multipinnata (Gong et al., 2015). Such phenomenon suggested the regional biota might be impacted by similar environmental drivers.

In summary, the high biodiversity levels found in southwestern China are rooted deeply in the Oligocene. The early tectonic events during the Oligocene drove the main lineage splits, while the fast uplifts of the Himalayas during the Miocene-Pliocene increased environmental heterogeneity and established substantial dispersal barriers (Meng et al., 2017; Xing and Ree, 2017; Yuan et al., 2019). Then, the Quaternary climate fluctuation led to distribution range contractions and expansions of the species, as well as the occurrence of Asian winter monsoons, and the dry season in winter and spring in southwestern China further restricted gene flow between the core YGP region and its periphery (Su et al., 2013; Ye et al., 2019). All these geological and climatic factors interacted during different timespans to shape the contemporary divergence pattern of the biota of the southeastern Himalaya biodiversity hotspot. Our phylogeographic study indicated the Oligocene tectonic induced divergence in Q. franchetii complex, which is a supplement to the review of Renner (2016) concerning the Paleogene events contributing to the species richness of the East Himalayan biodiversity hotspot.

Moreover, the population genetic structures inferred from nSSR and cpDNA markers of Q. franchetii complex were dissimilar. Notably, nSSRs mainly reflected the divergence between the populations in the core YGP region and the peripheral populations. The similar population genetic structure on nuclear genome was also reported in another sympatric species (Q. cocciferoides; Liu, 2019). All these evidences suggested the two species might underwent the similar selection pressure to trigger their divergence. The potential migration corridor analysis on the Q. franchetii complex showed that the populations in core YGP region maintained strong gene flows, but the marginal populations were mostly isolated since LGM (Figure 6). The rugged topography induced by the rapid YPG uplift during the late Neogene and (or) highly fragmented habitat of semi-humid evergreen forests in the peripheral areas around YGP might boost the allopatric divergence of these oaks. Nevertheless, another possibility of this pattern is that the quick evolution and possible backward evolution of SSR markers blurs the geological pattern. Further investigation using high throughput marker to illustrate the genetic structure at fine scale can provide a better understanding on the interplays between genetic diversity and environmental factors.

In contrast, the cpDNA data of Q. franchetii complex showed a much clearer phylogeographic structure, which has also been shown in other oaks, e.g., Q. delavayi (Xu et al., 2020), Q. cocciferoides (Liu, 2019), and Q. aquifoliodes (Du et al., 2017). Pollen- and seed-mediated gene flow have very different dispersal efficiencies in oaks (Du et al., 2017; Liu, 2019; Xu et al., 2020). In this study, nSSRs showed no significant differentiation of the populations in the East Red River, West Red River, and HDM regions, but cpDNA data showed a clear phylogeographic structure. This result suggested that these early tectonic activities during the Oligocene to early Miocene might have restricted seed-mediated gene flow in different regions, but only had limited impacts on pollen-induced gene flow.

In contrast, the niche modeling result suggested that the Q. franchetii complex populations are mainly located in the predicted suitable area. While the species complex experienced a southern contraction during the LGM, the distribution area at present and during the LGM largely overlapped in YGP. The quaternary glaciation had only minor impacts on its distribution. Likewise, the niche modeling results of Osteomeles schwerinae (Wang et al., 2015), Q. schottkyana (Jiang et al., 2016), and Q. kerrii (Jiang et al., 2018) from southwestern China show the similar pattern. Collectively, these studies suggest that central YGP region is an important refugia for species in semi-humid evergreen broadleaved forests in southwestern China.