Fine-Scale Spatial Patterns of the Genetic Diversity ofEuropean Beech (Fagus sylvatica L.) around a Mountainous Glacial Refugium in the SW Balkans

Abstract

Beech (Fagus sylvatica L.) is one of the most important forest trees in Europe and its current broad expanse over the continent is believed to be the outcome of the Holocene postglacial expansion of lineages originating from different glacial refugia. Studies using gene markers, pollen profiles and fossils suggest the main locations of glacial refugia in Southern Europe. In this study, we conduct a fine-scale genetic study on the beech populations surrounding the Almopia basin, an area that is said to have hosted one of the main European glacial refugia for several plant and animal species during the Pleistocene Epoch. We test the hypothesis of the existence of a local refugium in the study area for beech to understand the spatial genetic pattern of the putative refugial beech populations in the area and to investigate possible genetic connections between the local beech populations and the European expansion of the species. The genetic diversity of 100 sampled trees in 20 plots representing the expansion of beech in the area was studied using chloroplast and nuclear DNA microsatellites (cpSSR and nSSR, respectively). All three cpSSR regions were polymorphic, resulting in eight haplotypes, separated spatially in two distinct groups (one on the western and the other on the eastern part of the Almopia basin) that correspond to two different postglacial beech lineages. Furthermore, the cpSSR sequences of the eastern lineage are genetically identical to those of beech populations extending over central and northern Europe. The nSSR markers were highly polymorphic, and the trees studied were separated into two genetic groups that coincided with the cpSSR ones in locations where the topography is more pronounced. These results indicated that the Almopia region was indeed a major refugium for beech that possibly produced two main postglacial lineages for Europe, one of which is connected with the majority of beech populations growing on the continent. These southern refugial populations are important diversity centers that need to be the subject of special management and conservation.

1. Introduction

European beech (Fagus sylvatica L.) is one of the ecologically and commercially most important trees in Europe and forms extensive forests from the southern parts of the British Isles and Scandinavia in the north to Sicily and the southern Balkans in the south and from the Spanish Pyrenees in the west to Poland and Ukraine in the east [1,2]. In the southern part, beech forms island-type populations at higher altitudes on mountains [3,4,5,6]. The beech forests that currently cover the Mediterranean, such as on the Italian and the Balkan peninsulas, are considered to be adapted to warmer and drier conditions and harbour significant levels of genetic diversity, due to a complex history of migration and fragmentation that followed climatic changes during the Quaternary Period [7,8,9,10,11].

The current European expansion of beech is probably the greatest in history and is considered to be the outcome of the postglacial expansion of lineages from a number of southern glacial refugia after the Last Glacial Maximum [12,13,14,15,16,17,18,19,20]. These refugia were areas of relative ecological stability, with moisture and varied topography that ensured protection from extreme climatic conditions [11,21,22]. Southern glacial refugia sheltering native tree populations are considered important biodiversity hotspots and the identification of their locations has become a conservation priority [20,23,24]. Possible locations of glacial refugia especially for beech forests have been described in a number of paleobotanical and genetic studies on the European scale [12,13,14,15,16,17,18,19,25,26,27,28,29,30,31,32]. These studies and several more for specific countries [11,33,34,35] have suggested the existence of larger or smaller glacial refugia for beech in Greece, southern and northeastern Italy [36], Slovenia, northern Spain and southwestern France. The presence of a northern microrefugium in the Carpathian Mountains is debated [20]. While most authors agree about the existence of the southern European refugia, there are different opinions on their postglacial contribution to the current distribution of beech in Europe. Based on the fossil pollen maps of Huntley and Birks from 1983 [12], several authors support the hypothesis of a main beech lineage originating from the southern part of the Balkan Peninsula and spreading throughout Europe. According to this scenario, lineages originating from southern Italy were blocked by the Alps and contributed less to the European distribution of beech [15,26,28,29,37,38,39]. Magri et al. in 2006 [32] combined genetic data (nuclear and chloroplast DNA markers) and fossil data to suggest a different scenario: lineages deriving from refugia in the Balkans did not expand to the north during the Holocene Epoch. The recolonization of beech to Central and Western Europe occurred with lineages that originated from glacial refugia at the northern part of the Mediterranean region close to the Slovenian Alps and southern France. This hypothesis considers a complex pattern of postglacial migration in which beech lineages moved at different temporal rates and directions, influenced by local conditions [11,19,32,40,41].

