Ecological, genetic and evolutionary drivers of regional genetic differentiation in Arabidopsis thaliana

Background Disentangling the drivers of genetic differentiation is one of the cornerstones in evolution. This is because genetic diversity, and the way in which it is partitioned within and among populations across space, is an important asset for the ability of populations to adapt and persist in changing environments. We tested three major hypotheses accounting for genetic differentiation—isolation-by-distance (IBD), isolation-by-environment (IBE) and isolation-by-resistance (IBR)—in the annual plant Arabidopsis thaliana across the Iberian Peninsula, the region with the largest genomic diversity. To that end, we sampled, genotyped with genome-wide SNPs, and analyzed 1772 individuals from 278 populations distributed across the Iberian Peninsula. Results IBD, and to a lesser extent IBE, were the most important drivers of genetic differentiation in A. thaliana. In other words, dispersal limitation, genetic drift, and to a lesser extent local adaptation to environmental gradients, accounted for the within- and among-population distribution of genetic diversity. Analyses applied to the four Iberian genetic clusters, which represent the joint outcome of the long demographic and adaptive history of the species in the region, showed similar results except for one cluster, in which IBR (a function of landscape heterogeneity) was the most important driver of genetic differentiation. Using spatial hierarchical Bayesian models, we found that precipitation seasonality and topsoil pH chiefly accounted for the geographic distribution of genetic diversity in Iberian A. thaliana. Conclusions Overall, the interplay between the influence of precipitation seasonality on genetic diversity and the effect of restricted dispersal and genetic drift on genetic differentiation emerges as the major forces underlying the evolutionary trajectory of Iberian A. thaliana.

differentiation. Nonetheless, it must be noted that genetic differentiation is a spatially-explicit phenomenon, as genetic differentiation strictly depends on the genetic relationship that a given population has with its neighbors near and far.
The inherent spatial nature of genetic differentiation defines the theoretical and methodological framework of three models, which are not mutually exclusive, testing the major drivers of genetic differentiation: isolation-bydistance (IBD hereafter), isolation-by-environment (IBE hereafter) and isolation-by-resistance (IBR hereafter) models. In the classical IBD, genetic differentiation among populations exhibits a positive relationship with geographic distance [13][14][15][16][17][18]. In this case, dispersal limitation and genetic drift determine the greater genetic differentiation at larger distances. In fact, limited dispersal constrains gene flow among populations, which is not able to counteract the effect of genetic drift within populations [16][17][18]. Many types of organisms exhibit IBD [19], probably because unrestricted gene flow hardly occurs in nature. Nevertheless, quantifying the contribution of limited dispersal and genetic drift to genetic differentiation is not a straightforward task, as we largely ignore the actual extent of dispersal, the effective population sizes conditioning genetic drift, and the effects of historical factors shaping the genetic relationships among populations [20].
In contrast, the IBE model deals with the effects of environmental differences on genetic differentiation. IBE posits that gene exchange is strongest among populations located in similar environments, which would be mediated by environmental heterogeneity, the extent of local adaptation and spatial variation in gene flow across space [18,21,22]. Thus, genetic differentiation among populations increases with their environmental differentiation, independently of their geographic distance. IBE can arise due to multiple factors, such as biased dispersal due to preferences for particular environments, natural selection against maladapted immigrants, sexual selection against immigrants when they exhibit divergence in mating choices or sexual signals, and natural selection against hybrids when they show reduced fitness relative to non-hybrids [22].
Finally, the IBR model takes environmental heterogeneity across landscapes into account as a modulator of gene flow and its effects on genetic differentiation [23]. The IBR model predicts a positive relationship between genetic differentiation and resistance distance among populations [24][25][26]. The resistance distance between population pairs is a concept inspired in circuit theory, which considers the landscape features reducing the probability of dispersal and gene flow among populations [24][25][26][27]. Under IBR, the spatial structure of habitat suitability is of paramount importance to determine the least cost path between population pairs optimizing their connectivity and thus minimizing their resistance distance (24 and references therein). Given that the resistance distance between population pairs is a function of geographic distance, and that presence-background models estimate habitat suitability using environmental predictors, IBR inevitably conflates IBD and IBE [22].
Here, we tested these three hypotheses to identify the major drivers of genetic differentiation in the annual plant Arabidopsis thaliana across the Iberian Peninsula. This is the region of the species' distribution harboring the largest genomic diversity [28][29][30]. In addition, Iberian A. thaliana occurs in a wide array of natural environments practically across the whole region, spanning from seaside to sub-alpine locations [31][32][33]. These two features are likely the result of A. thaliana's history in the Iberian Peninsula [29,31,34], where the species long survived by developing adaptations to ample environmental heterogeneity over dramatic climatic oscillations. The occurrence of relict populations with an African origin [29,34] also supports such history of Iberian A. thaliana. In fact, habitat suitability of relict populations has been associated to more stable vegetation dynamics since the Last Glacial Maximum and during the Holocene in the Iberian Peninsula [35]. Overall, the geographic ubiquity, the large amount of genetic diversity, the broad variety of habitats occupied, and the long evolutionary history make Iberian A. thaliana an appropriate study system to disentangle the drivers of genetic differentiation at a regional scale.
Given the fact that a small sample size seriously reduces power and accuracy of spatial analyses [36], IBD, IBE and IBR were tested using 278 Iberian A. thaliana populations collected over a decade. About six individuals per population, totaling 1772 individuals, were genotyped with genome-wide, putatively neutral SNPs to estimate genetic diversity, differentiation and structure. We hypothesized that IBD and IBE largely accounted for genetic differentiation for two reasons. First, IBD is a common result in A. thaliana genetic studies-including the Iberian Peninsula [31][32][33]-due to dispersal limitation and high selffertilization rates. Second, Iberian A. thaliana shows significant adaptive variation in fitness-related traits across environmental gradients [30,32,37,38]. However, we ignore the importance of IBR and the joint contribution of IBD, IBE and IBR to the genetic differentiation of Iberian A. thaliana populations. For the sake of completeness, we also examined the geographical distribution of genetic diversity in Iberian A. thaliana. To this end, we employed a spatial hierarchical Bayesian model to identify regional hot and cold spots of genetic diversity and their potential environmental predictors. Overall, we stress the importance of identifying the drivers of genetic differentiation among populations, but also of pinpointing the forces that determine the amount of genetic diversity within populations, to understand the evolutionary dynamics of any organism.

