Climate change deeply shapes the geographic range of diverse species, reduces biodiversity, and affects the survival of numerous species^1,2. Understanding how species respond to climate change promotes species conservation and crop breeding.

Species distribution models (SDM) has been widely used to predict species range shift over time. However, SDM ignores evolutionary process. The evolutionary potential, namely, genetic diversity, affects the adaptation ability of natural populations to climate change³. Mutation facilitates the ability of species’ response to future climate change, but at the same time, could accumulate deleterious mutations and thus reduce fitness, which is called genetic load^4,5. Accordingly, we can combine two strategies to study how species respond to climate change based on population genomics. One is to screen adaptive mutations, which is well-known and could speed up adaptation^6,7,8. The other strategy is to evaluate genetic load (deleterious mutations) in the natural populations, in which high genetic load reduces fitness and adaptation potential⁹.

Many studies have revealed the genetic basis of local adaptation^6,7,8, built the relationship of adaptive mutations and present climatic factors, and predicted the responses to future climate change based on the relationship^10,11,12. Accordingly, one strategy of species conservation is to preserve functional genetic diversity associated with fitness¹³. However, adaptive mutations alone cannot determine the evolutionary fate of a population^3,14. At genomic level, beneficial mutations are rare, and most mutations are neutral or deleterious^5,15. The accumulation of deleterious mutations would reduce individual fitness, and accordingly, the potential for a population to respond to climate change would reduce if it had high genetic load^3,16. Therefore, study of genetic load could benefit the prediction of species’ response to climate change and species’ resilience^3,17. Nevertheless, most studies only focused on one aspect and ignored the other. Therefore, incorporating genetic load and other factors to predict climate change vulnerability is crucial for species conservation and crop breeding.

Knowing the driving forces of genetic load facilitates the assessment of population loss and species resilience. The accumulation of genetic load depends on selection coefficient, dominance coefficient and counts of deleterious mutations⁹. Given that selection coefficient and dominance coefficient estimation were still challenging, most studies focused on forces affecting deleterious allele counts. Deleterious allele counts was shaped by mutation, recombination, selection, genetic drift and GC biased gene conversion (gBGC)^18,19. Effective population size (N_e) was demonstrated to affect genetic load via changing the effectiveness of selection relative to genetic drift²⁰. Demographic history that reduces N_e, such as bottleneck effect^21,22,23 or founder effect²⁴, range expansion²⁵, introgression from low N_e population²⁶, mating system transition (from outcrossing to selfing)²⁷ and polyploidization²⁸, would induce relaxed purifying selection and increasing genetic load. gBGC mimic the effect of selection, lead to the fixation of G/C alleles, regardless of their fitness effect, thus possibly increasing the genetic load^18,29,30. Recombination affects genetic load via two ways. On the one hand, Hill-Robertson effect would lead to the accumulation of deleterious mutations in low-recombination regions^31,32. On the other hand, stronger gBGC strength in high recombination regions would increase genetic load^29,33. However, the relative contribution of each factor in a single species is not clear.

Arabidopsis thaliana have a human commensal non-relict group spreading beyond its native range and become naturalized worldwide, and several relict groups, which reflects its success as a colonizing species^34,35,36. It is generally interesting to reveal the mystery of its strong adaptation ability. The genetic basis of local adaptation has been explored in A. thaliana^6,7,8 and used to predict the fate of different populations in response to future climate change^11,12. In our study, we shifted our focus to genetic load to explore its variation and determinants along the genome and among diverse populations. We then integrated these findings with genetic offset and SDM to predict the vulnerability of different A. thaliana populations to climate change. Overall, our study highlights the driving forces of genetic load, and reveals how range shifts, adaptation and maladaptation jointly shaped the evolutionary fate of populations in response to climate change.

Identification and characterization of deleterious mutations in natural populations

To study the genetic load of A. thaliana populations, we leveraged 1114 globally distributed non-reference natural accessions that had published resequencing data^{35,37,38,39,40,41}. After SNP calling and quality control, we identified 8,678,999 homologous SNP sites with missing rate less than 10%. Using Arabidopsis lyrata and Capsella rubella as outgroups, we retained 2,958,899 SNP sites with determined ancestral and derived state for further analysis.

Because genetic load is hard to study based on fitness, here we measured with deleterious mutations derived from genomic sequences. We focused on SNPs in coding regions, and generated three datasets from the SNP datasets (Fig. 1A). The neutral datasets (fourfold degenerate sites) were used as control. Protein sequence conservation-based approaches were used to predict deleterious mutations from missense SNPs. For missense SNPs, we either grouped into tolerated missense SNPs (tnSNPs) and deleterious missense SNPs (dnSNPs) based on the prediction of PROVEAN or SIFT, or predicted the deleteriousness score of each SNP using LIST-S2, of which high score indicated high deleteriousness (Fig. 1A). Loss-of-function mutations (LoF), including stop-gain and splice site mutation, were included as the most deleterious SNP datasets (Fig. 1A).

