Genetic Diversity and Differentiation of Pedunculate Oak ( Quercus robur L.) Populations at the Southern Margin of Its Distribution Range—Implications for Conservation

: Understanding intraspeciﬁc genetic variation is one of the principal requirements for the evaluation of tree species capacity to cope with intensive climatic changes, as well as designing long-term conservation programs. Herein, we evaluated the genetic diversity and genetic structure of seven pedunculate oak (Quercus robur L.) populations, located at the southern margin of its distribution range on the Balkan Peninsula (Serbia). The objective of the study was to propose future in situ conservation measures aimed at protection of pedunculate oak adaptive and neutral genetic diversity at the species rear-edge. Genetic diversity and structure were estimated using twelve highly polymorphic simple sequence repeat (SSR) markers. The mean expected heterozygosity (H e ) was 0.769, allelic richness (A R ) 9.63, and private allelic richness (p AR ) 0.79, indicating high genetic diversity in the studied populations. Genetic differentiation among the populations was low (F st = 0.032). Structure analysis, the unweighted pair group method with arithmetic mean (UPGMA) showed the existence of two gene pools unrelated to the populations’ area of occurrence. Taking into consideration the results of the current study and previous conservation activities on the pedunculate oak in Serbia, as well as the importance of rear-edge populations in the long-term conservation of the species genetic diversity, we suggested establishing three additional gene conservation units for securing long-term sustainability of the species.


Introduction
Pedunculate oak (Quercus robur L.) is widely distributed across Europe, from the Northern Iberian Peninsula, South Italy, and the Balkan Peninsula and Turkey to Southern Scandinavia, and from Scotland to Asia through the Ural and Caucasus Mountains [1]. In Southern Europe, the Quercus species survived the last glacial period in micro-and macrorefugia, from which the species were spreading to the north and east depending on climate occurrences [2]. According to Temunović et al. [3], the margin of pedunculate oak continuous distribution in southeastern Europe is located in the lowlands of the Pannonian Basin and floodplains of large rivers. Furthermore, according to Muir et al. [4] results showed that although in Europe Q. robur and Q. petraea hybridize extensively, the two species are separate taxonomic units; and populations from Serbia were close to populations from Slovenia. In Serbia, the majority of pedunculate oak stands is situated in the northern part of the country, along with the Sava and Danube rivers, and in Central Serbia (e.g., separate populations) [5,6].
In the past decades, pedunculate oak forests in Europe have been subjected to numerous stress factors, such as low underground water level [7], drought stress [8], pest and diseases [9,10], and climate change [11]. As a consequence of their negative impact, trees' health status and growth rate have frequently deteriorated, sometimes resulting even in premature death [7,12]. Moreover, under the influence of intensive climate change, pedunculate oak forests may suffer significant damage from interacting causal agents, which have been shown to trigger a complex phenomenon of oak decline [13,14].
Genetic diversity is the basis of biodiversity, and it is very important for the adaptation, survival, and evolution of all living organisms in a changing climate [15,16]. Peripheral populations are defined as populations on the edge of species distribution that usually represent reservoirs of biodiversity [17], as well as a source of drought-resistant ecotypes [18]. For example, Erichsen et al. [19] reported notably higher genetic diversity of Fraxinus excelsior L. growing in the Hyrcanian forests (i.e., species rear-edge) than in European populations. In addition, Fagus sylvatica L. populations at the rear-edge showed higher resistance to xylem embolism in comparison to leading edge populations, indicating the potential of this species to adapt to different climatic conditions [20]. In this manner, peripheral populations at the rear-edge are responsible for long-term persistence of the species [21] and could save the species from extinction [22]. However, several studies have shown a negative effect of recent climate change on these populations. Indeed, growing under stronger selection pressures than core populations, rear-edge populations are frequently subjected to various biotic and abiotic stress factors which negatively affect their reproductive success and survival [23,24]. The problem becomes even more pronounced when these stresses act simultaneously or sequentially, chaining a large pool of negative factors. Therefore, peripheral populations have become the subject of conservation of genetic resources [25] and they were used in afforestation of other areas [26].
Molecular markers such as simple sequence repeats (SSRs) have a high rate of polymorphism, good reproducibility and codominant inheritance, making them ideal markers for studying genetic diversity [27]. Therefore, polymorphic SSR markers have proven to be a reliable tool for analyzing genetic diversity and population genetic structure in the case of many oak species [28][29][30]. Previous studies of the variability of pedunculate oak populations in Serbia were based on certain morphological and anatomical traits variation [5,6,31], failing to provide information on neutral genetic differences. Considering no consistent correlation between genetic variation and phenotypic variability [32], we evaluated the neutral genetic diversity of oak populations in Serbia, for the first time.
Bearing in mind the foregoing discussion, the aim of this paper was twofold: (1) to estimate the genetic diversity and structure of pedunculate oak populations in Serbia, and (2) to identify priority populations for conservation in order to preserve genetic diversity at the rear-edge of the species distribution range.

