Abundance and Phenotypic Diversity of the Medicinal Sideritis Scardica Griseb. in Relation to Floristic Composition of Its Habitat in Northern Greece

: Sideritis scardica (S. scardica) is an endemic medicinal species of the Central Balkan Peninsula. The aerial parts are traditionally used in folk medicine and, therefore, have been collected extensively from natural habitats. Overexploitation in combination with climate change has resulted in severely fragmented populations. In this context, the purpose of this study was to access the abundance and phenotypic diversity of S. scardica populations in relation to plant community structure and environmental and anthropogenic factors in six mountainous areas of Northern Greece. For this reason, the ﬂoristic composition and diversity was determined by accessing the number of plant species, number of individuals per plant species, and plant cover in each study area. In addition, the soil properties of the studied areas were determined and the phenotypic diversity of S. scardica populations was accessed through the imaging of leaf and inﬂorescence main characteristics. As a result, 141 plant species were identiﬁed in all studied areas, while the ﬂoristic composition clearly distinguished the North-Central from the North-Eastern studied areas. S. scardica was the predominant species in the habitats where the presence of forbs was favored, while a high presence of graminoid and shrub species in the study areas depressed its presence. A high coe ﬃ cient of variations was recorded among the six populations, varying from 12.2%–29.2% and 13.3%–43.1% for inﬂorescence and leaf traits.


Introduction
Sideritis scardica (S. scardica) is an endemic species of the Central Balkan Peninsula and its distribution is restricted to Bulgaria, Greece, North Macedonia, Albania, Serbia, and Turkey [1,2]. In Greece, the species is distributed mostly to the North and North-Eastern parts of the country. The plants of S. scardica are collected from natural habitats, as they are traditionally used in medicine (for herbal preparations, infusions, decoctions, etc.) for treating gastrointestinal disorders, sore throats, asthma, bronchitis, the common cold, etc. The reported bioactivities of hydro-alcoholic and other extracts-i.e., anti-inflammatory, neuroprotective, analgesic, antioxidant, antimicrobial, as well as its gastrointestinal protection properties-are attributed to a variety of secondary metabolites of S. scardica belonging to terpenes, di-and triterpenes, phenylpropanoids, and polymethoxylated flavones [3][4][5][6][7][8].
Due to the late reports on such bioactivities, S. scardica aroused considerable interest as a medicinal plant and its natural populations have been over-exploited over the past decade. This, in combination with abiotic (environmental factors) and biotic (plant community, herbivores, etc.) factors, as well as human activities.
The aim of the present study was to access the abundance and the phenotypic diversity of S. scardica populations in relation to floristic composition, environmental factors (slope, exposure, and soil), and anthropogenic management factors (overharvesting and grazing) in six mountainous sites of Northern Greece. A thorough understanding of its habitat will highlight the present and future threats and will guide its protection through sustainable management. This, in combination with the assessment of the phenotypic diversity of the species, could potentially substitute the basis for plant selection in terms of breeding purposes, and also sustainable cultivation as an alternative of ex situ preservation.  [23][24][25]. S. scardica is a cross-pollinated species with stems up to 40 cm high, branched and/or unbranched, opposite leaves covered by gray hairs, inflorescences with dense spikes, middle bracts 12-20 mm long, lemon yellow glandular corolla, and calyx tubular-campanulate [26]. The species is reported as diploid, with 2n = 32 [27], although in some materials, chromosomal complements of 2n = 32 + 0 2B and 2n = 30 have also been discovered [28]. In natural habitats, S. scardica is mostly reproduced through seeds, while in ideal environmental conditions (i.e., temperature, humidity, soil texture, etc.), they are also reproduced through offshoot clonal propagation.

