Phenotypic Diversity of Almond-Leaved Pear ( Pyrus spinosa Forssk.) along Eastern Adriatic Coast

: Almond-leaved pear ( Pyrus spinosa Forssk., Rosaceae) is a scientiﬁcally poorly researched and often overlooked Mediterranean species. It is an insect-pollinated and animal-dispersed spiny, deciduous shrub or a small tree, with high-quality wood and edible fruits. The aim of the study was to assess the phenotypic diversity of almond-leaved pear in the eastern Adriatic region. The examination of phenotypic diversity was based on a morphometric analysis of 17 populations using ten phenotypic traits of leaves. Varieties of multivariate statistical analyses were conducted to evaluate the within- and among-population diversity. In addition, the Mantel tests were used to test the correlations between geographic, environmental, and phenotypic differences among populations. High phenotypic variability was determined both among and within the studied populations. Leaf-size-related traits proved to be the most variable ones, in contrast to more uniform leaf shape traits. Furthermore, three groups of populations were detected using multivariate statistical analyses. The ﬁrst group included trees from northern- and southernmost populations characterized by high annual precipitation. However, the trees from the second and third group were highly overlapped without a clear geographical pattern. In addition, we revealed that both environmental and geographical interactions proved to be responsible for the patterns of phenotypic variation between almond-leaved pear populations, indicating signiﬁcant isolation by environment (IBE) and isolation by distance (IBD) patterns. Overall, our results provide useful information about phenotypic diversity of almond-leaved pear populations for further conservation, breeding, and afforestation programs.


Introduction
The genus Pyrus L., family Rosaceae, belongs to the subtribe Pyrinae that corresponds to the long-recognized subfamily Maloideae in which the fruit type is generally a pome [1,2]. The genus Pyrus is believed to have originated in Central Asia, the mountainous regions our knowledge, there has been no research on the almond-leaved pear population variability to date. Scientific research on almond-leaved pear is rare [43] and mostly oriented towards its usage as a rootstock in the Mediterranean area [19,30]. However, more attention should be given to this pear species because of its valuable wood and edible fruits, as well as its great role in ecosystems, since many species of mammals feed on its fruit, such as foxes, badgers, and martens [44,45]. Furthermore, the possibility of its application in the afforestation of burned and risk areas in the Mediterranean should be examined, as it is one of the few species that tolerate repeated passages of fire (personal observations).
In the present study, the diversity of almond-leaved pear populations along the eastern Adriatic coast, one of the Mediterranean hotspots of biodiversity, was examined on the leaf material from 17 populations. Our main objectives were (1) to determine the intraand interpopulation phenotypic diversity of almond-leaved pear and (2) to test whether the patterns of phenotypic divergence across the studied populations are better explained by geographical distances (isolation by distance-IBD) or by environment differentiation (isolation by environment-IBE). We hypothesized that (1) phenotypic variability was higher within populations than among them and (2) the influence of geographical and environmental factors has a positive relationship with phenotypic divergences of populations; therefore, we expected significant isolation by distance (IBD) and by environment (IBE).

Plant Material and Studied Leaf Traits
The material for this research was collected from 17 natural populations of almondleaved pear along the eastern Adriatic coast ( Figure 1, Table S1). This area is geomorphologically unique in the world and mostly karstic, with tectonically uplifted, rocky, and steep coastlines with a pronounced anthropogenic influence throughout history [46]. Belonging to the Mediterranean region, the covered area is characterized by hot drought summer periods with a cool and wet winter period [47].
In each of the studied populations, leaf morphometric material was collected from ten trees/shrubs. Leaves were collected only from short shoots, on the external, sunlit part of the tree's crown, as those are generally considered to be the most uniform ones [48]. From each tree, approximately 30 to 50 fully developed and undamaged leaves were collected during the vegetation season in 2020. The leaves were transported to the laboratory, and then pressed and fully dried for further morphometric analysis. Finally, 20 leaf samples were randomly selected per each tree/shrub and subjected to morphometric analysis. Vouchers of the studied populations were deposited in the herbarium at the Faculty of Forestry and Wood Technology of the University of Zagreb (DEND).
A total of 3400 leaves were scanned and measured using the WinFolia program [49] with a measurement accuracy of 0.1 mm. Ten phenotypic traits were measured on each leaf: leaf area (LA); perimeter (P); form coefficient (FC); leaf length (LL); maximum leaf width (MLW); leaf length, measured from the leaf base to the point of maximum leaf width (PMLW); leaf blade width at 90% of leaf blade length (LWT); angle closed by the main leaf vein and the line defined by the leaf blade base and a point on the leaf margin, at 10% (LA1) and 25% (LA2) of leaf blade length (LA1); and petiole length (PL). Finally, 34,000 simple data values were obtained.

