Genome-wide patterns of local adaptation in Western European Drosophila melanogaster natural populations

Signatures of spatially varying selection have been investigated both at the genomic and transcriptomic level in several organisms. In Drosophila melanogaster, the majority of these studies have analyzed North American and Australian populations, leading to the identification of several loci and traits under selection. However, several studies based mainly in North American populations showed evidence of admixture that likely contributed to the observed population differentiation patterns. Thus, disentangling demography from selection might be challenging when analyzing these populations. European populations could help identify loci under spatially varying selection provided that no recent admixture from African populations would have occurred. In this work, we individually sequence the genome of 42 European strains collected in populations from contrasting environments: Stockholm (Sweden) and Castellana Grotte (Southern Italy). We found low levels of population structure and no evidence of recent African admixture in these two populations. We thus look for patterns of spatially varying selection affecting individual genes and gene sets. Besides single nucleotide polymorphisms, we also investigated the role of transposable elements in local adaptation. We concluded that European populations are a good dataset to identify candidate loci under spatially varying selection. The analysis of the two populations sequenced in this work in the context of all the available D. melanogaster data allowed us to pinpoint genes and biological processes likely to be relevant for local adaptation. Identifying and analyzing populations with low levels of population structure and admixture should help to disentangle selective from non-selective forces underlying patterns of population differentiation in other species as well.

Because the effective population size of the X chromosome is ¾ of that of the autosomes, we expect reduced levels of polymorphisms on the X chromosome compared to the autosomes 47 . Indeed, we found that the average polymorphism on the X chromosome is reduced relative to autosomes both in the Swedish and Italian populations ( Fig. 2 and Table S3, Supporting Information). On the other hand, we did not find reduced polymorphism in the X vs the autosomes in Rwanda (Fig. 2 and Table S3, Supporting Information). This has previously been reported, and it was initially explained in terms of unequal sex ratio in African populations, where female population size was estimated to be 1.8 times larger than male population size 48 . However, re-analysis of this data suggested that differences in recombination rate between the X and the autosomes could not be excluded as an alternative explanation for the higher polymorphism on the X chromosome 49 . According to the African origin of D. melanogaster, diversity in non-African populations is expected to be a subset of that within Africa 50 . Consistently, we found that the Rwanda population showed the highest genomic diversity of the four populations analyzed (Fig. 2).
Global F ST values showed that the European populations are the less differentiated (Fig. 3). F ST between African and European populations is higher than between African and the North American population consistent with previous analysis reporting a secondary contact between African and North American strains [25][26][27] (Fig. 3).
Principal component analyses clearly separated strains by continent (Fig. 4). The two first principal components (PCs) explained a 39.3% of the observed variance. PC1, which explains 36.3% of the observed variance, separates African strains from non-African strains, and to a lesser extent American and European strains (Fig. 4). PC2, which captures 3% of the variance, separates American from European strains. Consistent with F ST estimates, European strains appeared closer to each other than American strains, and African strains (Fig. 4).
Finally, we also performed model-based global ancestry estimation (Fig. 5). The most likely hypothesis according to cross-validation errors, is a model with two or three ancestral populations ( Figure S3, Supporting information). The model with two ancestral populations (k = 2) clearly shows a cluster of African ancestry and a cluster of non-African ancestry, with the American population showing an average of 32% of African ancestry. The model with three ancestral populations (k = 3) separates strains by continent (Fig. 5). Increasing the number of ancestral populations does not separate the two European populations, suggesting that Swedish and Italian strains belong to a single cluster (Fig. 5).  Table S1, Supporting Information). Overall, as expected, our results are consistent with an out-of-Africa migration, as Rwandan samples are the most genetically diverse (Fig. 2). Our results are also consistent with population structure among continents (Figs 3 and 4), and with North American strains showing a considerable proportion of African admixture (Fig. 5). Finally, the two European population analyzed showed low levels of population structure (Figs 3 and 4) and no recent African admixture (Fig. 5) suggesting that they could be a good dataset to identify loci under spatially varying selection.

