Distribution of genetic diversity in Fagus sylvatica at the north-eastern edge of the natural range

An understanding of the genetic variation of the beech, especially at the edge of its natural distribution, is essential because of the change in natural distribution of the species resulting from changing climatic conditions. The main aim of the study was to determine the level of genetic diversity of European beech at the north-eastern edge of its natural range. The other aim was to check the genetic variation of beech from the two centres, the north and the south of Poland, which were identified in previous findings based on pollen analyses and phenotypic traits. The research material was the progeny of twelve beech provenances. The genetic structure of the populations was determined by ten highly variable microsatellite DNA loci. The results confirmed the high genetic diversity of beech at the north-eastern edge of its natural distribution, which infers the probability of their good adaptation to the changing climate and an extension of the range. Genetic analyses confirmed the existence of two genetic centres for beech in Poland. The populations from south-eastern Poland had a slightly higher diversity than the populations from the north-western area, which may indicate that the colonisation of Poland occurred by two routes. The results are important for creating the borders of the provenance regions and for limiting the transfer of seeds and seedlings. The choice of forest reproductive material, based on the knowledge of genetic diversity, is very important for the stability of future forests.


Introduction
The European beech (Fagus sylvatica L.) is one of the most widely found broadleaved tree species in southern, central and western Europe.Its range is from Sicily in the south to southern Scandinavia in the north and from northeastern Spain in the west to the Carpathian Mountains in the east (Hulten and Fries 1986).The north eastern edge of the natural distribution of beech is in Poland (Latałowa et al. 2004).The great adaptive potential of the European beech enables it to flourish in a wide range of natural conditions and so it plays an important economic and ecological role.In view of climate change, the study of genetic variation and differentiation of beech populations is crucial for determining the adaptive potential of the species.This is especially important for the populations situated at the edges of the natural distribution, which in light of the predicted climate change, will play a significant role in the dispersal of the species in the north-easterly direction (Matuszkiewicz 2001;Latałowa 2004).Identification of the genetic and demographic processes is essential for determining priorities in the management of genetic resources (Hampe and Petit 2005), including the legal regulations associated with the limits of transfer of forest reproductive material (FRM).The application of the knowledge gained will be of essential importance for the health and future of forests.
The European distribution model of the genetic diversity of trees illustrates the populations' evolutionary history (Leonardi and Menozzi 1995).The observed variation reflects the migration processes, including the inter-mixing of populations with a different evolutionary history.To correctly interpret the genetic diversity of the population of forest trees, inter alia, the number of refugial areas of the species and the routes of its post-glacial migration (Widmer and Lexer 2001) should be taken into account.The original model of the distribution of genetic diversity assumed that it was a result of the long-term isolation of populations during the last glacial period.The refugial areas, located in the south of Europe, were perceived as the places storing genetic richness, which was decreasing with the increasing distance from the distribution centres of a particular species, so called "southern richness and northern purity" (Hewitt 1999).Currently the refugial areas located in the south of Europe are considered as the sources of allele richness, but not as sources of genetic diversity (Widmer and Lexer 2001;Petit et al. 2003).According to new opinions in Central and Western Europe there could be an unrecognised local microrefugia, constituting, to date, an underestimated source of variation (Willis 2000;Mitka et al. 2014).The recent research shows that the history of post-glacial migration of many forest trees is much more complex than the explanation provided by previous hypotheses (Stewart and Lister 2001;Keir et al. 2011;De Lafontaine et al. 2013).Frequently, the process of species migration was controlled by the local conditions and therefore its dynamics in a region varied (Ralska-Jasiewiczowa et al. 2003;Magri et al. 2006).
For beech the most important period in its recolonisation history was the last interglacial era (Magri et al. 2006).Beech entered Poland in the subboreal period, first from the south and southwest into the Carpathian region -probably via the Beskid Niski Mountains and the Moravian Gate (Szafer 1932;Latałowa 1992).The second expansion route was from the northwest through the lakes of Western Pomerania, from where beech spread further to the east and the south (Ralska-Jasiewiczowa 1983;Giesecke et al. 2007;Magri et al. 2008).
As pollen data show there are clearly two centres of this species in Poland, one in the north-west and the second in the south-east (Latałowa et al. 2004).The differences between the capacity of beech from the northwest and southeast of Poland to adapt and grow have been confirmed by many years of research on interspecies variability in provenance experiments (Barzdajan and Rzeźnik 2002).
This paper presents the results of nuclear-microsatellite analyses of beech populations from the two Polish regions with a focus on the following questions: Is Fagus sylvatica in Poland structured into a northern and a southern genetic group?Is the genetic diversity of beech in the two regions significantly different and is this difference genetically justified?
Does Fagus sylvatica at the north-eastern edge of its natural range exhibit high or low genetic diversity?