A Neutral and deleterious mutation datasets generated in this study. B The site frequency spectrum (SFS) of neutral and deleterious mutations. DAF, derived allele frequency. 4-fold, fourfold degenerate sites; provean.t, tolerated missense SNPs (tnSNPs) predicted by PROVEAN; provean.d, deleterious missense SNPs (dnSNPs) predicted by PROVEAN; sift.t, tnSNPs predicted by SIFT; sift.d, dnSNPs predicted by SIFT. C The Spearman’s correlation (two-sided) between deleteriousness score of missense mutations and derived allele frequency (DAF). Missense mutations were binned into 1000 equal sized bins according to their deleteriousness scores. D Age distribution of neutral and deleterious mutations. E The Spearman’s correlation between deleteriousness score of missense mutations and mutation age. Missense mutations were binned into 1000 equal sized bins according to their deleteriousness scores. F Mutation spectrum of neutral and deleterious mutations. G The Spearman’s correlation (two-sided) between deleteriousness score of missense mutations and mutation spectrum. S > W, strong (G/C) to weak (A/T) mutations; W > S, weak to strong mutations; W > W/S-S, weak to weak or strong to strong mutations. Missense mutations were binned into 1000 equal sized bins according to their deleteriousness scores. H The Spearman’s correlation (two-sided) between genetic load proxy and fitness related traits (fruit number). Genetic load proxy was calculated as the derived allele counts of dnSNPs or sum of deleteriousness score at individual level, as suggested by previous study⁶². Derived allele counts of neutral SNPs and tnSNPs were used as control. Derived allele counts increased from left to right, and fruit number increased from bottom to top. Accessions were binned into ten equal sized bins according to their derived allele count. n = 26 accessions for each bin, the last bin has 28 accessions. Error bars indicate mean +/−95% confidence intervals. Source data are provided as a Source Data file.

Full size image

Consistent with the effects of purifying selection acting on deleterious mutations, the derived allele frequency (DAF) distribution suggested that LoF mutations and dnSNPs were more skewed toward low-frequency than neutral sites and tnSNPs, and the most deleterious LoF mutations had the rarest DAF (Fig. 1B). Similarly, the deleteriousness score of missense mutations were negatively correlated with DAF (Fig. 1C).

To reveal the age of deleterious mutations, we estimated their age using GEVA⁴², which could date variants based on the sequence divergence surrounding the variants. Due to the effect of purifying selection, deleterious mutations were eliminated quickly and were much younger than neutral sites and tnSNPs (Fig. 1D), which was consistent with that found in humans^42,43. In addition, the most deleterious LoF mutations were much younger than dnSNPs (Fig. 1D), and the deleteriousness score of missense mutations were negatively correlated with the age (Fig. 1E).

The mutation spectrum (base composition) suggested that deleterious mutations had higher proportion of strong (G/C; indicating the three hydrogen bonds between bases) to weak (A/T; denoting the two hydrogen bonds between bases) alleles than weak to strong alleles, especially for the most deleterious LoF mutations (Fig. 1F). Similarly, missense mutations with higher deleteriousness score tended to be strong to weak (S > W) alleles rather than weak to strong (W > S) alleles (Fig. 1G). Given the mutation spectrum of de novo mutation in A. thaliana mutation accumulation (MA) lines is biased towards A/T alleles^44,45, the enrichment of deleterious mutations in S > W mutations could be resulted from their relatively young age. While the enrichment of W > S mutations in neutral sites and tnSNPs probably reflect the effect of gBGC, which tend to fix G/C alleles and thus the older mutations would accumulate W > S alleles. gBGC is expected to cause a right-shifted site frequency spectrum of W > S mutations²⁹. Accordingly, we observed an excess of high-frequency W > S mutations (Supplementary Fig. 1A), indicating the impact of gBGC. In addition, high-recombination regions were enriched with W > S mutations (Supplementary Fig. 1B), which probably result from gBGC or purging of deleterious mutations in high-recombination regions. Although gBGC was assumed to be weak or absent in selfing species, it was directly observed or reported in A. thaliana in previous studies^30,46,47.

To explore the relationship between genome sequence derived genetic load proxy (the derived allele counts of deleterious mutations) and fitness consequences, we used published fruit number phenotype, which is cultivated and measured in greenhouse, and is a proxy for lifetime fitness in annual species⁴⁸. There was a generally negative correlation between genetic load proxy and fruit number (Fig. 1H), even after implementing a control for population structure (p = 0.00471 for LIST-S2 and p = 0.01 for LoF mutations, Supplementary Table 1). While the correlation between derived allele counts of neutral sites or tnSNPs and fruit number were not significant (Fig. 1H). Overall, we identified a deleterious mutation dataset and bridged the gap between sequence based genetic load proxy and fitness consequences, which was rarely demonstrated in previous studies of deleterious mutations in natural populations⁹.

