ISSUE 2/2021

Conservation of Quercus robur L. genetic resources in its south-eastern refugium using SSR marker system – a case study from Vojvodina province, Serbia

Article by Branislav Trudić, Evangelia Avramidou, Barbara Fussi, Charalambos Neophytou, Srdjan Stojnić, Andrej Pilipović | 05.07.2021 - 16:11

Schlüsselbegriffe: Stieleiche, DNA-Marker, genetische Vielfalt, Population

Abstract

Pedunculate oak, Quercus robur L., represents one of the most valuable tree species in Serbia and is mainly found in the Northern part of the country (Vojvodina province). These forests are among the most extensive pedunculate oak forests in the southern part of the distribution range in the Western Balkans. Serbian oak forests are under pressure due to climate change as well as affected by abiotic and biotic factors causing damage to stands since the 1950s. goal of current forest management is to ensure sustainability despite climate change and disturbances. The aim of this study was to estimate inter- and intra-population genetic variability of Q. robur using molecular markers (12 SSRs) and phenotypic traits to obtain the first genetic profile of the Pedunculate oak gene pool in Serbia. Our results indicated that individuals of pedunculate oak from Vojvodina province of Serbia showed higher expected heterozygosity (He) and number of migrants in comparison to other studies in Europe. The two sampled subpopulations of Q. robur had a highly diverse gene pool compared to other European populations. Suggestions for in-situ and ex-situ conservation were made in concordance with EUFORGEN’s technical guidelines on conservation and use of pedunculate oak in Europe.

Zusammenfassung

Die Stieleiche, Quercus robur L., ist eine der wertvollsten Waldbaumarten in Serbien, mit der Hauptverbreitung im Norden des Landes (Provinz Vojvodina). Auf dem West-Balkan gehören diese Wälder zu den ausgedehntesten im südlichen Teil des Verbreitungsgebiets dieser Art. Eichenwälder stehen durch den Klimawandel unter Stress und seit den 1950er-Jahren treten vermehrt biotische und abiotische Schäden an den Beständen auf. Ziel derzeitiger Waldbewirtschaftung ist es, die Nachhaltigkeit trotz Störungen und Klimawandel sicherzustellen. Ziel dieser Studie war es, die genetische Variabilität von Q. robur zwischen und innerhalb Populationen mithilfe molekularer Marker (12 SSRs) und phänotypischer Merkmale abzuschätzen, um das erste genetische Profil des Stieleichen Genpools in Serbien zu erstellen. Unsere Ergebnissen zeigten, dass Stieleiche Individuen aus der Provinz Vojvodina in Serbien im Vergleich zu anderen Studien in Europa eine höhere erwartete Heterozygotie (He) und Anzahl von Migranten aufweisen. Die beiden untersuchten Subpopulationen von Stieleiche weisen einen vielfältigeren Genpool im Vergleich zu anderen Populationen in Europa auf. Es werden Vorschläge für die In-situ- und Ex-situ-Generhaltung gemacht, die den technischen Richtlinien von EUFORGEN zur Erhaltung und Verwendung von Stieleiche in Europa folgen.

Introduction

Pedunculate oak (Quercus robur L., fam. Fagaceae) is an ecologically and economically important deciduous, noble hardwood species, which extends from western Asia to Europe (Trudić et al. 2013). Pedunculate oak forests represent one of the most valuable forests in the Republic of Serbia, especially in the complex alluvial hygrophilous forest types located in north-west area of Vojvodina Province (Banković et al. 2009). In recent years these ecosystems suffered from the oak decline syndrome, which is also found in several other European countries (Thomas 2008). Regardless of negative influences of various biotic and abiotic factors, it is evident that many pedunculate oak stands in the Republic of Serbia are characterized by a good phenotype, which was one reason for long history of collecting and exporting acorns to Western and Eastern Europe (Orlović et al. 2008). With the establishment of seed orchards, a large portion of the most valuable genetic material in natural populations of this species is preserved. Progeny tests, while still under development, could potentially represent another important assett for further selection of oak (Stojnić et al. 2014). A number of conservation measures have been applied over the last decades, new efforts are needed to keep pace of climate change in order to preserve these forests.

Forest genetic diversity is a key component to maintain species and ecosystem diversity under the threat of biotic, abiotic and climatic consequences (Fussi et al. 2016). Assessment of genetic diversity is further a fundamental step to accomplish the scope of Convention on Biological Diversity (CBD), which is to halt the ongoing erosion of biological variation of plant communities. Indeed, sustainable forest management practices affect diversity at various levels from ecosystem to species and from populations to individuals. Forest dynamics and the current situation shaped by past forest management practices must be reviewed and reconsidered in order to maintain future adaptation and secure resilience of forests (Kavaliauskas et al. 2018).

Flushing time of leaves determine the growing season length and is thus a major determinant of growth and survival (Dantec et al. 2015). Strong geographic patterns of variation in flushing times have been associated with elevation, latitude and longitude (Conover Duffy & Hice 2009, Kremer et al. 2014). Various studies on Quercus sp. showed that timing of flushing varies between trees from the same population (Pekeč et al. 2017) and has a genetic background (Perić et al. 2000, Alberto et al. 2011). For example, studies on sessile oak, Quercus petraea, showed considerable variability in flushing between populations (Vitasse et al. 2009a, b), due to both phenotypic plasticity (Vitasse et al. 2010) and genetic population differentiation (Vitasse et al. 2009c, Alberto et al. 2011).

Since the discovery of molecular markers for oak species (Steinkellner et al. 1997, Kampfer et al. 1998, Deguilloux et al. 2003, Guichoux et al. 2011), numerous studies have focused on evolutionary history and genetic structure comparisons between various oak species, especially considering Q. petraea (Streiff et al. 1997, Neophytou et al. 2010). A more comprehensive analysis from Neophytou et al. (2015) answered the question, whether human intervention played a significant role in the genetic diversity of three important oak species in Central Europe (Q. robur, Q. pubescens and Q. petraea). An important finding connecting bud burst and SSR markers was found by Scotti-Saintagne et al. (2004) in Q. robur. Specifically, locus QrZAG87 (which is situated in linkage group LG2 of Q. robur map), is associated with a QTL genome region responsible for bud burst (Scotti-Saintagne et al. 2004). In another study, Craciunesc et al. (2011) performed genetic diversity assessments based on isoenzymes and SSR markers on Q. robur samples from the Prejmer natural reserve in Romania, and found high genetic diversity and low inbreeding values. Evaluation of genetic diversity transmission from Q. robur adults to juvenile was also estimated from Vranckx et al. (2014) in an effort to examine if natural regeneration is an effective forest management practice to preserve genetic variation in small and low-density forest stands. A recent study for Q. robur in Croatia by Katičić et al. (2018) assessed the genetic diversity parameters based on SSR and chloroplast molecular markers in clonal seed orchards in order to examine the possibility to produce genetic improved material for future use. 

