A map of climate change-driven natural selection in Arabidopsis thaliana

Through the lens of evolution, climate change is an agent of directional selection that forces populations to change and adapt, or face extinction. Current assessments of the risks associated with climate change1,2, however, do not typically take into account that natural selection can dramatically impact the genetic makeup of populations3. We made use of extensive genome information in Arabidopsis thaliana and measured how rainfall-manipulation affected the fitness of 517 natural lines grown in Spain and Germany. This allowed us to directly infer selection at the genetic level4. Natural selection was particularly strong in the hot-dry Spanish location, killing 63% of lines and significantly changing the frequency of ∼5% of all genome-wide variants. A significant proportion of this selection over variants could be predicted from climate (mis)match between experimental sites and the geographic areas of where variants are found (R2=29-52%). Field-validated predictions across the species range indicated that Mediterranean and Western Siberia populations — at the edges of the species’ environmental limits — currently experience the strongest climate-driven selection, and Central Europeans the weakest. With rapidly increasing droughts and rising temperatures in Europe5, we forecast a wave of directional selection moving North, putting many native A thaliana populations at evolutionary risk.

Through the lens of evolu on, climate change is an agent of direc onal selec on that forces popula ons to change and adapt, or face ex nc on. Current assessments of the risks associated with climate change 1,2 , however, do not typically take into account that natural selec on can drama cally impact the gene c makeup of popula ons 3 . We made use of extensive genome informa on in Arabidopsis thaliana and measured how rainfall-manipula on affected the fitness of 517 natural lines grown in Spain and Germany. This allowed us to directly infer selec on at the gene c level 4 . Natural selec on was par cularly strong in the hot-dry Spanish loca on, killing 63% of lines and significantly changing the frequency of~5% of all genome-wide variants. A significant propor on of this selec on over variants could be predicted from climate (mis)match between experimental sites and the geographic areas of where variants are found (R 2 =29-52%).

