Introduction 

Hybridization between cork oak (Quercus suber L.) and holm oak (Q. ilex L.) has been known for a long time. The first written references to putative hybrid individuals date back to the mid-19th century ([5], [21], [3]). The ranges of cork and holm oak overlap in the western Mediterranean basin, where they form mixed stands on acid and decarbonated soils, the preferred ones for cork oak. Morphologically intermediate, presumably hybrid individuals, can be found in these stands, and are easily identifiable by their bark, and by some leaf intermediate traits, such as leaf thickness, which show characteristics of both parental species ([21]). These putative hybrid individuals occur at low frequency and are scattered among the pure cork and holm oaks. Introgression of genes from one species to another can have important adaptive and evolutionary consequences ([13], [30]). In the case of holm and cork oaks, it has been proposed that such process may have helped cork oak thrive in adverse environmental conditions (e.g., during the glacial pulses of the Pleistocene - [24]).

Molecular markers have confirmed the hybrid nature of phenotypically intermediate individuals and were used to assess the extent of hybridization between Q. suber and Q. ilex in natural mixed populations. Analysis of chloroplast DNA (cpDNA), which is maternally inherited in oaks, detected Q. suber populations in Southern France and Eastern Spain with a typical Q. ilex lineage cpDNA ([15], [23], [25]). Phenology of both species as well as post pollination barriers as those reported by Boavida et al. ([2]) favor hybridization with Q. suber acting as pollen donor, so that hybrids carry the ilex chloroplast. Chloroplast introgression pattern suggests hybrid trees carrying the ilex chloroplast recurrently backcrossed with Q. suber as pollen donor. Using isozymes ([11], [28]) or nuclear microsatellites ([34], [4]), previous studies confirmed the hybrid character of morphologically intermediate oak trees, and showed the current existence of hybridization. Burgarella et al. ([4]) estimated a low rate of ongoing introgression in cork oak over its range (~2%), a result confirmed by Lumaret & Jabbour-Zahab ([26] - ~4%). These studies suggested that introgression could occur bidirectionally. However, if the flowering phenology of first generation hybrids (F₁ hybrids) is more similar to Q. suber, as suggested by Perea García-Calvo ([29]), they would preferentially backcross with Q. suber, while mating with Q. ilex, especially when the latter species acts as pollen donor, would be much less frequent.

In the current work, we genotyped seedlings of morphologically putative hybrid trees, to identify backcrossing events with either Q. ilex or Q. suber as pollen donors, and to quantify the frequency of introgression towards each of the parental species. Our study provides a deeper insight in the hybridization and introgression processes which could be currently taking place between these two species.

  Material and methods 

Sampling site and acorn collection

The sampling site is located in Central Spain, (39° 59′ N, 05° 07′ W) in a mixed “dehesa”, i.e., highly anthropized open woodland. Quercus suber and Q. ilex are the dominant species, while there are some stands with comparatively few Q. faginea individuals. Four mature hybrid trees were identified by their phenotypes. These hybrid trees occurred in three distinct areas (IS, LG, and ZL) which differed in tree density and species composition. One hybrid was identified in IS, where Q. ilex and Q. suber occur at a density of 38 trees ha^-1 and where Q. suber is slightly more abundant. Two hybrids were found in LG, approximately 25 m apart and 740 m away from the IS hybrid. At LG, tree density was low (20 trees ha^-1) and dominated by Q. ilex. The fourth hybrid tree was found in ZL, where tree density was 123.5 trees ha^-1 , with more Q. suber trees (38.3 trees ha^-1 for Q. ilex and 85.2 trees ha^-1 for Q. suber), and approximately 2780 m away from the two LG hybrids and 3140 m away from the IS hybrid.

Acorns were collected in 2011, 2012 and 2014 on each of the four hybrids trees and also from four randomly selected cork oak and four holm oak trees located in IS (Tab. 1). Each year, all acorns from LG2 and ZL hybrid trees were collected, as well as approximately 1/3 of the acorns on the LG1 and IS hybrids.

All collected acorns were sown in peat-perlite substrate (3:1 proportion) in 3L containers in a greenhouse room kept at 22 °C, with a relative humidity of 50-60%. Seedlings were watered daily for 8 min. and five months after sowing, germination rate was estimated and leaf tissue collected from the surviving seedlings for DNA extraction.

