Population structure in Quercus suber L. revealed by nuclear microsatellite markers

Sampling and DNA extraction

Seventeen sampling sites were selected to broadly represent the entire distribution of Quercus suber around the Mediterranean. Leaf material from six sampling sites was collected in situ from natural stands (Portugal: Gerês, Serra da Estrela, Serra da Arrábida, Serra de Monchique, Serra do Buçaco, Serra de Sintra). For the remaining sites (Spain: Catalonia, Haza del Lino; Italy: Apulia, Lazio, Sardinia, Sicily; France: Corsica; Algeria: Forêt de Guerbès; Tunisia: Mekna ; Morocco: Taza, Kenitra), material was obtained from a cork oak provenance trial (FAIR I CT 95 0202) established in 1998 at Herdade Monte da Fava (Santiago do Cacém, Portugal; 8°7′W, 38°00′N) as part of the European Forest Genetic Resources Programme (EUFORGEN; Varela, 2003). Leaf material was obtained from 22 to 30 individuals from each location, to a total of 488 trees sampled, and kept at −80 °C. Leaves were ground manually using liquid nitrogen and genomic DNA was extracted with the DNeasy Plant Mini Kit (Qiagen, Valencia, CA, USA), according to the manufacturer’s protocol. To determine DNA concentration, quality and integrity, extracted DNA was analysed by gel electrophoresis and with an ND-1000 Nanodrop spectrophotometer. A map showing the 17 sampling sites analysed together with the distribution of Q. suber is presented in Fig. 1.

Microsatellite genotyping

A set of 13 simple sequence repeat (SSR) loci was used for analyses. Of these, ten are anonymous nuclear microsatellites (nuSSR): MSQ13, MSQ4 (Dow, Ashley & Howe, 1995); QpZAG9, QpZAG15, QpZAG36, QpZAG46, QpZAG110 (Steinkellner et al., 1997); QrZAG11, QrZAG7, QrZAG20 (Kampfer et al., 1998) and three are expressed sequence tag microsatellites (EST-SSR): QmOST1, QmD12, QmAJ1 (Ueno & Tsumura, 2008). Polymerase chain reaction (PCR) amplifications were performed in a volume of 15 µL reaction mixture containing 0.5 µL of DNA (50–100 ng), 1x PCR buffer (Promega, Madison, WI, USA), 1U Taq polymerase (Promega, Madison, WI, USA), 1.5(2.0) mM MgCl₂ , 0.12 mM dNTPs and 0.3(0.4) µM of each primer. Amplifications were performed under the following general conditions: an initial denaturation step at 94 °C for 5 min followed by 30 cycles consisting of denaturation at 94 °C for 30(60) s, annealing at 50(57) °C for 30 s, extension at 72 °C for 30(60)s and a final extension step at 72 °C for 10 min. PCR products were analysed by capillary electrophoresis using the ABI PRISM 310 automated sequencer. Genotypes were scored and visually controlled using GENEMAPPER v3.7 (Applied Biosystems, Inc., Waltham, MA, USA) The software MICRO-CHECKER v2.2.3 (Van Oosterhout et al., 2004) was used to identify and correct possible genotyping errors.

Genetic diversity analysis

Estimates of inter- and intra-population genetic diversity were calculated. The observed (H_o) and expected (H_e) heterozygosity, the W&C (Weir & Cockerham, 1984) inbreeding coefficient (F_IS ), fixation index (F_ST) and the number of observed alleles (Na) per locus were calculated using the GENEPOP package v 1.1.7 (Rousset, 2008) as implemented in R v 3.6.3 (R Core Team, 2013) using the Rstudio environment (R Studio Team, 2016).

F_IS per locus was estimated considering both the 17 sampling sites separately and the three clusters corresponding to K = 3 (see Results) where one cluster contained Lazio, Sardinia, Sicily, a second cluster contained Corsica, Mekna, Guerbès, Apulia and the third cluster contained all other sites. The frequency of null allele (fna) per locus was calculated using PopGenReport (Adamack & Gruber, 2014). Estimates of null allele frequency were obtained with both the method of Chakraborty et al. (1994), and the method of Brookfield (1996). Per-locus diversity statistics were calculated considering all 17 sampling sites as separate populations. Additionally, F_IS for each population was calculated using GENEPOP and the number of private alleles (Np) and allelic richness (Ar) per population were calculated with PopGenReport v 3.0.4. Pairwise F_ST across all populations was calculated for the 13 loci and the neutral loci data sets (see “Tests of non-neutrality”) with HIERFSTAT (Goudet, 2005), using the Nei 87 (Nei, 1987) genetic distance, and significance was tested with 1,000 bootstraps and a 95% confidence level. A dendrogram of parirwise F_ST for the neutral loci was constructed from 1,000 bootstrap replicates.

