Microsatellite analysis of genetic diversity in Czech populations of European beech (Fagus sylvatica L.)

We identified ten beech stands across the Czech Republic (Table 1, Figure 3) aiming at populations that are considered to be locally adapted and superior under local conditions. Each of the selected populations was represented either by 20 or 30 mature individuals according to the site characteristics, which turned into 250 tested individuals in total. Sampling was carried out in 2017. The minimal distance between sampled trees was 30 m to avoid sampling among closely related individuals (von Wühlisch, 2008). For each individual, we collected unburst buds and extracted the vascular cambium with a hole punch. Plant tissues were labeled appropriately and placed in plastic bags with silica gel. Then they were long-term stored in a freezer.

Initially, both the vascular cambium and leaf buds were tested for DNA extraction. Even though both types of tissues provided identical genetic profiles, we decided to extract the DNA from the buds, because of its higher purity and yield. From each individual, approximately 2 or 3 unsterilized fresh buds were cut into small pieces and placed into 2 mL safe-lock tubes (Eppendorf, Germany). Two 3 mm tungsten beads were added in each tube and frozen in liquid nitrogen before grinding (1 min, 30 Hz) using a mixer mill MM400 (Retsch, Haan, Germany). Genomic DNA was extracted with the Geneaid® Genomic DNA Mini Kit (Plant) according to the manufacturer's instructions with a minor modification in the lysis step (60 min of incubation phase). The DNA quantity and quality were measured with a NanoDrop® 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). The quality of the DNA was verified using 1.5% agarose gel electrophoresis. For subsequent Polymerase Chain Reaction (PCR) steps, we diluted all DNA samples to a working concentration of 20 ng/μL.

After amplification in single reactions, from preselected microsatellite loci, we excluded those which were difficult to evaluate or of low polymorphism (< 4 alleles), with poor amplification quality or low signal intensity. This led to a total number of 21 highly validated fluorescently labeled genomic markers – mfc5, mfc7 a mfc11 (Tanaka et al., 1999), FS1-03 (Pastorelli et al., 2003), sfc0036 (Asuka et al., 2004), concat14_A_0, EMILY_A_0, csolfagus_05, csolfagus_06, csolfagus_19, csolfagus_29, csolfagus_31, DE576_A_0, ERHBI_A_0, EEU75_A_0, DZ447_A_0 a DUKCT_A_0 (Lefèvre et al., 2012), Fagsyl_002929, Fagsyl_003849, Fagsyl_001018 and Fagsyl_003093 (Pluess & Määttänen, 2013) as presented in Table 2. The final loci selection was divided into two multiplexes (13- and 8-plex) according to their size ranges and PCR amplification compatibility. An issue of different annealing temperatures among selected primers was considerably decreased when using Type-it® Microsatellite PCR kit (Qiagen, Hilden, Germany) for PCR amplification. Both multiplexes showed clear patterns and consistent amplification profiles.

Amplification was performed in two multiplex polymerase chain reactions (Mastercycler® nexus, Eppendorf, Germany). The PCR mixture for multiplex 1 was performed in a total volume of 11.4 μL containing 1 μL template DNA (c = 20 ng/μL), 5.7 μL of Type-it® solution (Qiagen, Hilden, Germany), 0.26 μL sterile water and a total volume of 4.44 μL primer premix. The PCR for multiplex 2 was performed in a total volume of 6 μL, consisting of 1 μL template DNA (c = 20 ng/μL), 3 μL of Type-it® solution (Qiagen, Hilden, Germany), 0.14 μL of sterile water and 1.86 μL of primer premix. Concentrations and label dyes for each primer pair in the final primer premix are shown in Table 2. Both multiplexes had the same PCR conditions as follows: an initial incubation (5 min at 95 °C) was followed by 33 cycles consisting of the denaturation (45 s at 95 °C), annealing (45 s at 60 °C) and extension phase (60 s at 72 °C). Amplification was terminated with the final extension step (30 min at 72 °C). The fluorescently labelled PCR products, together with an internal size standard (GMC-GT500 LIZ) were processed on a 3500 Series Genetic Analyzer® (Applied Biosystems, Foster City, CA, USA).

Allele identification and genotyping were done using the GeneMarker® Fragment Analysis software (SoftGenetics, State College, PA, USA). All alleles were subsequently checked manually to minimize genotyping errors. Allele frequency analysis was performed using CERVUS v. 3.0.7 (Kalinowski et al., 2007) in order to identify the number of alleles per locus. Also, summary statistics were calculated for each locus, such as observed heterozygosity (H_o), expected heterozygosity (H_e), polymorphic information content (PIC) and fixation index (F_IS).

