Genetic Diversity in Relict and Fragmented Populations of Ulmus glabra Hudson in the Central System of the Iberian Peninsula

Ulmus glabra Hudson, or Wych elm, occurs as fragmented and relict natural populations in the Central System, which acts as a refugium in the Iberian Peninsula. Considering the importance of the Central System populations of U. glabra, the main objective was to assess their genetic diversity using nuclear microsatellite markers. A total of 360 different genotypes were detected in the 427 U. glabra individuals analyzed. Wych elm populations showed a highly significant genetic differentiation (24%; p = 0.0001). Of the 22 populations studied, population of Rozas de Puerto Real (ROZ) showed the highest values of effective number of alleles (2.803), mean Shannon’s diversity (1.047) and expected heterozygosity (0.590). Populations of ROZ and Mombeltrán (MOM) showed the highest values of observed heterozygosity (0.838 and 0.709, respectively), and highly negative values for inbreeding coefficient (−0.412 and −0.575, respectively). Also, most of putative hybrids (50 of 55) were observed in these two populations. Demographic analysis revealed signals for recent (four populations) and ancestral (fifteen populations) bottlenecks. Fragmented populations with diminishing number of individuals, along with anthropogenic intervention and Dutch elm disease (DED), are the main threats to U. glabra populations. From a future perspective, the information generated can be considered in the formulation of conservation strategies for U. glabra populations in the Central System.


Introduction
Wych elm (Ulmus glabra Hudson) is a deciduous tree distributed mainly in the Eurosiberian region. U. glabra is spread over central and northern Europe, reaching up to the Ural Mountains and further extending towards eastern and southern Europe (i.e., Mediterranean region) [1]. The Iberian Peninsula is the southwestern limit for Wych elm [2]. The continuity and abundance of U. glabra populations in the peninsula is greater in the northern strip and Pyrenean-Cantabric axis. The most meridional populations are separated by hundreds of kilometers [2]. Some of these populations belong to the Central System ( Figure 1). Due to its location, orientation and altitude, this Central System is System is important in modulating the general climate and acting as an island. These climatic and geographical characteristics have enabled the persistence of some areas with Eurosiberian conditions surrounded by zones with Mediterranean climate [3]. The persistence of these areas acting as refugia with a favorable climate, has enabled the survival of the Eurosiberian species till the present [4,5].  Table 1 for codes).
The presence of the genus Ulmus in the Iberian Peninsula dates back to the Miocene and also has been reported throughout the Pleistocene and Holocene [6]. Paleobotanic and paleoclimatic studies confirm a temperature increase over the Late Glacial and, also, in the Holocene [7]. Therefore, the Mediterranean climate has gradually expanded to the Iberian Peninsula. Consequently, Eurosiberian vegetation has been reduced both latitudinally and altitudinally [4,8,9]. In addition to the climatic factors, the development of human-related activities has drastically reduced the woody vegetation in these mountains [10,11], thus affecting the genus Ulmus [12].
Besides landscape, climatic and anthropogenic factors, natural interspecific hybridization and introgression among relatives of the genus Ulmus have had significant consequences on the genetic diversity, genetic structure and conservation of their natural populations [13][14][15]. In the Iberian Peninsula, three Ulmus species are considered natives: U. glabra, U. minor Mill. and U. laevis Pall. [16,17]. U. minor and U. glabra are genetically most close [18], so hybridization and introgression between the two species have been described in natural populations where Field elms and Wych elms coexist [15,19]. Both these Ulmus species are monoecious, self incompatible, wind pollinated and their seed dispersal is also by wind [1].  Table 1 for codes).
The presence of the genus Ulmus in the Iberian Peninsula dates back to the Miocene and also has been reported throughout the Pleistocene and Holocene [6]. Paleobotanic and paleoclimatic studies confirm a temperature increase over the Late Glacial and, also, in the Holocene [7]. Therefore, the Mediterranean climate has gradually expanded to the Iberian Peninsula. Consequently, Eurosiberian vegetation has been reduced both latitudinally and altitudinally [4,8,9]. In addition to the climatic factors, the development of human-related activities has drastically reduced the woody vegetation in these mountains [10,11], thus affecting the genus Ulmus [12].
Besides landscape, climatic and anthropogenic factors, natural interspecific hybridization and introgression among relatives of the genus Ulmus have had significant consequences on the genetic diversity, genetic structure and conservation of their natural populations [13][14][15]. In the Iberian Peninsula, three Ulmus species are considered natives: U. glabra, U. minor Mill. and U. laevis Pall. [16,17]. U. minor and U. glabra are genetically most close [18], so hybridization and introgression between the two species have been described in natural populations where Field elms and Wych elms coexist [15,19]. Both these Ulmus species are monoecious, self incompatible, wind pollinated and their seed dispersal is also by wind [1]. Thus, the available data about vegetation evolution together with the current distribution of U. glabra in the Iberian Peninsula emphasize the relict and peripheral distribution of Wych elm populations in the Central System. According to Rossignoli and Génova [2], these populations are not only relict but also have fragmented distribution. Collin et al. [20] have emphasized the need for special attention to marginal populations of U. glabra because of their rarity and adaptability to unique geographical conditions. These populations are also important for the natural heritage preservation, particularly in view of the fact that major threats to their survival were found [2]. The scattered distribution, small population size and impact of human activities, can increase the fragmentation of populations along with genetic erosion, mainly towards the southern limits of the species distribution, causing them to be more susceptible to stochastic, genetic, and demographic events [21].
Moreover, in the last decades, Dutch elm disease (DED) has reduced the number of trees in Wych elm populations [2,16]. Fifteen populations of U. glabra from the Central System have been previously described [2] as having low numbers of individuals in general and with a high contribution of young individuals combined with re-sprouting stumps. For the present study, these 15 populations and an additional seven populations have been included. Hence, 22 populations of U. glabra, located along a 600-kilometer stretch of the large mountainous region, which includes the Spanish provinces of Guadalajara, Segovia, Madrid, Ávila and Cáceres, and one population found in the Serra da Estrela in Portugal (Figure 1), have been analyzed in the present investigation.
For understanding and establishing adequate conservation management strategies in the natural populations of U. glabra in the Central System of Iberian Peninsula, characterization of genetic diversity is the first necessary step. To study genetic variability in populations of elm species, molecular markers such as Random Amplified Polymorphic DNAs (RAPDs) [22][23][24], Inter Simple Sequence Repeats (ISSRs) [23], and Amplified Fragment Length Polymorphisms (AFLPs) [15,16,25] have been used. AFLP markers have also been used to investigate interspecific hybridization between U. minor and U. glabra in Flanders [15].
In addition, numerous microsatellite markers have been developed for different Ulmus species, such as U. laevis [26], U. minor [27] and U. rubra Muhl. [28], and used for testing of cross-amplification in other elm species, including U. glabra. Many U. laevis populations located along Europe have been studied using these markers [17,29,30]. Genetic diversity of a few localized natural populations of U. minor (e.g., in Balearic Islands, [31]; in the Netherlands, [32]; in Croatia, [33]) and U. glabra (e.g., in Denmark, [34]) have also been studied using microsatellite markers. The most extensive analysis of European populations of these two species together, using restriction fragment-length polymorphisms in chloroplast DNA, AFLPs and nuclear microsatellites [16], included only five populations of Wych elm from Spain, of which one was from the Central System. Moreover, there are some studies using microsatellite markers focused on detecting hybridization and introgression between elm species, e.g., between invasive U. pumila L. and native U. rubra in USA [13,35] or native U. minor in Italy [14,36].
Keeping in view the importance of Central System populations of U. glabra in the Iberian Peninsula and simultaneous lack of studies on genetic diversity, interspecific hybridization and introgression, the major objectives of the present investigation were (i) to assess the level of the genetic diversity of U. glabra populations in the Central System; (ii) to analyze possible demographic processes (e.g., evidence of bottleneck events) they have undergone; and (iii) to detect possible introgression between U. glabra and U. minor, which can significantly influence the genetic diversity and genetic structure in Wych elm populations.