The aim of this study was to estimate inter- and intra-population genetic variability of Q. robur using genetic SSR molecular markers in order to describe the first genetic profile of the most valuable pedunculate oak gene pool in Serbia, which is found in the area of the Sava and Danube Rivers. Although a number of studies regarding its physiology, taxonomy and general biology (Kostić et al. 2019, Stojanović et al. 2015, Stojnić et al. 2019, Trudić et al. 2013, Vastag et al. 2020) were conducted in Serbia, this is the first genetic study using SSR molecular markers. The study has four objectives: 

1) to assess for the first time the genetic profile of pedunculate oak in Serbia; 

2) to compare the genetic variation of this gene pool with other European populations to estimate potential value of gene pool conservation; 

3) to compare late and early flushing individuals using genetic variation within subpopulations and 

4) to make recommendations for improvement of conservation practices in the southern part of the distribution range of pedunculate oak.

Material and methods

Trudic_Fig1.jpg

Figure 1: The map of two officially registered regions of pedunculate oak, Quercus robur, provenances in Serbia regulated by national guideline for establishing the regions of pedunculate oak provenances in Serbia: 91/2008-16. The brown color is the distribution range of the Q. robur in Europe (upper smaller picture) and Serbia (large picture), according to EUFORGEN (http://www.euforgen.org/species/quercus-robur/, accessed on 4th of November 2020). / Abbildung 1: Karte von zwei registrierten Herkunftsgebieten der Stieleiche, Quercus robur, in Serbien, geregelt durch die nationale Richtlinie zur Einrichtung der Herkunftsgebiete der Stieleiche in Serbien: 91/2008-16. Die braun markierte Fläche ist laut EUFORGEN (http://www.euforgen.org/species/quercus-robur/, vom November 2020) das Verbreitungsgebiet von Q. robur in Europa (kleines Bild rechts oben) und in Serbien (großes Bild).

In spring 2015, leaf samples from 83 individuals were collected from two subpopulations located in the Northern part of Serbia (i.e. Posavsko-podunavski region in autonomous province of Vojvodina), and managed by the Public Enterprise “Vojvodinašume”, Forest Estate “Sremska Mitrovica”, Forest Management Unit (FMU) “Morović” (Figure 1). The altitude of the management unit ranges from 81 to 83 m a.s.l. The geological substrate is composed of alluvial sand sediments forming mostly gley soils ranging from riparian black soil, black meadow soil to brown forest soil. Hydrological conditions are characterized by the absence of regular flooding and soil moisture is strongly dependent on fluctuations of groundwater table level. The investigated area lies in the Pannonian environmental zone PAN1 (Metzger et al. 2005). Average annual temperature for the investigated area is 10.8 °C with an average annual precipitation of 569.6 mm, of which more than 50% of precipitation during growing season. Data on ecological conditions of a study site (climate, hydrology and soils) were obtained from the forest management plans of the forest enterprise PE "Vojvodinašume".

42 samples from subpopulation 1, located in section 14 of management subunit 5a (equal number of late (var. tardiflora) and early (var. praecox) flushing individuals) and 41 (21 early flushing and 20 late flushing individuals) samples from neighboring subpopulation 2 located in sections 32 and 33 of management subunit 32f were collected. The distance between these two subpopulations is about 2 km. Phenological and morphometric traits assessment were used to evaluate individuals for upcoming genetic analysis: according to phenophase (early or late flushing) and physiological groups (vital or senescent). The main criterion in determining the affinity of trees to one of the physiological groups (phenotyping: vital or senescent) was the degree of canopy damage, which was assessed according to Dubravac et al. (2011); i.e. trees with canopy damage above 25% were considered significantly damaged trees, while trees with canopy damage up to 25% were considered undamaged or vital trees. Damage caused by pathogens was used as an additional criterion in the selection of trees. In both subpopulations, equal number of senescent and vital individuals was sampled. For all individual trees, morphometric measurements included diameter at breast height (DBH [cm]), absolute height of trees (H [m]) and height to the crown base (HC [m]). DBH was measured with a caliper (Haglöf, Sweden). H and HC were measured using “Vertex III” (Haglöf Sweden AB) device. Geographic coordinates for each individual were recorded by GPS in WGS84 (degrees, minutes, seconds) format for spatial analysis of genetic diversity. Samples were kept in 96% ethanol solution in 2 ml Eppendorf tubes and transported to the laboratory facilities of the Bavarian Office for Forest Genetics (AWG, Teisendorf, Germany). Total DNA extraction was performed with the ATMAB method (Dumolin et al. 1995). At the end, the DNA pellet was resolved in 50µl 1xTE buffer. DNA concentration and purity was measured by GeneQuant pro-spectrophotometer from Amersham Biosciences.

Trudic_Tab1.jpg

Table 1: Nuclear marker system consisted of twelve SSRs for Q. robur and Q. petraea. / Tabelle 1: Zwölf Kernmikrostellien-Genorte für Q. robur- und Q. petraea.

PCR conditions

In Table 1, our Quercus marker system and its characteristics are presented with followed sequences of all 12 used marker loci. Multiplex PCRs were performed with QIAGEN Multiplex PCR Kit (100) ® (for 100 x 50 µl multiplex PCR reactions: 2x QIAGEN Multiplex PCR Master Mix (providing a final concentration of 3 mM MgCl2, 3 x 0.85 ml), 5x Q-Solution (1 x 2.0 ml), RNase-Free Water (2 x 1.7 ml, with different concentrations of SSR which are presented in Table 1). PCR amplifications were performed in a Biometra Professional Termocycler (Analytik Jena, Jena, Germany) as follows: an initial step of 15 min at 95 °C, followed by 25 cycles for, each one including 30s at 94 °C for denaturation, 1 minute and 30s at 55/57 °C (depending on the multiplex set used; 55 °C was used for first four markers from Table 1 (MsQ13, QrZAG112, QrZAG96 and QrZAG11), while temperature for the other eight primers was 57 °C) for annealing, and 30s at 72 °C for elongation and a final extension of 30 min at 60 °C. PCR reactions were prepared according to manufacturer’s instructions for CEQ 8000 sequencer and results were analyzed using the accompanying software.

Statistical and Morphometric analysis

All statistical analyses for morphometric parameters were performed with Statistica programme, version 13 (TIBCO Software Inc. 2017). One-way analysis of variance (ANOVA) was conducted to determine the statistical significance of the main effects (flushing form and physiological status) to DBH, H and HC. All statistical analyses were performed at p≤0.05 level of significance.

SSR variation and genetic diversity within subpopulations

The automated binning process of all molecular marker combinations was carried out using the CEQ 8000 Fragment Analysis Software (Beckman Coulter Co.). The products of the binning process/scoring were presented as a table with an allele list in the same program and then exported into Microsoft Windows Excel 2007 program, being ready for subsequent statistical data analysis.

The number of alleles (Na), effective numbers of alleles (Ne), observed heterozygosity (Ho), expected heterozygosity (He), Shannon’s Information Index (I) were calculated using GenAlex 6.5.1b2 (Smouse et al. 2017). The inbreeding coefficient (FIS) was calculated for each locus with POPGENE (Yeh et al. 1999). Presence of null alleles, long allele dropout, and scoring errors were tested with Micro-Checker (Van Oosterhout et al. 2004). When null alleles were found, FreeNA (Chapuis and Estoup 2007) was used to infer the extent of bias imputed by their presence on FST values. Global FST values and frequencies of null alleles were computed by 10.000 iterations with FreeNa software after ENA correction (Chapuis and Estoup 2007). Furthermore, we conducted a Hardy-Weinberg test for heterozygote deficiency by Monte Carlo randomization (10.000) with ML-NULL software (Kalinowski 2005). If the p-value is small, there is a deficiency of heterozygotes, which may be because there are one or more null alleles at the locus. Allele richness (AR) and private allele richness (PAR; alleles which are unique to a particular population), were calculated using the hierarchical rarefaction method available in HP-RARE (Kalinowski and Taper 2006). 

Genetic divergence between subpopulations

Analysis of molecular variance (AMOVA), which assesses the hierarchical distribution of genetic variation among and within subpopulations for the twelve markers, was assessed with GenAlEx 6.5.1b2 software (Smouse et al. 2017). The tests were implemented using estimates of ΦST based on distances calculated from allelic data. Tests of significance were performed using 9999 permutations within the total dataset. Estimation of number of migrants (Nm) based on the private allele method was also implemented in GENEPOP (Raymond 1995). A model-based Bayesian clustering method was applied to the data to infer genetic structure and define the number of clusters (gene pools) in the dataset using the software STRUCTURE, version 2.3.2 (Evanno et al. 2005). Data was explored using the admixture model with allele frequencies correlated. Parameters were set as follows: 10 replicates of each simulation from K = 2 to K = 10 (K = predefined number of modelled clusters) with a burn in of 100,000 steps followed by 100,000 Markov chain—Monte Carlo iterations. Furthermore, K was computed with the Structure Harvester program (Earl 2012) by using the Evanno method (Evanno et al. 2005). In order to reveal differentiation patterns, we calculated population pairwise FST (Weir and Cockerham 1984) by using the Arlequin 3.1 software (Excoffier et al. 2005), applying 10.000 permutations with significance level 0.05.

Genetic divergence between early and late flushing individuals

To address whether there are differences in genetic diversity between early and late flushing individuals, we performed separate analysis within each category (Supplementary material 1). Individuals within each subpopulation (QRBMO & STEMO) were labelled according to their early (QRBE & STEMOE) or late flushing (QRBL & STEMOL). Specifically, number of alleles (Na), effective numbers of alleles (Ne), observed heterozygosity (Ho), expected heterozygosity (He), Shannon’s Information Index (I) were calculated using GenAlex 6.5.1b2 (Smouse et al. 2017). In order to test for correlation between genetic and geographic distance for each subpopulation of early and late flushing individuals in a population, a Mantel test based on 9999 permutations was performed using the matrices of pairwise genetic (GD) and geographical distances (GGD) in GenAlEx (Peakall and Smouse 2006). Moreover, using the same program, a Principal coordinate analysis (PCoA) based on genetic distances from SSR allele markers was performed for each subpopulation. Lastly, in order to study QrZAG87 which is the locus responsible for bud burst we calculated number of alleles (Na), effective numbers of alleles (Ne), observed heterozygosity (Ho), expected heterozygosity (He), Shannon’s Information Index (I) and number of private alleles for each subpopulation. To test if there is significant genetic differentiation between the four subpopulations for early and late flushing individuals we performed FST calculation for each pair of subpopulation with FreeNA software using ENA method (Chapuis and Estoup 2007).

Results

Morphometrics

Morphometric results are presented numerically in Table 2. Statistically significant differences in terms of H were observed between oak trees characterized by different flushing forms in the subpopulation 1, QRBMO (Figure 2a). For subpopulation 2, STEMO, significant differences were demonstrated between measured oak trees concerning HC (Figure 2b). Considering results presented as a single population, significant differences were recorded between oak tress concerning HC (Figure 2c), similarly to subpopulation 2.

Trudic_Fig2.jpg

Figure 2 (a-c): Statistically significant differences (p<0.05) were observed between Q. robur trees for a) different flushing forms in terms of absolute height of trees (H) in subpopulation 1, b) different physiological status in terms of height to the crown (HC) in subpopulation 2, and c) different physiological status in terms of height to the crown (HC) (after merging both subpopulations). / Abbildung 2 (a-c): Statistisch signifikante Unterschiede (p <0,05) zwischen Q. robur zeigten sich für a) unterschiedliche Austriebszeitpunkte hinsichtlich Baumhöhe (H) in Subpopulation 1, b) unterschiedlicher physiologischer Status hinsichtlich Kronenhöhe (HC) in Subpopulation 2 und c) unterschiedlicher physiologischer Status hinsichtlich Kronenhöhe (HC) (beide Subpopulationen zusammengeführt).