Determinants of genetic load variation along the genome

Although genetic load was much more extensively studied across populations, the variation, determinants and their relative contribution of genetic load along different genomic regions is largely unknown. In addition, parameters such as mutation rate and recombination rate were available at different genomic regions but not different populations. To study the determinants and their relative contribution to genetic load variation along genome, we split A. thaliana genome into multiple non-overlapping 50 kb windows, and estimated genetic load in each window. We scaled the derived allele count of deleterious mutations with that of fourfold degenerate sites, and used it as genetic load proxy, to exclude the confounding effects of linked selection, which makes the unscaled counts reflects polymorphisms rather than genetic load^49,50.

The potential evolutionary forces that shape the genetic load variation include mutation, recombination, selection, genetic drift and gBGC^18,19. Accordingly, mutation rate and recombination rate estimated from de novo mutations and crossover map could measure the strength of mutation and recombination^44,51, N_e measured with fourfold degenerate sites diversity (Ï€_s) could reflect the effectiveness of purifying selection²⁰, and the ratio of DAF of W > S mutations to that of S > W mutations (W > S/S > W) represent gBGC strength³⁰. We estimated these four parameters in each window, and used as predictor variables to evaluate their impact on genetic load. To this end, we used 1,986 windows with all predictors and genetic load proxies available.

Ï€_s was negatively correlated with genetic load, except for LoF mutations (Fig. 2A), suggesting low selection efficiency in low-diversity region. N_e varies along genome, and depends on recombination rate and density of selected sites⁵². Accordingly, recombination rate was positively correlated with Ï€_s and negatively correlated with genetic load (Supplementary Fig. 1C and Fig. 2B). Recombination reduces genetic load via Hill-Robertson effect, but at the same time increases genetic load via gBGC^29,31,32. The negative correlation between recombination rate and genetic load (except for LoF mutation) indicated that recombination affects genetic load mainly via Hill-Robertson effect rather than gBGC in A. thaliana. Similarly, several previous studies on other plants also found the negative correlation between recombination rate and genetic load^53,54,55,56. In contrast to the intuition that high-mutation region would accumulate more mutations, regardless of fitness effects, we found that mutation rate was negatively correlated with genetic load (Fig. 2C). Given the positive relationship between mutation and recombination (Supplementary Fig. 1D), it probably results from the purging of deleterious mutations in high recombination regions.

A Spearman’s correlation (two-sided) between π_s and genetic load proxies. π_s, the nucleotide diversity at fourfold degenerate sites, used as the proxy of effective population size. B Spearman’s correlation (two-sided) between recombination rate and genetic load proxies. C Spearman’s correlation (two-sided) between mutation rate and genetic load proxies. D Spearman’s correlation (two-sided) between W > S/S > W and genetic load proxies. W > S/S > W, ratio of derived allele frequency (DAF) of weak (A/T) to strong(G/C) (WS) alleles to the strong to weak (SW) alleles. E Variable importance of each predictor variable. Scaled permutation based variable importance from random forest model was used as the importance value. provean.d, deleterious missense SNPs (dnSNPs) predicted by PROVEAN; sift.d, dnSNPs predicted by SIFT; LIST-S2, missense mutations with deleteriousness score from LIST-S2. LoF, loss-of-function mutations. Error bars represent mean +/− 95% confidence intervals. Each predictor was grouped into 20 equal-sized bins (n = 100 windows for each bin, the last bin has 86 windows) in ascending order. Source data are provided as a Source Data file.

Full size image

In terms of the effect of gBGC, W > S/S > W was positively correlated with genetic load (Fig. 2D), similar to the reported deleterious mutation accumulation in high W > S/S > W regions³⁰. However, we found that W > S/S > W was higher in low-recombination rather than high-recombination regions as usually expected (Supplementary Fig. 1E). We also estimated gBGC strength using relative rate of fixation of W > S substitutions compared to W > W substitutions, as previously reported²⁹, and also found that it was higher in low-recombination regions (Supplementary Fig. 1F). Accordingly, GC content was higher in low-recombination regions (Supplementary Fig. 1G), and the reason remains unclear.

Based on the non-linear relationship between most predictor variables and genetic load, we fitted a random forest regression model on genetic load proxy, and evaluated the relative contribution of each predictor variable based on the permutation based variable importance. W > S/S > W was not included in the model, considering its unresolved relationship with recombination. The random forest model suggested that π_s was the most important contributor, and the second is recombination rate, while mutation rate has the least importance (Fig. 2E). Taken together, genomic regions with low π_s and low recombination rate would accumulate much more deleterious mutations.

