A Geometric Morphometrics Approach to the Study of Natural Variations and Hybrid Detection in Populations of Alnus incana (L.) Moench and Alnus rohlenae Vít, Douda and Mandák

Landmark-based geometric morphometrics (GM) was used to examine, for the first time, spontaneous hybridization between Alnus incana (L.) Moench and Alnus rohlenae Vít, Douda and Mandák, and to assess inter- and intrapopulation variability in leaf shape, leaf size and venation in natural populations in Serbia (Western Balkans). Two geographically distant (30 km) and two close (1.2 km) populations were selected to examine hybridization. The variability in leaf shapes was assessed by canonical variate analysis and linear discriminant analysis performed on the symmetric component of variation. Covariation between the symmetric component of shape variation and the number of pairs of secondary leaf veins was investigated with partial least squares analysis. Static allometry was examined for the first time in the genus Alnus Mill. A higher proportion of A. incana leaves was classified as A. rohlenae in geographically close populations, which is in accordance with the hypothesis about spontaneous hybridization. No single leaf of A. rohlenae was classified as A. incana, indicating that putative hybrids can only be found in grey alder populations. This study demonstrates that GM is a powerful tool for species delimitation and hybrid detection in the genus Alnus and it can be used for preliminary screening in hybrid zones.


Introduction
Studying the occurrence of gene flow between closely related species is of fundamental importance for inferring the patterns and mechanisms of speciation and hybridization [1][2][3].Hybridization is very important for biodiversity as it greatly contributes to the overall genetic variability in these species [4,5].However, it makes species delimitation more challenging.Spontaneous hybrids naturally occur in sympatric populations of the genus Alnus Mill., but only between very closely related species, namely, Alnus glutinosa × A. incana, A. serrulata × A. rugosa, A. sinuate × A. crispa and A. glutinosa × A. rubra [6].Among these, natural hybrids between A. glutinosa (L.) Gaertn.and A. incana (L.) Moench have been reported to exhibit superior growth rates due to a heterosis effect, characterized by bigger leaf size, increased protein production, improved mechanical properties of wood, better drought resistance and better resistance to Pythium rot [6,7].
Alnus glutinosa, A. incana and A. rohlenae Vít, Douda and Mandák belong to the subgenus Alnus, distinguished by stalked shoot buds and pistillate catkins closed over winter and are pollinated in late winter or early spring [8].Species of Alnus are monecious.Inflorescences (both staminate and pistillate) are formed in autumn, and flowering occurs early in spring, before leaf formation [9,10].Grey alder (A.incana) belongs to the circumpolar boreal floristic element and has a wide distribution in the Northern Hemisphere.It grows The canonical variate analysis (CVA) of the symmetric component revealed a distinct separation of species along CV1, accounting for 93.69% of the variation (Figure 1).Leaves of A. rohlenae, characterized by their circular-obovate shape, long petiole, narrow cuneate base and retuse apex, were distributed along the negative part of the CV1 axis.Leaves of A. incana, characterized by their ovate shape, short petiole, wide-cuneate base and acuminate apex, were distributed mainly along the positive part of the CV1 axis.However, some leaves collected from A. incana trees appeared in the negative part of CV1, resembling the shape of A. rohlenae.Intermediate leaf shapes were predominantly recorded in the population SA.Canonical variate 2, accounting for 3.77% of the variation, partially distinguished populations within species (Figure 1).The shape changes along CV2 were associated with leaf elongation and reflected integrated changes in the leaf's major and minor axes [37].The populations RM and SA were associated with longer and narrower leaves, whereas the populations PK and GR were characterized by leaves that were shorter and wider.Additionally, on the CVA plot, the geographically distant populations PK and GR appeared to be more separated than the populations RM and SA, where the presence of hybrids is assumed.
Pairwise comparisons showed that the differences between all populations were highly significant (Table 1).Based on the sum of interpopulation distances, population PK was the most divergent, while population RM was the least divergent.The linear discriminant analysis revealed that, in the geographically distant populations PK and GR, only 1.3% of leaves from four different individuals of A. incana were misclassified as A. rohlenae (Figure 2A).For the geographically close populations RM and SA, approximately 4% (i.e., 4.4% after cross-validation) of A. incana leaves were misclassified as A. rohlenae (Figure 2B, Supplementary Materials).This 4% of leaves came from nine different individuals.In addition, not a single leaf of A. rohlenae was classified as A. incana.