Trudic_Tab2.jpg

Table 2: Variability of diameter at breast height (cm), absolute height of trees (m) and height to the crown (m) between genotypes with different flushing form and physiological status. Underlined numbers indicate significant differences between groups determined by ANOVA at p≤0.05 level of significance. DBH – Diameter at breast height (cm); H – absolute height of trees (m); H – height to the crown (m). / Tabelle 2: Variabilität des Brusthöhendurchmessers (cm), der absoluten Baumhöhe (m) und der Kronenhöhe (m) zwischen Genotypen mit unterschiedlichem Austriebszeitpunkt und physiologischem Status. Unterstrichene Zahlen zeigen signifikante Unterschiede zwischen Gruppen an, die durch ANOVA bei einem Signifikanzniveau von p≤0,05 bestimmt wurden. DBH - Durchmesser in Brusthöhe (cm); H - absolute Höhe der Bäume (m); H - Höhe zur Krone (m).

Trudic_Tab3.jpg

Table 3: Genetic variation parameters: Na (total number of alleles), Ne (mean number of effective alleles), I (Shannon’s information index), observed heterozygosity (Ho), He (expected heterozygosity, Nei 1973), Fis (inbreeding coefficient), Global FST and Fnull (null allele frequency), are presented. Asterisk (*) indicates statistically significance and letter “a” indicates a small p-value which may be because there is one or more null alleles at the locus. / Tabelle 3: Parameter der genetischen Variation: Na (Gesamtzahl der Allele), Ne (mittlere Anzahl der effektiven Allele), I (Shannons Informationsindex), beobachtete Heterozygotie (Ho), He (erwartete Heterozygotie, Nei 1973), Fis (Inzuchtkoeffizient), Global FST und Fnull (Null-Allel-Frequenz); Das Sternchen (*) zeigt die statistische Signifikanz und Buchstabe “a” einen kleinen p-Wert an, was möglicherweise durch Null-Allele beeinflusst wird.