Determinants of genetic load variation among populations

To determine forces shaping genetic load variation among populations, we only consider N_e and gBGC, because it was difficult to estimate mutation rate and recombination rate variation among populations only based on population polymorphism data^57,58. The impact of N_e on genetic load variation has been extensively studied under different demographical scenarios in diverse organisms^21,22,23. In A. thaliana, a human commensal non-relict group originated near Balkans has experienced post-glacial range expansion and thus is a great model to study the effect of N_e on genetic load variation^36,59. We mainly focused on the east-west expansion across Eurasia, and excluded North America population, which was introduced recently and admixture among introduced lineages increased haplotype diversity and reduced genetic load⁶⁰.

The population history estimation using SMC++ suggested that populations at the expansion front, such as Yangtze River basin population (Yangtze), North Sweden population (N. Sweden) and Northwestern China and Central Asia population (NW. China & C. Asia) have undergone severe bottleneck effects and have smaller N_e compared to their origin population Balkans (Supplementary Fig. 2A). Consistently, the nucleotide diversity at fourfold degenerate sites (Ï€_s), which could reflect long-term N_e⁶¹, coincided with the population history estimation (Supplementary Fig. 2B). Therefore, Ï€_s was used as the N_e proxy in further analysis.

Using derived allele count of deleterious mutations as genetic load proxy, as suggested by a previous study⁶², both deleterious missense mutations and LoF mutations accumulated with elevated distance to the putative origin of non-relict expansion³⁶, and the most distant Yangtze River basin population accumulated the highest genetic load (Fig. 3A). Previous studies based on small sets of accessions also found that marginal populations accumulated more deleterious mutations^37,63. 1001 Genomes Project also found that populations in Sweden and Central Asia have more nonsynonymous mutations and LoF mutations³⁵. Because genetic load estimates probably differs using different load proxies^9,62, to validate the robustness of our results to different load proxies, we also used two other load proxies, ratio of derived allele count of deleterious mutations to that of neutral sites and R_A/B. R_A/B measures the relative derived allele counts of dnSNPs that occur more often in the expanded populations than that in Balkans⁵⁵. The results suggested that our conclusion was robust (Supplementary Fig. 2C, D).

A Spatial regression between genetic load and distance to putative origin. The regression model was built with the “spaMMâ€� (4.4.16) package of R⁹⁶. The shaded area was 95% confidence intervals around predictions while controlling for spatial effects. B Prediction of the origin of non-relict populations based on genetic load. Red denotes areas with top 1% correlation. Black star is the most likely origin. Color scale represents Spearman’s rank correlation (two-sided) between distance to each tested putative origin and genetic load. Red dots represent the non-relict accessions used in this study. C Pearson’s correlation (two sided) between genetic load and Ï€_s. Ï€_s, the nucleotide diversity at fourfold degenerate sites, used as the proxy of effective population size (N_e). D Pearson’s correlation (two sided) between GC biased gene conversion (gBGC) strength and Ï€_s, gBGC strength (W > S/S > W) was calculated as the derived allele frequency of W > S mutations to that of S > W mutations. Ï€_s used as the proxies of N_e. E Pearson’s correlation (two sided) between the ratio of weak to strong (W > S) load to strong to weak (S > W) load and Ï€_s. provean.d, deleterious missense SNPs (dnSNPs) predicted by PROVEAN; sift.d, dnSNPs predicted by SIFT; LIST-S2, missense mutations with deleteriousness score from LIST-S2. Source data are provided as a Source Data file.

Full size image

Because the origin of non-relict populations was putative based on polymorphism³⁶, here we predicted the origin of non-relict populations based on expansion load, by maximizing the pattern in Fig. 3A. Briefly, genetic load increases with range expansion, and areas with high positive correlation between expansion distance and genetic load are more likely to be the origin. The results suggested that Balkans is the most likely origin for non-relict expansion (Fig. 3B), similar to the prediction based on polymorphism³⁶.

In terms of the impact of N_e on genetic load, Ï€_s was negatively correlated with genetic load, and could explain 74–94% genetic load variation among populations, depending on deleterious mutation types (Fig. 3C). Since Yangtze River basin population appears as an outlier, we removed it to see if the correlation is still present. The results suggested that the correlation between Ï€_s and genetic load is robust, although to a lesser extent (Supplementary Fig. 2E). Consistently, previous study based on 80 accessions also found that N_e greatly affected deleterious mutation accumulation³⁷. Populations at the expansion front, which also have the smallest N_e, such as Yangtze River basin population, N. Sweden population and NW. China & C. Asia population, accumulated the highest genetic load (Fig. 3C). Consistently, fixed load was higher in populations with more severer bottlenecks (Supplementary Fig. 3), indicating the less efficient purifying selection and the stronger genetic drift in these populations.