Sampling of Plant Material
Leaf material of 427 individuals from 22 populations of U. glabra in the Central System was sampled for the present study (Table 1, Figure 1). Whenever possible, leaf material was collected from clearly-separated individuals to reduce the risk of sampling ramets. In addition, fifteen individuals of U. minor were sampled in three different Central System locations (Table 1). Also, since U. glabra populations have a fragmented distribution in the Central System and the number of estimated individuals/population (Nt) was extremely varied (see Table 1), sampling was done so as to reflect the distribution of the species in the Central System. Hence, a representative number of individuals (approximately between 20 and 60) which were mostly at a uniform separation within the population, were collected. For populations that had less than 20 estimated individuals/population (see Table 1), collection of all possible individuals was done, after taking care to avoid the collection of ramets. In such cases the number of individuals per population was less (less than 20). Leaf material of all individuals was stored at −80 • C until DNA extraction.

DNA Extraction and Microsatellite Genotyping
Total genomic DNA was extracted from frozen leaf material following the protocol supplied in the "NucleoSpin ® Plant II Kit" (Macherey-Nagel GmbH & Co. KG, Düren, Germany). Extracted DNA was quantified by visual comparison with known concentrations of lambda DNA on 1.2% agarose gels, and a working solution of DNA (approx. 10 ng/µL) was made.
One primer of each pair was labeled with a fluorescent dye