In Europe, genetic studies of beech populations have been used to infer the location of glacial refugia and to help reconstruct the postglacial recolonization routes of beech lineages [24,26,29,30,32]. These studies describe a similar diversity pattern, leading to an increase in the richness of genetic variants (alleles at the nuclear loci or cpDNA haplotypes) in southern and eastern Europe, in contrast to the more uniform situation in central and northern Europe. Haplotypes of the chloroplast genome have been proven especially useful for describing the postglacial movement of beech populations because of their uniparental inheritance. In a large-scale assessment of beech diversity in Europe from 2006, Magri et al. [32] described a specific pattern of cpDNA SSR-based haplotypes, with a single one being predominant in central and northern Europe, but a different one for Italy and a large diversity of haplotypes on the southern Balkan Peninsula. In Greece, this complex haplotype diversity pattern was confirmed by an extensive study of beech populations [34]. This study described at least four postglacial lineages for beech on the southern part of the Balkan Peninsula, possibly originating from local or distant glacial refugia. The typical European haplotype was abundant in the northern and western part of Greece, with the highest frequency observed in the mountains surrounding the valley of Almopia, close to the border with Northern Macedonia. This region is said to host refugial populations for other central European forest tree species, such as Abies alba Mill. [42,43,44], Quercus sp. [45,46], Castanea sativa Mill. [47,48,49] and Alnus glutinosa L. [50].

Most of the genetic studies mentioned above presented scenarios of postglacial history for beech over extensive areas (countries or continents) and therefore used geographically distant populations and a coarse sampling resolution. Some more recent studies have focused on the description of spatial, genetic and diversity patterns of populations connected with specific refugial areas on a smaller geographic scale. Such small-scale genetic studies of putative refugial areas in the mountains of Rodopi and Paggeo in Greece [22,51], southwestern France [11] and southern Italy [52] revealed high genetic diversity and complex spatial patterns, indicating that fine-scale studies are needed for a better understanding of the postglacial history of temperate plant species. The aim of this study is to investigate the fine-scale spatial patterns of genetic diversity of beech on the mountains surrounding the Almopia basin, using cpDNA haplotypes and nuclear microsatellite loci. We will further test the scenario that this mountainous area is a refugial area for beech and attempt to identify the postglacial lineages present in the area. Using these results, we aim to contribute to the discussion of the postglacial history of beech in the Balkans and the rest of Europe.

2. Materials and Methods

2.1. Area of the Study and Sampling

The region of Almopia lies in north-central Greece and comprises mainly an extensive plateau surrounded by an arc-shaped mountain range that effectively shields the area from the north, east and west. The surrounding mountains and their respective highest peaks include Voras (2524 m) to the west, Tzena (2182 m) and Pinovo (2154 m) to the north, and Paiko (1650 m) to the east. This mountain range is part of the easternmost branch of the Dinaric–Pindic high-mountain system that runs from the northwest Balkans to the southeast and forms an uninterrupted geomorphological connection with central Europe. The climate of the study area is humid-continental with a short summer [53]. The topography of the area is mostly rugged and diverse. Numerous mountain ridges descend to more or less steep slopes which are separated by broad or narrow valleys.

Under prevailing site conditions, beech forests occupy the mountain slopes between 600 and 1700 m a.s.l. (Figure 1). They grow especially well in deep fresh soils and find optimal development above 900 m. Based on the geographic distribution of the beech forests in the study area, 20 representative plots were chosen (Figure 1). Leaf samples were collected from five trees in each plot, with a minimum distance of 100 m between them. The sampled material was stored below −20 °C. DNA from four more beech samples from two locations in Switzerland (northeast in the Rhine valley and southwest in the Rhone valley) representing the main European cpSSR haplotypes [54,55] were used to compare with the Greek samples.

2.2. Laboratory Procedures

DNA was extracted from fresh leaves using the DNeasy^® Plant Kit (Qiagen, Hilden, Germany) following the instructions of the manufacturer. DNA quantity and quality were tested on 1% agarose gels after staining with ethidium bromide. For the chloroplast markers, amplification was conducted using three universal microsatellite primer pairs (cpSSRs), labeled with fluorescent dyes: ccmp4, ccmp7 and ccmp10 [56]. These markers were chosen because they have shown diversity and were broadly employed to describe cpDNA haplotypes and to study the postglacial history of beech both in Greece and in the rest of Europe [22,32,33,34,57]. The fragments derived from amplification were sequenced and published for Fagus sylvatica by Vettori et al. in 2004 [33] (GenBank accession numbers AY563559; AY563561; AY563558) and Papageorgiou et al. in 2017 [58] (GenBank accession numbers KY496670 to KY496688). The homology to these sequences was verified in DNA fragments from plant material from this study. Amplification followed Gailing and von Wuehlisch from 2004 [57], and the products were purified using the QIAquick Gel Extraction kit (Qiagen^®, Germany) following manufacturer’s specifications. Sanger sequencing reactions [59] were run on the ABI 3100 Genetic Analyser (Applied Biosystems) using the Big Dye^® Terminator v.3.1. Cycle Sequencing Kit (Applied Biosystems).