To explore the impact of gBGC on genetic load, we split genetic load into three types (W > S, S > W and W > W/S > S) based on base composition, and only the first two could be affected by gBGC. As a directional force, the strength of gBGC was counteracted by random drift. Therefore, N_e was expected to affect gBGC strength, which has been reported in human and mammals^64,65. Consistently, we found the fraction of W > S load was positively correlated with π_s, S > W load was negatively correlated with π_s, while that of W > W/S > S load was not correlated with π_s (Supplementary Fig. 4A–C). In addition, W > S/S > W was positively correlated with π_s (the correlation is robust to Yangtze outlier, Supplementary Fig. 4D), indicating that larger populations would have stronger gBGC strength (Fig. 3D). To explore if N_e contributes to gBGC load variation across A. thaliana populations, we used the ratio of W > S load to S > W load as gBGC load, and found that π_s was positively correlated with gBGC load (Fig. 3E), and the correlation is robust to Yangtze outlier (Supplementary Fig. 4E), suggesting that larger populations accumulated higher gBGC load. Taken together, N_e contributes to genetic load variation via regulating both the effectiveness of purifying selection and gBGC strength. Nevertheless, the purging of deleterious mutations by purifying selection overwhelms the fixation of deleterious mutations by gBGC in populations with large N_e, in terms of the overall negative correlation between π_s and genetic load.

Incorporating genetic load to predict climate change vulnerability

To predict the fate of different A. thaliana populations under future climate change, we incorporated SDM, genetic offset and genetic load. SDM ignores evolutionary process, genetic offset is a measure of maladaptation based on adaptive alleles. However, the fate of a species depends on more than just adaptive alleles. Populations with same genetic offset might differ in genetic load³. By integrating the three factors, we could capture the complexities of evolutionary process and more precisely evaluate the climate change vulnerability.

To estimate the change of habitat, using MaxEnt implemented with ENMeval, we modeled present distribution of A. thaliana and projected into future based on future climate (2040–2060 and 2080–2100) under SSP245 emissions scenario, which is a realistic future scenario that aligns with the current global policies and the level of emission reduction efforts⁶⁶. According to present distribution of A. thaliana, a suitability threshold 0.25 match the distribution of most accessions (Supplementary Fig. 5A), as suggested by a previous study⁶⁷. Consistent with previous studies^67,68, the modeling and projection suggested that future habitat suitability would become worse for southern region, while northern region would get much better (Fig. 4A and Supplementary Fig. 5B).

A Habitat suitability change (future-present) between future and present. Suitability was projected from species distribution models (SDM). B Genetic offset based on 2040–2060 SSP245 projections. We only showed regions with present suitability greater than 0.25, as the threshold matches the distribution of most accessions. C Genetic load projected with the spatial regression of genetic load and distance to putative origin. Genetic load was calculated using deleterious mutations predicted by PROVEAN, and normalized with the genetic load of putative origin. D Climate change vulnerability based on SDM. Vulnerability decreases from q1 to q4. E Climate change vulnerability based on genetic offset. Vulnerability decreases from q1 to q4. F Climate change vulnerability based on genetic load. Vulnerability decreases from q1 to q4. G Incorporating SDM, genetic offset and genetic load to predict the fate of different populations. We added the rank of each factor, low rank value indicates high risk, while high rank value indicates low risk. H Vulnerability of nine non-relict populations based on the integration of SDM, genetic offset and genetic load.

Full size image

To incorporate local adaptation into prediction framework, we identified climate associated adaptive SNPs, built the correlations between allele frequency of adaptive SNPs and 19 bioclimatic variables using gradient forest, and projected allele frequency into future based on future climate (2040–2060 and 2080–2100) under SSP245 emissions scenario. The allele frequency shift between future and present, we termed genetic offset, varied along space and time, with southern range, especially the eastern range reached higher genetic offset, indicating higher vulnerability in these regions (Fig. 4B and Supplementary Fig. 5C).

It is challenging to predict future genetic load based on climate change, because genetic load accumulates over many generations, while climate shifts in a few generations. Therefore, we only incorporate current genetic load to assess the vulnerability of different populations. Based on the spatial regression between distance to putative origin and genetic load (Fig. 3A), we projected genetic load on the map. We also utilized Kriging interpolation method to predict load across the landscape, which allow greater flexibility to account for some of the variability in load estimated. The projection indicated that the eastern range accumulated much higher genetic load, regardless of the deleterious mutation types and projection methods used (Fig. 4C and Supplementary Fig. 5D–G). Consequently, for the following analysis, we used only deleterious mutations predicted by PROVEAN to project genetic load based on spatial regression.

