Castilla et al. BMC Evolutionary Biology (2020) 20:71 https://doi.org/10.1186/s12862-020-01635-2 RESEARCH ARTICLE Open Access Ecological, genetic and evolutionary drivers of regional genetic differentiation in Arabidopsis thaliana Antonio R. Castilla1,2, Belén Méndez-Vigo3, Arnald Marcer4,5, Joaquín Martínez-Minaya6, David Conesa7, F. Xavier Picó2* and Carlos Alonso-Blanco3 Abstract 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. Keywords: Genetic diversity, Genetic structure, Iberian Peninsula, Nested maximum-likelihood population effect models, Precipitation seasonality, Spatial hierarchical Bayesian models Background causes, pace and the phenotypic effects of such changes— Genetic diversity is an important asset for the ability of constitute the raw material upon which natural selection populations to adapt and persist in changing environments eventually acts [9–11]. At any spatial scale, genetic diversity [1–8]. The spatio-temporal changes in genetic diversity con- typically becomes unevenly distributed across space [12], as stantly taking place in any population—regardless of the genetic diversity is determined by how genetic diversity is partitioned within and among populations across the distri- bution. In other words, the spatial distribution of genetic * Correspondence: xpico@ebd.csic.es 2 diversity depends on the extent of genetic differentiation Departamento de Ecología Integrativa, Estación Biológica de Doñana (EBD), Consejo Superior de Investigaciones Científicas (CSIC), Sevilla, Spain among populations whatever the sources of such Full list of author information is available at the end of the article © The Author(s). 2020 Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. Castilla et al. BMC Evolutionary Biology (2020) 20:71 Page 2 of 13 differentiation. Nonetheless, it must be noted that genetic therein). Given that the resistance distance between popu- differentiation is a spatially-explicit phenomenon, as genetic lation pairs is a function of geographic distance, and that differentiation strictly depends on the genetic relationship presence-background models estimate habitat suitability that a given population has with its neighbors near and far. using environmental predictors, IBR inevitably conflates The inherent spatial nature of genetic differentiation IBD and IBE [22]. defines the theoretical and methodological framework of Here, we tested these three hypotheses to identify the three models, which are not mutually exclusive, testing major drivers of genetic differentiation in the annual the major drivers of genetic differentiation: isolation-by- plant Arabidopsis thaliana across the Iberian Peninsula. distance (IBD hereafter), isolation-by-environment (IBE This is the region of the species’ distribution harboring hereafter) and isolation-by-resistance (IBR hereafter) the largest genomic diversity [28–30]. In addition, Iber- models. In the classical IBD, genetic differentiation ian A. thaliana occurs in a wide array of natural envi- among populations exhibits a positive relationship with ronments practically across the whole region, spanning geographic distance [13–18]. In this case, dispersal limi- from seaside to sub-alpine locations [31–33]. These two tation and genetic drift determine the greater genetic features are likely the result of A. thaliana’s history in differentiation at larger distances. In fact, limited disper- the Iberian Peninsula [29, 31, 34], where the species long sal constrains gene flow among populations, which is survived by developing adaptations to ample environ- not able to counteract the effect of genetic drift within mental heterogeneity over dramatic climatic oscillations. populations [16–18]. Many types of organisms exhibit The occurrence of relict populations with an African ori- IBD [19], probably because unrestricted gene flow hardly gin [29, 34] also supports such history of Iberian A. occurs in nature. Nevertheless, quantifying the contribu- thaliana. In fact, habitat suitability of relict populations tion of limited dispersal and genetic drift to genetic dif- has been associated to more stable vegetation dynamics ferentiation is not a straightforward task, as we largely since the Last Glacial Maximum and during the Holo- ignore the actual extent of dispersal, the effective popu- cene in the Iberian Peninsula [35]. Overall, the geo- lation sizes conditioning genetic drift, and the effects of graphic ubiquity, the large amount of genetic diversity, historical factors shaping the genetic relationships the broad variety of habitats occupied, and the long evo- among populations [20]. lutionary history make Iberian A. thaliana an appropri- In contrast, the IBE model deals with the effects of envir- ate study system to disentangle the drivers of genetic onmental differences on genetic differentiation. IBE posits differentiation at a regional scale. that gene exchange is strongest among populations located Given the fact that a small sample size seriously reduces in similar environments, which would be mediated by en- power and accuracy of spatial analyses [36], IBD, IBE and vironmental heterogeneity, the extent of local adaptation IBR were tested using 278 Iberian A. thaliana populations and spatial variation in gene flow across space [18, 21, 22]. collected over a decade. About six individuals per popula- Thus, genetic differentiation among populations increases tion, totaling 1772 individuals, were genotyped with with their environmental differentiation, independently of genome-wide, putatively neutral SNPs to estimate genetic their geographic distance. IBE can arise due to multiple fac- diversity, differentiation and structure. We hypothesized tors, such as biased dispersal due to preferences for particu- that IBD and IBE largely accounted for genetic differenti- lar environments, natural selection against maladapted ation for two reasons. First, IBD is a common result in A. immigrants, sexual selection against immigrants when they thaliana genetic studies—including the Iberian Peninsula exhibit divergence in mating choices or sexual signals, and [31–33]—due to dispersal limitation and high self- natural selection against hybrids when they show reduced fertilization rates. Second, Iberian A. thaliana shows sig- fitness relative to non-hybrids [22]. nificant adaptive variation in fitness-related traits across Finally, the IBR model takes environmental heterogen- environmental gradients [30, 32, 37, 38]. However, we ig- eity across landscapes into account as a modulator of gene nore the importance of IBR and the joint contribution of flow and its effects on genetic differentiation [23]. The IBD, IBE and IBR to the genetic differentiation of Iberian IBR model predicts a positive relationship between genetic A. thaliana populations. For the sake of completeness, we differentiation and resistance distance among populations also examined the geographical distribution of genetic di- [24–26]. The resistance distance between population pairs versity in Iberian A. thaliana. To this end, we employed a is a concept inspired in circuit theory, which considers the spatial hierarchical Bayesian model to identify regional hot landscape features reducing the probability of dispersal and cold spots of genetic diversity and their potential en- and gene flow among populations [24–27]. Under IBR, vironmental predictors. Overall, we stress the importance the spatial structure of habitat suitability is of paramount of identifying the drivers of genetic differentiation among importance to determine the least cost path between populations, but also of pinpointing the forces that deter- population pairs optimizing their connectivity and thus mine the amount of genetic diversity within populations, minimizing their resistance distance (24 and references to understand the evolutionary dynamics of any organism. Castilla et al. BMC Evolutionary Biology (2020) 20:71 Page 3 of 13 Results particular, two (separated by 0.6 km) and three popula- Genetic diversity tions (separated by 7.6–12.8 km) included the same non- We genotyped 240 genome-wide SNPs in 1172 individuals redundant multilocus genotypes. collected from 278 populations (Fig. 1a) to determine the To assess whether populations experienced major or genetic diversity of A. thaliana in the Iberian Peninsula. In sudden disturbances that could affect their genetic diver- this set of populations, the percentage of polymorphic loci sity and differentiation during the last decades, we re- (PL) ranged between 0 and 56.7% (mean ± SD = 24.4 ± trieved temporal series of publicly available 17.7%), the mean number of observed alleles per locus (na) orthophotographs from all populations. The analyses of between 0.95 and 1.56 (mean ± SD = 1.23 ± 0.18), and mean the habitats interpreted from digitalized orthophotographs gene diversity (HS) between 0 and 0.224 (mean ± SD = indicated that land-used changes and major disturbances 0.09 ± 0.07). HS exhibited a bimodal distribution clearly dif- (wildfires, landslides or the development of large infra- ferentiated in two groups of populations. On the one hand, structures) over the time considered did not affect land- 66 populations had very low HS values (54 of 0 and 12 be- scapes. In particular, the average change per habitat type tween 0.001 and 0.009; Fig. 1c). On the other hand, the between year intervals was concentrated around zero (Fig. remaining 212 populations had HS values distributed be- 1d), indicating that habitats barely changed over time. We tween 0.015 and 0.224 (Fig. 1c). The 66 populations with ignore, however, whether pathogen or diseases affected A. no genetic diversity did not show geographic or environ- thaliana populations, which could influence genetic diver- mental bias, as they occurred scattered across the region sity patterns. We have never observed noticeable cata- (Fig. 1a) and the altitude gradient (Fig. 1b). Average (± SD) strophic events caused by pests and diseases over 15-plus genetic differentiation, given by pairwise FST values, includ- years of experience sampling natural A. thaliana popula- ing all populations was of 0.644 ± 0.189 (FST = 0.538 ± 0.150 tions (C. Alonso-Blanco and F. X Picó, pers. obs.). without the 66 populations with no genetic diversity). Overall, we found 1613 non-redundant multilocus geno- Geographic distribution of genetic diversity types (NH) in the 1772 individuals, which showed an aver- To determine the geographic and environmental distri- age proportion of allelic differences between haplotype bution of HS in A. thaliana populations, we used a pairs of 0.30 ± 0.05. We only found five non-redundant spatial hierarchical Bayesian modeling to explain HS as a multilocus genotypes in different populations. In function of environmental variables. The model Fig. 1 Distribution of populations, genetic diversity and temporal habitat changes of Iberian Arabidopsis thaliana populations. a Geographic distribution of the 278 A. thaliana populations of study across the Iberian Peninsula. Red and blue dots represent populations with zero and non-zero genetic diversity (HS) values, respectively. b Frequency distribution of populations with zero and non-zero HS values as a function of altitude. Mean altitude for both groups of populations is almost identical, as indicated by dashed lines. c Frequency distribution of populations with zero and non- zero HS values. d Frequency distribution of the average percentage change between year intervals for each habitat type. Data from digitalized orthophotographs available from each population. For the sake of clarity, only one X-axis is shown, indicating the accumulation of populations with average percentage changes around zero. The map of Fig. 1a was obtained from the National Center for Geographic Information (CNIG) of Spain Castilla et al. BMC Evolutionary Biology (2020) 20:71 Page 4 of 13 considered simultaneously the populations with zero S1), and the low number of populations there with non- (N = 66 populations) and non-zero (N = 212 populations) zero HS values could be introducing some bias in the HS values. The best fitting model indicated that the oc- model. In contrast, the contribution of topsoil pH was currence probability of populations with non-zero HS negative for genetic diversity, indicating that populations values (populations with genetic diversity) was concen- with higher genetic diversity tended to occur in acidic trated towards the center and eastern areas of the Iber- soils, which is consistent with the preference of the spe- ian Peninsula (Fig. 2a). In contrast, northern, western cies for this sort of soils (C. Alonso-Blanco and F.X. and southern peripheral areas had higher odds for popu- Picó, pers. obs.). lations lacking genetic diversity (Fig. 2a). In addition, genetic diversity of A. thaliana was unevenly distributed Genetic structure across the Iberian Peninsula, with different nuclei with We analyzed the structure of the 278 A. thaliana popu- high diversity in large central and northern areas, as well lations in different genetic clusters with the 1613 non- as in a limited area in NE Spain (Fig. 2b). We did not de- redundant multilocus genotypes using Bayesian (Fig. 3a) tect any anomaly in the distribution of the mean effects and ordination (Fig. 3b) methods. Both analyses detected of the spatial component and its uncertainty (Fig. S2) four genetic clusters, three of them showing a strong that made us suspect that the spatial effects were affect- geographic structure (Fig. 3c): the northwestern cluster 1 ing the results. (NW-C1 hereafter), the northeastern cluster 2 (NE-C2 The environmental predictors with more influence on hereafter), and the southwestern cluster 4 (SW-C4 here- the spatial distribution of HS were precipitation season- after; Fig. 3b). In contrast, cluster 3, which corresponded ality (BIO15), precipitation of the warmest quarter to the group previously described as the relict A. thali- (BIO18) and topsoil pH. The contributions of precipita- ana lineage (relict-C3 hereafter), occurred scattered tion variables were negative for the occurrence of popu- across the Iberian Peninsula (Fig. 3c). These results were lations with no genetic diversity, as well as for the fully consistent with previous studies on Iberian A. thali- genetic diversity of populations with non-zero HS values ana based on just one accession per population [31, 32, (Table 1). Overall, these results indicated that areas with 34, 39]. higher precipitation seasonality, and to a lesser extent Since we used about six individuals per population, we higher precipitation in the warmest quarter, tended to could classify populations as homogeneous or heteroge- harbor populations with lower genetic diversity. In the neous based on the assignment of individuals to a single Iberian Peninsula, higher and lower precipitation season- or multiple genetic clusters, respectively. Most popula- ality characterizes xeric and mesic environments, re- tions (230 of 278) were homogeneous, with 118, 44, 35 spectively. The contribution of topsoil pH was positive and 33 belonging to genetic clusters NW-C1, NE-C2, for the occurrence probability of populations with zero relict-C3 and SW-C4, respectively (Fig. 3c). The or non-zero HS values, meaning that populations with remaining 48 heterogeneous populations were distrib- non-zero HS values tended to occur in areas with basic uted across the region, with an accumulation of them in soils. However, in the basic areas of the Iberian Penin- central and NE Spain (Fig. 3c). As expected, the most sula (E and SE Spain) the species is very rare (Fig. 1 and abundant genetic clusters, NW-C1 and NE-C2, exhibited Fig. 2 Geographic distribution of genetic diversity within Iberian Arabidopsis thaliana populations estimated by spatial hierarchical Bayesian modeling. a Distribution of populations with zero (N = 66; absence of genetic diversity) and non-zero (N = 212; presence of genetic diversity) genetic diversity (HS) values. Darker and lighter intensities indicate higher and lower odds for populations with non-zero and zero HS values, respectively. b Distribution of populations with HS values higher than zero (N = 212). Darker and lighter intensities indicate higher and lower HS values, respectively. In both cases, the uncertainty of spatial hierarchical Bayesian model is given as standard deviation units in small maps. Darker and lighter intensities indicate higher and lower uncertainty, respectively Castilla et al. BMC Evolutionary Biology (2020) 20:71 Page 5 of 13 a higher number of heterogeneous populations, because Fig. 4a). When environmental PCs were included separ- 46 of 48 heterogeneous populations were included in ately in the analyses, we found two top-ranked models these genetic clusters. exhibiting ΔAIC < 2 (Table 2) for SW-C4. However, model averaging discarded a relevant role of PC2 on the Drivers of genetic differentiation genetic differentiation among SW-C4 populations We applied nested maximum-likelihood population ef- (Fig. 4b). fect models (NMLPE) to evaluate simultaneously the ef- fects of the geographic (IBD), environmental (IBE) and Discussion resistance (IBR) variation on the genetic differentiation Genetic differentiation is a dynamic process because any among Iberian A. thaliana populations. IBD and IBE population is constantly under the effect of ecological, accounted for genetic differentiation at the entire Iberian genetic and evolutionary forces modifying the amount of Peninsula scale (Table 2A and Fig. 4a). In addition, we genetic diversity within and among populations. Disen- dissected IBE into three environmental PCs, which sub- tangling such forces accounting for genetic differenti- stantially contributed to the IBE (Table 2A and Fig. 4b). ation, which underlies major evolutionary processes These PCs explained 29.3, 23.8, and 22.2% of the envi- from local adaptation to speciation, has long promoted ronmental variance, respectively, which accounted for the development of a strong theoretical and methodo- different environmental gradients across the Iberian logical framework since practically the birth of the Peninsula. In particular, PC1 represented a gradient of Modern Synthesis. Here, we addressed this question by temperature and precipitation, pinpointing the negative testing three hypotheses (IBD, IBE and IBR) accounting relationship between the two variables in the Iberian for regional-scale genetic differentiation in Iberian Peninsula. PC2 only depicted the gradient of A. thaliana. temperature across the region. Finally, PC3 illustrated a The main result of our study, based on nested gradient of precipitation and pH in which soils with maximum-likelihood population effect models lower pH receive higher precipitation in the Iberian (NMLPE), indicated that genetic differentiation in A. Peninsula. thaliana was mostly accounted for by IBD, and to a We found that the single model including IBD and lesser extent by IBE (Fig. 4a). In other words, dispersal IBE was ranked as the top model for three genetic clus- limitation, genetic drift, and to a lesser extent local adap- ters (NW-C1, NE-C2, and relict-C3) with no models dif- tation, determined the distribution of A. thaliana’s gen- fering less than 2 in AIC (Table 2B and Fig. 4a). etic diversity within and among populations across its Although both IBD and IBE had significant effects on Iberian range. This result extends previous studies show- genetic differentiation, size effects were greater for IBD ing a marked IBD in A. thaliana in the Iberian Peninsula than for IBE (Fig. 4a). In addition, the three PC axes [30–33] and elsewhere [28, 34, 40–48]. It is widely ac- made substantial contributions to the IBE exhibited by cepted that IBD in A. thaliana results from the joint ef- the populations for the same three genetic clusters fects of dispersal limitation, self-fertilization, local (Table 2B and Fig. 4b). For SW-C4, a single model in- adaptation and demographic history. Our study contrib- cluding IBR was ranked as the top model (Table 2B and uted to disentangle the effect of some of these factors, Fig. 3 Genetic structure of the 278 Arabidopsis thaliana populations of study across the Iberian Peninsula depicting the four genetic clusters (NW- C1, NE-C2, relict-C3 and SW-C4). a Results from the Bayesian clustering method implemented in STRUCTURE. b Results from the Discriminant Analysis of Principal Components (DAPC). c Geographic distribution of homogeneous populations from each genetic cluster. Homogeneous populations (N = 230) were those with average membership proportions among individuals within populations greater than 0.3 for only one genetic cluster. Mixed or heterogeneous populations (N = 48) are also shown in grey. The map of Fig. 3c was obtained from the National Center for Geographic Information (CNIG) of Spain Castilla et al. BMC Evolutionary Biology (2020) 20:71 Page 6 of 13 Table 1 Coefficients of the best spatial hierarchical Bayesian model for the geographic distribution of genetic diversity of Iberian Arabidopsis thaliana HS ≈ 0 HS > 0 Variables β SD β SD Intercept 2.054 1.776 0.305 0.604 BIO15 −0.025 0.017 −0.023 0.006 BIO18 −0.009 0.004 −0.004 0.001 Topsoil pH 0.176 0.200 −0.187 0.067 Entries (β ± SD) are given for the degenerate distribution with point mass and zero (binary data including populations with HS values ≈ 0) and a conditional- to-presence continuous process (continuous data including populations with non-zero HS values). The best model had a WAIC value of − 365.43 and a LCPO value of − 0.31. BIO15 is the precipitation seasonality and BIO18 is the precipitation of the warmest quarter 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–49]. In contrast, IBR did not account for genetic differenti- ation at a regional scale (Fig. 4a), suggesting that land- scape heterogeneity was not relevant for genetic differentiation in Iberian A. thaliana. We believe that Fig. 4 Effects of isolation-by-distance (IBD), isolation-by-environment the lack of IBR found across the Iberian Peninsula is (IBE) and isolation-by-resistance (IBR) on genetic differentiation in likely due to the species’ cosmopolitan habit and a Arabidopsis thaliana. a Coefficients (± SD) of the top-ranked nested maximum-likelihood population effects models (NMLPE) testing the Table 2 Summary statistics for the top-ranked NMLPE models effect of IBD, IBE and IBR on genetic differentiation in A. thaliana. evaluating the effect of isolation-by-distance (IBD), isolation-by- b Model averaged coefficients (± SD) for three Principal Component environment (IBE) and isolation-by-resistance (IBR) on the (PC) axes accounting for more than 75% of the total variance. Model genetic differentiation among Iberian Arabidopsis thaliana averaging was conducted using the subsample of models exhibiting populations ΔAIC < 2 regarding the top-ranked model, if more than one. In all cases, model estimates for the analysis conducted for the entire Cluster Predictors AIC ΔAIC Weight Iberian Peninsula (IP) and the four genetic clusters (NW-C1, NE-C2, A – Model: IBD + IBE + IBR relict-C3, and SW-C4) are shown. Maps with the geographic IP IBD + IBE −128,716.6 0.00 1.00 distribution of populations used in each analysis are also given. Maps were obtained from the National Center for Geographic NW-C1 IBD + IBE −42,836.0 0.00 1.00 Information (CNIG) of Spain NE-C2 IBD + IBE −10,406.3 0.00 1.00 Relict-C3 IBD + IBE − 2938.5 0.00 1.00 distribution without major discontinuities (Fig. S1). SW-C4 IBR − 3753.4 0.00 0.79 Thus, most population pairs are connected by relatively high habitat suitability, diminishing resistance distances. B – Model: IBD + PC1 + PC2 + PC3 + IBR The demographic history, characterized by the presence IP IBD + PC1 + PC2 + PC3 −128,694.0 0.00 1.00 of the relict African lineage and other non-relict Eur- NW-C1 IBD + PC1 + PC2 + PC3 −42,856.2 0.00 1.00 asian lineages [28–31, 34, 35, 50], and the well-known NE-C2 IBD + PC1 + PC2 + PC3 −10,380.3 0.00 0.83 ability to adjust flowering time and seed dormancy to Relict-C3 IBD + PC1 + PC2 + PC3 − 2936.6 0.00 0.89 contrasting environments [38, 51–55], both account for SW-C4 IBR −3753.4 0.00 0.41 the cosmopolitan habit of Iberian A. thaliana. We also explored the effect of IBD, IBE and IBR for each IBR + PC2 − 3751.6 1.84 0.17 of the four genetic clusters, which represent the joint out- Entries are given for the entire Iberian Peninsula (IP) and the four genetic clusters detected using all populations available (N = 278). The two models come of the long demographic and adaptive history of the include IBD, IBE and IBR, and IBE given the principal analysis components species in the region [30–32, 34, 39, 56]. Overall, NW-C1, (PCs). Akaike information criterion (AIC), ΔAIC and model weight are reported NE-C2 and relict-C3 exhibited the same patterns than Castilla et al. BMC Evolutionary Biology (2020) 20:71 Page 7 of 13 those at the regional scale, i.e. IBD and IBE mostly thaliana populations rather isolated and locally adapted to accounted for genetic differentiation within each cluster particular environments can have very low genetic diver- (Fig. 4a). However, it is worth noting that the importance sity [67], a scenario that could also apply to some of our of IBE was not the same for these three clusters, since populations. Finally, sampling bias might also affect gen- NE-C2 and relict-C3 exhibited higher IBE than NW-C1 etic diversity [68], particularly because some populations (Fig. 4a), stressing probably the higher influence of local were sampled late in the reproductive season. adaptation for genetic differentiation in these clusters. In Whatever the probable combination of factors account- fact, NE-C2 is a cluster with pronounced altitudinal gradi- ing for the existence of A. thaliana populations with no or ents in NE Spain (from maritime to sub-alpine environ- very low genetic diversity, the spatial hierarchical Bayesian ments) along which the species is known to have adapted model allowed the analysis of the geographic distribution of by adjusting various life-cycle and physiological traits [57– genetic diversity in Iberian A. thaliana (Fig. 2), particularly 61]. In the case of relict-C3, the long evolutionary history to detect the environmental predictors of the geographic of this lineage has provided the means to survive and distribution of genetic diversity. Precipitation, and not adapt to diverse Iberian environments over the last millen- temperature, emerged as the major environmental factor nia [35]. On the contrary, that was not the case for SW- accounting for the split of populations with and without C4, whose genetic differentiation was accounted for by genetic diversity, as well as for genetic diversity (Table 1). IBR (Fig. 4a). Whereas SW-C4 is the genetic cluster with This result is in agreement with other studies on the distri- the most restricted geographic distribution (Fig. 3c), its bution of genetic diversity in plants, which indicated a trend populations are located across an area showing dramatic for a higher role of precipitation variables over those of changes in habitat suitability (Fig. S1). Most likely, this temperature [69–74]. In addition, precipitation seasonality, particularity accounted for the enormous weight of IBR in which determines climatic variation in the Iberian Penin- this cluster. sula spanning from drier Mediterranean to more humid At- Finally, we also conducted NMLPE using principal com- lantic climates, was the most important bioclimatic variable ponent (PC) axes to pinpoint the major environmental for the geographic distribution of genetic diversity in A. predictors underlying IBE. We found that the three PCs thaliana (Table 1). Previous studies also found that precipi- depicting gradients of temperature and/or precipitation tation seasonality was a good predictor of individual fitness across the region made substantial contributions to IBE in [33] and distribution of the four genetic clusters [39]. Albeit NW-C1, NE-C2 and relict-C3 Iberian clusters (Fig. 4b), precipitation seasonality may be hard to handle as a fixed probably due to the similar contributions of the three PCs factor, further experiments are needed to better understand to the environmental variation among these populations. how A. thaliana responds to this variable, which might be In addition, as IBE informs on the influence of local adap- a major evolutionary force in this species. tation on genetic differentiation, our results reinforce the view of the strong effects of environmental gradients on Conclusions local adaptation in Iberian A. thaliana. For example, pre- This study dissected the complexity of geographic, envir- vious studies showed that geographic variation in fitness- onmental and evolutionary factors contributing to the related traits in Iberian A. thaliana (flowering time and genetic differentiation among A. thaliana populations, seed dormancy) strongly co-varied as a function of while illustrating the power of dense collections to dis- temperature and to a lesser extent of precipitation [30, 32, entangle complex biological questions [8, 28, 30, 49, 63, 33, 38, 55, 62, 63]. In particular, A. thaliana adjusts its life 75, 76]. Beyond genetic differentiation patterns, we also cycle by advancing flowering time and increasing seed identified the environmental predictors accounting for dormancy as the environment becomes warmer, drier and genetic diversity within populations. Clearly, the pro- more seasonal [38, 55, 63]. cesses driving among-population genetic differentiation As genetic diversity and its within- and among- may not be the same than those determining the amount population distribution patterns represents the raw mater- of within-population genetic diversity. However, we need ial upon which genetic differentiation is estimated, we dis- to deal with the latter to understand the former. In this cuss some results dealing with genetic diversity of Iberian sense, recent resurrection studies showed that natural A. A. thaliana that are worth considering. For example, ana- thaliana populations exhibited substantial changes in lyses of the genetic diversity in this large number of popu- their genetic composition in just a decade [77, 78]. The lations showed that 24% of them had no or very low extent of temporal genetic variation in A. thaliana rep- genetic diversity (Fig. 1a–c). This is not exceptional as resents a reminder that populations are not static, and other authors also detected populations with practically that the geographic distribution of genetic diversity no genetic diversity in A. thaliana [40, 41, 49, 64, 65]. within and among-populations is changing too. Repeated Strong founder effects [66] and low migration rates [41] genomic scans over time on the same populations will could well account for this result. Furthermore, small A. provide real-time insights into the intensity and pace of Castilla et al. BMC Evolutionary Biology (2020) 20:71 Page 8 of 13 genetic variation within populations, which will increase among-population partitioning of genetic diversity in A. our understanding of the evolutionary processes shaping thaliana, which yielded interpretable results with four in- genetic differentiation in plant populations. dividuals per population [40]. Besides, other studies using more individuals per population and different markers, Methods albeit at smaller geographical scales, also found popula- Source populations tions with no or very low genetic diversity (see below), We sampled 278 natural populations of the annual plant indicating that is not unusual in A. thaliana [40, 41, 49, Arabidopsis thaliana (L.) Heyhn. (Brassicaceae) across the 64, 65, 77]. entire Iberian Peninsula (~ 800 × 700 km2; 36.00–43.48°N, We collected all seeds from wild populations. Arabi- 3.19ºE–9.30°W; Fig. 1a) during the decade of the 2000s. C. dopsis thaliana is a common plant species not catego- Alonso-Blanco and F.X. Picó identified and collected all rized as protected or endangered in any species list of the material. No voucher specimens of this material were the Convention on the Trade in Endangered Species of deposited in a publicly available herbarium. For each Wild Fauna and Flora. We carried out field sampling in population, we recorded its geographic coordinates and locations where no permission was required, except at altitude using a GPS (Garmin International, Inc., Olathe, Doñana National Park (permission issued by Estación KS, USA). We extracted environmental data from publicly Biológica de Doñana) and Sierra de Grazalema Natural available repositories, such as WorldClim v.2 [79]; https:// Park (permission issued by Red Andaluza de Jardines www.worldclim.org/; accessed 6 June 2019), The CORINE Botánicos de la Consejería de Medio Ambiente y Orde- Land Cover 2000 (https://land.copernicus.eu/pan-euro- nación del Territorio de la Junta de Andalucía). Based pean/corine-land-cover; accessed 6 June 2019), and The on the Royal Decree of the Spanish legislation (Real Soil Geographical Database from Eurasia v.4 (https:// Decreto 124/2017, de 24 de febrero; https://www.boe.es/ esdac.jrc.ec.europa.eu/tags/soil-geographical-database-eur- eli/es/rd/2017/02/24/124), the genetic resources in- asia; accessed 6 June 2019). cluded in this study fall within the definition of “exclu- Geographic distance among populations varied between sively for taxonomic purposes” as they were used for 1 and 1059 km (Fig. 1a) and altitude between 1 and 2662 scientific, educational or non-commercial purposes. m.a.s.l. (Fig. 1b). For this set of populations, annual mean temperature varied between 5.3 and 18.4 °C (mean ± SD = Historical factors 12.3 ± 2.7 °C), annual total precipitation between 216.2 Recent historical changes in land use or large-scale per- and 1778.8 mm (mean ± SD = 760.9 ± 289.5 mm), percent- turbations may dramatically affect the amount of genetic age of natural vegetation (i.e. forests, scrublands or grass- diversity in plant populations. Such changes can be a mis- lands) between 0 and 100% (mean ± SD = 61.3 ± 35.7%), leading factor when attempting to evaluate the joint effects and soil pH between 3.6 and 7.5 (mean ± SD = 5.7 ± 0.8). of IBD, IBE and IBR on genetic differentiation. This is be- Natural A. thaliana populations are made of patches of cause IBD, IBE and IBR do not contemplate perturbations individuals widely differing in size and density, which is im- producing sudden or erratic changes in genetic diversity portant when it comes to design a sampling scheme to that eventually may mask the patterns of interest. To study the spatio-temporal distribution of genetic diversity make sure that historical factors were not an issue in this in A. thaliana [49, 57, 65, 77]. For each population and study, we estimated the temporal landscape dynamics of whenever possible, we collected seeds from several individ- all study populations. We examined all aerial orthophoto- uals from different patches (separated 1–20 m from each graphs available for each population (spanning between other) representing well each study population. Every sam- 1945 and 2016), which were retrieved from different pub- pling year and a few months after sampling, we multiplied lic regional administrations in Spain [77] with a tool spe- field-collected seeds by the single seed descent method, in a cifically developed for this purpose. We excluded up to 15 glasshouse at the Centro Nacional de Biotecnología (CNB- populations because aerial orthophotographs were not CSIC) in Madrid. Bulked seeds were stored in dry, dark available. The final number of aerial orthophotographs conditions in cellophane bags at room temperature, storing was highly variable per population, ranging between 4 and conditions that can preserve A. thaliana seeds for years. 24 (mean ± SE = 8.6 ± 0.3 aerial orthophotographs) and In this study, each population included about six indi- covering between 9 and 71 years (37.1 ± 1.7 years between viduals, ranging between four and seven. Given the im- the first and the last orthophotograph). portance to work with a dense collection of natural To quantify temporal landscape dynamics, a circular populations, we selected this number as a compromise be- area (500 m radius) around the GPS coordinate was di- tween the number of populations and the number of indi- vided into a regular grid of 80 squares (100 m side) for viduals per population that could be handled to extract each aerial orthophotograph available per population. the genetic patterns of interest. This choice was also based Based on vegetation and land use, each square was catego- on a previous large-scale study exploring the within- and rized as forests, scrublands, dehesa (i.e. agro-silvicultural Castilla et al. BMC Evolutionary Biology (2020) 20:71 Page 9 of 13 ecosystems based on Mediterranean oak woodlands), bare ancestry membership proportions of each accession in soil, crops, urban, infrastructures, and water, when one of each population, we ran the algorithm 20 times for each these categories occupied at least more than half of the defined number of groups (K value) from 1 to 10. The square. Categories were assigned to each square by digital- number of distinct genetic clusters was determined by izing manually all squares from all aerial orthophoto- evaluating the differences between the data likelihood for graphs available (N = 2252 orthophotographs) with QGIS successive K values. The largest K value with significantly v.3.4 [80]. For each habitat type and population, we esti- higher likelihood than that of K-1 runs (two-sided P < mated the average change between year intervals available 0.005; Wilcoxon tests for related samples) gave the final K for each population to quantify landscape changes over number. This was supported by a high similarity among time. the ancestry membership matrices from different runs of the same K value (H′ = 0.99). We used CLUMPP v.1 [86] SNP genotyping to calculate the average symmetric similarity coefficient H A total of 1772 A. thaliana individuals from 278 popula- ′ among runs and the average matrix of ancestry member- tions were genotyped with 245 presumably neutral nu- ship proportions, derived from the 10 runs having the clear SNPs using the SNPlex technique (Applied highest likelihood. STRUCTURE simulations were con- Biosystems, Foster City, CA, USA) through the CEGEN ducted at the The Supercomputing Center of Galicia Genotyping Service (www.usc.es/cegen/). These (CESGA; http://www.cesga.es/). genome-wide SNPs are frequent polymorphisms in Cen- For the DAPC analysis, we assessed the number of tral Europe, the Iberian Peninsula, and worldwide collec- clusters using the ‘find.cluster’ function in adegenet, tions, which altogether minimized ascertainment bias which runs successive K-means clustering with an in- [31, 32, 34]. On average, there were about 49 SNPs per creasing number of clusters to determine the best sup- chromosome (range = 45–53 SNPs) located at approxi- ported number of genetic clusters using the Bayesian mately 0.5 Mb from each other (range = 0.11 Kb – 1.82 Information Criterion (BIC). The K value with the low- Mb). Four SNPs had percentages of missing data above est BIC value represents the optimal number of clusters, 25% and one was monomorphic. We discarded these five although BIC values may keep decreasing after the true SNPs from the analyses maintaining 240 SNPs. K value in case of genetic clines and hierarchical struc- ture [87]. Therefore, the rate of decrease in BIC values Genetic analyses was visually examined to identify values of K, after which For each A. thaliana population, we calculated the per- BIC values only decreased in a subtle manner [87]. The centage of polymorphic loci (PL), the mean number of ‘dapc’ function was used with the final K value, retaining observed alleles per locus (na) and mean gene diversity the axes of the Principal Component Analysis. (HS) using FSTAT v.2.9.3 [81]. In addition, we computed the percentage of differences among all pairs of non- Geographic distribution of genetic diversity redundant multilocus genotypes. We used the ‘pairwi- We modeled the effects of environmental heterogeneity se.WCfst’ function implemented in the R package on the spatial distribution of genetic diversity by using a hierfstat v.0.04–22 [82] to estimate genetic differen- modified version of a spatial hierarchical Bayesian model tiation as pairwise FST values according to Weir and recently used to model the effects of warming on the Cockerham [68], which represents the dependent vari- Iberian A. thaliana’s distribution [56]. The aim of this able to test IBD, IBE and IBR (see below). model was to visualize hot and cold spots of genetic di- The genetic structure of A. thaliana populations across versity across the Iberian Peninsula and to identify their the Iberian Peninsula was examined using two different environmental predictors. methods: the Bayesian clustering method implemented in We used a Bayesian approach to handle the particular- STRUCTURE v.2.3.3 [83, 84] and the ordination method ities of genetic diversity data, such as the semi-continuous represented by the Discriminant Analysis of Principal nature of the variable. In the case of Iberian A. thaliana, Components (DAPC) [85] available in the R package 66 of 278 populations had very low mean gene diversity adegenet v.2.1.1 [85]. In the case of STRUCTURE, model (HS) values (HS < 0.009), whereas the rest of populations settings included haploid non-redundant multilocus geno- exhibited non-zero HS values (Fig. 1c). To handle such types, correlated allele frequencies between populations complex distribution of HS, we developed a Hurdle-beta and a linkage model. We identified the number of identi- Bayesian model [88] to separate the zero structure (ab- cal multilocus genotypes using the ‘mlg.filter’ function im- sence of genetic diversity) from the non-zero structure plemented in adegenet. Each run consisted of 50,000 (presence of genetic diversity) of data. The Hurdle model burn-in MCMC iterations and 100,000 MCMC after- is defined as a finite mixture of two processes: a degener- burning repetitions for parameter estimation. To deter- ate distribution with point mass and zero, which deter- mine the K number of ancestral populations and the mines the presence or absence of data—in our case HS Castilla et al. BMC Evolutionary Biology (2020) 20:71 Page 10 of 13 values above or below the value of 0.009, respectively— 478 A. thaliana populations available to date (collected and a conditional-to-presence continuous process sup- between 2000 and 2019; Fig. S1), and the same environ- ported by the open interval between the lowest non-zero mental variables used in a previous habitat suitability HS value and 1. The first process (absence of genetic diver- model [39]. In particular, the model included topsoil pH sity) was modeled with a Bernoulli distribution, whilst the and eight bioclimatic variables, which were annual mean continuous component (presence of genetic diversity) was temperature (BIO1), mean diurnal temperature range modeled with a Beta distribution (see the full development (BIO2), isothermality (BIO3), temperature seasonality of the model in the Supporting Information section). The (BIO4), mean temperature of wettest quarter (BIO8), an- predictors of the model were the 19 bioclimatic variables nual precipitation (BIO12), precipitation seasonality available from WorldClim and topsoil pH available from (BIO15) and precipitation of the warmest quarter The Geographical Database from Eurasia. The Watanabe– (BIO18). We applied Maxent using default parameters Akaike information criterion (WAIC) [89]; was computed except for features using the hinge type, making it com- to determine the best models. The model included a sto- parable to a Generalized Additive Model [94]. We ob- chastic spatial effect (the spatial term) to remove the tained consistent results with previous habitat suitability spatial autocorrelation of data [56]. models conducted for Iberian A. thaliana with lower sam- ple sizes [39, 56]. Drivers of genetic differentiation Next, we transformed the resulting habitat suitability We estimated the ecological, genetic and evolutionary into resistance values as 1 minus the value of each pixel. drivers of genetic differentiation of Iberian A. thaliana Thus, greater pixel values represented lower occurrence as follows. IBD was based on geographic distance. How- probability of A. thaliana (higher resistance) and vice- ever, instead of using the Euclidian distance among versa. Using circuit theory [24, 25], we examined population pairs, we created a raster by assigning a 0.5 whether genetic differentiation among A. thaliana popu- value to all 1-km2 pixels [90, 91] of the Iberian Peninsula lations was accounted for by IBR. We used Circuits- map. We then calculated pairwise distances, known as cape v.4.0 [95] to calculate resistance distance matrices geographic resistance distance, among all A. thaliana assigning pixel values as resistance values. We calculated populations employing the new raster. We did not use resistance distance among all pairs of populations across the Euclidean distance to compute IBD because straight Iberian Peninsula, as well as among all pairs of popula- lines between population pairs may introduce some bias tions within each genetic cluster. when such straight lines traverses natural barriers for To quantify the relative contributions of IBD, IBE and the study organism (i.e. the coastline or open sea) [92]. IBR on the genetic differentiation (pairwise FST values) IBE was based on the 19 bioclimatic variables. In among A. thaliana populations, we fitted maximum- addition, we used The CORINE Land Cover to obtain likelihood population effect models (MLPE) [96]. Code the percentage of human-modified habitat (i.e. urban implementing the MLPE correlation structure within the areas, crops and semi-natural grasslands) within a 500 m R package nlme v. 3.1–143 comes from the corMLPE radius from the population GPS coordinate [39]. The package. The model used penalized least squares and a percentages of human-modified and natural habitat were residual covariance structure designed to account for the significantly negatively correlated (N = 278; r = − 0.95; non-independence of pairwise distances. Given the in- P < 0.001; Pearson’s correlation) and nearly summed to herent spatial dependence structure in pairwise compari- 100%. Finally, a topsoil pH layer was also used. Combin- sons between A. thaliana populations, all MLPE ing these 21 environmental variables, we conducted a considered spatial locations. To this end, we used a principal component analysis (PCA) with varimax rota- modification of the MLPE model incorporating the cor- tion using the software SPSS v.25 (IBM, Chicago, IL relation between pairwise measurements due to com- USA). Thus, we obtained the PC scores of the first three parison of populations and spatial locations (nested principal components (PC) for each population. We then MLPE or NMLPE) [91]. Since the independent variables used the PC scores to calculate environmental dissimi- (IBD, IBE and IBR) used in NMLPE had very different larity between populations as Euclidean distances de- units, we normalized them by subtracting the mean and rived with the ‘dist’ function in R. Furthermore, we also dividing by the standard deviation in all analyses. We calculated dissimilarity matrices among populations for did not linearize pairwise FST values. We conducted each PC axis separately, as they were interpretable in NMLPE on the complete dataset (N = 278 populations) terms of environmental gradients. as well as on the subset of those populations exhibiting To compute IBR, we first estimated habitat suitability HS different from zero (N = 212 populations with HS > for A. thaliana across the Iberian Peninsula with the 0.009). The results were consistent between the two maximum-entropy modeling technique implemented in datasets, showing that the 66 populations with no gen- the software Maxent 3.3.3 k [93, 94]. To do that, we used etic diversity were not biasing the results. Castilla et al. BMC Evolutionary Biology (2020) 20:71 Page 11 of 13 In addition, we ran NMLPE for each genetic cluster to Funding evaluate whether the effects of IBD, IBE and IBR on gen- This research has been funded by grant BIO2016–75754-P from the Agencia Estatal de Investigación (AEI) of Spain and the Fondo Europeo de Desarrollo etic differentiation differed among genetic clusters. This Regional (FEDER, UE) to CA-B. The funding body played no role in the design is because the number of genetic clusters represents the of the study and collection, analysis, and interpretation of data and in writing outcome of the demographic and evolutionary history of the manuscript. A. thaliana in the region. We considered each popula- Availability of data and materials tion to belong to a specific genetic cluster when the Data deposited in the Dryad repository: https://doi.org/10.5061/dryad. average value of the membership proportions of its indi- 47d7wm393 [98]. The tool to retrieve and examine aerial orthophotographs viduals was ≥0.3 for that specific cluster. Given the deposited in the Zenodo repository: https://zenodo.org/record/2552031. marked genetic structure of the study system, threshold Ethics approval and consent to participate values of 0.25 and 0.40 were consistent, but 0.3 provided The species sampled is not endangered or protected. a clearer classification of individuals into genetic clus- ters. The majority of populations were assigned to a sin- Consent for publication gle cluster (N = 230). We assigned populations reaching Not applicable. the minimum average membership proportion for more Competing interests than one cluster to multiple clusters. Overall, 47 popula- The authors declare that they have no competing interests. tions were assigned to two clusters, whereas only one population had average membership proportions higher Author details 1 Centre for Applied Ecology “Prof. Baeta Neves”, InBIO, School of Agriculture, than 0.3 for three clusters. University of Lisbon, Lisbon, Portugal. 2Departamento de Ecología Integrativa, We also ran NMLPE on the complete dataset and each Estación Biológica de Doñana (EBD), Consejo Superior de Investigaciones cluster separately, but replacing IBE by the PC axes. Thus, Científicas (CSIC), Sevilla, Spain. 3Departamento de Genética Molecular de Plantas, Centro Nacional de Biotecnología (CNB), Consejo Superior de we used IBD, PC1, PC2, PC3 and IBR as independent vari- Investigaciones Científicas (CSIC), Madrid, Spain. 4CREAF, Centre de Recerca ables. We used the Akaike information criterion (AIC) to Ecològica i Aplicacions Forestals, Bellaterra, E08193 Cerdanyola de Vallès, compare models containing all possible combinations of Catalonia, Spain. 5Universitat Autònoma de Barcelona, Bellaterra, E08193 Cerdanyola de Vallès, Catalonia, Spain. 6BCAM - Basque Center for Applied non-collinear predictions (r < 0.6), created with the Mathematics, Bilbao, Spain. 7Departament d’Estadística i Investigació ‘dredge’ function from the R package MuMIn v1.4 [97]. In Operativa, Universitat de València, Valencia, Spain. those cases where several models were top-ranked with Received: 13 February 2020 Accepted: 1 June 2020 ΔAIC < 2, we calculated the model averaged parameter es- timates (β) with standard errors for the explanatory vari- ables included in all top-ranked models also using MuMIn. References We considered effect sizes as significant when 95% confi- 1. Barrett RDH, Schluter D. Adaptation from standing genetic variation. Trends Ecol Evol. 2008;23:38–44. dence intervals did not overlap zero. 2. Jump AS, Marchant R, Peñuelas J. Environmental change and the option value of genetic diversity. Trends Plant Sci. 2009;14:51–8. Abbreviations 3. Hancock AM, Brachi B, Faure N, Horton MW, Jarymowycz LB, Sperone FG, AIC: Akaike information criterion; BIC: Bayesian information criterion; et al. Adaptation to climate across the Arabidopsis thaliana genome. DAPC: Discriminant analysis of principal components; IBD: Isolation-by- Science. 2011;334:83–6. distance; IBE: Isolation-by-environment; IBR: Isolation-by-resistance; 4. Wilczek AM, Cooper MD, Korves TM, Schmitt J. Lagging adaptation to warming MCMC: Monte Carlo Markov chain; NMLPE: Nested maximum-likelihood climate in Arabidopsis thaliana. P Natl Acad Sci USA. 2014;111:7906–13. population effect models; SNP: Single nucleotide polymorphism; 5. Matuszewski S, Hermisson J, Kopp M. Catch me if you can: adaptation from WAIC: Watanabe–Akaike information criterion standing genetic variation to a moving phenotypic optimum. Genetics. 2015;200:1255–74. 6. Fournier-Level A, Perry EO, Wang JA, Braun PT, Migneault A, Cooper MD, Acknowledgements et al. Predicting the evolutionary dynamics of seasonal adaptation to novel We thank Mercedes Ramiro for technical assistance. We are grateful to David climates in Arabidopsis thaliana. P Natl Acad Sci USA. 2016;113:E2812–21. Ragel, Marisol Muñoz and Rocío Gómez for digitalizing all orthophotographs 7. Exposito-Alonso M, Brennan AC, Alonso-Blanco C, Picó FX. Spatio-temporal used in this study. We are indebted to Joaquín Ortego, Rodolfo Jaffé, variation in fitness responses to contrasting environments in Arabidopsis Nathaniel Pope, and Bernardo Toledo for their insightful and constructive thaliana. Evolution. 2018;72:1570–86. comments on various aspects of the manuscript during its development. We 8. Exposito-Alonso M. 500 genomes field experiment team, Burbano HA, thank two anonymous reviewers and the editor for constructive comments Bossdorf O, Nielsen R, Weigel D. natural selection on the Arabidopsis on this manuscript. We particularly thank Nathaniel Pope for his help with thaliana genome in present and future climates. Nature. 2019;573:126–9. the implementation of NMLPE. The Supercomputing Center of Galicia 9. Fisher RA. The genetical theory of natural selection. Oxford: Oxford (CESGA) facilitated the analyses of genetic data. We acknowledge support of University Press; 1930. the publication fee by the CSIC Open Access Publication Support Initiative 10. Orr HA. The genetic theory of adaptation: a brief history. Nat Rev Genet. through its Unit of Information Resources for Research (URICI). 2005;6:119–27. 11. Hughes AL. The origin of adaptive phenotypes. P Natl Acad Sci USA. 2008; Authors’ contributions 105:13193–4. CA-B and FXP conceived the idea. BM-V, CA-B and FXP collected seeds, grew 12. Rauch EM, Bar-Yam Y. Theory predicts the uneven distribution of genetic plants and generated the SNP data. ARC and AM analyzed the data. AM de- diversity within species. Nature. 2004;431:449–52. veloped the habitat suitability model. JM-M and DC developed the Bayesian 13. Wright S. Isolation by distance. Genetics. 1943;28:114–38. model. FXP wrote the first version of the manuscript and all authors contrib- 14. Kimura M, Weiss GH. The stepping stone model of population structure and uted to it. the decrease of genetic correlation with distance. Genetics. 1964;49:561–76. Castilla et al. BMC Evolutionary Biology (2020) 20:71 Page 12 of 13 15. Malécot G. Heterozygosity and relationship in regularly subdivided 40. Bakker EG, Stahl EA, Toomajan C, Nordborg M, Kreitman M, Bergelson J. populations. Theor Pop Biol. 1975;8:212–41. Distribution of genetic variation within and among local populations of 16. Slatkin M. Isolation by distance in equilibrium and non-equilibrium Arabidopsis thaliana over its species range. Mol Ecol. 2006;15:1405–18. populations. Evolution. 1993;47:264–79. 41. Falahati-Anbaran M, Lundemo S, Stenøien HK. Seed dispersal in time can 17. Rousset F. Genetic differentiation and estimation of gene flow from F- counteract the effect of gene flow between natural populations of statistics under isolation by distance. Genetics. 1997;145:1219–28. Arabidopsis thaliana. New Phytol. 2014;202:1043–54. 18. Sexton JP, Hangartner SB, Hoffmann AA. Genetic isolation by environment 42. Sharbel TF, Haubold B, Mitchell-Olds T. Genetic isolation by distance in or distance: which pattern of gene flow is most common? Evolution. Arabidopsis thaliana: biogeography and postglacial colonization of Europe. 2014;68:1–15. Mol Ecol. 2000;9:2109–18. 19. Jenkins DG, Carey M, Czerniewska J, Fletcher J, Hether T, Jones A, et al. A 43. Nordborg M, Hu TT, Ishino Y, Jhaveri J, Toomajian C, Zheng H, et al. The meta-analysis of isolation by distance: relic or reference standard for pattern of polymorphism in Arabidopsis thaliana. PLoS Biol. 2005;3:e196. landscape genetics? Ecography. 2010;33:315–20. 44. Ostrowski M-F, David J, Santoni S, McKhann H, Reboud X, Le Corre V, et al. 20. Aguillon SM, Fitzpatrick JW, Bowman R, Schoech SJ, Clark AG, Coop G, et al. Evidence for a large-scale population structure among accessions of Deconstructing isolation-by-distance: the genomic consequences of limited Arabidopsis thaliana: possible causes and consequences for the distribution dispersal. PLoS Genet. 2017;13:e1006911. of linkage disequilibrium. Mol Ecol. 2006;15:1507–17. 21. Wang IJ, Summers K. Genetic structure is correlated with phenotypic 45. Beck JB, Schmuths K, Schaal BA. Native range genetic variation in divergence rather than geographic isolation in the highly polymorphic Arabidopsis thaliana is strongly geographically structured and reflects strawberry poison-dart frog. Mol Ecol. 2010;19:447–58. Pleistocene glacial dynamics. Mol Ecol. 2008;17:902–15. 22. Wang IJ, Bradburd GS. Isolation by environment. Mol Ecol. 2014;23:5649–62. 46. Lundemo S, Falahati-Anbaran M, Stenøien HK. Seed banks cause elevated 23. Manel S, Schwartz MK, Luikart G, Taberlet P. Landscape genetics: generation times and effective population sizes of Arabidopsis thaliana in combining landscape ecology and population genetics. Trends Ecol northern Europe. Mol Ecol. 2009;18:2798–811. Evol. 2003;18:189–97. 47. Lewandowska-Sabat AM, Fjellheim S, Rognli OA. Extremely low genetic 24. McRae BH. Isolation by resistance. Evolution. 2006;60:1551–61. variability and highly structured local populations of Arabidopsis thaliana at 25. McRae BH, Beier P. Circuit theory predicts gene flow in plant and animal higher latitudes. Mol Ecol. 2010;19:4753–64. populations. P Natl Acad Sci USA. 2007;104:19885–90. 48. Platt A, Horton M, Huang YS, Li Y, Anastasio AE, Mulyati NW, et al. The scale of 26. Spear SF, Balkenhol N, Fortin M-J, McRae BH, Scribner K. Use of resistance population structure in Arabidopsis thaliana. PLoS Genet. 2010;6:e1000843. surfaces for landscape genetic studies: considerations for parameterization 49. Bomblies K, Yant L, Laitinen RA, Kim ST, Hollister JD, Warthmann N, et al. and analysis. Mol Ecol. 2010;19:3576–91. Local-scale patterns of genetic variability, outcrossing, and spatial structure 27. Lundgren E, Ralph PL. Are populations like a circuit? Comparing isolation by in natural stands of Arabidopsis thaliana. PLoS Genet. 2010;6:e1000890. resistance to a new coalescent-based method. Mol Ecol Resour. 2019;19: 50. Lee C-R, Svardal H, Farlow A, Exposito-Alonso M, Ding W, Novikova P, et al. 1388–406. On the post-glacial spread of human commensal Arabidopsis thaliana. Nat 28. The 1001 Genomes Consortium. 1,135 Genomes reveal the global pattern Commun. 2017;8:14458. of polymorphism in Arabidopsis thaliana. Cell. 2016;166:481–91. 51. Donohue K. Niche construction through phenological plasticity: life history 29. Durvasula A, Fulgione A, Gutaker RM, Alacakaptan SI, Flood PJ, Neto C, et al. dynamics and ecological consequences. New Phytol. 2005;166:83–92. African genomes illuminate the early history and transition to selfing in 52. Debieu M, Tang C, Stich B, Sikosek T, Effgen S, Josephs E, et al. Co-variation Arabidopsis thaliana. P Natl Acad Sci USA. 2017;114:5213–8. between seed dormancy, growth rate and flowering time changes with 30. Tabas-Madrid D, Méndez-Vigo B, Arteaga N, Marcer A, Pascual-Montano A, latitude in Arabidopsis thaliana. PLoS One. 2013;8:e61075. Weigel D, et al. Genome-wide signatures of flowering adaptation to climate 53. Springthorpe V, Penfield S. Flowering time and seed dormancy control use temperature: regional analyses in a highly diverse native range of external coincidence to generate life history strategy. eLife. 2015;4:e05557. Arabidopsis thaliana. Plant Cell Environ. 2018;41:1806–20. 54. Burghardt LT, Edwards BR, Donohue K. Multiple paths to similar germination 31. Picó FX, Méndez-Vigo B, Martínez-Zapater JM, Alonso-Blanco C. Natural behavior in Arabidopsis thaliana. New Phytol. 2016;209:1301–12. genetic variation of Arabidopsis thaliana is geographically structured in the 55. Marcer A, Vidigal DS, James PMA, Fortin M-J, Méndez-Vigo B, Hilhorst HWM, Iberian Peninsula. Genetics. 2008;180:1009–21. et al. Temperature fine-tunes Mediterranean Arabidopsis thaliana life-cycle 32. Méndez-Vigo B, Picó FX, Ramiro M, Martínez-Zapater JM, Alonso-Blanco C. phenology geographically. Plant Biol. 2018;20:148–56. Altitudinal and climatic adaptation is mediated by flowering traits and FRI, 56. Martínez-Minaya J, Conesa D, Fortin M-J, Alonso-Blanco C, Picó FX, Marcer FLC and PHYC genes in Arabidopsis. Plant Physiol. 2011;157:1942–55. A. A hierarchical Bayesian Beta regression approach to study the effects of 33. Manzano-Piedras E, Marcer A, Alonso-Blanco C, Picó FX. Deciphering the geographical genetic structure and spatial autocorrelation on species adjustment between environment and life history in annuals: lessons from a distribution range shifts. Mol Ecol Resour. 2019;19:929–43. geographically-explicit approach in Arabidopsis thaliana. PLoS One. 2014;9: 57. Montesinos A, Tonsor SJ, Alonso-Blanco C, Picó FX. Demographic and e87836. genetic patterns of variation among populations of Arabidopsis thaliana 34. Brennan AC, Méndez-Vigo B, Haddioui A, Martínez-Zapater JM, Picó FX, from contrasting native environments. PLoS One. 2009;4:e7213. Alonso-Blanco C. The genetic structure of Arabidopsis thaliana in the South- 58. Montesinos-Navarro A, Wig J, Picó FX, Tonsor SJ. Arabidopsis thaliana Western Mediterranean range reveals a shared history between North Africa populations show clinal variation in a climatic gradient associated with and southern Europe. BMC Plant Biol. 2014;14:17. altitude. New Phytol. 2011;189:282–94. 35. Toledo B, Marcer A, Méndez-Vigo B, Alonso-Blanco C, Picó FX. An ecological 59. Montesinos-Navarro A, Picó FX, Tonsor SJ. Clinal variation in seed traits history of the relict genetic lineage of Arabidopsis thaliana. Environ Exp Bot. influencing life cycle timing in Arabidopsis thaliana. Evolution. 2012;66:3417–31. 2020;170:103800. 60. Wolfe MD, Tonsor SJ. Adaptation to spring heat and drought in 36. Legendre P, Dale MRT, Fortin M-J, Gurevitch J, Hohn M, Myers D. The northeastern Spanish Arabidopsis thaliana. New Phytol. 2014;201:323–34. consequences of spatial structure for the design and analysis of ecological 61. Zhang N, Belsterling B, Raszewski J, Tonsor SJ. Natural populations of field surveys. Ecography. 2002;25:601–15. Arabidopsis thaliana differ in seedling responses to high-temperature stress. 37. Méndez-Vigo B, Savic M, Ausín I, Ramiro M, Martín B, Picó FX, Alonso-Blanco AoB Plants. 2015;7:plv101. C. Environmental and genetic interactions reveal FLOWERING LOCUS C as a 62. Kronholm I, Picó FX, Alonso-Blanco C, Goudet J, de Meaux J. Genetic basis modulator of the natural variation for the plasticity of flowering in of adaptation in Arabidopsis thaliana: local adaptation at the seed dormancy Arabidopsis. Plant Cell Environ. 2016;39:282–94. QTL DOG1. Evolution. 2012;66:2287–302. 38. Vidigal DS, Marques ACSS, Willems LAJ, Buijs G, Méndez-Vigo B, Hilhorst 63. Martínez-Berdeja A, Stitzer MC, Taylor MA, Okada M, Ezcurra E, Runcie DE, HWM, et al. Altitudinal and climatic associations of seed dormancy and et al. Functional variants of DOG1 control seed chilling responses and flowering traits evidence adaptation of annual life cycle timing in variation in seasonal life-history strategies in Arabidopsis thaliana. P Natl Arabidopsis thaliana. Plant Cell Environ. 2016;39:1737–48. Acad Sci USA. 2020;117:2526–34. 39. Marcer A, Méndez-Vigo B, Alonso-Blanco C, Picó FX. Tackling intraspecific 64. Stenøien HK, Fenster CB, Tonteri A, Savolainen O. Genetic variability in genetic structure in distribution models better reflects species geographical natural populations of Arabidopsis thaliana in northern Europe. Mol Ecol. range. Ecol Evol. 2016;6:2084–97. 2005;14:137–48. Castilla et al. BMC Evolutionary Biology (2020) 20:71 Page 13 of 13 65. Gomaa NH, Montesinos-Navarro A, Alonso-Blanco C, Picó FX. Temporal 90. Castilla AR, Pope N, Jaffé R, Jha S. Elevation, not deforestation, promotes variation in genetic diversity and effective population size of Mediterranean genetic differentiation in a pioneer tropical tree. PLoS One. 2016;11: and subalpine Arabidopsis thaliana populations. Mol Ecol. 2011;20:3540–54. e0156694. 66. Mayr E. Systematics and the origin of species: from the viewpoint of a 91. Jaffé R, Veiga JC, Pope NS, Lanes ECM, Carvalho CS, Alves R, et al. Landscape zoologist. New York: Columbia University Press; 1942. genomics to the rescue of a tropical bee threatened by habitat loss and 67. Méndez-Vigo B, Gomaa NH, Alonso-Blanco C, Picó FX. Among- and within- climate change. Evol Appl. 2019;12:1164–77. population variation in flowering time of Iberian Arabidopsis thaliana estimated 92. Pereoglou F, Lindenmayer DB, MacGregor C, Ford F, Wood J, Banks SC. in field and glasshouse conditions. New Phytol. 2013;197:1332–43. Landscape genetics of an early successional specialist in a disturbance- 68. Weir BS, Cockerham CC. Estimating F-statistics for the analysis of population prone environment. Mol Ecol. 2013;22:1267–81. structure. Evolution. 1984;38:1358–70. 93. Phillips SJ, Anderson RP, Schapire RE. Maximum entropy modeling of 69. Manel S, Gugerli F, Thuiller W, Alvarez N, Legendre P, Holderegger R, et al. species geographic distributions. Ecol Model. 2006;190:231–59. Broad-scale adaptive genetic variation in alpine plants is driven by 94. Elith J, Phillips SJ, Hastie T, Dudík M, Chee YE, Yates CJ. A statistical temperature and precipitation. Mol Ecol. 2012;21:3729–38. explanation of MaxEnt for ecologists. Divers Distrib. 2011;17:43–57. 70. Cushman SA, Max T, Meneses N, Evans LM, Ferrier S, Honchak B, et al. 95. McRae BH, Dickson BG, Keitt TH, Shah VB. Using circuit theory to model Landscape genetic connectivity in a riparian foundation tree is jointly driven connectivity in ecology, evolution, and conservation. Ecology. 2008;89:2712–24. by climatic gradients and river networks. Ecol Appl. 2004;24:1000–14. 96. Clarke RT, Rothery P, Raybould AF. Confidence limits for regression 71. Mosca E, González-Martínez SC, Neale DB. Environmental versus relationships between distance matrices: estimating gene flow with geographical determinants of genetic structure in two subalpine conifers. distance. J Agr Biol Envir St. 2002;7:361–72. New Phytol. 2014;201:180–92. 97. Bartoń K. Package ‘MuMIn’. 2019. https://CRAN.R-project.org/package= 72. Mosca E, Gugerli F, Eckert AJ, Neale DB. Signatures of natural selection on MuMIn. Accessed 15 Oct 2019. Pinus cembra and P mugo along elevational gradients in the Alps. Tree 98. Picó FX, Castilla AR, Méndez-Vigo B, Marcer A, Martínez-Minaya J, Conesa D, Genet Genomes. 2016;12:9. Alonso-Blanco C. Data from: ecological, genetic and evolutionary drivers of 73. Bothwell HM, Cushman SA, Woolbright SA, Hersch-Green EI, Evans LM, regional genetic differentiation in Arabidopsis thaliana. Dryad Digital Whitham TG, et al. Conserving threatened riparian ecosystems in the Repository. 2020. https://doi.org/10.5061/dryad.47d7wm393. American west: precipitation gradients and river networks drive genetic connectivity and diversity in a foundation riparian tree (Populus angustifolia). Publisher’s Note Mol Ecol. 2017;26:5114–32. Springer Nature remains neutral with regard to jurisdictional claims in 74. Cuervo-Alarcon L, Arend M, Müller M, Sperisen C, Finkeldey R, Krutovsky KV. published maps and institutional affiliations. Genetic variation and signatures of natural selection in populations of European beech (Fagus sylvatica L.) along precipitation gradients. Tree Genet Genomes. 2018;14:84. 75. Lasky JR, Des Marais DL, McKay JK, Richards JH, Juenger TE, Keitt TH. Characterizing genomic variation of Arabidopsis thaliana: the roles of geography and climate. Mol Ecol. 2012;21:5512–29. 76. Frachon L, Bartoli C, Carrère S, Bouchez O, Chaubet A, Gautier M, et al. A genomic map of climate adaptation in Arabidopsis thaliana at a micro- geographic scale. Front Plant Sci. 2018;9:967. 77. Gómez R, Méndez-Vigo B, Marcer A, Alonso-Blanco C, Picó FX. Quantifying temporal change in plant population attributes: insights from a resurrection approach. AoB Plants. 2018;10:ply063. 78. Frachon L, Libourel C, Villoutreix R, Carrère S, Glorieux C, Huard-Chauveau C, et al. Intermediate degrees of synergistic pleiotropy drive adaptive evolution in ecological time. Nature Ecol Evol. 2017;1:1551–61. 79. Fick SE, Hijmans RJ. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. Int J Climatol. 2017;37:4302–15. 80. QGIS Development Team. QGIS Geographic information system. Open source geospatial foundation project. 2017. http://qgis.osgeo.org. Accessed 15 Mar 19. 81. Goudet J. FSTAT: a computer program to calculate F-statistics. J Hered. 1995;86:485–6. 82. Goudet J. HIERFSTAT, a package for R to compute and test hierarchical F- statistics. Mol Ecol Notes. 2005;5:184–6. 83. Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155:945–59. 84. Falush D, Stephens M, Pritchard JK. Inference of population structure using multilocus genotype data: linked loci and correlated allele frequencies. Genetics. 2003;164:1567–87. 85. Jombart T. Adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics. 2008;24:1403–5. 86. Jakobsson M, Rosenberg NA. CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics. 2007;23:1801–6. 87. Jombart T, Devillard S, Balloux F. Discriminant analysis of principal components: a new method for the analysis of genetically structured populations. BMC Genet. 2010;11:94. 88. Quiroz ZC, Prates MO, Rue H. A Bayesian approach to estimate the biomass of anchovies off the coast of Perú. Biometrics. 2015;71:208–17. 89. Watanabe S. Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. J Mach Learn Res. 2010;11:3571–94.