Morphological and molecular results from a geographical transect focusing on Quercus pubescens/Q. virgiliana ecological-altitudinal vicariance in peninsular Italy

Abstract Quercus pubescens and Q. virgiliana are the most cited taxa in the Italian deciduous oak forests floristic and phytosociological literature. According to some authors, Q. pubescens is typical of inland areas and higher altitudes whereas Q. virgiliana restricted to the coastal plain and the hilly belt. Seven pubescent oak populations distributed along an altitudinal transect in central Italy have been analyzed from a coenological, morphological and molecular point of view. The vegetation sampling was carried out using the phytosociological approach. Morphological variation of tree individuals was analyzed using 14 leaf traits, while leaf shape variation was investigated using the Procrustes ANOVA. Genetic analyses were carried out through twelve EST-SSRs markers. Results: only the leaf pubescence exhibited discriminating power among all the leaf morphological traits considered. Very low differences in the leaf-shape emerged from geometric morphometric analysis. Genetic analyses did not evidence statistically significant clusters. Bayesian analyses including data from genetically pure populations of Q. pubescens, Q. petraea and Q. frainetto assigned all the seven populations investigated to Q. pubescens. Neither the morphological nor the genetic results allowed to identify specimens attributable to Q. pubescens or Q. virgiliana nor to highlight a possible ecological-altitudinal vicariance between these two species.