Genetic diversity
We genotyped 240 genome-wide SNPs in 1172 individuals collected from 278 populations (Fig. 1a) to determine the genetic diversity of A. thaliana in the Iberian Peninsula. In this set of populations, the percentage of polymorphic loci (PL) ranged between 0 and 56.7% (mean ± SD = 24.4 ± 17.7%), the mean number of observed alleles per locus (n a ) between 0.95 and 1.56 (mean ± SD = 1.23 ± 0.18), and mean gene diversity (H S ) between 0 and 0.224 (mean ± SD = 0.09 ± 0.07). H S exhibited a bimodal distribution clearly differentiated in two groups of populations. On the one hand, 66 populations had very low H S values (54 of 0 and 12 between 0.001 and 0.009; Fig. 1c). On the other hand, the remaining 212 populations had H S values distributed between 0.015 and 0.224 (Fig. 1c). The 66 populations with no genetic diversity did not show geographic or environmental bias, as they occurred scattered across the region (Fig. 1a) and the altitude gradient (Fig. 1b). Average (± SD) genetic differentiation, given by pairwise F ST values, including all populations was of 0.644 ± 0.189 (F ST = 0.538 ± 0.150 without the 66 populations with no genetic diversity).
Overall, we found 1613 non-redundant multilocus genotypes (N H ) in the 1772 individuals, which showed an average proportion of allelic differences between haplotype pairs of 0.30 ± 0.05. We only found five non-redundant multilocus genotypes in different populations. In particular, two (separated by 0.6 km) and three populations (separated by 7.6-12.8 km) included the same nonredundant multilocus genotypes.
To assess whether populations experienced major or sudden disturbances that could affect their genetic diversity and differentiation during the last decades, we retrieved temporal series of publicly available orthophotographs from all populations. The analyses of the habitats interpreted from digitalized orthophotographs indicated that land-used changes and major disturbances (wildfires, landslides or the development of large infrastructures) over the time considered did not affect landscapes. In particular, the average change per habitat type between year intervals was concentrated around zero ( Fig.  1d), indicating that habitats barely changed over time. We ignore, however, whether pathogen or diseases affected A. thaliana populations, which could influence genetic diversity patterns. We have never observed noticeable catastrophic events caused by pests and diseases over 15-plus years of experience sampling natural A. thaliana populations (C. Alonso-Blanco and F. X Picó, pers. obs.).

