Genetic variability of indigenous (Quercus robur L.) and late flushing oak (Quercus robur L. subsp. slavonica (Gáyer) Mátyás) in adult stands compared with their natural regeneration

Slavonian oak (Quercus robur subsp. slavonica (Gáyer) Mátyás) is currently gaining interest in forestry due to forest restructuring in Germany caused by climate change. Slavonian oaks originating from Croatia have been introduced into Germany mainly in the Münsterland region of North Rhine-Westphalia since the second half of the nineteenth century. They are characterized by their late bud burst, long clear bole, stem straightness and faster height and diameter growth compared to indigenous oaks in Germany. In this study, the genetic differentiation of adult trees and their respective progeny of two Slavonian and two indigenous stands in Hamm-Westtünnen, was evaluated. Genetic diversity and structure were estimated using 23 nuclear simple sequence repeat (SSRs) and 5 maternally inherited chloroplast microsatellite markers (cpSSRs). The mean expected heterozygosity of 0.545 and allelic richness of 6.23 indicate high genetic diversity in the studied populations. The group of progenies (AR = 8.40, Ho = 0.524, He = 0.559, FIS = 0.064) shows similar levels of genetic variation as the adult stands (AR = 8.37, Ho = 0.513, He = 0.554, FIS = 0.075). The genetic differentiation between adult stands and progeny was low (FST = 0.013). Genetic assignment of individuals using STRUCTURE revealed that the studied populations were divided into two clusters. There was no evidence of extensive hybridization or gene flow between Slavonian and native populations, possibly due to the different timing of bud burst of the two taxa.


Introduction
In addition to forest fragmentation and urbanization, more frequent extreme weather events and the associated spread of pests and diseases, as well as the maladaptation of tree populations that have been translocated to new sites, threaten existing forests (Ivetić et al. 2016). Due to global warming and associated temperature extremes, forest trees will have to adapt to the major climate changes in the future. Thus, translocated populations also offer the possibility to study adaptation processes and to evaluate their suitability in the face of climate change. The adaptability of a plant is determined by its genetic material and by environmental factors at its growth site (Finkeldey and Hattemer 2010). A prerequisite for both this adaptation and for maintaining adaptive capacity for future generations is sufficient genetic variation in populations (Vornam et al. 2004;Gailing et al. 2008). The conservation of genetic variation plays an important role in forestry in the context of restructuring towards more climate-stable forests in order to sustain vital and productive stands. Therefore, the response of forest trees to climate change is determined by the interplay of genetic variation, adaptation, gene flow, and interspecific hybridization (Kremer 2010).
White oaks (section Quercus, family Fagaceae) are dominant and diverse forest tree species widely distributed across the Northern Hemisphere (Nixon 1993;Aldrich and Cavender-Bares 2011;Leroy et al. 2020). With regard to 1 3 climate change, it is assumed that oaks will play an increasingly important role in forests, as they have a high adaptation potential as well as drought tolerance and could thus respond positively to a deficit in the climatic water balance. Due to their major ecological as well as economic significance, and their propensity for hybridization and gene flow (e.g. , Burger 1975;Rushton 1993;Lepais et al. 2009;Kremer and Hipp 2020) oaks have been intensively studied, particularly with regard to their genetic diversity and especially their geographic structure (Dumolin et al. 1995;Curtu et al. 2007;Aldrich and Cavender-Bares 2011).
Pedunculate oak (Quercus robur L.) is one of the dominant oak species in Germany as well as most of Europe (Barreneche et al. 1998;Reif et al. 2016). Within the species Quercus robur, numerous taxonomic subspecies or ecological varieties have evolved over time in geographic sub-areas throughout its range, such as Quercus robur subsp. slavonica (Gáyer) Mátyás with distribution in Slavonia (eastern region of Croatia) . Slavonian oak is in demand as imported timber throughout Europe, mainly because of its very good qualities (Rieger 2018), and is currently gaining interest for forestry due to forest restructuring in Germany caused by climate change . Especially from a yield point of view, it has a great number of remarkable characteristics like high growth performance, straightness, a long clear bole as well as fine branches (Wachter 2001;Gailing et al. 2003). Slavonian oaks have been introduced into the western part of Germany, especially in the region around Münster, in the second half of the nineteenth century with the beginning of extensive seed trade through steam engines (Wachter 2001;Gailing et al. 2007a). Accordingly, the Slavonian oak stands are first generation stands in Germany.
Over the last two decades, the development and accessibility of several molecular markers enabled precise characterization of forest genetic resources as well as the study of hybridization and gene flow among oak species (Crăciunesc et al. 2013;Lind and Gailing 2013). For instance, nuclear microsatellites (nSSRs) can be used to study genetic differentiation and structures as well as gene flow between related species (e.g., Lexer et al. 2005;Lepais et al. 2009;Lind and Gailing 2013;Pérez-Pedraza et al. 2021), and chloroplast DNA markers (cpSSRs) for studying postglacial recolonization patterns (e.g., Petit et al. 2002b) in order to describe past episodes of introgression.
Here, we employed molecular markers (nSSRs and cpSSRs) to examine the genetic variation in adult stands and their natural regeneration and the presence of gene flow between populations of Quercus robur subsp. robur and Quercus robur subsp. slavonica. The present paper will address the following questions: (1) Is there any evidence of gene flow between the two closely related oak taxa, based on genetic assignment analysis comparing adult and offspring generations? (2) Is the haplotype composition of the natural regeneration in accordance with the one of the adult stands?
(3) Does the amount of genetic variation vary among varieties and their natural regeneration?

