1. Introduction

The presence of heterozygous individuals in a population is crucial for maintaining genetic diversity, which can positively affect fitness and adaptability to environmental changes (Reed and Frankham, 2003; Kremer et al., 2012). Heterozygosity excess in a population leads to higher genetic diversity, which may improve the population's ability to adapt to changing environmental conditions due to the presence of a reservoir of different alleles that may prove favourable under different selection pressures. However, in essentially all naturally outbreeding species heterozygosity excess can also have negative effects if accompanied by inbreeding depression, in which deleterious recessive alleles are more likely to accumulate and be expressed in closely related individuals. This can lead to reduced fitness (Wright, 1977; Charlesworth and Charlesworth, 1987; Thornhill, 1993; Frankham, 1995a; Falconer and Mackay, 1996), lower reproductive success, a general population decline (Balloux et al., 2004) and an increased risk of extinction (Frankel and Soulé, 1981; Frankham, 1995b; Newman and Pilson, 1997).

Inbreeding usually leads to a reduction in the proportion of heterozygous individuals in a population, as it increases the probability that two alleles in an individual come from the same ancestor and are thus identical by descent, which is referred to as the inbreeding coefficient (F) (Wright, 1969; Crow and Kimura, 1970). Therefore, the inbreeding coefficient of a population (F_IS) is a useful parameter for analyzing various aspects of the evolution of plant mating systems, e.g., the relationship between inbreeding history and certain traits, such as self-fertilization capacity or probability (Schultz and Willis, 1995). F_IS is highly dependent on the rate of asexual reproduction within a population (Stoeckel et al., 2006). For instance, asexuality limits the segregation of alleles, conserves ancestral heterozygosity through generations, and increases heterozygosity, since alleles of the same gene locus may independently accumulate mutations over generations, in a phenomenon also known as the Meselson effect (Judson and Normark, 1996; Stoeckel and Masson, 2014). The mutations only accumulate if they are (almost) neutral and there is no overdominance (i.e., there is an additive effect of the mutations on the phenotype). In addition, progressive or balanced selection can promote the overlapping of heterogeneous strains/clones (Haigh and Smith, 1974; Pavlidis and Alachiotis, 2017).

In general, polyploid individuals and populations retain a higher degree of heterozygosity and exhibit less inbreeding depression than their diploid ancestors and can therefore tolerate a higher degree of self-fertilization (Soltis and Soltis, 2000). Polyploidy is frequent in angiosperms, occurring in more than 50% of species (Weiss-Schneeweiss et al., 2013). Indeed, it is believed that environmental fluctuations can lead to adaptive responses via the formation of polyploids (Van de Peer et al., 2017), especially in connection with stress (Van de Peer et al., 2021). Polyploidy can improve tolerance to abiotic and biotic stress factors and boost disease resistance, which can have positive effects on plant growth and net production (Greer et al., 2018; Blonder et al., 2020; Tossi et al., 2022). The high levels of genetic diversity in polyploids can help individuals and populations to colonize regions and persist under strong climate variation (Dynesius and Jansson, 2000).

Populus tremuloides Michx. (quaking aspen) is one of the most widely distributed and ecologically important tree species in the Northern Hemisphere (Wang et al., 2016). The species distribution ranges from Alaska through Canada and the United States to central Mexico, covering a wide range of ecosystems and elevations (Little, 1971; Perala, 1990). Quaking aspen is an early successional tree that is adapted to natural disturbances (e.g., fire and disease); it is highly tolerant of environmental stress and is capable of colonizing new areas and rejuvenating existing stands (Swift and Ran, 2012; Goessen et al., 2022). These capacities are due to the reproductive flexibility of the species, which can reproduce sexually (through seeds and pollen suitable for wide-ranging wind dispersal) and asexually (through root suckers) (Barnes, 1966; Mitton and Grant, 1996), and to ploidy-level variation (Mock et al., 2012; Goessen et al., 2022). The quaking aspen, together with the European aspen, has the highest level of intraspecific genetic diversity ever recorded in a plant species (Cole, 2005; Callahan et al., 2013; Wang et al., 2016), which also contributes to its great adaptability.

However, the genetic diversity and structure of Mexican populations of Populus tremuloides, including the genetic signatures of their adaptation to a variety of environments, remain largely uncharacterized (Ouborg et al., 2010; Callahan et al., 2013). Previous studies used AFLP markers to identify significant large-scale spatial genetic structure in P. tremuloides (Quiñones-Pérez et al., 2014) and found that P. tremuloides had lower AFLP diversity than other Mexican tree species, as well as significantly positive associations between neighbor tree species diversity and genetic diversity (Simental-Rodríguez et al., 2014). In addition, one study reported that almost all P. tremuloides individuals in each of seven forest sample plots in the Mexican Sierra Madre Occidental were genetically different concluding that vegetative reproduction probably only plays a secondary role in quaking aspen in this region (Wehenkel et al., 2014). These findings are partially in accordance with those reported by Goessen et al. (2022), who investigated the entire Mexican genetic cluster by using single nucleotide polymorphism (SNP) markers. These researchers found that: (i) a large number of clonemates are located in Mexican aspen stands and (ii) several SNPs are under selection and impacted by temperature and precipitation variation across the four main genetic clusters in northeastern North America, northwestern North America, the western United States and Mexico.

