Genetic variation in Tertiary relics: The case of eastern‐Mediterranean Abies (Pinaceae)

Abstract The eastern‐Mediterranean Abies taxa, which include both widely distributed species and taxa with minuscule ranges, represent a good model to study the impacts of range size and fragmentation on the levels of genetic diversity and differentiation. To assess the patterns of genetic diversity and phylogenetic relationships among eastern‐Mediterranean Abies taxa, genetic variation was assessed by eight nuclear microsatellite loci in 52 populations of Abies taxa with a focus on those distributed in Turkey and the Caucasus. Both at the population and the taxon level, the subspecies or regional populations of Abies nordmanniana s.l. exhibited generally higher allelic richness, private allelic richness, and expected heterozygosity compared with Abies cilicica s.l. Results of both the structure analysis and distance‐based approaches showed a strong differentiation of the two A. cilicica subspecies from the rest as well as from each other, whereas the subspecies of A. nordmanniana were distinct but less differentiated. ABC simulations were run for a set of scenarios of phylogeny and past demographic changes. For A. ×olcayana, the simulation gave a poor support for the hypothesis of being a taxon resulting from a past hybridization, the same is true for Abies equi‐trojani: both they represent evolutionary branches of Abies bornmuelleriana.

The main aim of this work was assessing the patterns of genetic diversity and phylogenetic relationships among eastern-Mediterranean Abies taxa with a focus on A. nordmanniana s.l. and A. cilicica s.l. with a dense coverage that would allow better delineation of distribution ranges of each taxon, as they are disputed in some cases. In order to achieve this, populations covering major parts of the distribution ranges of A. nordmanniana and A. cilicica in Asia Minor and eastern Caucasus were sampled. Populations of the closely related European species A. alba and A. cephalonica at the eastern limit of their distribution were also included in order to set the results into larger geographical context of the eastern-Mediterranean area. Some taxa that we studied are not recognized by the state-of-the-art checklists and there is none or few information on their genetic structure, but F I G U R E 1 Scheme of speciation sequence in eastern-Mediterranean Abies based on Linares (2011) with modification they can be clearly geographically delineated and their names are widely used by local botanists; this is why they were of interest. The divergence of genetic lineages leading to recent A. cilicica, A. cephalonica, and A. alba occurred in the Miocene or even earlier (Linares, 2011;Palamarev, 1989). Therefore, we considered the reconstruction of the complete demographic history of the studied species complex beyond the reach of our possibilities with the type of data and the number of loci available. Instead, we focused on taxa within A. nordmanniana s.l. and A. cilicica s.l. In particular, we addressed the following questions: • Is A. × olcayana a hybrid between A. equi-trojani and A. bornmuelleriana (as suggested by Ata & Merev, 1987)?
• What is the realistic phylogenetic scenario for taxa or regional populations within A. nordmanniana s.l. and A. cilicica s.l.?