Plant material and DNA isolation
Plant material (leaves) of in total 270 trees from four neighboring stands (SLAV_161 A2, SLAV_161 A1, IND_160 B, IND_159 B) in Hüls (Hamm-Westtünnen, NRW) owned by Freiherr v. Boeselager was collected by shooting branch tips in early July 2019, Fig. 1, Table 1). Samples for SLAV_161 A2, SLAV_161 A1 and IND_160 B were randomly distributed across the plots, while for IND_159 B only the northern part of the plots were sampled (Fig. 1). Due to topographical and local conditions, as well as a very heterogeneous inventory picture, the northern part of 159 B was regarded as an independent stand from the rest of the entire stand. The stands are between 125 and 133 years old (Table 1). In addition, leaves of 369 seedlings from the natural regeneration were sampled in all stands at 7-12 sample points evenly distributed throughout the sampled part of each stand. Gailing et al. (2003) characterized the four stands based on 10 individuals at chloroplast -DNA markers and classified them at the tree level into early, intermediate, and late flushing according to Wachter (2001) (Table 1). The stands IND_159 B and IND_160 B showed an early bud burst and also the indigenous haplotype 1 (only one tree in IND_159 B showed haplotype 2), while the trees in the stand SLAV_161 A1 consistently sprouted later (Gailing et al. 2003). Stand SLAV_161 A1 was characterized by the Slavonian haplotype 2 (8 individuals) and haplotype 5 (two individuals, haplotype present in both Germany and Croatia). Stand SLAV_161 A2 showed a heterogeneous picture with early and late sprouting trees (intermediate) and was characterized by the haplotypes HP2 (4 individuals), HP5 (4 individuals) and HP10 (two individuals) (Table 1). Based on morphological and phenological classification of stands according to Wachter (2001) and genotyping with chloroplast markers by Gailing et al. (2003), two indigenous (IND_159 B, IND_160B) and two Slavonian (SLAV_161 A1, SLAV_161 A2) were identified.
The genomic DNA was extracted from in total 749 freshly collected leaves (approx. 1 cm 2 piece of leaf) using the DNeasy 96 Plant Kit from Qiagen (Hilden, Germany). Subsequently, the amount of DNA was checked on 1% agarose gels.
Thereof, 11 primer pairs were originally developed for Quercus rubra, of which 7 were derived from a Quercus rubra EST-library (https:// www. hardw oodge nomics. org/ Trans cript ome-assem bly/ 19630 23? tripal_ pane= group_ descr iption_ downl oad) and the other four markers were from Aldrich et al. (2002) and Sullivan et al. (2013). The other 12 markers were originally developed for Quercus robur by Durand et al. (2010). We chose these 23 markers because of the high transferability of EST-SSRs between Quercus species, and the availability of a large number of nSSRs and EST-SSRs developed for Q. robur, Q. petraea, and related species (Steinkellner et al. 1997;Barreneche et al. 1998;Tanaka et al. 1999;Durand et al. 2010) that have been used in previous studies on indigenous and Slavonian oaks (e.g., . The annotation of the sequences was obtained by searching the individual primer sequences in the respective contigs to identify the complete contig sequences for similarity searches against the UniProt Viridiplantae database (The UniProt Consortium 2021) using BLASTx (Basic Local Alignment Search Tool) (Altschul et al. 1990).
We used the same multiplex reactions for PCR as described in  with the exception that 3 additional primers (one in multiplex 2 and two in multiplex 3) were used.

Chloroplast microsatellites (cpSSRs)
The chloroplast microsatellites udt1, ucd4 and udt4 specifically developed for oaks (Deguilloux et al. 2003) (Table S. 1), were analyzed and amplified in all samples. Ucd4 and udt4 were used to differentiate between the main haplotypes (Table S. 2). In addition, udt1 was used to distinguish the southwestern European chloroplast haplotypes HP10-11 and HP12 from each other. Finally, the universal chloroplast markers ccmp2 and ccmp10 (Weising and Gardner 1999) were used for additional validation of 28 samples with haplotype 17 (Table S. 2).
The haplotypes HP5 and HP7-26 showed the same length at all used cpSSRs and could therefore not be distinguished (Gailing et al. 2007a, b) (Table S. 2). Consequently, we consider them as one haplotype HP5/7-26.
Each PCR reaction of the udt and ucd cpSSRs was also performed in a Biometra Thermal Cycler (MJ Research PTC 200, Analytik Jena, Germany) using the following touchdown program: An initial denaturation at 95 °C for 15 min, followed by 8 touchdown cycles at 94 °C for 1 min, 1 min at 53 °C (decreasing 1 °C each cycle) and 1 min at 72 °C, followed by 33 cycles at 94 °C for 1 min, annealing at 45 °C for 1 min and elongation at 72 °C for 1 min, and a final extension step at 72 °C for 20 min. The PCR profile for the universal cpSSRs (ccmps) consisted of an initial denaturation at 95 °C for 15 min, followed by 35 cycles at 94 °C for 1 min, a 1 min annealing step at 53 °C, an elongation step at 72 °C for 1 min and a final extension step at 72 °C for 10 min.
The forward primers used for PCR reactions are labelled with the fluorescence dyes 6-FAM or HEX (Table 4). The amplification products were diluted (1:200) and separated on an ABI 3130xl Genetic Analyzer (Applied Biosystems, Foster City, USA) in a multiplex analysis. Fragment sizes were determined with the GeneScan™ Rox-500 size marker using the software GeneMapper® version 4.1 (Applied Biosystems, Foster City, USA).

Statistical analysis of genetic diversity and differentiation based on nSSRs
The genetic variation indices observed heterozygosity (H o ) and expected heterozygosity (H e ), as well as the fixation index F ST and Hedrick's standardized G ST (G' ST (Hed)) for individual markers and across all markers were calculated in GenAlEx version 6.51b2 (Peakall and Smouse 2006;Smouse et al. 2017). With the software Fstat version 2.9.4 (Goudet 2003) allelic richness (A r ) (here: based on a minimum sample size of 48 individuals) as well as the inbreeding coefficients (F IS ) (Weir and Cockerham 1984) were determined. Differences in levels of A R , H o , H e and F IS between adults and progeny were analyzed using the 'comparison among groups of samples test' in Fstat computing 1000 permutations. In addition, significant deviations from zero of F IS were determined after standard Bonferroni correction (α = 0.05, p < 0.00022) implemented in the software Fstat (Goudet 2003). Linkage disequilibrium (LD) was calculated for each pair of loci in the 8 populations using Genepop version 4.7.2 (Rousset 2008) based on the following settings for the Markov chain: dememorization 10,000, batches 100 and iterations per batch 5000.
A dendrogram with UPGMA based on Nei's standard genetic distance D S (Nei 1972) was calculated using the software program Populations version 1.2.32 (Langella 1999), with bootstrap values based on 1000 permutations across loci. Applying the software TreeView version 1.6.6 (Page 1996), the dendrogram was visualized.
To detect signatures of recent genetic bottlenecks and to investigate whether the introduced seeds established from seeds representing a limited number of seed parents, the software BOTTLENECK version 1.2.02 (Piry et al. 1999) was used. Therefore, we performed a Wilcoxon signed-rank test (one-tailed) for heterozygosity excess for 'the two phasemutation model' (TPM) and 'the stepwise mutation model' (SMM) and a mode-shift analysis to test for a distortion in the allele frequency distribution.
Genotyping errors due to non-amplified alleles (null alleles) which can lead to overestimates of the inbreeding coefficient (F IS ) were identified using the Windows®-based software MicroChecker version 2.2.3 (Van Oosterhout et al. 2004). Furthermore, Arlequin version 3.5.2.2 (Excoffier and Lischer 2010) was run with 50,000 simulations of 100 demes per group with the infinite island model based on F ST in order to detect outliers which deviate significantly from the variation and differentiation expected under neutrality among all populations and between the groups of indigenous and Slavonian oaks.
An analysis of molecular variance (AMOVA) was used to analyze population differentiation and performed with GenAlEx using 9999 permutations.

