Dissecting Taxonomic Variants within Ulmus spp. Complex in Natural Forests with the Aid of Microsatellite and Morphometric Markers

: Spontaneous hybrids between the native elms (genus Ulmus L.) have been observed in the forests of Europe. Gene conservation raises questions regarding the genetic background for the complex morphology and taxonomy of elms. Our objective was to dissect morphological and genetic variation in the natural swamps of Ulmus species groups in Lithuanian forests with the aid of leaf morphology and microsatellite (SSR) markers. We sampled leaves from 189 elms at 26 locations to grasp the phenotypic diversity in variable natural habitats in Lithuanian forests. We assigned the elms into six taxonomic and genetics groups based on 31 leaf morphology parameters and tested the genetic differentiation between these six groups at six nuclear SSR loci by using Bayesian and genetic distance-based clustering. The genetic and leaf morphometric analyses of putative elm hybrid swamps indicated a low genetic exchange between U. laevis Pall. and the other Ulmus groups. The genetic and morphometric data supported the differentiation of U. glabra Huds. and U. glabra (female) × U. minor Mill. (male) spontaneous hybrids. In addition, the results of the genetic analysis also conﬁrmed the high level of genome sharing among U. minor and U. minor subsp. minor Richens., where leaf morphology failed to differentiate genetically discrete groups. For gene conservation, we would suggest considering separate gene conservation units selected based on leaf and stem morphology for U. laevis , U. glabra , U. glabra × minor , and the U. minor species complex.


