In the presence of population structure: From genomics to candidate genes underlying local adaptation

Abstract Understanding the genomic signatures, genes, and traits underlying local adaptation of organisms to heterogeneous environments is of central importance to the field evolutionary biology. To identify loci underlying local adaptation, models that combine allelic and environmental variation while controlling for the effects of population structure have emerged as the method of choice. Despite being evaluated in simulation studies, there has not been a thorough investigation of empirical evidence supporting local adaptation across these alleles. To evaluate these methods, we use 875 Arabidopsis thaliana Eurasian accessions and two mixed models (GEMMA and LFMM) to identify candidate SNPs underlying local adaptation to climate. Subsequently, to assess evidence of local adaptation and function among significant SNPs, we examine allele frequency differentiation and recent selection across Eurasian populations, in addition to their distribution along quantitative trait loci (QTL) explaining fitness variation between Italy and Sweden populations and cis‐regulatory/nonsynonymous sites showing significant selective constraint. Our results indicate that significant LFMM/GEMMA SNPs show low allele frequency differentiation and linkage disequilibrium across locally adapted Italy and Sweden populations, in addition to a poor association with fitness QTL peaks (highest logarithm of odds score). Furthermore, when examining derived allele frequencies across the Eurasian range, we find that these SNPs are enriched in low‐frequency variants that show very large climatic differentiation but low levels of linkage disequilibrium. These results suggest that their enrichment along putative functional sites most likely represents deleterious variation that is independent of local adaptation. Among all the genomic signatures examined, only SNPs showing high absolute allele frequency differentiation (AFD) and linkage disequilibrium (LD) between Italy and Sweden populations showed a strong association with fitness QTL peaks and were enriched along selectively constrained cis‐regulatory/nonsynonymous sites. Using these SNPs, we find strong evidence linking flowering time, freezing tolerance, and the abscisic‐acid pathway to local adaptation.

structure have emerged as the method of choice. Despite being evaluated in simulation studies, there has not been a thorough investigation of empirical evidence supporting local adaptation across these alleles. To evaluate these methods, we use 875 Arabidopsis thaliana Eurasian accessions and two mixed models (GEMMA and LFMM) to identify candidate SNPs underlying local adaptation to climate. Subsequently, to assess evidence of local adaptation and function among significant SNPs, we examine allele frequency differentiation and recent selection across Eurasian populations, in addition to their distribution along quantitative trait loci (QTL) explaining fitness variation between Italy and Sweden populations and cis-regulatory/nonsynonymous sites showing significant selective constraint. Our results indicate that significant LFMM/GEMMA SNPs show low allele frequency differentiation and linkage disequilibrium across locally adapted Italy and Sweden populations, in addition to a poor association with fitness QTL peaks (highest logarithm of odds score). Furthermore, when examining derived allele frequencies across the Eurasian range, we find that these SNPs are enriched in low-frequency variants that show very large climatic differentiation but low levels of linkage disequilibrium. These results suggest that their enrichment along putative functional sites most likely represents deleterious variation that is independent of local adaptation. Among all the genomic signatures examined, only SNPs showing high absolute allele frequency differentiation (AFD) and linkage disequilibrium (LD) between Italy and Sweden populations showed a strong association with fitness QTL peaks and were enriched along selectively constrained cis-regulatory/nonsynonymous sites. Using these SNPs, we find strong evidence linking flowering time, freezing tolerance, and the abscisic-acid pathway to local adaptation.
An important hurdle that GEA and other methods need to overcome, is disentangling the effects of selection from demographic history (Hoban et al., 2016;Lotterhos & Whitlock, 2014).
To examine the performance of population genomic methods (primarily F ST -based and GEA methods), there have been a multitude of simulation studies examining different demographic scenarios and selection regimes (De Mita et al., 2013;De Villemereuil, Frichot, Bazin, François, & Gaggiotti, 2014;Forester, Lasky, Wagner, & Urban, 2018;Lotterhos & Whitlock, 2014Perez-Figueroa, Garcia-Pereira, Saura, Rolan-Alvarez, & Caballero, 2010;Yoder & Tiffin, 2017). Some of the general findings of these studies are (a) that the underlying demographic history can have a significant impact on which approach (F ST or GEA) performs better (Lotterhos & Whitlock, 2015); (b) that under certain demographic scenarios the power of all methods can be low (De Mita et al., 2013;De Villemereuil et al., 2014); (c) the intensity/timing of selection, in addition to the number of loci involved, can have a significant impact on performance (De Villemereuil et al., 2014;Forester et al., 2018;Yoder & Tiffin, 2017); and (d) the sampling of individuals and markers can have a significant impact on the power to detect causative loci (De Mita et al., 2013;Yoder & Tiffin, 2017).