Plant Material and DNA Isolation
This study included seven natural pedunculate oak populations from Serbia (Table 1; Figure 1). In total, 210 adult trees (30 per population) were analyzed. Buds and flushing leaves were collected from trees located 50-100 m from each other. DNA was extracted from fresh plant material which was stored in deep freeze (−70.0 • C) before extraction. Total DNA extraction was carried out according to the DNeasy Plant Mini Kit protocol (Qiagen ® , Hilden, Germany). DNA concentration was measured with BioSpec-nano (Shimadzu, Kyoto, Japan) and standardized to 10 ng µL −1 .
For SSR analysis, reactions were performed in a total volume of 15 µL for each marker separately containing 10 ng µL −1 of template DNA, 1× reaction buffer (Promega GoTaq G2 Flexi, 5× reaction buffer with no magnesium), 2 mM of MgCl 2 (Promega) in the case of QrZAG11, 30, 87, 96, 101, 112, QpZAG110 and MSQ13, while 3 mM of MgCl 2 in the case of QpZAG9, QpZAG1/5, QrZAG15 and 20, 12 µM of each dNTP (Promega dNTP mix, 10 mM each), 0.5 µM each of the reverse and forward primer, except MSQ13 where 0.34 µM each (Applied Biosystems Custom DNA oligos) and 0.5 unit polymerase (Promega GoTaq G2 Flexi, 5U/µL). Fluorescent dyes compatible with the Applied Biosystems G5 Matrix were used on the forward primer 5 end. In order to multiplex the samples of different SSR loci by size and by dye, the following two mix combinations were set up: mix 1 contains the following six markers: QrZAG112 (NED), QrZAG96 (NED), QrZAG30 (VIC), QrZAG11 (VIC), QrZAG87 (PET), and QrZAG20(FAM)), while mix 2 contains the remaining six markers: QpZAG110 (NED), QpZAG1/5 (VIC), QpZAG9 (PET), MSQ13 (6FAM), QrZAG101 (NED), and QrZAG15 (6FAM). PCR amplifications were performed in a Veriti 96-well Thermal Cycler (Applied Biosystems, Waltham, Massachusetts, USA) in the case of QrZAG11, 96, 112, QpZAG110 with the following steps: an initial denaturation for 3 min at 94 • C, followed by 30 cycles of 30 s at 94 • C, 30 s at 50 • C, 90 s at 65 • C, and finished with a final extension step at 65 • C for 15 min. In the case of QrZAG30, 87 and 101, a short protocol was applied following the modification of Neophytou et al. [36] that included an initial denaturation step 8 min at 94 • C, followed by 23 cycles of 15 s at 94 • C, 15 s at the optimum annealing temperature (T a = 50 • C for QrZAG101 and T a = 56 • C for QrZAG30 and 87), 15 s at 72 • C, and 10 additional cycles with reduced denaturation temperature (89 • C), without a final elongation step. The original, long protocol was more appropriate in case of the remaining markers, so the following PCR programming was used: an initial denaturation for 3 min at 95 • C, followed by 35 cycles of 50 s at 95 • C, 50 s at the optimum annealing temperature (T a = 52 • C for MSQ13 and T a = 65 • C for QpZAG9 and T a = 55 • C for QpZAG1/5, QrZAG15 and 20), 105 s at 72 • C, finished with a final extension step at 72 • C for 10 min. SSR genotyping was performed on an ABI 3730 DNA Analyzer (Applied Biosystems, Waltham, MA, USA) by BIOMI Ltd. (Gödöllő, Hungary). Allele scoring was carried out using PeakScanner 1.0 software (Applied Biosystems, Waltham, MA, USA).