Six nuclear Expressed Sequence Tag (EST) microsatellite loci were used to genotype all individuals sampled (Table 1).

Three loci were originally developed for Fagus crenata Blume (FcC00464, FcC00927, FcC1877; [61]) and three for Quercus robur L. (GOT066, FIR065, FIR004; [63]). Primers were labeled with different fluorescent dyes and amplification procedures followed [64]. Fragments for nSSR markers were separated on the ABI 3100 Genetic Analyser (Applied Biosystems) using the internal size standard Gene Scan Standard 500 ROX (Applied Biosystems). Fragment sizes were determined using Genescan 3.7 and Genotyper 3.7 from Applied Biosystems.

2.3. Data Analysis

A cpSSR haplotype was assigned to each sampled tree based on a combination of the different length and sequences of the fragments of the three cpSSR loci. A minimum spanning haplotype network based on the distances between haplotypes in mutation steps was created using HapStar [65]. The frequencies of haplotypes were estimated for each plot. The number of haplotypes (Nh) and their effective number (Ne(h)) per plot were calculated using the GenAlEx 6.5 software [66]. An analysis of molecular variance (AMOVA) was carried out to separate total haplotypic diversity into two parts, among plots and within plots, using GenAlEx 6.5 software.

For each nSSR locus, fragments with different length were considered “alleles” and their frequencies within each plot were estimated with GenAlEx 6.5. The same software was used to calculate the number of alleles (Na), effective number of alleles (Ne), observed heterozygosity (Ho), expected heterozygosity (He) and inbreeding coefficient (F). Using the same software, two AMOVAs were performed to describe the distribution of the total diversity at three levels: among plots, within plots and within individuals. One AMOVA was based on genetic differentiation calculated solely by allele frequency (F_ST) and the other by considering “molecular” genetic distances between microsatellites as well (R_ST).

For both cpSSR and nSSR loci, each tree was assigned to a Bayesian panmictic group (cluster) using a MCMC simulation with the R package Geneland [67] and the non-spatial uncorrelated algorithm. After 10 simulations with 100,000 iterations each, the results of the run with the best mean posterior density were kept. All figures including maps were prepared using QGIS software [68].

3. Results

3.1. Diversity and Differentiation of cpSSR Haplotypes

All cpSSRs produced DNA fragments that were polymorphic in length and when combined resulted in eight beech cpDNA haplotypes (Table 2). The sequences of the cpDNA fragments demonstrated that the observed variation in length was due to different numbers of specific mononucleotide SSR repeats. No obvious fragment length homoplasy was detected. Three haplotypes were frequent (h01, h05 and h06) while the rest had low frequency. Four of the rare haplotypes were recorded for the first time (h03, h04, h07 and h08) in beech.

The minimum spanning network (Figure 2) showed that the two most frequent haplotypes (h01 and h06) were separated by two single mutation steps in different cpSSR fragments, while h05 and h02 were located in an intermediate position separated from the two frequent haplotypes by one mutation step. The rare haplotype h03 was possibly derived from a mutation of h01 and the other rare haplotypes (h08, h07 and h04) seemed to be related to the more frequent h06 and h05. This scenario was supported by the geographic distribution pattern of the cpSSR haplotypes (Figure 3). The most frequent h01 and the closely related h02, h03 and h04 were found exclusively on Mt. Voras on the western mountain range of the Almopia basin. On the eastern side of the study area (Mt. Tzena, Mt. Pinovo and Mt. Paiko), only haplotypes associated with the frequent h06 (h05, h07 and h08) were observed. An exception to this pattern was the dominance of h05 and h06 in the plots located in the western part of the basin, centrally on Mt. Voras, close to the Loutraki gorge. All five individuals from the Swiss population carried h06.