Study Area
The study was conducted in the mountainous areas of Northern Greece (prefecture of North-Central and North-Eastern Macedonia, Greece) (Figure 1), in which the study sites were identified in collaboration with local people and local forest services. Permission was asked from the Ministry of Environment and Energy in Greece so as to collect the samples legally for research activities. Six case study areas were selected, in which Sideritis scardica Griseb. (known as Greek mountain tea, Lamiaceae) is regarded as endemic; these areas were: Komnina (VerKom), Tetralofos (VerTet), Olympos (Oly), Falakro (Fal), Paggaio (Pag), and Menoikio (Men). The first two aforementioned areas are situated in Mountain Vermio, both of which belong to North-Central Greece in terms of the vegetation zone or phytosociological aspect, as well as the Olympos region [28,29]. The three final aforementioned areas belong to North-Eastern Greece.
All study areas, except for Komnina and Tetralofos, belong to Special Areas of Conservation (SAC) of the Natura 2000 network, or even to Special Protection Areas (SPA), in Northern Greece (Table 1) [30]. Moreover, all study areas are regarded as pseudo-alpine grasslands, apart from Komnina, which is regarded as a shrubland area. The habitat type 6170-"alpine and subalpine calcareous grasslands" [31] is the most dominant habitat type for Olympos, Paggaio, and Menoikio. The habitat type 62A0-"Eastern sub-Mediterranean dry grasslands (Scorzoneratalia villosae)" appeared in Falakro. S. scardica is a typical species of both habitat types. The case study areas lie at a range between 1047 m above sea level in Komnina to 2042 m in Olympos (Table 2). Initially, each area was visited, and continuous populations of the species were identified in collaboration with the local forestry services. The size of the area where these populations were identified ranged from 5 to 6 ha. Environmental factors, such as aspect and slope, were recorded, as well as plant cover, grazing intensity, and illegal collection of S. scardica. The forage utilization percentage (FUP) was estimated by the Ocular estimate-by-plot method [32] as an indicator of grazing intensity. Generally, the FUP was low to moderate, with less than 50% in all cases. There was no grazing in the areas of Komnina or Tetralofos ( Table 2). The forest service issues annual regulations for the collection of S. scardica. However, people collect it illegally for personal use and/or trade. The illegal collection of S. scardica was recorded by interviewing local people and forest service staff. Komnina and Olympos were the areas with the highest collection pressure, while there was no collection pressure in Falakro or Menoikio, which were relatively inaccessible sites and far from inhabited areas ( Table 2). As for mean annual temperature and mean annual precipitation, the data were downloaded by using the peer reviewed article of Karger et al. (2017) [33].   The case study areas lie at a range between 1047 m above sea level in Komnina to 2042 m in Olympos (Table 2). Initially, each area was visited, and continuous populations of the species were identified in collaboration with the local forestry services. The size of the area where these populations were identified ranged from 5 to 6 ha. Environmental factors, such as aspect and slope, were recorded, as well as plant cover, grazing intensity, and illegal collection of S. scardica. The forage utilization percentage (FUP) was estimated by the Ocular estimate-by-plot method [32] as an indicator of grazing intensity. Generally, the FUP was low to moderate, with less than 50% in all cases. There was no grazing in the areas of Komnina or Tetralofos ( Table 2). The forest service issues annual regulations for the collection of S. scardica. However, people collect it illegally for personal use and/or trade. The illegal collection of S. scardica was recorded by interviewing local people and forest service staff. Komnina and Olympos were the areas with the highest collection pressure, while there was no collection pressure in Falakro or Menoikio, which were relatively inaccessible sites and far from inhabited areas ( Table 2). As for mean annual temperature and mean annual precipitation, the data were downloaded by using the peer reviewed article of Karger et al. (2017) [33].

Analysis of Soil Characteristics
Six soil samples (from depth < 30 cm according to the root system of S. scardica) were taken from each study site (close to the quadrats for vegetation composition) and transferred in the laboratory. The samples were placed in trays, allowed to dry at an ambient temperature (32-35 • C), grinded in mortar, and sieved in order to obtain particles with the diameter of ≤ 2 mm. The soil pH was measured in a water-saturated soil sample (the ratio of distilled water:soil was 1:1). The soil organic matter content was recorded spectrophotometrically at 600 nm by the use of K 2 Cr 2 O 7 and H 2 SO 4 . KHC 8

H 4 O 4
Sustainability 2020, 12, 2542 6 of 16 was used as standard in order to generate a calibration curve. The content of boron (B) in the soil samples was determined by the azomethine-H method [34], using the double beam spectrophotometer HITACHI U-2900 IIO UV (Hitachi, Ltd., Tokyo, Japan). Total nitrogen (N) was determined through the Kjeldahl method with the Vapodest 50 system (Gerhardt, Königswinter, Germany). The quantitative determination of soil macronutrients was performed by an inductively coupled plasma optical emission spectrometry (ICP-OES) system (Perkin Elmer Optima 2100DV). The CaCO 3 soil was determined by the use of Bernards' apparatus after the collection and measurement of the released CO 2 , after the addition of hydrochloric acid on soil carbonates.