Tests of Hardy–Weinberg equilibirum and linkage disequilibrium

Deviations from Hardy–Weinberg Equilibrium (HWE) were estimated for the 13 loci and 17 sampling sites using GENEPOP, with the Monte Carlo based exact test for a null hypothesis of random union of gametes, and p-values obtained from a global test across samples using Fisher’s method (dememorization = 10,000, batches = 20, iterations = 10,000). Linkage disequilibrium (LD) across all pairs of loci and all 17 sampling sites was tested using GENEPOP with the default Monte Carlo based G-test for the null hypothesis of independence among genotypes at different loci, and the log likelihood ratio as test statistic (dememorization = 10,000, batches = 100, iterations = 5,000).

Population structure analysis

Structure patterns among sampled populations of Quercus suber were inferred using the clustering method implemented in STRUCTURE v 2.3.4 (Pritchard, Stephens & Donnelly, 2000). The program STRUCTURE was wrapped under Structure_threader v 1.2.2 (Pina-Martins et al., 2017). An analysis comprising 17 sampling sites and all 13 loci was run with values of K between one and 17, using the admixture model with priors for location information and correlated allele frequencies, a burnin period of 100.000 steps and a post-burning sampling of 500.000 steps, with 20 replicates. The 10 nuSSR loci and the three EST-SSR loci were also analysed separately, using the same parameters as for the analysis of the entire dataset. For each analysis, the best-fitting number of clusters (K) was estimated using the ΔK statistic (Evanno, Regnaut & Goudet, 2005). Finally, an analysis using the 11 “neutral” loci was also performed (see “Tests of non-neutrality”), with the same parameters as the first analysis.

Population structure and inference of K was further tested using the program MavericK v 1.0.4 (Verity & Nichols, 2016). The MavericK model differs from the model used in STRUCTURE as it does not include location information nor correlated allele frequencies. Analyses searched for values of K between 1 and 16 using an admixture model with a free alpha parameter of “1”, and comprised five runs of 10.000 iterations with 10% burnin, with thermodynamic integration set to 20 runs of 10.000 iterations and 20% burnin. MavericK was wrapped under Structure_threader v 1.2.2 (Pina-Martins et al., 2017). Analyses were run on all 13 loci, on the 10 nusSR and three EST-SSR separately, and on the 11 neutral loci data set.

Principal components analysis (PCA) and correspondence analysis (CA) were performed, respectively, using the methods implemented in the R packages pcaMethods v 1.84 (Stacklies et al., 2007) and ADEGENET v 2.1.3 (Jombart, 2008), on both the 13 loci and 11 neutral loci data sets.

Data files and scripts generated for analyses can be found at https://github.com/CoBiG2/Qsuber_mssats.

Tests of non-neutrality

Tests of non-neutrality were performed on the 13 SSR loci dataset using the program BAYESCAN v 2.1 (Foll & Gaggiotti, 2008) which assumes that a multinomial Dirichlet distribution can be used to model gene frequencies in neutrally structured populations (Excoffier, Hofer & Foll, 2009; Lotterhos & Whitlock, 2014). To avoid possible violations of the assumptions in the BAYESCAN model causing an excess of false positive outlier loci, and in order to maximise migration within rather than between clusters (Excoffier, Hofer & Foll, 2009), input data was rearranged to comprise two clusters, following the results from STRUCTURE analyses (see Results), rather than 17 groups corresponding to all sampling sites (e.g., Moore et al., 2014). Analyses were run with prior odds for the neutral model of 10, 100 and 1,000, and 20 pilot runs of 5,000 iterations followed by 10,000 iterations and a burnin length of 50,000, with a thinning interval of 10. To confirm the identification of candidate outlier loci, an additional analysis was run considering the results from MavericK analyses that estimated K = 3 (see Results), with one cluster comprising the sites Sardinia, Lazio, Sicily, a second cluster comprising the sites Corsica, Apulia, Guerbès, Mekna, and a third comprising the remainder of the sites.