Material and methods
We investigated the progenies of twelve beech provenances from northern and southern Poland (Table 1 and Fig. 1) in an experimental plot in Krynica (49°24´N, 20°55´E).The experimental plot was established in 1995 with provenances of beech from the limit of its natural range in Poland.
The collection of seeds to establish the project was made to achieve a minimum of 50 trees per stand.The detailed characteristics of the experimental plot and the provenances were presented in an earlier paper (Sabor and Żuchowska 2001).Leaves for DNA analysis were collected from approximately 46 trees per provenances.The trees were 22 years old.

Data analysis
The frequencies of null alleles (non-amplifiable fragments) per locus and per population were estimated using the Microchecker software (Van Oosterhout et al. 2004).The adjusted frequencies of null alleles according to Van Oosterhout et al. (2006) were used for further estimating the genetic variation.Genetic diversity within populations was estimated on the basis of the mean number of alleles per locus (N A ), the effective number of alleles (N E ), the observed and the expected heterozygosity (H O and H E ).To calculate these parameters the software GenAlEx 6.41 (Peakall and Smouse 2006) was used.FSTAT 2.9.3.2 software (Goudet 1995) was used to calculate the inbreeding coefficients (F IS ) for each locus.The significance was tested by performing 10 000 randomisations of alleles among individuals within the samples.The genetic differentiation patterns among populations were detected by calculating the populations pairwise F ST (Weir and Cockerham 1984).The statistical significance was tested by using 10 000 permutations of individuals between the compared populations.To quantify the distribution of total variance at different hierarchical levels, an analysis of the molecular variance (AMOVA) was carried out using GenAlEx 6.41 (Peakall and Smouse 2006).The same software was used for a Mantel test, to investigate the relationship between the pairwise population FST values and geographic distance.The correlation between the geographic coordinates (longitude, latitude) and the genetic diversity parameters were studied.A principal coordinate analysis (PCoA) was used to group populations according to their geographical origin based on pairwise F ST values (Hartl and Clark 1997).The Bayesian clustering method implemented in the software STRUCTURE 2.3.4 (Pritchard et al. 2000;Falush et al. 2003;Falush et al. 2007;Hubisz et al. 2009) was used to detect the genetic differentiation of populations.The multilocus SSR data of all individuals were used to calculate posteriori probabilities of membership and to define the K subpopulations based on their genetic similarity.The admixture model and correlated allele frequencies were used with no prior information.Five independent runs, each with 50 000 MCMC iterations after 10 000 burn-in periods, were carried out for the K set between 1 and 15.The ∆K statistics (Evanno et al. 2005) were used to detect the uppermost hierarchical level of the population structure, based on the rate of change between successive K values.The best estimate of K (∆K) was calculated using the web-based STRUCTURE HARVESTER program (Earl and von Holdt 2012).After identifying the optimal K by this method, samples were placed into the cluster for which they showed the highest percentage of contribution.