Genomic differentiation between Swedish and Italian populations.
To identify signatures of population differentiation between Swedish and Italian strains, we estimated pairwise F ST for all individual SNPs segregating at a minor allele frequency higher than 5% in these two populations (see Material and Methods). The empirical genome-wide distribution of F ST , which reflects primarily drift and other non-selective forces, had a maximum around zero and exponentially decays with a long tail of high F ST values ( Figure S4A, Supporting information). We found that the levels of differentiation were significantly heterogeneous among chromosomal arms (ANOVA p-value < 2e-16; Figure S4B, Supporting information). Therefore, we used the empirical F ST distribution of each chromosomal arm to identify candidate differentiated SNPs in each chromosome. SNPs falling in the top 5%, top 1% or top 0.05% tails of the distribution were considered as candidates ( Figure S4C, Supporting information).
To investigate the contribution of inversions to the observed population differentiation patterns, we first estimated the frequency of cosmopolitan inversions in Swedish and Italian populations (Table S4, Supporting Information). The only inversion present at a significant frequency was In(2 L)t: 18.8% in the Italian strains analysed and 31.5% in the Swedish strains analysed. In(2 R)NS was present at low frequencies in both populations, while In(3 R)P and In(3 L)P were either absent or present at low frequency (Table S4, Supporting Information). We thus only considered In(2 L)t for the rest of this work.  We then tested whether F ST values inside In(2 L)t common cosmopolitan inversion was different from those of the rest of the chromosomal arms. Additionally, we tested whether the proportion of candidate SNPs (top 5%, top 1% and top 0.05%) differed between the inversion and the rest of the chromosomal arm (Fig. 6A). The median F ST inside In(2 L)t inversion was significantly higher than the median F ST outside this inversion: 0.0080 vs 0.0068 (p-value = 0.0113). However, we did not find a significant accumulation of candidate differentiated loci inside the inversion (Fig. 6A).
We then analyzed the functional class of the differentiated SNPs. We first compared the distribution of non-differentiated and differentiated SNPs across the different functional categories and found that, as expected, small introns are enriched for non-differentiated SNPs and are depleted of differentiated SNPs (Fig. 6B). Our results also showed that core promoter, distal promoter, and intergenic regions are enriched for differentiated SNPs (Fig. 6B). Considering small introns as the background and controlling for allele frequency, chromosome, and In(2 L)t inversion status, we found that although other functional categories are enriched in differentiated SNPs, distal promoter, core promoter and intergenic regions are the three categories that showed a more significant enrichment of candidate SNPs (Table S5, Supporting Information).

