Sampling of endangered oak species Quercus austrocochinchinensis and Q. kerrii was granted and supported by National Forestry Bureau of China and Local National Nature Reserves.

In total, 57 and 36 individuals with typical traits of Q. austrocochinchinensis and Q. kerrii, respectively, were used and 15 morphological intermediates were included in this study. The 108 individuals were sampled in five populations, A–E (Figure 1 and Table 1). According to a previous investigation (Song et al., 2015), Q. austrocochinchinensis trees can only be found in populations D and E. Of these, E was considered as a pure Q. austrocochinchinensis population (Song et al., 2015) and D is located in the contact zone which contains both Q. kerrii and Q. austrocochinchinensis trees. Six sub-populations, D1–D6, were sampled within D population regions. Two putative Q. kerrii purebred populations were sampled from populations A and B. Population C is a putative hybrid population with morphological intermediates, but trees with typical traits of Q. austrocochinchinensis cannot be found in or close to region C.

Prior to the experiment, all individuals and their voucher specimens were carefully inspected and identified based on key morphological diagnostic traits, e.g., shape of cupule and trichomes on leaf abaxial surface. The results were used as species information for the studied samples in later analyses. Leaf tissues used for DNA extractions were collected from each individual and dried instantly using silica gel. All vouchers specimens of each tree were stored at the herbarium of the Shanghai Chenshan Botanical Garden (CSH).

In a previous study, we found that the two macro-morphological features leaf apex shape and leaf length-to-width ratio are the key diagnostic features that can be applied to identify Q. kerrii and Q. austrocochinchinensis and even to assign their hybrids (Song et al., 2015). To link the genetic pattern to the morphological features, we used ratios instead of absolute lengths to measure leaf morphological traits for the two parental species and the intermediates. The two ratios were leaf length-to-width ratio and leaf apex length-to-width ratio. The measuring method of the two ratios is shown in Figure 2. At least five leaves were measured from each voucher specimen to compute an average for the drawing of the scatter plot.

Total genomic DNA for each sampled individual was extracted from the silica gel-dried leaf tissue using the modified CTAB method (Huang et al., 2000). The DNA quality was checked by loading DNA on a 1.0% agarose gel, and the DNA concentration of each sample was measured using a TBS-380 Fluorometer (Turner BioSystems Inc., Sunnyvale, CA). The samples were screened for 11 SSR markers. Of these loci, Qk15874, Qk17139, Qk17611, and Qk20944 have been developed for Q. kerrii from transcriptome data (An et al., 2016), QmC00693, QpZAG9, QpZAG16, QpZAG36, QpZAG110, CG371 were from the non-coding region (Steinkellner et al., 1997; Ueno et al., 2008; Tong et al., 2012), and CR627959 was predicted in the coding region of Cys-3-His zinc finger protein in nuclear genome (Ueno and Tsumura, 2008). All of these loci show a relatively high degree of transferability within the genus Quercus and an adequate degree of polymorphism of the studied taxa. For the 11 SSR primer pairs, 5′ end of forward primers were labeled with fluorescent dye tags (6FAM or HEX or ROX) (Sangon, Shanghai, China). The PCR reactions were performed in 20 μl reaction volume containing 10 μl TIANGEN PCR Master Mix (TIANGEN, Beijing, China), 0.3 μl/L of each primer (10 mM), and 20 ng genomic DNA; PCR reactions were performed as follows: 5 min initial denaturation at 94°C, 35 cycles of 40 s at 94°C, 30 s at 55°C, 1 min elongation at 72°C, and 7 min extension at 72°C. Finally, Gene-Scan-500 LIZ size standard (Applied Biosystem™) was added to all the samples before loading on an automated sequencer ABI 3730 (Applied Biosystem™). The final step was performed by a private professional commercial lab (Shanghai Majorbio Bio-pharm Technology Co., Ltd, Shanghai, China).