Genetic Diversity
Genetic diversity parameters were calculated based on data from eleven microsatellite loci analyzed in 427 Wych elm samples from 22 populations. For each locus: allele frequencies, total number of alleles (N A ) and effective number of alleles (N E ), number of different genotypes observed (G O ), observed heterozygosity (H O ) and expected heterozygosity (H E ), and Shannon's information index (I) were calculated, using GenAlEx software version 6.5 [37,38]. Polymorphic information content per locus (PIC) was also calculated according to Bostein et al. [39]. Furthermore, the discrimination power (D; [40,41]) of microsatellite markers was estimated. This parameter is an estimate of the probability that two randomly-sampled trees could be distinguished by their microsatellite profiles, and was calculated for each locus as D = 1 − C, where C is the probability of coincidence, i.e., two samples match by chance at one locus (C = ΣP i 2 , and P i is the frequency of different genotypes observed at that locus). The discrimination power for all loci combined (m = 11) was calculated as D T = 1 − C T , where C T = ΠC m , and represents the probability of coincidence cumulative for all loci.
Considering the overall microsatellite loci analyzed, for each population the percentage of polymorphic loci (P) and G O were calculated. Also, other parameters, such as the mean over all loci of N E , I, H O , H E and inbreeding coefficient (F), were also calculated using the GenAlEx 6.5 software [37,38].
Relationships among 19 U. glabra populations (with at least six sampled individuals) were obtained using the UPGMA (unweighted pair group with arithmetic mean) method with NTSYS-pc software version 2.2 [42], based on Nei's standard genetic distances [43]. A co-phenetic matrix was derived from the genetic distance matrix to test the goodness-of-fit of the clusters by comparing the two matrices using the Mantel matrix correspondence test in the MxComp program of the NTSYS-pc package, using 10,000 random permutations. The robustness of the dendrogram was also tested using 10,000 bootstraps implemented with the Phylip software version 3.6 [44].
Additionally, in order to determine the significance of partitioning of genetic diversity among and within populations, an analysis of molecular variance (AMOVA; [45]) was conducted using the GenAlEx 6.5 software. Levels of significance of variance component estimates were computed by non-parametric permutational procedures using 10,000 random permutations. Genetic differentiation (F ST ) was also calculated over all populations in a group and then among all pairs of populations, and their significance was tested by a permutation procedure using 10,000 permutations.

Demographic Analysis
The effective population size (N e ) is a critical parameter to study the evolution and conservation in natural populations. As estimates of N e can be affected by recent population size changes, the test for evidence of recent bottlenecks in U. glabra populations was performed using the T 2 statistic implemented in the program BOTTLENECK version 1.2 [46,47]. This statistic represents an average over loci of the deviation of the actual H E from the H E expected from the number of alleles in the population assuming mutation-drif equilibrium in a constant size population. The program tests for excess of heterozygosity using several mutation models including the stepwise mutation model (SMM), infinite allele model (IAM), and a two-phase model (TPM). For microsatellite data, the TPM is suggested to be the most appropriate mutation model [47][48][49]. Thus, the TPM model was applied assuming 95% single-step mutations and 5% multiple-step mutations [47], and a variance among multiple steps of 26. Positive values of T 2 indicate bottleneck, whereas negative values are consistent with recent population expansions. Also, considering the number of polymorphic loci used in the present study, the Wilcoxon's signed rank test was chosen to determine the T 2 significance as recommended by Piry et al. [47]. The test was run with 10,000 iterations.
Additionally, the M-ratio values (calculated as a ratio of the number of alleles to the range in allele size at a microsatellite locus [50]) were calculated using Arlequin software version 3.5 [51]. M-ratio value was used to detect reductions in effective population size. In comparison to the BOTTLENECK methods that are tailored to detect recent bottlenecks, M-ratio reflects a population size decline over a longer timescale and may be more powerful for detecting ancestral bottlenecks [48,50]. Based on empirical evidence, Garza and Williamson [50] indicated that in any population study where seven or more loci are analyzed and if the population shows M-ratio below 0.68, then it may be assumed that they have undergone a genetic bottleneck [17]. Both T 2 bottleneck statistic and M-ratio method used to detect bottleneck events were carried out in only 15 U. glabra populations that were comprised of at least 10 sampled individuals.
Finally, the effective population size for each of the 10 U. glabra populations with at least 20 sampled individuals, was calculated using two different single-sample estimation methods: the linkage-disequilibrium information and approximate Bayesian computation. LDNe software version 1.31 [52] was used, which employs linkage-disequilibrium information among alleles at different loci caused by genetic drift in finite populations. This method corrects for biases associated with small sample sizes [52,53]. N e was estimated assuming random mating and excluding alleles with frequencies less than 0.05. This critical value has been shown to produce little bias from rare alleles while maintaining moderate precision [54]. The 95% confidence intervals for N e were obtained using a Jackknife method as recommended in Waples and Do [52].
Approximate Bayesian computations to estimate variance N e from summary statistics that are related to N e were performed using ONeSAMP software version 1.1 [55]. The user defines prior estimates of N e (lower and upper bounds) and ONeSAMP generates 50,000 simulated populations drawn randomly from the distribution of N e priors. Summary statistics are calculated and compared to the actual data and similar summary statistics are retained for use in generating an estimate of N e using weighted local regression [55]. Several lower and upper bounds were tested on the prior estimates for N e per population. Finally, 2-7000 for two populations (STR and IRU; see Table 1 for population codes) and 2-200 for the other eight populations tested, were used to estimate their N e along with 95% confidence intervals.