Biological processes involved in population differentiation.
To identify the biological processes likely to be involved in population differentiation, we first converted the SNP-centric F ST values into gene-centric Z ST scores while controlling for gene length bias (see Material and Methods and Daub et al. 45 ). Z ST scores are approximately normally distributed with positive values indicating high levels of population differentiation and negative values indicating low levels of population differentiation ( Figure S1, Supporting information). We obtained a Z ST score for 13,140 genes, covering 84% of all annotated D. melanogaster genes. We considered as candidate-differentiated genes the 657 genes with a Z ST score in the top 5% of the distribution (Table S6, Supporting Information).
We then performed GO enrichment analysis using an algorithm based on the classical Fisher's exact test (classic), and three additional more strict algorithms (elim, weight and weight01) that eliminate local similarities and dependencies between GO terms in order to reduce the false-positive rate (see Material and Methods). The candidate 657 genes were enriched in 43 GO terms according to at least one of the more strict algorithms, which we classified into broader categories to simplify their interpretation (Table S7, Supporting Information). The significant biological processes were related to response to xenobiotics, pigmentation, immunity, and developmental processes among others (Table S7, Supporting Information). We also checked whether there was overlap between the 43 GO terms identified in this and previous works, and we found that 6 of the 43 GO terms have been previously identified (Table S8, Supporting Information). Among the six overlapping GO terms, we found GO:0030707 and GO:0007297, which are related to ovarian follicle cell development and ovarian follicle cell migration, respectively. Interestingly, these processes are associated with diapause in D. melanogaster 51,52 , a complex phenothypic trait involved in the adaptation to temperate regions through overwintering 53,54 . Candidate genes and candidate SNPs involved in population differentiation. To identify the most likely candidate SNPs to play a role in local adaptation, we focused on the 22 genes that belong to the four GO terms that were consistently found to be overrepresented across algorithms: response to insecticide, cuticle pigmentation, neuropeptide signalling pathway, and cell wall macromolecule catabolic process (Table 1). For each gene in these four significant GO terms, we identified the SNP with the highest F ST (Table 1). We checked whether the candidate genes were located inside or nearby the breakpoints of the cosmopolitan inversion In(2 L)t. Only one of the 22 candidate genes, Duox, was located inside In(2 L)t inversion. However, this gene has been previously identified when analysing different natural populations, and also showed expression changes when comparing temperate with tropical populations 28,29,55 .
The strongest F ST signal in the genome corresponds to a coding change in the Acetylcholinesterase (Ace) gene that has been previously identified as conferring insecticide resistance in D. melanogaster 56 . Ace is a well-known insecticide resistance gene across insects 56,57 . In this gene, we also identified two additional coding mutations that are known to confer insecticide resistance (3 R: 9069408_C/G, Gly303Ala, and 3 R: 9063921_C/G, Gly406Ala with F ST values of 0.05 and 0.65, respectively) 56,58 . Besides these three coding mutations, we identified two new coding mutations with unknown functional impact (3 R: 9071797_T/C, Thr20Ala; 3 R: 9071847_A/G, Ile8Thr  (Table 1). This gene has also been repeatedly reported as a candidate gene in population differentiation studies across continents and it has been experimentally shown to confer resistance to insecticides 59-61 . Daborn et al. 60 identified a TE insertion upstream of Cyp6g1 that increased the expression of this gene resulting in resistance to DDT and neonicotinoids 59,60 . More recently, Battlay et al. 61 identified several SNPs in the Cyp6g1 gene associated with resistance to a different type of insecticide, the organophosphate azinphos-methyl, with the strongest association being from a SNP located in an intron. The mutation identified in this work is a coding mutation, and thus could represent an additional Cyp6g1 genetic variant involved in insecticide resistance (Table 1). Finally, we found that two other insecticide resistance related genes Cyp6a8 and Gr8a also have coding changes while Cyp6g2 has a mutation in the core promoter (Table 1). Regarding cuticle pigmentation and neuropeptide signalling, all changes were non-coding changes or synonymous changes, while cell wall macromolecule catabolic process contained non-coding or non-synonymous changes. This later GO term contains the gene with the second highest F ST in the genome: a SNP in the core promoter of LysC (Table 1). It is noteworthy that both in cuticle pigmentation (yellow-e3, yellow-e) and in cell wall macromolecule catabolic process (LysC, LysD, LysE, LysP, LysS) the differentiated genes are located next to each other, forming two physical clusters of genes that are co-differentiated and functionally related ( Table 1). The two clusters are located in regions with high local recombination rates, 2.69 cM/Mbp and 3.15 cM/Mbp respectively, and expand several kilobases, 11.8 kb and 17.8 kb respectively. Thus, although we cannot discard that some of the identified candidate SNPs are linked to the causative adaptive mutation, it could also be that several of them are adaptive (Table 1). Indeed, for LysP, LysE, and yellow-e there are independent lines of evidence suggesting that they are under selection 28,29,62,63 .
We also identified the genes in our dataset that have been reported as candidates in previous studies (Table S1, Supporting Information). The number of overlapping genes ranged from zero to 177: while no overlap was found between our candidates and the ones identified by Juneja et al. 63 when comparing two Australian populations, 177 out of 657 genes significantly overlapped between our study and the ones identified by Machado et al. 28 in the North American East coast (Fisher's exact test, p-value = 0.01904; Table S1, Supporting Information). Overall, 363 of the 657 genes (55%) identified in this study have already been identified in at least one previous study. Interestingly, 35 of these 363 genes, besides showing evidence of population differentiation in other continents, also showed differential gene expression when comparing tropical and temperate populations (Table S9, Supporting Information). These 35 genes, which are involved in several ecologically relevant biological processes, are thus good candidates for follow-up functional validation (Table S9, Supporting Information).
Gene set enrichment analysis. Classical GO enrichment analysis is powerful when several genes belonging to the same GO term show strong signals of differentiation that go beyond a threshold (top 5% in this work). However, it fails to detect more subtle trends in which most genes annotated to the same GO term might not 95 belong to the top 5% of the Z ST distribution but as a whole tend to have higher Z ST scores than expected by chance 64 . In order to detect those GO terms whose genes tend to occur toward the top of the Z ST distribution, we performed a Kolmogorov-Smirnov like test. A total of 31 GO terms were found to be significantly enriched by the three algorithms and 58 additional GO terms were found to be significantly enriched by at least one of the more strict algorithms. The list of 89 GO terms was classified into broader categories to simplify its interpretation (Table S10, Supporting Information). While some broader categories, such as developmental processes and genome organization, overlapped with the GO terms enriched in the top 5% differentiated genes, others such as circadian rhythm were only identified when performing gene set enrichment analysis (Table S10, Supporting Information). We found that 13 of the 89 GO terms were previously identified by other works (Table S11, Supporting Information).
Overall, the gene set enrichment analyses allowed us to identify additional biological processes potentially involved in differentiation between the two European populations analyzed.
Transposable elements also contribute to European population differentiation. Besides SNPs, we also investigated whether reference TE insertions, that is TEs annotated in the reference genome, showed patterns of differentiation between the Swedish and Italian populations. We first genotyped and estimated the frequency of TEs in these two European populations using the T-lex2 pipeline 65 . We then calculated F ST values for 161 TEs segregating at a minor allele frequency higher than 5% in the two populations (see Material and Methods). The empirical distribution of the F ST estimates for TEs was very similar to the SNPs F ST values distribution ( Figure S5, Supporting Information). Therefore, we used the SNPs critical F ST values of each chromosomal arm to determine the significantly differentiated TEs. We identified six TEs with significant population differentiation: four of them are present at higher frequencies in the Swedish population while the other two are present at higher frequencies in the Italian population (Table 2).
We identified the genes located nearby these six TEs (<1 kb), and analyzed their associated biological processes GO terms (Table 2). Interestingly, one of the genes nearby insertion FBti0020162, which is present at higher frequency in Swedish compared with Italian populations, is involved in cellular response to cold ( Table 2). Other ecologically relevant biological processes such as detection of chemical stimulus and response to starvation could be affected by FBti0019679 and FBti0019344 insertions respectively ( Table 2). Five of the eight genes nearby TEs with significant F ST between Swedish and Italian populations were previously identified in other works ( Table 2).

Discussion
In this work, we tested whether Western European populations could be a good dataset to disentangle the effects of demography and selection, and thus to identify the genetic basis of adaptation. These two populations were collected in locations with contrasting environments: Stockholm in Sweden, and Castellana Grotte, in Bari, Southern Italy (Fig. 1). We first explored the genetic diversity, population structure, and admixture patterns of these two populations in the context of a more global dataset including an African population from the ancestral range of the species, and a North American population 43 We cannot discard that the European populations analyzed in this work have African admixture from a different source African population. Kao et al. 27 reported some proportion of admixture in strains collected in France when a population from Cameroon was used as a source African population. However, Bergland et al. 25 analyzed admixture in North American and Australian populations using 19 different African populations, including the Gikongoro population used in this study. Admixture was detected independently of the donor population used in the analysis. If there was some degree of admixture that we failed to detect, the patterns of population differentiation identified for some of the European genetic variants could be the result of admixture rather than selection [25][26][27]67 . To tease apart selective from non-selective signatures of population differentiation, we analyzed our results in the context of all the previous available genome-wide datasets including both genomic and transcriptomic data. Our analyses allowed us to identify new genes likely to be involved in local adaptation, and to provide further evidence for previously identified candidate genes (Tables 1 and 2, and Table S6, Supporting Information). Thirty-five of the 657 genes identified in this study, besides showing evidence of population differentiation in other continents, also showed evidence of differential gene expression in temperate vs tropical populations (Table S10, Supporting Information). We argue that these genes, several of which are involved in ecologically relevant processes such as insecticide resistance and pigmentation, are good candidates to underlie local adaptation (Table S10, Supporting Information).
For some of the genes already known to play a role in local adaptation, we pinpointed new candidate genetic variants suggesting that different variants might be under selection in different populations. Adding European populations to the already available studies is thus enhancing our understanding of adaptive evolution. This is the case of Ace a well-known insecticide resistance gene for which, besides providing further evidence for three coding SNPs previously shown to confer resistance to insecticides, we identified two additional coding SNPs that could have also been targets of selection 56,58 . Another example is a coding SNP in the Cyp6g1 gene, which is also involved in insecticide resistance 60,61 . Besides Ace and Cyp6g1, two other genes functionally annotated as "response to insecticide" were present at higher frequencies in the Italian population. (Table 1). This is consistent with Italian and Swedish populations being subject to different insecticide pressures, which is in agreement with the characteristics of these two populations. While the Italian populations were collected in a vineyard in which pesticides were regularly used, the Swedish populations were collected in urban community gardens and thus are less likely to be subject to strong pesticide use.
Besides response to insecticide, other GO terms such as cuticle pigmentation, defense response, and cell wall macromolecular catabolic process previously related to environmental adaptation were also identified (Table S7, Supporting Information). Changes in four of the five genes involved in pigmentation were present at higher frequencies in the Italian population (Table 1). Variation in pigmentation in D. melanogaster natural populations has been related to different adaptive phenotypes: adaptation to cold environments, desiccation tolerance, and ultra-violet (UV) resistance (reviewed by True 68 ). The Swedish and Italian populations analyzed in this work differed in temperature, humidity, and UV exposure, and thus it is reasonable to find significant differentiation in SNPs related to pigmentation 69 . Defense response is also considered as an important trait affecting the ability of species to colonize new environments 70 . Besides the genes functionally annotated as defense response, three of the six genes functionally annotated as cell wall macromolecule catabolic process have been associated with anti-parasitic defense response 71 . Finally, the GO analysis also pinpointed six GO terms previously identified in clinal variation studies (   51,52 , an ecologically relevant trait putatively involved in the adaptation of D. melanogaster to temperate climates 53,54 . Besides analyzing individual genes likely to be involved in geographical adaptation, we also performed gene set enrichment analyses. This analysis complements the identification of mutations with strong effects by identifying biological pathways for which several small effect mutations have occurred 45 . We indeed found evidence for polygenic adaptation in several biological processes (Tables S10 and S11, Supporting Information). Some of them, such as circadian rhythm, were not identified with the outlier approach while others, such as detection of chemical stimulus, seem to be affected both by large effect and small effect mutations. These results suggest that in addition to identifying individual genes, gene set enrichment analyses should be performed to get a more complete picture of adaptation.
Finally, while previous studies focused on a single type of genetic variant, besides SNPs we also analyzed reference TE polymorphisms (Table 2). A total of 4% of the analyzed TEs showed significant patterns of population differentiation. Some of these TEs are located nearby genes involved in ecologically relevant biological processes such as cold adaptation or detection of chemical stimulus ( Table 2). Although we analyzed a much bigger dataset compared to previous studies (1,632 TEs vs 902 TEs), the number of TEs for which F ST could be estimated and as such the number of candidate differentiated TEs was relatively small [36][37][38] . However, our dataset is still incomplete as we focused on the TEs annotated in the reference genome, which is a North American strain 37 . We are thus missing all the TEs present in the two European populations that are not shared with the North American strain sequenced. De novo annotation of TEs in the two European populations analyzed in this work would probably result in the identification of more candidate TEs. Thus, further analyses are needed to better assess the role of TEs in local adaptation.
Overall, our results suggest that European populations are a good dataset to advance our knowledge on the genetic basis of environmental adaptation. By analyzing populations with low levels of population structure and no evidence of recent admixture, we were able to confirm and to de novo identify genomic targets of spatially variant selection, and to pinpoint genes and biological processes relevant for local adaptation. Understanding adaptation ultimately requires providing functional evidence linking the candidate loci to their relevant phenotypic effects. However, these studies are extremely challenging because they require the identification of the specific phenotype to be examined and the particular experimental conditions in which the phenotype should be tested 72 . Epistasis and pleiotropy further complicate the establishment of a genotype-phenotype relationship 73,74 . Thus, studies as the one reported here that contribute to the identification of a set of candidate genetic variants likely to be under selection are needed in order to pinpoint the most likely candidates for follow-up functional validation. Our study, however, was limited to only two European populations from contrasting environments collected at a single time point. The European continent is very diverse in terms of habitats and climatic areas, and as such, a more complete analysis of different European populations should help to identify candidate genes and biological processes under spatially varying selection. Besides spatial variation, analyzing seasonal variation is another venue of research that might be worth exploring in light of recent results obtained in North American populations 28,75 .
We conclude that adding European populations to the available genome-wide datasets of D. melanogaster populations allowed us to disentangle selective from non-selective signatures of population differentiation in this species. European populations did not show signatures of admixture and showed low levels of population structure and thus are a good dataset to evaluate the contribution of SNPs and TEs to adaptive evolution. The analysis of the genetic variants identified in this work pinpoint genes and biological processes that have been relevant for the adaptation of this species to the out-of-Africa environments. Furthermore, we identified both large-effect and small-effect mutations suggesting that polygenic adaptation has also been relevant in this colonization process.

Material and Methods
Sample collection and newly sequenced genomes. Genomes reported here are derived from two natural population samples collected in Stockholm (59.34°, 17.94°, Sweden) and in Castellana Grotte (40.90°, 17.16°, Bari, Southern Italy) in September and October 2011, respectively ( Fig. 1 and Table S2, Supporting Information). Isofemale lines were established from a single wild-caught female and maintained in the laboratory for at least twenty generations at room temperature, and therefore passively inbred, before sequencing. A total of 26 Swedish and 16 Italian strains were individually sequenced. For each strain, DNA was extracted from 30 adult females using Puregene Core Kit A (Qiagen, Germany) according to the manufacturer's instructions. Short-sequence libraries (~500 bp) were constructed and paired-end sequencing using Illumina High-Seq 2000 platform was performed (100 bp reads). Quality control was performed using FastQC 76 . Raw fastq files were deposited in the Sequence Read Archive (SRA) from the National Center for Biotechnology Information (NCBI) under the BioProject accession: PRJNA390275.
Previously sequenced genomes. A total of 24 inbred strains collected in Raleigh (North Carolina, USA) from the DGRP panel 44,77 , and haploid embryos derived from 17 inbred or isofemale strains collected in Gikongoro (Rwanda, Africa) from the DPGP2 panel 43 were also analyzed (Table S2, Supporting Information). All these strains were sequenced to a similar depth and quality as the European strains reported in this work (Table S2, Supporting Information). The 17 Rwandan strains are a subset of the 22 strains available for this population with no inferred cosmopolitan admixture (nine strains) or <3% admixture (eight strains). For all these strains, raw reads in fastq format were downloaded from the NCBI-SRA. The global cohort consisted of a total of 83 genomes (26 Swedish, 16 Italian, 24 North American, and 17 African) that were analyzed using the same pipeline ( Variant Calling and Joint Genotyping. We followed the two-step GATK variant calling pipeline as follows. First, haplotype calling was performed independently for each strain using GATK HaplotypeCaller with a minimum base quality of 31 77 . We checked the distribution of heterozygous SNPs and they were distributed along the genome without any apparent region being enriched for heterozygous SNPs. Second, joint genotyping of the cohort of 83 genomes was performed using GATK GenotypeGVCF, which aggregates the gVCF files obtained from GATK HaplotypeCaller into a multi-sample VCF file. A total of 6,820,800 variants were called, including indels. Out of 6,820,800 variants, 5,479,546 were annotated as SNPs. SNPs falling in low-recombination regions close to centromeres and telomeres were excluded. The coordinates corresponding to euchromatic regions used in our analyses were defined by the recombination rate calculator (RRC) estimates 82  Genetic diversity and population structure analysis in the four worldwide populations. Regions of known cosmopolitan admixture in the African strains were excluded from the initial dataset of 4,106,958 SNPs, yielding a total of 3,749,072 SNPs 43 . Genome-wide diversity levels (π) were estimated for non-overlapping windows of 10,000 SNPs using vcftools v.0.1.12. The average nucleotide diversity was calculated for the whole genome and separately for the autosomes and the X chromosome using an ad hoc R script. We calculated nucleotide diversity considering that we were sampling two alleles. Because it is possible that in each isofemale strain genome we are sampling greater than or less than two alleles per locus, we also estimated nucleotide diversity after sampling one allele per genome with a weighted probability based on allelic read counts using an ad hoc python script. Global weighted Weir & Cockerham 83 unbiased estimator of the fixation index (F ST ) was estimated between each pair of populations using vcftools 83,84 . Population structure was also assessed by principal component analysis using Eigensoft and model-based global ancestry estimation using Admixture. SNPs with a minor allele frequency <15% or a genotyping rate lower than 70% were excluded. Two African strains with a variant calling lower than 70% (RG4, RG18) were also excluded. The final set consisted on 68,390 high confidence SNPs. Model-based global ancestry was estimated using Admixture software 85 . The number of ancestral populations analyzed (k) ranged from 1 to 7. Cross validation error for each value of k was used to determine the most likely hypothesis. The ancestry coefficients for each strain were plotted using the ggplot2 v1.0.0 library loaded in R v-3.1.0.
Testing for population differentiation in Europe. The following additional filters were applied to the initial set of 4,106,958 SNPs: (i) only non-singleton biallelic SNPs were considered; (ii) A custom perl script was used to filter to N individual calls with coverage lower than 10x or higher than the 0.975 percentile of the coverage per base pair distribution for the strain considered (see Table S2, Supporting Information); (iii) SNPs with a genotyping rate lower than 70% were excluded. In our dataset, this threshold maximizes the number of positions that pass the filtering while excluding those positions with an excess of missing genotypes. (iv) SNPs segregating at a frequency lower than 5% in the European population were excluded.
A total of 1,123,227 SNPs passed our conservative filtering pipeline in the European populations and were used to calculate the pairwise F ST . Significance of the heterogeneity in population differentiation levels across chromosomal arms was determined by performing one-way ANOVA. In subsequent analyses, the empirical distribution of F ST in each chromosomal arm was used to identify outliers. Genomic differentiation within the In(2 L)t common cosmopolitan inversion (2 L:2225744.0.13154180) 43 was compared to the levels observed in the rest of the chromosomal arm by performing a Mann-Whitney U test. We also estimated the population frequency of common and rare cosmopolitan inversions in the 42 European strains analysed in this work using the panel of SNPs known to be associated with these inversions 34 .
Over/Under representation of candidate differentiated loci per functional category. Functional annotation of SNPs was performed using VariantAnnotation R/Bioconductor package 86 . Out of the initial 1,123,227 SNPs, 1,121,948 could be annotated into the following functional categories: intergenic, promoter (1,000 bp upstream of the transcription start site), core promoter (316 bp upstream of the transcription start site 87 ), 5′ UTR, coding, splice site, intron, and 3′ UTR using TxDb.Dmelanogaster.UCSC.dm3.ensGene R/ Bioconductor package v.2.14.0. Coding variants were further classified into synonymous and non-synonymous. In order to obtain a set of neutrally evolving polymorphisms, we extracted those SNPs falling within the first 8-30 base pairs of small introns (< = 65 bp) 88 . The rest of intronic SNPs, which mostly belong to large introns, were classified as non-small intronic SNPs.
Counts non-random associations between functional category and belonging to the 95% of the chromosomal arm F ST distribution, or to the top 5%, 1% or 0.5% tails, while controlling for allele frequency, chromosomal arm, and considering SNPs in small introns as background of our analyses. The same procedure was used to determine if there was an over-or under-representation of candidate differentiated loci specifically located inside In(2 L)t inversion with respect to the rest of each chromosomal arm.

Assignment of Z ST scores to genes.
A population differentiation score per gene was obtained while controlling for gene length bias. SNP-wise F ST values were converted to gene-wise F ST values by assigning to each gene the maximum F ST value among those SNPs found along its transcribed region and the 1 kb region upstream of the transcription start site. This has been shown to be a good representation of the F ST per gene 45 . Maximum F ST values per gene were normalized to a Z score using chromosome specific empirical F ST distributions. There was a positive correlation between gene length (and number of SNPs) and Z score (r = 0.37132; p-value = 2.2e-16) (see Fig. S1A, Supporting Information). In order to correct for this potential bias, genes were assigned to 16 bins containing all genes with a similar number of SNPs and showing non-significant correlation between number of SNPs and Z score within each bin. We followed the approach detailed in 45 to obtain a median based corrected Z score that will be referred to as Z ST score. As a result, the correlation between the number of SNPs per gene and Z ST was eliminated (r = 0.01121; p-value = 0.1981) (Fig. S1B, Supporting Information). A total of 657 genes in the top 5% tail of the empirical Z ST distribution were considered as candidate differentiated genes.
Functional enrichment analysis. Gene Ontology (GO) term enrichment analysis in the list of 657 candidate differentiated genes was performed using topGO v.2.16.0 89 . We also performed gene set enrichment analysis by applying the Kolmogorov-Smirnov like test, which computes enrichment based on gene Z ST scores using GO's biological process annotation as gene set definition. In addition to the conventional Fisher's exact test or the Kolmogorov-Smirnov test (classic), we used three other algorithms: elim, weight and weight01 to identify significantly enriched GO biological processes. While in the classic enrichment test each GO term is tested independently, in the other three algorithms the hierarchical structure of the GO is taken into account so that the enrichment score incorporates information about the whole GO topology. This eliminates local similarities and dependencies between GO terms and reduces the high false-positive rate that is usually attained with classic GO enrichment tests 89 . The elim method iteratively removes genes annotated to significant GO terms from their ancestors, following a bottom-up approach. This allows for the identification of more specific nodes in the GO graph. The Weight algorithm identifies the locally most significant terms in the GO graph (and not necessarily the most specific terms). Weight01 is a combination of weight and elim algorithms. We decided to look at the results yielded by elim, weight and weight01 in order to explore which areas in the GO graph contain enriched GO terms instead of finding the precise GO terms that are enriched. Only GO terms with 5 or more differentiated genes were considered. We considered that a GO term was significantly enriched when the p-value computed by the elim, weight or weight01 algorithms was smaller than 0.05. These methods are internally reducing the significance of most GO terms and practically disregarding many of them, reducing the probability of having false positives due to the multiple-testing problem. As such, the raw p-values obtained with this methods can be considered as corrected p-values 89 .

TE frequency estimations.
We used T-lex version 2 (T-lex2) to estimate presence/absence of TEs in each individual strain from the Italian and Sweden population 65

Data Availability Statement
The datasets generated during the current study are available in the Sequence Read Archive (SRA) repository, under accession number BioProject PRJNA390275.