Secondary Leaf Vein Variations
In the studied populations, the number of pairs of secondary leaf veins (NV) ranged from 7 to 15 for A. incana, aligning with previously reported values within the natural range of species distribution (Table 2).For A. rohlenae, a wider range and a lower limit were recorded compared to [12].The assumption that hybrids could be differentiated

Secondary Leaf Vein Variations
In the studied populations, the number of pairs of secondary leaf veins (NV) ranged from 7 to 15 for A. incana, aligning with previously reported values within the natural range of species distribution (Table 2).For A. rohlenae, a wider range and a lower limit were recorded compared to [12].The assumption that hybrids could be differentiated from parent species by an intermediate number of NV was also tested here.Leaves of A. incana that were misclassified as A. rohlenae by the linear discriminant analysis (LDA) based on shape were indeed characterized by intermediate NV, further indicating their hybrid origin (Table 3).Furthermore, significant differences in NV were recorded between species and among populations of A. incana (Table 3).The partial least squares (PLS) analyses were carried out to examine the relationship between leaf shape and NV for each species.The overall strength of association between blocks was significant for both species.For A. rohlenae, the correlation strength between leaf shape and NV (PLS1) was moderate (R 2 = 0.48, p < 0.05), and for A. incana, it was somewhat stronger (R 2 = 0.65, p < 0.001).Shape changes related to PLS1 are illustrated in Figure 3.For A. rohlenae (Figure 3A), a smaller NV was associated with a wider lamina in the lower part, a narrower lamina in the upper part, a longer petiole, and a retuse apex, whereas a larger NV was associated with a narrower lamina in the lower part, a wider lamina in the upper part, a shorter petiole, and a rounded apex.For A. incana (Figure 3B), shape changes along PLS1 coincided greatly with changes along CV2 (Figure 1): leaves with a smaller NV tended to have shorter and wider lamina, a longer petiole, and a slightly acuminate apex, while leaves with a larger NV were more elongated and narrower, with a shorter petiole and more pronounced acuminate apex.

Leaf Size Variations and Allometry
Alnus species did not differ in leaf centroid size (CS) (F = 1.34, p = 0.25), but the population SA had a significantly smaller CS compared to the other populations (Table 3).
Multivariate regression of the symmetric component on CS showed that shape changes were significantly correlated with changes in size in all the studied populations (p < 0.0001).Nevertheless, allometry accounted for only a small percentage of shape changes within each population (PK-4.490%,RM-3.014%,GR-6.606%,SA-6.017%).
The homogeneity of slopes test revealed that the studied populations exhibited different patterns of allometry (Z = 3.877, p < 0.01): the allometric slope for the population SA was steeper compared to those of the other populations (Figure 4).There was no significant difference in slopes between the populations PK, RM and GR, indicating that the trend of shape change with size was very similar in these populations.In all populations, with an increase in size, the leaves become shorter and wider with a shorter petiole (Figure 5).Major displacements of landmarks were associated with the beginning of the petiole (landmark 1), the widest part of the leaf blade (landmark pairs 4-6 and 8-12), and for the population SA, also with the blade apex (landmark 5) (see Figure 7 for the landmark definition).

inition).
For A. rohlenae, a weak positive correlation between NV and CS was determined in the population PK (R 2 = 0.113, p < 0.05).In the case of A. incana, positive correlations between NV and CS were recorded in both populations: GR (R 2 = 0.286, p < 0.001) and SA (R 2 = 0.316, p < 0.001).

Hybrids in Studied Populations
Considering various leaf morphological traits, hybrids between black and grey alder overlap with one or both parents or hold an intermediate position and appear to be closer to A. incana [14,16,18].Multivariate analysis of leaf morphological traits has shown that hybrids have an intermediate position and are closer to A. incana [14,15,19,20].It was reported by [39] that the germination of a hybrid seed was only successful when A. incana For A. rohlenae, a weak positive correlation between NV and CS was determined in the population PK (R 2 = 0.113, p < 0.05).In the case of A. incana, positive correlations between NV and CS were recorded in both populations: GR (R 2 = 0.286, p < 0.001) and SA (R 2 = 0.316, p < 0.001).

