Three sympatric oak species: Q. aliena, Q. dentata, and Q. variabilis were studied in their northern distribution, at Shangfang Mountain located in the southwest of Beijing Municipality, P. R. China (115°48′E, 39°39′N). Shangfang Mountain is a branch of the Taihang Mountains with fold-thrust geological structures. It has a mean altitude of 400 m and a summit of over 880 m. The mountain slopes range from 20 to 70°, forming very steep terrain with dozens of valleys. Within this region, hundreds of ancient Taoist and Buddhist temples scatter throughout the forest, with the oldest one being dated back to the Eastern Wei period (535 AD). The forest here is a well-preserved secondary forest in north China covering an area of approximately 340 hectares. Tree species of Shangfang Mountain mainly belong to Quercus, Pinus, and Platycladus with Quercus including mostly Q. aliena, Q. dentata, and Q. variabilis and a limited number of <20 individuals of Q. mongolica.

Leaf material was collected from seven 50 by 50 m sites in this region (Figure 1 and Supplementary Table 1). The sites were selected for each species on the basis of site accessibility and sample availability. In general, two leaves were collected from each mature tree of each of the three oak species provided that they were at least five meters apart and stored in silica gel for further analysis. In total, 432 individuals were collected including 207 Q. aliena (six sites), 68 Q. dentata (three sites), 156 Q. variabilis (five sites) and only one individual with intermediate morphology between Q. aliena and Q. dentata (at site TK). The latitude, longitude, and altitude of each individual sampled were recorded by GPS (Garmin, USA).

Genomic DNA was extracted from silica gel dried leaves using a plant genomic DNA kit (Tiangen Biotech, China) following the manufacturer's instructions. The quality of DNA was initially assessed by a 1.00% (w/v) agarose gel. Then, the quantity of DNA was measured by an ultramicro-spectrophotometer (Thermo Fisher, USA). Genomic DNA of each sample was adjusted to a final concentration of 20–30 ng/μL for subsequent use.

For the chloroplast DNA (cpDNA) study, 14 individuals from the seven sites (at least one individual of each oak species present at each site) were initially used to identify polymorphisms with 13 pairs of universal primers described in Shaw et al. (2005, 2007) and Du et al. (2017) (see Supplementary Table 2 for details). Five pairs of cpDNA primers (accD-psaI-75R, trnH-psbA, rpS12-rpL20, rpS16, and trnQ-trnS) that exhibited polymorphisms in this first scan were subsequently used to analyze all 432 individuals, out of which 343 individuals with high quality were taken forward for further investigation. Each polymerase chain reaction (PCR) amplification mixture contained 20 ng of genomic DNA, 67 mM Tris-HCL, 16 mM ammonium sulfate, 100 μM of each of the four dNTPs, 2 mM MgCl₂, 1 ng of BSA, 0.2% mercaptoethanol, 0.2 μM of each primer, and 0.2 units of Taq polymerase (Tiangen Biotech, China). The amplification was carried out using 1 cycle of 3 min at 94°C, 35 cycles of 30 s at 94°C, 45 s at 50 to 55°C, 45 s at 72°C, and 1 cycle of 10 min at 72°C as modified by Du et al. (2017). The PCR products were visualized on 2% agarose gels. Templates were sequenced in the reverse direction with five reverse primers. Sequencing was performed using BigDye Terminator V 3.1 (Applied Biosystems, USA) following the manufacturer's instructions. Samples were run on an ABI 3130 Genetic Analyzer (Applied Biosystems, USA). Sequences were checked using Chromas V 2.2 and only those of high quality giving single peaks were used in subsequent analysis. DNA sequence alignment was carried out using Clustal W in MEGA V6.0 (Tamura et al., 2013) with manual adjustments. Chloroplast sequences have been deposited in GenBank under the accession numbers of KY751544-KY751561.

