Between May 2014 and September 2017, we collected fresh leaf material of 419 individuals from 18 populations throughout the entire distribution of Q. chenii in China ( Figure 1A and Supplementary Table S1 ). In each population, eight to ten fresh leaves per tree were sampled for nine to 30 adult individuals at least 30 m apart from each other. Two species from Quercus section Quercus (Q. fabri and Q. aliena) and two species from genus Castanea (C. mollissima and C. henryi) were used as outgroups. Leaf tissues were quickly dried with silica gel and stored at room temperature in the laboratory. Spatially explicit information was recorded for each population using a handheld GPS unit (Magellan, USA). Voucher specimens of all individuals of Q. chenii and outgroups are deposited in the Herbarium of Nanjing Forestry University (HNFU) ( Supplementary Table S1 ).

Total genomic DNA was extracted from 30 mg leaf tissue of each individual using a Plant Genomic DNA Kit (Tiangen, Beijing, China). Four chloroplast intergenic spacers, atpB-rbcL, psbA-trnH, trnS(GCU)-trnG(UCC), and trnS(GCU)-trnT(GGU), were sequenced for six to ten individuals per population following protocols in Zhang et al. (2015). For primer details, please refer to Supplementary Table S2 .

All the 419 samples were genotyped at 14 nSSR loci, including quru-GA-0M07, quru-GA-1H14, quru-GA-1M17, and quru-GA-2G07 (Aldrich et al., 2002), ssrQrZAG4, ssrQrZAG7, ssrQrZAG11, ssrQrZAG59, ssrQrZAG96, and ssrQrZAG112 (Kampfer et al., 1998), ssrQpZAG15 and ssrQpZAG110 (Steinkellner et al., 1997), QM57-3M and QM67-3M1 (Isagi and Suhandono, 1997) ( Supplementary Table S3 ). The forward primers of each locus were labeled with fluorescent dyes 6-FAM or HEX (Applied Biosystems, USA). PCR reactions were performed following Zhang et al. (2018), with primer-specific temperature described in Supplementary Table S3 . PCR products were separated on an ABI3730xl automated Genetic Analyzer using ROX-500 as an internal standard (Applied Biosystems, USA). Allele sizes were determined manually using GENEMARKER 2.2.0 (Applied Biosystems, USA).