SSR variation and genetic diversity within subpopulations

A total number of 255 alleles were observed across the twelve SSR studied in the 83 individuals. The number of alleles (Na) per locus ranged from nine (MSQ13) to 35 (QrZAG65). The effective number of alleles varied from 1.615 (QpZAG96) to 22.774 (QpZAG65) with an average of 9.060 per locus. Average genetic diversity (He) within subpopulations for all SSR loci was 0.794. Shannon Information Index ranged from 0.930 (QrZAG96) to 3.296 (QrZAG65) with an average of 2.270 (Table 3). From Micro-checker ver.2.2.0.3 analysis, no evidence of large allele dropout, scoring of stutter peaks, and not amplified values was found for any locus. ML Null software showed a small p-value (0.000) for QpZAG104 indicating a deficiency of heterozygotes which may be because there are one or more null alleles at the locus. We further estimated null allele frequency with the FreeNa software (EM algorithm was used) and values ranged from 0.000 (for six loci QpZAG15, QpZAG110, QrZAG7, QrZAG87, QrZAG112 and QrZAG11) to 0.075 (QpZAG104). With respect to those results we retained all loci for further analysis due to the fact that according to a simulation study from Chapuis and Estoup (2007), ignoring the presence of null alleles within the observed range of frequency (5–8% on average across loci) will only slightly bias classical estimates of population differentiation (such as FST). Global FST was also computed with FreeNa software by using ENA correction and 10.000 bootstrap re-sampling over loci and values ranged from 0.000 (QpZAG110, QrZAG20 & QrZAG96) to 0.023 (QpZAG15). Fis ranged from -0.080 (QpZAG15) to 0.156 (QpZAG104) with an average of 0.011 across all loci. According to H-W tests from GENEPOP all loci presented a p-value <0.05 indicating no departure from H-W equilibrium (data not shown). 

Genetic variation within and between subpopulations

The effective number of alleles (Ne) was 7.893 and 9.192 for the two subpopulations (STEMO (subpopulation 1, 41 samples) and QRBMO (subpopulation 2, 42 samples, Table 4), respectively. The expected heterozygosity (He) was 0.776 and 0.799 while private allelic richness (PAR) was 3.730 and 5.120 respectively for the QRBMO and STEMO. AMOVA showed that the main proportion of the genetic variation was found within subpopulations (99.3%), rather than among subpopulations (0.7%). The number of detected migrants between all subpopulations after correction for size was eight (Nm = 7.55). The first two axes (principal coordinates) of the Principal Coordinate Analysis (Figure 3) explained 22.48% of the total variation. For population comparison FST was 0.003 with a p-value = 0.041 indicating significant differentiation between the two subpopulations.

Trudic_Tab4.jpg

Table 4: Genetic diversity parameters at the subpopulation level based on twelve loci: Ne (number of effective alleles), Ho (observed heterozygosity), He (expected heterozygosity), I (Shannon’s information index), AR (allelic richness), PAR (private allelic richness), FIS (inbreeding coefficient), FST (gene differentiation coefficient). / Tabelle 4: Genetische Diversitätsparameter auf Subpopulationsebene basierend auf zwölf SSR-Genorten: Ne (mittlere Anzahl der effektiven Allele), Ho (beobachtete Heterozygotie), He (erwartete Heterozygotie), I (Shannons Informationsindex), AR (allelische Vielfalt), PAR (private allelische Vielfalt)), FIS (Inzuchtkoeffizient), FST (Gendifferenzierungskoeffizient).

Trudic_Fig3.jpg

Figure 3: Principal coordinate analysis (PCoA) based on genetic distances from SSR allele markers explaining 22.48% of variation. The PCoA also explained 17.71 and 20.23% for QRBMO and STEMO respectively. / Abbildung 3: Hauptkoordinatenanalyse basierend auf der genetischen Distanz der SSR_Marker, die insgesamt 22.48% der Variation erklären. Die PCoA erklärte auch 17.71 und 20.23% für QRBMO bzw. STEMO.

Trudic_Fig4.jpg

Figure 4: Delta K value suggesting six genetic clusters using the method of Evanno et al. (2005) method. / Abbildung 4: Delta K Wert weist hin auf sechs genetischen Clustern basierend auf Methode nach Evanno et al. (2005).

Genetic divergence within early and late flushing individuals for each subpopulation