Results
All ten microsatellite loci used in this study were polymorphic and displayed a high number of alleles and a wide range of PCR products.A total of 132 alleles were detected.
The total number of alleles per locus ranged from 5 (FS3-04) to 25 (mfc5), with an average of 13.2 (Table 2).The highest number of alleles (mean value, N A ) and the highest effective number of alleles (N E ) were observed for the locus mfc5 with 16.1 and 8.09, respectively.The observed heterozygosity (H O ) ranged from 0.347 (Fs3-04) to 0.837 (csolf31), with an average of 0.698.The mean expected heterozygosity (H E ) was 0.701, and ranged from 0.346 (Fs3-04) to 0.866 (mfc5).In all analysed populations private alleles were found.The frequencies of these alleles were low, ranging between 0.01 and 0.13, with an average of 0.03.The values of the inbreeding coefficients (F IS ) were near zero in most cases.This indicates that there was no excess of homo-and heterozygotes in comparison to the expected values.However, significant and high F IS values were observed for two loci, indicating a deficiency of heterozygotes (Fs3-04 and mfc5).The genetic differentiation (F ST) between populations averaged 0.041.It varied from 0.030 to 0.087 at locus FS3-04, a very high value compared to the other nine loci.For the 12 provenances the mean number of alleles (N A ) varies from 7.7 in provenance 8 Jelenino to 9.2 in provenance 39 Bukowa with a mean of 8.3 (Table 3).The effective number of alleles (N E ) ranges from 3.52 (37 Marynin) to 4.88 (31 Pokrzywna).Generally, slightly higher values of N A and N E were found in the southern provenances.The observed heterozygosity (H O ) was the highest (0.750) in population 39 Bukowa and the lowest (0.673) in population 37 Marynin with the mean of 0.701.The expected heterozygosity (H E ) ranged from 0.665 (37 Marynin) to 0.730 (39 Bukowa).The mean expected heterozygosity was 0.704.Most provenances showed a lack of significant deviations from the expected Hardy-Weinberg frequencies.The mean value of the fixation index was 0.004.In most populations there were more homozygous genotypes than expected (F positive).Provenances 34 Jeleniów, 37 Marynin and 39 Bukowa, all from the south, showed significantly higher heterozygosity than expected.
The 'Analysis of Molecular Variance' (AMOVA), that allows the partitioning of variation among and within populations, indicated that variance was highest within populations (93%) (Table 4).A variation of about 5% of among populations occurred within regions, while only 2% variation was found between the regions.The Mantel test (Fig. 2) showed a positive and significant correlation coefficient between the geographic distance and pairwise F ST values (R = 0.721, p = 0.01).The existence of a distinct geographical structure of the studied populations was indicated.Significant, but negative, correlation with latitude was found only for the mean number of alleles N A (R = -0.812,p < 0.01).A longitudinal trend of genetic variation was not detected.
The results of the principal coordinate analysis (PCoA) based on pairwise F ST are presented in Fig. 3.The first two coordinates explain 47% of the total variation.Coordinate 1 indicates the separation of the northern and southern beech populations.One population from the south (31 Pokrzywna) was included in the group of northern populations.The spatial differentiation into a northern and a southern group was confirmed by the STRUCTURE analysis.
All loci were jointly used for calculating K subpopulations based on the membership proportions.The rate of change of data posterior probability (lnP(D)) between successive runs was the highest for K = 3 (Fig. 4).The clusters were not clearly distinguished.Admixture of clusters in the south and north provenances was high, but a trend of differentiation could be detected.The south populations are more represented by the green cluster and the north populations by the blue cluster.The average share of the green cluster in the southern provenances was 55%, whereas in  the north it was only 22%.For the northern populations the blue cluster averaged 43%, whereas in the south only 19%.The red cluster was characteristic for all provenances, averaging 34%.Two subpopulations were indicated by the STRUCTURE analysis, corresponding to the two groups of populations given by the PCoA analysis.The populations from the south form a clearly differentiated group.Only the 31 Pokrzywna is similar to the populations in the north.3).