| Plant material
A total number of 1,529 individuals from 52 indigenous populations were sampled; mostly in managed, naturally regenerated stands, less frequently in national parks, and nature reserves (Table S1;

| DNA extraction and genetic analysis
DNA extraction protocol according to Doyle and Doyle (1987) was used to extract DNA from approximately 50 mg of silica gel dried needles. Eleven nuclear microsatellite loci previously identified in various Abies species were used SF324, SF333, SFb04, SF1 (Cremer et al., 2006), ABF18 (Saito, Lian, Hogetsu, & Ide, 2005), NFF7, NFF3, NFH15, NFH3 (Hansen, Vendramin, Sebastiani, & Edwards, 2005), AB15 (Rasmussen, Andersen, Frauenfelder, & Kollmann, 2008), and AfSI16 (Josserand et al., 2006). Two multiplex and three singleplex reactions were done. The multiplex reactions were set up as 5 μl Amplification programs consisted of initial step at 94°C for 15 min., followed by 35 cycles of 30 s denaturation at 94°C, 90 s annealing at 58°C, and 90 s elongation at 72°C. A final elongation step was performed at 60°C for 20 min. 0.9 μl of multiplex A product with 0.1 μl F I G U R E 2 Results of the whole-dataset Structure analysis superimposed over the map of the eastern-Mediterranean. Charts represent inferred membership proportions of the studied populations. Distribution ranges of individual taxa (www.euforgen.org/species) are displayed in colors corresponding to the predominant Structure cluster of size standard GeneScan 500 LIZ (Applied Biosystems, Foster City, California) and 9 μl of formamide was loaded as a first batch for separation on capillary analyzer ABI 3130 (Applied Biosystems). Second batch consisted of the singleplex products (0.6 μl of each) together with 0.8 μl of multiplex B product, 0.1 μl of size standard, and 7.3 μl of formamide. GeneMapper 4.0 (Applied Biosystems) was used to analyze the raw data and produce genotypes.
Subsequently, the loci SFb04 and NFF7 were removed from further analyses. The locus SF324 was also discarded because it contains the same repeat segment as SF1.
To assess the patterns of gene diversity and allelic richness, expected heterozygosities of populations were calculated using arlequin and allelic richness and private allelic richness were calculated by the program hprare 1.1 (Kalinowski, 2005). Rarefaction to 22 and 140 gene copies was used to obtain allelic richness of individual populations and pooled groups (taxa), respectively. The populations Ann13, Aci4, and Aci5 were omitted from allelic richness analysis of individual populations because of low sample sizes, the same applies to A. alba and A. × olcayana for the analysis on the taxon level. Differences among taxa were tested by nonparametric Kruskal-Wallis test both on the subspecies and the species level; A. alba and A. cephalonica were excluded from this test.
The populations were checked for recent reduction in population sizes using the program Bottleneck v. 1.2.02 (Cornuet & Luikart, 1997), which relates the observed gene diversity with the heterozygosity predicted on the basis of the observed number of alleles under the assumption of mutation-drift equilibrium (Maruyama & Fuerst, 1985).
The distribution of the expected heterozygosity was obtained through simulating the coalescent process under two-phase model with the proportion of stepwise mutations of 95% and a variance among multistep mutations of 12%, as recommended by Piry, Luikart, and Cornuet (1999) for microsatellite loci. The probability of heterozygosity excess was tested using the Wilcoxon test (Cornuet & Luikart, 1997). In addition, the mode-shift approach applied in Bottleneck v. 1.2.02, based on search for transient distortions of allele frequency distributions induced by a reduction in population size, was used.
Several approaches to assess genetic differentiation of the studied populations were used. First approach was the program Structure 2.3 (Pritchard, Stephens, & Donnelly, 2000). It was run 16 times for each number of clusters ranging from 1 to 10 with 100,000 burn-in iterations and another 1,000,000 iterations without prior information on population structure. To determine the optimal number of clusters, the Structure harveSter script was used (version A.2 July 2014; Earl, 2012).
After assigning the populations to clusters inferred by Structure, the process was repeated second time for each cluster individually in order to find an additional genetic structure within each cluster. The same settings were used to run Structure, except for number of populations which ranged from 1 to 6, and again the optimal number of clusters was identified by Structure harveSter.
As a second approach, the Discriminant Analysis of Principal Components (DAPC) was used (Jombart, Devillard, & Balloux, 2010), implemented in the R package adeGenet 1.4.2 (Jombart, 2008). To more closely inspect the structure of the less distinct clusters, the most differentiated ones were gradually removed. In all the steps, a number of PCs corresponding to ~90% of the total variance were retained. Assignment of populations into groups followed the current taxonomic classification and geography.
Analysis of molecular variance (Excoffier, Smouse & Quattro, 1992) was carried out using the arlequin software separately for species and subspecies/taxa. In the latter case, the outcomes of the Structure analysis were considered, and the northern and southern populations of A. nordmanniana s.s. were considered as separate groups. The significance of variance components attributed to species/subspecies, populations, and individuals was tested using 99,999 random permutations.
Isolation by distance was tested according to Rousset (1997), re- (Slatkin, 1995) against logarithm of distance as recommended for a two-dimensional case. Significance of the relationship between genetic dissimilarity and distance was tested by Mantel test, the strength of IBD was quantified using a reduced major axis (RMA) regression, as applied in iBdwS v. 3.23 (Jensen, Bohonak, & Kelley, 2005). Confidence intervals of regression slopes were derived from 30,000 bootstraps over independent population pairs. Geographical distances were calculated by the program Geographical Distance Matrix Generator (Ersts, 2016). The Mantel tests were also performed separately for A. nordmanniana s.l., A. nordmanniana s.s. (whole subspecies and southern and northern populations separately), A. bornmuelleriana, A. cilicica s.l., A. cilicica s.s., and A. isaurica populations. The presence of barriers against gene flow was tested using the Monmonier's maximumdifferences algorithm implemented in the Barrier 2.2 software (Manni, Guérard, & Heyer, 2004), which is aimed at the identification of abrupt changes in genetic differences between pairs of populations in the geographical context, based on a network obtained by Delaunay triangulation. Again, Slatkin's linearized F ST (Slatkin, 1995) was used as a distance measure. Bootstrap support for the identified barriers was derived from 999 random resamplings. Barrier was run with a successively increasing number of groups (K), starting from 1, until new barriers with a bootstrap support of at least 50% were appearing.  (Cornuet et al., 2014). In all runs, 300,000 datasets were simulated for each scenario. A total of 1% of the simulated datasets most similar to the observed data were used to estimate the relative posterior. Simulations applied the Generalized Stepwise Mutation (GSM) model. As we have no information on the distribution of mutation rates of microsatellite loci in Abies, a mean of 2 × 10 −4 was used in analogy with Picea (Tsuda et al., 2016). One locus, showing excessively high frequency of single nucleotide insertions/deletions (AB15), was omitted from the analyses.
For the recalculation of the number of generations into a time scale, a reasonable estimate of generation time needed to be chosen.
In silver fir, isolated trees start flowering at 20 years of age and trees in a canopy at 60 years (Wolf, 2003), but this is for sure not the average age of offspring-producing trees in a natural forest ecosystem. Natural forests are characterized by a mosaic of developmental stages with cyclically changing standing stock and canopy closure. At the local culmination of standing stock, the height structure is typically one-layer with closed canopy, whereas the lack of light is strongly limiting for regeneration (in spite of abundant flowering and cone-bearing). In natural forests of silver fir (commonly growing in mixture with European beech, Norway spruce, and noble hardwoods), this stage usually occurs between 150 and 250 years of age of dominant trees (Korpeľ, 1995), before first canopy gaps enable survival of seedlings. A. nordmanniana s.l. grows in similar communities with similar stand structures (Mayer & Aksoy, 1986). On the other hand, firs in Asia Minor and the Caucasus may have experienced periods in their evolutionary history, when they spread by colonizing open areas; in that case, the generation turnover time must have been shorter. To account for this uncertainty and to use the analogy with A. alba, we counted with the generation turnover age of 150 ± 50 years (as a conservative choice) for the conversion. None of the populations exhibited a significant recent bottleneck, generally, heterozygosity deficiency was more prevalent than heterozygosity excess (Table S2). The mode-shift approach yielded identical outcome: even though a common graphical display of allele frequency distributions for 52 populations is difficult and the resulting picture is a bit confused, it is apparent that the distributions did not deviate from the L-shape typical for nonbottlenecked populations at mutation-drift equilibrium (Fig. S2). However, differentiation patterns indicate a strong fragmentation in some species: in spite of a range size comparable to the other studied Abies species and subspecies, A. cilicica subsp.

| RESULTS
isaurica was substantially more differentiated (F ST = 0.0562). In contrast, absence of differentiation within A. equi-trojani (F ST = −0.0027) indicates that the three studied samples were in fact drawn from the same population (Table 1).
Results of the Structure analysis strongly suggested the existence of six clusters (ΔK = 565; while the second highest ΔK was 11.7 for K = 4; Table S3, Fig. S3a). The first cluster comprised A. alba and Analysis of molecular variance showed that there is highly significant variation both at the species and the subspecies level, even though the respective variance components both represent around 5% of the total variation (   (Table S5, Fig. S6).
Estimates of effective population sizes (N e ) and divergence times are shown in

| Genetic variation
In contrast to continental Europe, the Mediterranean basin as a hotspot of species and genetic diversity (Cuttelod, García, Abdul Malak, Temple, & Katariya, 2008;Fady-Welterlen, 2005) was not as profoundly affected by Pleistocene glacial cycles as the more northern areas, and many taxa have persisted since the Tertiary (Palamarev, 1989). The degree of range fragmentation seems to be the main  (Alizoti, Fady, Prada, & Vendramin, 2011;Arbez, 1969;Mayer & Aksoy, 1986). The situation in the Caucasian part of the A. nordmanniana s.s. range (northern group) is not substantially different but the more rugged topography may contribute to isolation of local populations (Nakhutsrishvili, Zazanashvili, Batsatsashvili, & Montalvo Mancheno, 2015), and potentially to the loss of alleles by genetic drift, reflected in a lower allelic richness especially on the regional level.
The northern group is also located on the margin of the distribution of Mediterranean firs as such, the potential for genetic enrichment by gene flow from related species is very limited, also because of re-  Weir & Cockerham, 1984) prevent that the permanently small population size increases. In spite of this, A. equi-trojani retained relatively high allelic richness (at least at the population level) which, seen from the perspective of conservation, is a good news.
In accordance with the observation of Tayanç, Çengel, Kandemir, and Velioğlu (2012), levels of allelic richness and gene diversity in A. cilicica (notably in A. isaurica) were generally lower compared with the remaining taxa. The reasons may be associated with biogeography. The range of A. cilicica is more fragmented and the whole Taurus mountain range suffers from long-term forest degradation, especially at lower elevations (Gardner & Knees, 2013;Mayer & Aksoy, 1986).
Even though Sękiewicz et al. (2015) did not detect a difference in diversity between the two subspecies of A. cilicica, they reported a significant inbreeding and low effective population sizes in the studied populations.
Differences in biogeographical patterns are also reflected in the levels of differentiation. We observed a significant but very weak isolation by distance across all studied taxa. However, at the species/ subspecies level, A. bornmuelleriana was the sole taxon where significant isolation by distance was observed, which may, however, be respectively. Several of these clusters were also well-supported by the T A B L E 3 Posterior estimates of the parameters of the demographic inference based on the Approximate Bayesian Computation for the best-supported scenarios in different constellations of Mediterranean fir taxa or regional populations Posterior probabilities for the scenarios with the highest posterior probability of 10,000 sets of summary statistics most similar to the observed data through logistic regression. b Generations. c 1 generation = 150 ± 50 years.
genetic barriers found, especially in the case of A. cilicica and its subspecies. However, significant barriers also appeared within subspecies and regional populations, indicating strong effects of range fragmentation; this was mainly the case of A. isaurica growing in a harsh climate of southern Turkey. Obviously, geographical barriers and range discontinuities associated with the Pleistocene climatic fluctuations effectively prevented gene flow, which would otherwise contribute to homogenization of genetic structures, at least at neutral loci.

| Phylogeny
Even though the set of analyzed nSSR loci offers a quite poor coverage of the genome, the combination of distance-based (neighbornet, DAPC) and model-based Bayesian approaches (Structure, ABC) allowed getting an insight into phylogeny and demographic processes in eastern-Mediterranean firs.
All species of the section Abies are easily crossable with each other, as shown by artificial crossing experiments, which produced viable and mostly well-growing and fertile hybrids (Greguss & Paule, 1988 Ata and Merev (1987) and Klaehn and Winieski (1962). In both cases, the putative hybrid taxon was shown to represent a quite recent evolutionary branch of A. bornmuelleriana.
The ABC-based estimates of divergence times and effective population sizes should be considered a crude approximation. The problem of a small number of loci was already mentioned, and tandem-repeat loci are known to be more prone to homoplasy than the other types of markers. Moreover, we have no reliable information about mutation rates in nSSR loci in the genus Abies and relied on the study performed in another conifer genus (however, the mean mutation rate used by only two geographically very proximate silver fir populations from the Balkans were included in the analysis, which do not represent the gene pool of the species as such. This is also reflected in a very broad confidence interval of the divergence time in this case. Moreover, as suggested by Linares (2011), gene flow between fir populations on both sides of the present Bospor strait was probably maintained during most of the Pleistocene.
For the interpretation of ABC outcomes in terms of a temporal scale, a correct estimation of the average generation turnover time is crucial. In long-lived forest trees producing gametes for a large part of their lifespan, this is always associated with a big uncertainty, and true generation times may depend from species and ecological conditions: for example, for A. cilicica growing in forests with looser stand structure (Mayer & Aksoy, 1986), the average of 150 years may be an overestimation. If the proposed generation time is accepted, most speciation events occurred during the late Pliocene and Pleistocene. Thus, the estimated times of species or subspecies divergence fit well with the outline of the phylogeny of Abies in the eastern-Mediterranean space as drawn by Linares (2011) and Palamarev (1989).

| Implications for taxonomy and conservation
As already mentioned, the classification of A. cilicica has changed in the past, from section Piceaster to section Abies. However, these sections are both geographically and genetically close to each other anyway (Semerikova & Semerikov, 2014;Xiang et al., 2009Xiang et al., , 2015. Still, the detected high differentiation of A. cilicica from the rest of the studied species is in line with the fact that A. cilicica is the least consistently classified into the section Abies and confirmed earlier findings of Scaltsoyiannes et al. (1999) and Bergmann et al. (2013).
The taxonomical status of the three subspecies of A. nordmanniana s.l. is the most complicated issue in this study. Only A. equi-trojani is recognized by Euro+Med PlantBase (Euro+Med, 2006) and IUCN (Knees & Gardner, 2011) and the genus-level classifications usually deal only with A. nordmanniana as a whole species (Farjon & Rushforth, 1989;Liu, 1971;Xiang et al., 2009 (Bergmann et al., 2013;Liepelt et al., 2010;Linares, 2011;Scaltsoyiannes et al., 1999). Contrary to the former classification, our findings clearly differentiated A. bornmuelleriana from the nominate subspecies. A. equi-trojani was actually the group that was less clearly differentiated and tended to cluster under A. bornmuelleriana. Of course, A. bornmuelleriana is given as a synonym to A. equi-trojani in the former classification, but the ranks given by the IUCN Red List do not agree with our results. Our findings support the classification given by Coode and Cullen (1965): A. bornmuelleriana and A. equi-trojani are distinct groups, their differentiation is too low to support the classification as separate species. Nevertheless, from the point of view of conservation and forestry practice, it must be mentioned that all three taxa are treated as separate species in the OECD scheme for forest reproductive materials certification or in the European program for conservation of forest genetic resources Euforgen (Alizoti et al., 2011). There is not enough support for distinguishing A. ×olcayana as a separate taxon (whether hybridogenous or not) as suggested by Ata and Merev (1987).
The studied species are generally classified as Least Concern in the IUCN Red List (http://www.iucnredlist.org), except A. cilicica s.l., which is Near Threatened as a whole, although local populations (mainly outside Turkey) are considered Critically Endangered. The outcomes of our study generally support this view. Even though genetic variation at neutral loci does not directly imply high adaptive variation, the observed relatively high levels of allelic richness and gene diversity within the studied taxa, including endemic and quite isolated A. equi-trojani, are promising for their further survival. The only exception is A. isaurica; limited allelic richness and a strong differentiation from all remaining taxa make it a primary potential target of gene conservation efforts.

ACKNOWLEDGMENTS
Thanks are due to C. Ata, D. Güney, and Z. Yahyaoğlu, and the Leiba for providing the Abies nordmanniana samples from Russia and Georgia. Technical assistance of G. Baloghová is greatly appreciated.
We also thank to K. Willingham for linguistic correction. The study was supported by a grant of the Slovak Grant Agency for Science VEGA 1/0269/16.

CONFLICT OF INTERESTS
None declared.

AUTHOR CONTRIBUTIONS
LP designed and organized the study, MH, DG, and DK made the analyses and statistical treatment, SK, HS, ITu, and ITv designed sampling and provided materials, MH and DG wrote the first draft, all authors contributed to the final text.