Variation within Population
In order to describe variation among populations, effective numbers of alleles (N e ), observed heterozygosity (H o ), expected heterozygosity (H e ), and Shannon's Information Index (I) software POPGENE ver. [37] was used. The number of alleles in a genotype (allelic richness (A R ) and private allelic richness (p AR )) was calculated using HP-RARE software [38,39]. The inbreeding coefficient (F is ) was calculated for each locus with FSTAT software [40]. Frequency of null alleles was calculated with ML-Null software [41]. Microchecker software was used to test presence of null alleles, long allele dropout, and scoring errors [42].

Variation among Populations
In order to estimate among-population genetic variation, we calculated F st using FreeNA [43]. Means and significant values over the loci and populations were obtained by bootstrapping. GENEPOP was used to test genotypic disequilibrium between all pairs of loci, and Bonferroni correction was applied. The hierarchical distribution of genetic variation (AMOVA) [44] among and within the populations for twelve markers was assessed with GenAlEx software [45]. GENEPOP was also used to estimate the number of migrants (N m ) based on the private allele method. FreeNA software [43] was used to estimate Nei's unbiased genetic distance and significance among populations [46]. We also estimated the Cavall-Sforza and Edwards genetic distance for each pair of populations [43], with the ENA method (DC ENA ) and the INA method (DC INA ). The dendrogram among populations was produced with UPGMA based on Nei's genetic distance using POP-TREE2 software [47], and 1000 repetitions of bootstrapping were used.
STRUCTURE software was used to infer genetic structure and define the number of clusters (gene pool) in a dataset based on the Bayesian clustering method [48]. The most likely number of clusters (K) was defined by Structure Harvester [49]. For STRUCTURE analysis the following settings were used: burn-in 100000, MCMC 1000000, 20 repetitions per K. To detect structure among populations that are likely to be similar due to migration or shared ancestry, we used admixture model [50]. In order to test Cavalli-Sforza and Edward's chord, FreeNA software was used [43].

Variation and Genetic Diversity within Populations
In total, 253 alleles were obtained from the twelve SSR loci used in this study (Table 3). Several SSRs, QrZAG87, QrZAG30, QrZAG112, QrZAG101, and QrZAG11, contained more than 20 alleles. The number of alleles per locus (Na) ranged from a minimum of 11 (MSQ13) to a maximum of 36 (QrZAG30 and QrZAG11). Effective number of alleles (N e ) varied between 1.77 (QrZAG96) and 12.88 (QrZAG30), with an average value of 11.46 per locus. Expected heterozygosity (H e ) within population for all SSR loci studied was 0.77. F is varied between 0.021 (QrZAG87) and 0.453 (QrZAG112), with an average value of 0.209 alleles per locus. The genetic differentiation (F st ) was 0.032; the highest one (0.065) was for loci QrZAG96, whereas the lowest was for loci QrZAG20 and QpZAG1/5 (0.017). The mean frequency of null alleles (f n ) was 0.082 and ranged from 0.004 for locus QrZAG87, to 0.197 for locus QrZAG112 (Table 3).