Genetic diversity analysis

The 13 SSR loci comprised a total of 145 alleles (Table 1). The number of alleles per locus ranged between 5 (QpD12) and 24 (QpZag110), with an average of 11.2. H_e ranged between 0.187 (QpZag9) and 0.86 (QpZag110), while H_o ranged between 0.142 (QpZag9) and 0.777 (QpZag110). The average frequency of null alleles over all 13 loci, inferred with the method of Chakraborty et al. (1994), was 13%, whereas the value inferred with the method of Brookfield (1996), which takes into account the presence of null homozygotes, was 9%. F_ST per locus ranged between 0.013 (QpZag9) and 0.512 (MSQ13). The mean F_IS  was always positive and ranged between 0.019 (MSQ13) and 0.263 (QmAJ1) when considering the populations separately and between 0.07 (QrZag7) and 0.27 (QmAJ1) when considering K = 3 (Table 1).

Population diversity analysis

Ar ranged between 4.4 and 5.5 among the 17 sites, with the lowest values detected in populations Sicily, Gerês, Sardinia and Corsica, and the largest values found in populations Monchique, Sintra and Guerbès (Table 2). The average observed heterozigosity across 13 loci was lower than expected in all 17 populations. The number of private alleles was zero in five populations (Gerês, Haza del Lino, Kenitra, Apulia, Estrela) and was largest in populations Guerbès and Sintra (six and eight private alleles, respectively) (Table 2). F_IS  varied between 0,046 (Kenitra) and 0,183 (Lazio) (Table 2). Pairwise F_ST estimates among populations using the Nei 87 method on the neutral data set indicate that sites Lazio, Sicily and Sardinia show low differentiation (F_ST < 0.05) but differ from the remainder (F_ST = 0.1–0.24). Corsica also showed differentiation (F_ST = 0.1–0.21) from all other sites (Table S1).

Tests of Hardy–Weinberg equilibirum and linkage disequilibrium

Exact tests of Hardy–Weinberg Equilibrium (HWE) showed an overall departure in all populations (p < 0.05) and in all loci except QrZag7 (Table S2). Tests of linkage disequilibrium (LD) for all pairs of loci showed LD (p < 0.05) at eight out of 78 pairwise comparisons, when considering 17 sample sites (Table S3).

Tests of non-neutrality

Tests of non-neutrality using BAYESCAN and considering two clusters (K = 2) identified two candidate outlier loci, MSQ13 and QpZag110, for all prior odds (10, 100, 1,000) and with a false discovery rate of 5% (FDR = 0.05). The alpha parameter was positive for MSQ13 (alpha = 1.3091 for prior odds = 10) whereas for QpZag110 alpha was negative (alpha = −1.8945 for prior odds = 10). The additional analysis considering three clusters (K = 3) confirmed the identification of the two candidate outlier loci for prior odds of 10 and 100, and FDR =0.05. Plots for the tests of non-neutrality for K = 2 and K = 3 with prior odds of 10 and FDR = 0.05 are shown in Fig. S1.

Structure analysis

In all analyses using STRUCTURE (13 loci, 10 nuSSR loci, three EST-SSR loci, 11 neutral loci), the ΔK statistic indicated an optimal number of clusters of K = 2 (Fig. 2, Fig. S2). In the 11 neutral loci analysis (K = 2, ΔK = 123.52), Lazio, Sardinia and Sicily form a separate cluster, whereas populations Apulia, Mekna Guerbès and Corsica show admixture in all individuals. The admixture in Corsica is balanced, with the smallest parental contribution observed in a single individual corresponding to 39.5% (Fig. 2A).