To integrate SDM, genetic offset and genetic load into a prediction framework, we made a sensitivity analysis, by setting different thresholds of SDM suitability, genetic offset and genetic load. We only retained regions with present suitability greater than 0.25 to set different thresholds. For each factor, we used quartiles as the threshold, and divided each map into four vulnerability regions. The vulnerability map suggested that southern Spain, Central Asia and eastern range have higher vulnerability in terms of SDM (Fig. 4D and Supplementary Fig. 6A). Southern and eastern range have higher vulnerability in terms of genetic offset (Fig. 4E and Supplementary Fig. 6B). Northern and eastern range have higher vulnerability in terms of genetic load (Fig. 4F).

Taking these three factors together, we sequentially added the vulnerability rank, and predicted that the eastern range was the most vulnerable region, and populations at the range edge are more vulnerable to climate change (Fig. 4G and Supplementary Fig. 6C). In terms of the nine non-relict populations, Yangtze and NW. China & C. Asia population are the most vulnerable populations, Balkans and Spain are the second, N. Sweden and W. Europe are the third, while C. Europe, S. Sweden and Germany are the least concerned populations (Fig. 4H and Supplementary Fig. 6D).

Climate change threatens biodiversity. Understanding how species respond to climate change paves the way to species conservation and crop breeding. In this study, we incorporated SDM, local adaptation and genetic load to predict the fate of different A. thaliana populations in response to climate change. In particular, we comprehensively dissected the driving forces of genetic load variation, both along different genomic regions and among populations, and integrated genetic load to predict the fate of different populations. Our study showed that the integration of SDM, local adaptation and genetic load, instead of only local adaptation or SDM, could add valuable insight on predicting the fate of natural populations, which will be important for the crop breeding and species conservation.

It is common to estimate genetic load via counting deleterious mutations, but there is a gap between sequence derived load proxy and fitness consequence^9,69, although the negative correlation has been observed between deleterious mutation number and fitness related traits in crops^70,71,72. Here we filled this gap by observing the significant negative correlation between load proxy and fruit number. Nevertheless, it remains unclear how much fraction of the fitness variation could be ascribed to genetic load. Some of the variation in fitness could also be related to environment adapted alleles and their match or mismatch to the greenhouse environment. In addition, we only focused on coding SNPs, noncoding SNPs and structural variations could also contribute substantially to total genetic load and phenotypic variation^40,70,73. Therefore, future study incorporating all these variants are appreciated to improve our understanding of the role of genetic load in fitness and species conservation.

To decipher the driving forces of genetic load, we considered population genetic processes affecting genetic load and their relative contribution both along different genomic regions and among different populations, while most studies focused on only one or two of these processes and paid much attention to population level variation. At different genomic region level, we measured the relative contribution of each process and found that N_e was the major contributor. At population level, N_e could explain 74–94% genetic load variation among populations, depending on the deleterious mutation types. Previous studies conducted within species and among species also highlighted the important role of N_e on genetic load^22,24,25,74. Particularly, we found that the impact of N_e on genetic load not only depends on the effectiveness of purifying selection, but also gBGC strength, which was in contrast to the traditional view that N_e affected genetic load only via effectiveness of purifying selection. Therefore, large populations purge drift load and meanwhile accumulate gBGC load, but purging was much stronger than accumulating. The positive correlation between N_e and gBGC strength was also reported in human and mammals^64,65. However, the correlation between N_e and gBGC strength was not always monotonous in other species⁷⁵.

Genetic load accumulated with the increased distance to putative origin, similar to the study in human²⁵. Utilizing the relation between distance and genetic load, we speculated the origin of non-relict populations, which is similar to that predicted with population polymorphisms³⁶. Accordingly, we projected genetic load across the map based on distance to origin, and leveraged this information to predict population loss, based on the assumption that populations with higher genetic load have less potential to respond to climate change. Due to the contrast between long-term accumulation of genetic load and rapid shift of climate, predicting future genetic load based on climate change is still a challenge. Although the distance to niche centroid has been utilized to predict genetic load³, it is just a gross approximation.

Incorporating SDM, local adaptation and genetic load could provide comprehensive insights on predicting the climate change vulnerability of A. thaliana populations. Focusing on only one of them, one would be biased in making conservation strategies^13,14. In our study, predictions based solely on SDM or genetic offset might suggest that the N. Sweden population is safe, but when genetic load is considered, it may indicate an increased risk for this population. Nevertheless, except for SDM, local adaptation and genetic load, other factors, such as gene flow of adaptive or maladaptive alleles, the dispersal potential of populations, could also affect the response ability of populations to climate change³.

