Plant material and DNA sequencing

We combined previously studied (Denk & Grimm, 2010) and new material to develop a sampling design of 221 individuals completely covering the taxonomic breadth and distribution range ofQuercus sect.Cerris in western Eurasia (File S1); some individuals with intermediary morphology were labelled as presumed hybrids. 158 individuals were studied for the first time, and new plastid data were generated for the previously and newly studied individuals. DNAs were extracted from silica gel-dried leaf samples with the DNeasy Plant minikit, following the manufacturer’s instructions. ThetrnH-psbA intergenic spacer was amplified and sequenced followingSimeone et al. (2016). The nuclear ribosomal 5S intergenic spacer (5S-IGS) was amplified with the primer pair 5S14a and 5S15 (Volkov et al., 2001;Denk & Grimm, 2010). Individual PCR fragments were ligated into a pGEM-T easy vector (Promega). The ligation mixtures, purified with the Illustra GFX PCR DNA Purification kit (GE Healthcare), were used to transformE. coli strain XL1-Blue electroporation-competent cells (recA1, endA1, gyrA96, thi-1, hsdR17, supE44, relA1, lac, [F’ proAB, lacIqZ ΔM15, Tn10 (tetr)]). The positive clones, selected on LB/Ampicillin plates, were identified by colony PCR using the amplification primers. Five to ten recombinant clones per individual were sequenced with the vector-specific universal primers (SP6/T7) at LGC Genomics (Augsburg, Germany). The GenBanktrnH-psbA sequences of several East Asian members of sect.Ilex (Q.baloot,Q.floribunda,Q.phylliraeoides,Q.semecarpifolia,Q.baroni,Q.dolicholepis,Q.spinosa), East Asian species of sect.Cerris (Q.acutissima, Q. variabilis) and further sequence accessions of the investigated group were included in the analyses. In addition, all 5S IGS sequences ofQuercus sect.Cerris available on GenBank were included in the final dataset, and the sequences ofQ.baloot andQ.floribunda were used as outgroups, based onDenk & Grimm (2010); all GenBank accession numbers are reported inFile S1.

Data analyses

Eye-checked electropherograms were aligned in MEGA7 (Kumar, Stecher & Tamura, 2016). Highly dissimilar clone sequences showing no BLAST match with the targeted regions (Altschul et al., 1990) were filtered. Final multiple alignments were obtained with ClustalW 1.81 (Thompson, Higgins & Gibson, 1994) and checked by eye. The diversity of the investigated regions was evaluated with MEGA7 and DnaSP5.1 (Librado & Rozas, 2009). Median-joining (MJ) haplotype networks for thetrnH-psbA region were inferred with Network 4.6.1.1 (http://www.fluxus-engineering.com/), treating gaps as fifth state. The MJ algorithm was invoked with default parameters (equal weight of transversion/transition), in order to handle large datasets and multistate characters.

After removal of identical clones, the total 5S-IGS sequences were used to build a Maximum likelihood (ML) tree with RAxML v8.2 (Stamatakis, 2014) using the in-built GTR + Γ model with the ‘extended majority-rule consensus’ criterion as bootstopping option (Pattengale et al., 2009). To infer inter-individual relationships, we applied the approach described byGöker & Grimm (2008) that allows transformation of data matrices of ‘associates’ (here: cloned sequences) into ‘hosts’ (here: individuals). The programg2cef (available athttp://www.goeker.org/mg/distance/) was used to transform the primary character matrix (‘associates’, total cloned sequences) into a character consensus matrix of the individuals (‘hosts’) using an association file defining the list of clone sequences belonging to the same individual. The uncorrected pairwise distances of the primary character matrix (‘associates’, total cloned sequences) was calculated and used as input to the programpbc (Göker & Grimm, 2008). This program allows transforming the primary inter-clone pairwise distance matrix into inter-individual distances matrices using different flavours, of which the ‘Phylogenetic Bray-Curtis’ (PBC) transformation performed best in the original study that compared data sets with similar properties than our data set. Here, we applied three of the distance transformations tested byGöker & Grimm (2008), in addition to PBC distances (option -b) also the minimum (MIN; -i) and average (AVG; -a) inter-individual clonal distances. AVG, MIN and PBC distance matrices were generated setting different minimum number of associates per host (-m option);m = 4 (the number of cloned sequences obtained in most individuals) was then used to infer a phylogenetic network using the Neighbour-Net (NN) algorithm (Bryant & Moulton, 2004) implemented inSplitsTree4 (Huson & Bryant, 2006).

PlastidtrnH-psbA diversity and biogeography

After removing the 34-bp inversion, thetrnH-psbA marker showed pairwise uncorrectedp-distances ranging between zero and 0.008 (Table 2). The highest intra-specific distance (0.006) was found inQ.suber; four species showed similar values (0.002;Q. cerris,Q. ithaburensis,Q.trojana,Q. libani), while the marker variation in the remaining taxa converged to zero.

The total matrix was 503-bp characters long, including several indels (1–8 bp) and six polymorphic sites resulting in twelve haplotypes (labelled H1–H12) with a medium overall diversity (h = 0.515). Haplotype H1 hit 100% sequence identity with three non-representative individuals assigned to East Asian species of sect.Ilex in Genbank (haplotype list, occurrence and gene bank matches shown asFiles S1 andS2). It was the most common haplotype, occurring in 68.6% of individuals and all taxa exceptQ.brantii,Q.look, andQ. ithaburensis subsp.ithaburensis (henceforthQ. ithaburensis). Haplotypes H2, H5–H7 and H11 showed 100% sequence identity with Mediterranean members of sect.Ilex (Simeone et al., 2016;Vitelli et al., 2017). H2 is the second most frequent haplotype, found in 10.6% ofQ. cerris,Q. trojana,Q.  ithaburensis subsp.macrolepis (henceforthQ. macrolepis) samples from Turkey, the Balkans and Italy; H5–H7 were found inQ.brantii,Q.cerris andQ.macrolepis from Turkey, Iran, and Israel. They were all shared withQ.coccifera andQ.ilex of the Aegean ‘Cerris-Ilex’ lineage. Haplotype H11 found in Iberian samples ofQ.suber was shared withQ.ilex of the ‘Euro-Med’ lineage. Rare haplotypes restricted to a single species were H7 (one accession ofQ.cerris), H8 (three accessions ofQ.libani; new ‘Cerris-Ilex’ subtypes) and H11/H12 (eight accessions ofQ.suber; ‘Euro-Med’ types); all other haplotypes were shared by more than one species of sect.Cerris.Quercus cerris, the most widespread and ecologically diverse species of sect.Cerris, showed the highest number of haplotypes (eight), followed byQ. brantii (five) andQ. macrolepis (four). In the latter species, haplotype H6, exclusively found inQ.brantii, was also found in a suspected hybridQ.macrolepis ×Q. brantii (sample ml27). All samples ofQ. ithaburensis exhibited a single haplotype (H9). The geographically (more) restricted taxaQ. afares,Q.castaneifolia, Q.crenata andQ.trojana subsp.euboica (henceforth:Q.euboica) showed only the most frequent and widespread haplotype (H1).

In comparison (Table 2), the two East Asian members of sect.Cerris (Q.acutissima andQ.variabilis) displayed a higher variation at thetrnH-psbA locus, although mostly due to indels. A higher number of haplotypes was found in these species; none of them was shared with any species of sect.Ilex available in gene banks (identity range: 93–99% withQ.baroni,Q.dolicholepis and/orQ.spinosa), and only one haplotype was shared between the two species.

No shared parsimony informative characters (PICs) were found in the East Asian samples. In contrast, three PICs were exclusively shared by haplotypes H11 and H12 (Q. suber from Iberian Peninsula, North Morocco). One further PIC separated H9, including all individuals ofQ.ithaburensis, twoQ. look, two Israelian and one ItalianQ. cerris individuals. A single PIC also defined H4, including three co-occurringQ.trojana andQ.macrolepis accessions from the same locality in western Turkey, and another PIC was limited to twoQ.cerris andQ.libani accessions from southern Turkey, corresponding to H3.Table 3 shows that the highest mean intragroup divergence in the West Eurasian dataset was found inQ.suber andQ.libani. The haplotypes ofQ. suber andQ. look displayed the highest mean divergence from all other species.Quercus variabilis appeared more similar to its western Eurasian counterparts, whileQ.acutissima was highly distinct.

The haplotype network (Fig. 1) shows the general coherence of the ‘Cerris-Ilex’ lineage, which collects haplotypes typical of the western Eurasian members of sect.Cerris, clearly distinct from the haplotypes found in East Asian members of sect.Cerris and haplotypes H11–H12 (>5 mutations separating each lineage); this latter represents a unique, early diverged plastid lineage, most frequent in the western Mediterranean populations of sect.Ilex (Simeone et al., 2016;Vitelli et al., 2017). Based on the relative number of mutations (1–5) separating each haplotype, the ‘Cerris-Ilex’ lineage can be further subdivided into two groups: (L1) a group of potentially primitive (non-derived) haplotypes (H1–H4) including the most common haplotypes (H1, H2) and still close to haplotypes found in the north-easternmost species of sect.Ilex (Q. phylliraeoides), the plastid sister lineage of ‘Cerris-Ilex’ (Simeone et al., 2016); (L2) a group of derived haplotypes (H5–H10). Haplotypes not shared with sect.Ilex (H3–H4, H8–H10, H12) are derivates of the ‘Cerris-Ilex’ and ‘Euro-Med’ main types (H1–H2, H5-H7, H11), shared by both sections in the Aegean and the western Mediterranean regions. As shown inFigs. 2A–2B, plastid diversity is largely decoupled from species identity and related to geography; the least derived haplotypes within the ‘Cerris-Ilex’ lineage (H1 and H2) occur across the whole distribution range of the investigated group, except the Levant and the western Mediterranean (H2). All other haplotypes are more circumscribed, concentrated in Anatolia (H3–H5, H7–H8), the Levant (H6, H9–H10; the latter two showing single occurrences in Italy), and Iberian Peninsula + Morocco (H11-H12).

Nuclear 5S rDNA diversity and species phylogeny

In contrast totrnH-psbA, 5S-IGS sequence variation appeared generally correlated with the taxonomy of the studied individuals, and allowed inferences on potential reticulation and inter-species relationships within the western Eurasian members of sect.Cerris. The 5S-IGS clones varied greatly in sequence features and length (the multiple alignment of the cloned sequences can be viewed in the Online Supplementary Archive at the journal’s homepage). For instance, allQ. brantii clones displayed an intra-specific (ATTT)₁₋₇ simple sequence repeat (SSR) variation. In all the other species, this motif was either absent (replaced by a 5–12 bp long poly-T) or consisting of one to two repeats, with the exception of two clones of the suspected hybridQ. macrolepis ×Q. brantii (individual ml27) that showed four to five repetitions. Two clones ofQ. brantii (sample br02; C. Turkey) shared a four-bp insertion with several clones of sympatricQ. macrolepis individuals (ml20–ml22). ThreeQ. libani individuals (li02, 03, 04; S. and E. Turkey) displayed a long indel (ca. 100 bp) in (nearly) all clones (‘shortlibani variant’ cf.Denk & Grimm, 2010). The extended sample revealed that the ‘shortlibani variant’ is not exclusive toQ. libani but is rarely found also inQ. cerris (clones ce2104 and ce4704; Italy, W. Turkey) andQ. trojana (three clones of individual tj33, S. Turkey); the latter, however, is another suspected hybrid (Q. trojana xQ. libani).

Two other deletions were detected in the same region of the ‘shortlibani variant’. One (22 bp) was shared by single clones of twoQ.cerris individuals (ce18, ce22; S.W. and W. Turkey), four clones of aQ. look individual (lk2; Israel) and two clones ofQ. trojana (individual tj40; S. Turkey). The second (∼100 bp), largely overlapping with the deletion of the ‘short libani variant’, but beginning a few basepairs downstream, was shared by one clone ofQ. cerris (ce34; N. Turkey) and one clone ofQ. macrolepis (ml26; S. Greece). An 8-bp deletion occurred exclusively inQ.suber andQ.crenata, with the exception of single clones of samples su07, su09 (N.E. and S. Spain), su37 (Croatia), su53 (S. Italy), cr02 (C. Italy), two clones of sample cr04 (Slovenia) and cr06 (N.E. Italy), and three clones of sample cr05 (Croatia). The same deletion also occurred in two clones of sample tj08, aQ.suber ×Q.trojana cultivation hybrid. Further deletions (1–60 bp) were scattered along the alignment and found only in single individuals (e.g., it04, Israel; ml10, N.W. Greece). Finally, an 18-bp highly variable region was exclusively found in some clones of four co-occurringQ.trojana samples (tj03–05, tj16; S.C. Turkey).

The main diversity values of the investigated dataset are reported inTable 4. Identical 5S-IGS sequences typically occur in the same individual and species, and, to a lesser extent, in sympatric, different species (e.g.,Q.brantii,Q.cerris,Q.trojana,Q.look; see alsoFile S3). On the contrary,Q.afares,Q.castaneifolia,Q.libani, Q. ithaburensis, Q. macrolepis andQ. euboica showed high species coherence.Quercus suber andQ. macrolepis showed the highest number of intra-individual and intra-specifically shared clones, whereasQ. cerris, Q.  trojana andQ. ithaburensis displayed the highest levels of unique variants. No variants were shared betweenQ.trojana andQ. euboica;Q. suber andQ.crenata (but notQ. cerris) shared 69 identical sequences and are the genetically most similar taxon pair. The pairwise uncorrectedp-distance range of the total dataset was much higher than for the plastid marker (0–0.209), with highest values scored byQ.cerris andQ. trojana. The mean intra-specific molecular diversity estimated within sequence pairs (Table 5) was lowest in the two narrow endemicsQ. afares andQ. castaneifolia and highest inQ. brantii andQ. ithaburensis. Across the entire dataset,Q. brantii, Q. macrolepis andQ. ithaburensis were the most diverging taxa; the least divergent beingQ. afares andQ. castaneifolia. The mean divergence value betweenQ. macrolepis andQ.ithaburensis (0,0376), treated as subspecies ofQ. ithaburensis in current regional floras, was similar to values detected between these taxa and the other species (e.g.,Q. afares,Q. brantii,Q. libani). Likewise, the divergence recorded between the putative conspecificQ.trojana andQ. euboica (0.0266) was comparable to the estimates calculated between these and other taxa (e.g.,Q.afares,Q.cerris,Q.look,Q.libani). The putative hybrid taxonQ.crenata displayed the lowest divergence (0.0197) withQ.suber, one of the assumed parental species, and a slightly higher estimate (but similar to the values scored with other taxa, e.g.,Q.afares,Q.castaneifolia, Q. look) withQ.cerris (0.0266), the other putative parental species.

The clone-based ML tree rooted onQ.baloot andQ.floribunda (West-Asian members of sect.Ilex) showed four main topological features (grades/clades) generally coherent with taxonomy (Fig. 3, see alsoFile S4). These grades/clades collected to a large degree clones of (1)Q.crenata andQ. suber (resolved as proximal, weakly differentiated grade), (2)Q. brantii, Q. ithaburensis andQ. macrolepis (the most highly supported clade: BS_ML = 84), (3)Q.trojana (a large heterogeneous grade), and (4)Q.cerris (the distal, terminal, clade with diminishing support).Quercus libani clones (short and normal-length variants) were present in all clades/grades except grade 1. A moderately supported clade (BS = 63) including allQ.afares clones was placed as sister to the main clade including clades/grades 2–4;Q.castaneifolia clones were placed within grade 3. Clones ofQ. ithaburensis also occurred in grade 3,Q.brantii andQ. crenata in clade 4,Q.look andQ.euboica in grade 3 and 4. A few clones ofQ.cerris andQ.trojana occurred scattered across the tree (often in proximal positions).

Of the three clones sequenced from individual ml27, a suspectedQ. macrolepis ×brantii hybrid, one was identical to anotherQ. macrolepis clone (individual ml08) and the other two clustered together withQ.brantii. Likewise, three and two of the five clones sequenced in sample tj08, aQ.trojana ×Q.suber hybrid, clustered within the respective parental subtrees; the same applies to the five clones of the sample tj33, a supposedQ. libani ×Q.trojana hybrid. Conversely, all the three clones sequenced in sample tj02, another tree determined as possibleQ.libani ×Q.trojana hybrid, clustered withQ.  trojana.

The networks based on transformed 5S-IGS data (Fig. 4 based on AVG-transformed uncorrected distances;Fig. 5 based on PBC-transformed distance matrix; only individuals represented by more than four clones included) largely confirmed the earlier found intra- and inter-species relationships (Denk & Grimm, 2010; because of the amount of shared identical clones, the MIN-transformed networks are largely collapsed, but included in the Online Supporting Archive). Four clusters can be observed: Cluster 1, the ‘oriental’ lineage of sect.Cerris, is the least coherent cluster and equivalent to a grade in a corresponding outgroup-rooted (Q. baloot, Q. floribunda) tree. This lineage included, in the AVG network,Q.afares, four out of fiveQ.libani individuals and about half of theQ.trojana individuals (Fig. 4). Its counterpart, Cluster 2, the ‘occidental’ lineage, accomodated allQ.look, Q. euboica, the remainingQ.trojana andQ.libani samples, and all but oneQ.cerris individual (ce50;Figs. 4 and5). The PBC network (Fig. 5) reveals a more gradual shift between these two clusters, withQ. afares splitting off with two genetically similarQ. cerris andQ. libani individuals (ce50, li05). The reason for this is that the PBC transformation has a higher chance to capture evolutionary signals (Göker & Grimm, 2008). Cluster 3 includedQ.suber andQ. crenata. Here, the only difference is the boxyness inflicted by individuals cr05 and tj08. Cluster 4 included the ‘Vallonea’ (orAegilops) oaks,Q. brantii, Q.ithaburensis andQ. macrolepis, with twoQ. brantii individuals (br02 and br03, with diverging 5S-IGS features and variants;File S4) in proximal (br02 inFigs. 4 and5; br03 inFig. 4) or off-cluster (br03 inFig. 5) position. Thus, the basic structure of the AVG and PBC networks and the ML tree are equivalent, but they differ in placing the outgroup taxa, and the networks refine inter-species relationships.

The AVG network (Fig. 4) better captured putative hybrid and intrograded individuals. Strong ambiguous signals came from the hybridQ.trojana ×Q.suber (tj08), oneQ.crenata individual (cr05, terminals in the box-like structure connecting the ‘occidental’, cluster 2, andcrenata-suber lineage, cluster 3), and oneQ.ithaburensis individual (it03, terminal in the box-like structure between clusters 2 and 4). The placement of oneQ. libani individual (li01, inserted in cluster 2), with normal-long variants in the clone sample, and oneQ.cerris (ce50, in cluster 1 close to theQ. afares subgroup; cf.Fig. 5) does not follow the general trend. Long terminal edges indicative of unique individual clone samples (combinations) are found in each cluster. Besides the outgroupQ.baloot andQ.floribunda, these samples include individuals ofQ.cerris (ce29, 44),Q.trojana (tj03, 16, 24, 39, 45),Q.suber (su07, 29, 49),Q.brantii (br06),Q.ithaburensis (it04, 05), andQ. macrolepis (ml10). Some of these samples had unique deletions or highly divergent regions in their clones (e.g., it04, ml10, tj03, tj16; see above). The networks produced with individuals represented by ≥2, ≥3 (and ≥5) clones did not change this structure (File S5); they allowed inclusion of all individuals into the four clusters matching the general scheme and pinpointed a few other (possible) exceptions.Quercus castaneifolia (represented by two clones) and one sample ofQ.crenata (cr04; three clones) formed part of the ‘occidental’ lineage, cluster 2; one sample ofQ.suber (su09; three clones) was included in cluster 1, the ‘oriental’ lineage of Sect.Cerris (seeTable 4); sample br01 (three clones) was placed at the root of cluster 4, similarly to samples br02 and br03. The geographical distribution of the four clusters is shown inFig. 2C.

In contrast, the PBC network (Fig. 5) provided a better basis for inferring the evolution and differentiation (speciation) processes. The ‘oriental’ (cluster 1) and ‘occidental’ lineages are clearly connected and form a continuum, with the easternboundQ. trojana andQ. libani representing a diverged, differentiated pool from which the other species and westernQ. trojana derived. The western Mediterraneancrenata-suber lineage is clearly different and only linked to the main pool by occasional introgression or hybridisation with nearby members of the ‘occidental’ lineage (in nature:Q. cerris). The same holds even more for the ‘Vallonea’ oaks (Cluster 4), which appear to have split before the remainder of western EurasianCerris (but long-branch/-edge attraction with the extreme long-edged outgroup needs to be considered). A clear signal in the PBC network (Fig. 4) is the uniqueness ofQ. afares, a disjunct outpost of the putative ‘oriental’ lineage, genetically closely related to geographically very disjunct (C./S. Anatolian) individuals ofQ. cerris andQ. libani (cf.Fig. 3 showing aQ. afares subclade, andFile S4, same tree with clones labelled).

Molecular recognition of species and species diversity inQuercus sectionCerris

Widespread species such asQ. cerris and (to a lesser extent)Q.libani (Table 1) showed the highest plastid diversity in terms of number of haplotypes and parameters of molecular differentiation (Tables 2 and3);Q. suber diversity was inflated by the occurrence of few divergent haplotypes linked to—and possibly captured from—the ‘Euro-Med’ lineage of sect.Ilex.

The strikingly high haplotype richness ofQ. cerris, especially in the eastern part of this species’ range (cf.Bagnoli et al., 2016), mirrors the high morphological plasticity of this oak (many different varieties and subspecies are traditionally reported; for example, IOPI lists 30formae and 17 varieties) and its ecological adaptability. Likewise,Q. cerris also displayed a high nuclear (5S-IGS) diversity (Table 4,Fig. 3). This oak has the largest range and broadest climatic envelope (from perhumidCfa, Cfb via summer-dry warm temperate climates toBSk) and it is the only species of sect.Cerris naturalized on the British Islands (Cfb) and cultivated all over continental Europe (mostlyCfb, shelteredDfb). Indeed, establishment of a large range across the geologically and ecologically dynamic West Eurasian region might have provided many opportunities for diversification, isolation, drift, conservation of variants and eventual reticulation with sibling species.

Quercus suber, instead, displayed the lowest number of unique 5S-IGS variants (likeQ. macrolepis) with a low diversity (Tables 4 and5), which might indicate ongoing genetic erosion (possibly due to the species domestication for cork, tannins, wood and fruits exploitation). In view of the high species-coherence detected in the other conspecific samples, the few Iberian samples not following the general trend may indicate introgression ofQ.suber intoQ.ilex (individuals with ‘Euro-Med’ plastid haplotype H11) or reflect ancient reticulation and retention of ancestral signatures (individuals with the widespread, putatively ancestral H1 haplotype).

InQ. libani, the high haplotype diversity coincides with a moderate diversity at the 5S-IGS locus, characterized by interesting variation among the cloned sequences. Two individuals (li01, li05) show potential introgression (e.g., with sympatricQ.trojana samples tj05, tj35;Files S1 andS4). Alternatively, they might represent an ancestral line of diversification within the section (the short 5S-IGS variant so far only known forQ. libani, seeDenk & Grimm (2010), also occurs in clones of two non-sympatricQ.cerris individuals). Together withQ. brantii, Q. libani is the easternmost oak among the western EurasianCerris, and it is extremely variable at the morphological and ecological level (occurring in climates ranging fromCsa toDsb). For instance,Djavanchir-Khoie (1967) described up to 12 intra-specific taxa withinQ.libani, and introgression phenomena with other co-occurring oaks have been postulated (Menitsky, 2005;Khadivi-Khub et al., 2015). Likewise, in the nearby region of Iranian Kurdistan, this oak showed three distinct gene pools based on nuclear microsatellites (Khadivi-Khub et al., 2015).

Based on the phylogenetic reconstructions, the 5S-IGS sequence diversity detected inQ. trojana outmatched the high levels displayed byQ.cerris, rendering it the genetically least-coherent species of the section. Samples with highly diverse clones were detected in both species (see ‘Results’), but individuals ofQ.trojana can be found in two different clusters in the 5S-IGS network. Reflecting its more limited distribution and climatic niche (Csa, Csb), the haplotype diversity is substantially lower inQ. trojana s.l. (Q. trojana + Q. euboica) than inQ. cerris (Tables 2 and3). This finding indicates that the two species retain exceptionally high intra- and inter-individual variability, possibly conserving ancestral variants lost in more homogenized and/or geographically restricted species of sect.Cerris, and corresponds, at the plastid level, to the relative extension of their geographic (and ecological) ranges. Accordingly, the two endemic taxaQ.afares andQ. castaneifolia appeared the least diverse (Tables 4 and5), although more data are needed forQ.castaneifolia, which was here represented by only two samples. The same holds true for the two other taxa with narrow ranges in our dataset,Q.euboica andQ.look. Both are characterized by low levels of genetic diversity (Tables 4 and5); however, despite the low number of individuals investigated, our results allow first taxonomic inferences in both cases.

Quercus euboica appears genetically isolated fromQ.trojana, based on the number of unique 5S variants (Table 4) and the relative inter-taxa divergence (Table 5). This oak grows isolated fromQ.trojana on the Greek island of Euboea and differs morphologically by its coriaceous leaf texture and the conspicuous white tomentum of the abaxial leaf surface that is made up of stellate trichomes (T Denk, pers. obs., 2008 and 2017). In addition,Q. euboica is characterized by special edaphic conditions, growing on serpentine rocks. All these data indicate that the Euboean oak should be better considered as an independent species requiring special protection. Another (hairy) variant ofQ.trojana has been locally described at the south-eastern margin of the species’ range, in South-central Turkey (Q.trojana subsp.yaltirikii;Zielinski, Petrova & Tomaszewski, 2006). Some samples collected in the nearby area (tj03, 04, 05, 16) showed 5S-IGS clones with a unique, highly divergent motif, and grouped (mostly) in a specific sublcade (Fig. 3,File S4). However, more (morpho-ecological) data are needed to implement the description ofQ. trojana in this part of its range.

TheQ.look samples showed a distinct plastid-nuclear signature combination (Tables 2–4;File S4) linking it withQ. cerris (Figs. 3–5). In addition, this species showed the lowest mean estimate of evolutionary divergence of the nuclear 5S IGS together withQ. castaneifolia (Table 5). Although the taxonomic rank of this rare, enigmatic taxon cannot be yet established with certainty (a hybrid origin or a local diversification of an ancestral form ofQ.cerris seem equally probable), the two previous assessments of this oak as synonym ofQ.ithaburensis orQ.ithaburensis × Q. libani hybrid (Table 1) can be rejected. Additional investigations are required to evaluate if the morphology ofQ. look justifies its exclusion from the genetically and morphologically variableQ.cerris.

Finally, a distinct group with medium-high levels of plastid and nuclear diversity includesQ. brantii, Q. ithaburensis, andQ. macrolepis (Tables 2,4 and5; ‘Aegilops oaks’) with a range centred in theCsa climates of the central-eastern Mediterranean region, providing a low-land analogue to the situation inQ. trojana-euboica-libani. The Aegilops oaks are a highly specialized group, morphologically and ecologically well distinct from the other oaks in theCerris section (Menitsky, 2005). The detected genetic diversity at the plastid and, especially, the nuclear markers (Table 5) clearly indicates the genetic isolation from the rest of theCerris oaks and progressive inter-specific differentiation. Geographic, morphological and ecological differences are also evident inQ. ithaburensis andQ. macrolepis (Dufour-Dror & Ertas, 2002;Dufour-Dror & Ertas, 2004). On these grounds and considering the high inter-taxon 5S-IGS divergence supported by the different haplotypes (Tables 2–4), we suggest these two forms be treated as separate species (cf.Denk et al., 2017, appendix:http://dx.doi.org/10.1101/168146). Interestingly,Q.brantii appeared as the most diverse of the three taxa (Tables 4 and5), and displayed some 5S-IGS variants shared withQ.cerris andQ.suber (File S3), which might indicate the occurrence of ancestral traits.

Hybrid detection withinQuercus SectionCerris

Besides the deletion inQ. libani, some other sequence features, typical of other species (e.g., the SSR motif inQ.brantii, the deletion inQ.suber), confirmed the hybrid identity of a few samples included in our dataset (sample ml27, tj08 and tj33, supposed hybridsQ. macrolepis ×brantii,Q.trojana ×Q. suber andQ.trojana ×Q. libani, respectively). The haplotypes of these samples (ml27: H6, exclusive ofQ. brantii; tj08: H2, never found inQ.suber), also allowed identification of the maternal species. This finding confirms that these oaks can occasionally hybridize in sympatry, and evidence for such hybridization is found in the nuclear genome (see alsoFitzek et al., 2018).

Further instances of hybridization and/or introgression events could be inferred from common sequence features, inter-specifically shared variants (see ‘Results’ section,Table 4,Fig. 3,File S3), and the placement of individuals in the AVG Neighbour-Net (Fig. 4), mostly involving Anatolian samples ofQ.cerris,Q.brantii,Q.libani,Q.macrolepis andQ.trojana. However, in many cases the involved individuals grew hundreds to a few thousands of kilometres from each other. Based on their (relative) spatial proximity, introgressions could be suggested for samplesQ.brantii (br02; C. Anatolia) andQ.macrolepis (e.g., ml20–22; W. and S. Anatolia) sharing sequence features, and South AnatolianQ.libani (li01) andQ.trojana (e.g., tj05, tj35) sharing sequence features andtrnH-psbA haplotypes. Outside Anatolia, evidence of reticulation (shared variants) can be traced between IsraelianQ.cerris andQ.look (sample ce38, lk03) and between BalkanQ.cerris andQ.suber (sample ce43, Serbia; su37, Croatia). Introgression and past hybridization events between all these species or their precursors is a possible explanation. At the same time, retention of ancestral traits cannot be discarded, asQ.cerris (5S-IGS,trnH-psbA) andQ.trojana (5S-IGS) cover most variability found in sect.Cerris and are highly variable, especially in Anatolia.

In this context, the unresolved taxonomic status ofQ.crenata can be discussed in view of the present results (Figs. 3–5). The species is more closely related toQ. suber than toQ. cerris and part of the distinctQ. crenata-suber lineage. The relative distribution of the clones (in theQ. cerris- and in theQ. suber-dominated clades) and the large extent of variants shared withQ. suber can be interpreted as evidence of co-existing both (1)Q.cerris ×Q.suber F1 hybrids and (2) introgressive forms into eitherQ.cerris (North East Italy/Balkans) orQ.suber (Italian peninsula), in partial agreement withConte, Cotti & Cristofolini (2007). However, all forms would look phenotypically quite similar (intermediacy of habitus, leaf and bark shape betweenQ.suber andQ.cerris is traditionally used as a diagnostic character ofQ.crenata), which seems inconsistent with their presumable different genome composition. Clearly, the hybrid/introgressed phenotypes can be affected by several phenomena such as segregation, epistasis, heterosis, and maternal origin (Rieseberg & Ellstrand, 1993). At the same time, we note that the diagnostic traits used forQ.crenata occur in other species of sect.Cerris (corky bark:Q. afares,Q. variabilis, various forms ofQ. cerris;Menitsky, 2005; semi-evergreen habitus:Q.trojana andQ. libani;Yaltırık, 1984; crenate leaves: part of the morphological variation ofQ.cerris andQ.suber). Also,Q. cerris shares plastid and nuclear signatures with other species of sect.Cerris, including geographically isolated, morphologically distinct taxa such asQ. afares, Q. castaneifolia,Q. euboica (traditionally included inQ. trojana), andQ. look (traditionally included inQ. ithaburensis), which are part of the ‘occidental’ lineage (see ‘Results’). Besides occasional hybridizations (sample cr05;Fig. 4), the alternative explanation is thatQ.crenata represents a less-derived species possessing a limited gene pool within an autonomous evolutionary lineage includingQ.suber (Fig. 5). Being closer to the common root, it retained imprints of common origin, possibly ancient reticulation, with the (proto-)Q.cerris (‘occidental’) lineage, representing a geographic-evolutionary gradient (‘oriental’ lineage → ‘occidental’ lineage →crenata-suber lineage).Quercus crenata may then just represent the remainder of the ancestral form from whichQ. suber evolved rather than being the product of secondary contact betweenQ. cerris andQ. suber.

We also found no evidence to support the hybrid origin ofQ.afares (Q.suber xQ. canariensis; the latter is a member of sect.Quercus) as suggested byMir et al. (2006) based on cpDNA-RFLP and allozymes (cf.Welter et al., 2012;Mhamdi et al., 2013). 5S-IGS variants and plastid signatures of western Eurasian white oaks (‘roburoid oaks’ inDenk & Grimm, 2010; see alsoSimeone et al. 2016,Fig. 1) are very distinct fromCerris types and should be detectable unless the F1 hybrids, withQ. suber as a maternal parent, only backcrossed with the localQ. suber but notQ. canariensis. However, genetic exchange with localQ.suber can be excluded, since noQ.suber-typical 5S-IGS variants, or obvious phylogenetic linkage with theQ. crenata-suber lineage, were found inQ.afares. Ongoing next-generation target sequencing of the 5S-IGS region (producing several 10,000 5S-IGS sequences per sample/individual) showed, so far, no evidence for a clone-sampling artefact in the studied individuals ofQ. afares. The possibility of incomprehensive clone-sampling can thus be discarded. Analoguous toQ. crenata, the ancestry level ofQ. afares may explain earlier findings interpreted towards a hybrid origin: being much closer to the common ancestor of sect.Cerris thanQ. suber, this species may have retained (some) genetic imprints today found in members outside its section. This would explain also the association ofQ. afares with two other, geographically very distant, individuals of the ‘oriental’ lineage (ce50, li05, the onlyQ. libani without the ‘short’libani 5S-IGS variants but showing variants similar to AnatolianQ. trojana).

Aside from putative hybrid species and swarms (see also the ambiguous placement of sample it03 the 5S-IGS network), our data demonstrates a general permeability of species boundaries in members of Sect.Cerris, allowing occasional crosses. Indeed, our results demonstrate thatQ.trojana andQ.libani,Q. brantii andQ.macrolepis,Q. cerris andQ.suber, andQ. trojana andQ.suber are interfertile and hybridize in the wild. Further investigations (e.g., with metagenomics approaches, fine-scale geographic samplings) are needed to distinguish between ancient hybridization with subsequent incomplete lineage sorting and retention of ancestral traits, to clarify the status of several other samples that may represent both phenomena based on shared 5S-IGS variants, the occurrence of unique sequence features, length of terminal branches and odd-placing in the phylogenetic reconstructions. Adequately addressing these issues would be of great relevance to identify relict populations and/or past contact/hybrid zones, to assess the hybridization ability of species growing in sympatry, and further define the evolutionary history of theCerris oaks.

Taxonomic framework ofQuercus sectionCerris in western Eurasia

From the clusters identified by the 5S-IGS network based on the average inter-individual clone (Fig. 4) and PBC-transformed distances (Fig. 5), four major groups can be identified representing distinct evolutionary lineages and used as a framework.Figure 6 shows a scheme, a cactus-type branching silhouette (Podani, 2017;Morrison, 2018), based on the 5S-IGS andtrnH-psbA differentiation patterns, and with respect to the plastid tree provided inSimeone et al. (2016).

The first, most western lineage includesQ.suber andQ.crenata. Sample cr05 likely represents a F1 hybrid withQ.cerris. The same holds for tj08, aQ.trojana xQ.suber hybrid. A second lineage, tentatively termed the ‘occidentalis’ lineage includes the widespreadQ.cerris and the geographically restrictedQ.castaneifolia, Q. euboica andQ.look. The third lineage collects the Eastern Mediterranean ‘Vallonea’ oaksQ.brantii,Q. ithaburensis, andQ. macrolepis, with three potential outliers: br02, br03 and it03 (possible recent or ancient hybrids). The last, least coherent group, the ‘oriental’ lineage, includesQ.afares,Q.libani andQ. trojana. Aside from supposed hybrids and introgressed individuals (see section above), these groups almost perfectly match previous taxonomic observations (sect.Heterobalanus subsect.Suber + sect.Cerris subsect.Cerris andAegilops;Menitsky, 2005; sect.Suber + sect.Eucerris,Aegilops andErythrobalanus;Schwarz, 1936–1939), with the only exceptions ofQ.afares, accomodated by both authors within theQ.cerris-Q.castaneifolia group, andQ.crenata andQ.look that were not included in previous monographs.

Quercus trojana is the only species scattered over two clusters (cluster 1 and 2,Figs. 4 and5; see alsoDenk & Grimm, 2010), thus bridging between the ‘oriental’ and ‘occidental’ lineage. No geographic, haplotypic or subspecific relationships could explain this subdivision; it might therefore indicate the occurrence of two different, geographically overlapping but genetically isolated lineages within the species, possibly differentiated in the past (retained ancient polymorphism), especially considering the proximal positions of the sequences of samples tj24, 39 and 45 in the ML tree. BothQ.trojana (s.str.) sublineages occur in Italy and might correspond to the two main nuclear gene pools identified byCarabeo et al. (2017). Indeed, more ecological and molecular data are required to interpret this finding biologically.

Overall, the complex genetic differentiation patterns can only be explained by longer ongoing free genetic exchange or more recent common origin of the ‘oriental’ and ‘occidental’ lineages than in case of the other two lineages within western EurasianCerris: the corkish oaks (Q. crenata + Q. suber) and the ‘Vallonea’ (or Aegilops) oaks (Q.brantii,Q. ithaburensis, andQ. macrolepis).Quercus euboica possibly originated by geographical isolation from a proto-trojana still close to the proto-cerris (Fig. 6). Interestingly, all theQ.euboica samples occurred in theQ.cerris- dominated ‘crown’ clade 4 (Fig. 3), hence, are part of the ‘occidentalis’ lineage. Likewise, the microspeciesQ. afares, Q. castaneifolia, andQ. look as well as the eastern replacement ofQ. trojana, Q. libani, were isolated from the mastercerris-trojana genepool(s). It is impossible to provide absolute dates for the final isolation events. Comparison with intra- and inter-species 5S-IGS divergence in the other two main lineages of western Eurasian oaks (‘ilicoid’ oaks of sect.Ilex, 5S-IGS is species-diagnostic; ‘roburoid’ oaks of sectionsQuercus andPonticae, 5S-IGS is largely undiagnostic;Denk & Grimm, 2010) indicates that the original split possibly predates diversification in the ‘roburoid’ oaks. The final establishment of species in western Eurasian members of Sect.Cerris is however as young as in the ‘roburoid’ oaks, may be ongoing (Q. trojana), and younger than the main split within the ‘ilicoid’ oaks, i.e., ofQ. ilex fromQ. aucheri (+Q. coccifera).

Phylogeography ofQuercus sectionCerris in western Eurasia

The primitive (L1) and derived (L2) haplotype groups of ‘Cerris-Ilex’ lineage, and two haplotypes representing the ‘Euro-Med’ lineage characteristic for western populations of sect.Ilex (Simeone et al., 2016;Vitelli et al., 2017), describe the main evolutionary trajectories of the western Eurasian lineage of sect.Cerris and their contact with already established lineages of sect.Ilex. According to coalescent theory, the most frequent and widespread haplotypes, i.e., H1 and H2, are likely ancestral (Posada & Crandall, 2001;Fig. 1). The close relationship between haplotypes H1 and H2–found all across western Eurasia (Fig. 2B) and in all taxa except the south-eastern species group,Q.brantii,Q.look, andQ.ithaburensis—andQ. phylliraeoides (Fig. 5; cf.Simeone et al., 2016,Fig. 1), the only species of sectionIlex extending into Japan—points towards a north-eastern Asian origin of sect.Cerris and a westward migration of a large population into the Mediterranean region. The revised (see Introduction) fossil record ofCerris in (North-)East Asia (starting from early Oligocene) predates earliest records in western Eurasia (Oligocene/Miocene boundary) by ca. 10 Ma, thus rejecting the hypothesis that sect.Cerris evolved from the western stock of sect.Ilex populations with ‘WAHEA’ haplotypes (Denk & Grimm, 2010;Simeone et al., 2016). The effective population size of the early west-migratingCerris must have been (very) large in contrast to their East Asian siblings. The East Asian species of sect.Cerris are more heterogenous (Chen et al., 2012;Zhang et al., 2015), and differ much more profoundly fromQ. phylliraeoides (the JapaneseIlex oak), but also from the ‘WAHEA’ haplotypes of sect.Ilex andQ. baroni (Fig. 1), considered early diverged plastid lineages of subgenusCerris (‘Old World’ or mid-latitude clade; the earliest diverged plastid lineage being the western Mediterraenan ‘Euro-Med’ type found inQ. ilex;Fig. 5;Simeone et al., 2016). In this context, assessment of the 5S-IGS diversity in (East) Asian members of sectionsCerris andIlex would be needed, to check whether the Asian counterparts are equally coherent or more diverse than their western Eurasian relatives.

Once established in the Mediterranean region (H1, H2;Fig. 2B), local bottlenecks may have contributed to increased genetic drift in the plastome in the eastern part of the range. A likely trigger are the complex orogenies shaping modern-day Turkey and the Levant, areas with an increased haplotype diversity including the most derived ‘Cerris-Ilex’ haplotypes (Figs. 1–2). This, and the general west-east differentiation pattern (see alsoFigs. 4 and5), parallels the situation in sect.Ilex,Q. coccifera in particular (Vitelli et al., 2017). A notable difference to sect.Ilex is the lack of plastid structuring (and diversity) in the central and western Mediterranean region, indicating a rather recent, singular colonization by the master population, clearly not affected by Oligocene micro-plate tectonics as suggested forQ. suber byMagri et al. (2007).

The derived L2 ‘Cerris-Ilex’ haplogroup (H5–H10) starts in Anatolia and extends further east (Iran) and south (Levant). In addition to isolation during range establishment, specialization to drier climates (e.g., summer-dry Mediterranean climates:Csa, Csb, Dsb) can be considered as trigger for increased genetic drift, possibly linked to speciation. The Aegilops oaks,Q. brantii, Q. ithaburensis, andQ. macrolepis, a well-circumscribed group based on 5S-IGS differentiation (Figs. 4 and5) and morphology, are unique by showing only derived ‘Cerris-Ilex’ haplotypes.

A remarkable exception are two ItalianQ.cerris individuals showing the derived haplotypes H9/H10, which occur in locations more than 2,000 km apart from other individuals of this Levantine haplotype sublineage. H9/H10 derive from types found in Anatolia and eastwards (Figs. 1 and2). Long-distance seed dispersal is highly unlikely. The main animal vector for propagation of oaks are the jaybirds, which are sedentary birds, with a short evasion range (<50 km;Pesendorfer et al., 2016). Man-mediated dispersal (in historic times) could be a likely explanation, although we note that haplotypes shared by disjunct central Mediterranean and the Anatolian regions were also found inQ.ilex, and possibly reflect the remnants of a pre-Quaternary continuous range.

In this context, the genetic diversity detected in the ItalianQ.trojana populations (both at the nuclear and at the plastid level) and the very limited, amphi-Adriatic distribution of haplotype H2 inQ.macrolepis (Italy, Albania;Fig. 2,File S1) likely confirm that these oaks are native in Italy. Similar close intra-specific phylogeographic relationships have been detected in other plant species on both sides of the Adriatic Sea (Musacchio et al., 2006;Hilpold et al., 2014), including oaks (Lumaret et al., 2002;Fineschi et al., 2002;López de Heredia et al., 2007;Bagnoli et al., 2016). In this case also, the Apulian populations ofQ.trojana andQ.macrolepis can be interpreted as the remnants of a once continuous ancestral range (Simeone et al., 2016), or witness a colonization wave that was likely favoured by land connections between the Balkans and southeastern Italy during the Messinian salinity crisis and (or) the Pleistocene glaciations (Nieto Feliner, 2014).