Introduction
The pubescent white oak species (subgenus Quercus, section Quercus) are among the most important components of the south-European xerothermic deciduous oak forests (Mucina et al. 2016). In Italy, these are widespread from the coastal areas to the lower montane belt often as dominant species in final successional stages. Scientific interest in xerothermic deciduous oak forests has been increasing in recent years, especially in relation to topical issues such as reforestation policies, climate change and biodiversity safeguard. In fact, it is likely that xerothermic pubescent-oaks will progressively take over more mesophilous tree species if the increase in temperatures and drought periods will come true (Negi and Negi 2021;Petritan et al. 2021). Owing to their resistance to drought stress (Nardini and Pitt 1999;Saunier et al. 2018), xerothermic pubescent oaks could be an ecological answer for the arrangement of urban forest or rural reforestation (Pasta et al. 2016). Finally, xerothermic pubescent oak forests exhibit one of the highest specific richness among all temperate European deciduous broadleaf forests, both in the woody layers and in the undergrowth (COST Action BOTTOMS-UP 2020). However, whether it is for the protection and management of forest ecosystems or for the selection of tree species for environmental planning, a necessary pre-requisite would be to know the taxonomic identity of the pubescent oaks we are dealing with. Any assessment for ecological optimum, drought resistance, degree of resilience or of any other ecological parameter, to be exportable and replicable, cannot ignore a previous correct taxonomical classification. In fact, the uncertainty which involves identification processes in the pubescent oaks collective group remains a big problem for botanists (Wellstein and Spada 2015;Pasta et al. 2016). Italy hosts the locus classicus of various pubescent oak species, such as Quercus. amplifolia Guss., Q. congesta C.Presl, Q. dalechampii Ten., Q. ichnusae Mossa, Bacc. & Brullo, Q. leptobalana Guss., Q. virgiliana (Ten.) Ten., most of which are accepted in important Italian Floras and Checklists (Pignatti 1982;Brullo et al. 1999;Pignatti 2017Pignatti -2019, international institutions (e.g. The Plant List, IPNI, Euro + Med PlantBase) as well as in phytosociological syntheses (Brullo et al. 2001;Bacchetta et al. 2009;Biondi, Casavecchia, Beccarisi, et al. 2010). In the last two decades, several studies aimed at the taxonomic, and genetic circumscription of southern European pubescent oaks were carried out (Bordács et al. 2002;Bruschi and Grossoni 2004;Škvorc et al. 2005;Franjic et al. 2011;Fortini et al. 2009;Enescu et al. 2013;Curtu et al. 2015;Di Pietro et al. 2016;Di Pietro et al. 2021;Proietti et al. 2021). These studies highlighted the current difficulties in matching the rich-in-species taxonomic and nomenclatural classifications proposed by Floras and Checklists with a corresponding clustering based on morphological, and/or genetic characters. If a rich biometric morphological database associated to genetic information is being formed for southern Italy (Di Pietro et al. 2016, 2021, the same cannot be said for central Italy for which area there are only two genetic papers on Q. pubescens Willd. analyzed together with other white oaks (Bruschi et al. 2000;Fineschi et al. 2002). On the other hand, numerous are the studies concerning the phytosociology and syntaxonomy of the xerothermic pubescent oak forest communities in central Italy (Ubaldi et al. 1984;Blasi and Di Pietro 1998;Allegrezza et al. 2002;Blasi et al. 2004;Ballelli et al. 2006;Facioni et al. 2015). In some of these studies (those concerning the inner areas) the only pubescent oak species listed in the phytosociological tables is Quercus pubescens Willd., whereas in other studies (those concerning the thermophilous pubescent oak woods of the hilly belt) Quercus virgiliana (Ten.) Ten. is reported as the dominant oak species (Biondi et al. 2004;Rosati et al. 2005;Allegrezza et al. 2006;Biondi et al. 2008;Di Pietro and Misano 2009;Blasi and Biondi 2017). In Biondi, Casavecchia, Pesaresi, et al. (2010), an explicit proposal was advanced to replace the name of some associations validly published as Quercetum pubescentis with the epithet Quercetum virgilianae. In fact, there is not agreement on the effectiveness of Quercus virgiliana as a valid name (Wellstein and Spada 2015;Grossoni et al. 2021). The reasons for this lack of agreement are substantially two. The first concerns the inconsistencies between the protologue of Q. virgiliana made by Tenore (1835Tenore ( -1836 and the original material specimens deposited in various European Herbaria. The second is the lack of well-established diagnostic morphological traits useful to distinguish Q. virgiliana from Q. pubescens. In some classification frameworks Q. virgiliana is considered a distinct species based on some morphological characteristics such as the flattened apex of the leaf and the irregular edge of the cupule (Pignatti 2017(Pignatti -2019. In other classifications it is considered a doubtful species (Bussotti and Grossoni 1997;Govaerts and Frodin 1998;Conti et al. 2005;Pasta et al. 2016), or a synonym of Q. pubescens Willd. subsp. pubescens, (Govaerts and Frodin 1998;Bartolucci et al. 2018; Oaks of the world (http://oaks.of.the.world.free.fr); Euro + Med PlantBase (http://ww2.bgbm.org/EuroPlusMed)).
In the European taxonomic literature, there are several studies aimed to resolve the dualism between Quercus pubescens and Q. virgiliana but the results are often conflicting. For instance, in Croatia, Trinajstić (1974;2007) confirmed the taxonomical separation at species rank between Q. pubescens and Q. virgiliana based on a high number of morphological traits concerning leaf, flower, cupule and bark. In contrast Škvorc et al. (2005) andFranjić et al. (2006), for the same country, were more likely to consider only Q. pubescens due to the lack of an effective genetic differentiation between the latter and Q. virgiliana. In Turkey, Borazan and Babac (2003) virgiliana as a hybrid between Q. petraea (Matt.) Liebl. and Q. pubescens. Ballian et al. (2011) referred only to Q. pubescens in a study on chloroplast variability of the Q. pubescens complex in the western Balkans. Jerse and Batić (2007), Sofletea et al. (2011), andEnescu et al. (2013) for Romania likewise included Q. virgiliana in the range of variability of Q. pubescens based on morphological and molecular analyses.
As far as peninsular Italy is concerned, some authors (Brullo et al. 2001;Biondi, Casavecchia, Beccarisi, et al. 2010;Rosati et al. 2005;Filesi et al. 2010;Blasi and Biondi 2017) argue that the true Q. pubescens occur only in inland areas or at high altitude while in the coastal and hilly belt (up to the lower montane belt, as reported in Fascetti et al. 2013) it would be completely replaced by Q. virgiliana. To verify this hypothesis, we have analyzed different populations of pubescent oaks selected along an altitudinal gradient moving from the Tyrrhenian coastal side to the core of the Apennines in central Italy where the occurrence of the two species was reported in taxonomic and phytosociological literature (Moraldo et al. 1990;Brullo et al. 1999;Anzalone et al. 2010;Facioni et al. 2015). We integrated molecular data with morphological and morphometric leaf data, to identify possible gradients of leaf phenotypic variation associated to clusters of genetic diversity.
It is important to point out that the exact identification of the white oaks occurring in the European forest landscape has implications that are not limited to mere taxonomic issues. In fact, Quercus pubescens and Q. virgiliana are currently considered diagnostic species for Directive 92/43/EEC habitat 91AA* Eastern white oak woods for central Italy (Devillers et al. 1991;Biondi et al. 2009; Commission of the European Communities 2014).

Study area
The study area (Figure 1) is located in central Italy between the administrative regions of Lazio and Abruzzo (13.05723 Figure 1. map of the study area in central-italy and location of the pubescent oak sampling stands along an altitudinal transect. rectangles circumscribe hilly belt wood stands, pentagon the submontane belt stand and circles the montane belt stands. the site of mount Vairano (which will be discussed later in the paper) is also reported on the right side of the figure. E, 42.07215 N, NO vertex − 410910.405 E, 4571551.350 N, SE vertex). From a lithological viewpoint, this area derives from the fragmentation of the Lazio-Abruzzi Mesozoic limestone platform. Seven pubescent oak forest stands, distributed along an altitudinal transect ranging from the lower hilly belt of the Aurunci Mountains (166 m a.s.l.) to the montane belt of the Simbruini-Ernici Mountains in the core of the Apennines (1160 m a.s.l) were selected (Table 1).
The choice of the sites of collection was made taking care these forest sites covered a sufficiently wide altitudinal and bioclimatic range. To verify whether the populations of the higher altitudes corresponded to Q. pubescens and those at lower altitudes to Q. virgiliana, we rearranged these forest stands into three main groups corresponding to different altitudinal belts and bioclimatic regions, namely MSA, COR, ACU (hilly belt; Mediterranean Region), PDA (submontane belt; Transition Mediterranean Region), MAR, TOR, SMB (montane belt; Temperate Region).

Plant materials and vegetation
A total of 209 white oak trees were sampled. In each stand, 30 mature individuals were collected distant 20 m each other Viscosi et al. 2012). Voucher specimens were deposited in the herbarium of the University of Molise (IS) (Thiers 2015). The coenological features of each forest stand were investigated using the phytosociological approach (Braun-Blanquet 1964). Seven phytosociological relevés were performed in the seven collection sites. Five of these relevés (MSA, ACU, MAR, TOR, SMB) were already published in Facioni et al. (2015) whereas the other two (COR, PDA) were still unpublished. Syntaxonomy referred to Mucina et al. (2016), while species nomenclature followed Bartolucci et al. (2018).

Leaf traits assessment
Ten fully developed and healthy leaves per individual tree were analyzed for 209 individuals for a total of 2090 leaves. Leaves were scanned by placing the abaxial surface facing upwards on an Epson GT-15000 scanner (600 dpi resolution) and subsequently were measured with the ImageJ Tool instrument (Rasband 1997(Rasband -2007. Basing on sensu Kremer et al. (2002), 14 leaf traits considered as diagnostic in the identification of Quercus species were assessed (Table S1). The pubescence of the adaxial surface PU-AD, abaxial surface PU-AB and petiole PU-PI were assessed using a grading system applied on a standard area of 1 mm 2 using a stereoscope with 35× magnification. The datasets were subjected to i) univariate statistical analysis (normality tests, Kruskal-Wallis test for independent samples, Dunn multiple comparisons, p-values computed with the Monte Carlo method and Bonferroni correction of the significance level) using XLSTAT 2020.5.1 (Addinsoft 2021) and ii) multivariate statistical analysis using PAST 4.05 (Hammer et al. 2001).

Leaf geometric morphometric analysis
The scanned images of the 2090 leaves were used to obtain a 2D landmark configuration. Basing on Proietti et al. (2021), in total 13 landmark points (LMs) were digitized on each leaf with tpsDig2 2.31 program (Savriama and Klingenberg 2011;Rohlf 2015). TpsUtil 1.78 program (Rohlf 2019) was used to obtained NTS file. Landmark coordinate configurations were superimposed by translating them to the origin, scaling to unit centroid size, and rotating, to bring all of them into a standard position and orientation, by using the generalized Procrustes analysis (GPA) in MorphoJ program (Klingenberg 2011). Shape variation was so obtained, and the resulting Procrustes aligned shape coordinates were used as shape variables in subsequent analyses. The centroid size value (CS) of each specimen was saved for analyses of size variation in the stands. A two factor Procrustes ANOVA for bilateral shape symmetry, to evaluate the relative amount of variation of asymmetry for each stand was performed. For these analyses, two classifiers were used: the leaf was selected as the individual effect and the tree was used as the extra main effect. A one-way ANOVA on differences in centroid size within and among stands (averaged by tree) was performed with XLSTAT 2021.5.1 (Addinsoft 2021). To identify and examine the patterns of leaf shape variation, a principal component analysis (PCA) on the co-variance matrix of the symmetric component produced by the 2090 leaf configurations averaged by tree (n = 209) was performed. Leaf shapes were represented with wireframes according to the corresponding scale factor of PC1 axis. Both Procrustes ANOVA and PCA were carried out in MorphoJ program (Klingenberg 2011).
The GenAlEx software v. 6.503 (Peakall and Smouse 2012) was used to perform the main basic genetic analyses: mean number of alleles (N a ), observed heterozygosity (H o ), expected heterozygosity (H e ), fixation index (G ST ) per locus and population; inbreeding coefficient (F IS ) per locus; number of private alleles (N p ) and number of different alleles (N Al ) per population. An allelic distance matrix was used to perform a 999 permutations analysis of molecular variance (AMOVA). The correlations between pairwise co-dominant genotypic distance (LinGD) and geographical distance (GGD) were tested applying simple Mantel tests with 9999 permutations. F IS per population was calculated using FSTAT 2.9.4 (Goudet 1995). The number of individuals typed (Ni), number of alleles per locus (K), null allele frequencies (F null ), polymorphic information content (PIC) and deviations from Hardy-Weinberg Equilibrium (HW) were calculated by Cervus v3.0.7 (Marshall et al. 1998).
The allelic richness per population (A r ), suitable to assess the genetic diversity (Greenbaum et al. 2014), was calculated with rarefaction to the lowest sample size using the HP-Rare software v. 6 June 2006 (Kalinowski 2005).
STRUCTURE software was used to apply a Bayesian clustering method for inferring population structure based on the genotype data (Pritchard et al. 2000). Two Bayesian analyses were performed. A first analysis included only the seven populations of the altitudinal transect. A second analysis included the latter populations together with 90 oak individuals (Q. frainetto Ten.: 30; Q. petraea (Matt.) Liebl: 30; Q. pubescens Willd.: 30) from Mount Vairano (see Figure 1). This is a mountainous massif located in the neighbour Molise administrative region whose oak forests were investigated by Antonecchia et al. (2015) using the same genetic analysis we have used in the present work. All Mount Vairano individuals exhibited a Q value >90 so that we considered them as "pure individuals" according to Lepais et al. (2009). Admixture model (Alpha α) without prior information (geographical location or taxonomical classification) was set choosing the correlated allele frequency model (Lambda λ). The degree of admixture "Alpha" was fixed to be inferred for each wood stand and the starting value was set to 1/N (where N = 10). Based on the STRUCTURE manual, the λ value has been set to be inferred for each wood stand (starting from λ = 1). A burn-in period of 50.000 and Markov chain Monte Carlo (MCMC) simulations of 100.000 was used, with possible K from one to ten and twenty replications for each value of K. The delta-K (Evanno et al. 2005) have been obtained with STRUCTURE HARVESTER (Earl and Von Holdt 2012).

Coenological and phytosociological features
The forest stands investigated exhibited different floristic and coenological features (Table S3). MSA and COR stands, both occurring in the Ausoni-Aurunci Mountains area, are reported by Facioni et al. (2015) as dominated by Q. virgiliana, with high coverages of Acer monspessulanum L. and Fraxinus ornus L. in the lower tree layer and of Ruscus aculeatus L. in the shrub layer. In this latter also several evergreen Mediterranean vines were found, such as Asparagus acutifolius L., Clematis flammula L., Rosa sempervirens L, Smilax aspera L. MSA stand exhibited a marked dominance of Carpinus orientalis Mill. in the secondary tree layer. The ACU stand being located at slightly higher altitude than MSA and COR and in the inner side of the Liri river valley it was also characterized by a significant occurrence of Ostrya carpinifolia Scop. in the upper tree layer. From a syntaxonomic point of view, MSA, COR and ACU stands were classified by Facioni et al. (2015) as Pistacio terebinthi-Quercetum pubescentis (Blasi and Di Pietro 1998) Allegrezza et al. 2002 (Table S3, columns 1-3). To note that although referring to the validly published name Pistacio terebinthi-Quercetum pubescentis, the phytosociological tables published in Facioni et al. (2015) for these three sites reported Q. virgiliana (Ten.) Ten. as dominant oak species whereas Q. pubescens Willd. was reported in the original papers of Blasi and Di Pietro (1998) and Allegrezza et al. (2002) where this association was first described. In the montane belt of Simbruini-Ernici mountains the pubescent oak stands (MAR, SMB, TOR), were located within the south-facing slopes of the inland valleys on relatively shallow soils. These wood stands exhibited a low degree of canopy cover, with Quercus pubescens dominant, and Acer opalus subsp. obtusatum (Waldst. & Kit. ex Willd.) Gams and Ostrya carpinifolia Scop. as sub-dominant species. The shrub layer exhibited several species typical of the submontane or lower montane belt (Sorbus aria Wimm. ex Nyman, Cytisophyllum sessilifolium O. Lang, Citysus spinescens C.Presl, Juniperus deltoides R.P.Adams). The herb layer was found to be characterized by Sesleria autumnalis F.W.Schultz, Teucrium chamaedrys L., Anemone apennina L. These communities were classified by Facioni et al. (2015) as Cytiso sessilifolii-Quercetum pubescentis Blasi, Avena, and Feoli 1982 (Table S3, columns 5-7) which is considered the most "continental" Q. pubescens forest type in the central Apennines (Blasi et al. 1982;Blasi et al. 2004). The PDA stand has been here interpreted as a floristically impoverished form of the Cytiso-Quercetum pubescentis (Table S3, column 4).

Leaf trait results
The Kruskal-Wallis test (Table S4) showed that the leaf traits exhibiting a discriminant power were those linked to pubescence (PU-AD, PU-AB, PU-PI) and, in less extent, leaf dimension (A, P, LL, PL, LW). The three montane belt stands (MAR, SMB, TOR) showed higher values than the hilly belt stands (MSA, COR, ACU) in both the pubescence of the two sides of the leaf blade and of the petiole. Regarding the leaf dimensional traits, the montane belt stands exhibited leaves with a slightly smaller area (A), perimeter (P), lamina length (LL), lamina width (LW) and petiole length (PL). However, a continuous trend of dimensional decrease of the leaf dimensions linked to gradually increasing altitudinal values was not found. In fact, the highest values for all the dimensional leaf-traits were found in the PDA stand in the sub-montane belt. PDA stand is also the one showing the deepest lobed leaves (highest LDR) despite displaying high SW values too. LWR values (eccentricity of the leaf contour) were found to be similar in all stands. The leaves of the montane stands were on average more obovate (OB) than those of the hilly stands except for MSA (hilly group) that showed the highest OB value among all the stands. The other leaf traits did not show significant differences among the groups. The PCA with the partition superimposed ( Figure 2) showed a relative low degree of dissimilarity among the stands. The montane belt stands segregate together in the upper left quadrant and the variables more correlated to the first PCA axis (42.66%) are those pertaining to the leaf size values. The variables correlated to the second PCA axis (16.45%) are those concerning leaf pubescence where a separation between the hilly belt stands (bottom of the diagram) from the more pubescent oaks of the montane belt stands (upper part) is observable. The polygon including the trees of the submontane belt stand (PDA) occupies an intermediate position with respect to the second axis, whereas on the first axis it is located in the right part (high dimensional values of the leaf ).

Leaf geometric morphometric results
The Procrustes ANOVA results showed a higher leaf shape variation among all the leaves of all the trees of each single stand (LEAVES) than that observed among all the leaves of each single tree per stand (TREES) (p < 0.0001). This justified the data averaging per tree used for the principal component analysis (Table S5). Conversely, directional asymmetry (DA) and fluctuating asymmetry (FA) values were not statistically significant (p > 0.0001) except for TOR stand that showed a significant value only for directional asymmetry.
The first two components of the PCA scatterplot, aimed to examine the patterns of leaf shape variation, express 63.85% of the total variance (PC1 37.49%; PC2 26.36%). Observing the first component (PC1), a significant overlapping of the leaf shapes of all the trees of the seven stands was found. Wireframe diagrams were produced to graphically display the changes in the leaf shape along the PC1 and PC2. Moving along PC1, from the negative extreme (scale factor: −0.1) to the positive one (scale factor: +0.15), some changes in shape were found. These can be summarized in longer petiole (distance between LM1 and LM2), widening of the leaf in the basal region (distance between LM6 and LM11) and in the central region (distance between LM8 and LM13) and approaching of the basal lobes to those at the largest width of the leaf blade, in the basal and central region of the leaf. This trend of changing in the leaf shape was observed in all the stands. Moving along PC2, from the negative extreme (scale factor: −0.1) to the positive one (scale factor: +0.1), a slight lengthening of the basal part of the leaf (morphologically expressed into a more wedge shape) was observed. However, any defined clusters of trees sorted by stands emerged by the PCA scatterplot. The hilly belt stands trees were found to be completely mixed with those of the montane belt stands so that no ecological or altitudinal gradient was associable to the changing of the leaf shape neither along PC1 nor PC2 (Figure 3).
(EST-SSR variation and genetic diversity within populations) The alleles frequencies analysis per locus showed a mean of 11.83 (K) and a mean of 8.12 different alleles (N a ) over all populations. The summary of all alleles recorded showed a total of 142 alleles (ΣK) (   (Table 3).
The genetic distance dendrogram ( Figure S2) exhibited very low values of bootstrap at the nodes of the main two Figure 3. Pca scatterplot based on the first two principal axes (Pcs) calculated on the symmetric component generated for averaged data (leaf configurations averaged by tree of each stand. n = 209). extreme leaf shapes along Pc1 and Pc2 are displayed in the form of wireframe graphs, namely starting (average) shape in light blue and target shape in dark blue (scale factor of −0.10 and +0.15 for Pc1; −0.10 and +0.10 for Pc2). red: acu (dot), cor (triangle), msa (square); blue: mar (square), smB (triangle), tor (dot); light green: PDa. divisions of the dendrogram (bootstrap of 10 and 26 in the first and second division respectively). No separation in genetic clusters belonging to different altitudinal belts emerged. Some hilly belt stands were found to be more closely linked to montane belt stands than to other hilly belt stands (e.g. MSA and SMB). The AMOVA (Table 4) exhibited the most relevant genetic variation occurring "within individuals" (87%), followed by the variation "among individuals within populations" (11%) and a very low variation "among populations" (2%). All F-statistics showed a significant p-value (0.001).
STRUCTURE analysis carried out on the seven wood stands established K = 5 as the most probable number of groups of genetic diversity for simulations with K number ranging between 1 and 6 ( Figure S3). The five genetic diversity groups exhibited similar weights in all the wood stands investigated. None, out of the five groups established by STRUCTURE, exhibited individuals with Q value >90 ( Figure S4). The maximum Q value (83.35) was exhibited by an individual from SMB, whereas all the other individuals of the entire dataset exhibited Q values always lower than 63.
The second STRUCTURE analysis, which, in addition to the seven populations of the ecological-altitudinal gradient, showing in some cases a slightly higher value for Q. petraea (Q ranging between 56.98 and 78.85 occurring with a single individual in the sites ACU, COR and SMB) than for Q. pubescens. The contribution of Q. frainetto was never found to be dominant (Figure 4), the maximum Q value displayed by Q. frainetto in the whole dataset being 22.44 (an individual in SMB stand with Q value for Q. pubescens = 76.30).
The cluster analysis dendrogram obtained using the genetic dataset including Mount Vairano pure populations

Discussion
The Italian peninsula represented a crucial refugial site for the forest vegetation during the cold periods of the Quaternary and oaks played a major role in the arboreal component of the flora (Bennett et al. 1991;Brewer et al. 2002). In Italy deciduous oak forests occur on all types of substrates from the coastal plain to the montane belt (up to 1400-1500 m) where they give rise to a multitude of plant communities very different from each other for floristic composition and ecology (Blasi et al. 2004). To this wide pattern of oak forest vegetation corresponds an equally wide pattern of oak trees morphological traits variability which inspired generations of nineteenth century illustrious botanists, such as Presl, Tenore, Gussone, Lojacono, Borzì etc., to describe new taxa, most of which belonging to the collective group of Q. pubescens. Some of these taxa (e.g. Quercus virgiliana, Q. dalechampii, Q. congesta), were originally considered endemic to southern Italy, but subsequently were identified also in other south-European countries. It was probably due to their acceptance in the white oak lists published in the two main monographies on the European oaks carried out by Schwarz and Camus. Namely, Schwarz (1936-39) considered the three aforementioned taxa as valid species morphologically separated and identifiable from Q. pubescens (Q. virgiliana and Q. congesta) and from Q. petraea (Q. dalechampii). Subsequently in his systematic and taxonomical treatment of European oaks Schwarz (1964Schwarz ( , 1993 made a great reduction in the number of pubescent oak species until arriving to consider only Q. pubescens as a good species and the other taxa previously accepted by him as falling within Q. pubescens overall morphological and ecological range. Camus (1938), on the other hand, considered Q. congesta, Q. virgiliana and Q. amplifolia at the rank of varieties of Q. lanuginosa Lam. (that she considered having priority over Q. pubescens), while Q. dalechampii at the rank of subspecies of it. Thus, unlike Schwarz, Camus considered Q. dalechampii to be classified among the pubescent oaks and not to be included in the collective group of Q. petraea (for a more detailed dissertation on the taxonomic and nomenclatural vicissitudes of Q. dalechampii and Q. virgiliana, refer to Di Pietro et al. (2012) and Wellstein and Spada (2015). The most recent papers concerning the Italian white oaks Di Pietro et al. 2016, 2021Proietti et al. 2021) which were carried out through a multidisciplinary approach (micro and macro-morphology, morphometry, phytosociology and genetics) revealed that the high taxonomical splitting characterizing at present the pubescent oaks collective group would not correspond to a similar and well-identifiable differentiation based on diagnostic sets of genetic or morphological characters. The pubescent oak species that more than any other is currently considered to be vicariant of Quercus pubescens in the warm and xeric areas of the Italian peninsula and major islands is Q. virgiliana (Brullo and Marcenò 1985;Brullo et al. 1999;Rosati et al. 2005;Bacchetta et al. 2009;Brullo et al. 2009;Biondi, Casavecchia, Pesaresi, et al. 2010;Blasi and Biondi 2017). In order to verify the possible occurrence of a taxonomic/ecological dualism between Quercus pubescens and Quercus virgiliana, we selected three forest sampling sites from the hilly belt of the pre-Apennines and three from the lower montane belt of the core of the Apennines. In addition, we selected an intermediate sampling site in the submontane belt of the Apennines. As expected, the pubescent oak communities were found to be clearly different in floristic and coenological terms and belonging to two different associations, namely Pistacio terebinthi-Quercetum pubescentis for the hilly belt (with Q. virgiliana as guide species of the association) and Cytiso-Quercetum pubescentis for the montane (and submontane) belts (Facioni et al. 2015). This coenological and syntaxonomic distinction, however, was found not corresponding to an equally morphological or genetic distinction observable between the two dominant pubescent oak species of the communities that would have made it possible to assign the populations of the hilly belt to Q. virgiliana and those of the montane belt to Q. pubescens. The oak specimens of the montane belt stands exhibited a lower range of variability for most of the leaf morphological traits analyzed when compared to those of the of the hilly belt stands. This result was in some way expected, being the ecological niches available for Q. pubescens s.l. in the montane belt restricted to the shallow and rocky substrates of the south facing slopes, the north facing slopes and the flat areas being prerogative of forests dominated by Quercus cerris, Carpinus betulus L., Ostrya carpinifolia, Acer opalus subsp. obtusatum and Fagus sylvatica L. The only morphological trait that partially emerged as usable to distinguish the specimens of the hilly belt from those of the montane belt is the lower degree of pubescence of leaf and petiole. However, not all taxonomists agree in considering the degree of leaf pubescence as a diagnostic trait for the identification of the species in the collective group of Q. pubescens being this trait very variable and depending on the period of collection (Gellini et al. 1992;Bussotti and Grossoni 1997;Lamant 2013). However, the lower degree of hairiness we have found in the leaves of the hilly belt stands populations (putative Q. virgiliana) somehow agree with Tenore's original description of Quercus virgiliana (Tenore 1835(Tenore -1836. In this description, that was also accepted in Flora of Italy (Pignatti 1982) the abaxial leaf blade of Q. virgiliana is described as glabrous or subglabrous. On the other hand, this description disagrees with the diagnosis provided in Brullo et al. (1999) and subsequently in the new edition of Flora of Italy (Pignatti 2017(Pignatti -2019 where the abaxial leaf blade of Q. virgiliana is described as densely pubescent. The other morphological trait that would seem to bear a certain (weak) diagnostic power is the average size of the leaves as inferable by the slightly lower values of LL, P and A traits in the montane belt populations. However, these dimensional traits, like about all the other traits measured in this study, shows a wide range of variability in the investigated populations as demonstrated by the box-plot diagrams (Table S3). Moreover, although the montane belt oak populations show averagely smaller leaves than the hilly belt oak populations, it is not possible to establish a correlation between the decreasing leaf-size and the increasing altitude since the wood stand exhibiting the largest leaves (PDA), is the one located the intermediate altitude (submontane belt) ( Table 1).
The geometric morphometric analysis which should have provided information on the variability of the leaf shape minimizing the differences linked to the leaf size, did not show significant variations passing from hilly belt stands to montane belt ones. The PCA ( Figure 3) showed a high variability in leaf shapes both within and among stands, but no discrete gradients of leaf shape related to ecological or altitudinal factors.
The level of intra-population genetic diversity among the wood stands investigated shows an average value of alleles per population of 8.1 with an average allelic richness of 6.8. Observed and expected heterozygosity values are 0.698 and 0.722, respectively. These results are slightly higher than those evidenced for the pubescent oak population of southern Calabria, Sicily and Sardinia islands  and significantly higher than those regarding pubescent oak populations of Apulia region . It is interesting to note that although in the Apulian peninsula a similar study was applied on three times as many pubescent oak populations, the genetic diversity values at loci level were found to be lower than what was found in the present study. This emphasizes the importance of being a part of a continuous belt of white oaks forests which runs throughout the whole Apennines and Pre-Apennines compared with the situation of the Apulian peninsula where forest stands are scattered and isolated in a matrix of cultivated lands, steppes-like grasslands, and Mediterranean maquis (Biondi, Casavecchia, Beccarisi, et al. 2010;Di Pietro and Misano 2010;Terzi et al. 2010). It is possible that the forest fragmentation linked to the passage from a natural landscape to an almost totally agricultural one has determined a lowering of the Apulian inter-population geneflow and a consequent depletion of genetic diversity. In fact, Chybicki et al. (2012) and Memišević Hodžić et al. (2021) highlighted on oak populations isolated or located at the boundary of their distribution area displaying lower values in allelic richness and heterozygosity. There is a large literature dealing with evaluating the levels of genetic diversity in forests in relation to their degree of spatial fragmentation and the results do not always agree. Chung et al. (2002) and Craft and Ashley (2007) observed high levels of genetic diversity in highly fragmented oak forest habitats while Kittelson et al. (2009) reported an evident decreasing of it in similar environments.
GST values are rather low when compared with the results obtained for this factor in Sicily and Sardinia and similar to those obtained for of southern Calabria (Di . It is possible that despite the ecological-bioclimatic diversity between the forest stands investigated, their geographical proximity and the widespread presence of deciduous oak forests in the areas between one site and another allowed to maintain a high rate of pollen mediated geneflow and to keep the genetic diversity among populations relatively low. Pollen immigration from neighbour areas might counteract the effects of forest fragmentation and geographic isolation (Kramer et al. 2008) even across distances of more than 50-80 km (see Buschbom et al. 2011).
The Bayesian analysis classification, aimed at identifying the number of significant genotypes and to obtain a probabilistic estimate for each individual to belong to a pure species or to be an introgressive hybridization form gave unlikely results. The most probable value of K established through the Evanno method led to the surprising identification of five genotypes (K = 5) distributed more or less evenly within the 7 populations investigated. No tree individual was found displaying Q values >90 for a given genotype. The lack of pure individuals (or approximately pure) for a given genotype does not allow to easily interpret the data in a taxonomic key. The hypothesis of two vicariant species (Q. pubescens and Q. virgiliana) to be distinguished along an ecological-altitudinal gradient clearly does not match with K = 5. On the other hand, the hypothesis of populations composed entirely of individuals deriving from a massive hybridization and subsequent multi-species introgression among Q. pubescens, Q. virgiliana eventually including Q. petraea, Q. robur and Q. frainetto is rather unlikely too, although not completely excludable, considering the high panmictic power characterizing the European white oaks. Some authors (Guarino et al. 2015;Pasta et al. 2016) suggest that the great morphological variability observable today in what we identify as Q. pubescens may be the result of the introgressive hybridization deriving from repeated crossings of sympatric Mediterranean pubescent oaks (including Q. virgiliana) whose diagnostic characters and taxonomic boundaries would gradually fade over time. Several studies however, emphasized that although the concept of biological or morphological species is challenged in oaks due to their weak interspecific boundaries and high intraspecific variation, the percentages of hybrids and introgressive forms identified in mixed populations of sympatric white oaks are always relatively low compared to those of pure individuals (Valbuena-Carabaña et al. 2005;Curtu, Gailing, Leinemann, et al. 2007, Curtu et al. 2009Gailing and Zhang 2019). As pointed out in Lazic et al. (2021), oaks are syngameons acting as genomically coherent entities (even in sympatric interfertile populations) characterized by mechanisms for maintenance of species identity in the face of regular introgressive hybridization and gene flow. Thus, if we agree that "oak species are real" as stated in Hipp (2020) after a long series of supporting considerations, it is also easier to interpret the results of our study. In fact, the Bayesian analysis response is the one expected when a clustering request is carried out on a data-set that does not contain objects classifiable in different groups. When the same Bayesian analysis is carried out on a different dataset where in addition to the original seven pubescent oaks stands also includes the already established pure populations of white oaks from Mount Vairano (i.e. Q. pubescens, Q. petraea and Q. frainetto), the most probable value of genotypes drops to K = 3. The great majority of the individuals of the seven original stands exhibit a Q value >90. Based on these results, all the seven stands investigated from the hilly altitudinal belt to the montane one, have to be ascribed to a single pubescent oak species. If we accept the taxonomic identification provided in Antonecchia et al. (2015) for the pubescent oak specimens from Mount Vairano, the name of this single species is Q. pubescens. Curtu, Gailing, Leinemann, et al. (2007) revealed that the amount of gene flow between four white oak species living in sympatry (Q. frainetto, Q. petraea, Q. pubescens, Q. robur L.), although always identifiable, is nowhere near enough to hide the genetic (and taxonomic) identity of a species. The latter statement is somewhat confirmed in this paper observing the genetic autonomy between Q. frainetto, Q. petraea and Q. pubescens as emerging from STRUCTURE analysis and confirmed by UPGMA phylogenetic tree ( Figure 5). This tree clearly displayed that Q. frainetto and Q. petraea populations segregate (bootstrap values of 92 and 100 respectively) early and independently from the other oak populations. In contrast, the pure Q. pubescens population of Mount Vairano is included in a big sub-cluster of pubescent oaks containing also the seven original pubescent-oak stands of the ecological-altitudinal transect. The very low bootstrap value (5) at the base of this pubescent oak subcluster testify that the presence of further subclusters related to the possible occurrence of more than one pubescent oak species are unlikely and statistically not significant. To note, that in the phylogenetic dendrogram of Figure 5, Quercus frainetto forms a separate cluster (at higher level) from that composed of Q. petraea and Q. pubescens. Natural hybrids of Q. frainetto are not common although hybridizations of this species with other white oaks, are reported in biosystematics literature Curtu, Gailing, Leinemann, et al. 2007;. Certainly, the greater morpho-anatomical and phylogenetic similarities existing between Q. pubescens and Q. petraea allow these two species to produce hybrids with a considerably higher frequency than with Q. frainetto (Müller 1999;Curtu, Gailing, Leinemann, et al. 2007;Salvini et al. 2009;. This should be the reason for which a greater occurrence of inter-specific geneflow Q. pubescens x Q. petraea than Q. pubescens x Q. frainetto were found in our analyses. One might wonder why this happens even in MSA and COR stands where Q. petraea is not reported either in taxonomic or phytosociological literature (Q. petraea is absent from the entire pre-Apennine sector of southern Lazio) whereas woods dominated by Q. frainetto occurs adjacently to the pubescent oak stands investigated. At present, we do not have enough data to answer this question. It can be assumed, that the ascertained presence of Q. pubescens x Q. petraea introgressive forms in areas where Q. petraea is missing is not to be considered the result of natural crossings related to the current vegetational pattern. It is more likely that they are the result of natural crossings which took place within a different forest landscape which repeatedly occurred in central Italy in previous ages when the climatic conditions were different. It is probable that up to the first post-glacial stages the presence of Q. petraea and other forest species must certainly have been greater in the Tyrrhenian side of central Italy, this latter being a well-known glacial refugial area for forest vegetation (Follieri et al. 1988;Médail and Diadema 2009;Di Pasquale et al. 2020).

Conclusion
This study analyzed from different points of view (leaf morphology and geometric morphometry, genetic and phytosociology) different populations of pubescent oaks along an altitudinal gradient. The aim was to verify if in the Tyrrhenian side of the central Apennines there were two different taxa of pubescent oaks replacing each-other within an ecological and bioclimatic gradient as reported in some floristic and phytosociological literature. The result provided a negative response. The occurrence of a pubescent oak species typical of the inner areas or the high-altitude ones (i.e. Quercus pubescens Willd.) substituted by a thermophilous oak (i.e. Q. virgiliana (Ten.) Ten.) at lower altitudes did not emerge from either the genetic or morphological-morphometric analyses. Essentially, different ecological conditions give rise to pubescent oak forest communities different from each other in the floristic composition of the undergrowth but dominated by a single pubescent oak guide species. Waiting for a multidisciplinary study aimed to analyze and clarify the whole taxonomic-nomenclatural question of the pubescent oaks collective group at European level, this species can be provisionally here classified as Q. pubescens Willd. The high values of allelic richness and heterozygosity found in the investigated populations of central Italy bodes well for their persistence and adaptability to the climatic and environmental changes that Mediterranean and sub-Mediterranean oak forests will withdraw in the next future and whose negative effects they are already experiencing at present (Conte et al. 2019). Although it deals with a very low percentage of individuals, the Q. pubescens populations analyzed were found to be genetically more similar to Q. petraea than to Q. frainetto and this result seem to be in slight contrast with the large-scale infrageneric phylogenetic classifications recently published by Hipp et al. (2020). The preponderance of introgressive forms Q. pubescens-Q. petraea than Q. pubescens-Q. frainetto, although Q. petraea is currently absent in the Tyrrhenian sector of southern Lazio whereas Q. frainetto is rather common (Scoppola et al. 1990;Blasi et al. 2002;Copiz et al. 2006), suggests a past closer coexistence of Q. pubescens and Q. petraea populations in the study area probably due to the occurrence of glacial refugees.