The AFLP method was performed following the protocol by Vos et al. (1995), with minor modifications. Briefly, 500 ng of genomic DNA were double-digested using 5 U of EcoRI and 2 U of MseI (New England BioLabs). The digestion mixtures were incubated at 37°C for 3 h and the digested mixture was then incubated at 70°C for 10 min to denature the enzymes. Subsequently, 4 μL of digested DNA were added to 16 μL of ligation mix containing 2 U T4 DNA ligase (New England BioLabs), 5 pmol EcoRI, and 50 pmol MseI adaptor. The mixture was incubated at 10°C for 14 h and then denatured at 70°C for 10 min. The ligated DNA samples were diluted 5-fold with double sterile water. Pre-selective amplification reactions were carried out using EcoRI-A (5′-GACTGCGTACCAATTCA-3′) and MseI-C (5′-GATGAGTCCTGAGTAAC-3′) in a 50 μL volume containing 1.5 mmol/L MgCl₂, 200 μmol/L of each dNTP, 1.25 μmol/L of each primer, and 0.6 U rTaq DNA polymerase (Takara Biotechnology, Dalian, China) under the following cycle: 3 min at 72°C, 30 cycles of 30 s denaturing at 94°C, 30 s annealing at 56°C, 1 min extension at 72°C, and a final extension for 5 min at 72°C. After a 1:20 dilution of pre-selective PCR products, seven primer combinations were performed in selective amplification (EcoRI-ACG/MseI-CAC, EcoRI-AGG/MseI-CAA, EcoRI-AGG/MseI-CAC, EcoRI-AGG/MseI-CAG, EcoRI-AGG/MseI-CTA, EcoRI-AGG/MseI-CTC, EcoRI-AGG/MseI-CTT). The EcoRI primers were fluorescently labeled with 6-FAM. These primer pairs were chosen because they generated clear and fewer bands (thus decreasing the risk of fragment non-homology) with sufficient variability in preliminary tests. Selective PCRs were carried out in a 20 μL volume containing 2.5 mmol/L MgCl₂, 200 μmol/L of each dNTP, 1.25 μmol/L of each primer, and 0.2 U of rTaq DNA polymerase (Takara Biotechnology, Dalian, China) and under the following cycle: 3 min at 94°C, 9 cycles of 40 s at 94°C, 30 s at 65–57°C touchdown (reducing the temperature at 1°C per cycle), 15 min at 72°C, 20 cycles of 40 s at 94°C, 30 s at 56°C, 1.5 min at 72°C, and a final extension for 7 min at 60°C. The PCR products were 10-fold diluted and mixed with Gene-Scan-500 LIZ size standard (Applied Biosystem™); products from each primer combination were loaded separately on an automated sequencer ABI 3730 (Applied Biosystem™) by the same commercial service provider mentioned above.

Raw data of SSR and AFLP samples were collected and analyzed using GeneMarker®v2.2.0. The samples with low quality of size calibrations or peaks were excluded from allele calling. The allele sizes of SSR were read manually and MicroChecker v 2.2.3 (Van Oosterhout et al., 2004) was used to check for potential errors. For AFLP, the variable fragments in the size range 50–500 base pairs (bp) were manually scored as present (1) or absent (0). We only considered fragments with similar fluorescence profile and intensities across the samples to maximize the probability of homology.

A set of statistical tests on SSRs, including allelic richness (Na), allele frequency distribution, Shannon's Information Index (I), Observed and Expected Heterozygosity (H_O and H_E), were performed by GenAlEx version 6.5 (Peakall and Smouse, 2012), using all the individuals of Q. kerrii and Q. austrocochinchinensis. In addition, Unbiased Expected Heterozygosity (uH_E) was estimated as (2N/(2N-1)) * H_E. Fixation Index (F_IS), or inbreeding coefficient was estimated as (H_E–H_O)/H_E (Nei and Li, 1979). We used GenePop v 4.2 (Rousset, 2008) to test the departure from Hardy-Weinberg equilibrium (HWE) (heterozygote deficiency or excess) for each locus of the 11 SSR loci, and to test for homogeneity of alleles distributions between species. We also counted the number of private alleles for each species.

For AFLP, percentage of polymorphic loci, unbiased estimates of genetic diversity (H_j, analogous to H_E), and differentiation statistics were calculated using the AFLP-SURV v. 1.0 software (Vekemans et al., 2002). With this software, allelic frequencies at AFLP loci were calculated from the observed frequencies of fragments using the Bayesian approach proposed by Zhivotovsky (1999) for diploid species. A non-uniform prior distribution of allelic frequencies was assumed with its parameters derived from the observed distribution of fragment frequencies among loci. These allelic frequencies were used as the input for the analysis of genetic diversity within and between samples following the method described in Lynch and Milligan (1994).

The comparison of F_ST and Q_ST provides a basis to distinguish neutral from adaptive divergence (Leinonen et al., 2013). To investigate genetic differentiation between the two species and among populations, F_ST values were measured on AFLP and SSR data using AFLP-SURV v. 1.0 (Vekemans et al., 2002) and GenePop v. 4.2 (Rousset, 2008), respectively. The parameter Q_ST estimates the among-population proportion of the total additive genetic variance of a genetic quantitative trait. If Q_ST = F_ST, trait divergence among populations could have been driven only by genetic drift. If Q_ST > F_ST, the populations are likely to have been caused by directional selection, and if Q_ST < F_ST, there is evidence for uniform selection or stabilizing selection across the populations. Because this is not possible for the breeding designs for this study, we used P_ST, a proxy for Q_ST, to compare with F_ST (Leinonen et al., 2013). In this study, P_ST was calculated between all the Q. austrocochinchinensis and Q. kerrii individuals, using the equation of σGB2/(σGB2+2σGW2), where σGW2 and σGB2 are within- and among-population components of variance. Putative hybrids were excluded from F_ST and P_ST estimation.

In addition, F_ST for each locus was also estimated. For AFLP, allele frequencies of Q. austrocochinchinensis and Q. kerrii were calculated separately. We calculated F_ST values between 57 for Q. austrocochinchinensis and 36 for Q. kerrii individuals for all polymorphic loci via the formula F_ST = 1-H_S/H_T (Nei, 1973). The variable H_S represents average within-population heterozygosity and H_T represents expected heterozygosity for the total population; H_S and H_T were calculated using the following formula: H_S = 1/2(2p₁q₁+2p₂q₂) and H_T = 1/2(p₁+p₂)*(q₁+q₂), with q = 1-p. For SSR markers, F_ST for each loci was calculated using GenePop v. 4.2 (Rousset, 2008). Finally, the frequency distribution of F_ST values for both markers were plotted in a histogram.