Introgression between U. glabra and U. minor
To identify genetic groups within U. glabra populations and the relationships with U. minor genotypes, STRUCTURE software version 2.3.4 [56] was used. This Bayesian approach uses no a priori classification and assigns samples to K genetic clusters based on the allele frequencies at each locus. A first analysis was carried out on 427 U. glabra individuals using data from the eleven microsatellite markers, and the range of possible groups (K) tested was from 1 to 25, i.e., the total number of sampled populations plus three [57]. Then, to know the relationships between U. glabra and U. minor genotypes, a second analysis was performed by adding the 15 U. minor samples, using only data from the nine microsatellite loci that were amplified in both species, and wherein the number of groups tested ranged from 1 to 10 considering the results of previous analysis.
The estimate of the most likely number of genetic groups was performed following the procedure of Evanno et al. [57], which proposed the ad hoc statistic ∆K. Program settings used were the admixture ancestry and correlated allele frequencies models. The degree of admixture alpha was inferred from the data, and lambda, the parameter of the distribution of allelic frequencies, was set to 1 [57,58]. The program was run 20 independent times for each K value. In each run, a burn-in period of 10,000 iterations, and 100,000 post-burning MCMC (Markov chain Monte Carlo) simulations, were carried out. Finally, among the 20 runs performed for the optimal K value estimated, the run with the least negative log-likelihood value [59] was used to obtain the membership coefficients (q) for each individual in each of the K inferred groups.
In the second STRUCTURE analysis, admixed individuals could represent putative hybrids between U. glabra and U. minor, so an analytical tool more focused on hybrid detection was also applied. Bayesian statistical methods provided in the program NewHybrids version 1.1 beta [60] were used to estimate posterior probabilities of whether each analyzed individual belonged to parental classes (U. glabra and U. minor) or hybrid categories (i.e., F 1 , F 2 and backcrosses F 1 ). This was performed using prior information about individual assignments corresponding to genetic clusters previously detected by STRUCTURE analysis. The program was carried out using the default parameters for the six genotype class frequencies and ten independent runs with a burn-in period of 100,000 iterations followed by 600,000 MCMC simulations.

Genetic Diversity
An analysis of 11 microsatellite loci in 442 elm individuals (427 U. glabra and 15 U. minor) revealed a total of 101 alleles. Of the 101 alleles, 92 were detected in U. glabra (see Tables S1 and S2). Frequencies of 50 out of 92 alleles were below 0.05, whereas five alleles showed frequencies above 0.5 (Table S2).
Results of the genetic diversity analysis in U. glabra per locus are presented in Table 2 Table 2). The maximum and minimum values obtained for D were also for Ulmi1-165 (0.974) and Ulm2 (0.269) respectively, with a cumulative probability for all loci close to 1 ( Table 2). In general, microsatellite loci developed in U. minor (Ulmi1-SSR series) showed the highest values for N A , N E , G O , heterozygosity, PIC and D, and the lowest values in Ulm-SSR series, developed in U. laevis.  Microsatellite-based genotyping for the 11 combined SSR loci in 427 individuals of U. glabra revealed a high number (360) of distinct genotypes (Table S1). Only 25 of 360 genotypes were common to two or more individuals. The three Ulmi1-SSR loci alone showed a discrimination power of 0.99984. They could distinguish 301 different genotypes, which is about 84% of the total genotypes detected using the 11 SSR loci.
The results of genetic diversity analysis in populations of U. glabra are presented in Table 3.  (Table 3), when populations with sample size below 5 (CNT, BOC, and PAU) and population with one unique genotype (i.e., CER) were excluded from the analyses. Populations of ROZ and MOM showed the highest values of observed heterozygosity (H O = 0.838 and 0.709, respectively), which corresponded to highly negative values for inbreeding coefficient (F = −0.412 and −0.575, respectively; Table 3). In the present investigation, ten unique or private alleles were detected in seven of the 22 populations (NAV, MON, RAS, IRU, CAS, CVJ and MOM; see Table S1). The abundance of three of the unique alleles (Ulm2-111, Ulm3-164, UR173a-185 found in CAS, MON and RAS, respectively) was >50% (Table S1).

Population Relationships
The UPGMA dendrogram shows the relationships among 19 U. glabra populations with more than five sampled individuals ( Figure 2). The Mantel test revealed a good and significant co-phenetic correlation (r = 0.87; p = 0.0001), which indicates a good fit to the cluster analysis [42]. The dendrogram shows the populations grouped in two main clusters, separated at a genetic distance level of 0.48. Populations of ROZ and MOM group together in a cluster with a support of 66.5%, and in the other cluster are the remaining 17 populations with a support of 54.3% (Figure 2). Within the second cluster, CER, RAS and MON populations are distanced from the remaining 14 populations (Figure 2). On the other hand, the western populations (VIL, STR and BEN) group together, and especially VIL and STR are located in the same subgroup and are separated from the rest with a support of 72.4% (Figure 2).  Table 1 for codes), based on 11 microsatellite loci data using Nei's genetic distance. The numbers at nodes are support (>50%) by bootstrap analysis (10,000 replicates).
Results of AMOVA analysis revealed a highly significant genetic differentiation among the 19 U. glabra populations (24%; p = 0.0001). A major portion of variation was harbored within populations (76%). Genetic differentiation values among all population pairs were also calculated (Table S3), and all of them were highly significant (p < 0.0005). The highest FST value was observed between RAS and MOM populations (43.6%) and the lowest value was between CAS-CAN and RZA-CAN population pairs (8.5% in both cases).