Diversity Data
Plant cover, the number of plant species, and the number of individuals per plant species, were recorded in order to estimate vegetation composition and floristic diversity [35,36] in each region during the flowering period in July, 2018. Six quadrats (1mX1m) were randomly placed in each of the studied regions. The specimens were collected in order to be identified in the laboratory. Species identification was performed according to Tutin et al. (1968-80;1993) [37,38], Strid and Tan (1991) [28], Strid (1986) [29], and Pignatti (1997) [39]. Nomenclature followed Euro+Med (2006) [40]. A floristic bibliography referring to these regions is included in: "Mountain Flora of Greece" [21,22] and "Flora Europaea" [37,38]. The recorded plant species were grouped into graminoids (G), legumes (L), forbs (F), and woody (W) species.
Moreover, the species number, the diversity indices of Shannon-Wiener, Simpson, Margalef, equitability, evenness, Berger-Parker, and McIntosh [41,42], as well as a Jaccard index, were estimated for each population. The Jaccard index was estimated by the following formula: where: j = the number of species common to both sites, a = the number of species in site A, and b = the number of species in site B.

Evaluation of Phenotypic Characteristics
A phenotypic evaluation was conducted in twenty-four S. scardica individuals per population (for OLY, only 16 were found) at the stage of full blooming. The maximum stem length and the number of stems per plant were recorded on an individual basis in situ for the six study sites. All stems were collected separately for each individual, transferred in plastic containers to the laboratory, and hydrated for 24 h in water. The evaluation of the basal leaf and inflorescence morphological characteristics referring to area, perimeter, length, and width was conducted through Image Pro Plus Software (Media Cybernetics).

Statistical Analysis
All data concerning plant cover, species composition, floristic diversity, and phenotypic evaluation were analyzed statistically by using a one-way ANOVA in SPSS ver. 25 for Windows [43]. For the means comparisons, a Tukey's HSD (honestly significant difference) test and standard error (SE) were used at P ≤ 0.05 to establish significant differences. In cases where the number of replicates was not equal, a Welch test was also applied in order to evaluate the significance of the ANOVA analysis.
A detrended correspondence analysis (DCA) was used for the selection of the appropriate canonical ordination. As the length of the first DCA axis-which was indicative of the beta-diversity-was higher than 3 and lower than 4 SD (3.5149), the samples seemed to be rather heterogeneous, for which simple linear methods were not suitable. For that reason, a transformation-based principal component analysis (tb-PCA), and a transformation-based canonical redundancy analysis (tb-RDA) were performed using the Hellinger transformation. In both analyses, only the taxa that were presented to all plots in each site, and they had a percentage of more than 5% in species composition were included. Prior to the analyses, all the data were logarithmically transformed. The tb-RDA was carried out for the ordination of floristic composition constrained by the explanatory variables presented in Table 2. An automatic stepwise forward model building method based on the adjusted R2 and P-value was used for the identification of significant constrains (function ordi2step). While only the significant constrains were used during the calculation of the constrained tb-RDA, the rest of the constrains were fitted on the graph in a second step (function envfit). All the analyses was carried out by using the package Vegan (v2.5-6) of R.
An additional PCA analysis was conducted in order to study the patterns of variation in the data sets of phenotypic evaluation. The PCA output was used to construct biplots to visualize the distribution (and clustering) of S. scardica populations in relation to the morphometric traits. The biplots helped in the identification of the clusters of morphometric traits that could be related mostly to each population discrimination, thus, causing its clustering.