Twenty-seven nSSRs primers reported by Dow et al. (1995), Steinkellner et al. (1997), Kampfer et al. (1998), Ueno and Tsumura (2008), Ueno et al. (2008, 2009), and Durand et al. (2010), were tested for their ability to identify polymorphic SSR (see Supplementary Table 3 for details). Each pair of primers was initially amplified for three individuals (one for each oak species), and the PCR products were visualized on 2% agarose gels. The successfully amplified templates were sequenced to confirm the existence of repeated motifs and to determine the numbers of repeats. Each primer revealing SSR motifs was then amplified for 12 individuals (four individuals for each oak species) to test for polymorphism in simplex PCR reactions using the M13-tail technique (Schuelke, 2000). In brief, three primers were synthesized for each genotyping experiment: a 5′ M13-tailed forward primer, a reverse primer and a fluorescently labeled M13 primer carrying a FAM, HEX, TAMRA or ROX label (Sangon Biotech, China). The PCR reaction mixture contained 1 × Taq buffer, 0.2 mM dNTPs, 10–20 ng template DNA, 1.6 pmol of the reverse primer, 1.6 pmol of a single fluorescently labeled M13 primer, 0.4 pmol of the forward primer and 1 U Taq polymerase (Tiangen Biotech, China). The following conditions were used for PCR amplification: 8 min initial denaturation at 95°C and followed by 30 cycles of 30 s denaturation at 95°C, a 30 s annealing step at 56°C, a 30 s elongation step at 72°C, and a final extension step at 72°C for 8 min as modified by Xu et al. (2013). nSSRs with a satisfactory number of alleles (>5) and good amplification quality were selected. In total, 19 pairs of nSSRs primers (GOT021, PIE227, FIR026, QmC00716, POR017, QmC00196, QmC00963, GOT011, MSQ13, FIR015, PIE163, WAG066, DN950726, Qmc00932, DN950446, Qmc02052, PIE271, WAG068, and CcC00063; see details in Supplementary Table 3) were selected for genetic diversity analysis of all individuals. 0.5 μL samples of the PCR products obtained, combined with each of the four different fluorescently labeled primers, were added to 10 μL formamide and 0.5 μL of LIZ standard (Applied Biosystems, USA) and the mixture was analyzed using an ABI 3730 Prism Genetic Analyzer (Applied Biosystems, USA). Alleles were scored using the software GeneMarker V2.2 (Softgenetics, USA). Each genotype was checked by two readers, as described by Guichoux et al. (2011b). When an allele was found by only one reader, the inconsistency was classified as type A; when the two readers detected different alleles for the same sample, the inconsistency was classified as type B. The readers tried to determine a consensus genotype. If consensus result could not be achieved, the data were treated as missing.

A total of 3,344 bp from five cpDNA fragments were sequenced and a matrix of combined sequences was constructed. Twenty substitutions and six indels were detected in 343 individuals and were combined into seven chloroplast haplotypes (H1-H7) (Supplementary Table 5). H1 was the most common haplotype and was found in most Q. aliena (175 out of 178) and Q. dentata (47 out of 48), and in the single individual with intermediate morphology. H2 was the dominant haplotype in Q. variabilis, being found across all sites. H3 and H7 were restricted to Q. variabilis at site YSD. H4, H5, and H6 were only found at site XSL in Q. aliena, Q. dentate, and Q. variabilis, respectively (Table 1, Figure 1).

These haplotypes were clustered into two lineages separated by at least 19 mutations (Figure 1). H1, H4, and H5 formed a lineage containing the white oaks (Q. aliena and Q. dentata) and H2, H3, H6, and H7 formed the second lineage, that of Q. variabilis (Figure 1). Within the white oak lineage, the relationship between H1 and H5 was more distant with five mutation steps separating them; H4 differed from H1 by one mutation step. In contrast, in the Q. variabilis lineage, the four haplotypes were closely related: H3 and H6 were separated by one mutation from H2; H7 differed from H2 with two mutation steps.

The total haplotype diversity was high (H_T = 0.51) and more than ten times greater than the genetic diversity within populations (H_S = 0.03; Table 2). The N_ST value was significantly higher than the G_ST value (p < 0.05), suggesting the existence of phylogeography structure (Table 2). The cpDNA diversity of Q. variabilis was the highest among the three oak species, followed by that of Q. dentata and Q. aliena (Table 2). The results of AMOVA among the three species revealed that 99.83% of the variation could be attributed to among-species diversity and 0.17% to within-species diversity (Table 3). The variation between Q. aliena and Q. variabilis was similar to that between Q. dentata and Q. variabilis (Table 3) which was nearly fixed between species. In contrast, the variation between Q. aliena and Q. dentata was largely harbored within species (Table 3). In the intraspecies analysis, AMOVA results indicated a high level of genetic variation within sites while negligible variation among sites (negative estimates for Q. aliena and Q. dentata; 0.01 for Q. variabilis, Table 3).