Genotyping

DNA extraction was performed with a modification of the protocol by Doyle & Doyle ([8], [9]). We used eight polymorphic nuclear microsatellites (nSSRs) for genotyping: MSQ4, MSQ13 ([7]), QpZAG9, QpZAG15, QpZAG36, ([37]), QrZAG7, QrZAG11 and QrZAG20 ([18]). Microsatellites were amplified under standard PCR conditions ([34], [35]). Fluorescence labelled PCR products were analyzed in a 4300 LI-COR^® automated sequencer (LI-COR Biosciences, Lincoln, NE, USA). Allele sizes were determined with the SAGA Microsatellite Analysis Software (LI-COR Biosciences, Lincoln, NE, USA). The presence of null alleles was examined using ML-Null ([17]).

Nuclear admixture analyses

We carried out admixture analyses using the Bayesian clustering approach implemented in the program STRUCTURE ver. 2.3.4 ([32]). In this model, the posterior probability (q) provides an estimation of the proportion of an individual genome originating from each of K differentiated genetic pools. Values of K ranging from 1 to 5 were tested in order to find the optimal K value. A value of K=2 was identified as the optimal K value following the methodology described by Evanno et al. ([12]) and was used in all subsequent analyses. In order to improve q estimations ([4]), we included in the STRUCTURE analyses the genotypes of approximately 95 Q. ilex and 96 Q. suber parental trees from the same stand, previously analyzed by Soto et al. ([35]), along with 59.671 virtual individuals simulated with the software SimHyb ([36]). These simulated individuals included a majority of parental trees with some F₁ hybrids and backcrosses (Tab. S1 in Supplementary material).

Calculations were carried out under the admixture model assuming independent allele frequencies, as it has been done previously ([4]). A burn-in of 50.000 steps followed by 100.000 iterations were used, after verifying that results did not vary significantly across runs or with longer burn-in/iteration cycles. Assignation of seedlings to specific categories based on the q_s values (estimation of the proportion of the genome coming from Q. suber) used the following thresholds: q_s ≥ 0.85: suber; 0.85 < q_s ≤ 0.65: suberoid; 0.65 ≤ q_s < 0.35: hybrid; 0.35 ≤ q_s < 0.15: ilicioid; q_s ≤ 0.15: ilex. The software Colony ([16]) was used to infer selfing rates in the progeny using the adult trees as potential fathers and allowing up to one mismatch per seedling, and to detect potential crosses between hybrid trees.

The software R ([33]) was used to draw graphics, linear regressions between putative paternal q_s and q_so/q_sm (the ratio between the offspring and the maternal q_s), and to prepare a script to obtain all the possible paternal contributions for each offspring.

  Results and discussion 

Germination rates

A total of 671 acorns were collected from hybrid trees, and 53 and 49 acorns from Q. ilex and Q. suber, respectively (Tab. 1). Hybrid trees from LG2 and ZL produced fewer acorns than hybrid trees from LG1 and IS, and their germination rates were also lower (24% and 46%, respectively), while germination rates were above 70% for hybrids from LG1 and IS, as well as for the parental species trees (Fig. S1 in Supplementary material).

Assignation to specific categories

We did not detect null alleles for any locus. Autofecundation was ruled out, as the estimated selfing probability was <0.001 for every hybrid seedling. Moreover, mating between hybrid trees was also discarded after inspection of potential parents for each offspring estimated by Colony. The most probable number of differentiated genetic pools in the population was K = 2 (Fig. S2). The q_s values (estimation of the proportion of the genome coming from Q. suber) obtained for the hybrid trees were 0.8763 for IS, 0.5988 for LG1, 0.8680 for LG2 and 0.5489 for ZL. These results would classify the hybrid trees at LG1 and ZL as most likely first generation hybrids (F₁), and hybrid trees at IS and LG2 as more “advanced introgressed” individuals, at least second-generation backcrosses with Q. suber. These results are consistent with preferential backcrossing with Q. suber.