Floristic Composition and Diversity of the Studied Areas
The floristic composition significantly differed in the studied areas. The areas of Olympos and Falakro had a significantly higher percentage of graminoids and legumes in plant cover, but a lower percentage of forbs (Table 3) compared to the other areas. On the other hand, the areas of Tetralofos, Paggaio, and Menoikio had a lower percentage of graminoids, but higher of forbs. Finally, the lower percentage of legumes were recorded in the areas of Komnina and Menoikio, while the higher were in woody species and the S. scardica were in Olympos and Menoikio, respectively (Table 3). Regarding the number of species per functional group, only the number of forbs and woody species significantly differed among the studied areas. In particular, the highest number of forbs and woody species were recorded in Tetralofos and Olympos, respectively (Table 3). In detail, 141 species in total were recorded in the studied areas (supplementary Table S1). According to the PCA, the floristic composition clearly distinguished the studied areas, with the North-Central areas (Komnina, Tetralofos, and Olympos) on the left side of the PC1 axis and the North-Eastern areas (Falakro, Paggaio, and Menoikio) on the right (Figure 2). PC1 and PC2 accounted for 58.5% and 41.5%, respectively, of the total variation. The plant species formed four distinct clusters (Figure 2). The species in the green cluster characterized the vegetation in Komnina, while the species in the blue cluster characterized the North-Eastern areas and mainly the area of Menoikio. The species in the purple cluster mainly characterized the vegetation of Olympos, but they were also broadly distributed in both North-Central and Eastern areas. Finally, the species in the pink cluster, which included S. scardica, mainly characterized the species that were predominant in both the North-Central and North-Eastern areas. In accordance with the PCA, the highest Jaccard similarity index was among the North-Central areas and the Eastern areas (Table 4). Furthermore, the lowest Jaccard similarity index (   The tb-RDA analysis revealed that the floristic composition of the studied habitats was affected significantly by the altitude and the levels of precipitation. For instance, the presence of forbs was favored in habitats characterized from high precipitation levels in comparison with graminoids and woody species. Likewise, in higher altitudes, the presence of graminoids and shrubs was favored ( Figure 3). Grazing did not seem to be a statistically significant factor; however, in the studied areas that were grazed, it was observed that the forbs were mostly developed, in contrast with the un-grazed, where the presence of woody was favored (Figure 3). The presence of S. scardica was not a significant variable for the floristic composition. Nevertheless, a high presence of grasses, such as Festuca ovina and Sesleria sp., and shrubs in the study areas where S. scardica did not exist were observed, while, in contrast, the S. scardica was dominant in the habitats where the presence of forbs was favored. Eventually, in the study sites where S. scardica was present, no severe collection phenomena were recorded.
shrubs was favored (Figure 3). Grazing did not seem to be a statistically significant factor; however, in the studied areas that were grazed, it was observed that the forbs were mostly developed, in contrast with the un-grazed, where the presence of woody was favored (Figure 3). The presence of S. scardica was not a significant variable for the floristic composition. Nevertheless, a high presence of grasses, such as Festuca ovina and Sesleria sp., and shrubs in the study areas where S. scardica did not exist were observed, while, in contrast, the S. scardica was dominant in the habitats where the presence of forbs was favored. Eventually, in the study sites where S. scardica was present, no severe collection phenomena were recorded.    Table S2) presented no significant differences among the study areas.

Evaluation of the Phenotypic Diversity among the Populations in the Studied Areas
In the present study, an attempt was made to evaluate the phenotypic diversity concerning the morphometric traits of basal leaves and inflorescences of S. scardica individuals among the six study areas, and our findings are presented in Table 6. The maximum stem length of the plants varied between 11.4 cm and 38.2 cm among the six mountainous areas. The populations from Olympos (Oly) and Falakro (Fal) represented the shorter stem length compared to the rest of the populations, where a Sustainability 2020, 12, 2542 10 of 16 2 to 3 times higher stem length was recorded. The statistical analysis revealed significant differences on the maximum stem length of the plants, with the population from Komnina (VerKom) exhibiting a significantly greater maximum stem length. In contrast with the stem length, the population from Falakro (Fal) expressed a greater number of stems per individual, differing significantly from all the other investigated populations (VerTet, Oly, Pag, and Men), except for the one from Komnina (VerKom).
In regard to the morphometric traits of inflorescence, significant variations were observed among the six populations, with the area varying between 6.1-16.2 cm 2 , the perimeter from 13.4-28.0 cm, the length from 3.5-8.7 cm, and the width from 2.8-3.0 cm. According to the statistical evaluation of the morphometric traits, the population VerKom represented superior characteristics compared to the other populations, while these were considered almost 2-fold greater in area, perimeter, and length of inflorescence from the population Oly and Fal. Although the population of Fal (together with Oly) exhibited smaller inflorescence characteristics, they represented a significantly greater leaf area (3.6 cm 2 ), perimeter (13.1 cm), and length (6.3 cm) compared to the populations VerKom, VerTet, and Oly (except the leaf perimeter). A relatively high coefficient of variations was recorded among the six populations, varying from 12.2%-29.2% and 13.3%-43.1% for inflorescence and leaf traits, respectively.
The biplot of the PCA analysis, based on the leaf and inflorescence morphometric characteristics, explained 93.45% of the total variation, calculated from Jaccard's coefficient (Figure 4). According to the biplot of PCA, the six S. scardica populations represented a clear clustering in correlation to the phenotypic traits. Specifically, the populations Fal, Pag, and Men were grouped in the same cluster, and were highly positively correlated with leaf perimeter and length. The populations Oly, VerTet, and VerKom were classified on their own in different clusters, with the latter highly correlated with the perimeter and length of inflorescence. Nevertheless, the distribution of the individuals per population, according to the PCA analysis of the phenotypic traits (supplementary Figure S1), revealed highly variable populations, such as VerKom, Oly, Pag, and Fal. In contrast, the populations of VerTet and Men seemed to be more compact as their individuals clustered together closely. The biplot of the PCA analysis, based on the leaf and inflorescence morphometric characteristics, explained 93.45% of the total variation, calculated from Jaccard's coefficient ( Figure  4). According to the biplot of PCA, the six S. scardica populations represented a clear clustering in correlation to the phenotypic traits. Specifically, the populations Fal, Pag, and Men were grouped in the same cluster, and were highly positively correlated with leaf perimeter and length. The populations Oly, VerTet, and VerKom were classified on their own in different clusters, with the latter highly correlated with the perimeter and length of inflorescence. Nevertheless, the distribution of the individuals per population, according to the PCA analysis of the phenotypic traits (supplementary Figure S1), revealed highly variable populations, such as VerKom, Oly, Pag, and Fal. In contrast, the populations of VerTet and Men seemed to be more compact as their individuals clustered together closely.

