Genetic Diversity and Genetic structure of three sympatric Oak species in serbian Landscape of Outstanding Features "Kosmaj" Assessed by Nuclear Microsatellites

.)

Genetic diversity is a key requirement for the long-term survival of species on an evolutionary time scale (Pakkad et al. 2008, Fady et al. 2022) and is usually rather high in longlived forest tree species (Kremer et al. 2012, Rajora andZinck 2021).Forests play important roles in ecosystem functioning, species diversity maintenance, climate regulation, and soil and water conservation.Thus, the biodiversity and stability of the entire forest ecosystems are mainly dependent on the genetic diversity of forest trees (Sandurska et al. 2019), commonly shaped by pollen and seed dispersal (Dow and Ashley 1998), demographic changes and evolutionary processes, such as genetic drift and migration, which may account for the occurrence of genetic structure in forest tree populations (Fernández-M and Sork 2006).In forest trees, gene flow is particularly important, as it has the ability to both permit local genetic differentiation and introduce variability in isolated populations (Ducousso et al. 1993) as well as to erase present genetic structure over generations.
2 SEEFOR 14(2): early view Genetic diversity loss became particularly important in the context of climate change, relevant also for the production of well-adapted reproductive material used for reforestation (Ivetić et al. 2016, Katičić Bogdan et al. 2018, Popović et al. 2022).Information on genetic variability, genetic structure, and inbreeding is therefore required for the effective management of forest populations (Craft and Ashley 2007).
Widely known species richness of the Balkan Peninsula is a result of complex geological history, and interactions between populations, species, and ecosystems (Savić 2008).Among many forest tree species, oaks are one of the ecologically and economically most important genera in the northern hemisphere.A total of 12 oak species occur in the Balkans, including six species belonging to the Quercus (white oak) group (Nixon 1993).In the genus Quercus, several taxonomic groups are characterized by complex patterns of variation, leading to difficulties in the identification of taxa to the species level (Bruschi et al. 2000, Gömöry et al. 2001).
Resolving the species boundaries in oaks can be challenging due to the presence of natural hybrids and introgressive forms, especially in contact zones where two or more species occur (Curtu et al. 2007, Viscosi et al. 2009, Yücedağ and Gailing 2013, Lyu et al. 2018).One of such sites, where three oak species, namely Quercus petraea (Matt.)Liebl., Quercus pubescens Willd., and Quercus frainetto Ten., grow in sympatry, is found in Serbia, in the protected zone (protection degree II) within the Landscape of Outstanding Features "Kosmaj" at the Mt.Kosmaj.This site is geographically close to the metropolitan region and it is threatened by habitat loss and degradation, mostly due to the expansion of agricultural and urban areas, illegal cutting of forests for firewood and other wood products, and excessive grazing.These autochthonous oak populations are important both from scientific and economic point of view (Jovanović et al. 2023).Studies of genetic diversity of oaks in Serbia were conducted by several authors, using both chloroplast and nuclear DNA (e.g., Milovanović 2009, Šijačić-Nikolić et al. 2009, Ballian 2010, Neophytou et al. 2010, Kesić et al. 2021, Šijačić-Nikolić et al. 2021, Popović et al. 2022).However, at the region of the Landscape of Outstanding Features "Kosmaj" no such research has been conducted so far.
Due to the overexploitation of forests, together with poorly understood environmental changes, which have resulted in the decline of many oak species, oaks are involved in an increasing number of studies related to gene flow and genetic structure (e.g.Ducousso et al. 1993, Lepais et al. 2009, Fortini et al. 2015, Sandurska et al. 2019, Leroy et al. 2020).Gene flow levels in oaks are high and mostly related to the life history traits of its species, such as phenology, mating system, wind pollination, acorn production, etc. (Ducousso et al. 1993).Due to high levels of intraspecific diversity and hybridization between the species which cause the boundaries of species differentiation to be less distinctive, discerning some oak species at the molecular level and obtaining species-specific diagnostic markers can be challenging (Kelleher et al. 2005).Because of the highly informative and codominant nature, hypervariability, reproducibility, and the possibility of parallel amplification nSSRs are widely used for genotyping in plant and animal species.Although the usage of nSSRs may be challenging (e.g., Kerkez Janković et al. 2019), they have been proven suitable for assessing the levels of genetic diversity in oak species (Kesić et al. 2021).Furthermore, Muir et al. (2000) have demonstrated that nSSRs may be used for differentiation of phylogenetically distant species like Q. robur and Q. petraea despite the high levels of interspecific gene flow.For phylogenetically close species, genetic diversity in populations in hybrid zones or at the distributional margins of species, EST-SSR markers are useful (Ueno et al. 2008, Aizawa et al. 2018).
The aims of this study were to assess the genetic diversity and structure of Q. petraea, Q. pubescens, and Q. frainetto from the Landscape of Outstanding Features "Kosmaj" (Mt.Kosmaj, Serbia) using nuclear microsatellites.Levels of inter-and intraspecific genetic variability reveled in this study can lay the groundwork for setting the guidelines for the conservation of the available gene pool and prescribing measures for future forest management.