Hybrids in Studied Populations
Considering various leaf morphological traits, hybrids between black and grey alder overlap with one or both parents or hold an intermediate position and appear to be closer to A. incana [14,16,18].Multivariate analysis of leaf morphological traits has shown that hybrids have an intermediate position and are closer to A. incana [14,15,19,20].It was reported by [39] that the germination of a hybrid seed was only successful when A. incana was the maternal parent.When A. glutinosa was the maternal parent, the seed was sterile.Molecular marker analysis confirmed the predominance of A. incana alleles in A. × pubescens [15,21], indicating that hybrids tend to backcross to A. incana.In the studied populations in Serbia, A. incana was also identified as the mother species, because putative hybrids were found only in A. incana populations.Multivariate analyses of leaf shape supported our hypothesis regarding hybridization in geographically close populations: approximately 4% of A. incana leaves in the population SA were misclassified as A. rohlenae.This 4% of leaves belonged to nine different individuals.Given the small percentage of misclassified leaves, we assume that hybrids are primarily characterized by intermediate leaf shapes.However, a minor number of A. incana leaves (four to be exact) were also misclassified as A. rohlenae in the distant, control population GR, indicating that multivariate analyses of leaf shape cannot exclusively distinguish hybrids as well as some individuals with certain forms of phenotypic plasticity.Nevertheless, this method can be used as a preliminary screening tool for hybrids.
Previous studies featuring the determination of the ploidy level have not confirmed the existence of hybrids between A. incana and A. rohlenae.In mixed-ploidy populations in Serbia, [13] identified triploids and suggested that these are the result of a number of ongoing processes [17].They proposed that the triploid origin can be the outcome of hybridization between diploid A. glutinosa s. str.and tetraploid A. rohlenae.However, they did not explore the possibility of hybridization between diploid A. incana and tetraploid A. rohlenae.

Inter-and Intrapopulation Variability in Leaf Traits
Venation patterns and venation densities are variable traits influenced on an ontogenetic and evolutionary level by climatic and environmental parameters (e.g., temperature, water conditions, illumination, wind speed, nutrient status) [40,41].Vein plasticity is expressed within canopies and across environments for a given species, reflecting the gas and water exchange characteristics between the leaves and the atmosphere [41,42].Venation density is strongly positively correlated with the hydraulic conductivity of leaves [43].Leaves acclimated to higher irradiance, temperatures, nutrient supplies and lower water supply exhibit vein traits associated with increased leaf hydraulic conductance and greater drought tolerance [42].The findings of [41] imply that, compared to other leaf traits, the density of major veins (first-and second-order veins) is less sensitive to environment changes.This suggests that, at least to some degree, it is genetically determined for Quercus variabilis Blume.In Alnus species, NV displays less variability at the intrapopulation level and more variability at the interpopulation level with respect to other traits (e.g., leaf length, leaf width, leaf area, perimeter and petiole length) [19,38].Given these observations and the results from this study, NV appears to be a more conservative trait than leaf shape.This, together with the intermediate NV observed in putative hybrids, suggests that it is genetically determined to a great extent.Our findings further indicate that NV is not only species-specific but also population-specific for A. incana.Venation density may (e.g., Acer monspessulanum L.) or may not (e.g., Quercus petraea (Matt.)Liebl.)vary with leaf size [40], depending on the plant species.In the studied populations of A. incana, NV increases with increasing CS.This could be explained by biomechanical and physiological considerations: larger leaves typically require a relatively larger investment in the midrib and an extensive vein network to transport resources [37,44].
Leaf size is considered one of the most variable traits for many plant species [45][46][47], including Alnus spp.[19,38].It is responsive to a range of biotic and abiotic factors.Ref. [7] reported that the leaf surface of hybrids is larger than that of the parental trees due to a heterosis effect.However, in our study, putative hybrids tended to have smaller leaves, showing no evidence of heterosis.Smaller leaves were also recorded in one population of A. incana (SA).Furthermore, the allometric slope for the population SA was steeper, suggesting a higher rate of shape change relative to size in this population.Static allometries arise in response to variations in genetic or environmental regulators of size [48].In the population SA, the narrower and longer leaf shapes with higher NV and smaller CS likely represent a response to an increased evaporative demand in the stressful environment.Many studies have linked smaller leaf sizes to lower annual precipitation and humidity levels [49][50][51].However, according to the hydro-thermal coefficient (Table 4), all the studied populations grow in conditions with sufficient water supply, essential for the genus Alnus.More likely, the cause of the leaf size reduction in the population SA is the prolonged period of inundation in that area.The reduction in leaf transpiration surface enables plants to maintain high photosynthetic activity under flood stress conditions, as shown for Senna reticulata (Willd.)H.S. Irwin and Barneby [52], Populus angustifolia E. James and P. × jackii Sarg.[53].