To detect the outlier loci under selection of an AFLP dataset (781 loci), program BayeScan v. 2.1 (Foll and Gaggiotti, 2008; Fischer et al., 2011) was used to identify candidate loci under natural selection across all 10 populations. A threshold value for determining loci under selection was evaluated in accordance with Jeffreys (1961) interpretation, which is a logarithmic scale for model choice as follows: log₁₀PO > 0.5 (substantial), log₁₀ PO > 1.0 (strong), log₁₀ PO > 1.5 (very strong), and log₁₀PO > 2.0 (decisive support for accepting a model) (Fischer et al., 2011). We employed a threshold of log₁₀ PO > 2.0 for the rejection of the null hypothesis in each of the conducted tests. BayeScan analysis was conducted with a burn-in of 50,000 iterations, a thinning interval of 50, and a sample size of 5,000. The number of pilot runs was kept at 20, with a length of 50,000 each. The SSR loci were not used to calculate outliers because of their limited number.

To cluster individuals into genetically distinct groups, a Bayesian clustering approach was employed using STRUCTURE v. 2.3.4 (Pritchard et al., 2000; Falush et al., 2003), without consideration of sampling information. We adopted the admixture model with correlated allele frequencies (Lepais et al., 2009; Zalapa et al., 2009). No prior knowledge of the species was included in the analyzed data sets. To determine the optimal number of groups (K), we ran STRUCTURE, with K varying from 1 to 10 and with 10 runs for each K value. The ΔK was calculated using the mean log-likelihood for each K according to Evanno et al. (2005). For SSR and AFLP data sets, each run was performed for 100,000 Markov Chain Monte Carlo (MCMC) repetitions with a burn-in period of 50,000.

The admixture coefficient (q-value) generated from STRUCTURE was used to classify individuals into purebred and hybrids, using a threshold q-value of = 0.1, where samples with q-values < 0.1 or > 0.9 were classified as purebred and those with q-values between 0.1 and 0.9 as hybrids, including F₁ and backcrosses (Vähä and Primmer, 2006; Lepais et al., 2009). The F₁ hybrids result in q-values = 0.5, but the coefficient of backcrosses would be biased toward one of the parental species and produce q-values between 0.5 and 0.9 (Lepais et al., 2009). Taking errors into consideration, individuals with 0.6 < q-values < 0.9 were recognized as backcrosses.

As the SSR loci were tested to violate the HWE assumption (see result), we also used the program InStruct (Gao et al., 2007), an alternative software of STRUCTURE, to ensure that the results obtained from STRUCTURE were reliable. The HWE within loci for co-dominant markers is not compulsory for the model in InStruct. The optimal K (number of clusters) value was selected according to Deviance Information Criteria (DIC) and is presented in the result section. Mode 2 was selected to run 50,000 MCMC repetitions with 10,000 burn-in periods.

The NewHybrids software v. 1.1 beta, employing a Bayesian analysis (Anderson and Thompson, 2002; Anderson, 2008), was used as a second method to confirm the presence of hybrids in our dataset. It calculates the posterior probability that sampled individuals fall into one set of hybrid categories (Parent K, Parent A, F₁, F₂, BC to Parent K; BC to Parent A, thus covering parents and two generations of offspring). All the 11 SSR loci were analyzed using this software. For AFLP data, NewHybrids will assign the individuals more accurately, using the loci with high differentiation at interspecies level. Therefore, we filtered the AFLP loci according to F_ST. We applied three approaches to compare the results obtained by NewHybrids, using 100 loci with highest F_ST value at interspecies level, 279 loci (F_ST ≥ 0.1), and 450 AFLP loci (F_ST ≥ 0.03), respectively, for analysis. A burn-in period of 7,500 MCMC repetitions was defined and 10,000 iteration were run thereafter.

Principal coordinate analysis (PCoA) aims to visualize similarities or dissimilarities of individual data based on a distance matrix. This method is also an alternative algorithm without any assumptions about the population genetic model. For AFLP and SSR data, the pairwise Euclidian distance matrix was constructed, and the first two principal co-ordinates were visualized by GenAlEx v. 6.5 (Peakall and Smouse, 2012).

Morphological measurements showed that the two species could be distinguished using a combination of two ratio values, which were leaf length-to-width ratio (LR) and leaf apex length-to-width ratio (LAR) (Figure 2). We found that all Q. kerrii individuals were grouped together, with an average of LR = 2.45 (1.88–3.02) and LAR = 0.78 (0.23–1.36). Of the sampled Q. austrocochinchinensis individuals, only the samples from D2 and E were totally separated from Q. kerrii samples. They had larger LR and LAR values compared to other populations, with an average of LR = 3.71 (2.79–5.32) and LAR = 2.06 (1.50–3.20). Morphological variation was higher in Q. austrocochinchinensis individuals than in Q. kerrii specimens, with standard deviations (σ) for Q. austrocochinchinensis and Q. kerrii of 0.72 and 0.25 for LR and 0.64 and 0.25 for LAR, respectively.