Discussion
In this research the genetic diversity of twelve provenances of Fagus sylvatica from north-eastern part of its natural distribution range in Poland was studied.The analysis showed a clear genetic differentiation between provenances.The populations in the south-east had slightly higher variability (diversity, heterozygosity) than the populations in the north-west.In general the genetic diversity of the investigated populations was similar to populations of Fagus sylvatica from Central Europe, determined using the same marker type (nuclear microsatellites) Vornam et al. (2004): He = 0.765, Bilela et al. (2012): He = 0.740-0.792,Piotti et al. (2012): He = 0.728-0.810,and Dounavi et al. (2016) for EST microsatellites (He = 0.564-0.701)for German beech provenances.Lower values of genetic diversity (He = 0.617 resp.He = 0.614-0.632)for some German beech provenances have been reported Seifert (2011) and Rajendra et al. (2014).Differences between the observed and expected values of heterozygosity in the studied populations were small.According to Hartl and Clark (1997) this is a typical finding for a genetically stable population.
An equally low value of the fixation index was reported for beech by Lander et al. ( 2011) (F IS = 0.01-0.064)and Rajendra et al. (2014) (F IS = -0.003).However, generally higher F IS values have been reported in numerous isozyme studies for beech Comps et al. (1990): F IS = 0.065, Comps et al. (1991): F IS = 0.115, Leonardi and Menozzi (1995): F IS = 0.117, Harter et al. (2015): F IS = 0.02.The same applies to some nuclear microsatellite markers (Buiteveld et al. 2007: F IS = 0.224, Piotti et al. 2012: F IS = 0.139-0.233,Paffetti et al. 2012: F IS = 0.239).The positive value of fixation indices has been explained by a non-random micro-spatial structure due to mating preferences between neighbouring individuals (Comps et al., 1990;Merzeau et al. 1994;Shanjani et al. 2011).It is obvious that the values of the inbreeding coefficient depend, to a large extent, also on the selection of analysed loci.Thus, for example, microsatellite markers with a high frequency of zero alleles increase the F IS values (Buiteveld et al. 2007;Rajendra et al. 2014).In the reported study the mean F IS was close to zero, negative for the southern populations and slightly higher and positive for the northern populations.An increase of the inbreeding coefficient (F IS ) in the northern direction has also been reported by Comps et al. (2001) and Rajendra et al. (2014).
In the study most of the genetic diversity was attributable to the stands themselves.There was a small, but significant differentiation of 5% between stands.This low differentiation between stands is in accordance with other studies on Fagus sylvatica using both isoenzymes and DNA markers (Leonardi and Menozzi 1995;Larsen 1996;Konnert and Henkel 1997;Comps et al. 2001;Gömöry et al. 2003;Buiteveld et al. 2007;Kraj and Sztorc 2009).A low, but statistically important (p = 0.01), 2% genetic differentiation between regions was found, which reflects the geographical patterns of genetic variation.Genetic analyses confirmed the existence of two genetic centres for beech in Poland, one in the northwest and another in the southeast.Thus the reported analysis confirms previous findings based on pollen analyses and phenotypic traits (adaptation features) (Barzdajan and Rzeźnik 2002;Latałowa et al. 2004).A similar, significant geographic division of populations, at similarly low F ST values , was discovered earlier for beech in Italy (Leonardi and Menozzi 1995) and for beech from Austria, France, Germany, Italy and the Netherlands (Buiteveld et al. 2007).
Additionally, a continuous (clinal) increase in genetic diversity with increasing geographical separation was found, which confirms the work of Hampe and Petit (2005).The existence of a correlation between the genetic distance and geographical distance may provide evidence for a relationship between the genetic diversity and selective environmental gradients.This was the case in other studies using microsatellites with adaptive relevance (EST microsatellites) (Dounavi et al. 2016).However, it could also be a reflection of the hypothesis of genetic diversity based on different colonisation routes (Eckert et al. 2008).According to isoenzymes and DNA studies, a distinct variation in beech is found in the populations from various refugia, or from the various recolonisation routes after the ice age (Demesure et al. 1996;Gömöry et al. 1999;Comps et al. 2001;Magri et al. 2006).Originally it was assumed that beech in Central and Northern Europe arrived from the Balkan refugia (Huntley and Birks 1983;Taberlet et al. 1998).Currently the dominant opinion is that beech came to Central and Northern Europe from the refugia in southern France, the eastern Alps -Slovenia and Istria and probably from southern Moravia and southern Bohemia (Magri et al. 2006).The populations which survived the last glaciations in the Balkans and Italy were not the sources of today's Central European beech (Magri et al. 2008).A significant role in the development of the diversity of European beech was played by the areas where the migration routes from different refugia interacted during the postglacial recolonisation (Comps et al. 2001;Widmer and Lexer 2001;Petit et al. 2003).Moreover, the latest paleoecological and phylogeographical research indicates the former existence of potential northern refugial areas, which survived the period of glaciations due to favourable climatic-habitat conditions (Willis and van Andel 2004;Tzedakis et al. 2013).So far, the role of such refugia in the analysis of the genetic variation within Europe has been underestimated (Willis and van Andel 2004;Svenning and Skov 2007).Western Slovakia and Northeastern Hungary are regarded as potential cryptorefugia for beech (Stewart and Lister 2001).Genetic differences between the northern and southern beech provenances found in this study may confirm that colonisation of Poland occurred by two routes.The first led from the Czech Republic, Moravia and the plains of Hungary around 3000 years BC, entering Poland from the south, through the Sudetes and Beskids Mountains (Ralska-Jasiewiczowa 1983;Ralska-Jasiewiczowa et al. 2003;Latałowa et al. 2004).The second led from Germany and colonised the north-western part of Poland (Ralska-Jasiewiczowa et al. 2003;Giesecke et al. 2007), from where it continued to the east and the south (Ralska-Jasiewiczowa 1983; Latałowa et al. 2004;Magri et al. 2006).A higher level of beech genetic diversity in south-eastern Poland may be evidence for the colonisation of this part of the country from the area of the Czech-Moravian lowlands and Hungary.The longer migration route of the northern beech populations led to their genetic depletion (Demesure et al. 1996;Comps et al. 2001).The evidence of this is the observed tendency of the decreasing genetic multiplicity, given as mean number of alleles (N A ), from south to north, which closely corresponds to the reports of Comps et al. (2001) and Sułkowska (2010).
Climate change with increasing temperatures and dry periods will have an increasing impact on future beech stands (Fang and Lechowicz 2006).Results from modelling under different climate change scenarios predict a northward drift of the species range (Thuiller et al. 2006;Bellard et al. 2012;Falk and Hempelmann 2013) and the extinction of populations and a reduction of the range at the southern edge (Geßler et al. 2006;Dounavi et al. 2016).The climatic conditions in the north-eastern part of the present distribution range will become more favourable for beech.Actually beech is still in a phase of expansion as its range is showing further enlargement (Huntley and Birks 1983;Latałowa et al. 2004).Almost all prognoses state that the natural distribution range of the species will extend to southern and central Scandinavia, the Baltic countries and the larger part of Belarus (Thuiller et al. 2006;Bolte et al. 2010;Kramer et al. 2010).Analyses of climatic data suggest that the present conditions in north-eastern Poland should not be a limiting factor for the natural spread of beech into these areas (Matuszkiewicz 2001;Latałowa et al. 2004;Giesecke et al. 2007;Augustaitis et al. 2015).Due to the influence of climate change, the ecological, social and economic functions served by forests will depend on maintaining healthy, productive forest ecosystems, adjusted for the habitat (Chmura et al. 2010).In the near future the pool of species will probably not change drastically, but the frequency of their appearance, their growing performance, their seed production and resistance to external factors will change (Williams et al. 2012).Thus, the key issue is the implementation of management systems to fully maintain the diversity, while simultaneously maintaining the future stability of tree stands (Bolte et al. 2007).Therefore, forest management faces new challenges.The search for beech ecotypes having high adaptive potential in extreme growing conditions is essential to combat the predicted climate change, particularly the expected decrease of precipitation (Dounavi et al. 2016).An understanding of the genetic mechanisms which govern the physiological response of trees is very important.Research combining ecological and physiological approaches with assessments of tree growth in natural conditions is recommended (Geßler et al. 2006).The stands at the southern edge of the natural distribution range of the species will probably suffer drought stress (Netzer et al. 2016).Consequently, it is important for the genetic richness of these populations to be known, as it is assumed that greater genetic diversity provides greater possibilities for adaptation.A similar relationship exists at the northern edges of the natural distribution ranges.The populations at the northern leading edges, which in favourable conditions will be able to colonise new areas, should have a rich gene pool.All studied provenances showed a high genetic diversity, which implies the probability of good adaptation to the changing climate and the range being extended.The high adaptive capability of the Polish populations to stressful conditions compared to that of the German populations has already been confirmed by research (Rose et al. 2009).
An understanding of the genetic diversity of beech stands is also essential for activities associated with regeneration.Both natural and artificial regeneration are related to good seed years and good seed quality, but also with the management regime.Natural regeneration seems to permit an unrestricted flow of genetic information and thus the dynamic conservation of genetic multiplicity (Behm and Konnert 1990).For artificial regeneration, the choice of the appropriate FRM is of major importance.It is important to select provenances, not only on the basis of their adaptation value, but also on their genetic variability to provide the stability, needed for the generation of future forests.The present results complement the knowledge gained by multi-year observation of the variability of adaptation and growth features.They provide an insight into the genetic diversity and differentiation and demonstrate the high genetic variation of beech at the northern edge of its natural distribution.
The genetic differences between beeches from the two centres in Poland could result from their adaptation to local site conditions.The neutral molecular data could not be used as a surrogate for adaptive genetic information, but they provide partial insight into local adaptation or evolutionary potential (Kirk and Freeland 2011).The adaptation may lead to the genetic and phenotypic structuring of the population within the species, depending on the method of management over a few generations (Kramer et al. 2010).The local diversity of the beech population is frequently undervalued (Bolte et al. 2007) and without any doubt the maintenance of local diversity is important.However, to combat climate change a modest enrichment of populations from colder regions with provenances from more southern areas (Bilela et al. 2012) can be used to increase diversity and buffer against future climate uncertainty (Aitken and Bemmel 2015).To achieve this in a controlled way, the differences between the southern and northern provenances, indicated in the paper, should be maintained by the implementation of appropriate legal regulations.

Fig. 3 .
Fig. 3. Principal Coordinates Analysis (PCoA) of pairwise FST of twelve beech populations.Diamonds represent provenances from the north, crosses provenances from the south.

Fig. 4 .
Fig. 4. Results for STRUCTURE analysis: a) values of the statistic "Delta K" for all calculated K-values from 1-15 based on the combined dataset of nuclear microsatellites.The highest "Delta K" corresponds to the uppermost level of hierarchy among runs (Evanno et al. 2005).b) proportion of membership of each individual to three assumed subpopulations (K = 3).c) proportion of membership of each individual to three assumed subpopulations (K = 3).(Northern provenances: 3 to 19; Southern provenances 31 to 42.For location names see Table3).

Table 1 .
Description of investigated provenances.Geographical distribution of the studied provenances.The black circles indicate the location of the beech provenances, see Table1.The dotted line indicates the limit of the beech range in Poland.

Table 2 .
Summary statistics for the 10 microsatellite loci.
E -effective number of alleles, H O -observed hetrozygosity, H E -expected hetrozygosity, F IS -inbreeding coefficient.Significance of F IS values: *** p < 0.001

Table 3 .
Diversity measures in Fagus sylvatica.A -mean number of alleles, N E -effective number of alleles, H O -observed heterozygosity, H E -expected heterozygosity, F -fixation index.Significance of F values: * p < 0.01, number of private alleles and their frequency (Freq). N

Table 4 .
Hierarchical analysis of molecular variance (AMOVA) based on allelic distance matrix.
df -degree of freedom, SS -Sum of squares, MS -Mean squares, P -probability Silva Fennica vol.50 no.4 article id 1663 • Kempf et al. • Distribution of genetic diversity in Fagus sylvatica at…