Rare haplotypes were observed only in the southernmost plots of Mt. Voras (h02, h03, h04) and in plot 20 on Mt. Paiko (h07, h08). Genetic diversity differed largely within plots with the average number of haplotypes per plot being 2. While some plots were monomorphic, plots 2, 5, 13 and 20 carried more than two haplotypes and therefore demonstrated the highest diversity at cpSSR markers (Table 3). Differentiation among the plots was very high (Φpt = 0.794) since the largest part of the total genetic diversity of cpSSR haplotypes (79%) was among plots, as shown by the AMOVA, with diversity within plots being just 21% (Table 4). The Bayesian inference of the cpSSR haplotypes classified the trees into two genetic clusters corresponding to the groups seen in the haplotype network and their geographical coverage. Cluster 1 included trees that carried haplotypes belonging to the group of h01 located on the western side of the Almopia basin (Mt. Voras). Trees belonging to Cluster 2 carried h06, or one of the other haplotypes that belonged to this group, and were mainly found on the eastern part of the study area on Mt. Tzena, Mt. Pinovo, Mt. Paiko and Loutraki on the western side (Figure 4).

3.2. Diversity and Differentiation at nSSR Loci

A large number of different alleles (85) were scored in total for the six nuclear EST-SSR loci used, ranging from 9 alleles at GOT066 to 25 alleles at FcC1877 (Table 5). GOT066 was the least diverse marker across all plots (Na = 1.9, Ne = 1.329, Ho = 0.243, He = 0.211), while the highest diversity was observed at FcC1877 (Na = 6.35, Ne = 5.424, He = 0.806) and the highest observed heterozygosity at Fir065 (Ho = 0.818). The average inbreeding coefficient F was slightly positive for four markers, but reached negative values for Fir065 (F = −0.167) and GOT066 (F = −0.123).

Genetic diversity within the 20 plots of the study was very high, despite the small sample size for each plot (Table 3). The mean number of alleles (Na) over all loci ranged between 4 and 5 in all plots. The genetic diversity indices were highest in plot 17 (Ne = 4.968, He = 0.67) and lowest in plot 1 (Ne = 2.846, He = 0.557). The observed heterozygosity (Ho) ranged between 0.5 in plot 18 and 0.833 in plot 8. The mean inbreeding coefficient F value was −0.03 across all plots and most plots demonstrated similar values. However, the highest heterozygote deficit was observed in plot 18 (F = 0.226) and the largest heterozygote excess in plot 8 (F = −0.236). The percentage of total diversity explained as differentiation among plots in the AMOVA (Table 4) was very low, both based solely on differences in allele frequencies (F_ST = 0.023) and when distances between microsatellites were considered (R_ST = 0.006). In the first case, when only differences in allele frequencies were taken into account, the largest part of the total genetic diversity was found within individuals (87%), thus being the result of high heterozygosity. However, when the distances between microsatellites were considered, the largest proportion of diversity was found among individuals within plots (77%), an indication that in most heterozygous genotypes the two alleles were not very different in the number of microsatellite repeats.

The analysis of the Bayesian inference of the nSSR alleles resulted in two genetic clusters. Cluster 1 was most common among the trees of the study and is widely spread in almost all plots of the study. Cluster 2 was present in most plots as well, but with lower frequency than Cluster 1, with the exceptions of plot 16 on Mt. Paiko, where Cluster 2 was more abundant, and plots 7 and 8 on the western side near the Loutraki gorge on Mt. Voras, which were monomorphic for Cluster 2 (Figure 5).

4. Discussion

4.1. cpSSR Haplotypes and Postglacial Lineages

The frequent cpSSR haplotypes found in the beech populations of Almopia were connected with two postglacial lineages for beech in previous genetic studies that used the same markers. The most common haplotype, h01, which dominates the western part of the Almopia basin, had previously been found in central and northwestern Greece (North Pindos, Mt. Smolikas, Mt. Grammos and Mt. Voras) and is considered to be characteristic for the “Pindos” postglacial lineage [34]. This lineage has not been described outside of Greece, but a possible presence in other parts of the western Balkans and possibly Italy cannot be excluded. Common beech cpDNA haplotypes between Greece and Italy were suggested by Vettori et al. in 2004 [33], but gene flow between plant populations on the southern European peninsulas is considered unlikely [23].