Environmental Data
Climate data were obtained from the WorldClim 2 database with a spatial resolution close to a square km [50,51]. Prior to the correlation analysis between morphometric, geographic, and environmental data, the correlations among all 19 WorldClim bioclimatic variables for all studied populations were calculated to exclude the highly correlated ones [52]. The following bioclimatic variables were used in the statistical analysis: BIO1 (Annual Mean Temperature); BIO3 (Isothermality (BIO2/BIO7) ×100)); BIO4 (Temperature Seasonality (standard deviation ×100)); BIO5 (Max Temperature of Warmest Month); BIO8 (Mean Temperature of Wettest Quarter); BIO9 (Mean Temperature of Driest Quarter); BIO12 Forests 2021, 12, 1630 4 of 15 (Annual Precipitation); BIO15 (Precipitation Seasonality (Coefficient of Variation)); BIO17 (Precipitation of Driest Quarter); and BIO19 (Precipitation of Coldest Quarter). Finally, ten bioclim variables, solar radiation in June (SOLAR6), and altitude were selected to describe the ecological characteristics of the studied populations, for the principal component (PC) analysis and for the calculation of environmental distances (Table S1).

Environmental Data
Climate data were obtained from the WorldClim 2 database with a spatial resolution close to a square km [50,51]. Prior to the correlation analysis between morphometric, geographic, and environmental data, the correlations among all 19 WorldClim bioclimatic variables for all studied populations were calculated to exclude the highly correlated ones [52]. The following bioclimatic variables were used in the statistical analysis: BIO1 (Annual Mean Temperature); BIO3 (Isothermality (BIO2/BIO7) ×100)); BIO4 (Temperature Seasonality (standard deviation ×100)); BIO5 (Max Temperature of Warmest Month); BIO8 (Mean Temperature of Wettest Quarter); BIO9 (Mean Temperature of Driest Quarter); BIO12 (Annual Precipitation); BIO15 (Precipitation Seasonality (Coefficient of Variation)); BIO17 (Precipitation of Driest Quarter); and BIO19 (Precipitation of Coldest Quarter). Finally, ten bioclim variables, solar radiation in June (SOLAR6), and altitude were selected to describe the ecological characteristics of the studied populations, for the principal component (PC) analysis and for the calculation of environmental distances (Table S1).

Population Diversity and Phenotypic Traits
Descriptive statistics (arithmetic mean, standard deviation, minimum and maximum value, and coefficient of variation) were calculated for the particular trait for each

Population Diversity and Phenotypic Traits
Descriptive statistics (arithmetic mean, standard deviation, minimum and maximum value, and coefficient of variation) were calculated for the particular trait for each population in order to determine the range of their variation [53,54]. To detect the level of among-and within-population variability, hierarchical analysis of variance was used. The analyzed factors were populations and trees within populations (tree factor nested inside the population factor).

Correlations between Geographic, Environmental, and Morphometric Data
The simple Mantel tests were performed in order to evaluate the correlation between multicharacter differences among populations [55,56]. Dissimilarity matrices were calculated to test correlations between geographic (latitude and longitude), environmental (ten bioclim, altitude, and solar radiation in June), and phenotypic differentiation (all studied leaf variables). Environmental and morphometric distance matrices were assessed as the Euclidian distances between the population means for the first three factors of the principal component analysis. Geographic distances were calculated as the Euclidian distance between the population sites (latitude and longitude). The significance level was assessed after 10,000 permutations, and the Mantel tests were performed with the R package "Vegan" [57].