Number of alleles (Na), number of effective alleles (Ne), Shannon information Index (I), observed (Ho) and expected heterozygosity (He), Principal Coordinate analysis (PCA) for early and late flushing individuals are being presented for each of the two subpopulations in supplementary material. Moreover, analysis: a) for the locus QRZAG87, b) Cavalli-Sforza and Edwards (1967) genetic distance for each pair of populations using INA correction described in Chapuis and Estoup (2007) and c) pairwise Fst using ENA correction among four Quercus robur populations in Serbia, are also presented in the same file. The Mantel test results (data is provided in the Supplementary material) revealed correlation between genetic and geographic distances for early and late flushing individuals for STEMO population (R2 = 0.0256, p = 0.03) whereas no correlation was found for QRBMO early and late flushing individuals (R2 = 0.032, p = 0.15).

Quercus robur population structure and differentiation

A model-based clustering algorithm was used to identify subgroups with distinctive allele frequencies without prior information on population structure. Evanno et al. (2005) method identified K = 6 as the most likely number of clusters (Figure 4 and Figure 5).

Trudic_Fig5.jpg

Figure 5: STRUCTURE bar graphs of the 83 individuals at K = 6. Distribution of cluster memberships at population (upper figure) and individual levels in the two Q. robur populations, obtained using the STRUCTURE method (Pritchard et al. 2000). Number of clusters, K = 6. Each vertical bar in the histograms (lower figure) represents the proportion of cluster memberships in each of the 83 Q. robur individuals. / Abbildung 5: STRUCTURE Balken-Diagramm der 83 Individuen bei K = 6. Verteilung der Clustermitgliedschaften auf Population (obere Abbildung) und auf individueller Ebene in den beiden Q. robur-Populationen, erhalten nach der STRUCTURE-Methode (Pritchard et al. 2000). Anzahl der Cluster, K = 6. Jeder vertikale Balken in den Histogrammen (untere Abbildung) repräsentiert den Anteil der Cluster-Mitgliedschaften in jedem der 83 Q. robur Individuen.

Discussion

Trudic_Tab5a.jpg

Table 5a and 5b: Comparison of genetic diversity on SSR loci between our study, Katičić et al. 2018 (clonal seed orchard) and studies on natural Q. robur populations. N/P: number of trees/population, Ho: observed heterozygosity, He: expected heterozygosity, n.a: information not available (*: data obtained from Neophytou (2014)). / Tabelle 5a und 5b: Vergleich der genetischen Vielfalt an SSR-Genorten zwischen unserer Studie, der Studie von Katičić et al. 2018 (Samenplantage) und Studien, in natürlichen Q. robur-Populationen. N/P: Anzahl der Bäume/Population, Ho: beobachtete Heterozygotie, He: erwartete Heterozygotie, n.a.: Informationen nicht verfügbar (*: Daten wurden von Neophytou (2014) erhalten).

Oak populations present a high degree of pollen immigration rates (pollen may travel up to several tens of kilometers; Buschbom et al. 2011), and inbreeding levels are low (Neophytou et al. 2015, Vranckx et al. 2014); which explains the low genetic differentiation among populations found in several studies. In the current study, we also observed a similar pattern of low genetic differentiation between studied subpopulations, as FST was also low (FST = 0.003). This was also confirmed by the results of ANOVA, which showed that the majority of the variation was found within subpopulations (99%), rather than among subpopulations (1%).

The number of alleles per locus (Na) ranged from 9 to 35, with a mean value of 21.25 and was almost equal to Vranckx et al. (2014) (Na = 21.26) and higher from Bogdan et al. (2017) (Na = 20.25) where clonal seeds orchards were studied. Moreover, AR was significantly higher for both subpopulations (QRBMO: 14.4 & STEMO: 15.8) in comparison to Vranckx et al. (2014) study (adult: 11.5; seedling cohort: 10.7), indicating a rich genetic germplasm for the species. Bogdan et al. (2017) found lower value of AR for the three clonal orchards (12.67, 12.58 and 13.06, respectively) in Croatia, as well.

Trudic_Tab5b.jpg

Table 5a and 5b: Comparison of genetic diversity on SSR loci between our study, Katičić et al. 2018 (clonal seed orchard) and studies on natural Q. robur populations. N/P: number of trees/population, Ho: observed heterozygosity, He: expected heterozygosity, n.a: information not available (*: data obtained from Neophytou (2014)). / Tabelle 5a und 5b: Vergleich der genetischen Vielfalt an SSR-Genorten zwischen unserer Studie, der Studie von Katičić et al. 2018 (Samenplantage) und Studien, in natürlichen Q. robur-Populationen. N/P: Anzahl der Bäume/Population, Ho: beobachtete Heterozygotie, He: erwartete Heterozygotie, n.a.: Informationen nicht verfügbar (*: Daten wurden von Neophytou (2014) erhalten).

Results from the twelve SSR markers, showed that mean Ho and He values were smaller herein (0.779 and 0.794, respectively) in comparison to Vranckx et al. (2014) who studied four forest stands in Belgium, and Neophytou et al. (2010) who studied three populations from Greece, Bulgaria and Germany. This was expected, since we analyzed one population with two close neighboring subpopulations, where previously mentioned authors, have larger inter populational sample to compare. In more detail, the Tables 5a and 5b presents a comparison of Ho and He with other studies conducted in Q. robur genotypes. More specifically, marker QrZAG110 had lower values in comparison to other studies and together with QrZAG112 had Ho and He values almost equal, which was also the case in the study of Katičić et al. (2018). Interestingly, locus QrZAG96 displayed the smallest value of He, which is the locus linked to a QTL associated with leaf morphology and differentiate Q. robur and Q. petraea according to Scotti-Saintagne et al. (2004). Similar findings for this locus was found from Muir and Schloetterer (2005), Neophytou et al. (2010), Craciunesc et al. (2011) and Katičić et al. (2018) for Q. robur genotypes. Moreover, locus QrZAG87 (Scotti-Saintagne et al. 2004) which is situated in linkage group LG2 of Quercus robur map, in one of the genome regions coding for bud burst had the higher value compared to all other studies concerning Ho and He. This is probably the result of choosing in our sampling scheme, late (var. tardiflora) and early (var. praecox) flushing individuals. The same locus was nominated as one of the five loci which have higher intraspecific than interspecific provenance discriminating power according to Neophytou et al. (2010). Lastly, analysis for the locus QrZAG87 exhibited similar values for early and late flushing individuals. Further investigation is required in order to gain better understanding about the genetic component underlying bud burst in Q. robur.