Then, q_s values were also estimated for the progenies of the hybrids (Fig. 1). The high q_s values obtained for most of the hybrid offspring, particularly for IS, LG2 and ZL, indicate a high contribution of Q. suber to their genomes and therefore previous works could have classified these seedlings as backcrosses to Q. suber. Conversely, the offspring of LG1 shows a broader range of q_s, with an average value of 0.4512, closer to Q. ilex. However, assignation to specific classes based exclusively on q_s may be biased using this set of markers, especially for advanced introgressed individuals ([36]). Actually, more than 30% of the virtual introgressed individuals included in the simulation as control would be misassigned according to q_s. Following Soto et al. ([36]), we used the ratio q_so/q_sm to determine the most likely genetic pool of the pollen donor for each offspring (Fig. 2). These data support backcrossing with both parental species with many of the LG1 offspring having Q. ilex as father, and many offspring of the IS and ZL hybrid trees having Q. suber as father. Offspring of the LG2 hybrid tree had either Q. suber or Q. ilex as father. Fig. 3 shows the distribution of the q_s ratio throughout the studied years for IS and LG1, the trees that produced more acorns. Although the general trends are kept for the whole period, some differences among years can be noticed, particularly for LG1: in 2011, probably a masting year with high acorn production, the vast majority of the seedlings came from Q. ilex pollen (84.5%), while in 2014, with a much lower acorn production, the proportion of effective pollinations from Q. ilex decreased to 54.1%.

Fig. 1 - Boxplots of the q_s parameter across families for (a) all years combined and separately for (b) 2011, (c) 2012 , and (d) 2014. Sample sizes are presented in Tab. 1.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Fig. 2 - Histograms of the relative frequencies of the q_so/q_sm ratio between the offspring and the maternal q_s for Is (a), LG1 (b), LG2 (c), and ZL (d) families. The bars to the right of the vertical red line indicate pollination by individuals with a higher percentage of suber genomes. The bars to the left of the vertical red line indicate pollination by individuals with a higher percentage of ilex genomes.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Fig. 3 - Histograms for the relative frequencies of the ratio between the offspring and the maternal q_s for the families Is in 2011 (a), 2012 (b), 2014 (c), and LG1 in 2011 (d), 2012 (e), and 2014 (f). The bars to the right of the vertical red line indicate pollination by individuals with a higher percentage of suber genomes. The bars to the left of the vertical red line indicate pollination by individuals with a higher percentage of ilex genomes.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

To get additional support for the results on pollination direction we prepared an R script to compare each individual offspring and its mother’s genotypes, and to obtain all the possible paternal genetic contributions (haplotypes). We then estimated their q_s value using STRUCTURE. Fig. 4a shows a boxplot of the distribution of the q_s values assigned to all the possible paternal haplotypes for each mother tree, and Fig. 4b and Fig. 4c show these distribution through the studied years for IS and LG1.

Fig. 4 - Boxplots of the q_s parameter of the inferred putative fathers across families (a), and across cohorts for Is (b) and LG1 (c) families.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Finally, Fig. 5 shows the significant correlations above 0.80 for the estimated paternal q_s and the ratio q_so/q_sm, especially for IS (Fig. 5a) and LG1 (Fig. 5b), the trees with larger progenies. This result confirms that using the ratio q_so/q_sm is a good approach to identify the genetic pool of the pollen donor, provided the mother tree is known.

Fig. 5 - Correlation of the ratio between the offspring and the maternal q_s and the q_s parameter of the inferred putative fathers for Is (a), LG1 (b), LG2 (c), and ZL (d) families. Fitted linear model for each family are indicated.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

These results support bidirectional introgression, where hybrid trees backcross with both parental species. Bidirectional introgression was also suggested by genotypic data collected by Burgarella et al. ([4]) using the same set of nuclear microsatellite markers, on adult trees phenotypically classified as Q. suber and as Q. ilex sampled all across the overlapping distribution ranges of these two species. Lumaret & Jabbour-Zahab ([26]), using the same set of microsatellite markers, also reported the possibility of bidirectional introgression.

Each adult hybrid receives pollen from both parental species but is preferentially pollinated by either Q. ilex or Q. suber (Fig. 2). The major pollen donor for each hybrid tree seems to reflect the abundance of parental tree species in the area. For example, the hybrid tree LG1, located in La Laguna, is mostly pollinated by Q. ilex, the dominant species in that area (the number of seedlings from LG2 is too low to be conclusive); on the other hand, IS and ZL hybrids produced more seedlings from Q. suber pollen, the more abundant species in these areas. Therefore, the frequency of backcrossing with a parental species could be at least partially driven by pollen availability. Simulation studies using the oak complex Q. robur × Q. petraea have shown that the direction of introgression can be influenced by the composition of the pollen pool on individuals or stigmas which is itself affected by the spatial configuration of the parental species ([20], [19]). The direction of introgression depending on the abundance of the parental species in an area would be consistent with the findings in Eucalyptus ([14]), Senecio ([31]) and some mangroves ([39]).

Based on differential growth rate of pollen tubes, Boavida et al. ([2]) reported significantly higher success rate for the interspecific crosses with Q. suber acting as pollen donor rather than as female parent. However, the high number of hybrid offspring having Q. ilex as pollen donor, supports at least the partial breakage of this prezygotic incompatibility between introgressed hybrids and Q. ilex pollen. Notwithstanding, a potential effect of post-zygotic barriers on the observed gene flow directionality cannot be discarded. If some type of incompatibility existed between the two parental species, a high proportion of viable embryos could be expected for families of parent trees that have similar specific categories (estimated in our case through q_s). Interestingly, the majority of seedlings of the IS parent, classified as an “advanced introgressed” cork oak tree (q_s = 0.8763), had Q. suber as pollen donor. However, although both ZL and LG1 hybrids show intermediate q_s values (0.5489 and 0.5988, respectively), they show contrasting proportions of effective pollinations from Q. suber and Q. ilex. As previously stated, the Q. ilex neighborhood of LG1 could account for the abundance of LG1 × Q. ilex backcrosses. In the case of ZL an overabundance of backcrosses to Q. suber was scored, considering the specific composition of the surrounding trees. Although the number of seedlings from this tree is still low to draw definitive conclusions, this fact could be due to a phenological synchrony between this tree’s female flowers and pollination of the surrounding Q. suber ([29]), or to post-pollination processes ([38]), such as genomic incompatibilities between this hybrid and surrounding Q. ilex individuals. As shown for European white oaks, first-generation hybrids will also more likely mate with the more abundant species, leading to potential asymmetric introgression ([22]). Asymmetric introgression based on the abundance of the parental species has been described in other Mediterranean ([27]) and non-Mediterranean tree species ([1], [6], [10]).

  Conclusions 

These results support current bidirectional introgression with Q. suber and Q. ilex. Analysis of the progenies of four open pollinated hybrid trees from a natural mixed population reveals that each mother tree is predominantly fertilized by one of the two parental species, usually the most abundant in its vicinity. The frequency of backcrossing could at least be partially driven by the availability of pollen when female flowers are receptive. However, other post-pollination processes, such as genetic incompatibilities cannot be excluded. To clarify this point, controlled pollinations and monitoring of pollen tube development and abortion of embryos should be performed. Additionally, further investigations should take advantage of modern high-throughput sequencing methodologies that will allow the examination of larger proportions of the species’ genomes. Moreover, using these methodologies, it may be possible to identify the genome portions that determine hybridizing capability, and/or reproductive isolation in the Q. ilex - Q. suber complex.

  Acknowledgements 

AS and ULH conceived the idea, designed the experiments, collected the plant material and drafted the manuscript; HS and ULH genotyped the seedlings and adult trees; ULH, HS and AS performed the admixture and statistical analyses. All the authors have contributed to the final manuscript.

This work has been funded by the project AGL2015-67495-C2-2-R (Spanish Ministry of Economy and Competitiveness).

The authors thank the two anonymous reviewers for their helpful comments and suggestions.

  References

  Supplementary Material

Fig. S1 - Germination rates of the progenies of the hybrids, Q. ilex and Q. suber for 2011, 2012 and 2014.
Fig S2 - Most probable number of genetic groups (K) inferred with the test of Evanno et al. (2005).
Fig S3 - Barplots showing the relative contribution of Q. ilex and Q. suber to each offspring as the proporpotion of q_i and q_s. estimated with STRUCTURE.
Tab. S1 - Number of virtual individuals by the specific categories defined according to q_s intervals generated with SIMHYB (Soto et al. 2018).