Field-validated predic ons across the species range indicated that Mediterranean and Western
Siberia popula ons -at the edges of the species' environmental limits -currently experience the strongest climate-driven selec on, and Central Europeans the weakest. With rapidly increasing droughts and rising temperatures in Europe 5 , we forecast a wave of direc onal selec on moving North, pu ng many na ve A. thaliana popula ons at evolu onary risk.
To predict the future impact of climate change on biodiversity, the typical star ng point has been clima c tolerances inferred from the current species distribu ons. These tolerances are usually treated as sta c, and risks are assessed based on whether species' environmental niches will shrink 1,2 or shi faster than the species can migrate 1,6 . However, these approaches do not account for within-species gene c varia on, and for natural selec on causing species to gene cally change and adapt over me 3 . To predict the "evolu onary impact" of climate change on a species, i.e. how much gene c change is required for adapta on to climate change, we thus need to quan fy and model environment-driven natural selec on at the gene c level. Thanks to species-wide genome scans [7][8][9] , as well as genome associa ons with climate of origin [10][11][12][13][14] , we increasingly understand the genomic basis of past selec on and climate adapta on, which has been used to es mate future adapta on debt or "genomic vulnerability" 10,11 . Natural selec on, however, is only indirectly inferred in the types of analyses discussed above. The best way to directly quan fy selec on in a specific environment is provided by field experiments in which mul ple genotypes of a species are grown together in a common environment 15,16 . With such experiments, rela ve fitness can be directly associated with gene c varia on across popula ons 4,[17][18][19] . Ideally, one would carry out such field experiments at many different sites throughout the species range, but this is rarely prac cal 20 . Nevertheless, an emergent finding is that individuals are normally locally adapted and that local genotypes are o en posi vely selected over foreigners in their "home" environment, while nega vely selected in their "away" environments 21,22 . An intui ve conclusion is that it should be possible to derive a metric of natural selec on from the extent of change in local climate at an individual's home. Here we combine high-throughput associa ons of genome and current climate varia on with experimentally quan fied in situ natural selec on in the plant Arabidopsis thaliana . We exploit these associa ons to forecast natural selec on driven by future climate change, and how it impacts the genomic varia on of a species across its geographic range -what we interpret as a new metric of evolu onary risk of popula ons.
To study climate change-driven natural selec on in the annual plant A. thaliana , we In each experimental environment, we quan fied genome-wide selec on at the gene c level based on the difference in rela ve fitness of lines with the minor and the major allele at each genomic posi on (1,353,386 biallelic SNPs across 515 lines with high-quality genome informa on) ( Fig. 1 ). Our approach iden fies both causal variants, as well as many more variants that are in significant linkage disequilibrium (LD) with causal variants 24,25 -due to so-called background selec on or gene c hitchhiking. We use the term allelic selec on differen als ( , called total selec on by Thurman and Barre 4 ), to denote the realized selec on affec ng each SNP resul ng from the combina on of selec on ac ng directly on the focal variant, and the indirect effects due to selec on on causal SNPs that are in LD with the focal variant. Calcula ng allelic selec on differen als using Linear Models (LM-GEMMA, ref. 26, see Supplemental Methods V.3 ), we found a total of 421,962 SNPs with allelic selec on differen als below a 0.05 significance threshold ( Benjamini & Hochberg FDR correc on) in at least one of the eight environments (see Table S2 ) . Using more stringent Bonferroni correc on (<7x10 -7 ), we s ll detected 6,538 SNPs distributed throughout the genome, sugges ng that the polygenic model of natural selec on 27 prevails in this climate-manipula on experiment. These high numbers are not surprising, given that we expect to capture many SNPs that are only indirectly selected. Thinking about our experiment as studying a popula on of plants with mul ple genotypes, the change of allele frequencies in response to one genera on of selec on would be up to 10% in Spain and low precipita on, while it would not exceed 2% in the benign high-precipita on environment in Germany (see Supplemental Methods V , Fig. S9 ).
While variants inferred to be under posi ve or nega ve selec on a er Bonferroni-correc on were overall more likely to be located in intergenic regions than in genes (Fisher's Exact test Odds ra o [Odds]=1.11, P =7x10 -30 ), such variants were enriched for nonsynonymous muta ons (Odds=1.05, P =2x10 -4 ). The large number of variants affected by selec on implies a strong turnover of varia on across the en re genome as a response to the environment 28,29 , and a poten ally significant demographic decima on -what Haldane called "the demographic cost of natural selec on" 30 .
Changes in allele frequency are not only determined by the adap ve value of a variant but also the alleles it is linked to. We therefore improved the detec on of direct targets of selec on by correc ng for LD-driven effects 25,31 using Bayesian Sparse Linear Mixed Model associa ons with rela ve fitness (BSLMM-GEMMA, ref. 31), see Supplemental Methods V ). This analysis indicated that the frac on of the genome likely to be a target of selec on was only 8x10 -5 -5x10 -6 . This frac on was much smaller than what we had iden fied with significant allelic selec on differen als (2x10 -5 -0.001), confirming that selec on must be mostly indirect 4,32,33 . We studied whether the direc on and intensity of selec on was dependent on the experimental environment. Alleles that were posi vely selected under low precipita on tended to be nega vely selected under high precipita on, and vice versa, so called antagonis c pleiotropy 18  bio12: P <10 -5 ), but no differences in sweep likelihood ( P =0.9). Implemen ng genome-wide environmental niche models 10 (see Supplemental Methods VII ), we found that alleles selected in Germany and high precipita on were more likely to come from higher la tudes ( Fig. 1 C), while the opposite was true for alleles selected in Spain and low precipita on ( Fig. 1 D). In agreement, similarity in precipita on regime was a good predictor for alleles selected in Spain  We finally aimed to build an environmental model that can predict allelic selec on differen als based on the climate and diversity pa erns. We used as predictors the per-allele associa ons with climate of origin from clima c GWA as well as the signatures of past selec on at each SNP, , , and sweep likelihood, and their genome annota ons -all predictors were derived from public databases ( worldclim.org , 1001genomes.org , arabidopsis.org ). We used a regression with decision trees using Random Forests to build Genome-wide Environment Selec on ("GWES") models. Conceptually, GWES models are similar to concepts related to Environmental Niche Models (ENMs), but instead of training them with presence/absence data of a gene c variant 10,11 , we trained them with our measured allelic selec on differen als. This provided a means to predict whether alleles should increase/decrease in frequency in a certain climate, instead of merely an indica on of whether alleles are likely to be present, which is the indirect ENMs' version.
By training models jointly with experimental data from Spain and Germany and using cross-valida on, we confirmed that inferred and measured selec on differen als were correctly predicted, with a high correla on accuracy (0.56 < Pearson's r < 0.7) and explaining a large propor on of variance (R 2 = 29-52%) ( Fig. 3   a be er or worse posi on to respond to future climate than our global set of more diverse lines. We therefore looked for SNPs predicted to change most strongly in selec on by 2050 (fitness advantage or disadvantage changed over 5%), and evaluated whether the allele posi vely changing in selec on is locally present ( Fig. 3 G, Fig. S14 ). We found that most local alleles will become more nega vely selected; we therefore predict that the degree of local adapta on will decrease for many na ve popula ons ( Fig. 3 G), leading to an adapta on deficit.

Conclusion
The expected changes in climate during the 21 st century will threaten the survival of many species.
Because the distribu on of gene c diversity is so well characterized in A. thaliana , we have used it to address the challenge of predic ng the effects of climate-driven natural selec on genomic varia on across a species' range. Integra on of genome-climate associa ons with direct fitness observa ons allowed us to build models that directly predict selec on at the gene c level rather than mere probability of presence/absence of variants. This informa on enabled us to infer range-wide evolu onary vulnerability in the face of rapid climate change. The first two steps in our project, assembling a range-wide collec on and genome sequencing of a number of diverse lines, are in reach for many species of plants. A greater challenge is the genera on of fitness data, but this can be par ally solved by iden fying par cularly informa ve field sites -as we have done in our studyand by exploi ng the immense progress in field phenotyping at different scales 48,49 . Combining such observa ons with our new genome-wide environment modeling approach will help us to fully incorporate evolu on into predic ng the impacts of climate change on biodiversity.        Table S1 . The experiment resulted in observa ons from 23,154 pots.
Three measurements of fitness were produced: survival from seed to reproduc ve adult (propor on 0-1) and the average fecundity per reproduc ve adult (inflorescence skeleton lengths ranged from 18,400 to 1,622,000 pixels, which approximately corresponds to 1 to 6,127 seeds per plant).
Fecundity was only measured for plants with at least one fruit. We finally calculated an integrated life me fitness value by mul plying the survival propor on to adulthood with the total offspring produced. Data from only 515 accessions were used for subsequent analyses, because 2 accessions had insufficient genome informa on.

II. 1001 Genomes Project data
We

III. F st and selec ve sweep signatures from polymorphism data
We used the gene c groups previously defined for the same accessions 10

III.1 Geographic proxies of diversity metrics
In order to es mate and a proxy of Tajima's D at a regional scale, we used the 4 closest neighbouring accessions in our set (same pa erns were observed with three neighbours within a geographic area of 5°la tude-longitude radius), and computed the total number of polymorphisms P in the subset and the sum of all pairwise Hamming differences, H. Then we calculated , and D as: Where is the genome size, are all SNPs with full informa on that were used to count polymorphisms and distances, and are all SNPs of the genome matrix. In the denominators, 6 is the number of pairwise comparisons of four genomes, and 1.8666 is the harmonic number of 4.
Although D is normally divided by the standard error, we only wanted to rank our natural lines so we used the difference between and as a proxy of D.

IV. Heritability of fitness
To es mate how much variance in fitness is related to the genotypes of the lines, we used generalized linear mixed models using the R package MCMCglmm (ref. 53 ). We used fitness es mates per replicate and, appart from including the natural line ID, we controlled for block (growing tray) and posi on within the block (longitudinal, la tudinal, and the interac on). As this is a Bayesian approach, we used flat priors, we used 10,000 MCMC steps, a burnin of 10%, and confirmed that this was sufficient for convergence of the chain. For survival propor on we used a Binomial link, for number of seeds we used a Poisson link, and for the combined life me rela ve fitness we used a Gaussian link. The mode and 95% Highest Posterior Density of the posterior distribu on of each random effect were extracted ( Table S3 ).

V. Genome-Wide Associa on, selec on differen als, and direct selec on es mates
We .
The BSLMM model is also useful also used to calculate the propor on of variance explained ('chip heritability'). To do this, we used the last 1,000 samples of the MCMC chain and calculated the 95% Highest Posterior Density Interval (95% HPD), for which we report the median and the 2.5% and 97.5% percen les.
To illustrate the differences between a univariate (LM) and mul variate regressions (similar to BSLMM), we show an example of two SNP predictors, and . For mathema cal convenience we assume that the response variable fitness, , as well as the predictors, are mean centered and variance scaled. From the univariate approach, where the effect of a SNP is es mated marginally or independently, would be: This would be es mated separately for SNP one and two. In a mul variate regression framework, the regression coefficient, called condi onal or par al coefficient, , is corrected by the correla on between the two predictors, , which in the case of genotypes is called linkage disequilibrium, as in the form of: Thus we find an analogy between and and and from popula on gene cs (eq. 4 from 54 simplified as conceptual model),where a selec on differen al is dependent on the true selec on and an indirect term dependent on all other SNPs in the genome.

V.1.1 Across field experiments
We looked for gene c variants with a posi ve selec on differen al in one experimental environment that had a nega ve differen al in another (antagonis c pleiotropy). We also asked how o en a selected variant in one environment was neutral in another environment (condi onal neutrality).
Only for the purpose of two-environment comparisons, to calculate the odds ra o, we considered SNPs whose selec on differen al P -value was lower than 0.01 in one environment as condi onally neutral, while antagonis c pleiotropy SNPs were ones whose P -values in both environments were < 0.02 (because antagonis c pleiotropy requires two tests, one in each environment, the significance threshold should be 2 x P -value in each test) (see Fig. 2).

V.1.1 Across life history stages
Calcula ng allelic selec on differen als for survival and fecundity separately, we found no correla on between survival-only and fecundity-only es mates (r<0.07, Fig. S4-6 ), consistent with different stages of a plant being differen ally affected by environmentally imposed selec on 55 .

V.2 Intensity of selec on
The distribu on of absolute allelic selec on differen als, , has a shape resembling that of an exponen al. We calculated the expected rate using Maximum Likelihood op miza on in R, which can also can be approximated as the inverse of the mean: . For this, we use or the mean of as a metric of the overall intensity of selec on ( Fig. 1 B,   Fig. 3 D).

V.3 Important notes on popula on structure correc on in a wild species
The goal of genome-wide associa on studies in humans is the iden fica on of individual SNPs that are causal for traits such as disease suscep bility. It is therefore impera ve to penalize SNPs that are correlated with popula on structure, as the lack of controlled experiments can generate spurious associa ons between gene c variants more common in some human ancestry groups with regional measurement errors or different cultural and nutri onal environments 56

VII.2. Environmental Niche Models
Genome Environmental Niche Models (GEMs) were fit using decision trees with presence/absence of SNPs as response variable and the climate variables described in the previous sec on and la tude and longitude a predictors; as described 10  We used these models to predict from raster maps of the climate layers a probability between 0 and 1 that the alterna ve allele was in a map cell. We judge this as a more appropriate output than a discrete 0/1 outcomes, as some mes alleles were widespread or at intermediate frequencies in many regions and thus their environment niche was not strictly defined.

VII.3. Climate variability
To study spa al climate variability, for each Arabidopsis natural line, we extracted climate variables ( Table S8 ) in a 50 Km buffer where they were originally collected from and calculated the coefficient of varia on (CV) across grid cells.
To study temporal variability, we used climate data from 1958-2017 39 to calculate annual precipita on values for each popula on, from which we in turn derived the inter-annual CV..

VIII.1 The model
We used a decision tree approach with Random Forest using the R package randomForest (ref. 62,63) to predict the vector (n=1,353,386) of GWA results with rela ve fitness in one environment, which we call allelic selec on differen als s , from a 1,353,386 x 98 matrix of GWA associa ons with climate variables, ( Table S8 ,  . In the cases where we trained models with two environments, we also included the 2 x 98 climate variables at our field sta ons: . Because training a Random Forest with the full dataset would be computa onally expensive, we only trained with 10,000 observa ons (with smaller and larger SNP sets, we had determined that training with more than 10,000 observa ons did not improve predic ons). To test accuracy and bias we used a different set of 10,000 SNPs, divided into 100 bootstrap samples, and we report the intervals of the 95% bootstrap distribu on. The results presented in Fig. 3 were produced with 10,000 randomly drawn SNPs across the genome. To confirm that there was no confounding from non-independent samples in the training and tes ng SNPs, we repeated all analyses, training with 10,000 random SNPs from chromosome 1 and tes ng with 10,000 random SNPs from the four other chromosomes. There were no substan al changes in predictability.
Several combina ons of training and tes ng were performed to validate the predic ons of "unobserved" environments. We did four model fi ngs, each me using training observa ons from three environments only, and then used observa ons from a fourth environment for tes ng. The final model, used for produc on maps, was trained with observa ons of four experiments and tested with the same four environments, but different observa ons.

VIII.2 A simple visualiza on of environmental distance
To visualize more directly the rela onship between allelic selec on differen als at a loca on and the overall environment where the alleles are found, we calculated the distances between the field sta on's climate and the allele's home environment (as defined below). We started using the

VIII.3 Note on limita ons and interpreta ons
As in any predic ve exercise, our projec ons have limita ons (discussed below). We nevertheless firmly believe that they are indispensable to move forward in the field of forecas ng climate impacts.
Models such as ours are tremendously useful for subsequent experimental valida on (as we are currently doing through an experimental evolu on network: GrENE-net.org ) or with in situ observa ons collected as we move into the future (e.g. iNaturalist.org , iSpot.org ). This itera ve predic on ↔ valida on process will be key to advancing the complex field of predic ng the effects of climate change on biodiversity.
The limita ons, enumerated and discussed: E. Bet-hedging strategies such as seedbank demographic dynamics buffer allele frequency changes over me.

IX. Re-analysis of published data from common garden experiments
For a conceptual diagram of predictability valida on with external common garden datasets, see

IX.2 Fournier-Level et al. 2011
Fournier-Level and colleagues 36 germinated seeds in greenhouses, and two weeks a er germina on (established seedling stage), they transplanted seedlings to outdoor field sta ons where one plant was transplanted in one pot. They counted how many transplanted seedlings survived to reproduc on (par al survival propor on), and the number of fruits per plant (reproduc on, seed set). We again computed life me fitness as the product of par al survival and reproduc on.
We excluded the experiment in Finland in downstream analyses because only 58 natural lines were planted there in the original publica on 36 and because later we verified the imputa on accuracy was very low (Pearson's r<0.008).

IX.4 Sanity checks for imputa on and geographic predic ons
We carried out sanity checks to ensure that the imputed fitness from other experiments was not just an ar factual phenotype with the same structure as the rela onship matrix. This would mislead us to think there is predictability, as we would expect that allele selec on differen als calculated in such ar factual phenotype would depend on popula on structure and thus would likely be predictable from climate structure alone.
We shuffled the genotype iden es from Fournier-Level et al. and Manzano-Piedras et al.
with their fitness values. Then we repeated the GBLUP analysis with 50 rounds of shuffling and computed heritabili es and predic on accuracies. We confirmed that heritability with shuffled data was negligible (1x10 -9 < h 2 <1.6 -3 ) and so was the accuracy of imputa on (-0.047310 < r < 0.070380).
This indicated that in the absence of true heritable varia on, imputa on of fitness would be random and not dependent on the rela onship matrix.
We also were concerned that geographic predic ons could be driven by some underlying bias in our analyses, i.e. bias inherent to geographic sampling, popula on history of genotypes chosen, etc. In other words, we were concerned that the null expecta on of predictability would be non-zero. As before, we randomized fitness values with genotypes for all six datasets. Then, we repeated the GWA to es mate allelic selec on differen als (as Fig. 1), and trained different combina ons of GWES models to re-predict allelic selec on differen als at teach loca on based on climate (as Fig. 3). We confirmed that, differently from the analyses of real data presented in Fig. 3, there was no significant predictability ( Fig. S15 ).

IX.5 An explana on for "inverse predictability"
We no ced that using only our two experiments for model training, there was "inverse         Rela onships between rela ve fitness effect, rela ve fitness effect size, and P -values (calculated from GWA with rela ve fitness) and minor allele frequency of alleles in each environment.

Figure S10. Environmental distance and selec on differen als
As Fig. 2B, but for allele selec on differen als in Germany under high precipita on.

Figure S11. Future change in selec on for different climate change scenarios
Same as Fig. 3G, but for different climate change scenarios. The higher the predicted CO 2 emissions, the stronger the predicted increase in selec on intensity.

Figure S12. Field valida on conceptual chart
Conceptual workflow on field valida on procedure with data from published experimens.

Figure S13. Null expecta on of predictability
Same as Fig. 3, but with randomized fitness values associated to genotypes. We could not find any model combina on that had non-zero predictability (95% bootstrap confidence overlaps with zero).
This proof of concept indicates that the predictability we find must have a biological basis, in which the combina on of climate of origin for a gene c variant and the local climate allows to infer selec on over such a variant.  Figure S14. Change in selec on rela ve to local diversity Same as Fig. 3D, but coun ng the number of local alleles increasing or decreasing in selec on (total n= 10,752 SNPs ). Only changes with more than 5% advantage/disadvantage were considered (defined a posteriori from Bonferroni-significant alleles, which generated at least 5% effect in fitness).

SUPPLEMENTAL TABLES & DATASETS
SupplementaL tables are available in the online version of the paper with doi: xyz.

Table S1. List of accessions
The 1001 Genomes Project iden fier, country of origin, and la tude and longitude are reported.       Table S7.

Predictability of environmental models
A er training GWES models with a set of experiments, we inferred allelic selec on differen als on another set of experiments and compared those with the real allelic selec on differen als. We calculated Pearson's product-moment correla on r and % of variance explained R 2 using a regression. 95% confidence intervals were calculated with 100 bootstrap replicates.  into our 517 global accessions. We report heritability, Pearson's r between GBLUP predicted fitness and real fitness, and the significance of the correla on test.

Table S10. Correla on between inferred natural selec on intensity and other variables
Spearman's r between selec on intensity and diversity metrics or climate metrics is given.