Geographic distribution of genetic diversity
To determine the geographic and environmental distribution of H S in A. thaliana populations, we used a spatial hierarchical Bayesian modeling to explain H S as a function of environmental variables. The model H S values. The best fitting model indicated that the occurrence probability of populations with non-zero H S values (populations with genetic diversity) was concentrated towards the center and eastern areas of the Iberian Peninsula (Fig. 2a). In contrast, northern, western and southern peripheral areas had higher odds for populations lacking genetic diversity (Fig. 2a). In addition, genetic diversity of A. thaliana was unevenly distributed across the Iberian Peninsula, with different nuclei with high diversity in large central and northern areas, as well as in a limited area in NE Spain (Fig. 2b). We did not detect any anomaly in the distribution of the mean effects of the spatial component and its uncertainty (Fig. S2) that made us suspect that the spatial effects were affecting the results.
The environmental predictors with more influence on the spatial distribution of H S were precipitation seasonality (BIO15), precipitation of the warmest quarter (BIO18) and topsoil pH. The contributions of precipitation variables were negative for the occurrence of populations with no genetic diversity, as well as for the genetic diversity of populations with non-zero H S values (Table 1). Overall, these results indicated that areas with higher precipitation seasonality, and to a lesser extent higher precipitation in the warmest quarter, tended to harbor populations with lower genetic diversity. In the Iberian Peninsula, higher and lower precipitation seasonality characterizes xeric and mesic environments, respectively. The contribution of topsoil pH was positive for the occurrence probability of populations with zero or non-zero H S values, meaning that populations with non-zero H S values tended to occur in areas with basic soils. However, in the basic areas of the Iberian Peninsula (E and SE Spain) the species is very rare ( Fig. 1 and S1), and the low number of populations there with nonzero H S values could be introducing some bias in the model. In contrast, the contribution of topsoil pH was negative for genetic diversity, indicating that populations with higher genetic diversity tended to occur in acidic soils, which is consistent with the preference of the species for this sort of soils (C. Alonso-Blanco and F.X. Picó, pers. obs.).

Genetic structure
We analyzed the structure of the 278 A. thaliana populations in different genetic clusters with the 1613 nonredundant multilocus genotypes using Bayesian (Fig. 3a) and ordination (Fig. 3b) methods. Both analyses detected four genetic clusters, three of them showing a strong geographic structure (Fig. 3c): the northwestern cluster 1 (NW-C1 hereafter), the northeastern cluster 2 (NE-C2 hereafter), and the southwestern cluster 4 (SW-C4 hereafter; Fig. 3b). In contrast, cluster 3, which corresponded to the group previously described as the relict A. thaliana lineage (relict-C3 hereafter), occurred scattered across the Iberian Peninsula (Fig. 3c). These results were fully consistent with previous studies on Iberian A. thaliana based on just one accession per population [31,32,34,39].
Since we used about six individuals per population, we could classify populations as homogeneous or heterogeneous based on the assignment of individuals to a single or multiple genetic clusters, respectively. Most populations (230 of 278) were homogeneous, with 118, 44, 35 and 33 belonging to genetic clusters NW-C1, NE-C2, relict-C3 and SW-C4, respectively (Fig. 3c). The remaining 48 heterogeneous populations were distributed across the region, with an accumulation of them in central and NE Spain (Fig. 3c). As expected, the most abundant genetic clusters, NW-C1 and NE-C2, exhibited

Drivers of genetic differentiation
We applied nested maximum-likelihood population effect models (NMLPE) to evaluate simultaneously the effects of the geographic (IBD), environmental (IBE) and resistance (IBR) variation on the genetic differentiation among Iberian A. thaliana populations. IBD and IBE accounted for genetic differentiation at the entire Iberian Peninsula scale (Table 2A and Fig. 4a). In addition, we dissected IBE into three environmental PCs, which substantially contributed to the IBE (Table 2A and Fig. 4b).
These PCs explained 29.3, 23.8, and 22.2% of the environmental variance, respectively, which accounted for different environmental gradients across the Iberian Peninsula. In particular, PC1 represented a gradient of temperature and precipitation, pinpointing the negative relationship between the two variables in the Iberian Peninsula. PC2 only depicted the gradient of temperature across the region. Finally, PC3 illustrated a gradient of precipitation and pH in which soils with lower pH receive higher precipitation in the Iberian Peninsula.
We found that the single model including IBD and IBE was ranked as the top model for three genetic clusters (NW-C1, NE-C2, and relict-C3) with no models differing less than 2 in AIC (Table 2B and Fig. 4a). Although both IBD and IBE had significant effects on genetic differentiation, size effects were greater for IBD than for IBE (Fig. 4a). In addition, the three PC axes made substantial contributions to the IBE exhibited by the populations for the same three genetic clusters (Table 2B and Fig. 4b). For SW-C4, a single model including IBR was ranked as the top model (Table 2B and Fig. 4a). When environmental PCs were included separately in the analyses, we found two top-ranked models exhibiting ΔAIC < 2 ( Table 2) for SW-C4. However, model averaging discarded a relevant role of PC2 on the genetic differentiation among SW-C4 populations (Fig. 4b).