The optimal number of clusters inferred from the 13 loci, the nuSSR and the EST-SSR datasets, using the thermodynamic integration method implemented on MavericK, was K = 2 (Fig. S2). In the 11 neutral loci analysis, the estimated optimal number of clusters was K = 3. The Q-matrix plot for K = 3 (Fig. 2D) shows that Lazio, Sardinia and Sicily form a distinct cluster without pronounced admixture (except for a few individuals from Lazio and Sicily), and that Corsica forms another cluster with no pronounced admixture. Apulia, Mekna and Guerbès have pronounced admixture and a major genetic contribution from the cluster that includes Corsica. The remainder of the populations form the largest cluster and show admixture that is particularly pronounced in populations Arrábida and Kenitra (Fig. 2D). This same pattern is also seen in the STRUCTURE Q-matrix plot for K = 3 (Fig. 2B), whereas the STRUCTURE Q-matrix plot for K = 4 assigns Apulia, Mekna and Guerbès to the fourth cluster, while keeping Corsica isolated (Fig. 2C). The pairwise F_ST   dendrogram (Fig. 3) confirms the clustering of Lazio, Sicily and Sardinia and the separation of Corsica.

Overall, the PCA using 13 loci does not show a clear separation of clusters, but suggest that Lazio, Sardinia and Sicily form a distinct group, and that Corsica appears separated from the remainder (Fig. S3A), although PC1 and PC2 explain only 6.34% of the variance. The CA of the 13 loci data set separates Lazio, Sardinia and Sicily into one cluster, and isolates Corsica (Fig. S3B). PCA scatterplots of the 11 neutral loci analysis again show a cluster formed by Lazio, Sardinia and Sicily (Fig. 4A), but PC1 and PC2 explain only 6.85% of the variance. The CA scatterplots separate Lazio, Sardinia and Sicily from the remainder of the sites (Fig. 4B), and isolate Corsica and Guerbès.

Quercus suber nuclear SSR data show heterozygosity deficit and presence of null alleles

The utility of several nuclear SSR markers included in our analyses has been demonstrated in studies of Quercus suber diversity (López-Aljorna et al., 2007; Soto, Lorenzo & Gil, 2007; Ramírez-Valiente et al., 2010b), fine-scale structure (Soto, Lorenzo & Gil, 2007; Lorenzo et al., 2009), phenotypic prediction (Ramirez-Valiente et al., 2009) and introgression (Burgarella et al., 2009; López de Heredia, Sanchez & Soto, 2018). Here, we sampled a set of 13 nuclear SSR across 17 locations to infer genetic diversity and structure patterns across the entire natural range of the species. Our analyses of the 13 SSR loci and 17 population data set showed that the observed level of heterozygosity was lower than expected, for all loci, and F_IS was positive for all loci and samples. The observed presence of null alleles, whether true or apparent, is likely to have strongly contributed to the heterozygosity deficit and positive F_IS (Waples, 2015). Nevertheless, the frequency of null alleles inferred with the methods of Brookfield and Chakraborty (0.09 and 0.13, respectively) is in line with values commonly reported in the literature and is unlikely to cause a major bias in downstream population structure analyses (Dakin & Avise, 2004). Another possible explanation for low heterozygosity and positive F_IS is positive assortative mating, of which self-fertilization is an extreme case (Waples, 2015). Selfing is expected to occur to some extent in monoecious species, although Q. suber is known to have dichogamous flower development, which promotes outcrossing (Sobral et al., 2020).

Tests of deviation from HWE showed an overall departure that is not specific to particular loci or populations, and linkage disequilibrium (LD) was observed in a small number of different gene pairs. These observed deviations from HWE are not surprising given the presence of null alleles. Furthermore, the species does not conform to the expectations of the Hardy–Weinberg principle, that assumes discrete generations. Like most oaks, Quercus suber individuals have great longevity and both natural and managed stands are likely to contain individuals of very different ages. Wind pollination is also likely to promote mating among individuals from distant stands with different ages. In the presence of significant age structure, reproductive cycles do not correspond to random mating among adults of the same generation, but of different generations, each with different allele frequencies, meaning that offspring from different reproductive cycles will differ in their allele frequencies as well. Thus, samples containing mixed-age individuals or their progeny are composed, in practice, of different subpopulations, and both a deficiency of heterozygotes and mixture LD at pairs of loci are expected (Waples, 2015). Error caused by over-sampling of plants with particular phenotypes, for example if early- or late-developing plants are inadvertently preferred, may also cause deviations from HWE, even if the populations themselves are at Hardy–Weinberg equilibrium (Waples, 2015). However, an effect of sampling bias alone is unlikely to explain the HWE departure pattern seen here, which was observed across populations. As the deviations from HWE and sporadic LD in SSR loci are likely a result of reproductive traits of cork oak stands, as well as a consequence of the presence of null alleles, the retention of all 13 loci for downstream analyses was considered justified.