We employed a maximum entropy approach in MAXENT 3.4.1 (Phillips et al., 2018) to simulate the modern distribution of Q. chenii. Species occurrence records were mainly collected from the fieldwork, the literature, and the database of the Global Biodiversity Information Facility (GBIF, https://www.gbif.org/), the Chinese Virtual Herbarium (CVH, http://www.cvh.ac.cn/) and the Plant Photo Bank of China (PPBC, http://www.plantphoto.cn/). We filtered the data by removing duplicate records and retaining only one record among all locations falling within the same 2.5′ × 2.5′ grid. Finally, a total of 54 occurrence points were obtained for Q. chenii. The environmental layers of the eight bioclimatic variables ( Supplementary Table S6 ) were downloaded from the Worldclim database (http://www.worldclim.org/) with a resolution of 30″ for the present, the Mid Holocene (∼6 ka), the Last Glacial Maximum (LGM, ∼22 ka), and the Last Interglacial (LIG, ∼120–140 ka) under the Community Climate System Model 4 (CCSM4). All layers were clipped to the same spatial range (15°–40° N, 105°–140° E). The optimal settings of feature types (linear and quadratic) and regularization multiplier (value = 1) were selected through the R package ‘ENMeval’ (Muscarella et al., 2014). Model performance was evaluated using the areas under the receiver operating characteristic curve (AUC) produced by 10-fold cross-validation. The generated models were then projected onto the other three scenarios of past climates. Obtained layers were used to calculate the average, minimum, and 1-SD of occurrence probabilities across different historical periods in ARCMAP 10.1 (ESRI). The resulting layers were reclassified to show putative areas with moderate and high ecological stability (Chan et al., 2011), using 80% and 90% (for average), and 50% and 80% (for minimum and 1-SD) of the maximum values as thresholds. To assess the landscape connectivity of the cpDNA haplotype network, we applied SDMTOOLBOX (Brown, 2014) in ARCMAP 10.1 (ESRI) to convert species distribution models to ‘dispersal cost’ layers, and then to create a raster of the sum of least-cost corridors between populations that share haplotypes.

Using INEST, the frequency of null alleles was estimated to be lower than the threshold of 0.05 at each of the 14 loci across populations ( Table 1 ). Only three of the 91 locus pairs showed significant LD in two populations (quru-GA-1M17 × ssrQrZAG4 in QI, ssrQrZAG59 × ssrQrZAG4 in QI, and quru-GA-1M17 × ssrQpZAG15 in TH; P < 0.05 after Bonferroni correction). No consistent genotypic disequilibrium was found between any locus pairs across all populations, so all loci were used for further analyses. Significant deviation from HWE was only detected in one of the 18 populations (P < 0.05 after Bonferroni correction; Table 1 ).

At the population level, allelic richness (A _R) ranged from 3.849 to 5.901, and genetic diversity within populations (H _S) varied from 0.585 to 0.713 ( Table 1 ). Highland populations presented lower genetic diversity than lowland populations (P < 0.001 for both A _R and H _S; Table 2 ). Genetic differentiation among populations was significant over all loci (F _ST = 0.054, P < 0.001; G′_ST = 0.228, P < 0.001; Supplementary Table S3 ). Pairwise F _ST ranged from 0.015 to 0.124, with 152 of all the 153 pairs being significant (P < 0.05 after Bonferroni correction). The highest pairwise F _ST values were observed for JZ-ZN, corresponding to two montane habitats located at the northwest (the Dabie Mountains) and southeast (the Tianmu Mountains). A higher level of genetic differentiation was detected among highland populations (F _ST = 0.099, range of pairwise F _ST: 0.072–0.124) than among lowland populations (F _ST = 0.036; range of pairwise F _ST: 0.015–0.059) ( Table 2 and Figure 3 ). Genetic differentiation between these two groups was low but significant (F _CT = 0.004, P = 0.013), with only 0.43% of the variation partitioned among groups, and the most (94.33%) partitioned within populations ( Table 3 ).

Bayesian cluster analysis showed that the highest ΔK occurred at K = 2 and 3 ( Supplementary Figure S1 ). When K = 3, a significant decline in the probability of membership to genetic cluster III (Q _III) in each population was detected with increasing elevation (R = 0.711, P = 0.001) ( Supplementary Figure S2 ). According to this trend, we divided all populations into two groups: (1) Highland populations. Cluster III presented a much lower Q value ranging from 0.01 to 0.06. Four of those were dominated by genetic cluster I (Q _I: 0.64–0.94), and two were dominated by cluster II (Q _II: 0.93–0.94) ( Figure 1A , Supplementary Table S8 and Figure S3 ). (2) Lowland populations. Cluster III presented a higher Q value varying from 0.14 to 0.89. The other two clusters showed a Q value lower than or equal to 0.60 ( Figure 1A , Supplementary Table S8 and Figure S3 ). A higher level of genetic admixture was detected in lowland populations (mean of A _D = 0.578) than in highland populations (mean of A _D = 0.134) ( Table 2 ). The BARRIER analysis revealed multiple genetic boundaries that were well associated with persistent landscape barriers of East China. Among those, highland populations were isolated from lowland ones by the Dabie Mountains, the Jiuhua Mountains, the Baiji Mountains, and the Tianmu Mountains, with more than 70% bootstrap support. Lowland populations inside the Huaiyu Mountains, the Mufu Mountains, the Jiuling Mountains, and the Lu Mountains were also separated from the adjacent ones outside with more than 70% bootstrap support ( Figure 1B ).

ABC analyses clearly favoured the hypothesis that Q. chenii may have experienced a recent demographic expansion (scenario 4) within both high and low elevation regions ( Supplementary Table S9 ). All the summary statistics together with PCA plots indicated a good fit of this scenario to the observed data ( Supplementary Table S4 and Figure S4 ). Assuming an average generation time of 100 years, the recent expansion event was dated to 0.121 Ma (95% CI: 0.008–0.290 Ma) for the highland populations, with effective population size varying from 17,600 (95% CI: 868–243,000) to 526,000 (95% CI: 41,800–977,000). For the lowland populations, the recent expansion event was dated to 0.138 Ma (95% CI: 0.010–0.291 Ma), with effective population size varying from 12,800 (95% CI: 604–196,000) to 533,000 (95% CI: 46,000–977,000) ( Table 4 and Supplementary Figure S5 ). Both times corresponded to the Last Interglacial (∼0.14–0.12 Ma). Significant evidence of recent bottleneck was only detected in populations HS and TH under the two-phase mutation model (TPM) (P < 0.05; Supplementary Table S8 ).

Genetic diversity statistics (A _R and H _S) and genetic admixture index (A _D) were not associated with longitude or latitude (all P-values > 0.05), but correlated with elevation and climatic variables. The model that best explained genetic diversity contained elevation (all P-values < 0.01), with mean diurnal temperature range (bio2), temperature seasonality (bio4), and mean temperature of driest quarter (bio9) as covariates. The model that best explained A _D included only elevation (P < 0.01) ( Table 5 ). Significant negative correlations between elevation and these three genetic statistics were also revealed by simple linear regression (all P-values < 0.01, Figure 4 ). When using pure climatic predictors to describe these relationships, annual mean temperature (bio1) and temperature seasonality (bio4) remained in all the best-fit models. Significant positive associations of bio1 and bio4 were detected with A _R, H _S, and A _D (all P-values < 0.05). Other covariates, including mean diurnal temperature range (bio2) and mean temperature of wettest quarter (bio8), showed marginally significant (0.05 < P < 0.10) or insignificant (P > 0.10) positive correlations with genetic diversity ( Table 5 ). Multicollinearity was not detected in all the models as indicated by VIF < 2 for each variable ( Table 5 ).

The full RDA model (df = 16, F = 3.726, R _adj ² = 0.094, P = 0.001), and both partial RDA models corresponding to geography (df = 8, F = 3.591, R _adj ² = 0.046, P = 0.001) and climate (df = 8, F = 3.518, R _adj ² = 0.044, P = 0.001) were significant. The full model explained 12.91% of the total variance. The percentages of genetic variance attributed to geography and climate alone were 6.22% and 6.10%, respectively. Only 0.59% of the overall genetic variance was explained by the collinearity between geographic and climatic variables ( Supplementary Table S10 ). In the full model, all the eight spatial descriptors and the eight climatic predictors were significant (all P-values = 0.001; Supplementary Table S7 ). Forward stepwise selection also kept all the 16 variables in the best model based on the adjusted coefficient of determination (R _adj ²). Among those, three showed the highest scores on the first three axes, i.e., elevation on the RDA1, temperature seasonality (bio4) on the RDA2, and PCNM9 on the RDA3 ( Figure 5 and Supplementary Table S7 ).

The lengths of consensus sequences after alignment of atpB-rbcL, psbA-trnH, trnS(GCU)-trnG(UCC), trnS(GCU)-trnT(GGU), and concatenated cpDNA were 718, 639, 607, 902, and 2,866 bp, respectively ( Supplementary Table S2 ). Supplementary Table S11 Eighteen haplotypes identified in this study based on 16 nucleotide substitutions and six indels. Among those, only one (H1) was shared by seven populations, two (H6 and H18) were shared by two populations, and the other 15 were private to a single population. More than 70% populations were fixed for a single haplotype ( Figure 6 ). Such a geographic pattern of haplotypes resulted in a much lower average within-population genetic diversity (h _S = 0.126) compared with the total genetic diversity (h _T = 0.917) at the species level. Intraspecific differentiation at cpDNA markers was significant (F _ST = 0.887, P < 0.001), with 88.65% of the total genetic variation partitioned among populations, and 11.35% partitioned within populations ( Table 3 ). However, we did not detect any phylogeographic structure across the entire distribution (N _ST = 0.888 > G _ST = 0.863, P = 0.453). The genetic divergence between highlands and lowlands was also not significant (F _CT = 0.006, P > 0.05). Although the mismatch distribution for all haplotypes was bimodal (Supplementary Figure S6 ), both sum of squared deviation (0.017, P = 0.222) and Harpending’s Raggedness index (0.043, P = 0.140) did not reject the sudden expansion model. Following Qiu et al. (2009b), this statistical fit of the expansion model is here not taken as strong evidence of expansion. Tajima’s D test (Tajima’s D = −0.210, P = 0.475) and Fu’s F _S test (Fu’s F _S = −0.257, P = 0.548) also did not show evidence of extensive demographic expansion.

Both median-joining network and phylogenetic analysis grouped all the 18 haplotypes into three clades ( Figures 6B and 7 ). Clade A included eight haplotypes that were geographically scattered in the Mufu-Jiuling-Lu mountain region, and the Jiuhua-Huangshan-Tianmu-Baiji-Huaiyu mountain region. Clade B was characterized by a notable star-like pattern, with ancestral haplotype H1 and seven derived haplotypes widely distributed in both northern and southern regions. Clade C contained only two haplotypes that were separated from clade B by five to seven mutations. Among those, H16 was only detected in population NJ, and H14 was confined to population ZN in the East Tianmu Mountains ( Figure 6B ). The time to the most recent common ancestor of all the haplotypes was dated to the early Miocene (16.70 Ma, 95% HPD: 10.10–23.99 Ma). Clade A and clade B were estimated to have diverged at the end of the middle Miocene (11.63 Ma, 95% HPD: 6.78–17.73 Ma). The crown ages of clade A, B, and C were about 8.12 Ma (95% HPD: 4.31–13.13 Ma), 8.12 Ma (95% HPD: 4.06–12.99 Ma), and 6.21 Ma (95% HPD: 1.96–12.21 Ma), respectively ( Figure 7 ).

The mean AUC value (± SD) on the test dataset was 0.921 ± 0.031, indicating a good fit of ENM to the observed occurrence data of Q. chenii. Historical distributions for the LIG, LGM, and Mid Holocene are shown in Supplementary Figure S7 . An area with high ecological stability in the mountainous area of East China was identified by calculating the average, minimum, and 1-SD of occurrence probabilities across different periods ( Figures 8A–C ). A major north–south dispersal corridor along the plains and hills on the eastern side of the Poyang Lake Basin was detected for the LGM, the Mid Holocene, and the present. Two east–west dispersal corridors along the hills in the central Jiangxi Province and north of the Mufu Mountains were also identified for the LGM and the present, respectively ( Figures 8D–F ).

Our current study of Q. chenii indicated that a higher level of genetic differentiation occurred among highland populations than among lowland populations (P < 0.001, Table 2 ). A more detailed analysis using BARRIER identified multiple genetic boundaries that were well associated with the scattered distribution of mountain ranges in East China ( Figure 1B ). Such findings pinpointed that mountains as landscape barriers may have driven the population differentiation through long-term geographic isolation. The occurrence of multiple geographic barriers was also supported by the distribution pattern of cpDNA haplotypes ( Figure 6 ). Specifically, no shared haplotypes were detected between populations inside and outside the region surrounded by the Mufu Mountains, the Jiuling Mountains, and the Lu Mountains, and the region surrounded by the Jiuhua Mountains, the Huangshan Mountains, the Tianmu Mountains, the Baiji Mountains, and the Huaiyu Mountains. More strikingly, five of the seven haplotypes within these two mountainous regions were private to a single population. This pattern was more likely to result from long-term in situ diversification among isolated montane habitats. Considering the limited seed-mediated gene flow for oaks, differentiation may be enhanced by stochastic factors such as genetic drift (Zhang et al., 2013). Similar patterns with a geographic mosaic of cpDNA haplotypes have also been observed for several wind-pollinated tree species amid the complex topography in subtropical China, e.g., Fagus lucida (Zhang et al., 2013), Fagus longipetiolata (Zhang et al., 2013), and Juglans cathayensis (Bai et al., 2014), suggesting that the persistent landscape barriers have played an important role in shaping the pattern of long-term isolation for the extant plant species in this region.

In addition to geographic isolation, mountains may steepen climatic gradients, and further facilitate population differentiation through environmental isolation (Ohsawa and Ide, 2008; Sexton et al., 2016). For wind-pollinated trees like oaks, the pollen-mediated gene flow is not likely to be constrained by discontinuous habitats sharing similar environmental conditions (Wanderley et al., 2018), but can be disrupted by environmental dissimilarity (Sexton et al., 2014). Thus, a significant pattern of IBE would be expected for wind-pollinated plants, which has been observed for ash (Temunović et al., 2012), larch (Nishimura and Setoguchi, 2011), and also several oak species in East Asia (Q. liaotungensis, Yang et al., 2018) and North America (Q. lobata, Gugger et al., 2013; Q. engelmannii, Q. berberidifolia, Q. cornelius-mulleri, Riordan et al., 2016). Using redundancy analysis for nuclear microsatellite data, we revealed that both IBE and IBD explained similar amounts of constrained genetic variation across the species’ distribution, and elevation contributed the most to explaining the genetic variation constrained by RDA axis 1 ( Figure 5 and Supplementary Table S7 ), highlighting that environmental isolation associated with the complex topography in East China was also responsible for the local patterns of neutral microsatellite variation of Q. chenii.

The genetic differentiation along altitudinal gradients was also supported by Bayesian cluster analysis, which assigned a much lower proportion of genetic cluster III to the high elevation populations ( Figure 1A ). Multiple clues can be integrated to explain these findings. First, during the fieldwork, an elevation-related difference was observed regarding vegetation cover and topographic heterogeneity. Specifically, most individuals of lowland populations were found to be accompanied by short shrubs or sparse deciduous trees on plains and hills ( Figure 1C ), whereas most individuals of highland populations are scattered in thick forests, which are dominated by evergreen broad-leaved trees like Castanopsis eyrei and Quercus glauca, and have a continuous distribution across mountains and valleys ( Figure 1D ). Such a difference may have led to asymmetric pollen gene flow, and further reduced the genetic connectivity between highland and lowland populations (Gharehaghaji et al., 2017). Second, microclimate is also associated with elevation. Highland habitats are characterized by significantly lower temperature of spring, and a wetter environment surrounded by montane landscapes (Wang et al., 2010). These factors can influence the asynchronism of reproductive phenology, e.g., the timing of flowering and pollen release (Montoya-Pfeiffer et al., 2016), and then increase the effects of IBE through nonrandom mating due to climatic differences along altitudinal gradients. Finally, highland populations exhibited stronger differentiation and much lower within-population diversity than those in lowland areas ( Table 2 and Figure 3 ), suggesting that the species may have experienced genetic drift and long-term fragmentation among multiple island-like highland habitats, probably a consequence of historical colonization along altitudinal gradients and following ecological isolation (Qiu et al., 2009a; Wang et al., 2017).

Complementary to elevation, other environmental variables, such as temperature seasonality, annual mean temperature, and precipitation of warmest quarter, also contributed to the geographic patterns of climatically structured genetic variation significantly ( Supplementary Table S7 ). Previous studies have shown that both temperature and precipitation factors can affect the timing of reproductive phenophases of plants, including oaks (Matthews and Mazer, 2016; Gerst et al., 2017). Thus, nonrandom mating due to asynchronism of reproductive phenology may be a critical driver that shapes the pattern of environmental isolation for Q. chenii. More investigations on gene movement, natural selection, and local adaption to different ecological processes are expected to provide a better understanding of the pattern of genetic variation for Q. chenii.

Contrary to the role of mountains as geographic barriers, lowlands are more likely to be dispersal corridors for Q. chenii. Using least-cost path calculations, we confirmed that a major south–north dispersal corridor occurred along the plains and hills on the eastern side of the Poyang Lake Basin, and two east–west dispersal corridors occurred along the hills in the central Jiangxi Province and north of the Mufu Mountains ( Figures 8D–F ). Our investigations supported the conclusions of Fan et al. (2017), which showed that stronger population connectivity occurred in the eastern Poyang Lake Basin corridor for five plant species in subtropical China, and connections also occurred between the eastern and the western populations across the basin. However, the corridors detected in this study occupied a much lower elevation than those identified by Fan et al. (2017) (e.g., the Wuyi Mountains and the Yushan Mountains), implying that lowlands are important for the migration of plants that prefer a habitat at low elevation regions.

Furthermore, although the mismatch distribution analysis, Tajima’s D and Fu’s F _S tests did not show strong evidence of extensive demographic expansion, the star-like cpDNA phylogeny of clade B and the widespread distribution of the ancestral haplotype H1 suggested that lowland populations could have experienced local expansions during the interglacial periods. These results indicated that lowlands may have contributed to the post-glacial south-to-north range shifts of Q. chenii. The species may have had a southern refugium, probably at the locations of populations LC and LA that harbored the ancestral haplotype H1. During the interglacial periods, Q. chenii may have migrated northward to the region north of the Yangtze River along the dispersal corridors identified in this study.

The role of lowlands as dispersal corridors was also supported by the evidence inferred from nuclear microsatellite variation. The general linear model showed that elevation best explained the spatial pattern of genetic admixture across the species’ range ( Table 5 ). Compared with highland populations, lowland populations exhibited a higher degree of genetic admixture, but a lower level of population differentiation ( Table 2 ). These results suggested the landscape features in East China play different roles in shaping the patterns of gene flow within the species. Mountains as geographic barriers may have reduced gene exchanges among isolated habitats, and led to stronger patterns of both IBD and IBE. By contrast, lowlands were more likely to be dispersal corridors, which were critical for the post-glacial colonization, and contributed to the increased genetic admixture due to the arrival of immigrants originating from genetically differentiated populations (Ortego et al., 2015). To further support this hypothesis, we performed the ABC analysis to reconstruct the past demographic history of Q. chenii. The results clearly favoured the scenario that the lowland populations of Q. chenii may have experienced a recent expansion during the LIG, which is in agreement with the much larger potential distribution predicted by ENMs for this period ( Supplementary Figure S7A ). Overall, both cpDNA and nSSR markers showed evidence that lowlands were more likely to be dispersal corridors for Q. chenii, especially during the periods of range expansions. However, due to different effective population sizes and dispersal capabilities for nuclear and chloroplast genomes, they responded to topographic and climatic factors in different ways. Seed-mediated gene flow of oaks is more likely to be restricted by habitat fragmentation, and cannot counteract the effects of genetic drift (Petit et al., 1997). Thus, we detected shared haplotypes among lowland populations, but without much variation within populations. By contrast, pollen gene flow is much stronger than genetic drift, together with the higher migration rate, it would increase the genetic connectivity among regions, and also result in a higher level of both admixture and within-population diversity for lowland populations.

Accumulating evidence has shown that the tectonic–climatic interactions since the Miocene have significantly influenced the evolutionary dynamics of plant species in East Asia (e.g., Kou et al., 2015; Yu et al., 2017; Deng et al., 2018). Similar to prior studies, we conclude that the lineage divergence of Q. chenii has been shaped by the palaeoclimatic changes of the Neogene. The dating analyses of major cpDNA clades suggested that Q. chenii may have experienced the earliest divergence during the early Miocene (node E, 16.70 Ma, 95% HPD: 10.10–23.99 Ma, Figure 7 ). This age is more recent than the earliest unambiguous fossil records of Quercus section Cerris from the early Oligocene of Russia Far East (Denk et al., 2017), but close to the first occurrence of the reliable fruit fossils that were morphologically closely related to extant oak species from the section (i.e., Q. acutissima and Q. variabilis) in China (Song et al., 2000), implying a probable early Miocene origin of Q. chenii and its siblings. The estimated time here is also comparable with the crown ages of other tree species endemic to subtropical China, such as Cyclocarya paliurus (16.69 Ma, 95% HPD: 8.42–27.86 Ma, Kou et al., 2015), and dominant woody genera of subtropical evergreen broadleaved forests in East Asia, such as Machilus (15.8 Ma, 95% HPD: 9.9–21.6 Ma, Li et al., 2011), Castanopsis (17.7 Ma, 95% HPD: 11.5–24.7 Ma, Xing et al., 2014), and Michelia (19.2 Ma, 95% HPD: 13.2–25.9 Ma, Nie et al., 2008), suggesting that the initiation, or more likely, the intensification of the East Asian monsoon system (Sun and Wang, 2005; Licht et al., 2014), and subsequent increased humidification in the southern China (Guo et al., 2008) around the Oligocene–Miocene boundary may have triggered the lineage divergence of Q. chenii and other plant species of subtropical evergreen broadleaved forests in East Asia (Yu et al., 2017; Deng et al., 2018).

The time-calibrated phylogeny indicated that the divergence between clade C and D was likely to have occurred at the boundary of the middle and the late Miocene (11.63 Ma, Figure 7 ). Three major subclades A–C were estimated to have diverged during the late Miocene (8.12–6.21 Ma, Figure 7 ). Moreover, most haplotypes were found to have diverged during the late Pliocene to the early Pleistocene (3.86–1.60 Ma, Figure 7 ). Such findings were largely in agreement with the global climatic cooling after the Middle Miocene Climatic Optimum (17–15 Ma; Sun and Zhang, 2008), and the enhancement of East Asian monsoon intensity around 8–7 Ma and 3.5–1.6 Ma (An et al., 2001; Sun and Wang, 2005). Similar divergence history has also been reported for several woody taxa in East Asia, such as Cyclocarya paliurus (Kou et al., 2015), Quercus section Cyclobalanopsis (Deng et al., 2018), and Theaceae (Yu et al., 2017). More interestingly, a recent research also provided quantitative evidence to support the positive association of diversification rates with both global cooling and East Asian monsoon activity for Primulina, one of the most diverse herbaceous genera amid the complex karst landscapes in southern China (Kong et al., 2017). The combination of these findings pinpoints that the Neogene climate dynamics, especially the phased intensification of East Asian monsoons and global cooling of the Miocene, have stimulated the lineage divergence of plants in this region no matter whether they are herbs or long-lived trees.

Our analyses suggested that an area surrounded by multiple geographic barriers, i.e., the Tianmu Mountains, the Huangshan Mountains, the Baiji Mountains, and the Huaiyu Mountains, may have acted as a refugium for Q. chenii. This hypothesis was supported by three lines of evidence. First, ENM indicated that this area was characterized by moderate to high ecological stability throughout the LIG to the present ( Figures 8A–C ), implying that Q. chenii probably have survived in situ during the glacial periods of the Pleistocene. Second, a higher level of phylogenetic diversity was detected here. Although most populations were fixed for a single haplotype, all the three cpDNA haplotype clades, including several endemic haplotypes were concentrated, and populations dominated by the three genetic clusters identified by Bayesian cluster analysis were also confined to this area. Finally, macrofossils suggested that Q. chenii went extinct in central Japan along with the onset of glaciation in the early Pleistocene (Momohara, 2016), and probably retreated southward to the refugia on the East China Sea shelf (Li et al., 2016; Supplementary Figure S7B ), and then to the refugia in Mainland China. Palaeovegetation reconstruction of the LGM also supported that temperate deciduous trees may have confined to the low elevation regions of East China and extended to the offshore refugia on the continental shelf (Harrison et al., 2001). Thus, the combined evidence of long-term climatic stability, higher phylogenetic diversity, and a location southwest of Japan and close to East China sea land bridge suggested that this region may have served as a mainland refugium for Q. chenii. Recent studies also revealed that plants with a larger distribution may have survived in a refugium located at eastern subtropical China (e.g., Wang et al., 2017; Hohmann et al., 2018). Our inference strengthened this belief and confirmed that the tree species mainly distributed in East China may have persisted in situ, and more strikingly, have been isolated for a long term among multiple habitats within the mountainous area. However, caution is still needed, given that several haplotypes detected in this study were shared with a closely related species, Q. acutissima (Li et al., in press), although the same region may also have been a refugium for the latter, implying a probable scenario of interspecific introgression due to secondary contact or long-term co-occurrence in a shared refugium.