The putative hybrids of D5 within the hybrid zone had morphological values similar to Q. kerrii, with an average of LR = 2.52 (2.22–2.74) and LAR = 0.82 (0.47–1.13). However, the putative hybrids of population C had a sharper leaf apex with LAR = 1.54 (0.86–2.27), although the leaves were as broad as those of Q. kerrii with LR = 2.77 (2.35–3.47). In addition, there were seven Q. austrocochinchinensis individuals from D3 and D6 with leaf morphological traits similar to those of Q. kerrii.

In total, 151 alleles were obtained from the 11 SSR loci used in this study. Some SSRs were highly variable, e.g., CG371, QpZAG16, and QpZAG110, containing more than 20 alleles. These SSRs also had higher Observed Heterozygosity (H_O) and Expected Heterozygosity (H_E) (Table 2) than other loci. Although most frequent alleles were shared by Q. austrocochinchinensis and Q. kerrii, some loci showed greater variation in allele frequency (Figure S1). Mean Expected Heterozygosity (H_E) across all loci in Q. austrochichinensis was 0.706 and higher than that in Q. kerrii with 0.595. Species-specific alleles were found at several loci, especially rare alleles restricted to Q. austrocochinchinensis (Table 3). Frequency of the private alleles was 34.4 and 12.0% in Q. austrocochinchinensis and Q. kerrii, respectively. All SSR loci except Qk17611 significantly deviated from HWE (p < 0.01). Heterozygote deficiency was found at all loci except Qk17611 and QpZAG36 (p < 0.01). The estimation of genetic differentiation between Q. austrocochinchinensis and Q. kerrii was low across all SSR loci (F_ST = 0.138). The differentiation among population within species was also low for both species, with F_ST values of 0.130 and 0.042 for Q. austrocochinchinensis and Q. kerrii, respectively.

Application of the seven AFLP primer pairs to 108 individuals resulted in 859 loci, of which 781 were polymorphic. The levels of diversity within each species, either at the population (H_W) or the whole sample level (H_t), were very similar (Table 4). Genetic differentiation between Q. austrocochinchinensis and Q. kerrii across the 781 markers was low (F_ST = 0.095, p < 0.01) (Table 4). Similarly, the differentiation among populations within species was also low, with F_ST = 0.0246 (p < 0.01) for Q. austrocochinchinensis and F_ST = 0.0667 (p < 0.01) for Q. kerrii. The F_ST distribution of 781 AFLP alleles followed an L-shaped curve, which suggested that most of the alleles (471 of 781, 60.3%) had lower F_ST values (F_ST < 0.1). Only a few alleles exhibited higher F_ST values (Figure S2). The P_ST values for all individuals of the two species, except for the putative hybrids, were larger (0.338) than F_ST (0.095 for AFLP, 0.082 for SSR). By using BayeScan, 3 out of 781 AFLP loci (0.38%) were identified as outlier loci under directional selection at a threshold of log₁₀ PO > 2 (posterior probabilities higher than 0.99) (Figure S3).

The most likely number of clusters (K) of STRUCTURE analysis can be inferred by the peak of ΔK. We found that K = 2 was the optimal K value for both AFLP and SSR (Figure S4), representing two species, respectively. Following ten independent STRUCTURE runs with K = 2, individuals morphologically identified as Q. kerrii were assigned to one cluster with high probability, whereas those morphologically identified as Q. austrocochinchinensis were assigned to the other cluster with similarly high probability. Therefore, these two clusters were determined to represent Q. kerrii and Q. austrocochinchinensis, respectively. The STRUCTURE analysis results of SSR and AFLP are illustrated in Figures 3A,C, respectively, but the results are very different. Among individuals morphologically identified as Q. kerrii, both AFLP and SSR recognized populations A and B are Q. kerrii purebred population, but the results of the remaining populations are largely contradictory. For example, SSR results assigned most of the individuals in populations/subpopulations C, D1, D3, D4, D5, D6, E, and two individual of D2 to the Q. kerrii cluster, and the remaining individuals in D2 to another cluster, regardless of morphological identification. Based on the AFLP results, among the individuals identified as Q. kerrii in populations C, D1, D3, and D4, the mean proportion of Q. kerrii was 0.575 (0.412–0.941). The two morphological intermediate putative hybrid populations C and D5 contained a high proportion of Q. austrocochinchinensis germplasm (with a mean value of 0.901 [0.677–0.992]). Among the individuals identified as Q. austrocochinchinensis, in populations D1, D2, and E, the mean proportion of Q. austrocochinchinensis germplasm was 89.1% (with the lowest value, 21.7%, found in individual Ca44), but mean proportion was considerably lower in D3 with 34.5% (1.1–71.3%). Given that we considered individuals with q-values between 0.6 and 0.9 as backcrosses, we did not identify any backcrossed individuals based on SSRs. However, for AFLP, 25 out of 108 (23.1%) samples were recognized as backcrosses.

The software InStruct (Gao et al., 2007) was used to verify the reliability of SSR results generated by STRUCTURE. The result was highly consistent to the results obtained from STRUCTURE. The optimal K value selected by DIC was also 2, similar to the best K inferred by STRUCTURE. The D2 subpopulation with pure Q. austrocochinchinensis was still easily distinguished from the analysis (Figure S5).