Demographic Analysis
The results of demographic analysis are shown in Tables 4 and 5. T2 statistic using the TPM model implemented in BOTTLENECK program showed a significant (p < 0.05) excess of heterozygosity (positive values of T2) in four (NAV, ROZ, TIE and MOM) of the 15 populations studied ( Table 4), suggesting that these populations could have suffered a recent bottleneck. On the other hand, values of M-ratio (a recommended statistical parameter to detect a bottleneck lasting several generations followed by a demographic recovery) revealed bottleneck marks (based on M-ratio < 0.68; [50]) in all the analyzed populations ( Table 4).
The effective population size values for the ten Wych elm populations that were analyzed are shown in Table 5. The Ne was less than 20 in eight and seven of the ten analyzed populations using ONeSAMP and LDNe programs, respectively. The lowest Ne values were observed in ROZ population using both programs (Table 5).  Table 1 for codes), based on 11 microsatellite loci data using Nei's genetic distance. The numbers at nodes are support (>50%) by bootstrap analysis (10,000 replicates).
Results of AMOVA analysis revealed a highly significant genetic differentiation among the 19 U. glabra populations (24%; p = 0.0001). A major portion of variation was harbored within populations (76%). Genetic differentiation values among all population pairs were also calculated (Table S3), and all of them were highly significant (p < 0.0005). The highest F ST value was observed between RAS and MOM populations (43.6%) and the lowest value was between CAS-CAN and RZA-CAN population pairs (8.5% in both cases).

Demographic Analysis
The results of demographic analysis are shown in Tables 4 and 5. T 2 statistic using the TPM model implemented in BOTTLENECK program showed a significant (p < 0.05) excess of heterozygosity (positive values of T 2 ) in four (NAV, ROZ, TIE and MOM) of the 15 populations studied ( Table 4), suggesting that these populations could have suffered a recent bottleneck. On the other hand, values of M-ratio (a recommended statistical parameter to detect a bottleneck lasting several generations followed by a demographic recovery) revealed bottleneck marks (based on M-ratio < 0.68; [50]) in all the analyzed populations ( Table 4).
The effective population size values for the ten Wych elm populations that were analyzed are shown in Table 5. The N e was less than 20 in eight and seven of the ten analyzed populations using ONeSAMP and LDNe programs, respectively. The lowest N e values were observed in ROZ population using both programs (Table 5).

Introgression between U. glabra and U. minor
In order to analyze the existence of a possible introgression between U. glabra and U. minor, genetic data from nine microsatellite loci that amplified in both species were used. Seventy six alleles in 427 U. glabra individuals from 22 populations, and 28 alleles in 15 U. minor individuals from three locations in the Central System, were detected (see Table S2). Nineteen of the 28 alleles detected in U. minor were common to U. glabra also. Two of these alleles (Ulm3-182 and UR123-252) showed high frequency (≥0.5) in both species (Table S2). All U. glabra samples shared at least two alleles with individuals of U. minor (see Table S1). STRUCTURE analysis carried out on 427 U. glabra individuals using data from the 11 microsatellite markers, showed the modal value of the distribution of ∆K, i.e., the optimal number of groups estimated, for K = 2. Similar results were obtained when a second analysis was performed using data only from the nine microsatellite loci in 442 samples, 427 U. glabra and 15 U. minor ( Figure S1).
With respect to 427 U. glabra individuals, STRUCTURE analysis revealed the existence of two clearly-defined genetic clusters (A and B), based on coefficients q ≥ 0.9 for one of the two groups [61]. Thus, 368 individuals were assigned to genetic group A, 47 to group B, and 12 were considered as admixture samples, because they showed admixture values below 0.9 (Table 6 and Table S4). The 15 U. minor individuals were clustered with the 47 U. glabra from group B (Table S4). Table 6. Number and percentage of individuals assigned to clusters (A and B) or considered as admixture samples as defined by STRUCTURE analysis in each Ulmus glabra population studied. See Table 1 for population codes.