Analysis of population structure
Population structure was inferred using model-based cluster analysis implemented in the software STRU CTU RE version 2.3.4 (Pritchard et al. 2000) to identify possible subpopulations based on the microsatellite dataset. For this purpose, the admixture model was chosen because of its consideration of the ancestry of admixed individuals, and correlated allele frequencies were selected because of their potential to improve clustering in closely related populations, while using a burn-in period of 50,000, Markov chain Monte Carlo (MCMC) replicates of 100,000, and the LOCPRIOR model. However, we used additionally the default setting in STRU CTU RE 'admixed model without LOCPRIOR' to identify population structure solely based on genetic information. In addition, to determine the "best K" from the logarithmic results and the ΔK method (Evanno et al. 2005), the online software STRU CTU RE HARVESTER v. 0.6.94 (Earl and vonHoldt 2012) was used. The CLUMPAK software (Cluster Markov Packager Across K) was applied for post-processing and graphical representation of the results of the model-based population STRU CTU RE analysis (Kopelman et al. 2015).
Besides, clustering of populations was analyzed using a principal coordinate analysis (PCoA) with the pairwise, individual-by-individual (N x N) genetic distance implemented in GeneAlEx in order to find and plot the major patterns within this dataset Smouse 2006, 2012).

Chloroplast microsatellites (cpSSRs) and the frequency of cpDNA haplotypes
A total of 5 haplotypes could be distinguished in the studied populations corresponding to the haplotypes HP1, HP10-11, HP2, HP5/7-26, and HP17 (Petit et al. 2002b). Haplotype 1 is the most frequent haplotype in western Germany and dominant in Central Europe, but completely absent in the Balkan region (Bordács et al. 2002;König et al. 2002;Petit et al. 2002a, b). Its glacial refugium was most likely located in southern Italy (Petit et al. 2002a). The haplotypes HP10-11 show a center of distribution in southwestern and western Europe and had their glacial refugia on the Iberian Peninsula (Petit et al. 2002b). Moreover, these haplotype (HP10-11) are also apparently missing in the Balkan region (Petit et al. 2002a). In addition to these presumably indigenous haplotypes, HP2, HP5/7-26 and HP17 have a center of distribution in the Balkan region with presumed southern Italian origin (Bordács et al. 2002;Petit et al., 2002a, b). Although HP2, HP5/7-26, and HP17 are very frequent in the Balkan region (Petit et al. 2002a), none of them are restricted to this region, but only haplotype 5/7-26 occurs naturally in oak populations in both Germany and the Balkan region.
The haplotype composition of each of the four stands (two Slavonian and two indigenous) and their natural regeneration is summarized in Table 3 and visualized in Fig. 2.