The aim of this study was to analyze how inbreeding coefficient and ploidy level are associated with clonal richness, climate and soil traits in 91 marginal to small, isolated populations of Populus tremuloides throughout its entire distribution in Mexico (i.e., Baja California, Sierra Madre Oriental, Sierra Madre Occidental, and the Trans-Mexican Volcanic Belt), by using SNP markers derived from genotype by sequencing (GBS). We hypothesized that: (i) asexual reproduction will increase the expected allelic diversity within populations and will lead to more negative F_IS, i.e., extreme heterozygote excess (Stoeckel and Masson, 2014), and (ii) progressive selection against deleterious recessive alleles will result in less negative F_IS values in adaptive loci (Mitton (1989) in Stoeckel et al. (2006)). We also hypothesized that the rate of asexual propagation affects ploidy and that higher ploidy levels tend to occur in more extreme environmental conditions (Dynesius and Jansson, 2000; Goessen et al., 2022).

2. Materials and methods 2.1. Study area and sampling

Leaf samples were collected from 809 spatially georeferenced Populus tremuloides trees in 91 natural populations in Mexico. These populations were distributed in the Sierra de San Pedro Mártir (SSPM; Baja California), Sierra Madre Occidental (SMO; Durango, Chihuahua and Sonora), Sierra Madre Oriental (SMOr; Coahuila, Nuevo León and Tamaulipas) and the Trans-Mexican Volcanic Belt (Hidalgo, Michoacán and Queretaro, in central Mexico) (Fig. 1). Detailed information about the location of the sampled stands is provided in Table S1. The population cover (in hectares) was determined using Google Earth Pro v.6 (http://www.google.com/earth/index.html [Accessed 10 November 2023]) (see Table S2 for further details).

2.2. DNA extraction

DNA was extracted from dried leaves obtained in the field. The Nucleospin 96 Plant II kit (Macherey–Nagel, Bethlehem, PA) was used following the manufacturer's protocol for the centrifugation process, with modifications regarding the cell lysis step (PL2 buffer was heated at 65 ℃ for 2 h instead of 30 min). DNA concentrations were adjusted to 10 ng/μl before library preparation and randomization of samples.

2.3. Library preparation and sequencing

Libraries were prepared at the «Plateforme d'analyses génomiques, Institut de Biologie Intégrative et des Systèmes» (IBIS, Université Laval, Québec, Canada). Libraries were sequenced on an Illumina Novaseq 6000 S4 (1 lane) sequencing system (Centre d'expertise et de services Génome Québec in Montréal, Canada), and 150-bp paired-end were obtained in the FastQ format (as described in Poland et al. (2012), with some modifications following Goessen et al. (2022).

2.4. Assembly, quality filtering and SNP calling

DNA was digested using the restriction enzymes PstI, NsiI and MspI. Stacks v.2.4 was used to identify loci in the data and call genotypes (https://github.com/enormandeau; Catchen et al., 2013; Rochette et al., 2019). Cutadapt v.1.8.1 was used to remove adapters from raw sequences (Martin, 2011). The process_radtags module in R v.4.2.1 software (R Core Team, 2022) was used to demultiplex the libraries and perform quality trimming on reads of length 125 bp.

Sequenced reads were aligned to the P. tremuloides reference genome (v.1.1; "Potrs01b, " ~480 Mbp of sequence) available from the Populus Genome Integrative Explorer website (http://www.popgenie.org/; Sjödin et al., 2009; Sundell et al., 2015).

Chloroplast, mitochondrial and contaminant scaffolds were removed from the original reference genome. The unwanted scaffolds consisted of 214 sequences and 211, 552 nucleotides, which accounted for 0.06% of the genome. Reads were aligned to the reference genome on the basis of sequence similarity, determined using the Burrows-Wheeler alignment tool, BWA v.0.7.17 (Li and Durbin, 2009). SNPs were then called. The gstacks module was first used to assemble and merge paired-end contigs and call variant sites in the populations and genotypes in each sample (with the argument --max-clipped 0.1). The populations module was then run to filter the data (-p 200, -r 0.01, --vcf, --hwe, --smooth, --fasta-loci) and export it to variant call format (VCF).

For all subsequent analyses, VCFtools v.0.1.14 (Danecek et al., 2011) was used to generate a "final" SNP dataset, from which indels and problematic individuals with > 20% missing data were removed. In addition, SNPs that did not meet the following criteria were filtered out: (i) minimum allele count (MAC) of 2, (ii) minimal depth (minDP) of 10, and (iii) missing rate of higher than 15%.

2.5. Ploidy levels in samples

The FastPloidy script was used to detect ploidy variation in individuals resulting from filtering. FastPloidy, developed by Goessen et al. (2022), is available at https://github.com/RGoess/Ploidy_detection/blob/main/FastPloidy.R. The correct identification of ploidy levels by the FastPloidy script has already been verified (see also fig. 3 in Goessen et al., 2022) using flow cytometry and microsatellites on 105 individuals provided by Mock et al., (2012). These individuals, along with our material, were part of the dataset used by Goessen et al. (2022).

FastPloidy reads of VCF in R v.4.2.1 (R Core Team, 2022) allows: (i) extracting the allele depth for reference and alternative alleles (minimum depth of 18); (ii) dividing the depth of the reference allele by the total depth to obtain the allelic ratio of the reference allele and only retaining heterozygous SNPs (allelic ratio within 0.1–0.9); and (iii) assigning each SNP to a class depending on the value of the allelic ratio (A: 0.273–0.393, B: 0.44–0.56, C: 0.607–0.727, 'Other' would refer to a value anywhere outside the range). Triploids should have a higher allelic ratio, as most SNPs should be found within class A and C, while diploids should include most SNPs in class B, and thus, have a lower ratio. A comparison of FastPloidy and gbs2ploidy was performed to evaluate the efficiency of each (Gompert and Mock, 2017), and a slight advantage in accuracy for assigning ploidy by the FastPloidy package was observed (for more details see Goessen et al., 2022). The population cover (ha) and clonal richness, as well as the relative frequency of diploid and triploid individuals per population were then calculated in order to test the univariate association between polyploidy and bioclimatic and edaphic variables. Since there was probably only one tetraploid individual in the study (see also in Goessen et al., 2022), we merged it with the triploids for all analyses.

2.6. Assessment of clonal richness

The pairwise Jaccard similarity index (across all SNPs) (Jaccard, 1908) was used to assign an individual to a clone (genotype), and a similarity matrix was constructed using the vcf2Jaccard.py script developed by Rowe (2019) (available at: https://github.com/carol-rowe666/vcf2Jaccard.git). Consequently, samples from the same population with similarities ≥ 0.99 were assigned as belonging to the same clone (Blonder et al., 2021). In subsequent analyses to determine clonal richness within each study site (genotype), the clonal values were multiplied by a correction term (N∕(N1)) (Gregorius, 1978), where N refers to sample size per population.

2.7. Linkage disequilibrium

Natural selection, which supports adaptation to local environmental conditions, may increase the overall linkage disequilibrium (LD) caused by differences in allele frequencies between populations when alleles at different loci are favoured (Wright, 1940; Ohta, 1982a, 1982b; Agapow and Burt, 2001; Slatkin, 2008). LD analyses were computed with Tomahawk v.0.7.0 (Klarqvist, 2019; available from https://github.com/mklarqvist/tomahawk). The LD metrics (r² and D′) were calculated across all pairwise combinations of variant sites by using the following criteria: (i) P (Fisher's exact test/Chi-squared cutoff p-value) = 0.0001 and, (ii) r² (Pearson's R-squared minimum cut-off value = 0.0001).

2.8. Influence of putative adaptation on H_o, H_e and F_IS

To examine the influence of selection on F_IS in the natural populations of Populus tremuloides, putative outlier loci in the populations under study were identified using the multivariate method implemented in the pcadapt v.4.1.0 package (Luu et al., 2017) and Latent Factor Mixed Models 2.0 (LFMM2) (Caye et al., 2019) from the LEA package (Gain and François, 2021) of R v.4.2.1 (R Core Team, 2022).

Simulations showed that the pcadapt package compares favourably to other genome scanning software (BayeScan, hapflk, OutFLANK, sNMF), particularly in the presence of admixed individuals (data not shown). This method involves: (i) identification of outlier loci in admixed or continuous populations, (ii) determination of the population structure by principal component analysis (PCA) (rather than classifying individuals in admixed or continuous populations), and (iii) identification of SNPs under putative selection as those that are strongly correlated with population structure. The following criteria were applied: (i) exclusion of loci with global MAF < 0.05, (ii) an optimal K = 2 (according to a scree plot (Jackson, 1993; Fig. S1) and (iii) identification of outliers for different individuals by the significance of Mahalanobis distances for each SNP, using a Benjamini-Hochberg (HB) correction for a false discovery rate (FDR) of p < 0.05 (Thissen et al., 2002).

Since LFMM2 does not allow missing data, imputation was carried out with the impute (.) command of the R package "LEA" (Frichot and François, 2015). This imputation method replaces missing genotypes with predicted genotypes based on the population structure inference made in the snmf function of the package (Gain and François, 2021). We decided to run LFMM2 with two latent factors (according to a scree plot Fig. S1). The HB correction was used to adjust false-positive rates. SNPs with a corrected p value < 0.05 were retained.

The results of pcadapt and LFMM2 were processed using VCFtools to generate two VCF files, one for the putatively neutral SNPs and the other one for the putatively outlier SNPs. Only SNPs that were found to be putatively outlier SNPs by both methods and were also not identified as deleterious SNPs (more details in section 2.9.) were used further. Each dataset was analyzed separately to calculate values of observed heterozygosity (H_o), expected heterozygosity (H_e) and F_IS in the R package hierfstat (Goudet and Jombart, 2022), and to test univariate associations with bioclimatic and edaphic variables, population cover and clonal richness. Given that F_IS = 1 – H_o/H_e (Wright, 1969), F_IS and H_e are positively correlated when H_o is constant.

2.9. Gene annotation and association of deleterious mutations with heterozygous excess and clonal reproduction

Gene annotation was performed to categorize non-synonymous, synonymous and deleterious SNPs based on their location and impact using SnpEff v.5.2c and Sequences Ontology (SO) (http://www.sequenceontology.org/) for accurately annotating genetic sequences. The deleterious SNPs included the SO categories, in particular: (i) missense mutations with moderate impact, (ii) SNPs that cause the loss or disruption of the start codon of a gene (start_lost), (iii) SNPs that lead to the formation of a premature stop codon (stop_gained), (iv) SNPs that have both stop-gained and splice region effects (stop_gained & splice_region_variants), (v) SNPs that cause the loss of a stop codon (stop_lost), and finally (vi) SNPs that have both stop-loss and splice region effects (stop_lost & splice_region_variants). The SNP categories (ii) to (vi) are classified as high impact (Cingolani et al., 2012). To construct the input database for this gene annotation, we used the reference genome of Populus tremuloides, available at https://plantgenie.org/FTP?dir=Data%2FPlantGenIE%2FPopulus_tremuloides%2Fv1.1.

By comparing the accumulation of heterozygous deleterious SNPs due to heterozygous excess and clonal reproduction, these findings helped to decipher not only the causes but also the consequences of excessive heterozygosity.

2.10. Estimation of edaphic variables

For determination of 25 soil variables, four subsamples of 500 g of soil at a depth of 0–15 cm were randomly collected in each population. The subsamples were then mixed to produce a sample of 2 kg for each population. The techniques used to estimate the variables are described below. The concentrations of K⁺ (ppm), Ca²⁺ (ppm), Mg²⁺ (ppm), Na⁺ (ppm), Cu²⁺ (ppm), Fe³⁺ (ppm), Mn⁺² (ppm), Zn²⁺ (ppm), texture (relative proportion of sand, silt, and clay) and pH (CaCl₂ 0.01 M) in each soil sample were determined using the methods described by Castellanos et al. (1999). The concentration of P (kg/ha) was determined using the method detailed by Olsen et al. (1954). The relative organic matter content (%OM) was obtained by the León and Aguilar (1987) method, and that of nitrates (NO₃-) (kg/ha) by the Baker (1967) method. Electrical conductivity (CE) (dS/m) was determined according to Vázquez-Alarcón and Aguilar-Noh (2020). The hydraulic conductivity (HC) (cm/h) and saturation percentage (Sat%) were determined as described by Mualem (1976) and Herbert (1992), respectively. Finally, the relative proportions (%) of H⁺, Ca²⁺, Mg²⁺, K⁺, Na⁺ and other bases (%OB) were obtained by calculating the cation exchange capacity (meq 100 g soil) (CEC) using the ammonium acetate method (pH 8.5).

2.11. Estimation of bioclimatic variables

The Rehfeldt's (2006) climate model was used to estimate 22 bioclimatic variables in each sampling site. The model is based on the Thin Plate Spline (TPS) method (Hutchinson, 1991) and has previously been used for climate-change research in Mexico (Sáenz-Romero et al., 2010). The estimations of this model are based on the reports from more than 1700 Mexican weather stations for the period 1961–1990, estimating the standardized monthly mean values of minimum and maximum temperatures and precipitation (Table S1). The geographic coordinates (latitude, longitude and elevation) of each site were uploaded to a public database (available at http://forest.moscowfsl.wsu.edu/climate/[accessed 22 October 2023]) to obtain point estimates of the bioclimatic variables.

2.12. Statistical analysis

The non-parametric post-hoc Kruskal Wallis test (Tukey and Kramer method) (Kruskal and Wallis, 1952; Sachs, 2013) was used to test for any statistically significant differences in the medians of H_o, H_e and F_IS between SNP datasets (Table S3). The tests were performed in R v.4.2.1 (R Core Team, 2022), using the PMCMR package (Pohlert, 2014).

Once the population cover, relative frequency of diploids and triploids, clonal richness, number of heterozygous deleterious SNPs, bioclimatic and edaphic variables and H_o, H_e and F_IS were calculated for each sampling site, the degrees of association (measured by the Spearman's r_s) and the corresponding p-values were determined with the fuzzySim package (Barbosa, 2015).

To test the potential multivariate associations between both F_IS and ploidy level and the bioclimatic and edaphic variables, population cover and clonal richness, only the non-collinear variables and with an absolute value of the nonparametric Spearman correlation coefficient lower than 0.7 were selected for modelling eight machine learning algorithms (MLM), all including cross-validation (CV). The modelling methods used were: (i) Random Forest (rf), (ii) Regularized Random Forest (RRF), (iii) Tree Models from Genetic Algorithms (evtree), (iv) Multi-Layer Perceptron (mlpWeightDecay), (v) Bayesian Regularized Neural Networks (brnn), (vi) Model Averaged Neural Network (avNNet), (vii) Linear Regression (lm), and (viii) Neural Network (nnet). These methods were applied using the caret package together with the train function (Venables and Ripley, 1999; Williams et al., 2015; http://topepo.github.io/caret/index.html) in R software. The models were tested by 10-fold cross-validation, with 30 repetitions. The goodness-of-fit of the regression models was determined by the mean square error (MSE), root mean square error (RMSE) and the (pseudo) coefficient of determination (R²). Analysis was performed using the F_IS values separately for neutral, outlier and all SNPs.

Bioclimatic and edaphic variables, population cover, and clonal richness were tested regarding their conditional importance (i.e., considering all other variables studied) in predicting F_IS (all SNPs, outlier SNPs and neutral SNPs) and ploidy levels, considering the following criteria: (i) the multivariate association between F_IS and ploidy levels and the bioclimatic and edaphic variables, population cover and clonal richness; and (ii) the multivariate association between F_IS and ploidy levels and the bioclimatic and edaphic variables and the population cover. This procedure was performed using the cforest algorithm in the random forest implementation in the party package (Strobl et al., 2009) in R v.4.2.1 (R Core Team, 2022). The packages "ggpubr" (Kassambara, 2023) and "ggplot 2" (Wickham, 2016) were used in R v.4.2.1 to display associations graphically.

3. Results 3.1. SNP filtering, ploidy levels, clonal richness and pattern of linkage disequilibrium

After sequencing, SNP calling and filtering for quality using Stacks v.2.4, and script available in Stacks workflow, the raw VCF file contained 1, 556, 064 singleton SNPs with an average depth of 16.7 per SNP and 932 individuals. After SNP filtering, 36, 810 SNPs and 809 individuals from 91 populations were retained for subsequent analyses. The histogram of H_o is shown in Fig. 2.

FastPloidy script identified 718 diploids and 91 triploids. A total of 293 unique genotypes were identified from among the 809 individuals (Table S2).

The distribution pattern of linkage disequilibrium (LD) (r²) values between pairs of SNPs (in a total of 36, 810 SNPs) ranged from 0.0 to 1.0. The pairwise mean values of r² decreased rapidly with increasing physical distance ≥10 kbp. The overall mean LD value across the genome, 0.46, supports the occurrence of selection (Fig. S2).

3.2. Influence of putative adaptation on H_o, H_e and F_IS

Pcadapt and LFMM2 analyses identified 183 potential adaptive outlier SNPs (optimal K = 2 according to a skew plot; Fig. S1), excluding 37 of those outliers that were found to be deleterious (Fig. S4). The potential adaptive outlier SNPs without the detected deleterious SNPs (146 SNPs) were treated separately from the remaining (neutral) 36, 627 SNPs for estimating genetic diversity and performing correlations with the environment and with individual clonality (Tables S5, S6 and S7). Significance of the associations between SNPs and environmental factors across the genome are shown in a Manhattan plot (Fig. S3).

The non-parametric post-hoc Kruskal Wallis test (Tukey and Kramer method) revealed significant differences between the medians of the 36, 627 neutral SNPs and the 146 outlier SNPs for H_o (p < 0.0000001), H_e (p < 0.0000001), except for F_IS (p = 0.54). All H_e values using the outlier SNPs were much smaller than those using the neutral SNPs (Fig. 3).

Spearman's correlation showed that although H_o and H_e were not significantly correlated for neutral SNPs (r_s = 0.02, p = 0.83), a strong positive correlation was observed between H_o and H_e for the 146 outlier SNPs (r_s = 0.89, p < 0.0000001). There was also a positive, significant association between F_IS (both for all neutral SNPs and the 146 outlier SNPs) and clonal richness in the 91 populations of P. tremuloides (r_s = 0.88 and 0.79, p < 0.0000001). By contrast, ploidy level was negatively correlated (not significantly) with clonal richness (r_s = −0.10, p = 0.36). When all neutral SNPs were considered, clonal richness was negatively correlated with H_o (r_s = −0.23, p = 0.03) but positively correlated with H_e (r_s = 0.87, p < 0.000001). When we analysed only the outlier SNPs, clonal richness again was negatively correlated with H_o (r_s = −0.29, p = 0.01) and positively correlated with H_e, however, this positive correlation was not significant (r_s = 0.05, p = 0.65).

3.3. Associations of deleterious mutations with heterozygous excess and clonal reproduction

Out of the total 36, 810 filtered SNPs, 28, 998 were detected in gene regions, 3926 in intron regions, and 3886 in intergenic regions according to the gene annotation. In total, 5697 SNPs were classified as deleterious mutations (Fig. S4).

As the total heterozygous SNPs per individual increased, the total heterozygous deleterious SNPs per individual also increased (r_s = 0.88, p < 0.00000001), but the relative proportion of total heterozygous SNPs per individual significantly decreased (r_s = −0.23, p < 0.00000001). At the population level (with H_o per population remaining almost constant), there were no significant correlations between the relative proportion of total heterozygous deleterious SNPs and either clonal richness or F_IS. At the population level, the relative proportion of heterozygous deleterious SNPs in the total number of heterozygotic SNPs in clones varied from 0.11 to 0.14. This variation was slightly greater than in populations with higher clonal richness (Fig. 4).

3.4. Associations between F_IS, clonal richness, polyploidy and selected environmental variables

By selecting variables classified as non-collinear, according to the absolute value of the nonparametric Spearman correlation coefficient (ǀr_sǀ) < 0.7, population cover (Table S1), clonal richness (Table S2), five bioclimatic variables (Table S3) and eighteen edaphic variables (Table S4) were retained.

Using all filtered SNPs, Spearman's correlation showed statistically significant associations of H_o with relative proportion of other bases in cation exchange capacity (%) in the soil (r_s = −0.45, p = 0.000008), mean annual precipitation (mm) (r_s = −0.45, p = 0.000008), mean temperature in the warmest month (℃) (r_s = +0.41, p = 0.0001), degree-days below 0 ℃ (DD0) (r_s = +0.29, p = 0.0046) and relative frequency of triploids (r_s = +0.60, p < 0.000000001) at the population level (91 populations). Using the outlier SNPs, significant associations were also found between H_o and the above-mentioned variables (p < 0.05), except for H_o vs. DD0 (p = 0.44). There were also significant correlations of clonal richness with phosphorus concentration in the soil (ppm) (r_s = −0.35, p = 0.0006), clay proportion in the soil (%) (r_s = +0.31, p = 0.0032) and winter precipitation (mm) (r_s = −0.27, p = 0.0087) (Fig. 5). The relative frequency of triploids was significantly associated with mean temperature in the warmest month (r_s = +0.29, p = 0.0046), DD0 (r_s = +0.27, p = 0.0098) (Fig. 5) and winter precipitation (r_s = +0.26, p = 0.01). However, statistical associations of H_o, clonal richness or the relative proportion of triploids with geographic longitude or latitude were not observed (Fig. 1).

Using the non-collinear bioclimatic and edaphic variables, population cover and clonal richness, the Regularized Random Forest (RRF) algorithm proved to be the best for estimating the spatial pattern of F_IS (from the dataset with neutral SNPs) in Populus tremuloides populations, producing a pseudo R² of 0.95 and a RMSE of 0.08 (Table 1). The results obtained with the cforest random forest algorithm showed that clonal richness was strongly associated with F_IS (with 70.6% of all included variables), population cover (with 11.4%) and phosphorus concentration in the soil (with 7.7%). Using the 146 outlier SNPs, RRF was the best for estimating F_IS, with a pseudo R² of 0.65 and RMSE of 0.23) (Table 2). In this case, the cforest results obtained indicated that the variables most closely associated with F_IS were mean temperature of the warmest month (with 34.6%), relative proportion of Mg in the cation exchange capacity (11.1%), and Julian date of the first freezing date of autumn (8.4%). However, including the 37 deleterious SNPs together with the outlier SNPs resulted in an R² of 0.76 and RMSE of 0.19. Here, the cforest results were very similar to the findings using the neutral SNPs.

When using the bioclimatic and soil variables and population cover (i.e., excluding clonal richness), RRF proved to be the best model for estimating F_IS (from the dataset with neutral SNPs) in the populations, producing a pseudo R² of 0.35 and RMSE of 0.30 (Table S7). In this case, the cforest results obtained showed that the variables most closely associated with F_IS were population cover (40.6%), soil phosphorus content (18.5%), and winter precipitation (12.7%). Using the 146 outlier SNPs, RRF was found to be the best model for estimating F_IS, with a pseudo R² of 0.28 and RMSE of 0.34) (Table S8). Here, the cforest results obtained indicated that the variables most closely correlated with F_IS were mean temperature of the warmest month (with 27.6%), relative proportion of Mg in the cation exchange capacity (14.6%), and Julian date of the first freezing date of autumn (8.4%). Including the deleterious with the neutral outlier SNPs together gave an R² of 0.38 and RMSE of 0.30.

The Tree Models algorithm in Genetic Algorithms (evtree) was the best model for estimating the associations between ploidy level and bioclimatic and edaphic variables, population cover and clonal richness, yielding a pseudo R² of 0.26 and RMSE of 0.24 (Table 3). Using the outlier SNPs, Regularized Random Forest (RRF) had the best model for estimating F_IS, with a pseudo R² of 0.30 and RMSE of 0.33. The results obtained with cforest showed that the variables most closely (and positively) associated with ploidy level were winter precipitation (39.6%), soil phosphorus content (18.4%) and relative proportion of Mg²⁺ in CEC (16.8%).

4. Discussion

This study found that the values of the inbreeding coefficient (F_IS) in Mexican populations of Populus tremuloides were mostly negative (Table S2). Comparing F_IS levels between neutral and outlier SNPs in 91 populations of P. tremuloides showed that the median of putative outlier SNPs was not significantly larger (Fig. 3). It can be assumed that: (i) the negative F_IS can at first be attributed to a bottleneck effect in the populations studied, with a single allele retained for most loci, while the less common alleles represent more recent mutations (Cole, 2005); (ii) later processes of genetic drift led to accumulation of almost always recessive alleles derived from mutation and absence of recombination in small populations of trees with no gene flow between individuals of the same population or between populations; and (iii) natural selection forces subsequently acted to eliminate these recessive alleles and preserve rare beneficial alleles (Fig. 3), resulting in highly heterozygous populations and maintaining negative F_IS values (see the Meselson effect; Judson and Normark, 1996; Stoeckel and Masson, 2014). Genetic drift processes (claim (ii)) are supported by the strong positive correlation found between the total number of heterozygous SNPs per individual and heterozygous deleterious SNPs per individual. The negative relationship between the relative proportion of heterozygous deleterious SNPs per individual and the total number of heterozygous SNPs per individual (Fig. 4) are consistent with the action of natural selection forces (claim (iii)). Both conclusions correspond with the findings of low variation of the minor relative proportion of heterozygous deleterious SNPs independent of F_IS and clonal richness at the population level. These two results are in agreement if it is assumed that only a certain upper proportion of deleterious SNPs ensures the fitness of the populations and that only the fittest individuals could be found.

After constructing a population genetics model based on genotypic states in a finite population (less than 400 individuals) with mutations, Stoeckel and Masson (2014) concluded that reproduction by asexuality increased the expected allelic diversity within populations, as indicated by a more negative F_IS index, compared to similar, fully sexual populations. Asexuality fixes heritability and genetic drift at the genotypic level and preserves ancestral genetic states against drift, ultimately reducing allelic identities within individuals.

There is a lack of a significant correlation between H_o and H_e for neutral SNPs, in contrast to the strong positive and significant correlation between H_o and H_e for outlier SNPs. Moreover, H_o was nearly always larger than H_e for neutral and also for outlier SNPs, and H_e for outlier SNPs was much smaller than H_e for neutral SNPs (Fig. 3). The F_IS association of the outlier SNPs with the environmental variables studied was higher than those of the neutral SNPs with these variables. These findings suggest that these outlier SNPs found in Mexican populations of Populus tremuloides may be primarily subject to balancing selection, in which heterozygous genotypes have a fitness advantage over homozygous genotypes, resulting in multiple alleles being maintained in the population due to their adaptive advantages in different environments (Delph and Kelly, 2014). This contrasts with a whole-genome resequencing study from Canada and the USA by Wang et al. (2016), in which purifying and positive selection strongly affected the patterns of nucleotide polymorphism at linked neutral sites in P. tremuloides, P. tremula and P. trichocarpa.

By analyzing polymorphic enzyme loci, Jelinski and Cheliak (1992) also observed a substantial excess of heterozygotes in Populus tremuloides (F_IS of −0.10), attributing this result to the rather dry climate in the western range of the species, which limits seed germination, and results in predominantly asexual reproduction through the formation of adventitious shoots on the roots. Cheliak and Dancik (1982) also obtained a F_IS of −0.24 for seven locations in Alberta, Canada, indicating that asexual propagation could increase genetic diversity by facilitating the accumulation of mutations in offspring. By contrast, Latutrie et al. (2016) detected recent bottlenecks with an excess of heterozygotes in six populations distributed in the northwestern range of aspen distribution between Manitoba and Alaska.

Negative F_IS values have also recently been reported in other species of the genus Populus, e.g., Populus yunnanensis (F_IS = −0.65; Zhou et al., 2020) and Populus laurifolia (F_IS = −0.11; Wiehle et al., 2016), according to results based on the presence of high-frequency gene flow. On the other hand, after performing a mutation drift equilibrium test, Wu et al. (2020) reported significant heterozygosity excess (F_IS = −0.50) in Populus wulianensis, probably because of recently experienced bottleneck effects in the study populations. Similar results were observed in other forest tree species, e.g., Prunus avium L. (Stoeckel et al., 2006) and Tilia cordata Mill. (Erichsen et al., 2019).

While polyploidy is frequently observed in Populus tremuloides (Zhu et al., 1998; Mock et al., 2008), we reject the hypothesis that the ploidy level affected F_IS (r_s = −0.07, p = 0.53), which contradicts the claims made by Krieger and Keller (1998) and Ridout (2000). This level of ploidy has no significant effect on the ratio between H_o and H_e, which is crucial for F_IS (Wright, 1969). Although triploids are thought to reproduce almost exclusively by asexual methods (Mable, 2004; Baldwin and Husband, 2013; Chinone et al., 2014) clonal richness was not significantly lower (r_s = −0.10, p = 0.36).

A large number of clones was detected in this study (36% of the individuals). Similar results were reported by Goessen et al. (2022) for the entire Mexican aspen cluster (not at the population level) using a large part of our dataset. Our results provide evidence that asexual (clonal) reproduction is an indicator of extremely negative F_IS occurring within the study populations based on the strong association between these variables (r_s = 0.79–0.88, p < 0.0000001), as also discussed by Cheliak and Dancik (1982) and Jelinski and Cheliak (1992). Therefore, it can be inferred that clonal growth maintains or even increases the levels of heterozygosity and adaptive capacity by mutation over generations (Judson and Normark, 1996; Welch and Meselson, 2000; Delmotte et al., 2002).

In addition, we detected a clear influence of bioclimatic and edaphic variables on the prediction of F_IS and ploidy level of putatively neutral and outlier SNPs (Table 1, Table 2, S7 and S8). We can confirm that F_IS (from the dataset with all SNPs) was higher in populations with lower soil phosphorus content and less winter precipitation, but that triploids occurred in locations with higher winter precipitation. Therefore, our study, which only included 23 environmental variables, did not support the general expectation that polyploids should display a greater capacity than diploids to adapt to a wide range of conditions (Černá and Münzbergová, 2015). Our results contradict our hypothesis and the findings reported by Goessen et al. (2022), who analyzed Mexican populations as a single cluster. Thus, our results seem to be more consistent with those reported by Benson and Einspahr (1967) and Greer et al. (2018) as triploid individuals have more resource-oriented life histories, and may have a higher risk of induced mortality in drier environments (Dixon and De Wald, 2015; Blonder et al., 2021, 2022). Hence, the wide geographic range of the quaking aspen (Wang et al., 2016) implies that natural populations grow in diverse environmental conditions. This, in turn, leads to local adaptation (Savolainen et al., 2007; Neale and Kremer, 2011; Aitken and Whitlock, 2013). Thus, genotypes that originate from a particular habitat exhibit greater fitness in that specific habitat than genotypes from other habitats (Kawecki and Ebert, 2004), which indicates strong local adaptation. This highlights that adaptation processes rely on the emergence of advantageous mutations, which are then either fixed or maintained at an intermediate frequency through natural selection (Peischl and Kirkpatrick, 2012). This has also been observed in P. tremula, in which Wang et al. (2018) detected a locus (PtFT2) associated with an adaptive mutation from the post-glacial recolonization of northern Scandinavia and suggested that strong genetic drift at the front of the expansion range caused surfing of the adaptive allele in the newly colonized regions (Klopfstein et al., 2006; Excoffier and Ray, 2008).

In our study, the H_o for putative neutral and outlier SNPs, was positively and significantly correlated with clonal proliferation. This finding supports the claim of the possible presence of a Meselson effect, because accumulation of heterozygous loci was observed in populations with fewer genotypes and in which asexual reproduction was more common and, therefore, counteracted the inbreeding effect in small and isolated populations. However, the simultaneous reduction in genetic diversity (H_e) overwhelmed this positive effect and could probably have an overall negative impact on the survival of the species in Mexico from a genetic point of view (Leimu et al., 2006; Aitken and Whitlock, 2013). The Meselson effect is generally considered a strong indicator of long-term evolution under strict asexuality (Hartfield, 2016; Brandt et al., 2021); however, this phenomenon has also been shown to occur in relatively recent lineages, of less than 100, 000 years of age (Pellino et al., 2013). Some researchers have suggested that quaking aspen clones may be millions of years old and have survived several glacial cycles (Barnes, 1966; Kemperman and Barnes, 1976; Mitton and Grant, 1996).

Asexual reproduction in Mexican Populus tremuloides may play a contributing role in the ability of this species to survive under extreme ecological conditions (Goessen et al., 2022). For example, quaking aspen has developed physiological and ecological adaptations to survive and proliferate in diverse ecological environments, sometimes even extreme environments, such as the mountainous regions of the Sierra Madre Occidental in Mexico and areas experiencing severe drought (Ding et al., 2017; Rogers et al., 2020). Ding et al. (2017) considered the Sierra Madre mountain range in northeastern Mexico one of the few locations where quaking aspen populations show a moderate to high probability of surviving multiple glaciations.

5. Conclusions

This study contributes to a better understanding of the survival and adaptation to changing environmental conditions of natural clonal plant populations, especially of Populus tremuloides, one of the most widespread and ecologically important tree species in the Northern Hemisphere.

However, further research on the Mexican populations of Populus tremuloides is required, given the importance of obtaining information from natural populations located at the border of the general climatic niche of the species, where projections are not at all encouraging: (i) because global climate change is expected to transform the distribution of large numbers of forest trees (Noss, 2001), and (ii) it is estimated that about 26% of the current geographical distribution will no longer be suitable for quaking aspen by 2060, and that the loss of habitat will be particularly marked in the southern distribution range of the species, where the Mexican populations are distributed (Worrall et al., 2013). It has been suggested that these Mexican populations have developed strategies for surviving adverse conditions that are unique in the plant kingdom, in contrast to findings on populations of the same species in the southwestern USA (Crouch et al., 2023). It is important to continue exploring these adaptive mechanisms, as this would enable prediction of the impacts of climate change, which could be helpful for designing effective mitigation strategies.

Acknowledgments

We thank the Mexican Consejo Nacional de Humanidades, Ciencias y Tecnologías (CONAHCYT) for the financial support provided to the first author to carry out his training in the Institutional Doctoral Program in Agricultural and Forestry Sciences (PIDCAF-UJED) with Scholarship No. 334852 and financial support with agreement number CONACYT-FRQ-2016: 279459 for the project "Genome-wide scans for detecting adaptation to climate and soil in Populus tremuloides as the most widely distributed tree species in North América". Dr. Jesús M. Olivas-García assisted in the sampling in the state of Chihuahua, Mexico, and Katrin Groppe, Thünen Institute of Forest Genetics, Germany, provided excellent lab work. The Emerging Leaders of the Americas Program (ELAP) of the Government of Canada awarded a scholarship and the Institute of Integrative and Systems Biology (IBIS) of Laval University allowed the use of its campus and contributed to the training of the first author.

CRediT authorship contribution statement

Javier Hernández-Velasco: Writing – review & editing, Writing – original draft, Resources, Formal analysis, Data curation. José Ciro Hernández-Díaz: Writing – review & editing, Resources. Sergio Leonel Simental-Rodríguez: Writing – review & editing, Resources. Juan P. Jaramillo-Correa: Writing – review & editing, Resources. David S. Gernandt: Writing – review & editing, Resources. José Jesús Vargas-Hernández: Writing – review & editing, Resources. Ilga Porth: Writing – review & editing, Funding acquisition, Formal analysis, Data curation. Roos Goessen: Writing – review & editing, Formal analysis, Data curation. M. Socorro González-Elizondo: Writing – review & editing, Resources. Matthias Fladung: Writing – review & editing. Cuauhtémoc Sáenz-Romero: Writing – review & editing, Resources. José Guadalupe Martínez-Ávalos: Writing – review & editing, Resources. Artemio Carrillo-Parra: Writing – review & editing, Resources. Eduardo Mendoza-Maya: Writing – review & editing, Resources. Arnulfo Blanco-García: Writing – review & editing, Resources. Christian Wehenkel: Writing – review & editing, Writing – original draft, Resources, Funding acquisition, Formal analysis, Data curation.

Data availability

A filtered variant of the SNP file used for the population genomic analyses has been deposited in the Dryad repository under the accession doi: https://doi.org/10.5061/dryad.ttdz08m82. Scripts used for the analyses are available at https://github.com/JavierHdz88/Causes-of-heterozygosity-excess.

Declaration of competing interest

The authors have no competing interest to declare.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2024.12.006.