Introduction
Interspecific hybridization is a common feature in plants with vital evolutionary consequences [1][2][3]. There are variable opinions on the role of hybridization in plants. First of all, interspecific hybridization has often been considered a source for genetic and phenotypic novelties and a force for evolution, but hybridization can also cause genetic erosion, threaten species integrity, and lead to species extinction [2,4,5]. Hybridization can lead to "evolutionary innovation" through the generation of new genotypes, an increase in heterosis, a pool of constant genetic variation, and a decrease in genetic load [6][7][8][9][10]. Hybrid zones often represent regions with high genetic variation and unique combinations of alleles, where selection can be intense and evolution can be rapid [6,11]. Morphological data can provide information about evolutionary and environmental phenomena associated with hybridization, since hybrids are not always intermediate between parental species, but they often exhibit extremes and new traits [12]. The percentage of plants with distinctive or new traits increases with subsequent generations of hybrids, as does the percentage Microsatellite markers (simple sequence repeats-SSRs) are short, tandemly repeated DNA sequences that are valuable molecular tools for various eco-genetic studies due to the high degree of polymorphism, abundance, co-dominance, and easy transferability among the species and laboratories [37][38][39][40]. SSRs are very often in use for studying related trees species and their hybrids, but the development of new SSRs is an expensive and timeconsuming procedure. Thus, the transferability of nuclear microsatellite loci across species is important and depends on various factors, e.g., the level of hybridization, genome size, different breeding systems, and evolution history. [38,39]. The successful transferability of SSRs among species enables us to use them for related species [39]. There have been a number of studies of Ulmus species where the same microsatellite markers are in use, e.g., Zalapa et al. [41] developed 15 microsatellite loci in U. rubra and tested their crossamplification in U. pumila, which was successful-all 15 primers were amplified in both species; Zalapa et al. [42] used 13 nSSRs to study hybridization among U. rubra and an invasive tree-U. pumila; Nielsen and Kjaer [43] tested 22 nSSRs and selected seven to study U. laevis from Denmark, southern Sweden, and Finland; Venturas et al. [44] tested 19 nSSRs and selected nine to study U. laevis in Spain; Bertolasi et al. [45] studied gene flow between local U. minor and introduced U. pumila populations based on six microsatellite markers and found that species could hybridize when in sympatry; Buiteveld et al. [46] tested 23 nSSRs and selected ten to study U. minor in the Netherlands; and Martin del Puerto et al. [47] tested 22 nSSRs and selected 11 to study U. glabra in the Iberian Peninsula. Thus, a high number of studies have proven the successful transferability of SSRs among several Ulmus species.
There have been several studies [46,[48][49][50] on hybridization among Ulmus spp., e.g., Mittempergher and Porta [51] presented data on the cross-ability and rate of selfing derived from crossing trials among 11 elm species. Studies showed that barriers to hybridization among species were weak, with the success of several combinations dependent on malefemale interaction and the parental individual. An exceptional cross-ability barrier was found between U. laevis and the other Ulmus species. Zalapa et al. [42] identified a surprisingly large number of hybrids among an invasive U. pumila and local U. rubra hybrid individuals in the United States based on genetic analyses. Brunet et al. [52] studied the hybridization of Siberian elm (U. pumila) with native field elm (U. minor) in Italy. They used genetic markers to examine the extent of hybridization between these two species and to determine the pattern of introgression, and they found that hybrids between U. pumila and U. minor are quite common. Buiteveld et al. [46] used microsatellite markers to describe clonal diversity and structure, as well as to calculate genetic diversity parameters, in Dutch U. minor populations. At four locations, they found some individuals that might have been hybrids or at least not pure U. minor specimens based on STRUCTURE clustering analysis including parental species. Recently, Hirsch et al. [3] studied interspecific hybridization between the Siberian elm (U. pumila) and native elm species in the Midwestern United States, Italy, and Spain. They used a set of nuclear microsatellite markers and the program STRUCTURE to detect interspecific hybridization and determine the populations' genetic structure. DNA marker-based findings proved the presence of intraspecific hybridization. Hirsch et al. [3] supplemented evidence from previous studies and reported on intraspecific hybridization within Ulmus genus.
The aims of this work were (i) to elucidate the morphological differences between the critical groups of native elms, indicative for their taxonomic identity; (ii) to determine the morphological boundaries of taxa; and (iii) use the latter findings as the basis for the further study of species determination, hybrid identification, and genetic diversity based on nuclear microsatellite markers.

Objects and Material
Putative natural hybrids between smooth-leaved elm (U. minor subsp. minor), field elm (U. minor), wych elm (U. glabra), and European white elm (U. laevis), along with the individuals of pure species, were sampled for the laboratory examination of leaf morphology and DNA microsatellite analysis in natural mixed forests of Lithuania. This choice was based on our experience in exploring local elms [35,53]. The smooth-leaved elm (U. minor ssp. minor) is the most common type of field elm in continental Europe today [54,55]. In Europe, this subspecies has a more southerly distribution than U. glabra, and it is unknown in Denmark, Sweden, and Norway as a wild tree, though it is said to occur in the Baltic [56]. The smooth-leaved elm is identified by mature leaves that are smooth, glossy, bright green from above, and very unequal at the base [53] The two groups of hybrid elm trees were identified based on the visual examination of leaf morphology traits and stem morphotypes on site: (a) U. glabra × U. minor and (b) U. minor subsp. minor × U. glabra. The study sites were selected all over the country based on the presence of elm species and the elm hybrid swarms (Figure 1). At each site, study plots of 40 m radius were established, usually on moist fertile soils suitable for elm species. During the summers of 2018-2020 in each sample plot, 5-10 leaves were collected from long shoots at a height of 5-8 m of randomly selected mature elm trees (diameter > 15 cm) with aim of sampling approximately 30 trees for each of the four elm species and each of the two hybrid groups. Jeffers [27] suggested that the leaves on long shoots show far less sharp differentiation (in leaf length and width, petiole length, leaf base asymmetry, teeth number, teeth width, length, and depth) between different kinds of elm than those on short shoots. A minimum distance of 20 m was maintained among the sampled trees to avoid clones or close relatives.

Objects and Material
Putative natural hybrids between smooth-leaved elm (U. minor subsp. minor), field elm (U. minor), wych elm (U. glabra), and European white elm (U. laevis), along with the individuals of pure species, were sampled for the laboratory examination of leaf morphology and DNA microsatellite analysis in natural mixed forests of Lithuania. This choice was based on our experience in exploring local elms [35,53]. The smooth-leaved elm (U. minor ssp. minor) is the most common type of field elm in continental Europe today [54,55]. In Europe, this subspecies has a more southerly distribution than U. glabra, and it is unknown in Denmark, Sweden, and Norway as a wild tree, though it is said to occur in the Baltic [56]. The smooth-leaved elm is identified by mature leaves that are smooth, glossy, bright green from above, and very unequal at the base [53] The two groups of hybrid elm trees were identified based on the visual examination of leaf morphology traits and stem morphotypes on site: (a) U. glabra × U. minor and (b) U. minor subsp. minor × U. glabra. The study sites were selected all over the country based on the presence of elm species and the elm hybrid swarms (Figure 1). At each site, study plots of 40 m radius were established, usually on moist fertile soils suitable for elm species. During the summers of 2018-2020 in each sample plot, 5-10 leaves were collected from long shoots at a height of 5-8 m of randomly selected mature elm trees (diameter > 15 cm) with aim of sampling approximately 30 trees for each of the four elm species and each of the two hybrid groups. Jeffers [27] suggested that the leaves on long shoots show far less sharp differentiation (in leaf length and width, petiole length, leaf base asymmetry, teeth number, teeth width, length, and depth) between different kinds of elm than those on short shoots. A minimum distance of 20 m was maintained among the sampled trees to avoid clones or close relatives.

Statistical Analysis of Leaf Characteristics
One-way ANOVA was used to test the effect of species on leaf morphology traits with the XLSTAT2020 program. The traits with the highest F values from the ANOVA were used for multivariate principal component analysis (PCA) analysis with PC-ORD5 soft. Based on the performed statistical analysis, we identified the key leaf morphology traits for the discrimination among and identification of hybrids of elm species.

Molecular Data Analysis
Genetic diversity parameters were calculated for the six Ulmus spp. species groups (number of different alleles (Na), number of effective alleles (Ne), observed/expected/ unbiased expected heterozygosity, and fixation index (F)) based on six microsatellite loci using the GenAlEx 6.5 software [61]. Two loci were rejected from the genetic diversity analysis because they were found to be amplified in only one of the species (speciesidentification-suitable loci). An analysis of molecular variance (AMOVA) was performed using GenAlEx 6.5 [61] (for the significance test, we used 9999 random permutations). The estimation and visualization of private alleles were performed in the Poppr R package [62]. Allelic richness (Ar) was estimated with the FSTAT 2.9.3. software [63]. The software estimates allelic richness per locus, sample, and samples overall. Allelic richness is a measure of the number of alleles independent of sample size, thus allowing for comparison between different sample sizes among populations. The lowest number of samples (12) was used for rarefaction. Missing data estimation among the loci and target groups of Ulmus spp. performed and visualized by the R package poppr [62]. In addition, we used DAPC to examine the clustering of individuals: first with all six groups and later with five species groups (R package adegent 2.0.0 [64,65]). To test the associations among the species groups based on traditional Nei's 1978 genetic distances, we ran UPGMA cluster analysis with the R package poppr with 10,000 bootstrap replicates [66].
We used eight microsatellite loci for species differentiation and hybrid identification. First, the clustering algorithm in STRUCTURE 2.3.4 [67] was run to assess the genetic structure of the dataset of six Ulmus spp. groups for each K ranging from 1 to 12 using 100,000 Markov chain Monte Carlo iterations with a burn-in of 100,000 and 20 replicates per run. The admixture model was used and allowed for the correlation of allele frequen- cies among clusters. The approach by Evanno et al. [68] in STRUCTURE HARVESTER v0.6.94 [69] was used for selecting the most appropriate K clusters. The STRUCTURE analysis was performed for all samples from studied six groups (the group of U. laevis, the group of U. glabra and U. minor, and the groups of potential hybrids-133 individuals in total). After the clear identification of U. laevis form the other Ulmus spp. groups, European white elm was removed from further analysis. To improve the quality and accuracy of hybrid identification, we applied data filtering, and all individuals containing missing data in three or more loci were removed from further analysis. In total, 106 individuals remained for possible hybrid identification. Then, a second round of the clustering in STRUCTURE 2.3.4 was run [67] with the same parameters as above, only with K ranging from 1 to 10, including 106 individuals from five groups: UG-U. glabra; UM-U. minor; UMm-U. minor subsp. Minor; UMm × UG-U. minor subsp. Minor × U. glabra; and UG × UM-U. glabra × U. minor.
To facilitate the detection of putative interspecific hybrids in the populations, genetically pure individuals of the two respective parental species as reference populations were used. This method has been used in several studies [3,52] to identify intra-and interspecific hybridization in invasive Siberian elm and other Ulmus spp. Thus, sampled individuals were sorted according to leaf characteristics and sample locations into the most probably pure individuals of each species. Then, we used the program STRUCTURE, which uses the Bayesian clustering approach as the model-based clustering algorithm (version 2.3.4; [67]), to assign individuals to pure species or identify possible hybrids when pools were mixed. When two pure parental species were sampled as references, it was expected that the optimal value of K would consist of two genetic clusters (K = 2). This could be confirmed by testing values of K from one up to the number of populations in the respective groups using the STRUCTURE HARVESTER software [69], and we selected the optimum K following the method of Evanno et al. [68]. The program STRUCTURE generates an admixture coefficient (q) that represents the proportion of an individual's genotype that originates from each of the K genetic clusters. STRUCTURE can be run with the option ANCESTDIST, which computes the 95% posterior probability for each q value, equivalent to a 95% confidence interval. Following Blair and Hufbauer [10], individuals were classified as hybrids if their q value was <0.90. If an individual's proportion did not include one, introgression likely occurred [10]. In addition, species-specific alleles identified in the reference datasets could help to confirm the identification of hybrid individuals. Because we had two parental species (U. glabra and U. minor), we expected the optimal value of K to consist of two genetic clusters (K = 2). For each STRUCTURE analysis, we used an admixture model with 100,000 burn-in iterations, 100,000 Markov chain Monte Carlo repetitions, and 20 replicates at each level, and we allowed for the correlation of allele frequencies among clusters. The online software CLUMPAK was used to identify clustering modes and packaging population structure inferences across K, as well as to visualize the clustering [70].
Finally, we used the Bayesian algorithms provided in NewHybrids v.1.1 beta [71], which performed the independent classification of individuals as U. minor, U. glabra, or a hybrid based on their genotypic profiles, and it helped to further classify the hybrids into specific categories. We considered the following hybrid classes: first-(F1) and secondgeneration (F2) hybrids and first-(0_Bx) and second-generation (1_Bx) backcrosses. The NewHybrids algorithm was run with Jeffreys-like priors with 500,000 iterations following a 500,000-iteration burning. At the end, we combined the information obtained from Structure and NewHybrids to determine the specific hybrid class to which an individual tree was most likely to belong.

Leaf Morphology Variation
The ANOVA revealed significant species effects on all the leaf morphology traits except for two traits of Pl/I and GKPP (pubescens character of the corners of the main vein;  Table 2). The highest three F values were obtained for the leaf blade pubescens score (L), the secondary leaf vein branching score (L), and the leaf petiole lower half pubescens (KA).
The PCA analysis was performed for the morphology traits with the highest values of R 2 and F ( Table 2, rows 1-13). The results showed high correlation coefficients for the above-mentioned characteristics (Table 3). Table 3. Description of three main principal components (PCs) from the PCA on leaf traits with the highest R 2 from the ANOVA. Pearson's (R) and Kendall's correlations (tau) with the PCA ordination axes; n = 189. The PCA plot individual tree values against two major PCs accounted for 66.32%, 21.82%, and 9.66% of the total variance in the first, second, and third PC axes, respectively. Based on leaf morphology traits, pure elm species were divided into three groups: U. glabra, U. laevis, and U. minor subsp. minor with U. minor. The putative spontaneous hybrids between these elm species were located throughout the sampling sites (Figure 3a,c). A25 showed a strong statistical significance with Per and KP in the first, second, and third axes (Figure 3b,d).

Morphologic
For wych elm, the values of all investigated leaf morphology traits were high, except for the asymmetry at 25% blade (25K) length and the pubescens of the underside of the leaf blade (L) ( Figure S1). Though the field elm and smooth-leaved elm were found to occupy intermediate positions between wych elm and European white elm. Their L and underside of the petiole (KA) were found to usually be without pubescens. Smooth-leaved elm differed from field elm by a larger blade maximum width (PlMax) and secondary vein stretch (IS). The numerals of all characteristics of hybrids of wych elm and field elm (and smooth-leaved elm) varied from very small to very big. European white elm leaves were found to have no branches on the secondary veins, and they were the most asymmetrical.

Genetic Diversity
In total, 133 trees from six target tree groups were sampled and analyzed using six microsatellite loci. All six nuclear microsatellites were polymorphic and amplified 72 alleles in total. The number of alleles per locus varied from 7 at locus UR158 to 21 at locus Ulmi1165. Loci UR158 and Ulm6 were least polymorphic with the lowest expected heterozygosity (He), allelic richness (Ar), and number of effective alleles (Ne) (results not shown).
The mean number of alleles (Na) varied from 4.33 in U. laevis to 8.0 in U. glabra, with an overall average of Na = 6.89. The mean number of effective alleles (Ne) varied from 2.53 in U. laevis to 4.55 in U. glabra, with an overall average of Ne = 4.01. Allelic richness (Ar) varied from 3.6 in U. laevis to 6.68 in U. minor, with an overall average of Ar = 6.0. Expected heterozygosity (He) was high for all investigated species, except for U. laevis with He = 0.567. Most of the allelic diversity parameters were markedly lower in U. laevis in comparison with the other elm species (Table 4 and Figure S2). The deviation from random mating was markedly stronger in the putative hybrid groups than in the pure elm species (Fis from 0.21-0.18 for hybrids vs. Fis from −0.06 to 0.13 pure species; Table 4). In total, 19 private alleles were present and differently distributed over six elm species, with the highest number of private alleles observed in the U. laevis and U. glabra tree groups ( Figures S2 and S3).  For wych elm, the values of all investigated leaf morphology traits were high, except for the asymmetry at 25% blade (25K) length and the pubescens of the underside of the leaf blade (L) ( Figure S1). Though the field elm and smooth-leaved elm were found to occupy intermediate positions between wych elm and European white elm. Their L and underside of the petiole (KA) were found to usually be without pubescens. Smooth-leaved elm differed from field elm by a larger blade maximum width (PlMax) and secondary vein stretch (IS). The numerals of all characteristics of hybrids of wych elm and field elm (and smooth-leaved elm) varied from very small to very big. European white elm leaves were found to have no branches on the secondary veins, and they were the most asymmetrical.

Species Genetic Differentiation
We identified two microsatellite loci discriminating between U. laevis and the remaining Ulmus species: locus Ulm198 did not amplify in U. laevis, and locus Ulm19 amplified only in U. laevis, with some exception in U. minor ( Figure S4).
In agreement with the leaf morphology traits, the Bayesian clustering results based on the molecular data revealed a clear separation of U. laevis from the remaining elm species (Figure 4, upper plot). K = 2 was defined as the most likely number of genetic clusters (deltaK = 261.6; Figure S5 and Table S2). The discriminant analysis of principal components (DAPC) supported clear U. laevis genetic differentiation from the remaining elm species when plotted on the ordination axes ( Figure 4, lower plot).

Hybrid Identification
After excluding U. laevis and data filtering, we obtained a subset of 106 individuals for the identification of the parental species of U. minor and U. glabra, as well as their hybrids. Based on the leaf morphology, 28 of these trees were classified as U. glabra (UG), 19

Hybrid Identification
After excluding U. laevis and data filtering, we obtained a subset of 106 individuals for the identification of the parental species of U. minor and U. glabra, as well as their hybrids. Based on the leaf morphology, 28 of these trees were classified as U. glabra (UG), 19 were classified as U. minor (UM), 21 were classified as U. minor subsp. minor (UMm), 25 were classified as the U. glabra × U. minor (UG × UM) hybrid group, and 13 were classified as the U. minor subsp. minor × U. glabra (UMm × UG) hybrid group. The STRUCTURE clustering revealed two genetic clusters that best-explained the molecular genetic variation in this subset of 106 samples (delta K = 64.9; Figure S6 and Table S3). For two genetic clusters, STRUCTURE clustering clearly separated the group of U. glabra from U. minor and the putative hybrid groups ( Figure 5, upper plot). Furthermore, the putative UG × UM hybrids (UG as the female parent) contained stronger genetic associations with U. glabra than the UMm × UG hybrid group (UG as the male parent; Figure 5, upper plot). Based on this method, only a few individuals were assigned as hybrids, with equal membership to the UG and UM genetic clusters. These results indicated that leaf morphological traits may reliably distinguish the species of U. glabra and U. minor but fail to discriminate further within U. minor sensu lato. The leaf traits may be sensitive to the maternal-paternal identity of the hybrids.   [67]. In the CLUMPAK plot, everything is represented by a thin vertical line divided into K = 2 colored segments that represent the individual's estimated membership fractions in these two clusters. Black lines separate individuals from different species groups. The first cluster (orange) displays the clear separation of the U. glabra (UG) group, and the second cluster (blue) indicates groups of U. minor (UM) and U. minor subsp. minor (UMm). The lower plot illustrates the genetic structure of five Ulmus species groups based on two major principal components from the DAPC analysis (R package adegent 2.0.0) [64,65].
The DAPC analysis of the subset of five elm species groups separated three major genetic groups ( Figure 5, lower plot): (1) the U. glabra group (UG, rightmost group in the plot), (2) putative hybrids U. glabra × U. minor (intermediate group in the plot), and (3) a U. minor group containing the UM, Umm, and UMm × UG species groups (the leftmost group in the plot). These results suggested (a) genetic separation between U. glabra and U. minor groups and (b) confirmed the leaf traits as reliable indicators for putative hybridization between U. glabra and U. minor, with the latter species as the mother tree. However, DAPC could not confirm reliable leaf marker traits for identifying the hybrids with U. minor as the mother tree species.
The UPGMA clustering well-reflected the results of the Bayesian clustering and DAPC via the highly reliable separation of the pure species of U. laevis, U. glabra, the UG × UM hybrid complex, and the UM species complex ( Figure 6). However, it is worth noting the separation of the UMm × UG hybrid cluster from the UM × UMm cluster with a 97% bootstrap significance ( Figure 6).  [66]. The significance of branch nodes was tested with 10,000 bootstraps among loci (indicated by% of bootstraps separating a given branch).
Finally, the NewHybrids program identified nine F2 hybrids between UG × UM (with p ≥ 0.90) and an additional 22 individuals as putative F2 hybrids with a lower p value (p ≥ 0.50). In total, 11 hybrids were identified among 36 morphology-based putative hybrids. In addition, 20 individuals were identified as F2 hybrids among U. minor, U. minor subsp. minor, and U. glabra groups. In contrast to STRUCTURE clustering, where three hybrids were identified within the U. glabra (UG) group, the NewHybrids identified markedly more hybrids (10 hybrids). All the hybrids identified by the NewHybrids software were classified as second (F2) generation hybrids (with p ≥ 0.50).

Discussion
Three species of the Ulmus complex sensu lato naturally occur in Lithuania: U. glabra, U. minor, and U. laevis. Our country is at the northern margin of the natural range of these species. For U. minor, the northern boundary of the natural range coincides with 55° latitude and runs along the southern edge of the Baltic Sea through Lithuania towards the Urals. During field research, we did not find spontaneous hybrids between U. laevis and other Ulmus species based on morphology traits, thus agreeing with other studies on U. laevis hybridization. U. laevis belongs to the Blepharocarpus section of the genus, but the other two European elm species, U. glabra and U. minor, belong to the subgenus Ulmus. U. laevis does not easily hybridize with the other European elm species and is self-incompatible [51]. Townsend [72] determined that U. laevis is a diploid that exhibits reproductive isolation. The results of interspecific pollinations of U. americana with U. laevis have indicated that fertilization is prevented by a reproductive barrier that inhibits pollen tube growth at a stigmatic surface, regardless of belonging to the same Blepharocarpus section [73]. Meanwhile, we identified many hybrid individuals between U. minor and U. glabra. Hybrids of these species have also been found in many European countries, such as Spain, Belgium, England, France, and Slovenia [27,60,74]. The taxonomy of European elm trees  [66]. The significance of branch nodes was tested with 10,000 bootstraps among loci (indicated by% of bootstraps separating a given branch).
Finally, the NewHybrids program identified nine F2 hybrids between UG × UM (with p ≥ 0.90) and an additional 22 individuals as putative F2 hybrids with a lower p value (p ≥ 0.50). In total, 11 hybrids were identified among 36 morphology-based putative hybrids. In addition, 20 individuals were identified as F2 hybrids among U. minor, U. minor subsp. minor, and U. glabra groups. In contrast to STRUCTURE clustering, where three hybrids were identified within the U. glabra (UG) group, the NewHybrids identified markedly more hybrids (10 hybrids). All the hybrids identified by the NewHybrids software were classified as second (F2) generation hybrids (with p ≥ 0.50).

Discussion
Three species of the Ulmus complex sensu lato naturally occur in Lithuania: U. glabra, U. minor, and U. laevis. Our country is at the northern margin of the natural range of these species. For U. minor, the northern boundary of the natural range coincides with 55 • latitude and runs along the southern edge of the Baltic Sea through Lithuania towards the Urals. During field research, we did not find spontaneous hybrids between U. laevis and other Ulmus species based on morphology traits, thus agreeing with other studies on U. laevis hybridization. U. laevis belongs to the Blepharocarpus section of the genus, but the other two European elm species, U. glabra and U. minor, belong to the subgenus Ulmus. U. laevis does not easily hybridize with the other European elm species and is self-incompatible [51]. Townsend [72] determined that U. laevis is a diploid that exhibits reproductive isolation. The results of interspecific pollinations of U. americana with U. laevis have indicated that fertilization is prevented by a reproductive barrier that inhibits pollen tube growth at a stigmatic surface, regardless of belonging to the same Blepharocarpus section [73]. Meanwhile, we identified many hybrid individuals between U. minor and U. glabra. Hybrids of these species have also been found in many European countries, such as Spain, Belgium, England, France, and Slovenia [27,60,74]. The taxonomy of European elm trees has given rise to much debate, especially regarding the U. minor species complex and the nature and frequency of the hybridization of U. minor and U. glabra [75].
In our previous study [35], we used 14 leaf morphological parameters, and we found out that natural hybrids between field elm and wych elm (U. × hollandica) do occur in mixed forests of Central Lithuania. This study covered more elm species (we additionally studied field elm and European white elm) and more leaf morphological characteristics (we added pubescens of leaves). In this study, we confirmed that elms can easily be distinguished by leaf morphometric characteristics. The inclusion of leaf pubescence strength in morphological studies helps to better characterize individual elm species. According to our results, the U. laevis has a lot of pubescens on the main and secondary veins, as well as on the underside of the leaf. The upper half of the petiole is not always rich with pubescens, and the lower half is less frequent. Secondary veins are not branched. This elm is the most asymmetrical of elms. U. glabra secondary veins are usually all-branched, both the primary and secondary veins usually have a lot of pubescens, and the petiole has extremely many pubescens. U. minor and U. minor subsp. minor are closer to U. glabra because the leaf petioles and veins have much less pubescens than U. glabra. U. glabra and U. minor (and U. minor subsp. minor) hybrids have similar leaf traits to U. minor (or U. minor subsp. minor), but, in contrast to the pure species, the pubescens is absent along the edges of the main veins in their hybrids.
Examples of hybrid elm showed that their morphological and genetic boundaries do not fit, their genetic boundaries are much wider than the phenotypic, or vice versa. Morphological investigations were carried out to help to distinguish the species of elm from one another and to evaluate the abundance of hybrids, but morphological studies were not sufficient to investigate the degree of hybridization of the elm species.
We identified a set of eight nuclear microsatellite markers for the efficient study of genetic diversity and hybrid identification in the Ulmus species complex. Our findings were in good agreement with a number of studies where microsatellites or other DNA markers were used to study the hybridization in Ulmus species [3,42,45,46,52,76]. The results presented by L. Mittempergher and N. la Porta [51] on the cross-ability and rate of selfing from artificial mating trials among 11 elm species showed that hybridization barriers among Ulmus species were weak, with the success of several combinations dependent on male-female interactions and the parental individual. However, these studies reported that an exceptional cross-ability barrier was found between U. laevis and other Ulmus species. This finding may explain the significant genetic differentiation of U. laevis from other Ulmus species obtained in our study. In addition, this low interspecific crossing rate may have led to relatively lower genetic diversity and high frequency of private alleles in U. laevis in our study (Table 4 and Figure S2). A similar genetic diversity with a low mean number of alleles and a low heterozygosity (He) in U. laevis was observed at its northern distribution range in Denmark by Nielsen and Kjaer [43]. A significantly higher genetic diversity among other two parental Ulmus species and their possible hybrids was observed in our study. However, Martín del Puerto et al. [47] found a much lower genetic diversity (He) among the U. glabra populations in the Iberian Peninsula. The lower genetic diversity results of U. glabra in the Iberian Peninsula can be explained by small populations size and isolation. Overall, estimates of expected heterozygosity (He) in our Ulmus species groups (U. minor, U. minor subsp. minor, U. glabra, and their hybrids) were the same or slightly higher in comparison to other European field elm populations, e.g., Brunet et al. [52] found He = 0.59 in Italy, Fuentes-Utrilla et al. [77] found He = 0.333-0.592 in Balearic Islands, Bertolasi et al. [45] found He = 0.671 in Italy, Zebec et al. [48] found He = 0.418-0.642 in five natural field elm populations in Croatia, Buiteveld et al. [46] revealed a moderate genetic diversity (He = 0.483-0.628) in the Netherlands, and Collada et al. [60] found He = 0.49-0.73 in six Spanish populations. In addition, Venturas et al. [44] found He = 0.43-0.45 in two U. laevis populations in the Iberian Peninsula, and Whiteley et al. [59] found He = 0.08-0.74 in a group of U. laevis from Sweden (Öland). However, when comparing the results, we should consider the variation in sample size, geography, and number of loci used.
In our study, evidence for hybridization between U. glabra and U. minor was found from morphological and genetic backgrounds. Similarly, hybridization events have been reported in other studies on U. glabra and U. minor [35,74] and between U. minor and U. pumila, e.g., [3,42,45,46,52,76]. The lack of reproductive barriers between U. glabra and U. minor, their overlapping distribution, and their genetic proximity promote spontaneous hybridization [27,47,51,77]. Based on Bayesian clustering results, putative hybrids constituted 28.3% (30 out of 106 individuals), which was comparable with results of wych elm in Belgium, where Cox et al. [74] identified significant introgression (46%) in natural populations. In contrast, Martin del Puerto et al. [47] identified significantly less (13%) putative hybrids among elms in the Iberian Peninsula. The differences in the higher percentage of putative hybrids in Belgium may be explained by the higher abundance of both species U. glabra and U. minor at the same altitude [74]. In contrast, U. glabra and U. minor distributions are partly overlapping, partly divided, and limited to a certain altitude in the Iberian Peninsula, so hybridization is relatively weaker [47]. Thus, our results showing 28.3% of putative hybrids were more in line with the study of elms in Belgium [74].
Finally, based on results from NewHybrids, 31 individuals out of 106 in total were identified as possible F2 hybrids (with p ≥ 0.5). However, the threshold to consider individuals as putative hybrids in our study was lowered to p ≥ 0.5 in contrast to other studies, e.g., [74,78]. Additionally, in contrast to other studies e.g., [42,47], we did not identify a significant increase of genetic diversity due to hybridization between U. glabra and U. minor (data not shown).
The initial number of pure reference samples per parental species was not high in our study (e.g., 28 of U. glabra and 40 of U. minor and U. minor subsp. minor), which may have an effect on the genetic assignment into hybrid groups. However, a combination of molecular data and morphological characteristics showed high potential and could help the classification of the genus Ulmus, especially in hybrid individuals, e.g., [42,[45][46][47]52,76]. Our study showed the differentiation of U. laevis from other Ulmus species, as well as U. minor and U. glabra hybrid groups, based on both leaf morphology and genetic characteristics. In contrast, the classification of trees within the U. minor groups based on leaf morphology did not correspond well to the clustering based on eight microsatellite markers using genetic clustering. Thus, in our study, the morphological leaf characteristics typically used to identify elms in the field did not reliably distinguish taxonomy within the U. minor complex. Similar conclusions that leaf morphology-based clustering is not always congruent with clustering based on genetic markers in Ulmus spp.were reported by several authors, e.g., [45,46,52,74] Therefore, further autochthonous Ulmus species conservation measures should take a more detailed genetic examination within the U. minor complex into account to enable better species differentiation, which is the basis for in situ and ex situ conservation.
Our genetic diversity analysis showed markedly stronger deviation from random mating in the putative hybrid than in the pure species groups. This indicated differentially mating groups identifiable by a divergent leaf morphology in natural elm forests. As leaf morphology suggests, these groups could be the spontaneous hybrid formations within a mixture of elm taxonomic groups, especially with the U. glabra × U. minor hybrid group being separated by genetic clustering. There could be multiple generations of reciprocal hybrid mating that could obscure the leaf morphology-based taxonomic identification, especially within the U. minor species complex where the genetic data indicating stronger mating among the Ulumus spp. with no marked deviation from random mating.

Conclusions
In conclusion, the genetic and leaf morphology analyses of putative elm hybrid swamps indicate a low genetic exchange between U. laevis and the other species groups in the Ulmus species complex sensu lato. This also indicates a low probability of the contribution of U. laevis in forming spontaneous hybrids among the Ulmus species. However, the genetic data do support leaf morphology-based identification of U. glabra (female) × U. minor (male) spontaneous hybrids. There is a strong genome sharing among U. minor and U. minor spp. minor species. This supports a unified taxonomic reference for U. minor and U. minor spp. minor. For gene conservation, we suggest considering separate gene conservation units selected based on leaf and stem morphology for U. laevis, U. glabra, U. glabra × U. minor, and the U. minor species complex.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/f12060653/s1. Figure S1: The PCA ordination plots of elm trees given separately for speciesspecific leaf morphology traits. The symbol size indicates the relative size of the morphology traits in the entity. The minimum value (zero) is shown on an overlay as the smallest size for that symbol. Abbreviations and color of the triangles are shown in Figure 3; Figure S2: Distribution of genetic diversity among six sampled Ulmus spp. groups (Na-Mean no. of Different Alleles; Ne-Mean no. of Effective Alleles; Ar-Mean allelic richness (based on min. sample size of 12 diploid individuals.), Npriv-No. of Private Alleles; He-Expected Heterozygosity) (GenAlEx 6.5 [61]); Populations abbreviations in Table 4; Figure S3: Private alleles distribution among the studied six Ulmus spp. groups (133 individuals) (R package poppr [62]; Figure S4: Missing data among six target Ulmus spp. groups and among eight nSSR loci (R package poppr [62]); Figure S5: The results of Bayesian clustering (soft. STRUCTURE2.3.4 [67]) on the most likely number of genetic clusters within the studied six Ulmus spp. groups, indicated by the highest delta K value at K = 2 (STRUCTURE HARVESTER soft. [69]); Figure S6: The results of Bayesian clustering (soft. STRUCTURE2.3.4 [67]) on the most likely number of genetic clusters within the studied five Ulmus spp. groups (106 individuals), indicated by the highest delta K value at K = 2 (STRUCTURE HARVESTER soft. [69]); Table S1: List of nuclear microsatellite markers (nSSR's) used in our study;  [67]) on the most likely number of genetic clusters within the studied populations, indicated by the highest delta K value (STRUCTURE HARVESTER soft. [69]). Funding: The material presented in the article was collected during the long-term LAMMC research program "Sustainable Forestry and Global Change".

Data Availability Statement:
The data is available upon request by e-mail to the corresponding author.