The other common beech haplotype in our study was h06, found mainly on the eastern side of the mountains of the Almopia basin and in the Loutraki ravine on the western side. The three cpSSR fragments that comprise h06 in this study had identical sequences with the ones obtained from the trees from Switzerland. Furthermore, the same haplotype was described as being dominant in the main part of the European beech expansion and in parts of central and northern Greece [32,34,57]. This haplotype is characteristic for the postglacial “Europe” lineage that seems to have been the main contributor of the current expansion of beech on the continent, as explained in the introduction. Two mutation steps separated h01 from h06, while h05 differed in one mutation step from each one of the two most frequent haplotypes. This intermediate haplotype was less abundant than the previous two but was still frequent and present on both sides of the Almopia basin. In Greece, h05 was rare and found mainly in connection with h06; thus it was also considered to belong to the “Europe” lineage [34].

4.2. Spatial Genetic Pattern at cpSSR and nSSR Markers

The overall genetic diversity of beech in the Almopia basin was high in this study, both for cpSSR and nSSR markers. While a reliable comparison of genetic diversity levels within plots was not possible due to the small number of individual trees in each plot, clear trends of plot differentiation were observed. Previous genetic studies of cpSSR diversity at a fine scale in small geographical areas found similar levels of diversity and differentiation in beech [22,51]. At cpSSR markers, two different groups of haplotypes occupied two different sides of the Almopia basin, as explained previously. Plots were very different from each other, even in cases such as the Loutraki gorge where plots 7 and 8 were obviously different from the other plots on Mt. Voras, despite them being geographically close.

Overall differentiation at the nSSR loci was very low, which is an expected trend for the usually hyper-variable nuclear microsatellite markers [55,69]. However, cases of plots that differed greatly in their Bayesian cluster frequencies while being in close proximity to each other were observed. Interestingly, such a case was observed on the northern slope of the Loutraki gorge, where plots 7 and 8 are genetically different from most of the other nearby plots. This pattern is very similar to the one observed at the cpSSR markers. These findings are further supported by the AMOVA results of the nSSR loci, which reveal a strong phylogeographic signal for the differences between plots in the study area.

4.3. Several Smaller Refugia within a Larger Refugium

Both beech postglacial lineages of the Almopia basin described in this study appear to be part of the current expansion of the species F. sylvatica in Greece and European regions further north. No haplotype related to other beech species found on the southeastern part of the Balkans, F. orientalis Lipsky, was found despite the presence of characteristic haplotypes reported in other mountains of central and northern Greece [34]. Beech populations originating from both Almopia lineages “Pindos” and “Europe” coexist in a small geographical area in close proximity to each other. This coexistence of different postglacial lineages of beech in areas with high topographical diversity and intense relief has also been observed in other Greek mountains, such as Mt Paggaio [22], Mt Menoikio [70] and the mountains of eastern Rodopi [51]. The fine-scale spatial genetic pattern of the lineage specific cpSSR haplotypes indicates that both beech postglacial lineages found in the Almopia basin probably derived from the Holocene expansion from glacial refugia that would have existed in the region not far from the current beech populations. This scenario is initially supported by the complex landscape of the mountains surrounding the Almopia basin. Beech populations in Almopia spread on slopes over protected valleys and gorges, indicating a possible upward movement across a rather short distance from smaller glacial refugia to higher altitudes after the end of the last glacial period. It is broadly believed that glacial refugia of temperate plant species were climatically stable, generally located in moist, protected valleys or ravines that maintained plant populations during the recurrent ice advance [9,11,22,23]. Thus, the spatial distribution pattern of the cpSSR haplotypes and the complex topography of the region indicate the existence of several smaller, possibly isolated, glacial refugia along the mountains surrounding the Almopia basin. The studied beech populations appear to have originated from these refugia and their postglacial movement must have occurred upwards on the slopes of the mountains, over short distances.

The frequencies of the cpSSR haplotypes and the molecular distances among them further support the possible origin of these two postglacial beech lineages from glacial refugia that existed in the Almopia basin. The two most common haplotypes, h01 and h06, which are characteristic for the ‘’Pindos’’ and ‘’Europe’’ lineages, respectively, have a small molecular distance from each other, while two intermediate haplotypes were detected, one of which was rather frequent (h05). This indicated that these two lineages may have had a common origin from a glacial refugium in the Almopia basin. The latitudinal and altitudinal contraction and expansion of populations belonging to these lineages during subsequent glacial and interglacial cycles may have resulted in the formation of smaller isolated refugia within the broader initial glacial refugium. This pattern of smaller glacial refugia coexisting in isolation from each other in close proximity within larger refugial areas, separated by topographic genetic barriers belonging to a complex mountainous landscape, seemed to be common for the southern parts of the temperate European tree species expansion, the so called “rear edge” [9,22,24,71,72,73,74,75,76].