K E Y W O R D S
flowering time, population genomics, population structure, quantitative trait loci, selective constraint While simulation studies provide a controlled environment to evaluate methodologies, and determine the factors that can affect each method, we do not know how well they emulate real data and the complex demographic and selection forces underlying them.
To compensate for that, the current study compares various empirical evidence of local adaptation, recent selection, and function, to provide a general evaluation for each approach. As a study system, we use the species Arabidopsis thaliana, hereafter mentioned as Arabidopsis. Arabidopsis populations are found all across Eurasia (1001Genomes Consortium, 2016, traversing very different climates (Frachon et al., 2018;Lasky et al., 2012;Mojica et al., 2016;Monroe et al., 2018), and therefore have been thought to be locally adapted. Direct evidence of local adaptation has been observed between North Sweden and Central Italy populations (Ågren & Schemske, 2012;Fournier-Level et al., 2011), in which the fitness of genotypes was measured in a reciprocal transplant experiment at the native sites of each population (Ågren & Schemske, 2012); and recombinant inbred lines of these were used to map quantitative trait loci (QTL) explaining fitness variation (Ågren et al., 2013). Such QTL, can be used to evaluate population genomic methods by testing whether candidate SNPs are enriched within their confidence intervals. Furthemore, given that most of these QTL are of low resolution and span large genomic regions, a tight association between candidate SNPs and logarithm of the odds ratio (LOD) peaks should provide stronger support for an approach. In addition to a tight association to fitness QTL, when the underlying loci exhibit fitness tradeoffs (i.e., when one allele is advantageous in one environment but deleterious in another) we would expect significant allele frequency differentiation between locally adapted populations (Tiffin & Ross-Ibarra, 2014). Allele frequency differentiation is also expected under conditional neutrality (i.e., when one allele is advantageous in one environment but neutral in another) but to a lesser degree (Tiffin & Ross-Ibarra, 2014). GEAs are examined using a large sample of individuals across multiple populations and environments to identify candidate SNPs underlying local adaptation. If a large proportion of these SNPs includes causal or linked loci, then we would expect them to show significant allele frequency divergence between pairs of locally adapted populations and significant enrichment along fitness QTL peaks. Furthermore, if selection is relatively recent, candidate SNPs of both F ST -and GEA-based methods are expected to be within genomic regions that exhibit high linkage disequilibrium (Charlesworth, 2006;Kim & Nielsen, 2004) and a site frequency spectrum that is expected under a recent selective sweep (Kim & Stephan, 2002;Nielsen et al., 2005).
In conjunction with population genomic signatures of selection or significant associations with climate, genetic variation underlying local adaptation is expected to be enriched along sites that are functional and influence fitness. SNPs showing significant associations with climate were found to be significantly enriched among nonsynonymous, but also synonymous variation (Hancock, Brachi, et al., 2011;Lasky et al., 2012). The enrichment among synonymous variation (which largely evolve neutrally) may be the result of linkage disequilibrium due to neutral processes but also background selection (Charlesworth, Morgan, & Charlesworth, 1993) and/or hitchhiking (Gillespie, 2000). A stricter enrichment test will be one that controls for sequence conservation along coding and noncoding sites. Sites that are highly conserved among species (Haudry et al., 2013;Hupalo & Kern, 2013;Miller et al., 2007) are assumed to be under functional constraint and selectively important-that is, due to purifying selection the number of tolerated mutations is limited (Graur, 2016). Therefore, SNPs showing significant evidence of local adaptation across highly constraint sites are more likely to be true positives.
To provide an evaluation of the approaches frequently used to study the genetic basis of local adaptation, we include two methods that identify genotype-by-environment associations while accounting for population structure [genome-wide efficient mixed model (GEMMA; Zhou & Stephens, 2012) and the latent factor mixed model (LFMM; Caye et al., 2019)] and a method [BAYESCAN (Foll & Gaggiotti, 2008)] that identifies SNPs showing higher allele frequency differentiation than expected under various neutral models of evolution. The two GEA methods (GEMMA and LFMM) were applied across a set of 875 Eurasian accessions (1001Genomes Consortium, 2016 and four climate variables covering temperature, precipitation and photosynthetically active radiation (Lasky et al., 2012;Price et al., 2018). The 875 accessions excluded likely laboratory escapees or contaminants (Pisupati et al., 2017) and invasive lines that may reduce the signal of local adaptation (Lasky et al., 2012). BAYESCAN on the other hand was applied to a sample of accessions in North Sweden and South Italy.
After obtaining the results, we addressed the following ques- Addressing the above questions allowed us to assess these methods and consider the best approach to examine candidate genetic variation underlying 20 fitness QTL and evidence tying flowering time to local adaptation. Flowering time is a life-history trait that is thought to play a significant role in local adaptation to climate (Ågren et al., 2017); Dittmar, Oakley, Agren, & Schemske, 2014;Hall & Willis, 2006;Sandring & Agren, 2009;Verhoeven, Poorter, Nevo, & Biere, 2008), and whose genetic basis has been thoroughly studied (Salomé et al., 2011;Sasaki et al., 2015). To re-examine evidence

| Identifying SNPs showing significant associations with climate after accounting for population structure
To identify SNPs showing significant associations with climate, while accounting for population structure we used two prominent methods: GEMMA association (Zhou & Stephens, 2012) (Table S1) that excluded laboratory escapees/contaminants (Pisupati et al., 2017) and accessions from outside the native Eurasian and African range of A. thaliana that may have weaker patterns of local adaptation (Lasky et al., 2012).
After we filtered for biallelic SNPs with minor allele frequency >0.05, we tested association with home climate of ecotype and tested for potential confounding effects of population structure using the software GEMMA (Zhou & Stephens, 2012) with a missingness threshold of 0.05. For the linear mixed model option, we used the Wald test (default) to test the null hypothesis that the mean climate occupied by the two alleles is equal (Lasky et al., 2014). To apply LFMM, we first used the R function "prcomp," to perform a principal component analysis (PCA) and estimate structure in the genotypic data. "prcomp" was applied over 100 random samples of 20,000 polymorphic loci. After determining how many components (K) explained most of the genotypic variance, we applied LFMM ("lfmm_test") in R and calibrated p-values using the "gif" option. To estimate a false discovery rate for both GEMMA and LFMM, we used the "qvalue" function (Storey, 2002) implemented in R.

| Estimating allele frequency differentiation and recent selection across Arabidopsis thaliana populations
Allele frequency differentiation between 25 South Italy and 40 North Sweden accessions (Tables S2 and S3) was estimated using an F ST -based metric implemented in BAYESCAN (Foll & Gaggiotti, 2008) and a simple measure of absolute allele frequency differentiation (|f N.Sweden -f S.Italy |). |f N.Sweden -f S.Italy | was used as an alternative measure, because from prior experience  we noticed that significant F ST s were limited to SNPs in which alternative alleles were fixed between populations (>0.95).
In addition to allele frequency differentiation between Italy and Sweden populations, we also estimated derived allele frequencies across the set of 875 Eurasian accessions (DAF eurasia ) for SNPs that we could estimate the ancestral state. To infer an ancestral base probability, we used Phast suite's Prequel function (Hubisz, Pollard, & Siepel, 2011), and a MAF genome alignment of the species Arabidopsis thaliana (AT) (Berardini et al., 2015), Arabidopsis lyrata (AL) (Hu et al., 2011), Arabidopsis haleri (AH) (Briskine et al., 2017), Capsella rubella (CR) (Slotte et al., 2013), and Neslia paniculata (NP). Alignments were generated with LASTZ (Harris, 2007) and refined for high-confidence orthologs using the pipeline described by (Haudry et al., 2013). In order to remove reference bias, the A. thaliana genome was neutralized in the alignments before calling the ancestral base. The neutral tree used for ancestral infer- To estimate evidence of recent selection at the SNP level, we computed linkage disequilibrium (LD) using the coefficient of correlation (r 2 ). More specifically, using the package "PLINK" (Purcell et al., 2007) LD at a given SNP was estimated as the mean r 2 between it and neighboring SNPs within a 20-kb window (r 2 ). LD was measured using the 65 Italy-Sweden accessions and the 875 Eurasian accessions (r 2 eurasia ). As an additional measure of local adaptation across Italy and Sweden populations, we used SNPs that showed a high |f N.
Sweden -f S.Italy | > 0.70 and LD > 0.19 (hereafter referred to as AFD.LD SNPs). 0.70 and 0.19 represent the 95th percentiles of the respective genome-wide distributions.
For evidence of recent sweeps, we used previously calculated  composite likelihood ratios (CLRs) that were computed using Sweepfinder2 . We focused on CLRs in North Sweden, since preciously estimated CLR signals in other populations were very weak (Huber et al., 2014;Long et al., 2013;Price et al., 2018).

| Determining sites of functional importance
To narrow down SNPs to ones that are more likely to underlie differences in function/expression of protein-coding genes, we focused on cis-regulatory and nonsynonymous variation that was found along sites showing significant selective constraint. We regarded cis-regulatory SNPs as those found within 1 kb upstream from the transcriptional start site of a gene (Pass et al., 2017;Zou et al., 2011), unless these sites were found in genic regions of other genes (in which case they were excluded). To call nonsynonymous variation, we used biallelic sites, we used a publicly available python script (callSynNonSyn.py; archived at https ://github.com/kern-lab/), and gene models downloaded from the TAIR database (TAR10 genome release) (Berardini et al., 2015). To annotate regions showing significant selective constraint across the A. thaliana genome, we used phastCons scores (Siepel et al., 2005) derived using a nine-way alignment of Brassicaceae species from the study by (Haudry et al., 2013). We defined conserved regions as those with a score ≥0.8 over blocks of ≥10 nucleotides.

| Fitness and flowering time QTL underlying Italy and Sweden populations
Quantitative trait loci explaining fitness variation between natural A. thaliana Italy and Sweden populations were retrieved by the study of (Ågren et al., 2013). These 20 fitness QTL (Table S4) were assembled into 6 genetic trade-off QTL (Ågren et al., 2013); however, we treated them as independent given the very long genetic distances between fitness QTL peaks (Table S4). Furthermore, highconfidence QTL explaining flowering time variation were retrieved between these populations were retried from (Ågren et al., 2017) ( Table S5).

| Estimating enrichment of candidate SNPs along functional sites and fitness QTL
To test whether SNPs showing significant evidence of local adaptation are enriched along fitness QTL peaks and among cis-regulatory/ nonsynonymous variation at sites showing significant selective constraint, we used 1,000 circular permutations. For example, when examining whether SNPs showing significant correlations to climate are enriched within a given distance from fitness QTL peaks or cisregulatory/nonsynonymous sites we execute the following steps: 1. Shifting the q-values along the genome while keeping the genomic positions intact. During this step, we take the whole list of q-values ("qvalues") and shift them circularly from a random location ("random_loc") using the following R-code: shifted_qvalues < −qvalues[c(random_loc:total_no_SNPs, 1:( random_loc − 1))]. This shifts the tail-end q-values to the beginning of the list and vice versa.
2. We then choose significant SNPs using our threshold for significance (FDR/q-value <0.1) and estimate the proportion of SNPs within a given distance from a fitness QTL peak (e.g., 100 kb) or the proportion of SNPs at cis-regulatory/nonsynonymous sites showing significant functional constraint.

| Sliding window analysis of chromosomal variation SNPs showing evidence of local adaptation
To detect chromosomal regions with a high proportion of SNPs showing significant evidence of local adaptation, we used a sliding window approach. Specifically, for a window size of 20kb and a step size of 1kb we estimated the ratio of SNPs showing a specific requirement (e.g., |f N.Sweden -f S.Italy | > 0.70 and LD > 0.19) over the total number of SNPs within a 20-kb window.

| Flowering time estimates for A. thaliana Eurasian accessions and candidate flowering time genes
Estimates of flowering time for the 835 Eurasian A. thaliana accessions were downloaded from the study by the 1001 Genomes Consortium (1001Genomes Consortium, 2016

| Constructing rooted gene trees
To build neighbor-joining trees of genes showing significant local adaptation, we downloaded 1:1 orthologs between Arabidopsis thaliana and outgroups Arabidopsis lyrata and Capsella rubella from the Phytozome database (Goodstein et al., 2012) and after aligning the coding sequences with MAFFT (Katoh & Toh, 2008) we used MEGA (Tamura, Stecher, Peterson, Filipski, & Kumar, 2013) to build rooted gene trees.

| Climate-associated SNPs show low allele frequency differentiation and linkage disequilibrium across Italy and Sweden populations
To identify SNPs showing significant associations with climate while accounting for population structure, we applied GEMMA (Zhou & Stephens, 2012) (Ågren & Schemske, 2012)].
In addition to identifying SNPs showing significant associations with climate, we also applied BAYESCAN (Foll & Gaggiotti, 2008)  When examining significant SNPs (FDR < 0.1) after filtering out any significant associations with climate that did not segregate between Italy and Sweden populations, we find a very small overlap between methodologies (Figure 1b). The two GEA methods used (GEMMA and LFMM) showed a 25% overlap in predicted SNPs (Figure 1b), in addition to a negligible overlap with SNPs that showed significant allele frequency differentiation (F ST ) between Italy and Sweden populations according to BAYESCAN (Figure 1b). It is quite surprising that LFMM and GEMMA captured an extremely small proportion of SNPs that showed a significant F ST between Sweden and Italy populations (Figure 1b) given the strong adaptive divergence and evidence of genetic trade-offs underlying these populations (Ågren et al., 2013;Ågren & Schemske, 2012). This, however, may occur because BAYESCAN chooses a very different set of SNPs with a high F ST .
To further examine the evidence of recent selection and genetic differentiation captured by the three methods in Italy and Sweden, we looked at absolute allele frequency differentiation (AFD: |f N. Sweden -f S.Italy |) and linkage disequilibrium (LD) estimated as the SNP level (r 2 (20 kb)-the average r 2 between a SNP and its neighboring SNPs within 20 kb). These measures were examined across different false discovery rates (FDRs), in order to study the relation between these signatures of adaptation and FDR (Figure 1c  ( Figure 1c). A high |f N.Sweden -f S.Italy | was also observed across SNPs that showed a higher FDR (Figure 1c), but that significantly dropped at an FDR > 0.7, where 98.5% of the SNPs were found. On the contrary, |f N.Sweden -f S.Italy | associated with GEMMA and LFMM SNPs was low and near the genomic average across all FDRs (Figure 1c).
The same pattern was observed when examining LD, where climate-associated SNPs showed a very low LD that was repeated across all FDRs (Figure 1d), while SNPs showing a significant F ST according to BAYESCAN showed a significantly higher LD than the genome average (Figure 1d), and as expected, LD decreased as FDR increased. All in all, these results indicate that climate-associated SNPs capture very little evidence of local adaptation and recent selection across Italy and Sweden populations.

| Derived allele frequencies, linkage disequilibrium, and climatic differentiation, across Eurasian populations, significantly differ among candidate SNPs
The weak evidence of local adaptation observed across climate-associated SNPs may occur because they capture selection in other parts of the species range. To examine this possibility, we looked at their derived allele frequencies across the Eurasian range (DAF eurasia ), in association with linkage disequilibrium (r 2 eurasia (20 kb). Furthermore, we compared DAF eurasia and r 2 eurasia between climate-associated SNPs and SNPs that showed significant allele frequency differentiation and LD between Italy and Sweden populations. More specifically, the SNPs included were (a) all significant (FDR < 0.1) LFMM and GEMMA SNPs (i.e., unlike Figure 1b (Figure 2d). This is not surprising since GEA methods identify SNPs that explain more climate variation than the genome average. LD SNPs near QTL peaks was significantly (>95% percentile) higher than expected by chance (Figure 3c). The second largest proportion relative to expectation was observed across LFMM SNPs (Figure 3c). LFMM SNPs within 100 kb of fitness QTL peaks showed a much lower mean AFD (~0.36) and LD (~0.10) than AFD.LD SNPs, which by definition are high in AFD and LD (AFD > 0.70 & LD > 0.19). This is not surprising given the low mean AFD and LD across all significant LFMM SNPs (Figure 1c,d).

| High AFD and LD SNPs show a strong association with fitness QTL peaks and an enrichment at cis-regulatory and nonsynonymous sites showing significant selective constraint
SNPs that were significant according to BAYESCAN, showed a high AFD and LD (Figure 1c,d), but nonetheless they did not show a strong association (Appendix S4) and significant enrichment along fitness QTL, such as, high AFD.LD SNPs (Figure 3c). This may occur because significant BAYESCAN SNPs cannot capture all selection events underlying fitness QTL.
To examine where AFD.LD and BAYESCAN differed along genetic trade-off QTL that were built using 20 fitness QTL (Ågren et al., 2013), we scanned for regions showing a high proportion of SNPs that were common between the two sets (AFD. LD and BAYESCAN) and SNPs that were unique to the AFD.LD set. As shown in Appendix S5, genetic trade-off QTL that included regions that show high composite likelihood ratios for recent sweeps we see a high proportion of SNPs that are common to both approaches; on the other hand, along genetic trade-off QTL where evidence of recent sweeps is almost absent (Appenices S6 and S7) we only observe SNPs that were unique to the AFD.LD set. only along conserved nonsynonymous sites (Figure 3a). For methods that showed an enrichment at conserved sites (cis-regulatory and/or nonsynonymous sites), average LD in Eurasia (r 2 eurasia ) for BAYESCAN, LFMM, and AFD.LD SNPs was 0.27, 0.12, and 0.29, respectively. The average LD across LFMM SNPs was approximately the same as the genome-wide mean at conserved cis-regulatory and nonsynonymous sites (~0.12), while for BAYESCAN/AFD.LD SNPs was twice as high. Similarly, to the whole set of SNPs (Figure 2a), LFMM SNPs also showed a low mean DAF eurasia (~0.13).

| Genes underlying fitness QTL and controlling flowering time show significant evidence of local adaptation
Given the significant evidence of local adaptation and function underlying AFD.LD SNPs, we used them to detect potential genes that may underlie fitness QTL (note: We only focused on variants within 100 kb of their peaks), in addition to examining evidence of local adaptation underlying a list of ~170 genes affecting flowering time (Table S6). To further narrow down on SNPs that are more likely to underlie the fitness QTL examined, we only considered cis-regulatory/nonsynonymous variation at conserved sites that segregated between the parents used to derive the RILs (Ågren et al., 2013).
SNPs between the parental genomes were called in a previous study . Our analysis resulted in 24 genes within 100 kb of fitness QTL peaks and spanning three genetic trade-off QTL (2:2, 4:2, and 5:5) (Appendix S8). Many of these were involved in interesting biological processes such as response to different abiotic stress factors and the abscisic-acid signaling pathway which is important in abiotic stress response (Tuteja, 2007) (Appendix S8). Among these genes, two of them (AT4G33360 (FLDH) and AT4G33470 (HDA14)) showed strong expression G×E interactions [G×E interactions were identified in a previous study ] when Italy and Sweden plants were grown under cold acclimation conditions (4°C) for two weeks (Gehan et al., 2015). Interestingly, FLDH is a negative regulator of the abscisic-acid signaling pathway (Bhandari, Fitzpatrick, & Crowell, 2010). As shown in Figure 4, this gene was within a region of a genetic trade-off QTL that showed a high proportion of AFD.
LD SNPs. Furthermore, expression of FLDH under control and cold acclimation conditions was significantly lower in Sweden than Italy plants ( Figure 4).
To examine whether genes affecting flowering time (Table S6) show significant evidence of local adaptation, we tested whether the number of such genes containing cis-regulatory and/or nonsynonymous SNPs with a high AFD (Figure 5a), and a high AFD and LD (AFD.LD) (Figure 5b) was significantly higher than expected by F I G U R E 3 Testing for enrichment of candidate SNPs along QTL peaks explaining fitness variation between Italy and Sweden populations (Ågren et al., 2013) and cis-regulatory/nonsynonymous sites showing significant selective constraint across nine Brassicaceae species (Haudry et al., 2013). (a) Comparing the observed proportion of candidate SNPs (red lines) within 100 kb of fitness QTL peaks, to a distribution of expected proportions derived using circular permutations. Significance was set at the 95th percentile.  Table S5).
According to FlrT-5:4, the Sweden genotype was associated with longer flowering time in both Italy and Sweden (Ågren et al., 2017). In conjunction, with its overlap to a genetic trade-off QTL (Ågren et al., 2017), it indicates a possible role in fitness trade-offs. Studies have attributed flowering time variation within FlrT-5:4 to VIN3 (Ågren et al., 2017;1001Genomes Consortium, 2016. Although it may be an additional candidate, we did not find any significant genetic differentiation and selection along coding and cis-regulatory sites of VIN3.

| D ISCUSS I ON
In the quest to study the genetic basis of local adaptation using genome-wide associations with environment, linear mixed models have emerged as a powerful tool given their ability to account for population structure while testing for significant associations (Caye et al., 2019;Kang et al., 2010Kang et al., , 2008Yu et al., 2006;Zhou & Stephens, 2012). Although they provide a robust statistical framework, the current study shows that such approaches may significantly limit our ability to understand the polygenic basis of local adaptation.
Both GEA methods (GEMMA and LFMM) resulted in SNPs that showed poor associations with QTL explaining fitness variation, in addition to low genetic differentiation and evidence of recent selection across locally adapted populations. The poor performance of GEA methods when examining populations isolated by distance has also been shown when using simulations (Lotterhos & Whitlock, 2015). In fact, F ST -based approaches outperformed GEA methods in such instances (Lotterhos & Whitlock, 2015). In the current study, we show that SNPs exhibiting a significantly high F ST according to BAYESCAN, capture higher population genomic evidence of recent selection but fail to show a strong association with fitness QTL and a significant enrichment along regions with a high LOD score.
Using a more lenient FDR cutoff may capture some of the missing SNPs across fitness QTL that do not contain strong CLRs for recent sweeps (Appenices S6 and S7).
SNPs that show extreme AFD and LD are likely to contain many false positives, especially if no other information is considered.
These genomic signatures, however, show a promising future in identifying recent local adaptation within a statistical framework (Kemppainen et al., 2015). Such an approach can avoid several com-  (Bhandari et al., 2010), found within genetic trade-off QTL 4:2 (Ågren et al., 2013) and 100 kb from a fitness QTL peak (red and blue arrows represent QTL where the Sweden genotype had lower fitness in Italy and higher fitness in Sweden, respectively). The region including FLDH showed a high proportion of high AFD and LD (AFD.LD) SNPs, and furthermore, in Italy and Sweden plants FLDH showed strong expression GxE interactions under control (22°C) and cold acclimation conditions (4°C) for two weeks (FPKM: fragments per kilobase million) 2014; Lotterhos & Whitlock, 2015). An example of such populations is ones found in Northern and Southern Sweden, which follow a two-population island model (Huber et al., 2014 (Hancock, Brachi, et al., 2011;Hancock, Witonsky, et al., 2011;Lasky et al., 2012) or the parallel occurrence of low-frequency loss-of-function mutations showing significant associations with climate  has been interpreted as evidence of adaptation. Although such signals may represent instances of adaptation, they can also be explained by neutral evolution, in which relaxed selection across specific climates results in the enrichment of F I G U R E 5 Genes known to affect flowering time show significant evidence of local adaptation along putative functional sites. (a) The number of flowering time genes containing cis-regulatory/nonsynonymous SNPs showing a high AFD was significantly higher (>95th percentile) than expected by chance. Distribution of random numbers was derived using circular permutations. (b) The number of flowering time genes with a high AFD and LD cis-regulatory/nonsynonymous SNPs was also significantly higher than expected by chance. (c) PIF3 is a phytochrome interacting factor that has been found to affect flowering time (Oda et al., 2004) that was found underlying a region along chromosome 1 that showed the largest composite likelihood ratios (CLRs) for recent sweeps in Sweden, and windows with a high proportion of AFD.LD SNPs. A rooted phylogeny of the PIF3 coding region indicated that Eurasian accessions sharing the same allele as the Sweden parent (blue dot) show significantly higher flowering time than accessions sharing the same allele as the Italy parent (red dot). (d) COL5 is another gene that has been found to affect flowering time (Hassidim, Harir, Yakir, Kron, & Green, 2009) and in which Eurasian accessions show significant genetic differentiation and segregation in flowering time. This gene is also found within previously identified flowering time QTL (FlrT-5:4) (Ågren et al., 2017) in which the Sweden genotype was associated with longer flowering time in both Italy and Sweden independent loss-of-function mutations or nonsynonymous variation (Flowers, Hanzawa, Hall, Moore, & Purugganan, 2009;Zhen, Dhakal, & Ungerer, 2011;Zhen & Ungerer, 2008). In the context of local adaptation, these may represent instances of conditional neutrality, where in one environment expressing the gene has no significant impact on fitness. Some ways that could provide further support as to whether a recent loss-of-function mutation is adaptive are to compare LD or extended haplotype homozygosity (Sabeti et al., 2002) between individuals that have a loss-of-function mutation and ones that do not.
Among the genomic signatures of local adaptation examined, SNPs showing a high absolute allele frequency differentiation (AFD) and linkage disequilibrium (LD) between Italy and Sweden populations showed the strongest evidence of local adaptation and were enriched among nonsynonymous/cis-regulatory variation at sites showing significant selective constraint. Using these SNPs, we identified a list of candidate genes underlying fitness QTL. One of these was FLDH, a negative regulator of abscisic-acid signaling (Bhandari et al., 2010), that showed strong GxE interactions between Italy and Sweden plants under cold acclimation conditions. Abscisic-acid signaling is known to play an important role in abiotic stress response (Tuteja, 2007), with many studies supporting a key role in local adaptation to climate (Kalladan et al., 2017;Keller, Levsen, Olson, & Tiffin, 2012;Lasky et al., 2014;Ristova, Giovannetti, Metesch, & Busch, 2018). In addition to abscisic-acid signaling, our study provides further support for the important role of flowering time in local adaptation to climate. Among a list of genes that were experimentally shown to affect flowering time, we identified three genes (PIF3, FIO1, and COL5) that showed significant evidence of local adaptation and selective constraint along nonsynonymous sites. FIO1 was previously shown to contain SNPs that showed significant associations with flowering time among natural Swedish lines (Sasaki, Zhang, Atwell, Meng, & Nordborg, 2015) and COL5 was located within a QTL that explain flowering time variation among Sweden and Italy recombinant inbred lines (Ågren et al., 2017). Finally, PIF3, a transcription factor that interacts with phytochromes (Soy et al., 2012), has been implicated in multiple biological processes including early hypocotyl growth (Monte et al., 2004), photomorphogenesis (Dong et al., 2017), flowering time (Oda, Fujiwara, Kamada, Coupland, & Mizoguchi, 2004), and regulation of physiological responses to temperature (Jiang et al., 2017).
Interestingly, PIF3 was found within a large region that showed significant evidence of local adaptation. Regions of high divergence may involve a single causative variant, or a group of linked genes that interact with PIF3 and were under selection because they contributed to building an advantageous phenotype (Barton & Bengtsson, 1986;Yeaman & Whitlock, 2011). Although such "Islands of high divergence" can be the result of local adaptation, they can also be formed through nonadaptive processes (Pennisi, 2014), and therefore, caution should be exercised when drawing any conclusions.
When ignoring selective constraint, we identify a list of addition flowering time genes showing significant evidence of local adaptation along nonsynonymous/cis-regulatory sites. Genes such as SVP and MAF3 were previously associated with flowering time variation among natural Arabidopsis accessions (Caicedo, Richards, Ehrenreich, & Purugganan, 2009;Sasaki et al., 2015), and MAF3 showed strong allelic variation along a multivariate climate gradient (Lasky et al., 2012).
Although adaptation may involve sites that are not deeply rooted and/ or under strong selective constraint, including additional plant genomes when estimating sequence conservation across species may increase our power to detect selectively important regions. As shown by studies examining adaptation in species ranging from bacteria (Maddamsetti et al., 2017) to birds (Sackton et al., 2019), addressing selective constraint can improve our understanding of its genetic basis.
To sum up, the current study identifies candidate genes and life-history traits that may underlie adaptation of Arabidopsis populations to local environments. Furthermore, it shows that understanding organismal adaptation to local environments is a very complex venture, where several lines of evidence are needed to obtain a comprehensive and well-supported picture.

ACK N OWLED G M ENTS
We would like to thank the reviewers, Brenna R. Forester, and John McKay for suggestions and discussions that significantly improved the manuscript.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflict of interest.

AUTH O R CO NTR I B UTI O N S
Nicholas Price designed research, performed research, analyzed data, and wrote the paper. Lua Lopez and Adrian E. Platts performed research and reviewed the paper. Jesse R. Lasky reviewed the paper.

DATA AVA I L A B I L I T Y S TAT E M E N T
In our dryad submission, we included: