Genetic differentiation in functional traits among wild cherry (Prunus avium L.) half-sib lines

Understanding intra-specific variation in leaf functional traits is one of the key requirements for the evaluation of species adaptive capacity to ongoing climate change, as well as for designing long-term breeding and conservation strategies. Hence, data of 19 functional traits describing plant physiology, antioxidant properties, anatomy and morphology were determined on 1-year-old seedlings of wild cherry (Prunus avium L.) half-sib lines. The variability within and among half-sib lines, as well as the estimation of multi-trait association, were examined using analysis of variance (ANOVA) followed by Tukey's honestly significant difference test and multivariate analyses: principal component analysis (PCA), canonical discriminant analysis (CDA) and stepwise discriminant analysis (SDA). Pearson’s correlation coefficient was used to evaluate linear correlation between the study parameters. The results of the ANOVA showed the presence of statistically significant differences (P < 0.01) among half-sib lines for all study traits. The differences within half-sib lines, observed through the contribution of the examined sources of variation to the total variance (%), had higher impact on total variation in the majority of the examined traits. Pearson’s correlation analysis and PCA showed strong relationships between gas exchange in plants and leaf size and stomatal density, as well as between leaf biomass accumulation, intercellular CO2 concentration and parameters related to antioxidant capacity of plants. Likewise, the results of SDA indicate that transpiration and stomatal conductance contributed to the largest extent, to the discrimination of the wild cherry half-sib lines. In addition, PCA and CDA showed separation of the wild cherry half-sib lines along the first principal component and first canonical variable with regards to humidity of their original sites. Multiple adaptive differences between the wild cherry half-sib lines indicate high potential of the species to adapt rapidly to climate change. The existence of substantial genetic variability among the wild cherry half-sib lines highlights their potential as genetic resources for reforestation purposes and breeding programmes.

this species can be found under different climatic conditions and soil types (Welk et al. 2016). Wild cherry has a very important place in the wood processing industry, primarily for veneer production and furniture making. Likewise, it is widely used in agroforestry and as an ornamental tree in horticulture. The species has also been reported to promote biodiversity, and thus has significant ecological importance (Vuksanović et al. 2020).
Climate change has been projected to negatively affect forest ecosystems in Southeast Europe, causing changes in species bioclimatic niches and shifts in natural distribution (Stojanović et al. 2014;Pavlović et al. 2019). Wild cherry is drought sensitive due to the formation of a superficial root system. Moreover, as an early-flowering species, future increases in temperature may significantly affect flowering time, indirectly also affecting pollination biology, gene flow and population structure (Miller-Rushing et al. 2007).
The presence of substantial genetic variation at the intraspecific level may positively affect stress tolerance of populations and their adaptive potential (Pauls et al. 2013), thus contributing to faster response by the species to changing environmental conditions (Potter et al. 2017). Studying leaf functional traits variations may provide a good opportunity to evaluate intraspecific genetic diversity, as well as to understand the relationship between leaf structure and function, and the species adaptive responses to environmental factors. Moreover, such information is essential for designing breeding programmes and genetic resources conservation strategies. Application of multivariate statistical analyses for these purposes has been an efficient way in analyzing large datasets of measured phenotypic traits (Ganopoulos et al. 2018). For example, recent studies, performed on different broadleaved tree species grown in common garden experiments, showed that multivariate data analyses is a powerful tool in the assessment of a species genetic variability (Stojnić et al. 2016a;Torres-Ruiz et al. 2019;Vastag et al. 2019), as well as breeding for improved water-use efficiency .
Functional trait variations in wild cherry has been insufficiently investigated, despite their relevance for understanding the adaptation potential of a species, practical value for conservation of genetic resources, species genetic improvement and selection of suitable reproductive material. Indeed, although a number of studies have been carried out to assess the variability of wild cherry quantitative traits at different levels (e.g., between provenances, populations and half-sib lines), these researches have been restricted mainly to a single group of parameters and/or statistical methods Rakonjac et al. 2014;Popović and Kerkez 2016), thus providing limited information on multi-trait associations. Hence, the objective of this research was to evaluate and describe the diversity patterns between a broad spectrum of physiological, biochemical and morpho-anatomical characteristics of wild cherry half-sib lines, and to provide insight into the nature and magnitude of the relationship between the leaf functional traits examined. Moreover, the relationship between the half-sib lines and their adaptive response to the trial site climatic conditions, were studied in order to facilitate selection of appropriate, site-adapted reproductive material for future reforestation programs.

Plant material and experimental design
Seeds were collected in 2017 from 15 phenotypically dominant trees ranged across several sites with contrasting climatic conditions (Table 1; Fig. S1). Climate data (mean annual temperatures and annual sum of precipitation) was obtained from WorldClim 2 database for 1970 to 2000 (Fick and Hijmans 2017); altitudes were derived from Shuttle Radar Topography Mission (SRTM) digital elevation dataset (Van Zyl 2001). De Martonne Aridity index (I DM ) was calculated to evaluate the humidity of sites from which the seeds were collected, according to equation (Eq. 1): where P is the annual sum of precipitation (mm), and T stands for mean annual temperature (°C) (De Martonne 1926). Plants were grown in the nursery of the University of Novi Sad, Institute of Lowland Forestry and Environment, Serbia (45° 20′ N; 19° 51′ E; altitude 84 m). The information on climatic conditions and soil properties at the trial site was provided by Krstić et al. (2014). In August 2018, when seedlings were five months old (first growing season), gas exchange measurements and sampling of plant material for laboratory analyses (i.e., biochemical, anatomical and morphological analyses) were performed in 10 replicates per each half-sib line ( Table 2). All the leaves used for measurements/sampling were fully developed, with no symptoms of damages caused by pests or diseases.

Physiological analyses
ADC Bioscientific Ltd. LCpro + portable photosynthesis measuring system was used for measurements of net CO 2 assimilation rate (A), transpiration (E), stomatal conductance (gs), and intercellular CO 2 concentration (Ci). The irradiance and CO 2 concentration inside the leaf chamber were set on 1000 μmol m −2 s −1 and 350 μmol mol −1 , respectively, while the air temperature and air humidity were taken from the atmosphere. Intrinsic water-use efficiency (WUE) was (1) I DM = P∕ (T + 10) calculated as the ratio of net CO 2 assimilation rate and stomatal conductance (Hatfield and Dold 2019) (Eq. 2): (2) WUE = A ∕ gs ( mol mol −1 )

Biochemical analyses
Harvested leaf samples of different half-sib lines were frozen in liquid nitrogen, lyophilized and ground to a powder with a pestle and mortar. A 20 mg sample of lyophilized leaves was mixed and vortexed with 2 mL of 80% ethanol, and Petiole length PL mm the supernatant obtained after centrifugation at 10000 rpm for 20 min at 4°C was used for the determination of total antioxidant properties (ABTS test), ferric reducing power (FRAP) and phenolic content. Determination of total soluble protein was carried out using phosphate buffer extracts prepared by maceration of 0.1 g of fresh leaf with 2 mL of 0.1 M phosphate buffer (pH = 7.4), and the supernatant used for protein determination. Total soluble protein content (TSP) was measured according to Bradford (1976) using bovine serum albumin (BSA) as a protein standard. The protein concentration in the sample was expressed in mg BSA per gram of fresh weigh.
The amount of total phenolic content (TPC) was determined using the method of Kim et al. (2003). A 25 μL volume of the ethanol extracts, 125 μL of 0.1 M Folin-Ciocalteu reagent and 100 μL of sodium carbonate (7.5% Na 2 CO 3 ) were mixed and incubated for 5 min. Absorbance of the resulting dark blue colored liquid was measured at 760 nm using a spectrophotometric plate reader (MultiScan GO, Thermo Fisher Scientific, Waltham, MA, USA). The results were expressed as gallic acid equivalents (GAE) in milligrams per g of fresh weight of the initial sample.
Antioxidant activity of ethanolic extracts was estimated in terms of the ABTS + radical scavenging capacity following the procedure described by Miller and Rice-Evans (1997). A 30 μL sample was diluted to 270 μL of ABTS + solution in microplate wells. The absorbance was measured at 734 nm on a microplate reader (Thermo Fisher Scientific, Waltham, MA, USA). The radical scavenging capacity (RSC) was expressed as a percentage using the following equation (Eq. 3): where A control is the absorbance of the reagent, A sample the absorbance of the reagent in presence of the extract.
The FRAP test measures the ability of antioxidants to reduce the ferric 2.4.6-tripyridyl-s-triazine complex [Fe 3+ -(TPTZ) 2 ] 3− to the intensive blue-colored ferrous complex (Fe 2+ -(TPTZ) 2 ] 2− in an acidic medium. For successful reduction measurement, 20 µL of previously made plant ethanol extract was added to the mixture of 225 µL of FRAP reagent and 25 µL of water and shaken for 20 s, and the absorbance was recorded at 593 nm on a microplate reader (Thermo Fisher Scientific, model Multiscan GO, USA). Ascorbic acid solutions, with concentrations ranging up to 500 μmol, were used as standards for construction of the calibration curve. The results were expressed as ascorbic equivalents (AE) in milligrams per gram of lyophilized weight of the initial sample (Benzie and Strain 1996).

Leaf anatomical and morphological analyses
Stomatal density per mm 2 (SD), stomatal guard cell length (LS) and width (WS), and stomatal aperture length (SAL) and width (SAW) were determined from leaf prints made following the protocols of Kardel et al. (2010) and Stojnić et al. (2015a).
Morphological measurements were taken on the leaves from which stomatal imprints had been previously taken. Leaf length (LL), width (LW) and petiole length (PL) were measured using a scale ruler. Leaf area (LA) was determined with the ADC Bioscientific Ltd. AM300 Portable Leaf Area Meter on the leaves stripped of petioles. Afterwards, the same leaves were dried at 70°C for 72 h and dry mass (DM) determined. Leaf dry mass per unit area (LMA, mg cm −2 ) was calculated using the formula after De la Riva et al.

Statistical analyses
All statistical analyses were performed by STATISTICA software, version 13 (TIBCO Software Inc., 2017). For each dependent variable, Shapiro-Wilk's and Levene's tests were performed to analyze normality of distribution and homogeneity of variances, respectively (Sánchez-Gómez et al. 2013). Likewise, parameters expressed in percentages (ABTS) and parameters describing the number of stomata per unit leaf area (SD) were transformed by arcsin (arcsin √ x ) and square-root transformations ( √ x + 1 ), respectively, to meet normal distribution of frequencies and assumptions required for application of parametric statistics (Kovačević et al. 2020). The variability of the study parameters of the wild cherry half-sib lines was described by one-way analysis of variance (ANOVA) followed by Tukey's HSD (honestly significant difference) test, at a significance level of P = 0.05. Based on the outputs from ANOVA, the contribution of sources of variation to the total expected variance of measured parameters was evaluated. Coefficients of variation (CV %) were also determined as indicators of variability. The relationship between all parameters was assessed by Pearson's correlation coefficient. The significance of correlations was tested at P = 0.05.
Principal component analysis (PCA) and canonical discriminant analysis (CDA) were performed to estimate and summarize the relationships among the leaf phenotypic traits, as well as to examine patterns of variation between the half-sib lines. Means of the observed phenotypic traits at the level of genotype were used as entering data for PCA and CDA, by which principal components and canonical variables (roots) were extracted, respectively. Similarly, for the analysis of relationships between parameters, the loadings were calculated with principal components and canonical variables that had eigenvalues higher than 1 . In addition, the selected principal components were rotated by the varimax method to maximize variance of loadings within the same principal component. Parameters that had the highest loadings with a certain principal component or canonical variable are considered to be closely related and therefore assembled in the same group.
Within the procedure of the forward stepwise discriminant analysis (SDA), subsequent discriminant models were formed by the inclusion of those parameters that contributed to the discriminating ability of the model to the largest extend. The order in which the parameters were added, together with the coefficient of correct allocation of analogue subsequent models, were used to discuss the importance of the selected parameters in the discrimination of the half-sib lines (Kovačević 2014).

Results
The ANOVA and Tukey's HSD test for the measured quantitative traits of wild cherry half-sib lines are presented in Tables S1 and S2, respectively. Highly significant differences (P < 0.01) were observed between the half-sib lines for all traits. The highest F-values were recorded for gas exchange parameters, particularly for E (41.5) and Ci (21.3), as well as for TPC (17.6). The highest coefficients of variation between half-sib lines (CV B ) were observed for E, gs, WUE and LA (CV B > 20%), whereas the lowest CV B were found for leaf stomatal traits and LMA (CV B < 8%). In addition, Tukey's HSD test showed that characteristics E and gs had the greatest number of homologous groups (eight), thereby confirming their high variability, recorded previously by the F-test and coefficient of variation. The smallest number of homologous groups were found for traits related to stomatal guard cells and aperture sizes (SL, SW, SAL and SAW) (Table S1). Therefore, according to the F-test and coefficient of variation, E and gs were the traits by which the wild cherry half-sib lines differed the most (Table S2).
Besides the F-values and coefficients of variation, the variation patterns across the traits were analyzed by the contribution of within and between half-sib lines expected variation to the total expected variance (Fig. 1). The results showed that the greatest contribution to variability between the half-sib lines was for E, Ci, gs, WUE and TPC (more than 60%). Nevertheless, variability at the level of the halfsib lines was higher than between half-sib lines for the most of traits (12 out of 19). The highest impact of this source of variation was found for traits related to stomatal aperture size (SAL and SAW) and TSP (more than 80%).
The results of the Pearson's correlation analysis for wild cherry quantitative traits are in Table S3. Many of the relationships appeared to be significant, although they do not necessarily show high correlation coefficients (i.e., a correlation coefficient between 0.4 and 0.5 might be considered as Fig. 1 Contribution of the sources of variation to the total expected variance for parameters of the wild cherry half-sib lines. Abbreviations: Ci (intercellular CO 2 concentration), E (transpiration), gs (stomatal conductance), A (net photosynthesis), WUE (intrinsic wateruse efficiency), TPC (total phenolic content), FRAP (ferric reducing antioxidant power), TSP (total soluble protein content), ABTS (radi-cal scavenger capacity against ABTS), SD (stomatal density), SL (stomatal guard cell length), SW (stomatal guard cell width), SAL (stomatal aperture length), SAW (stomatal aperture width), LA (leaf area), LMA (leaf dry mass per unit of leaf area), LL (leaf length), LW (leaf width), PL (petiole length) a moderate relationship between the two traits). However, at the probability level of P < 0.05, statistically significant linear associations were revealed in approximately one fourth of the cases (40), suggesting a moderate multicollinearity within the set of parameters. Significant correlations were observed between E, gs and A, as well as between these parameters and LL and SD, emphasizing the considerable association of gas exchange with leaf size and stomatal density (Table S3). Likewise, SD was significantly correlated with the parameters of leaf size, such as: LA (0.59) and LL (0.73), and moderately with LW (0.47). In addition, Ci was significantly correlated with the most of the biochemical traits (ABTS, TPC and TSP), as well as with LMA.
PCA revealed which traits played an important role in phenotypic variation between the half-sib lines. According to the PCA, the first three principal components (PCs) accounted for 37.5%, 20.4% and 14.3% of the diversity between the half-sib lines, respectively (72.1% of the total variance) (Table S4). To analyze the relationship between the traits, the first four PCs were selected according to Keiser's rule (i.e., eigenvalue > 1) to be rotated by the varimax rotation. The quantitative traits that had the highest loading with certain PCs were assigned to the same group.
Considering the first four PCs, high values of maximum loadings (over 0.6) were detected after the rotation for almost all the traits, thus confirming the relationship between the traits based on the correlation matrix. The highest contribution on PC1 corresponded to parameters of gas exchange (E, gs, A, WUE), lamina size (LA, LL, LW), SD and FRAP values. Variables responsible for the differentiation along PC2 were Ci, LMA, and almost all biochemical parameters (ABTS, TPC and TSP), indicating a strong relationship between biomass accumulation and intercellular CO 2 levels with antioxidant capacity of the plants. Characters contributing to the differentiation along PC3 were those describing the stomata (SL, SW and SAL) and PL, confirming the negative association of petiole length with stomatal size. PC4 explained 8.7% of the total variance, and only the SAW had its highest loading with this principal component (Table 3).
Information on the relationship between the half-sib lines is presented at the 3D scatter plot that shows a genetic similarity in terms of the quantitative traits (Fig. 2). By starting from the negative to the positive values of PC1, wild cherry half-sib lines showed a gradual increase in parameters closely related to photosynthesis (A and gs), leaf blade size, and stomatal density, as well as a general decrease in E, WUE and FRAP. With regards to PC2, going from the negative towards positive values of the axis, the half-sib lines were characterized by higher values of LMA, TPC, TSP and ABTS, as well as lower Ci. Furthermore, the results showed that two groups related to the half-sib lines habitat humidity could be separated along the PC1. The first group, situated on the negative side of PC1, included half-sib lines originating from "humid" sites (28.0 ≤ I DM < 35.0; Baltas 2007) of Northern Serbia (CO1, CO2, JA11, JA13, JA14 and MO7), whereas the rest of the half-sib lines (except for LR35) originated from "very humid" sites (35.0 ≤ I DM ≤ 55.0; Baltas 2007) of Central, Western and Eastern Serbia were grouped on the positive side of PC1. In addition, the highest distance at PC1 was observed between the half-sib lines from the Bačevci locality (site with the highest value of I DM ) on the positive side, and Jamena and Čortanovci localities on the opposite side (sites with the lowest I DM values). Along the opposite sides of PC2, a grouping of the half-sib lines originating from Western Serbia (Bačevci site) and three half-sib lines from Central Serbia (Lomnička Reka), respectively, could be noticed.
A multivariate test to assess the influence of the traits under examination to the discrimination of the half-sib Table 3 Loadings between the parameters of the wild cherry halfsib lines and the first four principal components after their rotation by varimax method Bolded values represent the highest loading of the particular trait Abbreviations: Ci (intercellular CO2 concentration), E (transpiration), gs (stomatal conductance), A (net photosynthesis), WUE (intrinsic water-use efficiency), TPC (total phenolic content), FRAP (ferric reducing antioxidant power), TSP (total soluble protein content), ABTS (radical scavenger capacity against ABTS), SD (stomatal density), SL (stomatal guard cell length), SW (stomatal guard cell width), SAL (stomatal aperture length), SAW (stomatal aperture width), LA (leaf area), LMA (leaf dry mass per unit of leaf area), LL (leaf length), LW (leaf width), PL (petiole length) lines, showed a statistically significant contribution for the majority of characters except for certain morphological (LA, LL, LW and LMA) and biochemical (ABTS and TSP) parameters (Table S5). CDA was performed as a part of the analysis at the evaluation of the relationship between traits. Results of χ2 test statistics showed that the first nine canonical variables (roots) were significantly related to the variations between wild cherry half-sib lines. The first five had eigenvalues greater than 1 and were characterized by high canonical coefficients and low Wilks' lambda (Table S6). However, within the first three canonical variables that described more than 80% of the total variation, only E had its highest loading with the first canonical variable, no trait had its highest loading with the second canonical variable, and only TPC had its highest loading with the third canonical variable (Table S7). The relationship between the halfsib lines based on the first three canonical variables is shown in Fig. 3, and shows a similar grouping of half-sib lines as observed by PCA.
The forward SDA was performed in addition to determine the importance of the phenotypic traits contribution to discrimination of the half-sib lines. Fifteen traits characterized by significant F statistics were included in the discriminant model by the stepwise procedure (Table S8). The first parameter was E and the model with this trait alone achieved 40% of correct allocation, while the model with all parameters involved achieved 96% of correct allocation (Fig. 4).
The model based on the first six traits (E, gs, TPC, FRAP, Ci and WUE), included by SDA, was sufficient to achieve 77% of correct allocation. This result largely corresponds to the outputs from the analysis of variance components, since the contribution of expected variance between halfsib lines to the total variance was higher than 60% for all the traits except FRAP. The same traits had also the highest loadings with rotated PC1 and PC2, therefore suggesting a small contribution of traits from other groups (i.e., PC3 and PC4) to the differentiation of the half-sib lines. Finally, the traits E and TPC, which had their highest loadings with the first and third canonical variables, respectively, were among the first three variables to be included in the model by the forward stepwise procedure.

Variability of the half-sib lines
Variability between the half-sib lines was statistically significant for all traits, indicating a high genetic diversity of the wild cherry progeny, which has been already confirmed for certain biometrical (Petrokas and Pliûra 2014;Stanković-Neđić et al. 2018), and morphological characteristics (Temel 2018;Miljković et al. 2019). Although significant differences between the lines were observed for all parameters, the analysis of contribution of within and between the halfsib lines expected variances to the total expected variance showed that E, gs, Ci, WUE and TPC had a significant effect on divergence between the half-sibs. These results have been confirmed by SDA which showed that the same parameters, including FRAP, were the first six traits added to the model, thus indicating the importance of characters defining gas exchange and antioxidant capacity in distinction of the halfsibs. The differences with regards to different phenotypic traits are supposed to be the result of their different genetic constitutions, which develop as a consequence of the species evolutionary adaptability to divergent selective pressure across the distribution range (Stojnić et al. 2015b). Indeed, our study involved wild cherry half-sib lines originating from diverse habitats across Serbia; e.g., spanning from drier northern sites to the more humid, mountain sites of Central and Southern Serbia. The existence of substantial genetic diversity at the intraspecific level, as observed by our study, may positively affect the stress tolerance of wild cherry populations and their adaptive potential under changing climate (Pauls et al. 2013;Potter et al. 2017). Nevertheless, even though high genetic differences were observed between the half-sibs, the results show that most of the variation was harbored within the half-sib lines, i.e., variability at the level of half-sib lines had a bigger impact on total variation in the majority of the traits. These results are in agreement with a number of studies dealing with both quantitative traits (Popović and Kerkez 2016;Lobo et al. 2018;Miljković et al. 2019) and gene markers (Ganopoulos et al. 2011;Jarni et al. 2012) in wild cherry. For example, studying the preference of Myzus cerasi (Fabricius 1775) to half-sib lines of wild cherry from several populations, Poljaković-Pajnik et al. (2019) reported that the contribution of population variation was only 1.5% to the total expected variance of the damaged leaves percentage, whereas the contributions of half-sib lines within populations and genotypes within half-sib lines were significantly higher, 45.6% and 52.9%, respectively. High contribution from within half-sib line variance is possibly related to the mating system of wild cherry. The species is characterized by gametophytic self-incompatibility which is genetically controlled by multi-allelic S-locus, thus promoting out-crossing and preventing self-fertilization (Schueler et al. 2003;Vaughan et al. 2007). In other words, when the S-allele of the haploid genotype of the pollen grain matches one of the S-alleles of the diploid genotype of the stigmatic tissue, fertilization cannot occur, as the pollen tube development in the style is inhibited (Ganopoulos et al. 2012;De Rogatis et al. 2013). Nowadays, several researchers have shown evidence of a broad multiplicity of the S-locus in wild cherry natural populations; e.g., Schueler et al. (2006) have reported 15 S-alleles within a single population; Ganopoulos et al. (2012) found 16 S-alleles in five Greek populations; and Sharma et al. (2017) identified 18 different S-haplotypes across wild cherry populations in the Czech Republic. As the results of this mechanism, high genetic diversity is commonly observed within populations, whereas the differentiation between populations usually remains low to moderate (De Rogatis et al. 2013). Moreover, as entomophilous species, the pollen of wild cherry is spread by the insects often across large distances (Jolivet et al. 2012), whereas the seeds are efficiently dispersed by birds and mammals, thereby promoting an extensive gene flow.
Yet, PCA and CDA determined the pattern of variation among the studied half-sib lines. The results of PCA showed that they were grouped along PC1 according to the humidity at their original locations, indicating genetic variation among half-sibs induced by the climate of origin sites. The first group, situated on the negative side of PC1, included half-sib lines originating from Northern Serbia (CO1, CO2, JA11, JA13, JA14 and MO7), which might be considered "humid", according to the De Martonne climate classification (Baltas 2007). The rest of the half-sib lines, grouped on the positive side of PC1 (except LR35), are geographically separated and originate from regions south of the Sava and Danube Rivers which are characterized as "very humid" according to De Martonne Aridity Index (Baltas 2007). The same grouping was also observed using CDA.
Our results correspond to findings of other researchers who reported decreased physiological performance at the leaf level in plants originating from drier sites and grown under different watering regimes. The observed differences in physiological and morphological traits between two groups of half-sib lines may explain specific adaptation strategies to the environmental conditions of their sites of origin, indicating enhanced capacity of seedlings from drier sites to overcome prolonged water stress. For example, Zhang et al. (2004) demonstrated that dry climate Populus davidiana Dode. populations produced lower net photosynthesis and transpiration, and smaller leaf area under both watered and drought treatments, in comparison to populations originating from humid sites. Reduced gas exchange in plants from habitats with lower precipitation and higher air temperatures represent evolutionary mechanisms to conserve water usage under limited soil water availability. In addition, half-sib lines originating from drier sites of Northern Serbia were generally characterized by smaller leaf length and width, as well as smaller lamina size. Such characteristics are associated with a plant's water-saving strategy to reduce transpiration, as well as effective heat dissipation (Yan et al. 2016). Similarly, smaller LA is linked to a higher major vein density that might positively affect both the resistance of the leaf to hydraulic conductivity losses and specific conductivity (Scoffoni et al. 2011;Peguero-Pina et al. 2014).
Furthermore, the results show that SD was approximately 47% lower in the half-sib lines from drier sites, corresponding to higher FRAP values. In contrast to our results, increased SD in plants from drier sites has been frequently reported as an adaptation mechanism to minimize water loss (Pearce et al. 2006;Sun et al. 2014). However, some studies indicated high flexibility of SD in response to water status (Xu and Zhou 2008). Bussotti et al. (2005) concluded that increased SD in northern Fagus sylvatica L. stands in Italy, in comparison to southern stands, represented the adaptation to the physiological risk of drought stress. Moreover, studies have shown that drought tolerance of various crop species can be improved by reducing stomatal density (Hughes et al. 2017;Caine et al. 2019).
Previous studies have reported increased antioxidant activity in plants exposed to soil water deficits (Stojnić et al. 2016b;Sarker and Oba 2018), as well as populations adapted to water-deficit conditions (Visi-Rajczi et al. 2021). Accordingly, higher FRAP recorded by the wild cherry halfsib lines from localities with lower I DM possibly represent genetically-driven responses to drier climatic conditions of the sites of origin (Kučerová et al. 2018).

Relationship between leaf functional traits
The results of the PCA demonstrated a strong relationship between leaf structure and function, as revealed by the correlations among parameters with the highest loadings at PC1 and PC2. The highest contribution on PC1 corresponded to variables related to gas exchange (except Ci, which had the highest loading at PC2), leaf blade size, stomatal density and FRAP, an important indicator of plant antioxidant ability. By regulating water and CO 2 fluxes between leaves and the atmosphere, stomatal conductance has been recognized as the main driver of water relations and photosynthesis in plants (Urban et al. 2017). Therefore, there was a strong positive correlation between A and gs, while E was negatively correlated with both. Although A and E are regulated by gs and they mutually affect each other (Urban et al. 2017), the relationship between these parameters is not straightforward but depends on number of factors affecting stomatal aperture, subsequently modulating the contribution of stomatal regulation to photosynthesis (Kusumi et al. 2012;Stojnić et al. 2016b). Therefore, very often contrasting results have been reported on the relationship between these traits. Bojović et al. (2017), for example, showed that nature of the relationships between A and gs, as well as between E and gs were notably different in Quercus robur L. seedlings grown under drought and well-watered conditions. Based on Pearson's correlation analysis, a significant positive correlation was found between gs and SL (r = 0.53), confirming the results of Dow et al. (2014) who found evidence that stomatal conductance is affected by stomatal aperture characteristics. Similar results were also reported by Stojnić et al. (2019) for Q. robur half-sib progenies. Likewise, PCA showed a positive relationship between gs and SD, which corresponds to the results for different plant species (Pearce et al. 2006;Xu and Zhou 2008). In addition, there was a significant negative correlation between WUE on one side, and SD, gs and A, on another, suggesting that plants with lower stomatal density tend to have higher water-use efficiencies, as well as reduced CO 2 uptake, and consequently, lower assimilation rates (Drake et al. 2013). Our results are also in agreement with the findings of Stojnić et al. (2019) who demonstrated that variation of WUE in Q. robur half-sib lines was strongly dependent on both SD and FRAP.
The results in this study show a positive correlation between leaf size and gas exchange parameters, confirming the importance of a large surface area for the light interception and access to atmospheric CO 2 , essential for the process of photosynthesis (Perez et al. 2019). Several studies have highlighted the importance of plant architectural characteristics, including leaf angle, length and shape with respect to light absorption capacity and photosynthesis (Vos et al. 2010;Sarlikioti et al. 2011).
Regarding the traits with the highest loadings at PC2, there was a negative correlation between internal CO 2 concentration and plant antioxidant capacity, which is possibly related to the ability of plants to withstand oxidative and osmotic stresses that occur in plant cells simultaneously with drought stress. Under normal conditions, reactive oxygen species (ROS) are present in plant cells in low and controlled levels, as side products of oxidative respiration (Dumont and Rivoal 2019). However, in response to drought, plants produce abscisic acid (ABA) which triggers stomatal closure to reduce water loss and preserve hydraulic conductivity, thereby also decreasing the CO 2 available for photosynthesis (Sánchez-Rodríguez et al. 2010). Low internal CO 2 concentrations prevent the regeneration of the main electron acceptor NADP + , causing further imbalance between the generation of ROS and the use of electrons (Petridis et al. 2012). As the absorbed light energy cannot be completely used for photosynthesis, the misdirecting electron transport activities of chloroplasts initiate inevitable leakage of electrons to molecular oxygen, thereby overproducing ROS (Sharma et al. 2012). To avoid ROS-induced oxidative damage and to increase their antioxidant capacity, plants develop efficient enzymatic and non-enzymatic scavenging systems aimed to suppress negative effects towards oxidative injuries of the photosystem II (Laxa et al. 2019). For example, recent studies showed significant increases of ROS scavenging capacity evaluated by the ABTS test in different plant species exposed to soil water deficits and heavy metal contamination (Kebert et al. 2016;Klunklin and Savage 2017;Tani et al. 2019). Likewise, phenolic compounds are reported to have strong antioxidant effects, and commonly accumulate in plants subjected to water stress to prevent ROS formation (Sánchez-Rodríguez et al. 2010). Moreover, the increase of soluble proteins in leaves and branches has been described as part of the plants' strategies aimed at adjustment of the osmotic potential of cells and re-establishment of cell turgor, ensuring in that way the normal functioning of physiological processes under drought stress conditions (Li et al. 2010).
The results further demonstrate that LMA, although having the highest correlation with PC2 (0.625), had also moderate loading on PC1 (− 0.413). Likewise, Pearson's correlation analysis revealed moderate, but statistically significant correlation of LMA with E (r = 0.56) and WUE (r = 0.54). These results are in accordance with the previous findings that confirmed a strong relationship between LMA and plant photosynthetic and transpiration rates (Novriyanti et al. 2012;De la Riva et al. 2016), as well as LMA and the activity of certain antioxidant enzymes (Dai et al. 2017). Similarly, as the mesophyll structure greatly affects water and CO 2 diffusions inside leaves (Novriyanti et al. 2012;Muir et al. 2017), and consequently CO 2 assimilation rate and stomatal conductance (Quero et al. 2006;Hassiotou et al. 2010), several studies have shown significant correlation between LMA and WUE in different plant species (Park et al. 2016;Stojnić et al. 2019).

Conclusions
The existence of substantial genetic variability among and within wild cherry half-sib lines is a key requirement for this species rapid and effective genetic adaptation to ongoing climate change. However, in order to maintain high genetic adaptability, the measures of in situ and ex situ conservation have to be considered at the national level, also taking into account recommendations of programmes of international relevance, such as the European Forest Genetic Resources Programme (EUFORGEN).
Likewise, high genetic diversity of the wild cherry halfsib lines highlights their great potential as genetic resources for breeding programmes oriented toward increased resistance to various biotic and abiotic stress factors. Analysis of the contribution of the sources of variation to the total expected variance and multivariate statistical analyses (e.g., PCA and CDA) showed that half-sib lines were distinguished primarily on the basis of several leaf functional traits related to gas exchange capacity, leaf blade size, stomatal density and ferric reducing antioxidant power values. Additionally, the forward stepwise discriminant analysis identified traits (transpiration rate, stomatal conductance, total phenolic content, ferric reducing antioxidant power, intercellular CO 2 concentration and water-use efficiency) that affected these distinguishing features to the largest extend, indicating detrimental values of these parameters as reliable markers for the characterization of the wild cherry half-sib progeny.
Moreover, the success of reforestation, as well as longterm forest stability largely depend on the use of appropriate, habitat-adapted reproductive material. We found that halfsib lines originating from sites with lower annual precipitation and higher mean annual temperatures were characterized with adaptation mechanisms that prevented large water losses. Therefore, considering scenarios of climate change in South-east Europe, assisted migration should be considered in future reforestation programmes where these half-sib lines might be used to complement local seed sources in order to enhance their fitness.