Another argument that supports the existence of glacial refugia in the Almopia basin is the accumulation of rare cpSSR haplotypes in a rather restricted number of plots. More important, these rare haplotypes that coexist in the same location have a short genetic distance to each other and to any frequent haplotype close by, indicating that they are separated from each other by a few mutation steps and are phylogenetically related. These mutations probably occurred during long glacial phases in isolation in these specific locations. Based on the topography of the area and the occurrence of the rare haplotypes, we can identify the possible locations of glacial refugia for the two beech lineages in Almopia. These possible locations are in ravines on the southern tip of Mt. Voras for the “Pindos” lineage. Accordingly, the valley north of Loutraki and west of Loutra on Mt. Voras and the ravines on the north side of Almopia east of Mt. Tzena, as well as the southern part of Mt. Paiko can be considered possible locations for the refugial origin of the “Europe” lineage. On the north side of Mt. Voras and in plots 10–13, frequent haplotypes belonging to both lineages occur together, indicating a possible meeting point at higher altitudes.

Furthermore, since genetic drift is expected to have eliminated rare genetic variants in populations that migrated long distances from their refugia of origin, we can assume that current beech populations on the mountains around the Almopia basin are not located far from their original glacial refugia. In any other case, any haplotype would be abundant or non-existent. Thus, the existence and the spatial distribution of rare cpSSR haplotypes suggest that beech refugial populations performed altitudinal shifts to cope with changing climate during the interglacial periods and therefore managed to survive in the same broader refugial area during the Pleistocene glacial cycles and maintained high levels of genetic diversity. The scenario of refugia maintaining tree populations over a longer period of time, including several glacial cycles, has been suggested by several authors [9] and was confirmed by the pollen record in regions of northern Greece, such as Tenaghi Philippon [77] and Ioannina [23].

The results of the nSSR markers confirm the overall high genetic diversity for the beech populations in Almopia observed at the cpSSR markers. Despite the suggested long-term isolation of several smaller refugial populations and the small sample size of each plot (5 trees), the diversity within plots was very high in almost all cases. This may be the result of current gene flow through pollen, which connected populations originating from different small scale glacial refugia over short distances. It is worth noting that high differentiation between nearby populations was detected on slopes north of the Loutraki gorge at both cpSSR and nSSR loci. It seems that the pronounced mountainous topography of this area acted as a gene flow barrier for pollen transfer. A positive correlation between genetic variation at cpDNA haplotypes and nuclear SSRs in beech was also observed in Italian refugial populations [52].

4.4. A Possible Glacial Refugium for the Main European Lineage

Our study suggests that the postglacial “Europe” lineage found in the Almopia basin is the same, or closely related to, the main beech European lineage that dominates the central and northern parts of the continent [32]. This is supported by the identical sequences of the three cpSSR fragments of h06 with the ones observed in the Swiss sample and the existing sequences in the nucleotide databases. Although genetic differences between the Almopia and the European beech trees may exist at other sites of the cpDNA genome besides the three cpSSR fragments considered here, the fact that the main Europe-wide study of beech populations [32] used the same three fragments supports our suggestion. Following the results of this study, we assume that the glacial refugium of Almopia may have contributed to the current expansion of beech in Europe, possibly along with other refugia in the southwestern Balkans and further north.

Previous genetic studies considered the principal European beech refugium to be located in the southern part of the Dinarian range [29], a broader region immediately north of the Almopia basin. Several authors [37,38] support the hypothesis that is based on the fossil pollen maps of Huntley and Birks (1983) [12], suggesting that beech followed a north–west expansion from the SW Balkans along the Dinarian to the southern slopes of the Alps and from there to the rest of Europe. It seems that continuous mountain ranges do not act as geographical barriers to postglacial expansion of tree lineages; on the contrary, they assist this expansion because they allow altitudinal shifts while the climate changes [32]. This may be the case for the long mountainous continuum starting at the mountains in northern Greece and stretching over the Dinarians to the southeastern slopes of the Alps, where the location of another major glacial beech refugium is thought to be [19,32]. It is probable that different nearby or even distant glacial refugia were not differentiated genetically and may have hosted the same postglacial lineage during the LGM, due to migration events in older interglacial periods [29].