Genetic Variation between and within Populations
The effective number of alleles (N e ) ranged between 5.642 in the population OD and 6.429 in the population AP, with an average of effective 6.07 alleles per individual population. The expected heterozygosity (H e ) varied between 0.740 (KL) and 0.793 (KU). The overall average allelic richness (A R ) was 9.63; the highest one (10.47) was distributed in KU, whereas the lowest was in OD population (8.69). Private allelic richness (p AR ) varied between 0.41 (MO) and 1.12 (KU) ( Table 4). The genetic differentiation between populations was low (F st = 0.032, p > 0.05). AMOVA test proved that differentiation among populations was very low (2%), and 98% of the variation was within-population level, and statistically significant (p < 0.001, Table 5). Furthermore, the level of gene flow (N m ) was measured to be 9.82 individuals per generation between the populations. The results from Structure Harvester [49] revealed the existence of two pedunculate oak gene pools in Serbia (K = 2) (Figure 2).
The results of STRUCTURE revealed two pedunculate oak gene pools in Serbia The UPGMA dendrogram, based on Nei's [46] unbiased genetic distance matrix, was created to depict the genetic relationship among the seven studied populations applying MEGA software [53]. The populations were clustered into two distinct groups. Group I consisted of three populations, in which OD and KL are more similar to each other, and MO slightly separated from the previous two populations. The populations KU and KG belong to Group II together with the population in AP. This appears genetically differentiated from the other two, whereas Group II also included the most northern population BM that forms a completely separated branch in the same cluster ( Figure 4). The results from Structure Harvester [49] revealed the existence of two pedunc oak gene pools in Serbia (K = 2) (Figure 2).   The results from Structure Harvester [49] revealed the existence of two pedunculate oak gene pools in Serbia (K = 2) (Figure 2).    MEGA software [53]. The populations were clustered into two distinct groups. Group I consisted of three populations, in which OD and KL are more similar to each other, and MO slightly separated from the previous two populations. The populations KU and KG belong to Group II together with the population in AP. This appears genetically differentiated from the other two, whereas Group II also included the most northern population BM that forms a completely separated branch in the same cluster ( Figure 4).

Discussion
Genetic diversity has a key role for forest tree population adaptation under changing climate, as well as planning future conservation programs. We found that Na ranged between 11 (MSQ13) and 36 (QrZAG30 and QrZAG11), with a mean value of 29.99. We observed higher Na than the one that was observed by Vranckx et al. [54] (Na = 21.26) considering several mutual loci (QpZAG9, QpZAG15, QpZAG110, MSQ13, and QrZAG112). Likewise, the Na on the QrZAG11 (36) was notably higher than previously reported for