Discussion
In the present study, the floristic composition was analyzed in six mountainous areas that were habitats of S. scardica, a medicinal plant whose natural populations are endangered by overexploitation. The floristic composition clearly distinguished the North-Central areas (Komnina,

Discussion
In the present study, the floristic composition was analyzed in six mountainous areas that were habitats of S. scardica, a medicinal plant whose natural populations are endangered by over-exploitation. The floristic composition clearly distinguished the North-Central areas (Komnina, Tetralofos, and Olympos) and North-Eastern areas (Falakro, Paggaio, and Menoikio). S. scardica was in the group of plant species that were predominant in both North-Central and Eastern areas. However, its presence was very low in Komnina and Olympos in the areas with higher collection pressure. The main determinants of floristic composition in the studied sites were altitude and precipitation, while soil features did not affect it. It was the general assumption that vegetation spreading was regulated mainly by climate at regional scales [44], while soil type and biotic interactions were regulated at the site scale [45]. Similarly, Dainese et al. [46] referred that species composition in South Alpes mostly depended on temperature and precipitation.
According to the data of this survey, presence of S. scardica in the six study sites was higher in habitats where mostly forbs were composing the floristic profile, and lower in the areas where graminoids and shrubs were grown. This could be directly attributed to the competition phenomena between S. scardica and grass species, such as Festuca ovina and Sesleria sp., but also shrubs-i.e., Juniperus sp.-and for below (nutrients, water) and above ground (light) resources. In general, terrestrial plant species compete at individual basis interactions for nutrients and light, and one important mechanism for maintaining local coexistence is a trade-off in competitive ability for belowground and carbon via light [47][48][49]. However, preliminary field studies of our group revealed that although S. scardica plants have no strong demand for below ground resources, their development is highly depressed from competition, especially for light (unpublished data). As light availability is a dominant factor affecting seed set [50] through the photosynthesis assimilates, grasses or shrubs present in high density may also affect the seed-set success and maturation.
Although floristic composition seemed to directly affect the presence of S. scardica, livestock grazing could indirectly affect the species population density and communities through its effect on the relative abundance of different plant species [51,52], resulting in competition/'eutrophication' phenomena of herbaceous plants. It has been well-documented that herbivores prevent competitive exclusion through the increment of ground-level light, especially in productive systems [53,54]. Furthermore, intraspecific plant population density may substitute a significant factor that influences interactions between plant-pollinator, pollination, and plant reproductive success [55][56][57]. Hence, in habitats where herbivores modify plant communities (e.g., through browsing and trampling) this modification could indirectly affect the pollination efficiency and reproductive success of plant species that have escaped herbivores.
The studied areas were mostly grazed by free-ranging cattle. Cattle grazing starts at the beginning of spring, at lower altitude grasslands, where S. scardica does not exist. As the available forage is restricted in lowland during summer, they move to higher altitudes, where S. scardica is present. However, S. scardica is in vegetative growth and differentiates its inflorescence organs during spring and obviously receives competitiveness from the other highly dense plant species due to the absence of grazing at that time. On the other hand, grazing in subalpine Mediterranean ecosystems has been in constant decline over the last decades [58] and in some cases, such as Tetralofos in the present study, has been completely abandoned. The abandonment of grazing in these ecosystems results in shrub encroachment, the prevalence of competing species, and the reduction of diversity [58,59]. It should be noticed that this was not the case in Tetralofos, where high floristic diversity was recorded. However, the abandonment of grazing in Tetralofos is very recent and, also, there are references that the present floristic composition indicates the historical rather than the present land use [60,61]. Nevertheless, it seems that the abandonment of grazing in these ecosystems and its impact in floristic composition could be a threat for S. scardica populations.
Furthermore, it was observed that in the two least accessible study areas (Falakro and Menoikio), no anthropogenic management (collections) was reported and a higher presence of S. scardica was recorded. Such an observation was expected as, recently, an increment in illegal national collections of S. scardica was reported, which, together with false collection practices-i.e., harvesting all the plants via destructing its roots from the soil-resulted in severe fragmented populations and the unsustainable management of their natural habitats.
The phenotypic evaluation of S. scardica's populations from the six case study sites revealed a high morphological diversity of leaf and inflorescence traits, as observed from the coefficient of variations (higher than 15%). Similar to our data, Aneva and Zhelev (2019) [62] reported high phenotypic diversity (coefficient of variation (CV) higher than 20%) among six wild populations of S. scardica in Bulgaria. Several studies report the phytochemical profiling and bioactivities of Sideritis sp. herbal preparations, although there are no extensive studies on the phenotypic variation and/or plasticity phenomena of the species in correlation to the micro-climate conditions. However, such phenotypic variations seem to be common also in other Lamiaceae species-for instance, Teucrium orduini, Phlomis olivieri, and Lallemautia royleana [63][64][65]-and may also arise from the genetic and/or epigenetic diversity [20,66].

Conclusions
It could be concluded that the anthropogenic management factor, through the extensive collections of aerial parts, plays a pivotal role to the declined S. scardica populations. In addition, this study revealed the importance of generating a better knowledge on the natural habitats of endemic vulnerable species, such as S. scardica, as a variable floristic composition, and its plant community's structure directly affects the presence of S. scardica by depressing/favoring it. The interaction of S. scardica with the other functional groups of plant species identified within this study indicates that the species is not competitive, and its presence can be highly depressed from the presence of graminoids, especially in high density. The abandonment of grazing in some ecosystems has a direct impact on the floristic composition, while the absence of herbivores may indirectly affect the populations of S. scardica through favoring the eutrophication phenomena of competitive species in natural habitats. The high phenotypic variation of S. scardica populations observed in this survey revealed its putative evolutionary capacity to adapt to different environments in order to maintain their survival even under unfavored conditions. As a result of the aforementioned, the present study set the basic knowledge for the development of putative in-but also ex-situ conservation strategies for the sustainable management of the endemic species. Although this study focuses on a detailed description of the natural habitats that S. scardica is endemic, further analysis on the intra-and inter-population genetic and epigenetic diversity is performed by our group, combined with a phytochemical analysis of the six populations and a metagenomic analysis of soil samples of the case study sites. Further biological experiments on the evaluation of S. scardica's plasticity, but also the interactions and competition of the species with other plant commodities (i.e., grasses), would be of great importance.
Supplementary Materials: The following are available online at http://www.mdpi.com/2071-1050/12/6/2542/s1. Author Contributions: P.K.P. was involved in the plant species identification alongside the evaluation of floristic diversity, statistical analysis of the respective data, and drafting the manuscript. E.S. was involved in the phenotypic evaluation of the populations together with statistical analysis and interpretation of the data through PCAs, drafting, and editing the manuscript. E.A. contributed to the phenotypic evaluation and drafting the manuscript. E.M.A. was involved in the conceptualization, supervision of the present research, interpretation of the data, writing, and editing the manuscript as well. All authors have read and agreed to the published version of the manuscript.