Our findings initially seem to contradict the scenario proposed by Magri et al. in 2006 [32] and later supported by Hatziskakis et al. in 2009 [34], that the current European expansion of beech derives mainly from glacial refugia on the eastern side of the Alps and that the populations in the southern Balkans derived from an expansion from these refugia to the south. Since the existence of a glacial refugium for beech in the Almopia basin—and possibly in other similar regions in the southern Balkans—is considered possible, it is highly improbable that beech arrived in the region from the north. We considered that the existence of several non-differentiated glacial refugia, on both sides of the Dinarians was far more probable and that the movement of the main European lineage probably took place during an earlier interglacial period from Almopia to the north, not vice versa. The secondary glacial refugium of the Alps possibly originated in the SW Balkans and may have indeed made the most significant contribution to the current geographic distribution of beech on Europe, as suggested by Magri et al. in 2006 [32]. In this case, and based on our results, the Almopia basin is a strong candidate for the source of the major European beech lineage. It is worth mentioning that a rare population of the tertiary relict pine Pinus peuce Griseb. survives in this area, as do several boreal and central European floristic elements, such as the southernmost population of silver fir (Abies alba Mill.) [43] and many companion beech and fir forest species of northern origin [78]. Finally, the Loutraki gorge on Mt. Voras has been suggested as a major glacial refugium for several species of fauna [79,80].

The Almopia basin is one of the numerous southern European locations that made an important contribution to the evolutionary history of temperate tree species and is a significant diversity center for them because it harbors the “rear-edge” populations that survived through the climatic fluctuations of the Quaternary Period, or for even a longer time [8]. However, these southern European populations will probably suffer most from the climate change that is currently on the way and are at risk of becoming extinct [8,81,82]. Further studies should focus on glacial populations of temperate tree species in Europe and elsewhere and their findings should be connected with management and conservation efforts.

Author Contributions

Conceptualization, A.C.P. and R.F.; samples collection, O.T.; methodology, A.C.P.; software, A.C.P., O.T.; validation, A.C.P., O.T.; formal analysis, A.C.P., O.T.; data curation, A.C.P., O.T.; writing—original draft preparation, A.C.P., O.T., G.K., L.L., O.G. and R.F.; writing—review and editing, A.C.P., O.T., G.K., L.L., O.G. and R.F.; supervision, A.C.P. and R.F.; project administration, R.F. All authors have read and agreed to the published version of the manuscript.

Funding

The research received no external funding.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Acknowledgments

We would like to thank Alexandra Dolynska and Apostolos Manolis for their assistance in the laboratory. We acknowledge the support of the Open Access Publication Funds of the Georg-August University of Göttingen. We would also like to thank the reviewers for comments on the paper.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

cpSSRs	chloroplast Single-Sequence Repeats
nSSRs	nuclear Single-Sequence Repeats
EST	Expressed Sequence Tag
AMOVA	Analysis of Molecular Variance
MCMC	Markov chain Monte Carlo
QGIS	Quantum geographic information system
LGM	Last Glacial Maximum
SW	Southwest
SE	Standard Error

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Map of the Almopia basin with the 20 sampling plots of this study. Light-gray shading shows mountains and dark-gray shading the expansion of beech.

Figure 2. Minimum spanning network of the cpSSR haplotypes. Each line denotes a mutation step. The numbers refer to the absolute frequency of each haplotype in the study area.

Figure 3. Map of the geographical distribution of the cpSSR beech haplotypes in the 20 plots of this study.

Figure 4. Map of the average Bayesian cluster membership probabilities of the cpSSR beech haplotypes in the 20 plots of the study.

Figure 5. Map of the average Bayesian cluster membership probabilities of the nSSR gene loci for beech in the 20 plots of the study.

Table 1

Characteristics of 6 simple sequence repeat (SSR) markers [60,61,62].

Locus(Cluster ID)	Motif Type/Repeat Number	Forward Primer Sequence (5′–3′)	Reverse Primer Sequence (5′–3′)
GOT066	(GAA)10	TCCCTAGATGATGGGGATGA	TTTTACGTCGGCCAACTTTT
FIR004	(CT)18	TCTCTCTCAGGGCAGCTTCT	AACCAAACTCAGATCCAGATTCA
FIR065	(CTT)9	ATTCCCATGCATCAAAATCC	TCCTTCAGTTTGAGAGCTCCTT
FcC00464	(AG)3GGCTCCATTACCTGGTCCT(AG)16ATTGAAGGGCACCATTACTAGTCCT(AG)14	GCACGAGGGTGAGAGTGAGAGTA	GCCCGTGTGTGAACATGTAGGAC
FcC00927	(TC)3CCAAAAAAAAAA(AC)3(TC)16	CCTGCCGCTCAATTATTAGTCCC	ACGAGTCAGCTGGAAAACACGAG
FcC01877	(TC)15	TCATAAACCCAAAATGCCTTCCC	CATTTGCTCCAAGAAGTCGTCGT