Discussion
Genetic differentiation is a dynamic process because any population is constantly under the effect of ecological, genetic and evolutionary forces modifying the amount of genetic diversity within and among populations. Disentangling such forces accounting for genetic differentiation, which underlies major evolutionary processes from local adaptation to speciation, has long promoted the development of a strong theoretical and methodological framework since practically the birth of the Modern Synthesis. Here, we addressed this question by testing three hypotheses (IBD, IBE and IBR) accounting for regional-scale genetic differentiation in Iberian A. thaliana.
The main result of our study, based on nested maximum-likelihood population effect models (NMLPE), indicated that genetic differentiation in A. thaliana was mostly accounted for by IBD, and to a lesser extent by IBE (Fig. 4a). In other words, dispersal limitation, genetic drift, and to a lesser extent local adaptation, determined the distribution of A. thaliana's genetic diversity within and among populations across its Iberian range. This result extends previous studies showing a marked IBD in A. thaliana in the Iberian Peninsula [30][31][32][33] and elsewhere [28,34,[40][41][42][43][44][45][46][47][48]. It is widely accepted that IBD in A. thaliana results from the joint effects of dispersal limitation, self-fertilization, local adaptation and demographic history. Our study contributed to disentangle the effect of some of these factors, suggesting that dispersal limitation and genetic drift are probably more efficient than local adaptation in shaping genetic differentiation patterns in A. thaliana. We found additional support to this conclusion in the low number of populations located at short distances sharing the same multilocus genotypes, stressing the limited natural A. thaliana's dispersal ability [41,[47][48][49].
In contrast, IBR did not account for genetic differentiation at a regional scale (Fig. 4a), suggesting that landscape heterogeneity was not relevant for genetic differentiation in Iberian A. thaliana. We believe that the lack of IBR found across the Iberian Peninsula is likely due to the species' cosmopolitan habit and a distribution without major discontinuities (Fig. S1). Thus, most population pairs are connected by relatively high habitat suitability, diminishing resistance distances. The demographic history, characterized by the presence of the relict African lineage and other non-relict Eurasian lineages [28-31, 34, 35, 50], and the well-known ability to adjust flowering time and seed dormancy to contrasting environments [38,[51][52][53][54][55], both account for the cosmopolitan habit of Iberian A. thaliana.
We also explored the effect of IBD, IBE and IBR for each of the four genetic clusters, which represent the joint outcome of the long demographic and adaptive history of the species in the region [30-32, 34, 39, 56]. Overall, NW-C1, NE-C2 and relict-C3 exhibited the same patterns than   those at the regional scale, i.e. IBD and IBE mostly accounted for genetic differentiation within each cluster (Fig. 4a). However, it is worth noting that the importance of IBE was not the same for these three clusters, since NE-C2 and relict-C3 exhibited higher IBE than NW-C1 (Fig. 4a), stressing probably the higher influence of local adaptation for genetic differentiation in these clusters. In fact, NE-C2 is a cluster with pronounced altitudinal gradients in NE Spain (from maritime to sub-alpine environments) along which the species is known to have adapted by adjusting various life-cycle and physiological traits [57][58][59][60][61]. In the case of relict-C3, the long evolutionary history of this lineage has provided the means to survive and adapt to diverse Iberian environments over the last millennia [35]. On the contrary, that was not the case for SW-C4, whose genetic differentiation was accounted for by IBR (Fig. 4a). Whereas SW-C4 is the genetic cluster with the most restricted geographic distribution (Fig. 3c), its populations are located across an area showing dramatic changes in habitat suitability (Fig. S1). Most likely, this particularity accounted for the enormous weight of IBR in this cluster. Finally, we also conducted NMLPE using principal component (PC) axes to pinpoint the major environmental predictors underlying IBE. We found that the three PCs depicting gradients of temperature and/or precipitation across the region made substantial contributions to IBE in NW-C1, NE-C2 and relict-C3 Iberian clusters (Fig. 4b), probably due to the similar contributions of the three PCs to the environmental variation among these populations. In addition, as IBE informs on the influence of local adaptation on genetic differentiation, our results reinforce the view of the strong effects of environmental gradients on local adaptation in Iberian A. thaliana. For example, previous studies showed that geographic variation in fitnessrelated traits in Iberian A. thaliana (flowering time and seed dormancy) strongly co-varied as a function of temperature and to a lesser extent of precipitation [30,32,33,38,55,62,63]. In particular, A. thaliana adjusts its life cycle by advancing flowering time and increasing seed dormancy as the environment becomes warmer, drier and more seasonal [38,55,63].
As genetic diversity and its within-and amongpopulation distribution patterns represents the raw material upon which genetic differentiation is estimated, we discuss some results dealing with genetic diversity of Iberian A. thaliana that are worth considering. For example, analyses of the genetic diversity in this large number of populations showed that 24% of them had no or very low genetic diversity (Fig. 1a-c). This is not exceptional as other authors also detected populations with practically no genetic diversity in A. thaliana [40,41,49,64,65]. Strong founder effects [66] and low migration rates [41] could well account for this result. Furthermore, small A. thaliana populations rather isolated and locally adapted to particular environments can have very low genetic diversity [67], a scenario that could also apply to some of our populations. Finally, sampling bias might also affect genetic diversity [68], particularly because some populations were sampled late in the reproductive season.
Whatever the probable combination of factors accounting for the existence of A. thaliana populations with no or very low genetic diversity, the spatial hierarchical Bayesian model allowed the analysis of the geographic distribution of genetic diversity in Iberian A. thaliana (Fig. 2), particularly to detect the environmental predictors of the geographic distribution of genetic diversity. Precipitation, and not temperature, emerged as the major environmental factor accounting for the split of populations with and without genetic diversity, as well as for genetic diversity (Table 1). This result is in agreement with other studies on the distribution of genetic diversity in plants, which indicated a trend for a higher role of precipitation variables over those of temperature [69][70][71][72][73][74]. In addition, precipitation seasonality, which determines climatic variation in the Iberian Peninsula spanning from drier Mediterranean to more humid Atlantic climates, was the most important bioclimatic variable for the geographic distribution of genetic diversity in A. thaliana (Table 1). Previous studies also found that precipitation seasonality was a good predictor of individual fitness [33] and distribution of the four genetic clusters [39]. Albeit precipitation seasonality may be hard to handle as a fixed factor, further experiments are needed to better understand how A. thaliana responds to this variable, which might be a major evolutionary force in this species.

Conclusions
This study dissected the complexity of geographic, environmental and evolutionary factors contributing to the genetic differentiation among A. thaliana populations, while illustrating the power of dense collections to disentangle complex biological questions [8,28,30,49,63,75,76]. Beyond genetic differentiation patterns, we also identified the environmental predictors accounting for genetic diversity within populations. Clearly, the processes driving among-population genetic differentiation may not be the same than those determining the amount of within-population genetic diversity. However, we need to deal with the latter to understand the former. In this sense, recent resurrection studies showed that natural A. thaliana populations exhibited substantial changes in their genetic composition in just a decade [77,78]. The extent of temporal genetic variation in A. thaliana represents a reminder that populations are not static, and that the geographic distribution of genetic diversity within and among-populations is changing too. Repeated genomic scans over time on the same populations will provide real-time insights into the intensity and pace of genetic variation within populations, which will increase our understanding of the evolutionary processes shaping genetic differentiation in plant populations.
Natural A. thaliana populations are made of patches of individuals widely differing in size and density, which is important when it comes to design a sampling scheme to study the spatio-temporal distribution of genetic diversity in A. thaliana [49,57,65,77]. For each population and whenever possible, we collected seeds from several individuals from different patches (separated 1-20 m from each other) representing well each study population. Every sampling year and a few months after sampling, we multiplied field-collected seeds by the single seed descent method, in a glasshouse at the Centro Nacional de Biotecnología (CNB-CSIC) in Madrid. Bulked seeds were stored in dry, dark conditions in cellophane bags at room temperature, storing conditions that can preserve A. thaliana seeds for years.
In this study, each population included about six individuals, ranging between four and seven. Given the importance to work with a dense collection of natural populations, we selected this number as a compromise between the number of populations and the number of individuals per population that could be handled to extract the genetic patterns of interest. This choice was also based on a previous large-scale study exploring the within-and among-population partitioning of genetic diversity in A. thaliana, which yielded interpretable results with four individuals per population [40]. Besides, other studies using more individuals per population and different markers, albeit at smaller geographical scales, also found populations with no or very low genetic diversity (see below), indicating that is not unusual in A. thaliana [40,41,49,64,65,77].
We collected all seeds from wild populations. Arabidopsis thaliana is a common plant species not categorized as protected or endangered in any species list of the Convention on the Trade in Endangered Species of Wild Fauna and Flora. We carried out field sampling in locations where no permission was required, except at Doñana National Park (permission issued by Estación Biológica de Doñana) and Sierra de Grazalema Natural Park (permission issued by Red Andaluza de Jardines Botánicos de la Consejería de Medio Ambiente y Ordenación del Territorio de la Junta de Andalucía). Based on the Royal Decree of the Spanish legislation (Real Decreto 124/2017, de 24 de febrero; https://www.boe.es/ eli/es/rd/2017/02/24/124), the genetic resources included in this study fall within the definition of "exclusively for taxonomic purposes" as they were used for scientific, educational or non-commercial purposes.

Historical factors
Recent historical changes in land use or large-scale perturbations may dramatically affect the amount of genetic diversity in plant populations. Such changes can be a misleading factor when attempting to evaluate the joint effects of IBD, IBE and IBR on genetic differentiation. This is because IBD, IBE and IBR do not contemplate perturbations producing sudden or erratic changes in genetic diversity that eventually may mask the patterns of interest. To make sure that historical factors were not an issue in this study, we estimated the temporal landscape dynamics of all study populations. We examined all aerial orthophotographs available for each population (spanning between 1945 and 2016), which were retrieved from different public regional administrations in Spain [77] with a tool specifically developed for this purpose. We excluded up to 15 populations because aerial orthophotographs were not available. The final number of aerial orthophotographs was highly variable per population, ranging between 4 and 24 (mean ± SE = 8.6 ± 0.3 aerial orthophotographs) and covering between 9 and 71 years (37.1 ± 1.7 years between the first and the last orthophotograph).
To quantify temporal landscape dynamics, a circular area (500 m radius) around the GPS coordinate was divided into a regular grid of 80 squares (100 m side) for each aerial orthophotograph available per population. Based on vegetation and land use, each square was categorized as forests, scrublands, dehesa (i.e. agro-silvicultural ecosystems based on Mediterranean oak woodlands), bare soil, crops, urban, infrastructures, and water, when one of these categories occupied at least more than half of the square. Categories were assigned to each square by digitalizing manually all squares from all aerial orthophotographs available (N = 2252 orthophotographs) with QGIS v.3.4 [80]. For each habitat type and population, we estimated the average change between year intervals available for each population to quantify landscape changes over time.

SNP genotyping
A total of 1772 A. thaliana individuals from 278 populations were genotyped with 245 presumably neutral nuclear SNPs using the SNPlex technique (Applied Biosystems, Foster City, CA, USA) through the CEGEN Genotyping Service (www.usc.es/cegen/). These genome-wide SNPs are frequent polymorphisms in Central Europe, the Iberian Peninsula, and worldwide collections, which altogether minimized ascertainment bias [31,32,34]. On average, there were about 49 SNPs per chromosome (range = 45-53 SNPs) located at approximately 0.5 Mb from each other (range = 0.11 Kb -1.82 Mb). Four SNPs had percentages of missing data above 25% and one was monomorphic. We discarded these five SNPs from the analyses maintaining 240 SNPs.

Genetic analyses
For each A. thaliana population, we calculated the percentage of polymorphic loci (PL), the mean number of observed alleles per locus (n a ) and mean gene diversity (H S ) using FSTAT v.2.9.3 [81]. In addition, we computed the percentage of differences among all pairs of nonredundant multilocus genotypes. We used the 'pairwise.WCfst' function implemented in the R package hierfstat v.0.04-22 [82] to estimate genetic differentiation as pairwise F ST values according to Weir and Cockerham [68], which represents the dependent variable to test IBD, IBE and IBR (see below).
The genetic structure of A. thaliana populations across the Iberian Peninsula was examined using two different methods: the Bayesian clustering method implemented in STRUCTURE v.2.3.3 [83,84] and the ordination method represented by the Discriminant Analysis of Principal Components (DAPC) [85] available in the R package adegenet v.2.1.1 [85]. In the case of STRUCTURE, model settings included haploid non-redundant multilocus genotypes, correlated allele frequencies between populations and a linkage model. We identified the number of identical multilocus genotypes using the 'mlg.filter' function implemented in adegenet. Each run consisted of 50,000 burn-in MCMC iterations and 100,000 MCMC afterburning repetitions for parameter estimation. To determine the K number of ancestral populations and the ancestry membership proportions of each accession in each population, we ran the algorithm 20 times for each defined number of groups (K value) from 1 to 10. The number of distinct genetic clusters was determined by evaluating the differences between the data likelihood for successive K values. The largest K value with significantly higher likelihood than that of K-1 runs (two-sided P < 0.005; Wilcoxon tests for related samples) gave the final K number. This was supported by a high similarity among the ancestry membership matrices from different runs of the same K value (H′ = 0.99). We used CLUMPP v.1 [86] to calculate the average symmetric similarity coefficient H ′ among runs and the average matrix of ancestry membership proportions, derived from the 10 runs having the highest likelihood. STRUCTURE simulations were conducted at the The Supercomputing Center of Galicia (CESGA; http://www.cesga.es/).
For the DAPC analysis, we assessed the number of clusters using the 'find.cluster' function in adegenet, which runs successive K-means clustering with an increasing number of clusters to determine the best supported number of genetic clusters using the Bayesian Information Criterion (BIC). The K value with the lowest BIC value represents the optimal number of clusters, although BIC values may keep decreasing after the true K value in case of genetic clines and hierarchical structure [87]. Therefore, the rate of decrease in BIC values was visually examined to identify values of K, after which BIC values only decreased in a subtle manner [87]. The 'dapc' function was used with the final K value, retaining the axes of the Principal Component Analysis.

Geographic distribution of genetic diversity
We modeled the effects of environmental heterogeneity on the spatial distribution of genetic diversity by using a modified version of a spatial hierarchical Bayesian model recently used to model the effects of warming on the Iberian A. thaliana's distribution [56]. The aim of this model was to visualize hot and cold spots of genetic diversity across the Iberian Peninsula and to identify their environmental predictors.
We used a Bayesian approach to handle the particularities of genetic diversity data, such as the semi-continuous nature of the variable. In the case of Iberian A. thaliana, 66 of 278 populations had very low mean gene diversity (H S ) values (H S < 0.009), whereas the rest of populations exhibited non-zero H S values (Fig. 1c). To handle such complex distribution of H S , we developed a Hurdle-beta Bayesian model [88] to separate the zero structure (absence of genetic diversity) from the non-zero structure (presence of genetic diversity) of data. The Hurdle model is defined as a finite mixture of two processes: a degenerate distribution with point mass and zero, which determines the presence or absence of data-in our case H S values above or below the value of 0.009, respectivelyand a conditional-to-presence continuous process supported by the open interval between the lowest non-zero H S value and 1. The first process (absence of genetic diversity) was modeled with a Bernoulli distribution, whilst the continuous component (presence of genetic diversity) was modeled with a Beta distribution (see the full development of the model in the Supporting Information section). The predictors of the model were the 19 bioclimatic variables available from WorldClim and topsoil pH available from The Geographical Database from Eurasia. The Watanabe-Akaike information criterion (WAIC) [89]; was computed to determine the best models. The model included a stochastic spatial effect (the spatial term) to remove the spatial autocorrelation of data [56].

Drivers of genetic differentiation
We estimated the ecological, genetic and evolutionary drivers of genetic differentiation of Iberian A. thaliana as follows. IBD was based on geographic distance. However, instead of using the Euclidian distance among population pairs, we created a raster by assigning a 0.5 value to all 1-km 2 pixels [90,91] of the Iberian Peninsula map. We then calculated pairwise distances, known as geographic resistance distance, among all A. thaliana populations employing the new raster. We did not use the Euclidean distance to compute IBD because straight lines between population pairs may introduce some bias when such straight lines traverses natural barriers for the study organism (i.e. the coastline or open sea) [92].
IBE was based on the 19 bioclimatic variables. In addition, we used The CORINE Land Cover to obtain the percentage of human-modified habitat (i.e. urban areas, crops and semi-natural grasslands) within a 500 m radius from the population GPS coordinate [39]. The percentages of human-modified and natural habitat were significantly negatively correlated (N = 278; r = − 0.95; P < 0.001; Pearson's correlation) and nearly summed to 100%. Finally, a topsoil pH layer was also used. Combining these 21 environmental variables, we conducted a principal component analysis (PCA) with varimax rotation using the software SPSS v.25 (IBM, Chicago, IL USA). Thus, we obtained the PC scores of the first three principal components (PC) for each population. We then used the PC scores to calculate environmental dissimilarity between populations as Euclidean distances derived with the 'dist' function in R. Furthermore, we also calculated dissimilarity matrices among populations for each PC axis separately, as they were interpretable in terms of environmental gradients.
To compute IBR, we first estimated habitat suitability for A. thaliana across the Iberian Peninsula with the maximum-entropy modeling technique implemented in the software Maxent 3.3.3 k [93,94]. To do that, we used 478 A. thaliana populations available to date (collected between 2000 and 2019; Fig. S1), and the same environmental variables used in a previous habitat suitability model [39]. In particular, the model included topsoil pH and eight bioclimatic variables, which were annual mean temperature (BIO1), mean diurnal temperature range (BIO2), isothermality (BIO3), temperature seasonality (BIO4), mean temperature of wettest quarter (BIO8), annual precipitation (BIO12), precipitation seasonality (BIO15) and precipitation of the warmest quarter (BIO18). We applied Maxent using default parameters except for features using the hinge type, making it comparable to a Generalized Additive Model [94]. We obtained consistent results with previous habitat suitability models conducted for Iberian A. thaliana with lower sample sizes [39,56].
Next, we transformed the resulting habitat suitability into resistance values as 1 minus the value of each pixel. Thus, greater pixel values represented lower occurrence probability of A. thaliana (higher resistance) and viceversa. Using circuit theory [24,25], we examined whether genetic differentiation among A. thaliana populations was accounted for by IBR. We used Circuitscape v.4.0 [95] to calculate resistance distance matrices assigning pixel values as resistance values. We calculated resistance distance among all pairs of populations across Iberian Peninsula, as well as among all pairs of populations within each genetic cluster.
To quantify the relative contributions of IBD, IBE and IBR on the genetic differentiation (pairwise F ST values) among A. thaliana populations, we fitted maximumlikelihood population effect models (MLPE) [96]. Code implementing the MLPE correlation structure within the R package nlme v. 3.1-143 comes from the corMLPE package. The model used penalized least squares and a residual covariance structure designed to account for the non-independence of pairwise distances. Given the inherent spatial dependence structure in pairwise comparisons between A. thaliana populations, all MLPE considered spatial locations. To this end, we used a modification of the MLPE model incorporating the correlation between pairwise measurements due to comparison of populations and spatial locations (nested MLPE or NMLPE) [91]. Since the independent variables (IBD, IBE and IBR) used in NMLPE had very different units, we normalized them by subtracting the mean and dividing by the standard deviation in all analyses. We did not linearize pairwise F ST values. We conducted NMLPE on the complete dataset (N = 278 populations) as well as on the subset of those populations exhibiting H S different from zero (N = 212 populations with H S > 0.009). The results were consistent between the two datasets, showing that the 66 populations with no genetic diversity were not biasing the results.
In addition, we ran NMLPE for each genetic cluster to evaluate whether the effects of IBD, IBE and IBR on genetic differentiation differed among genetic clusters. This is because the number of genetic clusters represents the outcome of the demographic and evolutionary history of A. thaliana in the region. We considered each population to belong to a specific genetic cluster when the average value of the membership proportions of its individuals was ≥0.3 for that specific cluster. Given the marked genetic structure of the study system, threshold values of 0.25 and 0.40 were consistent, but 0.3 provided a clearer classification of individuals into genetic clusters. The majority of populations were assigned to a single cluster (N = 230). We assigned populations reaching the minimum average membership proportion for more than one cluster to multiple clusters. Overall, 47 populations were assigned to two clusters, whereas only one population had average membership proportions higher than 0.3 for three clusters.
We also ran NMLPE on the complete dataset and each cluster separately, but replacing IBE by the PC axes. Thus, we used IBD, PC1, PC2, PC3 and IBR as independent variables. We used the Akaike information criterion (AIC) to compare models containing all possible combinations of non-collinear predictions (r < 0.6), created with the 'dredge' function from the R package MuMIn v1.4 [97]. In those cases where several models were top-ranked with ΔAIC < 2, we calculated the model averaged parameter estimates (β) with standard errors for the explanatory variables included in all top-ranked models also using MuMIn. We considered effect sizes as significant when 95% confidence intervals did not overlap zero.