Plant Material
The studied populations of Alnus spp.are located on Golija Mountain (Golija-Studenica Biosphere Reserve, UNESCO).Plant material was collected in autumn 2017 from mature trees of similar age according to [55].Two populations of A. rohlenae were previously classified as A. glutinosa s. str.[56], due to a lack of data on A. rohlenae's exact distribution.Morphologically, the two species are very similar, making determination without estimating ploidy level challenging.Following the publication on the exact distribution of the tetraploid A. rohlenae [13] and the private consultation with Professor Dr. Dmitar Lakušić from the University of Belgrade, these populations were correctly reclassified.The sample size and geographic locations of the studied populations are described in Table 4 and Figure 6.

Digitalization of Leaves
A total of 1315 leaves were sampled from the middle part of the branches.A priori assumptions about a hybrid origin were not made in order to reduce subjectivity.Analyzing several leaves from the same individual is considered the best method for establishing hybridization, because it is known that a single branch can exhibit leaf shapes characteristic of both species [14,19], as well as variable intermediate shapes [18,20,21].Each leaf was scanned with the abaxial surface facing up with an Epson Perfection V370 scanner at 600 dpi resolution.The scanned images were used to obtain 2D landmark configurations.In our previous study [56], we developed a landmark configuration to describe the leaf shape of Alnus species (Figure 7).The digitalization of specimens was carried out in tpsUtil version 1.70 ( [57], downloaded on 9 September 2016) and tpsDig version 2.26 ([58], downloaded on 9 September 2016).
In addition, the number of pairs of secondary leaf veins (NV) was measured for each leaf.Previous studies based on the traditional morphometric approach have identified NV as the most distinctive feature of hybrids, noting that it varies within a very narrow, intermediate range [7,14,[18][19][20].This trait was used to compare with the results obtained from shape and size analysis.

Digitalization of Leaves
A total of 1315 leaves were sampled from the middle part of the branches.A priori assumptions about a hybrid origin were not made in order to reduce subjectivity.Analyzing several leaves from the same individual is considered the best method for establishing hybridization, because it is known that a single branch can exhibit leaf shapes characteristic of both species [14,19], as well as variable intermediate shapes [18,20,21].Each leaf was scanned with the abaxial surface facing up with an Epson Perfection V370 scanner at 600 dpi resolution.The scanned images were used to obtain 2D landmark configurations.In our previous study [56], we developed a landmark configuration to describe the leaf shape of Alnus species (Figure 7).The digitalization of specimens was carried out in tpsUtil version 1.70 ( [57], downloaded on 9 September 2016) and tpsDig version 2.26 ([58], downloaded on 9 September 2016).

Assessment of Measurement Error
Measurement error can be quantified using Procrustes ANOVA [59][60][61].To estimate the measurement error associated with digitalization, a pilot study was carried out with 400 specimens (5 trees per population, 10 leaves per tree and 2 replicates per leaf).A full Procrustes fit and Procrustes ANOVA were performed for each population separately in MorphoJ ( [62], 2011, downloaded on 27 March 2017).The results showed that for each population, the digitizing error was small relative to the smallest level of variation (individual-by-side interaction) (Table 5).Therefore, for the original dataset of 1315 specimens, replicate digitalization was not necessary.
Fluctuating asymmetry, estimated by the individual-by-side interaction, was significant in all the studied populations.Fluctuating asymmetry represents the difference of each individual s asymmetry from the average asymmetry in the whole population and reflects the combined effects of phenotypic plasticity and developmental instability [63].
Given that this study s subject is interspecific hybridization, we focused on the symmetric component of variation, which reflects genetic variation and phenotypic plasticity.In addition, the number of pairs of secondary leaf veins (NV) was measured for each leaf.Previous studies based on the traditional morphometric approach have identified NV as the most distinctive feature of hybrids, noting that it varies within a very narrow, intermediate range [7,14,[18][19][20].This trait was used to compare with the results obtained from shape and size analysis.

Assessment of Measurement Error
Measurement error can be quantified using Procrustes ANOVA [59][60][61].To estimate the measurement error associated with digitalization, a pilot study was carried out with 400 specimens (5 trees per population, 10 leaves per tree and 2 replicates per leaf).A full Procrustes fit and Procrustes ANOVA were performed for each population separately in MorphoJ ( [62], 2011, downloaded on 27 March 2017).The results showed that for each population, the digitizing error was small relative to the smallest level of variation (individual-by-side interaction) (Table 5).Therefore, for the original dataset of 1315 specimens, replicate digitalization was not necessary.Fluctuating asymmetry, estimated by the individual-by-side interaction, was significant in all the studied populations.Fluctuating asymmetry represents the difference of each individual's asymmetry from the average asymmetry in the whole population and reflects the combined effects of phenotypic plasticity and developmental instability [63].Given that this study's subject is interspecific hybridization, we focused on the symmetric component of variation, which reflects genetic variation and phenotypic plasticity.

Data Analysis
The full Procrustes fit of 1315 leaf configurations from pooled populations was performed in MorphoJ [62], taking into account that leaves have object symmetry.The main trends in shape variation were captured by the symmetric component of variation, which accounted for 80.77%.
Inter-and intraspecific variability in leaf shapes were assessed by canonical variate analysis (CVA) performed on the covariance matrix of the symmetric component.The delimitation of the species was further assessed with linear discriminant analysis (LDA).
Differences between species and populations in NV and centroid size (CS) were assessed with nested ANOVA in Minitab 17 statistical Software (2010, Minitab Inc., State College, PA, USA).Individual trees were assigned as random factors and nested in populations (fixed factors).The Pearson product-moment correlation coefficient was used to examine the relationship between NV and CS.
Partial least squares (PLS) analysis was used to investigate the covariation between the symmetric component of shape variation as block 1 and NV as block 2. Prior to analysis, data were averaged per tree and NV values were square root-transformed.
CVA, LDA and PLS analysis were carried out in MorphoJ.
Static allometry was examined with the multivariate regression of the symmetric component of shape variation on CS for every population separately.The homogeneity of allometric slopes between populations was tested using the R-package GEOMORPH (version 3.0.5,[64]) in R statistical software (version 3.2.3., downloaded on 30 January 2016) with the 'procD.allometry'function.Allometric slopes are displayed using a predicted shape regression which calculates the predicted values of a regression of shape on size and plots the first principal component scores of these predicted values as an allometric trend [65].

Figure 2 .
Figure 2. Linear discriminant analysis of leaf shape: histograms of linear discriminant values of leaves from geographically distant populations PK and GR (A) and geographically close populations RM and SA (B).Population abbreviations: PK-Prilički kiseljak, RM-Rimski most, GR-Golijska reka, SA-Sastavci.

Figure 2 .
Figure 2. Linear discriminant analysis of leaf shape: histograms of linear discriminant values of leaves from geographically distant populations PK and GR (A) and geographically close populations RM and SA (B).Population abbreviations: PK-Prilički kiseljak, RM-Rimski most, GR-Golijska reka, SA-Sastavci.

Figure 3 .
Figure 3. Leaf shape changes (black line) related to the number of pairs of Consensus leaf shapes (grey line) for A. rohlenae (A) and A. incana (B).

Figure 3 .
Figure 3. Leaf shape changes (black line) related to the number of pairs of secondary leaf veins (NV).Consensus leaf shapes (grey line) for A. rohlenae (A) and A. incana (B).

Figure 6 .
Figure 6.Geographical locations of the studied populations in the Moravica River basin.

Figure 6 .
Figure 6.Geographical locations of the studied populations in the Moravica River basin.Voucher specimens have been deposited in the Herbarium of the Department of Ecology at the Institute for Biological Research, University of Belgrade, Serbia.

Table 1 .
Mahalanobis and Procrustes distances calculated from canonical variate analysis symmetric component of 1315 Procrustes-aligned 15-landmark configurations of leaves.

Table 1 .
Mahalanobis and Procrustes distances calculated from canonical variate analysis of the symmetric component of 1315 Procrustes-aligned 15-landmark configurations of leaves.

Table 2 .
Literature data for number of pairs of secondary leaf veins (NV) within the natural range of species distribution.
a-d homogenous groups obtained by Tukey and Bonferroni pairwise co dence level.