Table 2

Description and absolute frequencies of cpSSR haplotypes for beech in the Almopia basin.

	DNA Fragment Length (bp)
cpSSR Haplotype	ccmp4	ccmp7	ccmp10	Individuals
h01	116	144	115	42
h02	115	144	115	3
h03	117	144	115	2
h04	116	144	117	1
h05	116	145	115	14
h06	115	145	115	35
h07	116	145	116	1
h08	115	145	116	2
Total				100

Table 3

Measures of genetic diversity at cpSSR and nSSR markers for the 20 plots of the study.

		cpSSR	nSSR
Plot	N	Nh	Ne(h)	Na	Ne	Ho	He	F
1	5	1	1.000	4.000	2.846	0.567	0.557	−0.007
2	5	4	3.571	4.333	3.512	0.633	0.590	−0.051
3	5	2	1.471	4.000	3.291	0.567	0.523	−0.102
4	5	2	1.471	4.500	3.579	0.667	0.637	−0.076
5	5	3	2.273	4.833	3.624	0.700	0.663	−0.085
6	5	2	1.471	4.833	3.802	0.667	0.607	−0.078
7	5	2	1.471	3.667	2.958	0.633	0.583	−0.096
8	5	2	1.923	5.000	3.889	0.833	0.687	−0.236
9	5	2	1.471	5.000	4.121	0.633	0.623	−0.001
10	5	2	1.471	4.833	3.822	0.533	0.617	0.072
11	5	1	1.000	4.833	3.757	0.533	0.593	0.045
12	5	1	1.000	4.667	3.172	0.533	0.597	0.054
13	5	3	2.778	5.000	3.830	0.733	0.673	−0.098
14	5	2	1.471	4.667	3.840	0.633	0.680	0.037
15	5	1	1.000	4.333	2.959	0.600	0.567	−0.106
16	5	2	1.923	4.833	3.869	0.667	0.663	0.076
17	5	1	1.000	4.667	4.968	0.733	0.670	−0.097
18	5	1	1.000	5.000	4.036	0.500	0.643	0.226
19	5	2	1.471	4.333	3.172	0.667	0.640	−0.081
20	5	4	3.571	4.500	3.904	0.583	0.633	0.018
Mean		2	1.657	4.642	3.648	0.631	0.622	−0.030
SE		0.201	0.172	0.174	0.166	0.025	0.022	0.023

N—number of individuals, Nh—number of haplotypes, Ne(h)—effective number of haplotypes, Na—number of alleles, Ne—effective number of alleles, Ho—observed heterozygosity, He—expected heterozygosity, F—inbreeding coefficient, SE—standard error.

Table 4

AMOVA classification of genetic diversity based on cpSSR haplotypes, F-statistics of nSSR markers using allele distances and R-statistics of nSSR markers using microsatellite distances.

Percentage of Variation	cpSSR Haplotypes	nSSR FST	nSSR RST
Among plots (%)	79	2	1
Among trees (%)	21	11	77
Within trees (%)	-	87	22
	Φpt = 0.794 (p = 0.010)	FST = 0.023 (p = 0.001)	RST = 0.006 (p = 0.001)

Table 5

Average genetic diversity indices over all plots for each nSSR marker.

	N	Na	Ne	Ho	He	F	Total N
Fir065	4.950	4.700	3.760	0.818	0.709	−0.167	9
GOT066	4.950	1.900	1.329	0.243	0.211	−0.123	3
Fir004	4.950	3.450	2.210	0.505	0.490	−0.056	14
FcC00464	5.000	5.900	4.878	0.750	0.781	0.045	17
FcC00927	5.000	5.550	4.285	0.680	0.737	0.071	17
FcC1877	5.000	6.350	5.424	0.790	0.806	0.022	25
Mean	4.975	4.642	3.648	0.631	0.622	−0.030	14.167
SE	0.014	0.174	0.166	0.025	0.022	0.023	3.081

N—number of individuals, Na—number of alleles, Ne—effective number of alleles, Ho—observed heterozygosity, He—expected heterozygosity, F—inbreeding coefficient, F_ST —differentiation, SE—standard error.