Analysis of SSR data using NewHybrids did not obtain results in agreement with the morphological identity of Q. kerii and Q. austrocochinchinensis individuals examined. Most individuals previously identified as Q. austrocochinchinensis, based on morphological characters from subpopulation D2, were still assigned to Q. austrocochinchinensis, with high posterior probabilities (>0.974), except for one individual (Ca21). The remaining individuals identified as Q. austrocochinchinensis in D1, D3, D6, and E were assigned to either Q. kerrii or backcrosses to Q. kerrii, as well as the morphological intermediates individuals in Population C and D5. Among the 36 individuals identified as Q. kerrii, most were assigned to Q. kerrii, with mediate to high probabilities, two (CK10, CK232) were assigned to backcrosses to Q. kerrii, with probabilities of 0.935 and 0.802, respectively (Figure 3B).

The assignment of individuals to certain genotypes using AFLP data better fits the morphological identification, but is still very messy. The results obtained from 100, 279, and 450 loci of AFLP were similar in terms of assigning the individuals of Q. austrocochinchinensis and Q. kerrii from the D region, but the results were different in terms of assigning the individuals in Q. kerrii purebred populations A and B, as more backcrosses to Q. kerrii were found when adding more loci with low divergence of F_ST to the analysis (Figure 3, Figures S5B,C). To accurately estimate hybrid classes, alleles with large divergence of F_ST are required, or misclassification of backcrosses and purebred parental individuals are very likely to happen (Vähä and Primmer, 2006). Therefore, the NewHybrids results obtained from 100 AFLP loci with highest F_ST value were more reliable and sensitive on assigning the genealogical classes.

Based on morphological features, 57 individuals were identified as Q. austrocochinchinensis in populations D and E, of which 30 were assigned to pure Q. austrocochinchinensis, one was assigned to F₂ hybrids, six to backcrosses to Q. austrocochinchinensis, and four to backcrosses to Q. kerrii with high probabilities. The remaining genotypes were mixtures of certain amounts of Q. kerrii, backcrosses to both parental species, and F₂ (Figure 3D). Among the 15 morphological intermediate individuals of populations C and D5, seven were assigned to Q. austrocochinchinensis purebreds, seven to backcrosses to Q. austrocochinchinensis with high probabilities, and one was an admixture of Q. austrocochinchinensis and its backcrosses. The morphologically identified individuals of Q. kerrii in A and B populations were predominantly assigned to pure Q. kerrii. The remaining genotypes were mixtures, with few admixtures of certain amounts of Q. kerrii and backcrosses to Q. kerrii. In D3 and D4, no purebred parental individuals existed, three individuals were assigned to F₂ with high probability. The genotypes the remaining individuals were either backcrosses to both parental species or a wide spectrum of admixtures of the backcrosses to both parental species, F₁, F₂, and Q. kerrii purebreds. Of the six individuals of Q. kerrii in subpopulations D3 and D4, three were assigned to backcrosses to Q. austrocochinchinensis, one to backcross to Q. kerrii, and two were identified as F₂ (Figure 3D).

The PCoA results of the AFLP and the SSR data are shown in Figures 4A,B, respectively. The PCoA results were in agreement with the results of the STRUCTURE analysis. Individuals of the two species were mostly separated for both AFLP and SSR data, which indicated that interspecific differentiation was stronger than intraspecific differentiation.

For the AFLP data, most Q. austrocochinchensis individuals in D2, D6, and E showed little mixing with Q. kerrii and were grouped together (left dashed line circle in Figure 4A). The Q. kerrii purebred populations A and B were grouped unambiguously (right dashed line circle in Figure 4A). Populations D1 and D3, which contained both species, were grouped between the two dashed line circles that represented each species. The morphological intermediates in C and D5 were grouped together in the middle position and biased toward Q. austrocochinchensis. For SSR analysis, most Q. austrocochinchinensis individuals in D2 were grouped together with smaller PC1 values (dashed line circle in Figure 4B). Other samples did not separate into distinct clusters. Still, the result of the PCoA seems to be more reliable, as it revealed more groups than that inferred by STRUCTURE, which only recognized two groups. Except for D2, the putative pure parental species and hybrids were mostly grouped together, but the resolution was not as high as the AFLP data.

Quercus austrocochinchinensis is a rare species with only two known distribution sites in China: Xi-Shuang-Ban-Na (XSBN) Nature Reserve in Yunnan province and Ba-Wang-Ling (BWL) Nature Reserve in Hainan Province (Huang et al., 1999; Song et al., 2015). Trees with intermediate morphological form between the two parental species have been found previously in XSBN, suggesting ongoing natural hybridization (Song et al., 2015). However, morphological intermediacy is not invariably associated with hybrids, as it can be a result of hybridization or phenotypic plasticity of the species in these areas (Rieseberg et al., 1993). Therefore, one major objective of this study is to confirm the possible hybridization and introgression in the two species.