Population
Ns For performing NewHybrids analysis, the 368 U. glabra individuals which had clustered into group A (by STRUCTURE analysis) and the 15 U. minor individuals which had clustered into group B (by STRUCTURE analysis), were treated as reference samples for U. glabra and U. minor respectively. Similar to STRUCTURE analysis, all individuals with a posterior assignment probability p ≥ 0.9 [61], were placed in one of the six predefined genotypic classes (Table S4). Based on NewHybrids analysis, 372 individuals were assigned to parental class of U. glabra, 22 (all from ROZ population and with the same genotype; see Table S1) together with 15 U. minor samples were assigned to U. minor class, and 26 to F 1 hybrid class (16 from ROZ population, 7 from MOM, 2 from CVJ and 1 from IRU). The remaining seven samples (5 from ROZ population, 1 from IRU and 1 from CAN) showed values of posterior assignment probability lower than 0.9 and were not assigned to any parental or hybrid class (Table S4). The 55 individuals (22 of ROZ + 26 of F 1 hybrid class + 7 unassigned) which were not classified into U. glabra parental class following NewHybrids analysis, showed the highest number of heterozygotic loci (7 to 10; see Table S1) and these individuals also shared 5 to 10 alleles with U. minor (see Table S1).

Genetic Diversity and Population Relationship
An analysis of microsatellite diversity values (Ns = 427, H O = 0.477 and H E = 0.566; see Table 2) of populations of U. glabra from the Central System yielded lower values than a similar study in the Suserup Forest of Denmark (Ns = 29, H O = 0.759 and H E = 0.774; [34]). However, as each study analyses a different number of samples, one should exercise caution while comparing the results. In the present analysis as well as in the study by Nielsen and Kjaer [34], Ulmi1-165 was the most polymorphic and informative locus, suggesting that Ulmi1-165 is an appropriate locus to study genetic diversity of U. glabra.
Populations with the highest number of individuals (IRU, Nt > 6500; STR, Nt > 5000) showed high diversity indexes (Table 3). However, this trend cannot be generalized (e.g., populations MOM and CAN which have small population size also show high diversity values). The location of IRU and STR in protected areas (Natural Reserve of Valley of Iruelas and ZEPA-Zona de Especial Protección para las Aves-of Sierra de Gata, respectively) has probably helped to maintain the high number of individuals. Also, most of the IRU and STR individuals are young with very limited rates of sexual reproduction [2]. According to Rossignoli and Génova [2], in the last few decades, DED has affected the populations of Wych elm from the Central System. As a consequence of this, many ramets are formed by sprouting from the stumps of the old trees, resulting in an increase of young individuals in a population. Thus, DED may be the reason behind the high number of young individuals, as the older elms are preferred by bark beetles and therefore perish due to disease [62].
Population CER had only six young individuals in nature, located in an area of approximately 20 m 2 , which showed the same genotype. This probably indicates that they may be sprouts of one adult tree, which does not exist anymore. In general, the presence of sprouting was detected in all Ulmus populations analyzed in this study. Wych elm does not show root sprouts but sprouts on trunks of young trees have been described by Cox et al. [15]. Vegetative reproduction by sprouting helps to diminish the effects of fragmentation and genetic drift, thus avoiding allelic loss [63] and contributing to clonal restoration of U. glabra populations on disturbed sites [34].
In general, signs of population isolation and low gene flow were detected throughout the Central System (F ST = 0.239) compared to that reported for some other wind-pollinated temperate forest species: F ST = 0.049 in Quercus garryana Douglas ex Hook [64], F ST = 0.052 in Populus euphratica Olivier [65], Φ ST = 0.013 in Fraxinus excelsior L. [66], and F ST = 0.074 in Castanea sativa Mill. [67]. In the present study, high and significant F ST values among population pairs indicated both high genetic differentiation and low average diversity within populations, probably due to small population size in most of them. Pollen flow up to 1000 m was estimated in U. glabra [34], so geographic distances between elm populations appear to have favored population differentiation observed in the Central System.
The UPGMA tree showed a main cluster of 17 populations and another of two populations. Populations MOM and ROZ group together and form a separate cluster, probably because both these populations are similar in that they contain a high number of putative hybrids (38/46 in ROZ and 7/10 in MOM) with U. minor and also low F ST amongst them. Within the majority cluster, MON and RAS populations are separated from the rest, probably because they are geographically isolated populations. Also, these populations showed high frequency of private alleles. The presence of private alleles has been associated with low gene flow and, therefore, population isolation [68]. The clustering results are related to geographic isolation and suggest a decrease in gene flow because of larger geographical distances [69]. In addition, the presence of BEN in the majority cluster suggests its genetic relatedness with other Iberian population in the Central System. This indicates U. glabra in Portugal (BEN) is a part of natural populations, as proposed by Monteiro-Henriquez et al. [70].