Tests of non-neutrality revealed two outlier loci that are possibly under selection. Locus MSQ13, which had the highest F_ST, had a positive value for the alpha parameter on the BAYESCAN analysis, which could indicate diversifying selection. Locus QpZag110, which had the highest observed heterozygosity and highest number of alleles, had a negative value of alpha, which suggests purifying selection. Outlier detection may result in a high false positive rate when the assumptions of the tests are not met (García-Verdugo et al., 2015). This could be the case of the Dirichlet distribution model in BAYESCAN, which assumes neutrally structured populations (Beaumont, 2005) and is not adequate for species that have undergone range expansion or that display hierarchical population structure (Lotterhos & Whitlock, 2014). Indeed, an exploratory analysis considering all 17 sampling sites as distinct populations resulted in the detection of a very high number of non-neutral loci, which contradicted the expectation that most SSR loci are neutral. Detection of the two outlier loci was attained when the original 17 populations were grouped under two or three clusters. As the detection of the same two outlier loci was consistent across different neutrality prior odds and in the analyses of two and three clusters, the two loci were considered non-neutral, and their exclusion from analyses was deemed justified.

Quercus suber is structured into three population clusters along a longitudinal gradient

The nuclear SSR markers chosen in this study proved to be an effective means for detecting population structure in Quercus suber, despite being in much lower number than SNPs . The relative low cost per sample allowed for the inclusion of more samples per location (22 to 30) compared to earlier studies using SNPs, that typically sampled around six individuals per site (Pina-Martins et al., 2019; Vanhove et al., 2021), thus decreasing the probability of biases due to limited sample size.

Analyses of the complete data set, with STRUCTURE and Maverick, and of the neutral data set, with STRUCTURE, show a separation of Quercus suber populations into two clusters (K = 2) that roughly correspond to an East-West separation. This structural pattern is different to those of Pina-Martins et al. (2019) and Vanhove et al. (2021), based on SNP data. In both those studies, the East-West divide involved a different set of sampling sites and the eastern cluster comprised other sites besides Lazio, Sicily and Sardinia. In contrast, our analysis of the neutral data set using Maverick (K = 3) shows a very distinct and homogeneous eastern cluster comprising the central italian population of Lazio and the island populations of Sicily and Sardinia, and a third cluster including Corsica and other sites with pronounced admixture (Apulia, Guerbès and Mekna). The placement of Corsica in a distinct third cluster was also evident in Vanhove et al. (2021).

The hypothesis of Quercus suber being structured in three, rather than two clusters, agrees with the pairwise F_ST calculations, which indicate that populations Lazio, Sicily and Sardinia have significant differentiation from all other populations but not among them, and that population Corsica is distinct from the remainder. Results from PCA and CA analyses also distinguish the group formed by Lazio, Sicily and Sardinia, and separate Corsica from the remainder, often showing a degree of affinity between the latter and populations Apulia, Guerbès and Mekna, which may be a reflection of the admixture pattern for these populations suggested by the Q-plot for K = 3. The identification of population Corsica as belonging to a separate cluster, distinct from the eastern and western clusters of Q. suber populations, was also hypothesized in the work of Vanhove et al. (2021), with a SNPs dataset. In their work, although the inferred number of K, using bayesian approaches, was also K = 2, a shared coancestry matrix assigned individuals from Corsica to a third cluster that is possibly a sub-cluster of the western group. This sub-cluster included populations from southern France (Landes), which are missing from our analyses.

