Temperature, rainfall and wind variables underlie environmental adaptation in natural populations of Drosophila melanogaster.

While several studies in a diverse set of species have shed light on the genes underlying adaptation, our knowledge on the selective pressures that explain the observed patterns lags behind. Drosophila melanogaster is a valuable organism to study environmental adaptation because this species originated in Southern Africa and has recently expanded worldwide, and also because it has a functionally well-annotated genome. In this study, we aimed to decipher which environmental variables are relevant for adaptation of D. melanogaster natural populations in Europe and North America. We analysed 36 whole-genome pool-seq samples of D. melanogaster natural populations collected in 20 European and 11 North American locations. We used the BayPass software to identify single nucleotide polymorphisms (SNPs) and transposable elements (TEs) showing signature of adaptive differentiation across populations, as well as significant associations with 59 environmental variables related to temperature, rainfall, evaporation, solar radiation, wind, daylight hours, and soil type. We found that in addition to temperature and rainfall, wind related variables are also relevant for D. melanogaster environmental adaptation. Interestingly, 23%-51% of the genes that showed significant associations with environmental variables were not found overly differentiated across populations. In addition to SNPs, we also identified 10 reference transposable element insertions associated with environmental variables. Our results showed that genome-environment association analysis can identify adaptive genetic variants that are undetected by population differentiation analysis while also allowing the identification of candidate environmental drivers of adaptation.

INTRODUCTION

Understanding how organisms adapt to different environments is a major goal in evolutionary biology (Hoban et al., 2016; Nelson et al., 2019). The genetic basis of adaptative traits has been studied in several organisms, such as lactase persistence (Tishkoff et al., 2007) and skin colour in humans (Norton et al., 2007), and dark colour in the peppered moth Biston Betularia (Van't Hof et al., 2016) among many others. Genome‐wide studies aimed at elucidating the genetic basis of environmental adaptation have also been conducted in several species such as plants (Flood & Hancock, 2017), bacteria (Gorter et al., 2016) and Drosophila (Rech et al., 2019). However, knowledge on the specific environmental variables driving these adaptations lags behind.

In the past few years, the availability of whole genome sequences as well as the development of different analytical tools, have facilitated the performance of genome‐environment association (GEA) analyses. GEA analyses are useful approaches to identify the genetic variants and the environmental factors that are involved in the adaptive processes. Combining outcomes from GEA analysis with classical population genome‐wide selection scans, such as those based on differentiation statistics, may help to link the genetic variants underlying local adaptation with their environmental drivers (Ahrens et al., 2018; Hoban et al., 2016). These analyses have already been applied to several species including Plant, Chordata, and to a lesser extent Arthropoda, Mollusca, Cnidaria, Echinodermata, and Nematoda (Ahrens et al., 2018). However, there are still important limitations behind GEA analyses. One of the main drawbacks is the difficulty in distinguishing the patterns associated with demographic processes from those that are the consequence of selection (reviewed in Rellstab et al., 2015). A second limitation is related to the choice of environmental variables to include in the analysis. Prior selection of the most relevant variables for any particular GEA analysis is complicated since some previous knowledge about which variables may be relevant in the adaptation process is needed. Indeed, most environmental variables used in GEA studies are related to temperature and precipitation, while other variables such as solar radiation, daylight hours, evaporation, and wind, that could also play a role in adaptation are not widely used. Solar radiation, and more specifically UV‐B radiation, could be relevant as DNA damage responses are known to play a role in adaptation of several species such as birds, insects or fungi (Körner, 2007; Svetec et al., 2016; Wu et al., 2019; Zhou et al., 2020). Daylight hours is related to the circadian rhythm, which for example is known to play a role in Drosophila behavioural adaptation to high latitudes (Helfrich‐Förster et al., 2020). Evaporation is involved in organism thermoregulation and response to desiccation stress (Ferveur et al., 2018; Rajpurohit et al., 2018; Smit et al., 2018). Finally, wind direction is involved in plant adaptation by modifying pollen flow and therefore changing the spatial genetic structure (Balkenhol et al., 2017; Gardiner et al., 2016; Wang et al., 2016) and in the case of insects, antennae and specifically the Johnston's organ, are directly involved in neuron response to wind (Fuller et al., 2014; Patella & Wilson, 2018).

In addition, climate variables such as temperature and precipitation can be highly correlated (Lotterhos et al., 2018). Relationship between explanatory variables, i.e., multicollinearity, compromises the results of multivariate regression analysis (Kim, 2019). Multicollinearity could yield unreliable regression parameter estimation, magnitude and sign of regression, which impedes the assessment of the relative importance of the explanatory variables (Sokal & Rohlf, 2013). This problem may be overcome by using synthetic variables obtained via principal component analysis (PCA) of the environmental variables of interest. However, using PCs based on climate variables may lead to a limited interpretation of the environment drivers of selection. The PCs will represent the environmental variables that covary the most, but this may not coincide with the combination of variables that drive divergent selection and local adaptation (Houle et al., 2002; Lotterhos et al., 2018).

Recently developed software such as the BayPass package (Gautier, 2015), have overcome some of the limitations of the GEA analyses mentioned above. On one hand, this software identifies those genetic variants with statistically different allele frequencies between populations and those associated with environmental variables, while taking into account the covariance between population allele frequencies due to, for instance, the joint demographic history of the samples analysed (Gautier, 2015). On the other hand, this software includes different modules based on different models: a single‐covariate regression model where the association is estimated for each covariate, and a multiple‐covariate regression model where the association is estimated for several covariates assumed to be orthogonal (Gautier, 2015).

Drosophila melanogaster is a valuable model organism to study environmental adaptation. This species originated from Southern Africa and has recently expanded worldwide colonizing a wide range of environmental conditions (12,000–19,000 years ago; Arguello et al., 2019; Pool et al., 2012; Sprengelmeyer et al., 2020). In addition, this species offers many key advantages as it has a small and well‐annotated genome which facilitates the identification of putatively adaptive loci (Mohr et al., 2014), as well as a short lifecycle implying many generations in short periods of time (15 generations per year in nature; Pool, 2015). Past studies carried out with North American and Australian D. melanogaster populations have already shown clinal and seasonal genetic patterns suggesting that this species could be a good model to study environmental adaptation (Bergland et al., 2014, 2016; Fabian et al., 2012; Hoffmann & Weeks, 2007; Kolaczkowski et al., 2011; Machado et al., 2019). Indeed, this species has already been studied in other continents such as Europe, where clinal patterns and correlations between genetic variants and environmental variables have also been identified (Kapun et al., 2020; Lerat et al., 2019).

In this study, we combined genome scans for adaptive differentiation and whole‐genome GEA analysis using pool‐seq data available for 36 samples of D. melanogaster, representative of the genetic diversity across the European continent (n = 20 locations) and across a latitudinal cline in eastern North America (n = 11 locations). We focused on these two continents because they have an approximately similar range of climatic conditions (mostly temperate climates) and they were both recently colonized, which allow us to focus on short‐term evolutionary events. Our threefold aims were to characterize: (i) to which extent environmental variables contributed to adaptive differentiation in D. melanogaster; (ii) which climatic variables, namely temperature, rainfall, evaporation, solar radiation, wind, daylight hours and soil type, may be contributing to this environmental adaptation; and (iii) to which extent the observed signals were parallel across two different geographic areas: Europe and North America.

MATERIALS AND METHODS

Data sets

European pool‐sequencing samples were obtained from the 2014 DrosEU data set (Kapun et al., 2020). We discarded 18 out of the 48 samples available, for which Tajima's D was very low (Tajima's D < –0.2; Kapun et al., 2020). For some locations, samples were collected several times across 2014. When several samples were available for the same season, summer or fall, we only included the earliest collected sample in the analysis. Thus, overall, we analysed 25 samples from 20 different locations (Figure 1). To perform the analysis, we created three data sets including only one sample per location (Table S1): Europe (20 samples), Europe Summer (14 samples), and Europe Fall (10 samples). Note that for the Europe data set, when samples of both seasons were available for the same location, we only included the summer sample. Average sequencing coverage among samples ranged from 25 to 190X (Table S1). VCFs are available at http://hdl.handle.net/10261/180630.