Discussion
Genetic diversity has a key role for forest tree population adaptation under changing climate, as well as planning future conservation programs. We found that N a ranged between 11 (MSQ13) and 36 (QrZAG30 and QrZAG11), with a mean value of 29.99. We observed higher N a than the one that was observed by Vranckx et al. [54] (N a = 21.26) considering several mutual loci (QpZAG9, QpZAG15, QpZAG110, MSQ13, and QrZAG112). Likewise, the N a on the QrZAG11 (36) was notably higher than previously reported for this species (i.e., Craciunesc et al. [55] and Pohjanmies et al. [56] found 18 and 20 alleles for QrZAG11, respectively).
Expected heterozygosity (H e ) across the studied pedunculate oak populations was similar to other populations at the south-east margin. For example, marker QrZAG87 had similar value of H e (0.87) to those observed by Katičić-Bogdan et al. [57] (0.87), Morić [58] (0.88), and Muir and Schlötterer [59] (0.89). Based on QTL associated studies [60], locus QrZAG87 is considered to be linked to genome regions responsible for bud burst timing. Therefore, considerably high H e value of this locus is probably the consequence of sampling procedure, since the material of both late (var. tardiflora) and early (var. praecox) flushing individuals were collected. According to Neophytou et al. [36], the same locus has higher contribution to intraspecific than interspecific provenance discrimination. Considering locus QrZAG96, the H e observed in our study (0.40) was similar to that reported in Croatia (0.42) [58] and higher in comparison to southern population in Greece (0.33) and Bulgaria (0.39) [36]. In contrast, marker QpZAG110 had smaller H e (0.62) in comparison to Greek (0.74) and Bulgarian (0.84) populations [36]. Marker QrZAG112 had equal H e value (0.81) in our and other studies at the southern margin [36,55,57,58].
The mean A R in our populations (9.63) was lower than A R measured in Belgium (11.5) considering four mutual loci (QpZAG 9, QrZAG15, QrZAG112, MSQ 13) [54]. Moreover, according to Katičić-Bogdan et al. [57], the mean A R in Croatian pedunculate oak populations (12.77) was higher than in the Serbian (9.63), considering eight mutual loci (QpZAG9, QrZAG11, QrZAG30, QrZAG87, QrZAG96, QrZAG101, QrZAG112, and MSQ13). Considering our results, lower A R in comparison to those observed by the aforementioned authors could be explained by the founder effect impact [61]. It has been empirically shown that A R is more sensitive to the founder effects than H e , meaning that the loss of rare alleles is likely to reduce A R to the larger extent than H e (a rare allele which is lost will probably not reduce H e as much as A R ) [62,63]. In our case, it is possible that the measures on stands regeneration were not synchronized with the pedunculate oak seed years (as is the case with modern shelterwood system); therefore, the studied populations originate from the seeds collected from a limited number of mother trees.
The observed value of F st (0.032) in Serbian pedunculate oak populations was lower than the ones previously reported for Bosnia and Herzegovina (0.051) [64] and Romania (0.052) [55], respectively. On the other hand, genetic differentiation was higher in Serbian populations than in other central European populations (0.024) [65]. Generally, peripheral populations have a tendency to show higher genetic differentiation than populations from the center of the species distribution [55,56,64]. For example, Jiménez et al. [66], reported higher genetic variation in marginal than central populations of Quercus suber L. Likewise, analysis of molecular variance (AMOVA) showed greater variability within (98%) than among populations (2%), which is commonly observed in geographically widespread tree species and populations, and mostly attributed to natural selection, different ecological conditions, and anthropogenic factor effects [6].
The, UMPGA and STRUCTURE analyses revealed that the studied populations were grouped into two different gene pools. Both gene pools consisted of geographically distinct populations; i.e., the first gene pool included populations MO, OD, and KL, whereas the second group consisted of geographically scattered populations KU, AP, KG (southernmost population), and BM (the northernmost). Low correlation between geographical and genetic distances could be explained by the anthropogenic impact on the current genetic structure of pedunculate oak, as well as human-assisted regeneration of these forests. Namely, situated on the border between the Austro-Hungarian and Ottoman empires, the northern part of Serbia (formerly part of Austro-Hungarian Monarchy; until 1918) has been subjected to turbulent historical occurrences over the centuries that frequently lead to intensive deforestation, and in some cases even forest devastation [67]. The first interim guideline for forest management (K.u.K.Őstereichischie-Ungarische Monarchie) was enacted in 1754, after which the legal regulations related to forest management changed several times depending on historical circumstances during 18th and 19th centuries [68]. The beginnings of modern management of pedunculate oak forests are linked to 1860, when the Forest Law of Austro-Hungarian Monarchy (passed in 1852) was extended to this region [69]. Nevertheless, although historical data on forest management in this area date back to the middle of the 18th century, there is no clear evidence on the origin of acorns used for regeneration of pedunculate oak forests. Likewise, although there is no evidence that populations were coppiced in the past, we assumed that coppicing was the simplest form of usage and forest reproduction, especially around settlements. This assumption can be supported by the fact that due to the excessive use of forests around settlements, the first Decree Law on forests has passed on 1755 to regulate the future management and protection of the forests in this region [69]. In contrast, certain records highlight that pedunculate oak forests in the area were largely preserved until second part of XIX century, primarily due to low population density, existence of impassable swamps on the large areas and lack of industrial wood processing, etc. (see [69]).
Even though some authors stated that the seed of Slavonian pedunculate oak (Quercus robur subsp. slavonica) was broadly used across Austro-Hungarian Monarchy [69,70], certain records also stressed that the afforestation of the nearby areas was done using seed originating from the Czech Republic and Slovakia [71]. In addition, the populations KL and KG are 97 and 90 years old, respectively, meaning that they were regenerated after World War I, when the Austro-Hungarian Empire collapsed, and the Kingdom of Yugoslavia was proclaimed (e.g., Northern part of Serbia was annexed to Serbia). Thus, it is very likely that these populations are established from the acorns collected in some of already existing populations, especially considering the fact that genetically close populations (OD to KL and KU to KG) are also the oldest ones, and theoretically could produce seed in that time (see [72,73]). Another reason to consider human-assisted regeneration in shaping pedunculate oak genetic differentiation also lies in the tendency of oak periodical yield; e.g., according to Gradečki-Poštenjak et al. [74] abundant yield could be expected every 4 years. Thus, due to the problematic acorn sourcing outside mast years, a common practice even today is that seed for afforestation is provided from other refugia, which are sometimes tens or even hundreds of kilometers away. For example, all pedunculate oak forests in Northern Serbia are managed by the single public enterprise for forest management and belong to the same provenance region (see Figure 1), thereby the production, market, and use of seed within this region is in accordance with national laws on forests and forest reproductive material, which are fully compliant with EU Directive 1999/105/EC.

Implications for Conservation
The most common approach in conservation of genetic diversity of forest trees is establishment of in situ dynamic conservation units. The prime objective of dynamic genetic conservation is not only to capture the current genetic diversity, but also to support continued evolutionary processes within populations [75]. Considering climate change as the major challenge for long-term conservation of forest genetic resources, European Forest Genetic Resources Programme (EUFORGEN) made pan-European recommendations aimed primarily at preserving adaptability of tree populations in the future, emphasizing the following activities: (1) identification of most vulnerable tree populations and species, (2) conservation of genetic diversity in marginal populations, and (3) conservation of adaptive diversity [76].
Our study revealed that the studied pedunculate oak populations were characterized by a high level of genetic diversity (e.g., expected heterozygosity), which is likely for Southern European tree populations as they are located in areas of glacial refugials [77]. According to Reed and Frankham [78], heterozygosity is tightly related to the overall population fitness, and therefore important for its capacity to adapt to changing environmental conditions. However, as found in Southern Europe, near the edge of the species environmental limits, these populations are supposed to be negatively threatened by climate change, imposing the need for implementation of measures aimed at conservation of their genetic diversity. So far, only one gene conservation unit (GCU) of pedunculate oak has been established in Serbia (population MO), in the framework of the European Forest Genetic Resources Programme (EUFORGEN). However, our results showed that this population handles the lowest p AR , as well as quite low A R and H e values (e.g., only the populations OD and KL have lower AR and He, respectively), imposing the need to set aside additional conservation units. Considering the importance of A R and p AR for reliable diversity assessment and conservation purposes [30,65,79,80], we used these parameters as the most relevant criteria for prioritization of populations for conservation. The population that stood out by the highest values of H e (0.793), A R (10.47), and p AR (1.12) was population KU; therefore, we suggest the establishment of an additional GCU to cover its genetic variability.
Furthermore, to ensure conservation of sizable amount of adaptive diversity throughout species distribution range, EUFORGEN recommended that at least one GCU per country and per environmental zone should be identified [76]. Moreover, in the case of marginal populations, it is recommended to duplicate the number of GCUs within each "country x environmental zone". Bearing this in mind, we suggest establishing one GCU to cover the southernmost population KG, which is the only one situated southern of the Danube and Sava rivers, in the zone characterized by cool and moist climate (the Central Serbia region). Likewise, although the rest of the studied populations lie in the region with cool and dry climate, one more GCU should be established to cover the northernmost population BM for two reasons: (1) Establishing more GCUs along the transect from the North to the South will ensure the conservation of the species' spatial genetic diversity, (2) the population BM belongs to Special Nature Reserve "Upper Danube", which was officially designated UNESCO Biosphere Reserve in 2017 under the name "Bačko Podunavlje" [81], and will be part of the future transboundary Mura-Drava-Danube Biosphere Reserve.

Conclusions
In this study, twelve SSR markers were used to estimate the genetic diversity and genetic structure of seven pedunculate oak populations in Serbia, which are situated at the rear edge of the species distribution range. The pedunculate oak populations studied herein showed high genetic diversity, which was, as expected, more pronounced at the intrapopulation level than between different populations. Considering the results, we suggest a conservation strategy based on the establishment of additional in situ conservation units across the species distribution range in Serbia to cover all environmental zones, as well as populations characterized by exceptional genetic diversity. The evaluation of genetic diversity and structure of Q. robur populations at the rear edge allowed us to identify populations that are important for conservation. In addition, this genetic information could greatly improve knowledge of the species genetic resources which is available along the species distribution.  Institutional Review Board Statement: "Not applicable" for studies not involving humans or animals.
Informed Consent Statement: "Not applicable" for studies not involving humans.

Data Availability Statement:
The data presented in this study are available from the corresponing author upon request.