Given the mixed results obtained here with the different methods for estimation of K, it is necessary to consider all evidence in order to ascertain which hypothesis is more likely. Recent studies have shown that the ΔK parameter of Evanno, Regnaut & Goudet (2005) has a tendency to estimate K = 2, and strongly recommend against using ΔK alone (Janes et al., 2017). Furthermore, the TI method, which makes use of several closely related MCMC chains to infer K, has been suggested to have a superior performance in this task compared to other methods, and in fact has been shown to outperform ΔK when estimating K from large data sets (Verity & Nichols, 2016). Considering the evidence from pairwise Fst estimates and from PCA and CA analyses, which point to a more complex scenario than a simple separation of Quercus suber populations into two groups, and the fact that estimates of K = 2 with the ΔK method should be regarded with caution, we posit that a structural pattern of Q. suber containing more than two clusters, as estimated by the TI integration method from the neutral loci dataset, is a valid hypothesis.

The possible existence of more than two genetic clusters in Quercus suber suggests differentiation caused by isolation in refugia, with posterior expansion and secondary contact between isolated populations. Earlier studies using chloroplast and AFLP markers had already suggested the existence of several glacial refugia and a division of Q. suber into three genetic clusters, albeit with different compositions (Lumaret et al., 2005; López de Heredia et al., 2007a; López de Heredia et al., 2007b). The detection of high levels of private alleles may indicate putative refugia, while low genetic diversity indicates newly colonized areas (de Dato et al., 2020). The highest number of private alleles encountered in our data set belong to populations Guerbès and Sintra, while populations Sardinia, Mekna, Monchique and Arrábida also presented an above average number of private alleles. This result is compatible with a scenario of contraction of Q. suber into multiple southern refugia during the last glacial maximum, not only in central southern Mediterranean, but also in southwestern Iberian regions (Lumaret et al., 2005; López de Heredia et al., 2007b; Vessella, Simeone & Schirone, 2015). Thus, the hypothesis that genetic differentiation occurred prior to a recent northwards expansion of Q. suber seems conceivable.

The observation of an incomplete East-West structure in Quercus suber has been interpreted as either the result of a balance between gene flow and local adaptation or the consequence of differential hybridisation, while recent expansion from refugia would only be detected in chloroplast lineages (Pina-Martins et al., 2019). Considering the type of marker used here, essentially neutrally behaving genomic regions which at most segregate passively along with regions under selection, it can be argued that local adaptation is unlikely to account for the observed structure pattern. On the contrary, the more complex scenario involving three genetic clusters was only revealed when loci under selection were excluded from the analysis. This veiling of structure due to selection seems to result from two combined forces. The first is the overall homogenizing effect of balancing selection, where the same alleles are kept in similar frequencies in different populations. The second is directional selection, which promotes identical allele frequencies under similar ecological/climatic circumstances, but differentiates populations in contrasting ecological circumstances, consequently favouring the two cluster scenario. Hybridisation between Q. suber and both Q. ilex and Q. cerris is indeed known to occur in the western and eastern ranges of the species’ distribution, respectively (Lumaret et al., 2005; Magri et al., 2007; Lumaret & Jabbour-Zahab, 2009; López de Heredia, Sanchez & Soto, 2018; López de Heredia et al., 2020). Introgression from related Quercus species may have contributed to the survival of Q. suber through glaciation periods (López de Heredia, Vázquez & Soto, 2017) and may also have influenced the number of private alleles of certain populations. However, it is unclear how post-glacial and contemporary introgression from these two species alone could explain the persistence of three hypothetical Q. suber genetic clusters. Furthermore, the two species known to hybridise with Q. suber are also likely to have been affected by recent contraction-expansion cycles (Guzmán et al., 2015; Bagnoli et al., 2016) and it is in fact possible that hybridisation between species occurred already within shared refugia (Bagnoli et al., 2016). Thus, differentiation of Q. suber populations in multiple glacial refugia prior to a northwards expansion of the species is perhaps the main cause of the tripartite structure observed with nuclear SSR loci, even if introgression from other species during and after range contraction occurred. While its role in differentiation is not completely clear, local adaptation should nevertheless be expected to occur in Q. suber populations, given the latitudinal and altitudinal variation across the distribution of the species. It is also plausible that human activity and migration after the last glacial cycle and to present time may have influenced diversification in Q. suber by facilitating dispersal and gene flow. SSR markers do not always reflect genome-wide genetic diversity in natural populations (Väli et al., 2008), but the identification of distinct genetic clusters and genetically diverse populations based on these 13 nuclear SSR markers should nevertheless provide additional evidence on where new allelic variants and possibly adaptive variation may be encountered.