Two nSSRs loci (QmC00196 and WAG066) presenting null alleles detected by Micro-Checker and 18 individuals showing poor amplifications were excluded from further analysis. In total, genetic diversity of 414 individuals was estimated based on 17 loci. All 17 nSSRs loci were polymorphic with the mean number of different alleles per locus (N_a) varying from 9.18 to 12.29 (Table 2). In contrast to the results of cpDNA analysis, nSSRs data indicated that the genetic diversity of Q. variabilis was the lowest among the three oak species (H_e = 0.63; Table 2). The observed heterozygosity (H_o) ranged from 0.59 (Q. variabilis) to 0.70 (Q. dentata) (Table 2). Within species, the N_a of Q. variabilis ranged from 5.06 (site TL) to 7.59 (site XSL). The expected heterozygosity (H_e) ranged from 0.55 (site TZF) to 0.64 (site XSL), and the H_o from 0.53 (site TZF) to 0.62 (site XSL) (Table 1). The N_a of Q. aliena ranged from 7.18 (site TZF) to 9.47 (site XSL). The highest H_e was 0.68 at site ZL and the highest H_o was 0.69 at sites ZL and XSL (Table 1). Both N_a (9.53) and H_e (0.71) for Q. dentata had their highest values at site ZL and the H_o ranged from 0.68 (site XSL) to 0.77 (site TZF) (Table 1). The results of the overall AMOVA suggested that most of the variation occurred within species (79.4%; F_ST = 0.21), and only 20.6% of the variation was due to interspecies differences (Table 3). AMOVA results of the three species pairs (Q. aliena and Q. dentata; Q. aliena and Q. variabilis; Q. dentata and Q. variabilis) revealed a high level of genetic variation within species (Table 3). The intraspecies analysis indicated a high level of genetic variation within sites (Table 3).

The likelihood of classifying the data increased greatly from K = 1 to K = 3 and then increased slightly from K = 3 to K = 4, where it reached a plateau (Supplementary Figure 2). The ΔK (ΔK = 4157.04 when K = 2, Supplementary Figure 2) indicated that two was the optimal number of clusters (K = 2). However, the statistics also gave some support for K = 3 (Supplementary Figure 2). When K equaled to 2, one cluster corresponded to Q. variabilis (red), and the other one included Q. aliena and Q. dentata (green). When K equaled to 3 (ΔK = 1280.58), each species was classified into a distinct cluster with red, green and blue representing Q. variabilis, Q. aliena and Q. dentata, respectively (Figure 2). The PCoA result showed three distinct clusters, corresponding to Q. aliena, Q. dentata, and Q. variabilis. The first axis separated Q. variabilis from a cluster of Q. aliena and Q. dentata, and the second axis separated Q. dentata from Q. aliena and Q. variabilis (Figure 3).

If the threshold value of Q was set to 0.8, when K = 2, no individual was identified as a hybrid (Figure 2); when K = 3, seven individuals were hybrids between Q. aliena and Q. dentata with Q-values ranging from 0.24 to 0.46. These included the individual with intermediate morphology at site TK, one morphologically inferred Q. aliena at TL, two (one Q. aliena and one Q. dentata) at XSL and three (two Q. aliena and one Q. dentata) at site ZL. Hence, hybrids accounted for 1.7% of the total samples with only one individual being inferred to be an F1 hybrid (Q-value was 0.464, at site ZL) (Figures 2, 4).

If the Q-value was set to 0.9, when K = 2, only one individual at site TL was identified as a hybrid between Q. aliena and Q. variabilis, with a Q-value of 0.163. When K = 3, 14 individuals, including the individual with intermediate morphology at site TK, were considered to be hybrids. Among them, one hybrid between Q. aliena and Q. variabilis and six hybrids between Q. aliena and Q. dentata were identified with Q-values from 0.1 to 0.2 and another seven hybrids between these two species were detected with Q-values from 0.2 to 0.5; only one of these seven was considered to be an F1 hybrid (Q-value of 0.46, at site ZL) (Figures 2, 4).

A comparison of different models calculated by Migrate-n based on Bayesian factor values showed that model 2, with three population sizes and two migration rates (from Q. aliena to Q. dentata and from Q. dentata to Q. aliena) was the best model for evaluating the historical gene flow (Chen et al., 2017) among the three oak species. This model verified the occurrence of bidirectional historical gene flow between Q. aliena and Q. dentata (Table 4, Q. aliena→Q. dentata, M = 56.3; Q. dentata→Q. aliena, M = 57.0).

The analysis using cpDNA reveals low intraspecies diversity and high interspecies diversity in the oak complex on Shangfang Mountain, whereas both inter- and intraspecific diversity revealed by nSSRs data are high (Tables 1, 2). Sharing of haplotype H1 between Q. aliena and Q. dentata (Table 1) suggests very recent speciation and incomplete lineage sorting or introgression of H1 from one species to another. Sharing of cpDNA haplotypes has also been reported in other sympatric distributed Quercus species (Whittemore and Schaal, 1991; Dumolin-Lapègue et al., 1999; Kanno et al., 2004; Okaura et al., 2007). The inter- and intraspecific diversity revealed by nSSRs data (Tables 1, 2) are consistent with other oak species studied on a fine scale (Antonecchia et al., 2015).