MAtErIALs AND MEtHODs study site and Plant Material
A total of 187 trees found in the protected area of Mt.Kosmaj (44°28'17.68"N,20°34'32.04"E)were selected for genotyping -65 of Q. pubescens, 60 of Q. petraea, and 62 of Q. frainetto (Figure 1).The protection zone occupies 3514.50 ha of a forest complex surrounded by predominantly agricultural land and rural settlements.At the lower altitudes (250-400 m a.s.l.) dominant forest type is Quercetum frainetto-cerris Rudski, above 400 m a.s.l.Querco-Carpinetum serbicum Rudski and from 500 to 626 m a.s.l.Fagetum montanum Rudski.At Mt. Kosmaj oaks are abundantly present in sympatry, which can be seen from forest typology previously mentioned, but are negatively impacted by the lack of natural regeneration and the coppice origin.
Samples for molecular genetics analyses comprised of young leaves which were collected in May 2022.The selection of trees was based on phenotypic characteristics (trunk and crown appearance) and general health status (absence of entomological and phytopathological damage).From each selected tree up to five normally developed, healthy leaves were sampled, herbarized, dried in silica gel, and kept in a freezer prior to the analyses.Thus, only adult trees, distant at least 50 m from each other, were sampled.

DNA Extraction, PCr Amplification and Fragment sizing
Samples were processed in Biotechnology Laboratory at the Faculty of Forestry, University of Belgrade (Serbia).For each selected tree up to 20 mg of dried leaves were homogenized with TissueLyser II (Qiagen, Valencia, CA, USA) and used for extraction of the total genomic DNA with peqGOLD Plant DNA Mini Kit (PEQLAB).Genomic DNA was quantified and assessed for purity utilizing NanoVue (GE Healthcare Europe, Freiburg, Germany).DNA solutions were diluted to working concentrations of 50 ng•µl -1 .
For genotyping, 20 nuclear microsatellites were selected, 14 of which proved to be very informative and reliable after testing on a panel of 8 samples (Table 1).For the parallel amplification Type-it Microsatellite PCR Kit (Qiagen) was used.Microsatellite loci were grouped into two mixes: ОМ1 -PIE239, FIR004, QrZAG90, QrZAG108, MSQ13, GOT004, QrZAG87, QpZAG104, QrZAG11, QrZAG103, QrZAG102, and ОМ2 -QpZAG36, QrZAG101, MAQ4, PIE242, QrZAG20, QpZAG1/2, QpZAG58, QrZAG7, QpZAG110.Forward (F) primers from each primer pair of all selected loci were labelled with one out of the four different fluorescent dyes from the Dye set G5 (DS-33, Applied Biosystems, USA) to enable fragment sizing on an automated DNA fragment analyser.PCRs were carried out using Multigene Opti Max (Labnet International, Inc.).The PCR amplifications were performed as follows: initial denaturation at 95°C for 5 min; then 28 cycles of denaturing for 30 s at 95°C, annealing for 90 s at a 60°С, and extension for 30 s at 72°C; with a final extension at 60°C for 30 min.The PCR products were visualized on 1% agarose gels.PCR amplification products were separated commercially by Center for Forensic and Applied Molecular Genetics at the Faculty of Biology, University of Belgrade via capillary electrophoresis using 96-capillary 3730xl DNA Analyzer automated sequencer (Applied Biosystems, Inc. USA).The lengths of PCR amplification products were assessed using GeneMapper (Applied Biosystems Inc., Foster City, USA), with the GeneScanTM-600LIZTM Size Standard (Applied Biosystems).

Data Analysis
The standard parameters of genetic diversity: number of alleles (A), number of private alleles (PA), average number of alleles (Na), effective number of alleles (Ne), observed heterozygosity (H O ), expected heterozygosity (H E ) and coefficient of inbreeding (F), were determined using the GenAlEx 6.5 software package.The effective population size (Ne) was assessed using the NeEstimator software package (Do et al. 2014) with the rejection of alleles whose frequency was ≤0.02.Allelic richness was estimated using the HP-Rare 1.0 software package (Kalinowski 2005).Genetic structure was determined using two approaches: (1) obtained allele frequencies and genotypes generated for 14 molecular markers were used to determine the parameters of genetic differentiation expressed through the F ST values, and PCoA analysis using the program package GenAlEx 6.5; obtained F ST values were also used to assess levels of gene flow, Nm, according to the formula Nm = [(1 / F ST ) -1] / 4; (2) Bayesian method implemented in the STRUCTURE 2.3.4 software package (Pritchard et al. 2000) was used to determine the optimal number of genetic groups (K) using the ΔK Evanno model (Evanno et al. 2005), and to determine the potential substructure of populations.Monte Carlo Markov Chain (MCMC) simulation had burn-in and run lengths of 700,000 iterations each.Ten independent analyzes for each assumed group K = 1-6 were performed.The model of correlated (dependent) frequencies of alleles was used as the allele frequency model, and the admixture model was used as the model of individual kinship.4 SEEFOR 14(2): early view rEsULts

Genetic Diversity
The summary of genetic diversity parameters obtained upon analyzing variability at 14 microsatellite loci in the three oak species is shown in Table 2, while the summary of genetic diversity parameters per species and in overall sample is shown in Table 3.A total of 314 alleles were detected in overall sample.The smallest number of alleles in Q. petraea was observed at the PIE239 locus (6), in Q. pubescens at QrZAG108 (6), and in Q. frainetto at QrZAG87 (7).The largest number of alleles for all tree species was observed at the QrZAG90 locus (A = 36, 35, and 31; Ае = 24.69, 20.64, 22.58 in Q. petraea, Q. pubescens and Q. frainetto, respectively).A total of 249 alleles were detected in Q. petraea, 238 alleles in Q. pubescens, and 213 alleles in Q. frainetto (Table 3).The average number of alleles per locus in the entire sample was 16.67 (SE = 1.07), the average effective number of alleles per locus 8.35 (SE = 0.87), and the average allelic richness obtained by the rarefaction method for 98 gene copies was 16.34.The number of private alleles was 29 in Q. petraea, 34 in Q. pubescens and 18 in Q. frainetto.The effective population size was 305.5 (190.2;719.2) in Q. pubescens, 484.3 (254.8; 3374.5) in Q. petraea, and ∞ (1859.9;∞) in Q. frainetto.The observed heterozygosity (Ho) ranged from 0.731 (SE = 0.051) (Q.pubescens) to 0.753 (SE = 0.053) (Q.petraea), with an average value of 0.745 (SE = 0.030).The expected heterozygosity (He) ranged from 0.794 (SE = 0.033) (Q.petraea) to 0.834 (SE = 0.027) (Q.pubescens), with an average value of 0.817 (SE = 0.019).A low but statistically significant excess of homozygotes was detected in the populations of all tested oak species.

Genetic Differentiation
To better understand the relationships between the three sympatric oak species, the gene flow (Nm) and genetic differentiation coefficient (F ST ) were estimated for all pairs of species (Table 4).The genetic differentiation between species was low but statistically significant (P ≤ 0.05), and ranged from 0.032 among Q.pubescens and Q. petraea, to 0.047 among Q.petraea and Q. frainetto.Consequently, the highest gene flow (Nm = 7.563) was observed between Q. petraea and Q. pubescens.The lowest gene flow (Nm = 5.069) was observed between Q. petraea and Q. frainetto.
Principal coordinate analysis (PCoA) results, obtained by summarizing the genetic distances between genotypes within each of the populations showed that Q. frainetto separated from Q. pubescens and Q. petraea along the first principal coordinate, while the separation of Q. pubsecens and Q. petraea was observed along the second principal coordinate (Figure 2).The first principal coordinate (Coord.1) explained 63.32% of the variability, and the second principal coordinate (Coord.2) 36.68% of the variability, suggesting high reliability of the obtained results.PCoA results obtained by summarizing the genetic distances between genotypes within each of the individuals of the three species (Figure 3) showed the same grouping pattern as in population analysis (Figure 2) despite the small percentage of the total variability explained by the first two principal coordinates (12.27%,Coord. 1 = 6.84%,Coord. 2 = 5.43%).

Genetic structure
The optimal number of genetic groups in the study sample, obtained by the ΔК Evanno model, was four (Figure 4).However, a higher level of hierarchical structure, with six genetic groups, was observed as well.
Results of the STRUCTURE analysis with four (K = 4) and six (K = 6) genetic groups in Q. petraea, Q. pubescens, and Q. frainetto are shown in Figures 5 and 6.
The STRUCTURE analysis showed that Q. frainetto represents a coherent and distinct genetic group that was not substructured, i.e., in which a large number of individuals had a high proportion of assignment to the unique gene pool.In contrast to Q. frainetto, the populations of Q. petraea and Q. pubescens were substructured, i.e., comprised individuals strongly assigned (qi > 0.80) to distinct gene pools.The population of Q. petraea in the area of Mt.Kosmaj consisted of two separate genetic groups.The substructure has also been observed in the population of Q. pubescens (Figure 6).The results of the STRUCTURE analysis under the assumption of six genetic groups (K = 6) showed that Q. pubescens population comprised individuals strongly assigned (qi > 0.80) to three distinct gene pools.8 SEEFOR 14(2): early view

DIsCUssION
For the conservation of forest genetic resources and sustained forest management, knowledge of the levels and distribution of genetic variation is crucial, especially if it is related to genetic processes taking place during reproduction and stand establishment, i.e. during developmental stages that are under high anthropogenic influence (Liesebach and Zaspel 2004).However, complex patterns of genetic diversity and characteristic phenotypic plasticity of oaks hampered the efforts to manage their gene pools and accurately assess the conservation status of most of the oak taxa (Sullivan et al. 2013).While oaks are capable of high rates of local adaptation, their long generation times and immobility make them vulnerable to rapid environmental changes (Yan et al. 2019).The high levels of genetic diversity favor rapid adaptation, and the result of this local adaptation is the development of a phenotype that optimizes the response to environmental pressures, and is associated with the highest fitness (Kremer 2010).
Q. petraea, Q. frainetto and Q. pubescens belong to the Quercus (white oak) group (Nixon 1993) and their populations are rather common in Serbia, in the Western Balkans, where they occasionally occur in sympatry.This is the case with the site at Mt. Kosmaj, which acquired the status of the Landscape of Outstanding Features, due to its species richness and diversity, in addition to cultural, historical, and geological features.Q. petraea, Q. frainetto and Q. pubescens are autochthonous in this area (Stajić et al. 2019), and have not been studied to date at the molecular level.We analyzed 160 individuals of these three species with 14 nuclear microsatellites in order to assess the levels of genetic diversity in their populations and their genetic structure, needed for the formulation of conservation strategies and management practice.
The parameters of genetic diversity obtained in our study are consistent with previous reports in which oak species were analyzed with the same type of molecular markers (e.g., Muir and Schloetterer 2005, Salvini et al. 2009, Neophytou et al. 2010, Alberto et al. 2010, Katičić Bogdan et al. 2018).The values of expected heterozygosity for the three oak species were also recorded by different authors across Europe -e.g., ranging from 0.781 to 0.815 in Q. petraea from central Europe and Balkan Peninsula (Neophytou et al. 2010), from 0.251 to 0.890 in Q. pubescens populations from Italy (Di Pietro et al. 2020), or from 0.701 to 0.929 in Q. pubescens and from 0.181 to 0.922 in Q. frainetto from Romania (Curtu et al. 2011).As our study obtained the values of expected heterozygosity of 0.824 in Q. petraea, 0.834 in Q. pubescens, and 0.794 in Q. frainetto, it can be concluded that the obtained values are in accordance with the values obtained in previous studies.Thus, we found that all three analyzed species are characterized by rather high levels of genetic diversity, which indicates good prospects for their long-term survival, especially in conditions of changing climate and habitat degradation.
Despite living in sympatry, all three species of oaks have retained their genetic integrity, which is a rather important finding relevant for the conservation and management practice.Low but statistically significant genetic differentiation, expressed via F ST values, was observed (0.047 among Q.petraea and Q. frainetto, 0.039 among Q.pubescens and Q. frainetto, and 0.032 among Q.petraea and Q. pubescens).Furthermore, Q. frainetto was clearly separated from the other two oak species in PCoA analysis along the first coordinate, while Q. petraea and Q. pubescens were clearly separated in the PCoA analysis, along the second coordinate.Also, in the case of Q. frainetto, almost all individuals were strongly assigned to one gene pool, while individuals belonging to Q. petraea and Q. pubescens were strongly assigned to two or three distinct gene pools, indicating complex substructure of their populations.Population substructure is commonly associated with limitations to the gene flow (Sork 2016), and introgression, which is rather common in oaks, known for interspecific hybridization (e.g., Curtu et al. 2007, Salvini et al. 2009, Neophytou et al. 2010, Ortego and Bonal 2010).The highest gene flow and the lowest genetic differentiation were observed between Q. petraea and Q. pubescens, and the lowest gene flow and the highest genetic differentiation between Q. petraea and Q. frainetto.These results were expected because lower genetic differentiation and higher levels of gene flow are commonly observed among species that are more closely related, such as Q. petraea and Q. pubescens, than those that are more distantly related, such as Q. petraea and Q. frainetto (Curtu et al. 2007).A similar pattern has been previously observed in Italian (Salvini et al. 2009, Fortini et al. 2015) and Romanian (Curtu et al. 2007, Curtu et al. 2011) populations of oak species that are more or less related.The substructure in Q. petraea and Q. pubescens could exist due to higher gene flow and potential presence of introgressive forms in both species.
An important finding is rather high effective population size in all three oak species, which was the highest in Q. frainetto (∞, 1859.9;∞) and the lowest in Q. pubescens (305.5, 190.2; 719.2).This finding suggests a rather high number of parent's contribution to the formation of the next generation, which is important for the maintenance of high levels of genetic diversity in next generations.Nevertheless, low but statistically significant inbreeding was observed in all examined oak populations in the study area.The observed excess of homozygotes, however, most likely reflects population substructure, i.e., Wahlund effect (Wahlund 1928).It is well-known that the variation of allelic frequencies among subpopulations may create a heterozygote deficiency at the scale of the whole population (Bacillieri et al. 1994).This is supported by the outcomes of the STRUCTURE analysis which revealed population substructure in populations of Q. petraea and Q. pubescens, which comprised individuals strongly assigned to two or three distinct gene pools.It is worth mentioning that the lowest F value was observed in the population of Q. frainetto (F = 0.065, SE = 0.041), for which population substructure was not found in the STRUCTURE analysis.

Figure 1 .
Figure 1.Distribution of the selected oak trees used in the study in the Landscape of Outstanding Features "Kosmaj", with the representation of the protection zones (II and III).

Figure 1 .
Figure 1.Distribution of the selected oak trees used in the study in the Landscape of Outstanding Features "Kosmaj", with the representation of the protection zones (II and III).

Figure 5 .
Figure 5.The results of the STRUCTURE analysis under the assumption of four genetic groups (K = 4).The proportion of each cluster group for each individual is shown by the color code.

Figure 6 .
Figure 6.The results of the STRUCTURE analysis under the assumption of six genetic groups (K = 6).The proportion of each cluster group for each individual is shown by the color code.

table 1 .
A list of 14 informative nuclear microsatellites used in the study, with primer sequences for their amplification, repeat motif, allele range available in literature, and references.
https://www.seefor.euSEEFOR 14(2): early view 5 N -population size; A -number of alleles; Ae -effective number of alleles; H O -observed heterozygosity; H E -expected heterozygosity; F -coefficient of inbreeding; SE -standard error N -population size; A -number of alleles; PA -number of private alleles; Na -average number of alleles per locus; Ae -average effective number of alleles per locus; Ar98 -allelic richness according to the rarefaction method for 98 gene copies; Ne -effective population size; 95% CI -95% confidence intervals; H O -observed heterozygosity; H E -expected heterozygosity; F -coefficient of inbreeding; SE -standard error

table 2 .
Summary of the genetic diversity parameters based on genotyping with 14 nuclear microsatellite loci in the three oak species.

table 3 .
Standard genetic diversity parameters in populations of the three oak species.

table 4 .
The F ST values and Nm between the pairs of populations of Q. petraea, Q. pubescens and Q. frainetto.