Overall, our study revealed that there were fewer Q. austrocochinchinensis purebreds than previously expected based on morphological diagnostic traits. Although AFLP and SSR have different distinguishing power on the genotypes of examined samples, both markers revealed the presence of hybridization and introgression between the two species. The SSR data indicated that subpopulation D2 of Q. austrocochinchinensis has a unique germplasm composition and might be the only existing purebred. Further STRUCTURE analysis identified that only one individual (Ca78), located at contact zone subpopulation D1, is F₁. NewHybrids analysis is more sensitive to assign the genotype to different genetic categories and detected a general presence of backcrosses to Q. kerrii, both in morphologically identified Q. kerrii and Q. austrocochinchinensis and the intermediate populations, but no F₁, F₂, or backcrosses to Q. austrocochinchinensis. The SSR results indicate unidirectional introgression from Q. kerrii to Q. austrocochinchinensis.

For the AFLP data set, both Bayesian clustering approaches used (implemented in STRUCTURE and NewHybrids) detected an unexpected high number of backcrosses and hybrid genotypes. The threshold values and loci used on assigning individuals to different genetic categories are different in STRUCTURE and NewHybrids. Therefore, the two methods provided different percentages on F₁ and F₂ hybrids and backcrosses. STRUCTURE inferred the existence of F₁ hybrids and backcrosses, whereas NewHybrids inferred absences of F₁, but predominant backcrosses. STURECTURE is more efficient to evaluate the presence of hybrids in wild populations, whereas NewHybrids algorithm explicitly searches for hybrid and parental classes with assumption of two parental classes, which generally showed higher assignments accuracy than STRUCTURE (Marie et al., 2011). The NewHybrids result suggested that the formation of first-generation hybrids is less likely to occur than the interbreeding of hybrids with purebreds or with other hybrids. The occurrence of F₂ hybrids and the predominance of first-generation backcrosses to both parental species also reflect recent hybridization between the two species. Meanwhile, the genotypes of some individuals fell between F₂ hybrids, the first-generation backcrosses, pure Q. kerrii, and pure Q. austrocochinchinensis, with NewHybrids having no category available to assign them. Those individuals might in reality be second- or later-generation backcrosses or hybrids. This phenomenon is quite prominent in subpopulations D3 and D4, as no purebred individuals exist and the individuals are all hybrids, with varying percentages of backcrossing and parental types. Such evidence indicates that the two locations (D3 and D4) contain historical and ongoing gene flow between the two species.

The intermediate individuals in populations C and D5 are morphologically similar to Q. kerrii, mainly in terms of the shallow cupule, persistent trichomes on the leaf abaxial surface, and a relatively thick bark. However, their sharp leaf apex and leaf margin teeth more resemble Q. austrochochinensis (Song et al., 2015). The STRUCTURE results of the AFLP and SSR data on individuals of C and D5 were contradictory. The AFLP result indicated that individuals of C and D5 were Q. austrocochinchinensis purebred and backcrosses to Q. austrocochinchinensis, but the SSR suggested that all individuals in C and D5 were Q. kerrii purebreds. In the PCoA analysis, morphologically intermediate individuals from C and D5 were grouped together and located between two parental purebreds (Figure 4). Interestingly, C and D5 are geographically distant. Population D5 is located between regions with Q. kerrii and Q. austrocochinchinensis, but population C does not have any individuals with the typical features of Q. austrocochinchinensis; in addition, all the individuals are young trees, as the area is almost entirely occupied by farming land. The AFLP data revealed a high percentage germplasm of Q. austrocochinchinensis in C and D5. A similar situation was also found in population E and subpopulations D1 and D6, indicating that they are genetically “swamped” and that this “swamping” occurred recently, as a high percentage of germplasms of Q. austrocochinchinensis can still be detected.

A hybrid swarm is often characterized by a wide spectrum of phenotypic variation, the existence of backcrosses, and high genetic variation (Cockayne and Allan, 1926; Keim et al., 1989). Gathering all the evidences from molecular markers and morphology, there is incidence that a hybrid swarm had been established at the contact region of Q. kerrii and Q. austrocochinchinensis in XSBN Nature Reserve, especially at populations D3 and D4. The trees of Q. austrocochinchinensis in the adjacent regions (locations C, D1, D3–D6, and parts of D2) might have been already genetically swamped. The same situation might also have occurred in population E in Hainan, as SSR indicated that no Q. austrocochinchinensis purebreds exist, and the AFLP inferred co-existence of Q. austrocochinchinensis purebreds and its backcrosses.

It is worth noting that NewHybrids detected that the genotypes of few individuals contain different admixture levels of backcrosses to Q. kerrii and Q. kerrii purebreds in the two Q. kerrii purebred populations A and B (e.g., CK10, CK27, and CK30), which are distant from all the known populations of Q. austrocochinchinensis. Interspecific gene flow is a widespread and ongoing process among oaks, especially in species with close genetic relationship (Coart et al., 2002; Burgarella et al., 2009; Lepais et al., 2009; Moran et al., 2012). There is strong evidence that Neogene climatic changes had little impact on plant distribution in tropical and subtropical Asia (An, 2000; Su et al., 2013; Jacques et al., 2014), although evergreen oaks and other Fagaceae still experienced range shifts in this region (Xu et al., 2015; Jiang et al., 2016; Sun et al., 2016). A recent biogeographical study has indicated that the tropical zone could have extended further north in the geological past than it does today, e.g., the line 20°30′N was the northern biogeographical boundary of the tropical zone in south and southeastern China during the mid-Holocene (Zhu, 2013). The germplasm of Q. austrocochinchinensis recovered in Q. kerrii purebred populations probably reflects the historical gene flow between the two species and a once wider distribution of Q. austrocochinchinensis at this geological time in Indo-China. Future work will need to include more populations of both species, using both maternal markers and high throughput SNP markers to explore the genetic structure and couple the niche modeling to estimate the historical population size; such an approach could provide a better understanding on the genetic patterns of the two oak species.