Population Structure
Population differentiation was identified using multivariate statistical methods [58]. The K-means method was applied to detect phenotypic structure and define the number of K-groups that best explained the phenotypic variation of populations [59][60][61][62]. If the proportion of a specific population was equal to or higher than 0.7, that population was assumed to belong to one cluster, and if it was lower than 0.7, that population was considered to be of mixed origin. A dendrogram of the closest Euclidean distances on the basis of the unweighted pair-group method using arithmetic means (UPGMA) was constructed to check the structure between the studied populations. In addition, population structure was assessed using principal component (PC) analysis across all individuals and all studied traits. The input data in multivariate statistical methods were previously standardized, i.e., standardization of traits to zero mean and unit standard deviation was performed prior to each multivariate analysis. All statistical analyses were performed using the software packages STATISTICA version 13 [63] and R v.3.4.3 [64].

Population Diversity and Phenotypic Traits
The results of the conducted descriptive statistical analysis are shown for each studied population (Table S2), and for all populations together (Table 1). Of all measured leaf traits, leaf area (LA) and petiole length (PL) proved to be the most variable ones, with a coefficient of variation above 40%. Contrarily, the leaf trait with the lowest variability was the form coefficient (FC), with a coefficient of variation under 20%. In general, leaf traits related to its form (FC, LA1, and LA2) showed less variability compared to those relating to leaf size. The most variable was the southernmost population Konavle, while no population stands out as the least variable one. The population Žminj was distinguished by having the highest values for four out of the ten measured leaf traits (LA, FC, MPW, and LWT). On the other hand, population Hvar has the smallest value for the majority of measured variables (LA, P, LL, MLW, PMLW, LWT, and PL). The highest and the smallest values for traits considering leaf angles (LA1 and LA2) were noted for populations Krka and Muć, respectively. The already mentioned population Hvar, along with other southern island and coastal hinterland populations, was found to have generally smaller leaves. On the other hand, northern populations located in Istria and far southern coastal populations tend to have rather larger leaves. Statistically significant differences among and within populations were confirmed for all studied leaf traits at a significant level p < 0.001. The percent of variation was significantly higher among the trees within populations compared to the one among populations ( Figure 2). However, the component of error accounted for the greatest part of total variation for the majority of measured leaf traits. Statistically significant differences among and within populations were confirmed for all studied leaf traits at a significant level p < 0.001. The percent of variation was significantly higher among the trees within populations compared to the one among populations ( Figure 2). However, the component of error accounted for the greatest part of total variation for the majority of measured leaf traits. Figure 2. Partitioning of total variance by hierarchical level for studied leaf phenotypic traits. Acronyms for leaf morphometric traits as in Table 1.

Climate Differences among Sampling Sites
The ten bioclim variables, solar radiation in June, and altitude varied across the sampling sites ( Figure 1B). The principal components analysis revealed distinct climates based on the annual mean temperature (BIO1), seasonality of precipitation (BIO15), precipitation of driest quarter (BIO17), solar radiation in June (SOLAR6), and mean temperature of driest quarter (BIO9), along the first principal component. The second principal component was in negative correlation with the annual precipitation (BIO12) and precipitation of the coldest quarter (BIO19). The first two components accounted for 64.15% of the total variance (Table S3). These two principal components clearly distinguished the two northernmost and two southernmost populations that share high annual precipitation values from other studied populations. Acronyms for leaf morphometric traits as in Table 1.

Climate Differences among Sampling Sites
The ten bioclim variables, solar radiation in June, and altitude varied across the sampling sites ( Figure 1B). The principal components analysis revealed distinct climates based on the annual mean temperature (BIO1), seasonality of precipitation (BIO15), precipitation of driest quarter (BIO17), solar radiation in June (SOLAR6), and mean temperature of driest quarter (BIO9), along the first principal component. The second principal component was in negative correlation with the annual precipitation (BIO12) and precipitation of the coldest quarter (BIO19). The first two components accounted for 64.15% of the total variance (Table S3). These two principal components clearly distinguished the two northernmost and two southernmost populations that share high annual precipitation values from other studied populations.