For cpDNA, AMOVA results show high genetic differentiation among the three oak species and between the two sections of oak species. Variation between the combination of Q. variabilis and either Q. aliena or Q. dentata is completely fixed between species whereas most of the variation between Q. aliena and Q. dentata is within species (Table 3). For nSSRs, The AMOVA results show the majority of variation occurs within species (Table 3). The inconsistencies between the AMOVA results for cpDNA and nSSRs markers for all species, Q. aliena & Q. variabilis and Q. dentata & Q. variabilis agree with those reported by Petit et al. (2005), which can be explained by the different modes of inheritance of cpDNA and nuclear DNA. cpDNA is maternally inherited dispersing limited gene flow through seeds while nuclear DNA is bi-parentally inherited and can mediate extensive gene flow via both seeds and pollen. However, the AMOVA results for Q. aliena and Q. dentata indicate that interspecific differences are lower in cpDNA (Mean F_ST = 0.03 for cpDNA; F_ST = 0.11 for nSSRs) which may due to large proportion of individuals sharing haplotype H1 between the two species. Our results indicate chloroplast haplotypes are less species-specific than nuclear markers, probably because chloroplast DNA is more susceptible to drift than is nuclear DNA (Whittemore and Schaal, 1991; Dumolin-Lapègue et al., 1999; Petit and Excoffier, 2009).

For each species the intraspecific variation within sites is very high, with almost no genetic differentiation among sites within species for either cpDNA or nSSR markers (Table 3). These results indicate that even the rugged mountain ranges in this region have not impeded extensive gene flow mediated by seeds and pollen among sites which has helped to maintain the species integrity. Human modifications of the landscape such as the building of ancient temples, caves and bridges may have fragmented the forest, decreasing habitat size and increasing isolation between habitats and populations (Young et al., 1996; Su et al., 2003). However, for outcrossing and wind-pollinated tree species such as oaks, the fragmented spatial organization may not impede the genetic exchange via pollen among sites (e.g., reviewed by Herrera-Arroyo et al., 2013), as demonstrated here.

The genetic assignment by STRUCTURE using nSSRs shows that two clusters are optimal, reconfirming the taxonomic classifications of these three oak species. The two white oaks were classified into one cluster as they are more closely related to each other than to Q. variabilis. Given that K = 3 is more biologically representative of the oak species sampled, we also discuss the admixture results for three clusters. With three genetic clusters, Q. aliena, Q. dentata, and Q. variabilis are all clearly separated (Figure 2). Low frequency of bi-directional gene flow between Q. aliena and Q. dentata were visually detected in STRUCTURE bar plots (Figure 2). These results are also confirmed by the PCoA results (Figure 3). The existence of reproductive barrier(s) between these two closely related species may be one of the explanations for the limited gene flow between them. At our study site, the flowering periods of Q. aliena and Q. dentata are asynchronous but they occur in fairly rapid succession. Moreover, the relative abundance of parental species and pollen discrimination can contribute to pre-pollination barriers, limiting interspecific pollen transfer, since abundant species are proportionally less involved in hybridization than minority species (Lepais et al., 2009, 2013; Salvini et al., 2009). The relative abundance of Q. aliena is more than twice of that of Q. dentata in our study sites, so interspecific gene flow may be more frequent in the core ranges of Q. dentata (Figure 2 shows more hybrids at site ZL than at site XSL). However, because neither seeds nor seedlings have been sampled in this study, we cannot rule out the existence of postzygotic isolation mechanisms such as seed abortion, hybrid lethality, weakness and sterility.

It should be noted, however, the high number of immigrants per generation 4Nm (calculated by θM, higher than one) detected by Migrate-n in Q. aliena and Q. dentata, which suggests the existence of relatively high level of symmetric historical gene flow between them (Table 4). These results are consistent with genetic differentiation based on cpDNA with the mean F_ST being 0.03 and 0.11 for cpDNA and nSSRs, respectively (Table 3). This will have to be borne in mind when attempting to interpret the species boundaries.

For either K = 2 or K = 3, only one individual was identified as a hybrid between Q. aliena and Q. variabilis (Q = 0.16) and no hybrids between Q. dentata and Q. variabilis were detected (Figure 2), a finding which suggests that there is almost no gene flow either between Q. aliena and Q. variabilis or between Q. dentata and Q. variabilis. Nearly no gene flow was detected between Q. variabilis and either Q. aliena or Q. dentata by Migrate-n (Table 4). Combining the cpDNA and nSSRs results, clear species boundaries between the two sections of oak on Shangfang Mountain were detected.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.