Results from structure analysis set K = 6 according to Evanno et al. (2005) indicating no significant clustering between the subpopulations. Taking into account the estimated number of migrants which was eight and the fact that the two subpopulations are near in distance (approx. 2 km distance from each other) this result was expected. Buschbom et al. (2011) also pointed out in a Q. robur isolated stand that effective long-distance gene flow highly contributes to genetic diversity and the impact of long-distance pollen immigration from very diverse pollen sources on genetic diversity is an important factor associated with the survival and the diversification of stands.

Moreover, results from comparison of late and early flushing individuals (supplementary material document) did not show to follow any clear pattern. Mean observed and expected heterozygosity was quite similar between early and late flushing individuals for both subpopulations (QRMBO and STEMO). Mantel test was not significant for QRBMO whereas was for STEMO subpopulation. An important finding was that according to genetic distances and Fst values the two subpopulations differ significantly. 

According to EUFORGEN’s technical guidelines on conservation and use of pedunculate (and sessile) oaks, in situ conservation methods should be generally preferred. If natural regeneration methods are not sufficient, an adapted and specified ex situ conservation programme including a controlled autochthonous reproductive material system (e.g. clonal seed orchards) should be used as well to preserve the above subpopulations (Ducousso and Bordacs 2004).

Conclusions and recommendations for future conservation and management measures

According to the results found herein, individuals of pedunculate oak from Vojvodina province of Serbia showed high He values and a high number of migrants in comparison to previously published studies in Europe. These two stands of Q. robur have significant gene richness that should be conserved firstly in situ and afterwards, in ex situ scheme.

In their latest publication, Stojnić et al. (2014) gave also specific recommendations for future pedunculate oak breeding and conservation programmes in AP Vojvodina region of Serbia. Some of those recommendations include the setup of experiments that will focus on genetic, phenological, chemical and morphometric analysis of pedunculate oak in natural conditions, greenhouses and in vitro culture. Stojnić et al. (2014) also recommended creation of mapping oak populations in Serbia in order to monitor the development and inheritance of quantitative traits that are regulated by multiple genes. In conclusion, a strategic plan for conservation of Q. robur populations (in situ & ex situ) in Serbia will ensure the survival of the species under the ongoing climatic changes.

Author contributions

B.T., E.A. and B.F. conceived the ideas, designed the study and coordinated the research activities. B.T., A.P. and S.S. collected the leaves samples in the field and performed morphometric measurement and observational phenotypic candidacy of all individual trees. B.T. and B.F. performed laboratory analyses. E.A. and S.S. conducted entire statistical analysis; C.N supported this work with his valuable inputs and comments. B.T. and E.A. wrote the first version of the manuscript. C.N. was involved in shaping the manuscript and facilitating its development until its drafting version. Regarding to her immense contribution to the development of study and paper, B.T. and E.A. are declaring shared first authorship. All other authors contributed critically to the drafts and gave final approval for publication.

Acknowledgments

This study was conducted under the project “Spruce and Oak Genetic Structure” (SOGeneS), supported by FP7 project Trees 4 Future—Transnational Access.

Conflicts of Interest

The authors declare no conflicts of interest.

References

Alberto, F., Bouffier, L., Louvet, J.M., Lamy, J.B., Delzon, S. & Kremer, A. (2011) Adaptive responses for seed and leaf phenology in natural populations of sessile oak along an altitudinal gradient. Journal of Evolutionary Biology, 24, 1442–1454.

Banković, S., Medarević, M., Pantić, D., Petrović, N. (2009) Nacionalna inventura šuma Republike Srbije – šumski fond Republike Srbije. Ministarstvo poljoprivrede, šumarstva i vodoprivrede, Uprava za šume, Beograd. str. 244. 

(PDF) Conservation of pedunculate oak (Quercus robur L.) genetic resources at the territory of Public Enterprise "Vojvodinašume". Available from: https://www.researchgate.net/publication/301661824_Conservation_of_pedunculate_oak_Quercus_robur_L_genetic_resources_at_the_territory_of_Public_Enterprise_Vojvodinasume [accessed Sep 07 2020].

Belkhir, K., Borsa, P., Chikhi, L., Raufaste, N., Bonhomme, F. (2004) GENETIX 4.05, Windows TM software for population genetics Laboratoire génome, populations, interactions, CNRS UMR 5000

Bogdan, S., Ivanković, M., Temunović, M., Morić, M., Franjić, J., & Bogdan, I. K. (2017) Adaptive genetic variability and differentiation of Croatian and Austrian Quercus robur L. populations at a drought prone field trial. Annals of Forest Research, 60(1), 33-46.

Buschbom, J., Yanbaev, Y., Degen, B. (2011) Efficient long-distance gene flow into an isolated relict oak stand. Journal of Heredity, 102(4), 464-472.

Cavalli-Sforza, L. L. & Edwards, A. W. (1967) Phylogenetic analysis. Models and estimation procedures. American journal of human genetics, 19(3 Pt 1), 233.

Chapuis, M-P. & Estoup, A. (2007) Microsatellite null alleles and estimation of population differentiation. Molecular biology and evolution, 24:621-631.

Conover, D.O., Duffy, T.A. & Hice, L.A. (2009) The covariance between genetic and environmental influences across ecological gradients: reassessing the evolutionary significance of countergradient and cogradient variation. Annals of the New York Academy of Sciences, 1168, 100–129.