More importantly, systematic integration of all these factors is still a challenge. A previous study introduced the FOLDS model, which could integrate multiple factors and apply empirical thresholds to different factors³. Nevertheless, empirical thresholds were still much ambiguous, and determining the weight of each process is difficult. Here we assigned risk levels for each factor based on quartiles, and summed them to assess the relative risk of each population. Although in its early stages, it can be readily used in conservation and management, with the flexibility to incorporate additional processes. Future study systematically integrated all these factors would enable more precise prediction of the evolutionary fate of populations. Overall, our study highlighted the importance of integration of SDM, local adaptation and genetic load on species conservation and crop breeding.

Datasets and population structure

The whole genome resequencing data of 1,114 globally distributed A. thaliana non-reference accessions were retrieved from previous studies^{35,37,38,39,40,41}. The 1114 non-reference accessions were composed of one relict group and ten non-relict populations⁴⁰.

Paired end reads of the 1114 non-reference accession were first aligned to the reference genome using BWA (0.6.2)⁷⁶. Following alignment, the “rmdupâ€� function of Samtools (0.1.18)⁷⁷ was employed to eliminate duplicate reads resulting from library preparation or sequencing. SNPs were called using Genome Analysis Toolkit (GATK 2.5)⁷⁸. Additional filtering of raw SNPs was conducted using the “VariantFiltrationâ€�, applying criteria such as quality score (Q)  ≥  30, mapping quality (MQ)  ≥  20, quality-by-depth ratio (QD)  ≥  10, ReadPosRankSum ≥ −8.0, depth coverage (DP)  ≥  3, strand bias probability (FS)  ≤  10.0, and having no more than three SNPs within a 10 bp window. Only homozygous SNPs with missing rate less than 10% were used for further analysis. Allele state of each SNP site was determined based on the whole genome alignment of Col-0 (TAIR10) to Arabidopsis lyrata (MN47, v1)⁷⁹ and Capsella rubella (MTE, v1)⁸⁰ using LASTZ (1.04.00)⁸¹. Ancestral alleles were defined as those that matched in both A. lyrata and C. rubella, while derived alleles were those where A. lyrata and C. rubella had identical states, but A. thaliana did not match them.

Deleterious mutation identification

SNPs were annotated with SnpEff(v4.3t)⁸². PROVEAN (1.1.5)⁸³, SIFT 4G⁸⁴ and LIST-S2 (1.10)⁸⁵ were used to predict deleterious missense SNPs. For each SNP site, the derived allele was submitted for prediction. NCBI nonredundant protein database⁸⁶, Uniref90⁸⁷ and UniProt TrEMBL/SwissProt were used as the reference database of PROVEAN, SIFT 4G and LIST-S2, respectively, as recommended by these softwares. In PROVEAN prediction, missense SNPs with score ≤ −2.5 were defined as deleterious SNPs (dnSNPs), while SNPs with score > −2.5 were defined as tolerated SNPs (tnSNPs). In SIFT 4G prediction, missense SNPs with score <0.05 were defined as deleterious SNPs (dnSNPs), while SNPs with score ≥0.05 were defined as tolerated SNPs (tnSNPs). In LIST-S2 prediction, the score ranges from 0 to 1, higher score indicates higher deleteriousness. Loss of function (LoF) mutations (including stop-gain and splice site mutation) were obtained from our previous study⁴⁰.

SNP age estimation

For each non-singleton SNP site with determined ancestral state and passed our quality and missing rate control, GEVA (v1beta)⁴² was used to estimate age distribution. Although in human, a species level parameter of N_e and mutation rate was used⁴², different populations have varied N_e and mutation rate. In A. thaliana, N_e and mutation rate were estimated to be about 300,000 and 7e−9 per site per generation in previous studies^37,44,45.

To test the impact of parameter setting on robustness of age estimation, we randomly selected 1000 SNP sites from the genome, and estimated their age with gradient N_e (from 100,000 to 500,000) and mutation rate (from 5e−9 to 9e−9) setting that fluctuating around the species level parameter. When we evaluated the impact of N_e, mutation rate was set as 7e−9, and for the impact of mutation rate, N_e was set as 300,000. The evaluation suggested that age estimation was robust to parameter setting, based on their high correlation under different parameter setting (Supplementary Fig. 7). Therefore, we used the species level parameter “–Ne 300000â€� and “–mut 7e-9â€� to estimate SNP age. The mean of the composite posterior distribution under joint model was used as the age estimates for a given SNP site.

Genetic load proxy calculation

For population level comparison, three genetic load proxies were calculated, derived allele count of deleterious mutations, ratio of derived allele count of deleterious mutations to that of neutral sites and R_A/B (a relative measure of dnSNPs frequency enhancement). For PROVEAN and SIFT 4G predictions, derived allele count was calculated as the derived allele count of dnSNPs. For LIST-S2 predictions, derived allele count was calculated as the sum of deleteriousness scores. R_A/B was calculated as previously described⁵⁵ between the putative origin Balkans and other non-relict populations to measure the relative derived allele counts of dnSNPs that occur more often in the expanded populations than in Balkans.