Correlations between Geographic, Environmental, and Morphometric Data
In order to determine whether the observed phenotypic variability was caused by geographical (IBD) or environmental distances (IBE) between the studied populations, the Mantel tests were performed. The results identified significant correlations between the phenotypic, geographical, and environmental distance matrices ( Figure 3). Correlations were higher between phenotypic and environmental distance matrices (r = 0.4134, p = 0.002), and slightly smaller but still statistically significant between phenotypic and geographical distance matrices (r = 0.2516, p = 0.029). Therefore, leaf phenotypic traits in our populations show greater dependence on environmental than geographic factors, which explains why some geographically distant populations, such as Pula in the north and southernmost Konavle, share great similarities, despite being more than 440 km apart. phenotypic, geographical, and environmental distance matrices (Figure 3). Correlations were higher between phenotypic and environmental distance matrices (r = 0.4134, p = 0.002), and slightly smaller but still statistically significant between phenotypic and geographical distance matrices (r = 0.2516, p = 0.029). Therefore, leaf phenotypic traits in our populations show greater dependence on environmental than geographic factors, which explains why some geographically distant populations, such as Pula in the north and southernmost Konavle, share great similarities, despite being more than 440 km apart.

Population Structure
The structure of the 17 almond-leaved populations was illustrated by the K-means clustering method, which revealed the most probable division of populations into three groups ( Figure 1A To additionally analyze the relationships between the studied populations, a dendrogram was constructed using the UPGMA method of cluster analysis on the closest Euclidean distances between populations ( Figure 4). The results are compliant with the results of the K-means clustering method, and the studied populations were divided into three major clusters. The first cluster (A) that stands out the most was formed by the northernmost populations Škropeti and Žminj. The second extensive cluster (B) was formed by populations Blato na Cetini, Nin, Konavle, Pula, Slano, Sinj, Vir, and Muć. Populations Biokovo, Brač, Pelješac, Hvar, Drniš, Krka, and Obrovac formed the third cluster (C). Populations Biokovo and Brač, with almost identical altitude, were shown to be the most similar ones, followed by populations Škropeti and Žminj, which share similar altitudes and

Population Structure
The structure of the 17 almond-leaved populations was illustrated by the K-means clustering method, which revealed the most probable division of populations into three groups ( Figure 1A To additionally analyze the relationships between the studied populations, a dendrogram was constructed using the UPGMA method of cluster analysis on the closest Euclidean distances between populations ( Figure 4). The results are compliant with the results of the K-means clustering method, and the studied populations were divided into three major clusters. The first cluster (A) that stands out the most was formed by the northernmost populations Škropeti and Žminj. The second extensive cluster (B) was formed by populations Blato na Cetini, Nin, Konavle, Pula, Slano, Sinj, Vir, and Muć. Populations Biokovo, Brač, Pelješac, Hvar, Drniš, Krka, and Obrovac formed the third cluster (C). Populations Biokovo and Brač, with almost identical altitude, were shown to be the most similar ones, followed by populations Škropeti and Žminj, which share similar altitudes and almost identical precipitation values. Interestingly, populations Pula and Konavle, which are located at the opposite ends of the study area, showed great similarity, corresponding to similar environmental conditions. PC analysis showed that 90.0% of the total variation is explained by the first two principal components ( Table 2). The first principal component participates in the overall variance with 54.75% and is in highly negative correlation with the leaf size variables. In other words, individuals with larger leaves and longer petioles are grouped on the left side of the diagram, while those with smaller leaves and shorter petioles are grouped on the right side of the diagram ( Figure 5). The second principal component accounts for 35.25% of the overall variance and is negatively correlated with leaf shape variables, i.e., leaf angles (LA1 and LA2), and form coefficient (FC). A significant overlap in PC diagram was observed between the studied populations. PC analysis showed that 90.0% of the total variation is explained by the first two principal components ( Table 2). The first principal component participates in the overall variance with 54.75% and is in highly negative correlation with the leaf size variables. In other words, individuals with larger leaves and longer petioles are grouped on the left side of the diagram, while those with smaller leaves and shorter petioles are grouped on the right side of the diagram ( Figure 5). The second principal component accounts for 35.25% of the overall variance and is negatively correlated with leaf shape variables, i.e., leaf angles (LA1 and LA2), and form coefficient (FC). A significant overlap in PC diagram was observed between the studied populations.

Population Diversity and Phenotypic Traits
According to available literature, almond-leaved pear is described to have 2.5-7 cm long and 1-3 cm wide leaves with 1-2 cm long petioles [15,17,21]. Although closer to the lower limit, our results are in accordance with these descriptions, with the leaf length and width ranges between 1.3-7.3 cm (average 3.4 cm) and 0.4-3.3 cm (average 1.3 cm), and petiole length between 0.1-3.6 cm (average 1.1 cm). Comparing the morphometric leaf

Population Diversity and Phenotypic Traits
According to available literature, almond-leaved pear is described to have 2.5-7 cm long and 1-3 cm wide leaves with 1-2 cm long petioles [15,17,21]. Although closer to the lower limit, our results are in accordance with these descriptions, with the leaf length and width ranges between 1.3-7.3 cm (average 3.4 cm) and 0.4-3.3 cm (average 1.3 cm), and petiole length between 0.1-3.6 cm (average 1.1 cm). Comparing the morphometric leaf characteristics between P. spinosa and other South-East European and Western Asian narrow-leaved pear species (P. elaeagnifolia, P. nivalis and P. syriaca), almond-leaved pear stands out as the one with the smallest and narrowest leaves [9,15,21,65]. However, all of the above-mentioned species have overlapping phenotypes with regard to leaf dimensions and shape, and the similarity between them is indisputable, especially considering that all species are described as highly variable. This was also confirmed in the study conducted by Korotkova et al. [26], where P. spinosa samples were divided into two separate clades. The same authors stated that the taxon concept of P. spinosa as used in Flora Europaea [66] might not define a natural entity. Although our study was not aimed at resolving the phylogenetic relationships between narrowed-leaved pear species from southern Europe and western Asia, it is quite evident that those species are not easily differentiated because of high phenotypic plasticity and gradual trait transitions.
Our research showed leaf area and petiole length as the most variable leaf traits. Such results are in agreement with previously conducted research on other Pyrus [67], Prunus L. [68][69][70][71], and Malus Mill. species [72]. Considering that leaf area plays a major role in effective light capture, water balance, and temperature regulation [73], and the petiole has an important role in the adjustment of foliage and its inclination angles for optimal light capture [74], such results are expected. In addition, our results confirmed that the leaf shape is genetically fixed in a particular type [75,76]. The least variable leaf traits in this study were those related to the shape of the leaf lamina (FC) and leaf blade base (leaf angles LA1 and LA2).
Although almond-leaved pear forms small and isolated populations, they were highly diverse. Based on the analysis of variance, statistically significant differences were confirmed for all measured leaf traits on inter-and intrapopulation levels. In general, intrapopulation variability was higher than interpopulation variability for all of the measured leaf traits, which indicates high gene flow between populations.
Compared to the studies conducted on other fruit species, such as Prunus avium (L.) L. (CV = 10.6-38.2%) [77], Sorbus domestica L. (CV 7.2-30.7%) [78], and S. torminalis (L.) Crantz (12.7-29.3%) [48], our results exhibit a significantly higher variability of leaf traits. Such results are noticeable even compared to other Pyrus species, such as P. mamorensis Trab. (CV 15.3-30.2%) [79], P. communis L. (CV 22.3-26.5%), P. pyrifolia (Burm.f.) Nakai (3.8-15.8%), and P. syriaca (11.3-41.5%) [80]. High phenotypic variability of almond-leaved pear leaves probably resulted from a combination of multiple indicators related to the diversity of environmental factors and genotypic variability [81]. The ability of a particular genotype to exhibit different phenotypic characteristics under diverse environmental conditions is usually described as phenotypic plasticity [82,83]. Phenotypic plasticity enables species to survive in a wide range of environmental conditions and reduces the risk of species loss due to climate change [84]. Such adaptation to changing environmental conditions is crucial for a species' persistence [85] and could have a decisive role in the future global warming context [84]. Since heterogeneous species, such as pears, respond to climate change within their natural range, they are directly influenced by phenomena such as local adaptation, intra-specific diversity, and phenotypic plasticity [84], which contributes to genetic variation between populations [86]. This might explain the existing variability between our populations, as they grow under different bioclimatic conditions and therefore try to adapt in order to fully utilize available resources.

Correlations between Geographic, Environmental, and Morphometric Data
A significant relationship was found between both geographical and phenotypic distance matrices (isolation by distance-IBD) and environmental and phenotypic distance matrices (isolation by environment-IBE). The isolation by distance (IBD) model indicates that genetic differentiation between populations increases with geographical distances [87]. On the other hand, isolation by environment (IBE) explains genetic differentiation through environmental differences between populations, where phenotypic and environmental distances are positively correlated, independent of geographical distance [88]. As a consequence, populations often diverge in ways that are crucial to their interaction with the landscape, including dispersal patterns, habitat preference, and adaptation to different environmental conditions, all of which may influence patterns of gene flow. Today, ecologically driven diversification is progressively being acknowledged as the main driver of population diversification [89,90]. Our results show that the predominant pattern of phenotypic differentiation among P. spinosa populations is isolation by environment (IBE), that is, environmental effects play an important role in shaping population diversity and manifesting phenotypic differences between environmentally heterogeneous populations. Overall, our results suggest that reduction in leaf size is a likely outcome of almond-leaved pear growth in an environment with low precipitation during the vegetation season accompanied by higher temperatures and solar radiation. However, our results confirmed that phenotypic differentiation is correlated with both geographical distances and environmental differences between studied populations, that these adaptive and neutral processes are not mutually exclusive, and that natural populations are expected to experience them in combination [91,92].

Population Structure
The structure of researched populations followed IBD and IBE patterns to some extent. We obtained three clusters of populations. Interestingly, the first cluster consisted of two northernmost and partially two southernmost populations that share high annual precipitation values. The other two clusters of populations were highly intermixed, and they did not appear to follow any clear geographical or environmental pattern. Almond-leaved pear is an insect-pollinated and animal-dispersed species, which outlines the assumption that seeds are dispersed in any possible direction [93]. As there are no major geographical barriers between these populations, animals and birds that feed on almond-leaved pear fruit spread genetic material into other populations, even those on islands, causing mixed ancestry of populations. In addition, we cannot exclude the possibility that, despite the heterogeneous Mediterranean landscape and the low-density populations of P. spinosa in the study area, most of the adult trees formed part of an extensive network of pollen flow spanning discontinuous, widely spaced, open bush associations. This could possibly explain the similarity between geographically close populations. In any case, we assume that the selective pressure of the environment is strong for almond-leaved pear and likely explains the differences between the studied populations along the precipitation gradient. We suggest that the adaptability of P. spinosa to local climate affects its phenotypic traits.

Conclusions and Practical Implications
Almond-leaved pear showed great phenotypic diversity within and among natural populations distributed along the eastern Adriatic coast. The greatest diversity was found among phenotypic traits related to leaf size and water use efficiency, such as leaf area and petiole length, while mostly genetically driven leaf shape indicators were the least variable. In addition, we revealed that both geographical and environmental interactions play an important role in the patterns of phenotypic variation between almond-leaved pear populations. The predominant pattern of phenotypic differentiation among P. spinosa populations was isolation by environment (IBE), indicating that populations from similar environmental conditions exhibit similar phenotypes. This has reflected on population structure, where environmentally similar populations, as well as geographically close populations, were classified into three clusters, with relatively heterogeneous ancestry. The problematic taxonomy of this species requires comprehensive molecular and morphological research throughout its natural range. For now, this question remains open, which leaves room for future research.
The almond-leaved pear is of great importance for biodiversity in the xerophytic and thermophilic habitats of the Mediterranean area. Like other pears, the almond-leaved pear also has high quality wood, which can be used for the carving of small objects. However, on more favorable sites, it can reach bigger dimensions, so the wood can be used for furniture or as high-quality firewood. It has edible fruits that can be dried and used for tea or processed as desired. Various birds and mammals, such as martens, badgers, and foxes, feed on its fruit, which contributes to the normal functioning of the ecosystem. In addition, its currently most important use is as a rootstock for grafting various pear cultivars, especially in the Mediterranean area with calcareous soils, where drought resistance is a highly desirable trait. When collecting samples for this study, it was noticed that the species is highly fire-resistant and acts as a passive pyrophyte, indicating the potential use in afforestation of burned and other areas exposed to such risks. Consequently, the almond-leaved pear can be a potentially useful species in forestry in the Mediterranean area, especially taking into account the increasingly pronounced climate changes to which the area is highly exposed. Decreasing amounts of precipitation and increasing temperatures will force foresters to resort to new solutions in afforestation and forest management. In addition, by selection and crossbreeding, cultivars with desirable fruit traits could be obtained and could be commercially grown in such areas. It is also useful to mention its ornamental value, especially during autumn when fruits are ripening.