Demographic Analysis
Demographic analyses using Wilcoxon's test indicated bottlenecks in fewer populations (4 of 15; NAV, ROZ, TIE and MOM) compared to M-ratio values, where the 15 analyzed populations revealed bottlenecks. A similar trend of detection of bottlenecks by the two statistical methods was observed in U. laevis populations by Fuentes-Utrilla et al. [17]. Low M-values in the 15 populations indicated an ancestral and extended bottleneck [50]. Bottlenecks described in Spanish populations of U. laevis were attributed to Holocene expansion of more xeric habitats in the Iberian Peninsula [17] or glacial periods in Pleistocene [30]. It is probable that the same factors have affected the Wych elm, contributing to an ancestral and progressive loss of trees in the Central System, reducing the interpopulation gene flow (shown by high genetic differentiation) and leading to fragmented populations.
Four of the studied populations (NAV, ROZ, TIE and MOM) showed recent bottleneck signals as detected by T 2 statistic and Wilcoxon's test. In contrast to the M-value, T 2 bottleneck statistic is more effective in detecting bottlenecks in populations that have experienced a recent and severe reduction in the number of individuals [47].
In most populations of U. glabra from the Central System, clear signs of anthropogenic activity were observed, which can be attributed to different forestry land-uses (i.e., wood removals, grazing, etc.). In the case of NAV (now a managed pinewood forest), there is a history of successive and high wood removals [71]. Thus, it is possible that the Ulmus stands from NAV population also have been affected. According to the National Plan of Aerial Orthophotography (PNOA, © Instituto Geográfico Nacional; http://pnoa.ign.es/), in TIE population also a similar process of wood removal was observed in the 1970's. Therefore, it is possible that forest land-uses could be the reason for the existence of a recent bottleneck in TIE and NAV populations. In such cases, bottlenecks can further increase the interpopulation differentiation that existed prior to the anthropogenic intervention [72].
In the MOM and ROZ populations, only two genotypes in 10 samples and 13 genotypes in 46 samples, respectively, were observed (i.e., <30% of the sampled individuals). This may be a consequence of high rates of sprouting in these populations. The recent bottlenecks observed in these populations may be due to samples being the ramets of the same clone which the software recognized as a bottleneck. Also, repeated genotypes due to vegetative reproduction cause a reduction in N e of populations. Another effect of asexual reproduction is the absence of segregation, leading to maintenance of some of the allelic combinations in heterozygotes, which resulted in negative F values [73].
Furthermore, a very low N e (except in IRU), was observed in populations of the Central System. A reduction of N e of population is caused by factors leading to genetic drift, and such populations with low N e have a greater probability of extinction [74]. Wych elm is a self-incompatible species [34], and therefore even more susceptible to habitat fragmentation and reduction in population sizes. Thus, population fragmentation in Wych elm in the Central System is evident, with a scarce number of populations and a limited number of trees.

Introgression between U. glabra and U. minor
In the present study of populations of the Central System, introgression between U. glabra and U. minor was detected. Such introgression has also been observed in other European populations of Wych elm [15,75]. The genetic proximity between U. glabra and U. minor [18], and the lack of reproductive barriers between them [15], have favored natural hybridization. These two factors (i.e., genetic relatedness and absence of reproductive barriers) are probably responsible for the large number of alleles shared between both species, some of them in high frequency.
Putative hybrids constituted 13% (55 of 427 individuals) and were detected in five of 22 populations of the Central System (ROZ, IRU, CVJ, MOM and CAN). Although significant introgression was observed, it was less compared to other studies, e.g., in natural populations of Wych elm in northern Belgium (46%; [15]). The attributed reason is that in Belgium, U. glabra is more abundant and coexists with U. minor at the same altitude [15]. However, in the Central System the presence of U. glabra is limited and its populations can be found at an average altitude of 1150 m reaching <1800 m [76], while U. minor does not usually surpass 1200 m [77], thus introgression in the Central System is more restricted.
Besides co-existence of the two species at the same altitudes, introgression is influenced by anthropogenic activities as well [78][79][80]. Traditionally, U. minor is the more cultivated species and is frequently naturalized [76]. Populations of ROZ and MOM are located close to individuals of U. minor and are exposed to anthropogenic activities. These factors have probably resulted in maximum hybridization (38/46 individuals in ROZ and 7/10 individuals in MOM). The presence of large numbers of putative hybrids in ROZ and MOM is responsible for its high observed heterozygosity and estimated diversity values. Also, the inbreeding coefficient of the MOM and ROZ populations was highly negative, indicating excess of heterozygosity. Similar results were obtained by Zalapa et al. [13], where the extent of hybridization in naturalized populations of U. pumila influenced the genetic diversity.
Throughout the study area, 26 individuals were identified as possible F 1 hybrids by NewHybrids analysis. These samples were heterozygous for the majority of loci (9 or 10 of 11) analyzed, as expected in hybrid individuals. Zalapa et al. [13] have shown a similar increase of genetic diversity due to hybridization between U. pumila and U. rubra. In addition, the putative hybrid individuals detected in the present study shared numerous alleles with U. minor, some of which were present in high frequency in U. minor and uncommon in Wych elm, which indicated asymmetric introgression of U. minor into U. glabra. This could happen because at < 1200 m of altitude in the Central System, U. minor is in abundance compared to U. glabra, thus facilitating asymmetric gene flow. Due to large differences in abundance of parental species, asymmetric gene flow has been observed in other studies [81,82].
The NewHybrids program also classified 22 individuals of U. glabra from ROZ population, which presented the same genotype, into U. minor parental class. The lack of sufficient numbers of species-specific exclusive alleles can complicate the classification of individuals when using the NewHybrids software [60]. In the present investigation, the 15 samples of U. minor analyzed (which presented six different genotypes and eight exclusive alleles) could be insufficient to assign correctly some hybrid individuals. Morphologically, the 22 samples did not differ much from the rest of U. glabra samples collected. The morphological characteristics have been shown to be questionable in the classification of the genus Ulmus, especially in hybrid individuals [13,14,23,32,36]. These samples showed a high number of alleles shared with U. minor (9 of 18 possible alleles), and were heterozygous for most loci (10 of 11) analyzed. Consequently, these 22 individuals probably have undergone a high level of introgression with U. minor.
As in U. minor, root suckering was also observed in introgressed individuals from ROZ and MOM populations, which increases the number of introgressed samples with the same genotype in these populations. Similar observations were made by [15].
In the present study, despite extensive sampling, no F 2 individuals were detected using the NewHybrids program. Cox et al. [15] also did not find F 2 offspring of F 1 hybrids among elms sampled in northern Belgium. One possible explanation for failing to detect F 2 hybrids could be that an increase of homozygous incompatibility in F 2 individuals would lead to loss of fitness [83]. The other possibilities are the existence of pollen-stigma incompatibilities in F 1 hybrid, which would contain S-alleles from the two parental sources [35], or because of outbreeding depression [84].

Conclusions and Future Prospects
Populations of ROZ and MOM showed lowest numbers of genotypes, high heterozygosity, lowest genetic differentiation between the two populations, and highly negative values of inbreeding coefficient. This points towards the existence of introgression and propagation through vegetative means (as indicated by high H O , low G O and highly negative F) and that the populations are genetically not very different (low F ST ). The estimated population size is much lower (ROZ = 46, MOM = 18) compared to IRU (>6500) and STR (>5000). This could probably be due to the location of the IRU and STR populations in protected and undisturbed areas. Populations IRU and STR not only maintained a high number of individuals but also represented a very high number of genotypes (>95% of the sampled individuals), indicating low frequency of clonal propagation. A distinct feature of ROZ and MOM populations was that maximum numbers of putative hybrids (with U. minor) were detected, which enhanced their genetic diversity indexes. Population of RAS on the other hand did not show any putative hybrids and revealed lowest diversity indexes. However, a high number of genotypes with private alleles was detected (>87% of the analyzed individuals) in RAS population compared to ROZ and MOM (<30% of the sampled individuals), and high genetic differentiation between RAS-ROZ and RAS-MOM population pairs was observed. Thus, the RAS population appears to be more isolated with respect to interpopulation gene flow. Such information can play an important role in the formulation of appropriate conservation strategies for Ulmus populations of the Central System.
In general, populations of U. glabra analyzed in the Central System showed recent and/or ancestral bottlenecks and low effective population size. As U. glabra is a self-incompatible species, has bottlenecks and reduced N e , the chances of population fragmentation and decrease in number of individuals, and probability of extinction increases. At this point, it is also important to note that existence of clonal propagation in U. glabra populations possibly contributes to reducing the effects of fragmentation and genetic drift and hence allelic loss. Also, existence of natural hybridization with U. minor influences genetic diversity and may elevate the number of alleles and effective alleles in U. glabra populations.
To summarize, the fragmented and relict natural populations of U. glabra of the Central System of the Iberian Peninsula are marked by moderate genetic diversity, high genetic differentiation and ancestral as well as extended bottlenecks. The small number of fragmented populations with diminishing number of individuals in most populations, along with anthropogenic intervention and the biotic stress of DED are the main threats to Wych elm in the Central System. Further, there is no continuity of U. glabra populations of the Central System with the populations of the rest of the peninsula and no strategies for conservation of these relict populations. This is especially true in the province of Avila, where most of the populations have aggregated. Therefore, the present information generated from this study (i.e., high genetic differentiation and lack of continuity of Central System populations with the rest of the peninsula) emphasizes the need for formulating regional conservation strategies to secure the future of the Central System populations.
Supplementary Materials: The following are available online at www.mdpi.com/1999-4907/8/5/143/s1, Figure S1: Estimation of the optimal number of genetic groups (K) based on the rate of change in the statistic ∆K between successive K values [57], for K = 1 to 25 (using data from 11 microsatellite markers analyzed in 427 U. glabra individuals; a) and K = 1 to 10 (using data from the nine microsatellite loci that were amplified in both species, 427 U. glabra individuals and 15 U. minor samples; b)); Table S1: Allele sizes in base pairs at each of eleven microsatellite loci analyzed in 427 U. glabra individuals from 22 populations and 15 U. minor individuals from three locations sampled in the Central System; Table S2: Allele sizes (in base pairs) and frequencies at each of eleven microsatellite loci analyzed in 427 U. glabra and 15 U. minor individuals sampled from the Central System; Table S3: Pairwise F ST (below the diagonal) and p (above the diagonal) values among 19 U. glabra populations (see Table 1 for codes) based on eleven microsatellite loci data using the AMOVA method. Significance was obtained on 10,000 random permutations; Table S4: Results of Structure and NewHybrids analysis obtained in 442 elm individuals (427 U. glabra and 15 U. minor) studied.