FIGURE 1

Drosophila melanogaster samples used in this study. Samples were collected across Europe and North America east coast (see Table S1) in four main climate zones and seven subclimate zones (depending on precipitation and level of heat), according to the Köppen‐Geiger climate distribution (Kottek et al., 2006) [Colour figure can be viewed at wileyonlinelibrary.com]

Eleven North American pool‐sequencing samples collected from 2003 to 2014 were obtained from Machado et al., (2019) sampled in eleven different locations in the North American East coast (Figure 1; Table S1). We focused in these samples because clinality has been detected in previous studies (Bergland et al., 2014; Fabian et al., 2012; Figure 1; Table S1). VCFs are available at https://datadryad.org/stash/share/rHMqJSiXuGX12eBYyPvKE_Ng1b‐FMTrLLnmegosbQ74.

In addition to SNPs, we also included in our analysis transposable element (TE) insertions. We estimated population frequencies for 1,630 euchromatic reference TE insertions (Rech et al., 2019) using T‐lex3 (Bogaerts‐Márquez et al., 2019). For each data set, we only analysed those insertions that were polymorphic in at least one of the populations (403 TE insertions; Table S2).

Environmental variables

We downloaded environmental data from four different sources: WorldClim (Fick & Hijmans, 2017), Copernicus (Hersbach et al., 2020), the US Naval Observatory (https://www.usno.navy.mil/USNO/astronomical‐applications/data‐services/data‐services) and The Astronomical Data Portal UK Hydrographic Office (http://astro.ukho.gov.uk/). From WorldClim, we used the 19 Bioclimatic variables, which are derived from the monthly temperature and rainfall values from the 1970–2000 time range. We used the R package raster (v.2.6–7) for downloading this data. In addition, we also used year‐specific environmental variables from the year previous to the collection date of each sample. For the year‐specific variables, we used ERA5 "hourly data on single levels from 1979 to present" database from Copernicus, to obtain data on temperature (2 m temperature), rainfall (Total precipitation), evaporation (Evaporation), solar radiation (Clear‐sky direct solar radiation at surface), wind (10 m u‐component of wind and 10 m v‐component of wind), and soil type. This data was downloaded as GRIB files and parsed using ecCodes package (v.2.8.2). Finally, daylight hours were obtained from the US Naval Observatory for the European data sets, and from The Astronomical Data Portal UK Hydrographic Office for the North American data set. For the year‐specific variables, similar variables as the ones in WorldClim were constructed with the data from the year previous to the collection date using different python scripts (v.2.7.12). In total, we analysed 59 environmental variables related to temperature (11), rainfall (8), evaporation (11), solar radiation (7), wind (14), daylight hours (7) and soil type (1) (Table S3). We tested whether those variables that are common between WorldClim and Copernicus (i.e., temperature and rainfall) but that were obtained from different time ranges (average data from 1970 to 2000 and year‐specific data from the year before the collection date) were significantly correlated using a Spearman correlation test (ρ > 0.8) (Table S3A). Most of the eight rainfall related variables were not correlated (six in Europe, six in Europe Summer, five in Europe Fall, and eight in North America), while most of the 11 variables related to temperature in Europe and NA data sets were correlated (seven in Europe, eight in Europe Summer, nine in Europe Fall, and seven in NA; Tables S3A–B). For variables that were not correlated, we included the two values of the same environmental variable in the analysis. When the correlation was significant, we performed the analysis using the variable corresponding to the WorldClim database.

In order to study the correlation among environmental variables, we performed Spearman rank correlation test for each pair of variables in each data set using the R function cor.test() and R package corrplot (v.0.84) (Table S3C). We considered and reported as strong correlations those with Spearman's ρ > 0.8. We found that solar radiation and daylight hours in Europe and in North America were highly correlated. In North America, but not in the European data sets, we also found that temperature, solar radiation, and daylight hours variables were highly correlated (Figure S1A–D). We also performed a PCA of the environmental variables for each of the four data sets analysed using stats (v.3.5.2) and plotnScree function in nFactors (v.2.4.1) R packages (Table S3D1–4). We found that temperature and solar radiation explained most of the variation of the PC1 both in Europe and Europe Summer data sets. In the Europe Fall data set solar radiation and daylight hours explain most of the PC1 variation while temperature and rainfall explain most of the PC1 variation in North America (Table S3D1–4).

Whole‐genome scans for adaptive differentiation

Whole genome‐scans for adaptive differentiation were performed using BayPass (v.2.1) (Gautier, 2015). The model underlying BayPass accounts for the correlation structure among allele frequencies and allows identifying putative genetic variants subjected to adaptive differentiation based on the XtX statistic (Günther & Coop, 2013). This method, as demonstrated in Gautier (2015), has several critical advantages: (i) it explicitly and efficiently accounts for the confounding factors of the shared demographic history (via the covariance matrix ‐Ω); (ii) it makes no simplifying assumptions about the underlying demographic model; and (iii) it can explicitly model pool‐seq data (via a binomial likelihood under the so‐called poolseqmode) to account for the extra‐variance introduced when sequencing pools of DNAs that are not individually identifiable (which basically prevents from distinguishing reads that are identical because they were obtained from the same sequence or from two distinct but identical sequences) (see Gautier, 2015 for a more detailed explanation). The genotyping input files contain the read count data (reference and alternative) per site and per sample. For SNPs, this information was obtained from the VCF files, while for the transposable elements (TEs) the information was obtained using T‐lex: the absent read count information was used as the alternative read count, and the present read count information as the reference read count (Bogaerts‐Márquez et al., 2019; Fiston‐Lavier et al., 2015). To generate the input files for the North America samples, VCFs were parsed using the poolfstat (v.1.1.1) R package (Hivert et al., 2018). For the three European data sets, VCFs were parsed using python scripts, and bash and awk command lines. TE frequencies were added to the data using python scripts. Invariant and nonbiallelic positions were removed from each data set. We ran BayPass for autosomes (2L, 2R, 3L and 3R) and X chromosome separately because the autosome and X‐linked variants have different haploid sample sizes as samples were obtained from male flies, and more importantly autosomes and X‐chromosome have different demographic histories (e.g., Clemente et al., 2018).

We ran BayPass under the core model for the computation of the XtX genetic differentiation statistics for each data set separately (Europe, North America, Europe Summer and Europe Fall). As the number of SNPs in the autosomes in the four data sets was large (Table 1), we used a subsampling strategy as in Gautier et al. (2018), dividing each data set into 50 subdata sets containing only one SNP every 50 SNPs. We run the 50 subdata sets in parallel, thus all the SNPs available were used for the analysis. This strategy allows a more efficient analysis as it requires less computational time because each of the 50 pseudo‐independent files are run in parallel. Note that the SNPs and TEs in each pseudo‐independent file have low level of background linkage disequilibrium (LD). Three independent runs (using the option ‐seed) were performed for each data set. To check for consistency the results of the three independent runs were evaluated using the Forstner and Moonen Distance (FMD) (Förstner & Moonen, 2003) between pairs of covariance matrices (Ω) with the R function fmd.dist() (provided within the BayPass package). We compared the covariance matrices among the 50 different subdata sets, and among the three different seeds runs. We found that FMD was low for all comparisons (Table S4A). Consistency was additionally tested comparing the posterior means of the two parameters α and β of the Beta distribution of the estimated population allele frequencies (Table S4B).

TABLE 1

Summary of the four data sets used in this analysis: three European and one North American data set

Data set	No of populations	Autosomes		X chromosome		Total
		SNPs	TEs	SNPs	TEs	SNPs	TEs
Europe	20	2,846,701	249	119,228	53	2,965,929	302
North America	11	2,147,276	280	291,632	64	2,438,908	344
Europe Fall	10	2,663,700	227	115,147	49	2,778,847	276
Europe Summer	14	2,725,176	222	117,332	42	2,842,508	264

Prior to further analysis, we removed SNPs and TEs with very low allele frequency (MAF <0.01) based on the mean of the posterior distribution of the frequency of the reference allele across populations for each site (included in column M_P of the BayPass output summary_pi_xtx.out).

To obtain a calibrated estimator of the XtX statistic, we relied on the XtX* estimator recently described in (Olazcuaga et al., 2020). We further derived bilateral p‐values assuming that XtX* follows a Chi‐square distribution with npop degrees of freedom under the null hypothesis of neutral differentiation (Olazcuaga et al., 2020). To control for multiple testing, we estimated the associated q‐value with the R package qvalue (v 3.9) (Storey & Tibshirani, 2003). For further analysis, we focused on the SNPs and TEs with the most highly significant XtX* scores (top 0.05% and q‐value ≤0.05).

Genome‐environment association analysis

The so‐called BayPass STD model extends the previous analysis to evaluate association of the genetic variant allele frequencies with population‐specific covariates. We ran this model with the environmental variables previously described for each data set, and default options except for the ‐scalecov option that was used to scale each covariate. As we did with the previous model, we run the four data sets independently, as well as autosomes and X chromosome separately. For the autosome data sets, we also used the subsampling strategy mentioned above. Three independent runs were performed for each data set (using the option ‐seed).

The support for association between the genetic variants and the environmental variables was assessed using a Bayes factor (BF) measured in deciban (dB) units and estimated with an importance sampling algorithm from the MCMC samples (Coop et al., 2010; Gautier, 2015). More specifically, we used as an estimate the median BF among the three independent MCMC runs. We discarded SNPs and TEs present at very low allele frequency (MAF <0.01) in both observed and simulated data (see below). We considered a BF threshold of 20 dB (i.e., “decisive evidence” according to the Jeffreys’ rule [Jeffreys, 1961]) as evidence for association between an environmental variable and a TE and an even more stringent threshold of 30 dB for SNPs to limit the number of false positives (SNPs being far more numerous than TEs). We evaluated the false positive rates (FPR) associated with these thresholds based on the analysis of pseudo‐observed data sets (PODs) generated using the same parameters as for the observed data sets according to the approach described in Gautier (2015). Briefly, the rationale of this approach is to provide an empirical null distribution of the BF statistic, i.e., neutrally evolving SNPs are simulated under the generative model parameterized with the matrix Ω, which is estimated on the real data to summarize the joint demographic history of the populations (and to capture its effect on the neutral covariance structure on population allele frequencies). The estimated FPR for the 20 BF and 30 BF thresholds ranged from 0% to 2.40% and 0% to 0.33%, respectively for the association tests with the different covariates (Table S5). When a SNP or a TE was significantly associated with more than one environmental variable, we considered as the primary variable the one with the relative highest BF compared to the 99.9% of the POD distribution (Table S5), which usually coincides with the absolute highest BF value.

Analysis of candidate genes

We identified the genes where the significant SNPs were located using bedtools intersect (v.2.27.1) and the D. melanogaster FlyBase reference genome annotations v6.12 and v5.50 for the European and the North American data sets, respectively (Thurmond et al., 2019). We also identified significant SNPs located in gene regulatory region (<1 kb upstream of genes) (Hoskins et al., 2011) using SnpEff (v.4.3) (BDGP5.75 data base for North American and BDGP6.86 for European data sets). For TEs, we also used FlyBase annotations to check whether they were located <1 kb upstream or downstream of a gene, inside a gene, or in intergenic regions.

We performed a gene ontology (GO) enrichment analysis of the candidate genes using DAVID (v.6.8). A cluster was considered to be significant when the enrichment score was higher than 1.3 (Huang et al., 2009).

To test if there was a significant overlap of candidate genes between Europe and North America (624,069 shared SNPs corresponding to 15,944 genes), we use the SuperExactTest R package (v.1.0.7).

RESULTS

Development and signalling underlie population differentiation in Europe and North America

To characterize the patterns of genetic differentiation in European and in North American D. melanogaster natural populations, we ran the BayPass core model in two data sets containing 20 and 11 populations, respectively (Figure 1, Table 1 and Table S1). Samples were collected from seven climate types distributed in four climate zones (Figure 1, Table S1; Kottek et al., 2006). For some European populations, we had samples collected in summer and fall. Thus, in addition to the whole European data set, we also analysed the summer (Europe Summer data set) and the fall (Europe Fall data set) samples separately (Kapun et al., 2020).

We first analysed the distribution of the SNPs that showed significant population differentiation patterns across chromosomes (Figure 2). We tested whether any of the four main cosmopolitan inversions described in D. melanogaster (In(2L)t, In(2R)NS, In(3L)P, and In(3R)P) were enriched for significantly differentiated SNPs by comparing the SNPs located inside each insertion with the rest of the chromosome (Corbett‐Detig & Hartl, 2012; Figure 2, Table S6A). In the European and European Fall data sets, inversion In(2L)t was enriched for significantly differentiated SNPs (Fisher exact test p‐value <0.001). This inversion shows a strong frequency gradient in European populations ranging from 2% to 50% (Table S6E–F; Kapun et al., 2020). On the other hand, in the North American data set inversion In(3R)P was enriched for significantly differentiated SNPs (Fisher's exact test p‐value <0.001), which also shows a strong frequency gradient in North American populations ranging from 0% to 41% (Tables S6E–F; Kapun et al., 2016; Kapun & Flatt, 2019). Our results are consistent with previous analyses that found that these two inversions show latitudinal and/or seasonal clinal patterns mainly in Australia and North America (Kapun et al., 2016; Kapun & Flatt, 2019).

FIGURE 2

Genome‐wide distribution of the XtX* values (population differentiation) associated with single nucleotide polymorphisms (SNPs) and transposable elements (TEs) in the four data sets. Nonsignificant SNPs and TEs are plotted in grey and black, respectively. Significant SNPs located inside and outside of inversions are plotted in dark blue and blue, respectively. Significant TEs are plotted in red. Genes for the five most significant SNPs for each data set are highlighted [Colour figure can be viewed at wileyonlinelibrary.com]

We considered genes with at least one significant SNP located in the gene body region or in their 1 kb upstream regions to be candidates for adaptive differentiation, i.e., to be subjected to selection (Table S7; see Materials and Methods; Hoskins et al., 2011). Overall, we identified 1,300 candidate genes. Among our candidates, we found genes previously known to play a role in adaptation, such as cpo, involved in reproductive dormancy (Cogni et al., 2014; Schmidt et al., 2008), sgg involved in circadian rhythm (Rand et al., 2010; Wolf et al., 2007), mth involved in longevity and stress response (Schmidt et al., 2000), and Ace, involved mainly in insecticide resistance (Fournier et al., 1992; Menozzi et al., 2004). We also found other interesting genes, which have been previously reported as candidates in North America but not in Europe, such as obst‐F, which has been suggested to be involved in longitudinal clinality (Table 2; Campo et al., 2013). obst‐F is involved in the cuticle formation, which acts as a barrier between the fly and the environment protecting it from insecticides, solar radiation, and desiccation (Balabanidou et al., 2018; Behr & Hoch, 2005; Ferveur et al., 2018; Rajpurohit et al., 2018). Indeed, most of the SNPs with the highest differentiation score were located in genes that are candidate for several stress responses, such as oxidative and starvation stress response, and behavioural phenotypes (Table 2).

TABLE 2

Candidate genes showing the most significant population differentiation patterns

Gene	SNP location	Data set	XtX*	Phenotype
BORCS6/klar	Gene body/Upstream	NA	89.63	–/Alcohol, Starvation
Gale	Upstream	NA	88.85	Aggressiveness; Diapause; Immunity; Starvation
Ca‐alpha1T	Gene body	NA	88.40	Olfactory
tok	Gene body	NA	83.40	Circadian; Starvation
lncRNA:CR43314	Upstream	NA	77.60	–
CG6951	Upstream	EuS	171.38	Alcohol; Oxidative
Argk	Upstream	EuF	117.21	Immunity; Starvation
capu	Gene body	EuF	105.60	Alcohol, Circadian behavior, Oxidative, Xenobiotic
ed	Gene body	EuF	104.07	Olfactory, Oxidative
CG7102	Gene body	EuF	116.21	–
Ace	Gene body	Eu	197.83	Diapause, Insecticide resistance, Olfactory, Starvation
	Gene body	EuS	166.55
obst‐F	Gene body	Eu	213.90	–
	Gene body	EuS	194.73
CG17233	Gene body	Eu	228.63	–
	Gene body	EuS	172.06
Clc	Upstream/ Gene body	Eu	192.62	Starvation
		EuS	171.38
RapGAP1	Gene body	Eu	211.90	Aggressiveness
	Gene body	EuF	110.64

For each data set, top 5 genes with SNPs located in the gene body or upstream region (< 1kb) with the highest significant XtX* values and their associated phenotype (see Table S12).

Abbreviations: Eu, Europe; EuF, Europe Fall; EuS, Europe Summer; NA, North America.

To identify which biological processes underlay the population differentiation in the four data sets analysed, we performed a GO term enrichment analysis (Table 3 and Table S8). Both in Europe and in North America, the most significant clusters were related to development and signalling, suggesting that similar biological processes have been involved in adaptation in the two continents. Signalling was the most enriched cluster in the Europe Fall data set, while development and morphogenesis were enriched both in the Europe Summer and Europe Fall data sets (Table 3). These results are similar to previous analysis performed in D. melanogaster: development and morphogenesis have been reported in population differentiation studies in different continents such as Europe, Australia, and North America (Fabian et al., 2012; Kolaczkowski et al., 2011; Mateo et al., 2018; Reinhardt et al., 2014). Note that excluding the SNPs that are located inside inversions led to very similar enriched biological processes GO terms (Table S8B).

TABLE 3

GO enrichment analysis of candidate genes for local adaptation

Data set	Significant SNPs	Significant genes	GO enrichment terms
Europe	719	410	Neuron development; eye development; signalling; organ morphogenesis; growth
North America	1,164	583	Response to stimulus; organ development; regulation of growth; nervous system development; localization and transport
Europe Summer	752	396	Learning/memory; eye development; neuron development; sensory perception of pain; organ morphogenesis
Europe Fall	821	412	Signalling; localization/transport; organ morphogenesis; neuron development; heart morphogenesis

For each data set, the number of genes and significant single nucleotide polymorphisms (SNPs), located in the gene body and upstream region (<1 kb), and the top five most enriched GO terms (significance >1.3). The significance of the SNPs was determined based on the empirical distribution of the calibrated XtX* values (top 0.05%), which corresponds to q‐value thresholds of 7.56e‐10 in Europe, 3.70e‐07 in Europe Summer, 1.10e‐05 in Europe Fall, 9.90e‐06 in North America for autosomes; and to q‐value thresholds of 1.49e‐05 in Europe, 0.0003 in Europe Summer, 0.000441 in Europe Fall and 0.03 North America for X chromosome.

Finally, we also tested whether there was a significant overlap between the candidate genes for local adaptation found in Europe and in North America. We found 55 significant genes overlapping in the two continents (SuperExactTest p‐value <0.05; Table S9A, see Materials and Methods). Among these 55 genes, we found already known genes such as cpo and Ace, as well as other genes previously identified in North American clinal studies such as Cow, involved in neuromuscular junction development (Kopke et al., 2020) or dpy, involved in wing and trachea development (Wilkin et al., 2000) (Table S7; Table S9B). We performed a GO enrichment analysis with these 55 overlapping genes among continents, and the main clusters were related to regulation, signalling and response to stimulus, and development (Table S9C).

Temperature, rainfall, and wind are the most contributing variables in the genome‐environment association analyses

To identify the environmental variables that are relevant for adaptation in D. melanogaster natural populations, we looked for significant associations between SNPs frequencies and several environmental variables using the BayPass standard model (see Materials and Methods). We analysed 59 environmental variables related to temperature, rainfall, evaporation, solar radiation, wind, daylight hours, and soil type (Table S3 and S5, see also Materials and Methods), and we found significant associations between at least one of these variables and 748 genes (Table S10).

For all data sets, temperature was the variable associated with the highest number of genes, followed by wind in Europe and Europe Fall, and rainfall in the North America and Europe Summer data sets (Table 4). Significant SNPs located in some of these genes were associated with more than one variable as expected from the correlation found between some of the covariates (Figure S1). For instance, in North America most of the SNPs associated with solar radiation variables were also associated with Temperature variables (84/103), which is consistent with the correlation found between these variables (Figure 3b). Note that, the majority of SNPs that were associated with wind were not associated with any other environmental variables (Figure 3a,d), which is consistent with the lack of significant correlation between wind and other environmental variables (Figure S1). On the other hand, in Europe Summer the majority of SNPs associated with evaporation were also associated with temperature, although we did not find a strong correlation between evaporation and temperature variables (Figure 3c, Figure S1D). However, there are studies reporting similar responses to cold and desiccation stress (Sinclair et al., 2007). Among the top five genes with the highest BF scores in the four data sets, we found several associations with Annual mean temperature and Annual mean solar radiation (Table 5).

TABLE 4

Summary of results obtained for environmental association

	Europe	NA	Europe Summer	Europe Fall
Temperature	143	217	87	33
Wind	118	83	17	20
Rainfall	29	152	62	13
Evaporation	36	79	51	7
Solar radiation	45	103	19	3
Soil	–	4	–	–
Daylight hours	18	52	8	1
Total	296	382	155	64

Number of genes with significant SNPs (BF ≥ 30) located in the gene body or upstream region (< 1 kb) for each type of environmental variable. In bold, three top type of environmental variables with more genes for each data set.

FIGURE 3

Overlap between genes with single nucleotide polymorphisms (SNPs) significantly associated with environmental variables. For each data set, genes with SNPs in the gene body and upstream region (<1 kb) significantly associated with the three groups of environmental variables with more genes associated with them are depicted [Colour figure can be viewed at wileyonlinelibrary.com]

TABLE 5

Candidate genes associated with environmental variables

Gene Name	SNP location	Data set	Strongest association variable	BF	Phenotype
Ace	Gene body	Eu	Annual mean temperature	84.89	Diapause, insecticide resistance, olfactory, starvation
Sap‐r	Gene body	Eu	Annual mean solar radiation	62.74	Starvation
obst‐F	Gene body	Eu	Annual mean temperature	72.05	‐
apn	Gene body	Eu	Mean temperature of warmest quarter	60.59	‐
Ptr	Gene body	Eu	Mean evaporation of warmest quarter	70.75	‐
tok	Gene body	NA	Annual mean solar radiation/Solar rad mean diurnal range	64.34	Circadian, starvation
Mhc	Gene body	NA	Solar rad mean diurnal range	61.76	‐
CG13705	Gene body	NA	Temperature seasonality	58.18	‐
Abd‐B	Gene body	NA	Annual mean solar radiation	51.11	Alcohol, dessication, pigmentation
Gale	Upstream	NA	Wind mean diurnal range	50.57	Aggressiveness, diapause, immunity, starvation
Ace	Gene body	EuS	Annual mean temperature	74.38	Diapause, olfactory, starvation
obst‐F	Gene body	EuS	Max temperature of warmest month	68.08	‐
CG7290/ CG7298	Gene body/Upstream	EuS	Max temperature of warmest month	58.37	‐ /Hypoxia, immunity, xenobiotic
CG10257	Upstream	EuS	Precipitation of driest quarter	58.48	Xenobiotic
Ptr	Gene body	EuS	Mean evaporation of warmest quarter	62.11	‐
capu	Gene body	EuF	Wind variability index	49.59	Alcohol, circadian behaviour, oxidative, xenobiotic
Kek2	Gene body	EuF	Wind variability index	44.49	‐
Argk	Upstream	EuF	Temperature seasonality	65.35	Immunity, starvation
RFeSP	Gene body	EuF	Wind seasonality	48.51	Hypoxia
CG43750	Gene body	EuF	Wind seasonality	42,35	‐

For each data set, the top five genes with significant single nucleotide polymorphisms (SNPs) located in the gene body and upstream region (<1 kb) with the highest significant Bayes factor (BF) scores, the environmental variable with the strongest association and their associated phenotype (see Table S12). All significant genes can be found in Table S10.

Abbreviations: EuA, Europe; EuF, Europe Fall; EuS, Europe Summer; NA, North America.

We also tested whether candidate genes with SNPs associated with environmental variables were enriched inside cosmopolitan inversions (Figure 2; Corbett‐Detig & Hartl, 2012). We found an enrichment of significant SNPs in the In(2L)t inversion in the Europe Fall data set and in the In(3R)P inversion in the North America data set (Fisher's exact test p‐value <0.001; Table S6B). The In(2L)t inversion in the Europe Fall data set was enriched for SNPs associated with temperature variables while the In(3R)P in the North America data set was enriched for SNPs associated with rainfall, solar radiation, and wind variables (Table S6C–D). In previous studies In(2 l)t and In(3R)P were correlated with climatic factors varying latitudinally in North America, specifically with temperature and rainfall (Kapun et al., 2016). Our analysis suggests that other climatic factors such as wind may also be correlated with inversions.

Finally, we also found a significant overlap between the genes with SNPs significantly associated with environmental variables identified in the North American and the European data sets (SuperExactTest p‐value <0.05; see Materials and Methods; Table S9 and S10). Among the 32 significantly overlapping genes, fipi is involved in the Drosophila courtship song (Fedotov et al., 2018) and was associated with variables related to wind in Europe and North America. Courtship song, as well as wind have been shown to activate neurons which are related to the antennal and mechanosensory motor center in the central brain in D. melanogaster (Matsuo & Kamikouchi, 2013; Yorozu et al., 2009).

23% to 51% of the genes significantly associated with environmental variables did not show adaptive differentiation across populations

We found that, across data sets, only 12% to 37% of the genes that showed patterns of population differentiation (XtX*) were associated with at least one environmental variable (Table S7). Indeed, most of the genes that showed the highest association with environmental variables, such as Ace and obst‐F, were also among the top candidates for significant population differentiation (Tables 2 and 5). Another example is Gale, which in North America was associated with a wind variable, and has been related with aggressiveness and diapause responses as well as with immunity and starvation stresses (Tables 2 and 5; Clark et al., 2013; Edwards et al., 2006; Fukuyama et al., 2013; Grönke et al., 2005; Harbison et al., 2005; Kučerová et al., 2016; Shorter et al., 2015; Zhao et al., 2016).

On the other hand, we found that 23% to 51% of the genes that showed a significant association with at least one environmental variable, did not show population differentiation patterns (Table S10). Among these genes, RFeSP has one of the highest association scores with Wind seasonality in the Europe Fall data set (Table 5 and Table S10). This gene encodes Rieske iron sulphur proteins which are highly conserved functional constituents of energy‐transducing respiratory complexes (Gontijo et al., 2011).

Ten transposable elements insertions are associated with environmental variables

In addition to SNPs, we also analysed the population differentiation patterns and correlations with environmental variables for TE insertions (Table S11A). We found that nine out of the 403 TE insertions showed patterns of population differentiation (XtX*) in at least one of the data sets analysed; however, we did not find overlap between continents (Figure 2, Table 6, Table S11A–B). Four of these TEs have previously been identified as candidate adaptive TEs (Table 6).

TABLE 6

Significant transposable element insertions found in the population differentiation analysis

Transposable element	Family	Location	Gene	Data set	Evidence of selection
FBti0019112	pogo	First intron	lilli	Eu, EuF	iHS, H12, nSL (Rech et al., 2019)
FBti0019164	X‐element	First intron	CG9932	Eu	Population differentiation (González et al., 2008)
FBti0019144	Rt1b	First intron	CG44153	EuF	Population differentiation (González et al., 2008)
FBti0019276	S‐element	Second intron	Adf1	EuF	CSTV (Lerat et al., 2019)
FBti0019862	G6	432 bp downstream	Tif‐IA	Eu, EuS	–
FBti0020056	BS	507 bp downstream	bin	NA	–
FBti0020306	hopper	Third intron/first intron	atms/CG44098	NA	–
FBti0060187	G2	First intron	Syn1	NA	–
FBti0061428	hobo	52 bp upstream/529 bp downstream	CG31809/CG6012	NA	–

TE insertions were considered significant if their associated XtX* values were above the top 1% of the empirical distribution of XtX* values, and q‐value <0.05. When the TE insertion is located in intergenic regions, genes located nearby are reported (Table S11A).

Abbreviations: CSTV, correlation with spatiotemporal variables; Eu, Europe; EuF, Europe Fall; EuS, Europe Summer; NA, North America.

In addition, we identified six TE insertions in the Europe data set that showed significant associations with different environmental variables: four of them showing the highest association with temperature variables, one with evaporation, and one with rainfall (Table 7 and Table S11C). Three of these insertions also showed significant patterns of population differentiation (FBti0019112, FBti0019164 and FBti0019862; Table 6). FBti0019112 showed the highest BF value and was associated with the Minimum temperature of the coldest month variable (Table 7). This insertion is located in an intron of the lilli gene, which is mainly involved in cell identity and growth, and plays a role in retinal development (Distefano et al., 2012; Wittwer et al., 2001). This result is interesting given that other studies in Drosophila showed the impact of temperature in eye development genes (Del Bel et al., 2018). In addition, llili has been suggested to have a role in local adaptation, as it was recently reported to be part of a strong outlier region in a study comparing D. melanogaster collected in wilderness areas and collected in the nearby of towns in southern‐central Africa (Sprengelmeyer et al., 2020). FBti0018880 showed the second strongest association also with temperature, Isothermality (temperature variability index), and has been reported to play a role in oxidative stress response (Guio et al., 2014). Other studies performed in D. melanogaster showed that metabolites involved in oxidative stress are altered by selection in cold tolerance (Koštál et al., 2016; Williams et al., 2014).

TABLE 7

Significant candidate TE insertions associated with environmental variables (BF ≥ 20)

Transposable element	Environmental variable	Significant XtX*	BF	Data set
FBti0018880	Isothermality	No	30.53	Eu
FBti0019112	Min temperature of coldest month	Yes	43.38	Eu
FBti0019164	Temperature Annual range	Yes	24.79	Eu
FBti0019165	Evaporation Mean diurnal range	No	20.52	Eu
FBti0020057	Precipitation Seasonality	No	22.12	Eu
FBti0019862	Isothermality	Yes	23.69	Eu
FBti0061428	Annual mean wind	Yes	43.95	NA
FBti0020086	Solar radiation variability index	No	26.83	NA
FBti0020306	Precipitation of wettest quarter	Yes	26.02	NA
FBti0019318	Mean evaporation of coldest quarter	No	28.57	NA

The environmental variable with highest score is reported (Table S11).

Abbreviations: Eu, Europe; NA, North America.

We also identified four TE insertions in the North America data that showed significant associations with wind, solar radiation, rainfall and evaporation (Table 7; Table S11A,C). Two of these insertions also showed patterns of population differentiation (FBti0061428 and FBti0020306; Table 6).

DISCUSSION

In this work, we aimed to identify the main environmental drivers of adaptation of D. melanogaster natural populations in a large continental geographical scale. To accomplish this, we used GEA analysis on a large set of population samples (up to 20 in Europe and 11 in North America) representative of different environments and considering a wide‐range of environmental covariates to capture this variability. We found that in addition to temperature and rainfall, wind related variables appear to be also relevant for D. melanogaster environmental adaptation. Temperature and rainfall are the most widely used variables in GEA analysis in several species including D. melanogaster (Božičević et al., 2016; Cavedon et al., 2019; Gao et al., 2017; Hopley & Byrne, 2019; Kapun et al., 2020; Leroy et al., 2020; Mayol et al., 2020; Pina‐Martins et al., 2019; Todesco et al., 2019). Our results show that the majority of genes associated with environmental variables were associated with a temperature‐related one (400/748), while the number of genes associated with rainfall was smaller (241/748) (Table S10). These results are consistent with similar GEA analysis performed previously in D. melanogaster (Kapun et al., 2020). Moreover, among the 748 candidate genes associated with environmental variables in our study, 226 were associated with a wind‐related variable being the third variable group with most associations, and far from the following ones (evaporation and solar radiation with 153/748 genes each) (Table S10). Wind‐related variables have been studied mainly in plant adaptation (Balkenhol et al., 2017), and are often included as part of PCs where individual wind effect cannot be properly measured (Gao et al., 2017). In other species, wind has also been reported to be involved in desiccation stress and thermoregulation (Baig & Tranquillini, 1980; Ortega et al., 2017). In Drosophila species, including in D. melanogaster, it has already been suggested that wind might be relevant for adaptation (Fuller et al., 2014; Patella & Wilson, 2018). The effect of wind variables in Drosophila could be related to the Johnston's organ, which is the largest mechanosensory organ in Drosophila. This organ is involved in a variety of behaviours such as hearing, touch, vestibular sensing, proprioception and wind sensing (Patella & Wilson, 2018). In addition, Fuller et al. (2014) reported how flies regulate flight speed according to the information from their visual system and their antennae, and how they can overcome the effect of sudden wind disturbances. To the best of our knowledge, our analysis is the first that identifies genome‐wide variants associated with wind‐related variables in D. melanogaster. Although temperature, rainfall, and wind seem to be important drivers of adaptation, we still lack information about other variables directly related to them and that may be actually underlying adaptive processes. Further analysis testing the direct effect of these three variables on the genetic variation should be performed to obtain conclusive evidence.

We found that an important proportion of the genes showing signals of adaptive differentiation did not show associations with any of the environmental variables studied (>60%). As we addressed previously, it is difficult to know a priori the variables that may be relevant for adaptation, so, it could be that we are not including in our analysis the environmental variables responsible for the adaptive processes in which these genes are involved. For example, 185 of the 1,300 genes showing population differentiation patterns are candidates for xenobiotic stress response (Rech et al., 2019). For the majority of these genes, we did not find any association with environmental variables. Thus, including variables related to pollution might help explain the population differentiation in some of these genes. It could also be that although the relevant environmental variables are included in the analysis, our samples do not allow us to capture the whole range of the variation of these environmental variables making the GEA analyses less powerful. Alternatively, population differentiation patterns in some of these genes might be due to selective pressures not related to the environment. We also found that between 23% and 51% of the candidate genes showed association with an environmental variable but did not show significant population differentiation. Thus, GEA analyses not only identifies the relevant environmental variables, but also allows to identify genetic variants involved in environmental adaptation that cannot be detected through population differentiation analysis, as expected if they result in too subtle changes in allele frequencies across populations (Gautier, 2015).

Our work also aimed at evaluating to which extent our observed signal of adaptation were consistent across the European and North American continents. We found that 55 genes showing patterns of population differentiation, and 32 genes showing association with at least one environmental variable, overlapped in these two continents. Out of these 32 genes, 14 were associated with a different environmental variable group in each continent, and 12 were associated with different variables in the same group. This suggests that although there is a pattern of parallel adaptation, there may be different environmental pressures which may drive adaptation for the same genes.

We also assessed whether samples collected in European populations in summer and fall differed in their association with environmental variables. Recent studies have shown the role of temperature in seasonal variation (Machado et al., 2019). We found that in the summer and in the fall data sets the majority of genes were associated with temperature and rainfall (Figure 3c, d, Table 4). However, there was a substantial proportion of candidate genes associated with evaporation in Europe Summer but not in Europe Fall (33% vs. 11%; Table 4). On the other hand, there were more candidate genes associated with wind variables in the Europe Fall than in the Europe Summer data set (31% vs. 11%; Table 4). These results suggest that different environmental variables, evaporation and wind, might play a role across seasons. However, temporal data series from several years should be analysed to confirm these results.

Finally, we identified four TE insertions showing significant population differentiation patterns, five TE insertions associated with an environmental variable, and five insertions showing both. We described as candidates for the first time three of these TE insertions, FBti0019862, FBti0020306 and FBti0061428, which are associated with environmental variables and showed significant population differentiation patterns, as well as FBti0019164 only reported in González et al. (2008) and FBti0019112 which has shown previous evidence of selection (Rech et al., 2019). This analysis is, however, limited as we could only investigate those TE insertions present in the reference genome and that were polymorphic in our samples. We suggest that both reference and nonreference insertions should be included in future analysis in order to get a comprehensive picture of the role of TE insertions in environmental adaptation.

Overall, we identified temperature, rainfall and wind as environmental variables which may play a critical role in environmental adaptive processes in D. melanogaster. In addition to increasing the number of populations and of TE insertions analysed, we further suggest that performing GEA analysis in populations collected across time should inform us about how the role of environmental variables changes through time and contributes to the dynamics of genetic variation across populations and to the maintenance of adaptive variants. Extending this analysis to other continents should also further enhance our understanding of the role of environmental variables in adaptive evolution.

AUTHOR CONTRIBUTIONS

Designed research: Josefa González. Performed research: María Bogaerts‐Márquez, Sara Guirao‐Rico and Mathieu Gautier. Contributed new reagents or analytical tools: María Bogaerts‐Márquez. Analysed data: María Bogaerts‐Márquez and Josefa González. Drafted the manuscript: María Bogaerts‐Márquez and Josefa González. Edited and approved the manuscript: María Bogaerts‐Márquez, Sara Guirao‐Rico, Mathieu Gautier and Josefa González.

Fig S1

Click here for additional data file.

Table S1‐S12

Click here for additional data file.