Introduction

The majority of Southern European forest communities are characterized by the dominance of species belonging to the genus Quercus L. The genus includes about 435 species throughout the world, but this figure is probably an underestimate (Denk et al. 2017). In Europe, there are about 34 species which are distributed across the whole continent. The most recent classifications based on molecular studies (Manos et al. 2008; Denk and Grimm 2010; Hubert et al. 2014; Denk et al. 2017; Hipp et al. 2019) have supported the hypothesis that the genus Quercus is divided into two monophyletic groups, which correspond to two different subgenera: subgen. Cerris and subgen. Quercus, the latter composed of five sections. The European white oaks belong to the section Quercus and originated from two waves of Palaeogene radiation from north America across both the north-Atlantic land bridge and the Bering land bridge (Denk et al. 2017).

The total number of oak species occurring in Italy is not known exactly, because the classification has been subject to continuous change over time. It ranges between nine and fifteen, depending on flora or checklist. This discrepancy concerning the number of accepted oak species names is caused, almost exclusively, by the descriptive and interpretative nebulosity that still characterizes the species complex of Q. pubescens (Table 1).

Table 1 List of the taxa belonging to the species complex of Quercus pubescens occurring in Italy on the basis of Pignatti (1982); Pignatti et al. (2017–19, Bartolucci et al. (2018), the most widely-respected classifications regarding the genus Quercus and the available National floras and/or checklists of neighbouring Mediterranean countries

The classification of oaks is based on a variable number of morphological traits of leaves and acorns, such as the degree of pubescence, its persistence on twigs and leaves, the dimensional traits of the leaf, the numbers, shapes and depths of the leaf lobes, length of the petiole, type of the margin and the surface of the cupule, among others. However, while the classification system enables reasonably clear separation among well-known white oak species (e.g. Q. faginea, Q. frainetto, Q. pyrenaica, Q. robur) within particularly critical groups (the species complexes of Q. petraea and especially that of Q. pubescens), taxonomists are facing considerable difficulties caused by more or less continuous patterns of morphological variation (see Bakis and Babaç 2014). This phenomenon, more prominent in the oaks than in other genera of European trees, is generally ascribed to a complex pattern of introgression and hybridization and to convergent morphological evolutionary processes in oaks (Oh and Manos 2008; Kremer et al. 2012; Denk et al. 2017). Such hybridization events occurred mostly in geographical areas such as southern Italy, where the thermophilous woody flora took refuge during cold periods of the Quaternary (Follieri et al. 1989; Petit et al. 2006; Magri et al. 2007; Leroy et al. 2019b).

Even if we exclude all binomials described from southern Italy and Sicily, which are now generally considered synonyms for other oak names, such as Quercus apennina Lam., Q. cupaniana Guss., Q. minae Lojac., Q. sicula Borzi ex Lojac., Q. tenoreana Borzì (cf. Gussone 1844; Borzì 1905, 1911; Lojacono-Pojero 1907, 1913–15), a high number of pubescent oak species is still reported in the most recent floras and checklists. According to Brullo et al. (1999a) and Pignatti et al. (2017–2019) the Euro-Mediterranean/SE-European Q. pubescens Willd. is absent from southern Italy where it is replaced by other deciduous thermophilous oak species with a steno-Mediterranean or Euro-Mediterranean distribution, for example, Q. congesta C. Presl., Q. dalechampii Ten. Q. virgiliana (Ten.) Ten., or southern Italian endemic distribution, for example Q. leptobalana Guss. and Q. amplifolia Guss. However, these taxa are not treated in the same manner by oak taxonomists in the rest of Europe. In fact, Q. pubescens is considered as a valid binomial throughout Europe, except for Spain, where it is, for nomenclatural reasons, substituted by the name Q. humilis Mill., of which the name Q. pubescens is considered a later synonym (Amaral Franco 1990). Concerning the other taxa, some of these are accepted in many European floras and checklists (e.g. Q. congesta and Q. dalechampii), while others (e.g. Q. virgiliana) are reported in some Southeastern-European floras, and Q. amplifolia and Q. leptobalana are considered endemic to southern Italy.

In addition to the lack of agreement on the number of white oak species, there are persistent cases of ‘misinterpretation’ that destabilize the Italian and European deciduous oak classification system. A typical example is Q. dalechampii, which although unequivocally typified and definitively assigned to the pubescent oaks group (Di Pietro et al. 2012), is in some cases still classified as belonging to the Q. petraea species complex. This is evident when looking at distribution maps that mix together records of pubescent and glabrous specimens under the name Q. dalechampii (see Euro+Med 2006−2014).

Several papers aiming to establish the possible occurrence of other pubescent oak species in addition to Q. pubescens Willd. have been published in the last decade for Southern Europe. Some of them were primarily based on the analysis of morphological characters, but also involved geometric morphometric analyses (Viscosi and Fortini 2011), micro-morphological leaf traits (Scareli-Santos et al. 2007; Fortini et al. 2009; Panahi et al. 2012; Scareli-Santos et al. 2013), fractal analysis of leaf morphology (Musarella et al. 2018), and molecular data (Curtu et al. 2007a, b, 2015). As far as southern Italy is concerned, two recent papers (Di Pietro et al. 2016, 2020) analysed 24 south-eastern Italian populations of pubescent oaks from a macro-morphological and molecular standpoint.

In the present study, we analysed the pubescent oak populations of Sicily and southern Calabria, which are crucial areas for the diversity of deciduous Mediterranean oaks. Sicily and southern Calabria host the highest number of pubescent oak species currently reported for the Italian Peninsula (Pignatti et al. 2017–2019; Bartolucci et al. 2018), the majority of their loci classici (Brundu et al. 2017; Peruzzi et al. 2019) and the highest number of pubescent oak forest phytosociological associations (Blasi et al. 2004). Furthermore, southern Italy provided refuge areas, both primary and secondary, during glacial periods which constitute an important reservoir of genetic diversity that has been preserved over time despite the strong anthropogenic pressure over the last three millennia (Médail and Diadema 2009). It should be noted that the reference to phytosociology is relevant to the taxonomic treatment if we consider that according to some authors (e.g. Guarino et al. 2015) the high coenological and floristic diversity found in the pubescent oak associations of southern Italy can be fully explained only by admitting the occurrence of ‘a blend of different species, partially blurred by systematic hybridization and introgression’ (ibid.).

The aim of this paper is to verify whether this high diversity, which involves several pubescent oak taxa, is supported by (i) the results of multivariate statistical analyses and (ii) sets of unambiguous and identifiable diagnostic morphometric characters.

Material and methods

Sampling

Thirteen pubescent oak populations were sampled throughout Sicily and southern Calabria. According to the floristic and phytosociological literature, these populations included the following taxa: Q. amplifolia, Q. congesta, Q. dalechampii, Q. leptobalana and Q. virgiliana. The collection sites were selected from the entire deciduous oak woodland bioclimatic belt, from sea level to the lower montane belt (Gianguzzi et al. 2016), and different types of geological substrates (Fig. 1 and Table 2). Where possible, priority was given to collection from populations located near to the loci classici of the species in question and to those already described from a syntaxonomic point of view in published phytosociological tables (Table 3). Leaves and fruits were collected from 27 to 33 adult trees per site in autumn 2017 and 2018. Up to three branches were selected randomly in the middle-upper part of the crown. Voucher specimens are stored in the Herbarium of the University of Molise (IS).

Fig. 1
figure 1

Spatial distribution of the thirteen pubescent oak populations sampled in the study area.

Table 2 Geographical location and syntaxonomic references of the Quercus populations sampled (ST). Data concerning dominant taxon, elevation and type of substrate of each site of collection were drawn from Brullo et al. (2009) and Gianguzzi et al. (2016)
Table 3 Description of leaf morphological characters

Ten leaves per tree specimen were pressed dried and scanned using an Epson GT-15000 scanner with a resolution of 600 dpi, with the abaxial surface facing upwards, and measured with ImageJ (Rasband 1997–2007).

In total, seventeen morphological characters were assessed (Table 3 and Online Resource 1): Seven primary characters (leaf area, leaf perimeter, leaf lamina length, petiole length, width of the major lobe, sinus width, widest point) were measured as described in Fortini et al. (2015). Five characters were derived characters (leaf compactness, obversity, lamina/petiole ratio, lobe depth ratio, lobe width ratio). One character (number of lobes on the left side) was numerical (counted), and one other (scoring of the basal shape of the lamina) was evaluated following Kremer et al. (2002). Finally, three characters (abaxial and adaxial laminar pubescence and twig pubescence) were observed on voucher specimens following Kissling (1980).

Three fruits (acorns) per tree specimen were measured using a digital electronic calliper MAURER 93110. In total, eleven fruit characters were assessed from voucher specimens (Table 4 and Online Resource 1), including four dimensional characters (cupule length, cupule width, acorn length, acorn width), two derived characters (cupule width/cupule length and cupule length/acorn length) and five categorical (dimensionless) characters: regularity of cupule edge, scale shape, gibbosity (tuberculate bases of cupule scales), cupule internal hairiness and type of exterior of cupule edge.

Table 4 Description of fruit morphological characters

In total, 3,887 leaves and 1,047 fruits were collected from 391 tree specimens. Since not all sampled tree specimens bore fruit, the final and complete matrix subjected to statistical analysis was composed of 333 specimens.

Because the dichotomous keys of Flora D’Italia (Pignatti et al. 2017–2019) takes into consideration the appearance of the bark ribs, this character was also observed for all the sampled oak specimens to aid identification (see Online Resource 2).

Data analysis

Partition of the total variance of leaf morphology was evaluated through the two-levels nested analysis of variance (nested ANOVA) which gave rise to results readable at different hierarchical level: (a) inter-population (variance among different populations), (b) intra-population (variance among the trees of a single population) and (c) individual (variance among the leaves of a single tree). A spreadsheet suitable for this purpose (McDonald 2009) was modified to analyse more than fifty subgroups, with up to 200 observations per subgroup, and used subsequently to perform the nested ANOVA (with Satterthwaite approximation, because of the unequal sample sizes of the subgroups) on the entire data set of 3,887 leaves and fourteen leaf variables.

The matrix of 333 tree specimens and 28 leaf and fruit variables was analysed. To identify possible diversity groups, an exploratory analysis was performed on standardized data through Hierarchical clustering with PAST v. 4.03 using Euclidean distance, Ward’s minimum variance method (Hammer et al. 2001); Hierarchical clustering with JASP v. 0.13 (estimation method ‘moment’, model optimized with respect to the silhouette value (Jasp Team 2020).); fuzzy k-means clustering with XLSTAT vers. 2020.3.1 (coefficient of fuzziness = 1.1, method: cosine, 20 repetitions, not weighted, initialization method: K||); k-means clustering with JASP vers. 0.13 (estimation method ‘moment’, model optimized with respect to the silhouette value).

A principal component analysis (PCA) was performed with PAST and the results showed with the superimposition of the groups highlighted by the cluster analyses.

Furthermore, using XLSTAT 2020.1.1 (ADDINSOFT 2020), the consistency of the groups identified by the four cluster analyses was tested with the Kruskal–Wallis test (P calculated with the Montecarlo method, 1000 resamplings) and then by the multiple pairwise comparisons (Dunn’s procedures with Bonferroni correction).

Finally, a linear discriminant classification (LDC, estimation method ‘moment’) using JASP was performed.

To confirm the result of the a-priori classification based on the guide species of the oak communities investigated, an expert classification of all the tree specimens collected was also performed based on the analytical key for oak identification published in the new edition of Flora d’Italia (Pignatti et al. 2017–2019). In some cases of particularly critical specimens, it was necessary to do some forcing to complete the identification process. When the incompatibility among diagnostic characters turned out to be too pronounced and would have required an excessive degree of subjectivity to complete the identification process, the specimen was classified as ‘preliminarily indeterminate’. When we found a lower degree of incompatibility among diagnostic characters, we decided to choose one of the two options proposed by the identification key arbitrarily. Specimens identified in this way were subsequently classified adding the epithet ‘doubtful’ (e.g. Q. virgiliana doubtful, Q. congesta doubtful, etc.).

Our expert identifications of the specimens were compared with what was expected for the sampling sites based on the phytosociological tables of the associations already published for those areas or for strictly neighbouring areas. All the pubescent oaks occurring at each collection site were collected by randomized sampling and not selectively directed toward the collection of one or other putative species. In this context, the expected values of the various oak species were obtained by making a correlation between the number of specimens we collected and the values of frequency and coverage of the oak species reported in the phytosociological tables of the pubescent oak forest associations published for those areas.

Results

The nested ANOVA showed that the morphological character variability observed among the leaves of an individual tree (ranging between 42.40% and 89.73%) accounted for the highest portion of the total variability. This was followed by the variability observed among the leaves of the tree specimens belonging to the same population (9.37%–46.73%) and by that observed among the leaves of the specimens belonging to different populations (0.83%–17.72%; Online Resource 3).

The hierarchical cluster analysis of the matrix based on leaf and fruit characters produced a dendrogram showing two main groups of oak specimens (Online Resource 4a). The occurrence of two groups was also revealed by the other three cluster analyses performed (Online Resource 4b–e).

By making a pairwise comparison between the different partitions obtained by the cluster analyses on the basis of the presence of the tree specimens in the two groups and regardless of the order of the groups, an agreement always higher than 80%, was found (Online Resource 4f,g).

The relatively low degree of dissimilarity between the two main groups identified by the cluster analysis may be observed in the PCA diagram with partition superimposed (Fig. 2). The two clusters exhibited a significant overlap when considering the grouping based on the hierarchical classifications whereas the superimposion of k-means and the fuzzy k-means clustering showed a clearer separation along PCA axis I (20.03% of the variance) (Online Resource 4h). The variables correlated to the first PCA axis are those pertaining to leaf dimensions (Area, Per, LL, LW, WP), while the variables correlated to the second PCA axis are those concerning the dimensions of the fruit (CW, AL, AW). Summary tables from Kruskall–Wallis test (α = 0.05) showed that the separation of the two clusters was clearer using the fuzzy k-means and the k-means clustering (Online Resource 4i) for which the only insignificant variables were Ob, LWR, CL, RCE, TCE and Co, LWR, CL, RCE, TCE, respectively. The linear discriminant classification LDC attributed to the fuzzy k-means clustering the best value of the accuracy test (Online Resource 4l–m) and for this reason we decided to accept the two groups identified by fuzzy k-means clustering for further analysis and discussion.

Eleven populations out of the thirteen investigated contained specimens occurring in the both main two clusters identified by the fuzzy k-means clustering (Fig. 3). Populations ST11 and ST13 were the only populations where all specimens belonged to the same main clusters (Group 1). In four populations (ST01, ST04, ST07 and ST09), one of the clusters (Group 2) was dominant, comprising more than 70% of the specimens. All the other oak populations investigated were mixed with less than 70% of the specimens in one or the other main cluster.

From the box-plot diagrams, produced from the Kruskall-Wallis test (α = 0.05), and the multiple pairwise comparisons using Dunn's procedure, it was found that the separation of the two main clusters of the fuzzy k-means clustering was mostly related to leaf dimensions (Area, LL, LW, Per, PL, SW and WP), lobe depth ratio (LDR), petiole ratio (PR), number of lobes (NL), basal shape of the lamina (BSL), compactness (Co), and leaf pubescence (AB-PU, AD-PU and TW-PU). Among fruit characters, cupule width (CW), acorn length and width (AL and AW), cupule form (CW/CL), cupule internal hairiness (CIH), and cupule length/acorn length ratio (CL/AL), scale shape (SS), tuberculate bases of the cup scales (G), and cupule internal hairiness (CIH) differentiate slightly between the two groups (Online Resource 5).

The expert identification of the 333 samples, achieved through rigorous use of the analytical keys of the most recent edition of Flora d’Italia (Pignatti et al. 2017–2019), proved tortuous. In fact, we found impossible to assign thirty individual trees out of the 333 sampled whereas 86 specimens were classified as doubtful.

The labelling of the dendrogram with the names derived from our expert identification did not show the separation into groups containing specimens belonging to the same putative species (Online Resource 4a). Also in the PCA there is a marked dispersion of the specimens belonging to the same putative species (Fig. 2a).

Fig. 2
figure 2

Principal component analysis (PCA) based on the 28 morphological traits. The first axis explains 20.028% of the total variance, while the second explains 12.114%. a – The two groups obtained by the fuzzy k-means clustering are emphasized by the convex hulls (Group 1 on the left and Group 2 on the right). The oak taxa listed comprise both specimens that were positively identified and specimens which remained indeterminate (dots – Quercus dalechampii, Q. congesta, Q. ichnusae, Q. leptobalana, Q. pubescens, Q. virgiliana; Squares – doubtful specimens). b – Projection of the original axes (variables). c – The thirteen populations sampled emphasized by convex hulls: ST01 (red), ST02 (fuchsia), ST03 (grey), ST04 (orange), ST05 (dark green), ST06 (black), ST07 (blue), ST08 (yellow), ST09 (purple), ST10 (light blue), ST11 (green), ST12 (sky blue), ST13 (dark golden).

The comparison of the number of specimens identified and the number of specimens expected on the basis of the published phytosociological papers (Table 5) showed that for Q. congesta (36 specimens identified with certainty and 85 with doubtful specimens included) and Q. dalechampii (87/101) the number of specimens identified (when the doubtful specimens are included) more or less correspond to what was expected (64 and 93 specimens, respectively). More than half of the doubtful specimens in the whole dataset belong to Q. congesta, which is the only taxon that showed a higher number of doubtful specimens than identified specimens. However, more than half (21/36) of the specimens identified with certainty as Q. congesta came from sites where this species was not expected while this ratio was significantly lower for Q. dalechampii (21/87). Two out of the three collection sites where unexpected Q. dalechampii specimens were found (ST6 and ST2, both referring to the association Oleo-Quercetum virgilianae), showed Q. dalechampii as the most abundant species. Only twelve specimens were identified with certainty as Q. virgiliana (14 including doubtful specimens) whereas the number expected was 110. The number of specimens identified with certainty as Q. leptobalana was 39 (56 including doubtful specimens). Only fourteen specimens out of the 39 Q. leptobalana specimens identified were collected in the two sites reported as Quercetum leptobalanae (4 specimens in ST1 and 10 in ST11), while the remaining specimens came from other collection sites in Sicily and southern Calabria at which no Q. leptobalana was expected. No specimens were identified as Q. amplifolia, despite 38 specimens being expected, whereas we identified 35 specimens of Q. pubescens and seven specimens of Q. ichnusae, although these two taxa are currently not reported from the study area (cf. Brullo et al. 1999a, b; Pignatti et al. 2017–2019).

Table 5 Comparison between the number of specimens identified in the thirteen collection sites on the basis of the diagnostic keys published in (Pignatti et al. 2017–2019) and the number of specimens expected (in bold) on the basis of frequency and cover values reported in the phytosociological tables published in Brullo and Marcenò (1985) – (A) and Brullo et al. (2001) – (B). For each taxon (columns) both the number of specimens identified with certainty (normal font) and that of doubtful specimens (italic font) are reported. In the last two columns the number of indeterminate specimens (grey column) and the total number of specimens collected are shown. (Q. lept = Quercetum leptobalanae; Oleo-Qv = Oleo-Quercetum virgilianae; Fe-Q.con = Festuco heterophyllae-Quercetum congestae; Ar-Q.con = Arabido turritae-Quercetum congestae; Er-Qv = Erico arboreae-Quercetum virgilianae; Er-Qc = Erico arboreae-Quercetum congestae).
Fig. 3
figure 3

Percentage of distribution of pubescent oak specimens in the two groups identified by fuzzy K-means clustering for each of the thirteen populations sampled.

Discussion

Over recent years a major international effort has been made to revise the entire Quercus genus, which, worldwide, had reached the number of 500–650 species (Govaerts and Frodin 1998). In the taxonomical framework of Southern Europe the Q. pubescens species complex remains an unresolved tangle. Papers published in Southern Europe in the last fifteen years regarding the morphological variability of the European white oaks have been moving in more or less the same direction: that a critical evaluation of the many oak taxa is needed regardless of whether they were included in the national floras or checklists as accepted species or not (Franjic et al. 2006; Enescu et al. 2013; Wellstein and Spada 2015; Di Pietro et al. 2016, 2020).

The concept of species in Quercus

Does it still make sense to classify oak populations to different species when it has been known for at least forty years that the genus Quercus is subject to interspecific hybridization? Actually, a definitive answer to this question has never been provided and is well beyond the scope of this paper. We incline to think that the answer is yes. The topic is still widely debated and a ‘measured scepticism about oak species is not uncommon among botanists even today’ (Hipp 2015). Already in the first half of the 1970s some authors (e.g. Van Valen 1976; Burger 1975; Hardin 1975) had argued that if gene flow among different oak species might exceed gene flow between populations of a single species, then the concept of species in oaks needed to be significantly rearranged. Accordingly, these authors proposed that oak species should be considered as ecological-taxonomic entities rather than reproductively defined ones. More recently, the availability of increasingly refined molecular markers has led to this issue being reconsidered from a different perspective. Certainly, hybridization may have been one of the factors underlying the palaeogenic radiation of oaks and their postglacial spread northwards. Indeed, recent population genomic studies have defined hybridization as a microevolutionary process able to reinforce adaptation and enhance migration, thus improving the ability of individuals to colonize new habitats (Cavender-Bares 2019; Leroy et al. 2019a, b; Kremer and Hipp 2019). However, most researchers in the field argue that hybridization does not lead to the fragmentation of the integrity of species if even a small part of the genome is able to maintain species barriers (Eaton et al. 2015). Particular genes (e.g. those under selection for drought-tolerance or photoperiodic control of growth) were found to be highly contributing to the maintenance of genetic divergence between ecologically distinct species, suggesting that selection manages to cope with gene-flow (Lind-Riehl et al. 2014; Oney-Birol et al. 2018). Much of the variation in oaks that has been interpreted as evidence of hybridization in the past can more properly be considered as morphological variation within highly variable species attributable to intraspecific phenotypic plasticity (i.e. the ability of a species to produce different phenotypes with varying environmental conditions). Furthermore, it has been demonstrated that ‘gene flow among species is normally overwhelmed by gene flow within species’ (Gerber et al. 2014). These and other arguments led Hipp et al. (2015) ‘to say confidently against the backdrop of the last 20 years of oak research, that oak species are real’.

Several recent works using genome-scale data have been successful in resolving relationships throughout the Quercus phylogenetic tree, from the definition of the two major clades worldwide (Denk et al. 2017; Hipp et al. 2015, 2018, 2019; Crowl et al. 2019), up to the fine-scale infrasectional level (Olalde et al. 2002). However, as underlined in Denk et al. (2017), the major challenge for future studies is to establish the molecular and morphological circumscription at species rank and their biogeographic and ecological characterization. In the last two decades important results have been achieved in identifying the occurrence (or the absence) of species boundaries between taxonomically critical groups of sympatric Southern-European white oaks using molecular markers (Bruschi et al. 2000; Jedináková-Schmidtová et al. 2004; Franjic et al. 2006; Fortini et al. 2009; Curtu et al. 2011; Enescu et al. 2013; Crăciunescu et al. 2015; Di Pietro et al. 2020). In our opinion, the necessary further step should be to link discontinuities highlighted at the genomic level to a well-defined and unambiguous set of diagnostic morphological, ecological and coenological characters to circumscribe oak species. Also, proper species names should be found, carefully applying the rules of the International Code of Nomenclature (Turland et al. 2018). As a consequence, a name linked to a type specimen is also linked to a certain geographical location (locus classicus) and possibly (but with a significantly lower degree of certainty) to a larger geographical area. Unfortunately, correspondence between provenance and nomenclatural identity of the collected material tends sometimes to be disregarded, leading inevitably to mistaken or deviated conclusions. This paper, even though restricted to a relatively small geographical area, represents an essential step in this direction.

The current state of the art

What has emerged from the most recent taxonomic investigations of European white oaks on the basis of morphological characters (micro- and macro-) is that the current ‘state of art’ of identification of pubescent oak specimens does not allow unequivocal identification at the species rank. Once one arrives at the objective dichotomy separating ‘pubescent-twig oaks’ (Q. pubescens species complex) and ‘glabrous-twig oaks’ (Q. petraea species complex), the subsequent steps in taxonomic discrimination become markedly subjective (for both groups) and therefore not universally shared. As regards Italy, we are fairly sure about the identification of the four main species of white oaks, namely Q. frainetto, Q. petraea, Q. robur and Q. pubescens. As regards the first three of these taxa, there have never been serious doubts that they had to be accepted as valid species and only a few others considered morphological variations within them to deserve the rank of variety or, at most, subspecies (see Q. petraea subsp. austrotyrrhena Brullo, Guarino and Siracusa or Q. robur subsp. brutia (Ten.) O. Schwarz). Only in the Quercus petraea group there is no general agreement on the validity of some taxa (e.g. Q. petraea subsp. polycarpa and Q. petraea subsp. iberica) and on the choice of the name to refer to what was previously called Q. dalechampii (cf. Bock and Tison 2013, Di Pietro et al. 2012; Raab-Straube and Raus 2013; Kučera 2018). This is in no way comparable, however, with the taxonomic and nomenclatural chaos in which the group of Q. pubescens is immersed. Indeed, the well-known morphological variability of this group has often been used to advance proposals to split it in a multitude of narrowly defined taxa. The majority of these ‘narrow endemic eco-morphotypes’ (see Pasta et al. 2016) are generally considered synonyms of Q. pubescens. In some rare cases the opposite path occurs, that is, populations attributed preliminarily to Q. pubescens that have subsequently acquired taxonomic autonomy. An example may be Q. kotschyana O. Schwarz, an endemic species of Lebanon previously included in Q. pubescens as Q. pubescens subsp. pinnatifida Gmelin (Stephan and Teeny 2017; Hipp et al. 2019) or considered an hybrid between Q. pubescens and Q. cerris under the name Quercus ×baenitzii (see Camus 1936–54 and http://oaks.of.the.world.free.fr/quercus_baenitzii.htm), or Q. pubescens subsp. crispata (already accepted in Schwarz 1993).

Among the Southern European countries, Italy is certainly the one where the splitter tendency is more deeply rooted. Accordingly, several Italian authors keep on referring several pubescent oak taxa other than Q. pubescens in their floristic, taxonomic and phytosociological contributions (cf. Brullo and Marcenò 1985; Brullo et al. 1999b, 2001; Blasi et al. 2004; Brullo et al. 2009; Casavecchia et al. 2017). On the other hand, the group of ‘lumper taxonomists’ tend to reduce the number of white oaks occurring in the Italian peninsula to four species only.

In our opinion, the problem is partly conceptual and partly methodological. The wide morphological variability exists not only in the pubescent oaks group as a whole, but also within populations of the same putative pubescent oak species. The ‘splitters’ consider unlikely that such a wide variability could be displayed by a single species. On the other hand, there is evidence that when collections include a large number of populations and a large number of specimens per population, the variability tends to follow more or less continuous pattern in nearly all the diagnostic morphological characters. This gradual morphological pattern was clearly demonstrated in Di Pietro et al. (2016) in the Apulia region, who covered the entire bioclimatic and lithological variability of the study area and collected an adequate number of populations and specimens per population. In this study, no groups of morphological and molecular diversity and no combination of diagnostic characters were found to identify the putative species reported to dominate the study sites, such as Q. dalechampii, Q. pubescens and Q. virgiliana (Biondi et al. 2004; Di Pietro and Misano 2009). Indeed, a more recent study using fractal analysis applied to the leaves of Calabrian deciduous oaks (Musarella et al. 2018) was also unable to distinguish different species inside the Q. pubescens group.

In the present research on the Calabrian and Sicilian oaks, we decided to further improve the sampling protocol by increasing the number of specimens per site and to include the loci classici of some of the species sampled. No consistent groups were found at the morphological level whereas, according to the published phytosociological studies, we should have found Q. amplifolia, Q. congesta, Q. dalechampii Q. leptobalana and Q. virgiliana. Taking the high frequency and cover values of these species in the original phytosociological tables into account, it would have been unlikely to avoid collecting a sample of one or more of these species given the randomized sampling protocol we adopted. The distribution of specimens into the two main clusters identified by cluster analyses was found to be rather random and unrelated to ecological or geographical characters. For example, specimens we collected from a putative thermo-Mediterranean and xerothermic Q. virgiliana–Q. amplifolia community on limestone were found to group together with specimens of a putative Q. congesta–Q. dalechampii community from the montane belt of the Etna Volcano where neither Q. virgiliana nor Q. amplifolia were reported. As a further example, we would have expected that specimens collected from the only two sites known to host Quercetum leptobalanae (cf. Brullo 1984; Brullo and Marcenò 1985; Gianguzzi and La Mantia 2004) would have grouped together. Instead, all the specimens from the site ‘Piano Zucchi’ (locus classicus of the Quercetum leptobalanae) clustered in Group 1 whereas two thirds of the specimens of Quercetum leptobalanae from ‘Bosco Ficuzza’ clustered in Group 2.

The morphological characters that were diagnostic for the two clusters deriving from fuzzy k-means analysis seem to suggest the existence of two morphotypes. One (Group 2) is characterized by greater dimensional characteristics of leaves (as regards the leaf lamina, petiole length, depth of the lobes) and bigger acorns, while the second (Group 1) has smaller and more compact leaves and smaller acorns. However, observing the box plots (Online Resource 5), the situation is much more complex and all the characters show a wide variation range. From a taxonomical point of view, the two morphotypes are difficult to match to any of the putative oak species that were the object of this study, neither in respect of their original descriptions nor of those proposed in analytical keys (Brullo et al. 1999a, b; Pignatti et al. 2017–2019). In fact, many of the morphological characters we found to discriminate between the two major clusters are dimensional characters. These are known to exhibit a wide range of variation depending on the site environmental conditions and on seasonal variations of the climate factors and may change from one year to another (cf. Škvorc et al. 2005; Bonito et al. 2011; Martiník et al. 2014). Theoretically, the use of many morphological characters in statistical analyses could have caused excessive background noise caused by characters that are non-diagnostic. This possibility is, however, largely refuted by the recent taxonomic literature, which tends to use the greatest number of available characters to discriminate at the interspecific and intraspecific levels. Moreover, it has been shown that if populations of other white oak species, such as Q. petraea and Q. frainetto, are included in statistical analysis together with populations of Q. pubescens s.l., the specimens belonging to the two former separate quite clearly from the undifferentiated bulk of specimens belonging to Q. pubescens s.l. (Viscosi and Fortini 2011; Fortini et al. 2015).

Aware of the difficulties in identifying and classifying white oaks, we also opted to perform our expert identification of the specimens using the analytical keys proposed in Pignatti et al. (2017–2019) and to compare results of our classification with cluster analyses and the PCA (Online resource 4a,h). We expected the identified species to correspond to some ecological factors (e.g. bioclimatic belt, type of substrate) or at least to the species reported in the published phytosociological tables. These expectations were largely unmet (Table 5). For example, only twelve specimens were ascribable with certainty to Quercus virgiliana whereas we would have expected the number to be not less than 100. A total of 39 specimens were attributable to Q. leptobalana but only fourteen specimens (36%) came from the two sites of Quercetum leptobalanae (4 specimens in ST1 and 10 in ST11), while the remaining specimens were distributed in other collection sites. Such results would be new sites for Q. leptobalana in Sicily and first records of this taxon for continental Europe. In contrast, no specimens classifiable as Q. amplifolia were found even though, according to the phytosociological literature, this species should occur in at least 60% of the communities we sampled. Finally, we found specimens attributable to Q. pubescens (35) and Q. ichnusae (7) although the former species is considered to be absent in southern Italy and the latter endemic to Sardinia Island (Mossa et al. 1999).

It is possible that in southern Italy and Sicily the morphological variability among pubescent oaks is greater than in other countries, and this fact led botanists from the early nineteenth century (Presl 1822, 1826; Tenore 1830, Tenore 1835–36; Gussone 1844; Borzì 1905, 1911; Lojacono-Pojero 1907, 1913–15) to the present days (Brullo et al. 1999a, b; Giardina et al. 2007; Pignatti et al. 2017–2019) to consider the local pubescent oak populations as a multi-species complex. In opposite, we suggest that this high morphological variability (if really higher than in similar other areas at all) may be explained by biogeographic history processes and hybridization. The Apennine Peninsula is universally considered one of the most important refuge sites for the thermophilous oak forest during the Quaternary cold periods (Follieri et al. 1989; Hewitt 1999; Allen et al. 2000; Petit et al. 2002).

In particular, the southern Apennines and the major islands suffered only in a very moderate way the effects of the glaciations, compared with the central and northern Apennines (see Acquafredda and Palmentola 1986; Di Pietro et al. 2017), so extinction of genotypes and populations was minimized. As suggested by Nieto Feliner (2011), adverse periods of climatic oscillations may have given rise to a low spatial scale dynamism with species movement ‘across a wide possibility of patches and altitudes with more source diversity at genotype and population level’. Southern Italy oak species had thus higher chance to be in close contact and possibly to produce hybrids than in more northern areas. At this point, the question raises which oaks are more likely to have hybridized. Our opinion is that the most likely hybrids are those between a pubescent oak (e.g. Quercus pubescens) and other white oaks today unquestionably recognized as valid species (Q. petraea, Q. robur and Q. frainetto). Indeed, Fortini et al. (2015) identified presence of genomes of Q. frainetto, Q. petraea and Q. pubescens (in variable quantities) in several oak specimens collected in mixed oak woods of southern Italy, showing intermediate morphological characters.

The results of the present research suggest that all the populations we collected can be assigned to a single species of pubescent oak characterized by high morphological variability and presumably a wide ecological amplitude, too. A name that would correspond perfectly to the identity of this morphologically and ecologically highly variable species is, in fact, Q. pubescens Willd. This statement, despite being in contradiction with the taxonomic and phytosociological literature concerning southern Italy cited above, is actually in accordance with what has already been published for the whole Southern Europe: it is impossible to trace reliable species boundaries within the complex of pubescent oaks (cf. Franjic et al. 2006; Enescu et al. 2013; Wellstein and Spada 2015; Di Pietro et al. 2016, 2020). More information useful for unravelling the taxonomic-nomenclatural tangle of oaks can be derived from comparative studies combining morphological and molecular analyses. For example, interesting results on the association between one genomic region and leaf traits were recently obtained in a study on two interfertile and partially sympatric red oak species of North America (Gailing et al. 2018). As regards the Italian white oaks, the preliminary results of the genetic analysis carried out on the same set of populations and specimens as used in this paper (Di Pietro et al. submitted) show significant concordance between molecular and morphological results in terms of an absence of genetic groups of variation. Similar comparative studies have already been published for other Italian areas (Fortini et al. 2015; Di Pietro et al. 2016, 2020) and have also suggested that occurrence of more than one species of pubescent oaks in those study areas was unlikely.

Conclusions

Sicily and southern Calabria form a very important phytogeographical district, which is located in the centre of the Mediterranean Basin and is notoriously populated by a high number of endemics, relict and rare species (Brullo et al. 2011; Sciandrello et al. 2015; Spampinato et al. 2018). As far as the genus Quercus is concerned, this area is known to host five pubescent oak species, four of which were typified there. These species, in addition to all being accepted in the latest edition of the Flora of Italy and partly in Flora Europaea and in other European National floras and checklists, were used as name-giving species in several phytosociological syntaxa and considered diagnostic for some Forest Habitats included in the 92/43/EEC Directive (Biondi et al. 2009; European Commission 2013).

The morphological analysis of specimens from thirteen pubescent oak populations of Sicily and southern Calabria performed both without a priori identification and through expert identification based on the published analytical keys, led to results that were incompatible with each other and with published data on pubescent oak species occurring in the study area and their distributional and ecological ranges. It seems more probable that the morphological variability can be assigned to a single species that is characterized by a wide phenotypic and ecological amplitude.