Random forest modeling for determinants of genetic load variation along the genome

For comparison along different genomic regions, the reference genome was divided into non-overlapping 50 kb windows, and genetic load proxy was calculated in each window. The ratio of derived allele counts of deleterious mutations to that of fourfold degenerate sites, instead of derived allele count of deleterious mutations, was used as genetic load proxy, considering the varied strength of linked selection along genome.

Mutation rate in each window was calculated from the de novo mutations from 107 mutation accumulation lines⁴⁴. Recombination rate at each window was estimated from a crossover map in 17 F₂ populations derived from 18 A. thaliana lines⁵¹. N_e of each window was calculated as the nucleotide diversity of fourfold degenerate sites (Ï€_s)⁵². The strength of GC-biased gene conversion (gBGC) in each window was measured with the ratio of derived allele frequency of weak-to-strong (W > S) mutations to that of strong-to-weak (S > W) mutations (W > S/S > W), as previously suggested³⁰.

To analyze the impact of four predictors (mutation rate, recombination rate, π_s and W > S/S > W) on genetic load, 1986 windows with all these parameters were used. Each predictor was ranked and split into 20 equal sized bins, and Spearman’s rank correlation test was performed to test the relationship between these predictors and genetic load proxy.

Random forest model implemented with R package “caretâ€� (6.0.94) (method = “rangerâ€�) was fitted on genetic load proxies of each window to explore the relative contribution of mutation rate, recombination rate, and Ï€_s to genetic load variation along the genome⁸⁸. A fivefold cross-validation procedure was used to evaluate the accuracy of model, by using 70% of the datasets as training data, and the remained 30% as testing data. The scaled permutation based variable importance for each predictor was obtained from the random forest model “varImpâ€� function and used to measure their contribution.

Population history and N
_e estimation

SMC++(v1.15.4)⁸⁹ was used to estimate the historical population size in each population. As A. thaliana is a selfing species, for each population, 20 pseudo-diploid genotypes were generated by randomly selecting 40 individuals. Mutation rate was set to 7e−9 based on previous studies^44,45, and the timepoints was set to “100 100000â€�. Composite likelihood was formed by create five distinct datasets, by varying the identity of the distinguished individual and treating the remaining samples as “undistinguishedâ€�.

Predicting the origin of non-relict expansion based on expansion load

A map with 10 arc-minutes geographic resolution was downloaded from worldclim v2, and the range was limited to 20 W to 160 E, and 0 to 69 N, based on the distribution of non-relicts. Each coordinate on the map was assumed to be the origin, and distances from the putative origin to non-relict accessions were calculated with ‘rdist.earth’ function in the “fields� (14.1) package of R. The Spearman’s rank correlation between distance and genetic load was calculated for each putative origin, and plotted on the map. The top 1% region with the highest correlation was defined as the putative origin area, and coordinates with the highest correlation was assigned as the most likely origin.

Species distribution modeling

To model the species distribution of A. thaliana, occurrence data from Global Biodiversity Information Facility (GBIF) and accessions used in this study were used. The occurrence data was further thinned to one sample per 50 km using the “poThin� function of “RSpatial� (0.3.0) package of R to reduce sampling bias.

19 bioclimate variables of present and future (2040–2060 and 2080–2100 under SSP245 emissions scenario) were downloaded from worldclim v2. To achieve a balance between reducing multicollinearity and retaining meaningful variability, we selected Bio3, Bio6, Bio7, Bio8, Bio10, Bio15, Bio16 and Bio17 for modeling. These variables were chosen based on their relatively low correlation (Pearson’s correlation coefficients among variables <0.75) and hypothesized importance, as previously described⁶⁷. MaxEnt was used to model species distribution⁹⁰. “ENMevalâ€� (0.3.1) package of R⁹¹ was used to optimize parameters of MaxEnt, using “randomkfoldâ€� for cross validation. Model with the lowest Akaike information criterion (AIC) value was used to build and project species distribution under present and future climate.

Genetic offset calculation

“pcadaptâ€� (4.3.3) package⁹², latent factor mixed models (LFMM) model from “LEAâ€� (2.8.0) package⁹³ and omegaplus (3.0.3)⁹⁴ were used to scan adaptive SNPs (SNPs with top 5% score), respectively. SNPs identified as adaptive SNPs by all software were regarded as the adaptive SNP datasets. Random SNP datasets was generated by randomly sampling the same number of SNPs as the adaptive SNP datasets. The correlation between allele frequency and 19 bioclimate variables from worldclim v2 was built, and genetic offset, the Euclidean distance between present and future values, in future climate (2040–2060 and 2080–2100) under SSP245 emission scenarios was projected with the “gradientForestâ€� (0.1.32) package⁹⁵.

Statistical analyses

All statistical analyses were performed in R.