Previous leaf anatomical work has suggested that population E is a Q. austrocochinchinensis purebred population with very distinct features compared to Q. kerrii, such as narrow leaves, sharp leaf apex, and margin tips (Song et al., 2015). However, the AFLP data showed evidences of hybridization in population E. Across all surveyed populations, only D2 was identified as a pure Q. austrocochinchinensis subpopulation by both molecular markers. It is situated in the core region of the Xi-Shan-Ban-Na Nature Reserve with thick woods and geographically isolated from the Q. kerrii population, which likely reduces the chances of gene flow and hybridization. Population D2 has a much higher forest canopy density than that of other populations, including population E.

According to our field observation, Q. austrocochinchinensis and Q. kerrii have different habitat preferences, as the former is likely to grow in closed and moist forests, but the latter tends to grow on open slopes or by roadsides (Figure S6). The fine-scale closed and shady habitats with high humidity might favor the growth of Q. austrocochinchinensis, but are also crucial to maintain its genetic distinction. Subpopulations D3, D4, and D5 are located in the buffer area of XSBN Nature Reserve, where human activities are intensified. These regions were inferred as the location of the hybrid swarm of the two species based on our study. Opening habitat as a result of deforestation will likely favor colonization by Q. kerrii and facilitate pollen invasion. Our study also showed that introgression of the two oaks occurred at disturbed habitats, and pure Q. austrocochinchinensis was found at core protection areas with minimum disturbance. Habitat selection is intimately tied to niche differentiation and coexistence in plant communities (Bazzaz, 1991). Although there are no studies demonstrating the habitat preference of these two species, our data suggests that fine-scale heterogeneous habitats may play an important role in shaping the genetic structure of the two species in areas where they co-exist.

The results of STRUCTURE and NewHybrids showed discrepancies between AFLP and SSR markers (Figure 3). For many samples, AFLP and SSR provided different classifications of species from the same set of samples. Such similar results have also been detected in other plant groups, e.g., Abies ziyuanensis (Tang et al., 2008) and the arctic-alpine genus Draba (Skrede et al., 2009). This is probably caused by the limited genetic differentiation between the two species, especially in a fine sampling scale. If the two species are not fully differentiated, they will have a number of incompletely differentiated alleles, which may reduce the diagnostic power to distinguish one species from another. In the cases of species with limited genetic differentiation, sampling large numbers of loci across the genome is required when using a molecular approach. The AFLP method has advantages over the SSR method in sampling loci numbers; consequently, AFLP generally reveals higher polymorphism than SSR does (Varshney et al., 2007; Sun et al., 2008; Skrede et al., 2009); it also has a higher assignment success and solution compared to SSR loci (Zeng et al., 2010). Our results also demonstrated that the AFLP data was more consistent with morphology in STRUCTURE and NewHybrids analysis, and the two species separated much better based on using AFLP in PCoA than by using SSR markers. Therefore, we speculate that the results generated from AFLP are more reliable than those obtained from SSR. Our results also demonstrated that when studying hybridization between two genetically closely related oaks, high throughput markers e.g., AFLP or RAD, are more appropriate.

In another aspect, the discrepancies between the two molecular markers might be due to different selective pressures. Similar to the results found by Scotti-Saintagne et al. (2004), the F_ST values of 781 AFLP alleles fit an L-shaped curve. Most alleles (60.3%) have low F_ST values (F_ST < 0.1). The results of BayeScan also suggested that only three loci (0.38%) are potentially under selection among the 10 populations. Therefore, most of the AFLP markers are selectively neutral. However, SSR markers may have been subjected to more selective pressure. Morgante et al. (2002) mentioned that SSR repeats tend to occur at transcribed regions of the genome. Selective pressures acting on coding regions are higher than in non-coding regions. In turn, such pressure may lead to genetic differentiation and as a result, SSR loci may have higher genetic differentiation levels than AFLP. We particularly used SSR loci that had been developed from coding sequences, e.g., Qk15874, Qk17139, Qk17611, Qk20944, and CR627959.

Although SSR provided less information than AFLP, it is still a valuable method to compare different marker types. For a pair of species undergoing hybridization, the impermeable genomic regions that are under high selective pressure often serve as a way to maintain species integrity, and thus, these regions accumulate genetic divergence, whereas the regions with low selective pressure are permeable to introgression and show decreased differentiation (Wu, 2001). Because of the higher selective pressures, SSR is more representative of regions with high species integrity, while AFLP is more likely to reflect the effects of gene flow. It is important to note that SSR identified more Q. kerrii than AFLP in the STRUCTURE analysis. Moreover, in the PCoA analysis for AFLP, the morphological intermediate populations C and D5 were clustered with Q. austrocochinchinensis purebreds. However, for SSR, they were closer to Q. kerrii purebreds instead. As discussed above, we suggest that although the trees in the swarm were influenced by the gene flow of Q. austrocochinchinensis, the alleles of Q. kerrii were preferred in these hybrids. If this is the case, it indicates the high genetic extinction risk of Q. austrocochinchinensis.

Our study revealed the discrepancies between the results inferred by morphological traits and molecular markers. For example, although morphological intermediates only make up a small proportion of samples in the hybrid zone (15 of 61), the molecular method detected more hybrids than expected. Moreover, individuals with one parental morphological characteristic may be identified as another parental purebred, e.g., a part of trees from population D3. There may be two explanations for these discrepancies. First, trees with clear parental characteristics could in fact be F₁ hybrids or backcrosses, and the morphological phenotypes may be a result of some of the loci having dominant effects on controlling a certain morphological trait. Moreover, hybrids that are in the second or later generation may return to homozygosis at specific loci and thus resemble one parent only (Stebbins, 1950; Rieseberg et al., 1993). However, this offers only a partial explanation, because of the low chance that all loci controlling morphological diagnostic characteristics return to homozygosis at the same time.

Another possibility of the inconsistency is that the phenotype with one of the parental phenotype may be under positive selection in natural habitats, as described previously (Hipp and Weber, 2008). As discussed above, the two species appear to have different habitat preferences. The two populations with typical Q. austrocochinchinensis in appearance were all located in the national nature reserve, where canopy density is usually high. The stellate and fused fascicular trichomes and dense stomata, as well as the thick leaf cuticles are important features in seasonal dry areas with strong sun light (Mediterranean-type) climate (Tattini et al., 2007). The glabrous leaves with sharp, spiny teeth were generally an adaptation trait to humid and shady environments (Sun et al., 2003), which may explain why the trees are distributed in dense canopy areas of the Xi-Shuang-Ban-Na Nature Reserve of Yunnan Province and Ba-Wang-Ling Nature Reserves likely show the phenotypes of typical Q. austrocochinchinensis, although introgression is common, as this phenotype might have a selective advantage over the hairy, obtuse toothed phenotype in ravine habitats where the environment is usually humid and shaded. Meanwhile, we also could not rule out the possibility that the morphological variation in Q. kerrii might encompass the typical morphology of Q. austrocochinchinensis. In this case, the morphological features used to distinguish Q. austrocochinchinensis from Q. kerrii need to be further clarified. Common garden experiment coupling the high throughput genotyping on the different populations of the two taxa and their hybrids are need to reveal how and why traits evolved in response to nature selection and local adaptation in the future.

Our data indicates ongoing and historical introgression in the two Q. austrocochinchinensis populations. Although the introgression directions inferred by SSR and AFLP are not the same (as SSR indicated unidirectional introgression from Q. kerrii to Q. austrocochinchinensis, but bidirectional introgression was inferred by AFLP). BayeScan analysis indicated most AFLP loci were selection neutral, but almost all the SSR loci violated HWE assumption. Therefore, AFLP result is more reliable to infer the gene flow between the two species. Recently simulation model study demonstrated when outercrossing rate between the two parental species is different, hybridization can facilitate invasions of the species with high outcrossing rate, even without enhancing local adaption (Mesgaran et al., 2016). Considering the fast shrinking and degradation of habits favoring Q. austrocochinchinensis at middle to low elevation of tropical Asia, the ongoing hybrid swarm and the predominance of Q. kerrii at co-occur region, Q. austrocochinchinensis faces critical extinction risks both in terms of lost habitat and genetic assimilation. Furthermore, our results indicate that only subpopulation D2 maintains its genetic integrity, as the most pure individuals of Q. austrocochinchinensis are still purebreds. The other subpopulations in XSBN Nature Reserve and adjacent areas formed a hybrid swarm. Once introgression is detected in an endangered species, protective intervention measures should be adopted immediately to prevent species integrity loss. Otherwise, all subpopulations in the reserve could become hybrid swarms.

Firstly, in situ conservation of species genetic resource can be realized by establishing nature reserves. According to our results, the most pure Q. austrocochinchinensis population D2 has the highest forest canopy density and lowest human disturbance. Although it is difficult to prove the correlation between ecological integrity and hybrid degree in this case study, we hypothesize that it is important for in situ conservation to maintain the existing ecological balance of the habitat of population D2. To avoid further asymmetrical introgressions, forest landscape restoration is also essential for the hybrid Q. austrocochinchinensis populations. Secondly, considering the abundance of Q. kerrii and the efficient pollen dispersal abilities of oaks, ex situ conservation of pure Q. austrocochinchinensis should also be considered. Although we cannot claim that Q. austrocochinchinensis homozygous individuals still exist, molecular markers and leaf morphological features all identified D2 as the most pure Q. austrocochinchinensis population. Therefore, we suggest that the D2 population area should be designated as the priority conservation zone.

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.