Variation within and among stands at nuclear microsatellite markers
Inbreeding coefficients across all markers were not significantly different from zero in any population.  Table 4.
The mean values for effective number of alleles, allelic richness and inbreeding coefficients (not significant after Bonferroni correction (5%)) were lower in the natural regeneration (A R : 6.02, H e : 0.546, F IS : 0.045) compared to adult stands (A R : 6.43, H e : 0.545, F IS : 0.068). However, the observed and expected heterozygosity was slightly lower in the group of adult stands (H o : 0.513, H e : 0.545) than in the natural regeneration (H o : 0.524, H e : 0.546).
The mean observed heterozygosity per locus ranged from 0.092 at locus 2P24-0.818 at locus GOT040, and the mean expected heterozygosity per locus ranged from 0.088 at locus 2P24 -0.879 at locus OC11 (Table S. 3).
The pooled sample of naturally established seedlings across all stands showed no significantly higher (exact test, 1000 permutations, p = 0.045) allelic richness (A R ) in adult trees (6.43) than in the seedling cohort (6.02; Table 4). Higher levels of A R in adults compared with seedlings were detected for each population separately. In addition to the mean A R across populations, observed heterozygosity (H o ), expected heterozygosity (H e ) and inbreeding coefficient (F IS ) were also not different between the adult and progeny cohort (exact test, 1000 permutations, p > 0.05) ( Table 2).
The genetic differentiation (F ST ) between the 8 populations (adult and natural regeneration) was relatively low at most loci ranging from 0.009 for FS_C2791 -0.174 for Qr0332. Overall, the mean differentiation among all populations was 0.022 (G' ST (Table S. 4). In addition, Table S. 5 shows also the pairwise F ST values between adult and progeny separately for Slavonian and indigenous stands. The low differentiation between adult and offspring generations as well as between indigenous and Slavonian oaks and results of genetic assignment analyses reveal similar genetic structures in both generations and no evidence for extensive gene flow between Slavonian and indigenous stands.
Further, the outlier test with Arlequin showed that there are no markers with signature of selection among populations and between clearly indigenous (IND_160_B) and Slavonian (SLAV_161 A1) stands ( Fig. S. 1).
The percentage of all pairs of loci in LD was 0.79% in the natural regeneration SLAV_NR_161 A2, as well as 0.40% in the natural regenerations SLAV_NR_161 A1, IND_NR_159 B and IND_NR_160 B (mean over all populations: 0.20%). All adult populations had no loci in LD.
Null alleles were detected in each of the 10 populations at least at one of the markers OC11, 3D15, FIR028, FIR035, FIR043, GOT040, PIE040 (highest null allele frequencies in IND_NR_159 B: 0.1716, in SLAV_NR_161 A1: 0.1848 and IND_159 B: 0.1716), PIE125, VIT023, VIT107, FS_C2660 and FS_C2791 (highest null allele frequency in SLAV_161 A1: 0.1615) (Table S. 6). However, in all populations null alleles only occurred with a frequency between 0.51% and 2.51% over all loci (Table S. 6). Furthermore, null alleles across all populations per locus ranged from 0.006 at locus PIE125 and GOT040 to 0.092 at locus PIE040 (Table S. 3). The comparison of null alleles between marker types showed that two of the four gSSRs (OC11, 3D15) showed null alleles with a low frequency (< 0.040), null allele frequencies of EST-SSRs developed for Q. robur and Q. petraea ranged from 0.006 at locus PIE125 to 0.092 at locus PIE040 and null allele frequencies of EST-SSRs developed for Q. rubra ranged from 0.013 at locus Qr0057 and FS_C2660 to 0.016 at locus FS_C2791 (Table S. 3). Across all loci, EST-SSRs developed for Q. robur and Q. petraea have the highest null allele frequency (0.021), followed by gSSRs (0.012) and EST-SSRs developed for Q. rubra (0.006).
The results of the program BOTTLENECK mode-shift analysis showed that the allele frequencies of all populations followed a normal L-shaped distribution as expected for mutation-drift equilibrium (Table S. 7). This may indicate that the introduced seeds did not originate from only a limited number of seed parents. The AMOVA overall populations showed that 98% of the molecular variance was within populations (91% within and 7% between Individuals) and only 2% among populations.

Population structure
Using principal coordinate analysis (PCoA), the basic patterns among the populations are presented in Fig. 3. Here, the PCoA results do not show a clear differentiation between the indigenous and the Slavonian populations. Even though there is a separation along axis 1, the natural regeneration of IND_159 B groups with Slavonian stands (Fig. 3, see also Fig. 4). In comparison, PCoA based on chloroplast microsatellite markers shows that each natural regeneration is grouped with their adult stand (Fig. 5).
Using the Evanno method (Evanno et al. 2005) and STRU CTU RE HARVESTER (Earl and vonHoldt 2012) an optimal and significant K = 2 was determined ( Fig. S. 2, Fig. S.  Fig. 4. In addition, Table 5 shows the proportion of membership of each sampled population in the two clusters, but there is no specific cluster for each origin, since both clusters occur in both origins. A proportion of ancestry > 0.5 or < 0.5 in cluster 1 is characteristic for indigenous and Slavonian stands, respectively. In addition, Fig. 6 combines the STRU CTU RE results with the haplotypes and shows that the Slavonian haplotypes (HP2, HP5/7-26, HP17) are differentiated from non-Slavonian HP1, HP10-11 at nSSRs (Fig. 6). Furthermore, a table in the Online-Resource 2 lists a direct comparison of cpSSR haplotype vs. STRU CTU RE-based proportion at tree level for each population. However, we found no putative F 1 hybrids. Because of the low differentiation between Slavonian and indigenous oaks at nSSRs, the identification of F 1 hybrids is challenging. While adult stand and natural regeneration of stand IND_160 B are clearly assigned to the "indigenous" cluster, seedlings and adult trees of IND_159B  show a variable assignment which is also reflected in the PCoA analysis and cpSSR haplotype composition (see above, Fig. 3).
To show the genetic relationship between the 8 studied populations (4 adult stands and their natural regeneration), an UPGMA dendrogram was also prepared in the appendix, based on the unbiased genetic distance matrix of Nei (1972) (Fig. S. 4). In accordance with PCoA (Fig. 3), it shows that the Slavonian and native stands differ only slightly at nuclear markers.

Genetic diversity of adults and natural regeneration
To our knowledge, nuclear microsatellites are among the most precise tools for genotype determination and are important genetic markers widely used in population genetic analysis in oaks. According to Spence et al. (2021), microsatellites are an accessible method to gain an understanding of genetic diversity and structure.
A previous study ) conducted on a larger number of populations enabled the assessment of genetic differentiation between native and Slavonian adult stands using the same nuclear microsatellite marker set (except for 3 additional markers used in the present study). Accordingly, our dataset also shows high and similar genetic diversity at the microsatellite loci examined (Table 4 and Online Resource 1) across all populations with expected heterozygosity (H e ) of 0.545 and allelic richness (A R ) of 6.23. Such a relatively high genetic variation is common for outcrossing and wind pollinated woody species such as oaks (e.g., Shi et al. 2017: H e = 0.707, A R = 7.79; Spence et al. 2021: H e = 0.73, A R = 5.29). Furthermore, the overall pattern of genetic variation we observed at nuclear microsatellite loci is typical of longlived and outcrossing species such as trees with relatively high genetic diversity (A R = 6.23, H e = 0.545), a low level of inbreeding within (F IS = 0.059) and a low differentiation among populations (F ST = 0.022) (Hamrick et al. 1992). According to Kremer and Petit (1993) such a low species differentiation is typically for European white oaks.
No clear differences between the genetic variation indices of the natural regeneration and the adult trees were found in the stands of Hamm-Westtünnen in North-Rhine Westphalia (Table 4). Other studies, however, showed higher genetic diversity estimates, but also no clear differences between adult and progeny of pedunculate oak stands. For instance, based on 10 genomic microsatellite loci, Vranckx et al. (2014) reported for the adult trees a mean number of A R of 11.5, a mean H o of 0.798, a mean H e of 0.821, and for the progeny a mean A R of 10.7, a mean H o of 0.804, and a mean H e of 0.817. Based on 19 nuclear microsatellite markers (12 genomic and 7 EST-SSRs), also Sandurska et al. (2017) reported slightly higher values (adult: A R = 17.823, H o = 0.784, H e = 0.821; seedlings: A R = 17.128, H o = 0.759, H e = 0.803) as compared to the present study. However, in assessing diversity indices from different studies, it should be mentioned that the values are not directly comparable, as different microsatellite marker sets have been used by different  . 6 Differentiation of haplotypes based on nSSRs (admixed model with LOCPRIOR) displayed using CLUMPAK (Kopelman et al. 2015).
Populations separated by black vertical lines. More orange = indigenous haplotype. blue = Slavonian haplotype researchers. For instance, higher genetic variation could be due to the generally higher variability at genomic SSRs compared with EST-SSRs (Rungis et al. 2004). The estimates of F IS were positive and show significant deviations from HWE at some loci (Online-Resource 1). Positive F IS values indicate inbreeding only if they are significant at most loci, which was not the case here (only at 2 loci per population, Online-Resource 1). Otherwise, high F IS values at only certain loci may indicate selection for these loci, selection for genes associated with these loci, or the occurrence of null alleles (Katičić et al. 2018). However, due to the very low null allele frequency (0.51%-2.71%) on average across loci, little impact on population genetic structures and differentiation (F ST ) can be assumed based on Oddou-Muratorio et al. (2009) and Chapuis und Estoup (2007).
Positive and higher F IS in the Slavonian adult and offspring cohorts as compared to the indigenous stands could be due to the Wahlund effect and non-random mating (e.g., inbreeding, or preferential mating within groups of different geographic origin, differences in flowering time). However, no clear effects of genetic drift and inbreeding on the genetic diversity of the progeny could be observed, so non-random mating is rather unlikely.
Reproductive isolation of late flushing Slavonian oaks from indigenous stands has likely resulted in the observed pattern of genetic differentiation between Slavonian and indigenous seedling cohorts. Here, it is important to note that the flushing date was assessed a proxy for flowering time, two traits which are highly correlated in oaks (e.g., Chesnoiu et al. 2009). Low coefficients of pairwise fixation levels show also low difference between the two generations (F ST = 0.013) of the studied stands, thus providing no evidence for fixation or genetic drift. In addition, the overall effect of genetic drift on allelic richness (A R ) was also rather limited, and even when all populations were pooled, allelic richness was not significantly different between adult and seedling cohorts (Table 4).
Likewise, the analysis of molecular variance for the two generations shows greater variability within populations (98%) than among populations (2%), which is commonly observed in pedunculate oaks (Gömöry et al. 2001;Neophytou 2015;Kesić et al. 2021). Furthermore, a comparison between Slavonian (A R = 6.26, The low genetic differentiation between generations and similar levels of expected heterozygosity and allelic richness in the seedlings of Slavonian stands indicated that natural regeneration is efficient in these stands, possibly due to the low density of herbivores/deer. According to Hamrick (2004), sustaining high genetic diversity increases the likelihood of having well-adapted genotypes to ensure the long-term viability of forest tree populations.

Haplotype composition of adult stands and their natural regeneration
Since chloroplast markers (cpSSRs) in angiosperms are usually maternally inherited (Dumolin et al. 1995) and natural regeneration share the same chloroplast information as the maternal tree, they show high uniformity within stands but marked differentiation among them, which is especially true for barochoric tree species such as oaks (Petit et al. 2005). Thus, chloroplast markers are particularly suitable in combination with biparentally inherited microsatellite markers to classify and detect admixture of natural regeneration of Slavonian and indigenous oaks. The five cpSSRs used in this study have been shown to be suitable and useful in previous studies to identify different haplotypes (Gailing et al. 2007a(Gailing et al. , b, 2009. The populations analyzed in this study were first examined for haplotype composition in 2003 by Gailing et al. (2003) on 10 representative individuals per population. Thereby, haplotypes 1-4 found in Gailing et al. (2003), correspond to haplotypes 2, 5, 10 and 1 described in Petit et al. (2002a). We found similar haplotype composition for the SLAV_161 A1, SLAV_161 A2, IND_159 B, and IND_160 B stands as described in Gailing et al. (2003), except that HP1 and HP17 were additionally found in individual trees in SLAV_161 A2 and HP2, HP5/7-26, and HP17 in IND_159 B (Table 3, Fig. 2). Since haplotype 5/7-26 also occurs naturally in Germany, it is not possible to infer Slavonian reproductive material from cpSSR markers alone. Therefore, the combination of both marker types (maternally inherited cpSSRs and biparental inherited nSSRs) was useful to detect this admixture and to assign haplotypes, which occur naturally in both Germany and Slavonia, to their geographic origin.
The natural regeneration shows a similar to equal haplotype composition compared to the respective adult stands (Table 3). This can be also seen in the PCoA based on cpSSRs (Fig. 5), which assigns each natural regeneration to its respective old-growth stand with low genetic distances. Only the relative frequencies of haplotypes differ slightly between adults and seedlings (Fig. 2), but the most frequent haplotype remains the same in both generations (Table 3).

Genetic assignment of populations
We used the Bayesian model-based clustering method implemented in the program STRU CTU RE (version 2.3.4, Pritchard et al. 2000) to assign individuals to K populations based on 23 nuclear microsatellites. The STRU CTU RE analysis (Fig. 4, Table 5) revealed that the studied populations were grouped into two clusters, Slavonian (a higher proportion of blue) and indigenous oaks (a higher proportion of orange). However, it is difficult to identify two distinct groups (Slavonian and Native) based on the PCoA (Fig. 3), since the natural regeneration of each adult stand is shifted substantially against their maternal stands in one direction. On the one hand, a possible reason for the shift in natural regeneration IND_NR_159 B and IND_NR_160 B along axis 1 could be gene flow between native stands (e.g., the surrounding native stands that were not sampled, cf. Figure 1). On the other hand, the shifts in the natural regeneration of the Slavonian stands could suggest a (non-random) gene flow (pollen influx) from another Slavonian stand (Table S. 4, 5, Fig. S.4). However, a PCoA with all Slavonian and native populations in North Rhine-Westphalia from  and the populations investigated in this study also revealed two groups in which the Slavonian populations (SLAV) could be assigned to the Slavonian cluster of other Slavonian populations in North Rhine-Westphalia ( Fig. S. 5). Therefore, the nuclear marker set is well suited to distinguish at the population level (but not at the tree level) between Slavonian and native oak in both old-growth stands and their progeny and to identify mixed stands (see stand IND_159 B). In addition to these nuclear marker results, the inclusion of cpDNA haplotypes shows that Slavonian and indigenous oak stands are not always pure stands (Fig. 2). Thus, the Slavonian stand SLAV_161 A2 also showed 18% of the non-Slavonian haplotypes 1 and 10-11, and IND_159 B showed 25% Slavonian old-growth trees (HP5/7-26, 5, 17) (Fig. 2). We demonstrated that chloroplast markers are suitable to identify individual trees (tree level) with non-Slavonian haplotypes in Slavonian old-growth stands and vice versa.
Comparison of STRU CTU RE results with cpSSR haplotypes for haplotype 2 showed that within the parent population 85 trees had HP2 and of these 73 were also assigned to the Slavonian cluster (ancestry coefficient > 0.6) based on nSSRs (Online-Resource 2, Fig. 6). Thus, 86% of the parent trees with haplotype 2 were also assigned to the Slavonian variety based on nSSRs.
Furthermore, low genetic differentiation between Slavonian and indigenous oaks did not allow to unambiguously distinguish between varieties at the tree level (Online-Resource 2). The different timing of bud burst (Table 1), as described for these stands in Gailing et al. (2003) suggests that flowering times between taxa show little to no overlap (due to synchronization of flowering with flushing, Chesnoiu et al. 2009), thus restricting gene flow between Slavonian and indigenous stands. In Germany, Slavonian oaks show up to three weeks later bud burst compared to indigenous oak stands (Wachter 2001;Gailing et al. 2003).

Practical applicability
For sustainable forest management, it is important to provide forest enterprises with reproductive material that complies with both economic objectives, such as increment and wood quality, and ecological objectives, such as genetic diversity and local adaptation (Konrad et al. 2012). In this context, Slavonian oak is currently gaining interest as a seed source for forestry not only because of its very good economic characteristics (wood quality, fast height growth, long clear stem), but also with regard to climate change-induced forest conversion due to its broad site amplitude, possibly better adaption and late bud burst and associated lower susceptibility to the oak leaf roller (Schirmer 2017;Rieger 2018;. Currently, an alternative form of sustainable forest management is practiced, based on certification systems such as certification of the origin of reproductive material, which provides the forest owner with detailed information about the origin of the planting material. The microsatellite marker set (cpSSRs and nSSRs) used in this study can reliably be used for identification and certification of Slavonian and native forest reproductive material. For this purpose, a previous study  already indicated that this nuclear microsatellite marker set is suitable for identifying mixed stands, detecting the origin of reproductive material, and distinguishing between Slavonian and native oak stands.
Furthermore, the marker set is also suitable for establishment and test of admixture in seed orchards. Thus, Katičić et al. (2018) based on other nSSRs and the same cpSSRs were able to investigate clonal seed orchards in Croatia for their genetic diversity and differentiation as well as spatial genetic structure and transfer of forest reproductive material between clonal seed orchards.
In a further study, however, it would be necessary to investigate the extent to which gene flow occurs between 1 3 intermediate and early flushing stands. Especially seed stands must be chosen to avoid outside pollen influx. The advantage in this case could be that plantations or seed orchard with the particularly late-flushing haplotype 2 may not require geographic isolation from conspecific indigenous stands, as they are reproductively isolated from early-flushing neighboring stands due to different flowering times.

Conclusion and future perspectives
Our study shows that the genetic information of old-growth stands is almost completely passed on to the next generation. Minor changes in the genetic structure from adult to offspring generation suggest the maintenance of the same level of genetic variation in the progeny. Since genetic diversity has a strong impact on the adaptability of populations to changing environmental conditions and survival, the main objective for effective conservation of forest genetic resources is to achieve the highest possible level of genetic variability (Andonovski and Velkovski 2019). In addition, the genetic variability of tree populations ensures the stability and sustainability of forest ecosystems. Moreover, limited gene flow was revealed between the native and Slavonian oaks, likely due to differences in bud burst and flowering time between the two oak taxa, especially between the very late-flushing (predominantly haplotype 2) and early-flushing (predominantly haplotype 1) oaks.
In conclusion, seed production stands and seed orchards of late flushing Slavonian oak would require little geographic isolation from early flushing indigenous oak stands.
Author contributions KB wrote the main manuscript text and prepared all figures and tables. OG revised it critically for important intellectual content. All authors reviewed the manuscript.
Funding Open Access funding enabled and organized by Projekt DEAL.

Competing interests
The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.