For the analysis of population genetic parameters, we used software GenAIEx v. 6.5 (Peakall & Smouse, 2006, 2012). To group populations according to their pairwise distances in a distance matrix, a Principal Coordinate Analysis (PCoA) was carried out. The distribution of total variance at different hierarchical levels was quantified with an Analysis of Molecular Variance (AMOVA). We also evaluated the discrimination power of assembled multiplexes. Therefore, the Probability of Identity (PI) for unrelated and related (PI_sibs) individuals was calculated, which represents a measure for the statistical power of individual identification with the chosen marker combination (Hartl & Clark, 2007).

All samples (n = 250) were genotyped on 21 SSR loci grouped into two multiplexes (Table 3). The number of genotyped individuals on each locus (N) varies because some of the samples did not show any apparent alleles on a particular locus. We did not detect vegetative copies among the samples. We observed the most significant dropout of samples (56 individuals) on locus Fagsyl_003093.

The total amount of 208 alleles was detected. The number of alleles per locus (k) ranged from 4 to 22 with a mean value of 9.905 ± 4.098. The highest count of alleles per locus (22) was found at mfc5 where we also observed significant deviation from Hardy-Weinberg equilibrium. The expected heterozygosity (H_e) ranged from 0.436 (csolfagus_29) to 0.906 (mfc5), with an average of 0.685 ± 0.134. The observed heterozygosity (H_o) ranged from 0.216 (Fagsyl_003093) to 0.859 (csolfagus_06) with an average of 0.625 ± 0.183. The mean polymorphic information content (PIC) across all loci was 0.649 ± 0.146. Most loci did not show a significant deviation from the expected Hardy-Weinberg frequencies, except 5 loci (mfc11, FS1-03, mfc5, Fagsyl_001018, Fagsyl_003093). We also detected 6 loci with a significant occurrence of null alleles (α = 0.05, marked bold in Table 3).

Private alleles were found in all analyzed populations except populations 5 and 10. The frequencies of rare alleles were generally low and ranged between 0.017 and 0.077, with an average of 0.033 ± 0.019. The highest amount of private alleles was observed in populations 1 and 9. Private alleles with the highest frequencies were found in population 2.

In order to evaluate the discriminatory power of used multiplexes, the Probability of Identity (PI, GenAlEx v. 6.5) was carried out. Figure 1 shows the decreasing probability of identity with an increasing number of markers used (e.g., the effect of cumulative markers).

Figure 1

In our study, we also evaluated a modified version PI_sib for estimating the probability of identity when related individuals are included in the sample. When using 5 randomly selected loci, the probability of statistically significant distinction was for unrelated individuals PI = 7.2 × 10⁻⁴ and for related individuals PI_sib = 0.04. With 8 loci used the probability decreases to PI = 3.13 × 10⁻⁷ and PI_sib = 2.2 × 10⁻². With 15 loci the probability further declines to PI = 4.97 × 10⁻¹⁵ and PI_sib = 3.78 × 10⁻⁶. For an entire set of 21 loci, the observed PI for unrelated individuals was PI = 1.85 × 10⁻²⁰ and for related individuals PI_sib = 2.77 × 10⁻⁸. These values confirm a high discriminatory power of used multiplexes even when family structures are expected among samples and the discriminatory power of loci is therefore significantly lower.

The parameter of observed heterozygosity (H_o) among populations ranged from 0.595 (Pop. 3) to 0.654 (Pop. 9) with the mean value of 0.625 ± 0.014. Expected heterozygosity (H_e) ranged from 0.651 (Pop. 6) to 0.678 (Pop. 8) with the mean value of 0.664 ± 0.010.

The Principal Coordinate Analysis was performed based on genetic distances between the populations (e.g. the pairwise population matrix of Nei's genetic distance, Table 5). The result of PCoA is presented in Figure 2. The first two coordinates explain 35% of the total variation. PCoA sorted the populations in three groups, which can be simplified as south-west (Pops. 5, 4, 2, 3 and 7), north (Pops. 1 and 10), and southeast (Pops. 6, 9 and 8). Population 9 was the most distant from the others. The approximate regionalization of the populations is shown in Figure 3. However, in the scale of the Czech Republic, the geographical distance cannot be seen as the main driver of the genetic variation between the tested indigenous populations. AMOVA indicated that 98% of total molecular variance occurs within populations. Only 2%of the variation can be observed on an interpopulation level. It corresponds with the previous findings that a significant part of the genetic variability of forest trees is attributable to the intrapopulation level.

Figure 2

Figure 3