Convention on Biological Diversity, 1992, UN (https://www.cbd.int/doc/legal/cbd-en.pdf, access date: 25.11.2020.)

Craciunesc, I., Ciocîrlan, E., Sofletea, N., Curtu, A. (2011) Genetic diversity of pedunculate oak (Quercus robur L.) in Prejmer natural reserve Bulletin of the Transilvania University of Brasov Forestry, Wood Industry, Agricultural Food Engineering Series, II 4:15.

Dantec, C. F., Ducasse, H., Capdevielle, X., Fabreguettes, O., Delzon, S., Desprez‐Loustau, M. L. (2015) Escape of spring frost and disease through phenological variations in oak populations along elevation gradients. Journal of Ecology, 103(4), 1044-1056.

Deguilloux, MF., Dumolin‐Lapègue, S., Gielly, L., Grivet, D., Petit, R. (2003) A set of primers for the amplification of chloroplast microsatellites in Quercus. Molecular Ecology Notes, 3:24-27.

Dow, BD., Ashley, MV., Howe, HR. (1995) Characterization of highly variable (GA/CT)n microsatellites in the bur oak, Quercus macrocarpa. Theor Appl Gen, 91:137-141.

Dubravac, T., Dekanić, S., Roth, V. (2011) Dinamika oštećenosti i struktura krošanja stabala hrasta lužnjaka u šumskim zajednicama na gredi i u nizi – rezultati motrenja na trajnim pokusnim plohama. Šumarski list – posebni broj: 74-89.  

(PDF) Conservation of pedunculate oak (Quercus robur L.) genetic resources at the territory of Public Enterprise "Vojvodinašume". Available from: https://www.researchgate.net/publication/301661824_Conservation_of_pedunculate_oak_Quercus_robur_L_genetic_resources_at_the_territory_of_Public_Enterprise_Vojvodinasume [accessed Sep 07 2020].

Ducousso, A. and S. Bordacs. (2004) EUFORGEN Technical Guidelines for genetic conservation and use for pedunculate and sessile oaks (Quercus robur and Q. petraea). International Plant Genetic Resources Institute, Rome, Italy. 6 pages.

Dumolin, S., Demesure, B., & Petit, R. J. (1995) Inheritance of chloroplast and mitochondrial genomes in pedunculate oak investigated with an efficient PCR method. Theoretical and applied genetics, 91(8), 1253-1256.

Earl, DA. (2012) STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conservation genetics resources, 4:359-361.

Evanno, G., Regnaut, S., Goudet, J. (2005) Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Molecular ecology, 14:2611-2620.

Excoffier, L., Laval,. G, Schneider, S. (2005) Arlequin (version 3.0): an integrated software package for population genetics data analysis. Evolutionary bioinformatics, 1:117693430500100003.

Fussi, B., Westergren, M., Aravanopoulos, F., Baier, R., Kavaliauskas, D., Finzgar, D., Alizoti, P., Bozic, G., Avramidou, E., Konnert, M. and Kraigher, H. (2016) Forest genetic monitoring: an overview of concepts and definitions. Environmental Monitoring and Assessment, 188(8), p.493. 

Guichoux, E., Lagache, L., Wagner, S., Léger, P., Petit, R. (2011) Two highly validated multiplexes (12‐plex and 8‐plex) for species delimitation and parentage analysis in oaks (Quercus spp.). Molecular Ecology Resources, 11:578-585.

Kalinowski, ST. (2005) Hp‐rare 1.0: a computer program for performing rarefaction on measures of allelic richness. Molecular ecology notes, 5:187-189.

Kalinowski, ST., Taper, ML. (2006) Maximum likelihood estimation of the frequency of null alleles at microsatellite loci. Conservation Genetics, 7:991-995.

Kampfer, S., Lexer, C., Glossl, J., & Steinkellner, H. (1998) Brief report characterization of (GA) n microsatellite loci from Quercus robur. Hereditas, 129(183), 1-86.

Katičić Bogdan, I., Kajba, D., Šatović, Z., Schüler, S., Bogdan, S. (2018) Genetic diversity of Pedunculate oak (Quercus robur L.) in clonal seed orchards in Croatia, assessed by nuclear and chloroplast microsatellites. South-east European forestry, 9:29-46.

Kavaliauskas, D., Fussi, B., Westergren, M., Aravanopoulos, F., Finzgar, D., Baier, R., Alizoti, P., Bozic, G., Avramidou, E., Konnert, M. and Kraigher, H. (2018) The interplay between forest management practices, genetic monitoring, and other long-term monitoring systems. Forests, 9(3), p.133. 

Kostić, S., Levanič, T., Orlović, S., Matović, B., & Stojanović, D. B. (2019) Pendunctulate and Turkey Oaks Radial Increment and Stable Carbon Isotope Response to Climate Conditions through Time. Topola, (204), 29-35.

Kremer, A., Potts, B.M. & Delzon, S. (2014) Genetic divergence in forest trees: understanding the consequences of climate change. Functional Ecology, 28, 22–36.

Lepais, O., Léger, V., & Gerber, S. (2006) Short note: high throughput microsatellite genotyping in oak species. Silvae Genetica, 55(1-6), 238-240.

Metzger, M.J., Bunce, R.G.H., Jongman, R.H.G., Mücher, C.A. and Watkins, J.W. (2005) A climatic stratification of the environment of Europe. Global Ecology and Biogeography, 14: 549-563. https://doi.org/10.1111/j.1466-822X.2005.00190.x

Morić, M. (2016) Genetic diversity of pedunculated oak (Quercus robur L.) in field trials whit progeny from selected seeds stands (Doctoral dissertation, Šumarski fakultet, Sveučilište u Zagrebu).

Muir, G., Schloetterer, C. (2005) Evidence for shared ancestral polymorphism rather than recurrent gene flow at microsatellite loci differentiating two hybridizing oaks (Quercus spp.). Molecular ecology, 14:549-561.

Neophytou, C., Aravanopoulos, FA., Fink, S., Dounavi, A. (2010) Detecting interspecific and geographic differentiation patterns in two interfertile oak species (Quercus petraea (Matt.) Liebl. and Q. robur L.) using small sets of microsatellite markers. Forest Ecology and Management, 259:2026-2035.

Neophytou, C., Gärtner, SM., Vargas-Gaete, R., Michiels, H-G. (2015) Genetic variation of Central European oaks: shaped by evolutionary factors and human intervention? Tree Genetics & Genomes, 11:79.

Neophytou, Charalambos (2014) Bayesian clustering analyses for genetic assignment and study of hybridization in oaks: effects of asymmetric phylogenies and asymmetric sampling schemes, Dryad, Dataset, https://doi.org/10.5061/dryad.b64b4 and GeneAlex was used to estimate parameters Ho & He for Q. robur.

Orlović, S., D. Šimunovački, Z. Đorđević, A. Pilipović and N. Radosavljević (2008) The preservation of the gene pool and the production of seeds of oak (Quercus robur L.), Monograph: 250 Years of Forestry Ravni Srem, JP Vojvodinašume (Sr).

Peakall, R., Smouse, PE. (2006) GENALEX 6: genetic analysis in Excel. Population genetic software for teaching and research. Molecular ecology notes, 6:288-295.

Pekeč, S., Orlović, S., Katanić, M., Stojnić, S., Drekić, M. (2017) Fenološka osmatranja hrasta kitnjaka (Quercus petrea Matt./Liebl.) i hrasta lužnjaka (Quercus robur L.) na području Vojvodine. Topola, 199/200:11-12.

Perić, S., Gračan, J., & Dalbelo-Bašić, B. (2000) Flushing variability of pedunculate oak (Quercus robur L.) in the provenance experiment in Croatia. Glasnik za šumske pokuse, 37, 395-412.

Pilipović, A., Drekić M., Stojnić S., Nikolić N., Trudić B., Milović M., Poljaković-Pajnik L., Borišev M., Orlović S. (2020). Physiological Responses of Two Pedunculate Oak (Quercus robur L.) Families to Combined Stress Conditions - Drought and Herbivore Attack. Šumarski list (In press). 

Raymond, M. (1995) Population genetics software for exact test and ecumenicism. J Heredity, 86:248-249.

Saintagne, C., Bodénès, C., Barreneche, T., Pot, D., Plomion, C., & Kremer, A. (2004) Distribution of genomic regions differentiating oak species assessed by QTL detection. Heredity, 92(1), 20-30.

Scotti-Saintagne, C., Mariette, S., Porth, I., Goicoechea, P.G., Barreneche, T., Bodénes, C., Burg, K. and Kremer, A. (2004) Genome scanning for interspecific differentiation between two closely related oak species [Quercus robur L. and Q. petraea (Matt.) Liebl.]. Genetics, 168(3), pp.1615-1626.

Smouse, PE., Bank, SC,, Peakall, R. (2017) Converting quadratic entropy to diversity: Both animals and alleles are diverse, but some are more diverse than others. PloS one, 12:e0185499.

Steinkellner, H., Lexer, C., Turetschek, E., & Glössl, J. (1997) Conservation of (GA)n microsatellite loci between Quercus species. Molecular Ecology, 6(12), 1189-1194.

Stojanović, D. B., Levanič, T., Matović, B., & Orlović, S. (2015) Growth decrease and mortality of oak floodplain forests as a response to change of water regime and climate. European Journal of Forest Research, 134(3), 555-567.

Stojnić, S., Kovačević, B., Kebert, M., Vaštag, E., Bojović, M., Neđić, M. S., & Orlovic, S. (2019) The use of physiological, biochemical and morpho-anatomical traits in tree breeding for improved water-use efficiency of Quercus robur L. Forest systems, 28(3), 5.

Stojnić, S., Trudić, B., Galović, V., Šimunovački, Đ., Đorđević, B., Rađević, V., & Orlović, S. (2014) Conservation of pedunculate oak (Quercus robur L.): Genetic resources at the territory of public enterprise 'Vojvodinašume'. Topola, (193-194), 47-71.

Streiff, R., Labbé, T., Bacilieri, R., Steinkellner, H., Glössl, J., Kremer, A. (1998) Within‐population genetic structure in Quercus robur L. and Quercus petraea (Matt.) Liebl. assessed with isozymes and microsatellites. Molecular Ecology, 7:317-328.

Thomas, F. M. (2008) Recent advances in cause-effect research on oak decline in Europe. CAB Reviews: Perspectives in Agriculture, Veterinary Science, Nutrition and Natural Resources, 3(37), 1-12.

Trudić, B., Galović, V., Orlović, S., Pap, P., & Pekeč, S. (2013) A strategy for the identification of a canditate gene for drought induced stress in pedunculate oak (Quercus robur L.(Q. pedunculata Ehrh.)), Fagaceae. Bulgarian journal of agricultural science, 19(2), 338-346.

Van Oosterhout, C., Hutchinson, WF., Wills, DP., Shipley, P. (2004) MICRO‐CHECKER: software for identifying and correcting genotyping errors in microsatellite data. Molecular Ecology Notes, 4:535-538.

Vastag, E., Cocozza, C., Orlović, S., Kesić, L., Kresoja, M., Stojnić, S. (2020) Half-sib lines of pedunculate oak (Quercus robur L.) respond differently to drought through biometrical, anatomical and physiological traits. Forests, 11, 153. doi: https://doi.org/10.3390/f11020153

Vitasse, Y., Bresson, C.C., Kremer, A., Michalet, R. & Delzon, S. (2010) Quantifying phenological plasticity to temperature in two temperate tree species. Functional Ecology, 24, 1211–1218.

Vitasse, Y., Delzon, S., Bresson, C.C., Michalet, R. & Kremer, A. (2009c) Altitudinal differentiation in growth and phenology among populations of temperate-zone tree species growing in a common garden. Canadian Journal of Forest Research, 39, 1259–1269.

Vitasse, Y., Delzon, S., Dufrêne, E., Pontailler, J. Y., Louvet, J. M., Kremer, A., & Michalet, R. (2009a) Leaf phenology sensitivity to temperature in European trees: Do within-species populations exhibit similar responses?. Agricultural and forest meteorology, 149(5), 735-744.

Vitasse, Y., Porté, A. J., Kremer, A., Michalet, R., & Delzon, S. (2009b) Responses of canopy duration to temperature changes in four temperate tree species: relative contributions of spring and autumn leaf phenology. Oecologia, 161(1), 187-198.

Vranckx, G., Jacquemyn, H., Mergeay, J., Cox, K., Kint, V., Muys, B., Honnay, O. (2014) Transmission of genetic variation from the adult generation to naturally established seedling cohorts in small forest stands of pedunculate oak (Quercus robur L.). Forest Ecology and Management, 312:19-27.

Weir, BS., Cockerham, CC. (1984) Estimating F-statistics for the analysis of population structure evolution: 1358-1370.

Yeh, FC., Yang, R., Boyle, T., Ye, Z., Mao, JX. (1999) POPGENE, version 1.32: the user friendly software for population genetic analysis Molecular Biology and Biotechnology Centre, University of Alberta, Edmonton, AB, Canada.