Genomic Analysis of European Drosophila melanogaster Populations Reveals Longitudinal Structure, Continent-Wide Selection, and Previously Unknown DNA Viruses Martin Kapun ,*†¶§,1,2,3 Maite G. Barr on,†,1,4 Fabian Staubach,‡,1,5 Darren J. Obbard,1,6 R. Axel W. Wiberg ,1,7,8 Jorge Vieira,1,9,10 Clement Goubert,1,11,12 Omar Rota-Stabelli,1,13 Maaria Kankare ,‡,1,14 Marıa Bogaerts-Marquez,1,4 Annabelle Haudry,1,11 Lena Waidele,1,5 Iryna Kozeretska,1,15,16 Elena G. Pasyukova,1,17 Volker Loeschcke,1,18 Marta Pascual,1,19,20 Cristina P. Vieira,1,9,10 Svitlana Serga,1,15 Catherine Montchamp-Moreau,1,21 Jessica Abbott,1,22 Patricia Gibert,1,11 Damiano Porcelli,1,23 Nico Posnien ,1,24 Alejandro Sanchez-Gracia,1,19,20 Sonja Grath,1,25 Elio Sucena,1,26,27 Alan O. Bergland,‡,1,28 Maria Pilar Garcia Guerreiro,1,29 Banu Sebnem Onder,1,30 Eliza Argyridou,1,25 Lain Guio,1,4 Mads Fristrup Schou,1,18,22 Bart Deplancke,1,31 Cristina Vieira,1,11 Michael G. Ritchie ,1,7 Bas J. Zwaan ,1,32 Eran Tauber,1,33,34 Dorcas J. Orengo,1,19,20 Eva Puerma,1,19,20 Montserrat Aguade,1,19,20 Paul Schmidt ,‡,1,35 John Parsch,1,25 Andrea J. Betancourt,1,36 Thomas Flatt ,*†‡1,2,3 and Josefa Gonzalez*†‡1,4 1 The European Drosophila Population Genomics Consortium (DrosEU) 2 Department of Ecology and Evolution, University of Lausanne, Lausanne, Switzerland 3 Department of Biology, University of Fribourg, Fribourg, Switzerland 4 Institute of Evolutionary Biology, CSIC-Universitat Pompeu Fabra, Barcelona, Spain 5 Department of Evolutionary Biology and Ecology, University of Freiburg, Freiburg, Germany 6 Institute of Evolutionary Biology, University of Edinburgh, Edinburgh, United Kingdom 7 Centre for Biological Diversity, School of Biology, University of St. Andrews, St Andrews, Scotland 8 Department of Environmental Sciences, Zoological Institute, University of Basel, Basel, Switzerland 9 Instituto de Biologia Molecular e Celular (IBMC), University of Porto, Porto, Portugal 10 Instituto de Investigaç~ao e Inovaç~ao em Sa ude (I3S), University of Porto, Porto, Portugal 11 Laboratoire de Biometrie et Biologie Evolutive UMR 5558, CNRS, Universite Lyon 1, Universite de Lyon, Villeurbanne, France 12 Department of Molecular Biology and Genetics, Cornell University, Ithaca, NY 13 Research and Innovation Centre, Fondazione Edmund Mach, San Michele all’ Adige, Italy 14 Department of Biological and Environmental Science, University of Jyv€askyl€a, Jyv€askyl€a, Finland Article 15 General and Medical Genetics Department, Taras Shevchenko National University of Kyiv, Kyiv, Ukraine 16 State Institution National Antarctic Scientific Center of Ministry of Education and Science of Ukraine, Kyiv, Ukraine 17 Laboratory of Genome Variation, Institute of Molecular Genetics of RAS, Moscow, Russia 18 Department of Bioscience—Genetics, Ecology and Evolution, Aarhus University, Aarhus C, Denmark 19 Departament de Genètica, Microbiologia i Estadıstica, Facultat de Biologia, Universitat de Barcelona, Barcelona, Spain 20 Institut de Recerca de la Biodiversitat (IRBio), Universitat de Barcelona, Barcelona, Spain 21 Universite Paris-Saclay, CNRS, IRD, UMR Evolution, Genomes, Comportement et Ecologie, 91198, Gif-sur-Yvette, France 22 Section for Evolutionary Ecology, Department of Biology, Lund University, Lund, Sweden 23 Department of Animal and Plant Sciences, Sheffield, United Kingdom 24 Johann-Friedrich-Blumenbach-Institut für Zoologie und Anthropologie, Universit€at Göttingen, Göttingen, Germany 25 Division of Evolutionary Biology, Faculty of Biology, Ludwig-Maximilians-Universit€at München, Planegg, Germany 26 Instituto Gulbenkian de Ci^encia, Oeiras, Portugal 27 Departamento de Biologia Animal, Faculdade de Ci^encias da Universidade de Lisboa, Lisboa, Portugal 28 Department of Biology, University of Virginia, Charlottesville, VA 29 Departament de Genètica i Microbiologia, Universitat Autonoma de Barcelona, Barcelona, Spain 30 Department of Biology, Faculty of Science, Hacettepe University, Ankara, Turkey 31 Institute of Bio-engineering, School of Life Sciences, EPFL, Lausanne, Switzerland ß The Author(s) 2020. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution. This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com Open Access Mol. Biol. Evol. 37(9):2661–2678 doi:10.1093/molbev/msaa120 Advance Access publication May 15, 2020 2661 Kapun et al. . doi:10.1093/molbev/msaa120 MBE 32 Laboratory of Genetics, Department of Plant Sciences, Wageningen University, Wageningen, Netherlands 33 Department of Evolutionary and Environmental Biology, University of Haifa, Haifa, Israel 34 Institute of Evolution, University of Haifa, Haifa, Israel 35 Department of Biology, University of Pennsylvania, Philadelphia, PA 36 Department of Evolution, Ecology, and Behaviour, University of Liverpool, Liverpool, United Kingdom † These authors contributed equally to this work. ‡ Drosophila Real Time Evolution (Dros-RTEC) Consortium. ¶ § Present address: Department of Evolutionary Biology and Environmental Sciences, University of Zürich, Zürich, Switzerland Present address: Division of Cell and Developmental Biology, Medical University of Vienna, Vienna, Austria *Corresponding authors: E-mails: martin.kapun@uzh.ch; thomas.flatt@unifr.ch; josefa.gonzalez@ibe.upf-csic.es. Associate editor: Daniel Falush Abstract Genetic variation is the fuel of evolution, with standing genetic variation especially important for short-term evolution and local adaptation. To date, studies of spatiotemporal patterns of genetic variation in natural populations have been challenging, as comprehensive sampling is logistically difficult, and sequencing of entire populations costly. Here, we address these issues using a collaborative approach, sequencing 48 pooled population samples from 32 locations, and perform the first continent-wide genomic analysis of genetic variation in European Drosophila melanogaster. Our anal- yses uncover longitudinal population structure, provide evidence for continent-wide selective sweeps, identify candidate genes for local climate adaptation, and document clines in chromosomal inversion and transposable element frequencies. We also characterize variation among populations in the composition of the fly microbiome, and identify five new DNA viruses in our samples. Key words: population genomics, adaptation, demography, selection, clines, SNPs, structural variants. Introduction wild samples, and a well-developed context for population Understanding processes that influence genetic variation in genomic analysis (Powell 1997; Keller 2007; Hales et al. 2015). natural populations is fundamental to understanding the Importantly, this species is studied by an extensive interna- process of evolution (Dobzhansky 1970; Lewontin 1974; tional research community, with a long history of developing Kreitman 1983; Kimura 1984; Hudson et al. 1987; shared resources (Hales et al. 2015; Bilder and Irvine 2017; McDonald and Kreitman 1991; Adrian and Comeron 2013). Haudry et al. 2020). Until recently, technological constraints have limited studies Our study complements and extends previous studies of of natural genetic variation to small regions of the genome genetic variation in D. melanogaster, both from its native and small numbers of individuals. With the development of range in Southeastern Africa and from its world-wide expan- population genomics, we can now analyze patterns of sion as a human commensal. The expansion into Europe is genome-wide genetic variation for large numbers of individ- thought to have occurred 4,100–19,000 years ago and into uals, with samples structured across space and time. As a North America and Australia in the last few centuries (David result, we have new insight into the evolutionary dynamics and Capy 1988; Lachaise et al. 1988; Li and Stephan 2006; of genetic variation in natural populations (Begun et al. 2007; Keller 2007; Kapopoulou, Kapun, et al. 2018; Arguello et al. Hohenlohe et al. 2010; Cheng et al. 2012; Pool et al. 2012; 2019; Sprengelmeyer et al. 2020). The colonization of novel Harpur et al. 2014; Zanini et al. 2015). But, despite this tech- habitats and climate zones on multiple continents makes nological progress, extensive large-scale sampling and genome D. melanogaster especially useful for studying parallel local sequencing of populations remains prohibitively expensive adaptation, with previous studies finding pervasive latitudinal and too labor-intensive for most individual research groups. clines in allele frequencies (Schmidt and Paaby 2008; Turner Here, we present the first comprehensive, continent-wide et al. 2008; Kolaczkowski et al. 2011; Fabian et al. 2012; genomic analysis of genetic variation of European Drosophila Bergland et al. 2014; Kapun, Fabian, et al. 2016; Machado melanogaster, based on 48 pool-sequencing samples from 32 et al. 2016), structural variants such as chromosomal inver- localities collected in 2014 (fig. 1) by the European Drosophila sions (reviewed in Kapun and Flatt 2019), transposable ele- Population Genomics Consortium (DrosEU; https://droseu. ments (TEs) (Boussy et al. 1998; Gonzalez et al. 2008, 2010), net). Drosophila melanogaster offers several advantages for and complex phenotypes (de Jong and Bochdanovits 2003; genomic studies of evolution in space and time. It boasts a Schmidt and Paaby 2008; Schmidt et al. 2008; Kapun, relatively small genome, a broad geographic range, a multi- Schmidt, et al. 2016; Behrman et al. 2018), especially along voltine life history which allows sampling across generations the North American and Australian east coasts. In addition to on short timescales, simple standard techniques for collecting parallel local adaptation, these latitudinal clines are, however, 2662 Genomic Analysis of European D. melanogaster Populations . doi:10.1093/molbev/msaa120 MBE 65 Summer based on SNPs, is provided in Supplementary Material Fall online (see supplementary results and table S1, Summer + Fall Supplementary Material online). For each sample, we es- 60 timated genome-wide levels of p, Watterson’s h, and Tajima’s D (corrected for pooling; Futschik and 55 Schlötterer 2010; Kofler et al. 2011). In brief, patterns of genetic variability and Tajima’s D were largely consistent Latitude 50 with what has been previously observed on other conti- nents (Fabian et al. 2012; Langley et al. 2012; Lack et al. 45 2015, 2016), and genetic diversity across the genome varies mainly with recombination rate (Langley et al. 2012). We also found little spatiotemporal variation 40 among European populations in overall levels of sequence variability (table 2). 35 Below we focus on the identification of selective sweeps, -10 0 10 20 30 40 previously unknown longitudinal population structure across the European continent, patterns of local adaptation and Longitude clines, and microbiota. FIG. 1. The geographic distribution of population samples. Locations of all samples in the 2014 DrosEU data set. The color of the circles Several Genomic Regions Show Signatures of indicates the sampling season for each location: ten of the 32 loca- Continent-Wide Selective Sweeps tions were sampled at least twice, once in summer and once in fall To identify genomic regions that have likely undergone selec- (see table 1 and supplementary table S1, Supplementary Material tive sweeps in European populations of D. melanogaster, we online). Note that some of the 12 Ukrainian locations overlap in used Pool-hmm (Boitard et al. 2013; see supplementary table the map. S2A, Supplementary Material online), which identifies candi- date sweep regions via distortions in the allele frequency also affected by admixture with flies from Africa and Europe spectrum. We ran Pool-hmm independently for each sample (Caracristi and Schlötterer 2003; Yukilevich and True 2008a, and identified several genomic regions that coincide with 2008b; Duchen et al. 2013; Kao et al. 2015; Bergland et al. previously identified, well-supported sweeps in the proximity 2016). of Hen1 (Kolaczkowski et al. 2011), Cyp6g1 (Daborn et al. In contrast, the population genomics of D. melanogaster 2002), wapl (Beisswanger et al. 2006), and around the chime- on the European continent remains largely unstudied ric gene CR18217 (Rogers and Hartl 2012), among others (Bozicevic et al. 2016; Pool et al. 2016; Mateo et al. 2018). (supplementary table S2B, Supplementary Material online). Because Eurasia was the first continent colonized by These regions also showed local reductions in p and Tajima’s D. melanogaster as they migrated out of Africa, we sought D, consistent with selective sweeps (fig. 2 and supplementary to understand how this species has adapted to new habitats figs. S1 and S2, Supplementary Material online). The putative and climate zones in Europe, where it has been established sweep regions that we identified in the European populations the longest (David and Capy 1988; Lachaise et al. 1988). We included 145 of the 232 genes previously identified using Pool- analyze our data at three levels: 1) variation at single- hmm in an Austrian population (Boitard et al. 2012; supple- nucleotide polymorphisms (SNPs) in nuclear and mitochon- mentary table S2C, Supplementary Material online). We also drial (mtDNA) genomes (5.5106 SNPs in total); 2) struc- identified other regions which have not previously been de- tural variation, including TE insertions and chromosomal scribed as targets of selective sweeps (supplementary table inversion polymorphisms; and 3) variation in the microbiota S2A, Supplementary Material online). Of the regions analyzed, associated with flies, including bacteria, fungi, protists, and 64 showed signatures of selection across all European popu- viruses. lations (supplementary table S2D, Supplementary Material online). Of these, 52 were located in the 10% of regions Results and Discussion with the lowest values of Tajima’s D (SuperExactTest; As part of the DrosEU consortium, we collected 48 population P < 0.001). These may represent continent-wide sweeps samples of D. melanogaster from 32 geographical locations that predate the colonization of Europe (Beisswanger et al. across Europe in 2014 (table 1 and fig. 1). We performed 2006) or which have recently swept across the majority of pooled sequencing (Pool-Seq) of all 48 samples, with an av- European populations (supplementary table S2D, erage autosomal coverage 50 (supplementary table S1, Supplementary Material online). Indeed, 43 of the 64 genes Supplementary Material online). Of the 32 locations, ten (67%) that showed signatures of selection across all European were sampled at least once in summer and once in fall populations were located in regions with reduced Tajima’s D (fig. 1), allowing a preliminary analysis of seasonal change in (lowest 10%) in African populations, suggesting that selective allele frequencies on a genome-wide scale. sweeps in these genes might predate the out-of-Africa expan- A description of the basic patterns of genetic variation sion (supplementary table S2D, Supplementary Material of these European D. melanogaster population samples, online). 2663 Kapun et al. . doi:10.1093/molbev/msaa120 MBE Table 1. Sample Information for All Populations in the DrosEU Data Set. Number Lat Lon Alt ID Country Location Coll. Date ID ( ) ( ) (m) Season n Coll. Name AT_Mau_14_01 Austria Mauternbach 2014-07-20 1 48.38 15.56 572 S 80 Andrea J. Betancourt AT_Mau_14_02 Austria Mauternbach 2014-10-19 2 48.38 15.56 572 F 80 Andrea J. Betancourt TR_Yes_14_03 Turkey Yesiloz 2014-08-31 3 40.23 32.26 680 S 80 Banu Sebnem Onder TR_Yes_14_04 Turkey Yesiloz 2014-10-23 4 40.23 32.26 680 F 80 Banu Sebnem Onder FR_Vil_14_05 France Viltain 2014-08-18 5 48.75 2.16 153 S 80 Catherine Montchamp- Moreau FR_Vil_14_07 France Viltain 2014-10-27 7 48.75 2.16 153 F 80 Catherine Montchamp- Moreau FR_Got_14_08 France Gotheron 2014-07-08 8 44.98 4.93 181 S 80 Cristina Vieira UK_She_14_09 United Kingdom Sheffield 2014-08-25 9 53.39 21.52 100 S 80 Damiano Porcelli UK_Sou_14_10 United Kingdom South Queensferry 2014-07-14 10 55.97 23.35 19 S 80 Darren Obbard CY_Nic_14_11 Cyprus Nicosia 2014-08-10 11 35.07 33.32 263 S 80 Eliza Argyridou UK_Mar_14_12 United Kingdom Market Harborough 2014-10-20 12 52.48 20.92 80 F 80 Eran Tauber UK_Lut_14_13 United Kingdom Lutterworth 2014-10-20 13 52.43 21.10 126 F 80 Eran Tauber DE_Bro_14_14 Germany Broggingen 2014-06-26 14 48.22 7.82 173 S 80 Fabian Staubach DE_Bro_14_15 Germany Broggingen 2014-10-15 15 48.22 7.82 173 F 80 Fabian Staubach UA_Yal_14_16 Ukraine Yalta 2014-06-20 16 44.50 34.17 72 S 80 Iryna Kozeretska UA_Yal_14_18 Ukraine Yalta 2014-08-27 18 44.50 34.17 72 S 80 Iryna Kozeretska UA_Ode_14_19 Ukraine Odesa 2014-07-03 19 46.44 30.77 54 S 80 Iryna Kozeretska UA_Ode_14_20 Ukraine Odesa 2014-07-22 20 46.44 30.77 54 S 80 Iryna Kozeretska UA_Ode_14_21 Ukraine Odesa 2014-08-29 21 46.44 30.77 54 S 80 Iryna Kozeretska UA_Ode_14_22 Ukraine Odesa 2014-10-10 22 46.44 30.77 54 F 80 Iryna Kozeretska UA_Kyi_14_23 Ukraine Kyiv 2014-08-09 23 50.34 30.49 179 S 80 Iryna Kozeretska UA_Kyi_14_24 Ukraine Kyiv 2014-09-08 24 50.34 30.49 179 F 80 Iryna Kozeretska UA_Var_14_25 Ukraine Varva 2014-08-18 25 50.48 32.71 125 S 80 Oleksandra Protsenko UA_Pyr_14_26 Ukraine Pyriatyn 2014-08-20 26 50.25 32.52 114 S 80 Oleksandra Protsenko UA_Dro_14_27 Ukraine Drogobych 2014-08-24 27 49.33 23.50 275 S 80 Iryna Kozeretska UA_Cho_14_28 Ukraine Chornobyl 2014-09-13 28 51.37 30.14 121 F 80 Iryna Kozeretska UA_Cho_14_29 Ukraine Chornobyl Yaniv 2014-09-13 29 51.39 30.07 121 F 80 Iryna Kozeretska SE_Lun_14_30 Sweden Lund 2014-07-31 30 55.69 13.20 51 S 80 Jessica Abbott DE_Mun_14_31 Germany Munich 2014-06-19 31 48.18 11.61 520 S 80 John Parsch DE_Mun_14_32 Germany Munich 2014-09-03 32 48.18 11.61 520 F 80 John Parsch PT_Rec_14_33 Portugal Recarei 2014-09-26 33 41.15 28.41 175 F 80 Jorge Vieira ES_Gim_14_34 Spain Gimenells (Lleida) 2014-10-20 34 41.62 0.62 173 F 80 Lain Guio ES_Gim_14_35 Spain Gimenells (Lleida) 2014-08-13 35 41.62 0.62 173 S 80 Lain Guio FI_Aka_14_36 Finland Akaa 2014-07-25 36 61.10 23.52 88 S 80 Maaria Kankare FI_Aka_14_37 Finland Akaa 2014-08-27 37 61.10 23.52 88 S 80 Maaria Kankare FI_Ves_14_38 Finland Vesanto 2014-07-26 38 62.55 26.24 121 S 66 Maaria Kankare DK_Kar_14_39 Denmark Karensminde 2014-09-01 39 55.95 10.21 15 F 80 Mads Fristrup Schou DK_Kar_14_41 Denmark Karensminde 2014-11-25 41 55.95 10.21 15 F 80 Mads Fristrup Schou CH_Cha_14_42 Switzerland Chalet a Gobet 2014-07-24 42 46.57 6.70 872 S 80 Martin Kapun CH_Cha_14_43 Switzerland Chalet a Gobet 2014-10-05 43 46.57 6.70 872 F 80 Martin Kapun AT_See_14_44 Austria Seeboden 2014-08-17 44 46.81 13.51 591 S 80 Martin Kapun UA_Kha_14_45 Ukraine Kharkiv 2014-07-26 45 49.82 36.05 141 S 80 Svitlana Serga UA_Kha_14_46 Ukraine Kharkiv 2014-09-14 46 49.82 36.05 141 F 80 Svitlana Serga UA_Cho_14_47 Ukraine Chornobyl Applegarden 2014-09-13 47 51.27 30.22 121 F 80 Svitlana Serga UA_Cho_14_48 Ukraine Chornobyl Polisske 2014-09-13 48 51.28 29.39 121 F 70 Svitlana Serga UA_Kyi_14_49 Ukraine Kyiv 2014-10-11 49 50.34 30.49 179 F 80 Svitlana Serga UA_Uma_14_50 Ukraine Uman 2014-10-01 50 48.75 30.21 214 F 80 Svitlana Serga RU_Val_14_51 Russia Valday 2014-08-17 51 57.98 33.24 217 S 80 Elena Pasyukova NOTE.—Origin, collection date, season, and sample size (number of chromosomes: n) of the 48 samples in the DrosEU 2014 data set. Additional information can be found in supplementary table S1, Supplementary Material online. We then asked if there was any indication of selective sweeps “response to stimulus,” “transport,” and “nervous system devel- particular to a certain habitat. To this end, we classified the opment.” For cold climates, it showed enrichment for “vitamin populations according to the Köppen–Geiger climate classifi- and cofactor metabolic processes” (supplementary table S2E, cation (Peel et al. 2007) and identified several putative sweeps Supplementary Material online). There was no enrichment of exclusive to arid, temperate, and cold regions (supplementary any GO category for sweeps associated with arid regions. table S2A, Supplementary Material online). To shed light on Thus, we identified several new candidate selective sweeps potential phenotypes affected by selective sweeps, we per- in European populations of D. melanogaster, many of which formed a gene ontology (GO) analysis. For temperate climates, occur in the majority of European populations and which this analysis showed enrichment for functions such as 2664 Genomic Analysis of European D. melanogaster Populations . doi:10.1093/molbev/msaa120 MBE Table 2. Clinality of Genetic Variation and Population Structure. Factor Latitude Longitude Altitude Season Moran’s I p(X) 4.11* 1.62 15.23*** 1.65 0.86 p(Aut) 0.91 2.54 27.18*** 0.16 20.86 h(X) 2.65 1.31 15.54*** 2.22 0.24 h(Aut) 0.48 1.44 13.66*** 0.37 21.13 D(X) 0.02 0.38 5.93* 3.26 22.08 D(Aut) 0.09 0.76 5.33* 0.71 21.45 PC1 0.63 118.08***(***) 3.64 0.75 4.2*** PC2 4.69* 7.15* 11.77** 1.68 20.32 PC3 0.39 0.23 19.91*** 0.28 1.38 NOTE.—Effects of geographic variables and/or seasonality on genome-wide average levels of diversity (p, h, and Tajima’s D; top rows) and on the first three axes of a PCA based on allele frequencies at neutrally evolving sites (bottom rows). The values represent F ratios from general linear models. Italic type indicates F ratios that are significant after Bonferroni correction (adjusted a0 ¼0.0055). Asterisks in parentheses indicate significance when accounting for spatial autocorrelation by spatial error models. These models were only calculated when Moran’s I test, as shown in the last column, was significant. *P < 0.05; **P < 0.01; ***P < 0.001. merit future study, using sequencing of individual flies and geographic distance, consistent with isolation by distance functional genetic experiments. (IBD; Mantel test; P < 0.001; R2¼0.12, max. FST  0.045; fig. 3A). Second, a principal components analysis (PCA) on European Populations Are Structured along an East– allele frequencies showed that the three most important PC West Gradient axes explain >25% of the total variance (PC1: 16.71%, PC2: We next investigated whether patterns of genetic differenti- 5.83%, PC3: 4.6%, eigenvalues ¼ 159.8, 55.7, and 44, respec- ation might be due to demographic substructuring. Overall, tively; fig. 3B). The first axis, PC1, was strongly correlated with pairwise differentiation as measured by FST was relatively low, longitude (F1,42 ¼ 118.08, P < 0.001; table 2). Again, this pat- particularly for the autosomes (autosomal FST 0.013–0.059; X- tern is consistent with IBD, as the European continent chromosome FST: 0.043–0.076; Mann–Whitney U test; extends further in longitude than latitude. We repeated the P < 0.001; supplementary table S1, Supplementary Material above PCA using SNPs in 4-fold degenerate sites, as these are online). The X chromosome is expected to have higher FST also assumed to be relatively unaffected by selection (Akashi than the autosomes, given its relatively smaller effective pop- 1995; Halligan and Keightley 2006; supplementary fig. S3, ulation size (Mann–Whitney U test; P < 0.001; Hutter et al. Supplementary Material online), and found highly consistent 2007). One population, from Sheffield (UK), was unusually results. differentiated from the others (average pairwise FST ¼ 0.027; Because there was a significant spatial autocorrelation be- SE¼ 0.00043 vs. FST ¼ 0.04; SE¼ 0.00055 for comparisons tween samples (as indicated by Moran’s test on residuals without this population and with this population only; sup- from linear regressions with PC1; P < 0.001; table 2), we re- plementary table S1, Supplementary Material online). peated the analysis with an explicit spatial error model; the Including this sample in the analysis could potentially lead association between PC1 and longitude remained significant. to exaggerated patterns of geographic differentiation, as it is To a lesser extent PC2 was likewise correlated with longitude both highly differentiated and the furthest west. We therefore (F1,42 ¼ 7.15, P < 0.05), but also with altitude (F1,42 ¼ 11.77, excluded it from the following analyses of geographic differ- P < 0.01) and latitude (F1,42 ¼ 4.69, P < 0.05; table 2). Similar entiation, as this approach is conservative. (For details, see to PC2, PC3 was strongly correlated with altitude (F1,42 ¼ Supplementary Material online; including or excluding this 19.91, P < 0.001; table 2). We also examined these data for population did not qualitatively change our results and their signatures of genetic differentiation between samples col- interpretation.) lected at different times of the year. For the data set as a Despite low overall levels of among-population differenti- whole, no major PC axes were correlated with season, indi- ation, we found that European populations exhibit clear ev- cating that there were no strong differences in allele frequen- idence of geographic substructuring. For this analysis, we cies shared between all our summer and fall samples (P > 0.05 focused on SNPs located within short introns, with a length for all analyses; table 2). For the ten locations sampled in both 60 bp and which most likely reflect neutral population summer and fall, we performed separate PC analyses for sum- structure (Haddrill et al. 2005; Singh et al. 2009; Parsch et al. mer and fall. Summer and fall values of PC1 (adjusted R2: 0.98; 2010; Clemente and Vogl 2012; Lawrie et al. 2013). We re- P < 0.001), PC2 (R2: 0.74; P < 0.001) and PC3 (R2: 0.81; stricted our analyses to polymorphisms in regions of high P < 0.001) were strongly correlated across seasons. This indi- recombination (r > 3 cM/Mb; Comeron et al. 2012) and to cates a high degree of seasonal stability in local genetic SNPs at least 1 Mb away from the breakpoints of common variation. inversions (and excluding the inversion bodies themselves), Next, we attempted to determine if populations could be resulting in 4,034 SNPs used for demographic analysis. statistically classified into clusters of similar populations. We found two signatures of geographic differentiation us- Using hierarchical model fitting based on the first four PC ing these putatively neutral SNPs. First, we identified a weak axes from the PCA mentioned above, we found two distinct but significant correlation between pairwise FST and clusters (fig. 3B) separated along PC1, supporting the notion 2665 Kapun et al. . doi:10.1093/molbev/msaa120 MBE A -1.5 -1.0 -0.5 0.0 0.5 D 12000000 12100000 12200000 12300000 exp CR44306 CG8888 CG18343 wash CG8854 CG42700 ana3 PI31 otk2 CG13186 SkpB CG33964 CR44313 CR45271 CG8839 ADD1 Mppe jeb CG13177 CR45276 CG13168 CG8850 CG30047 CG30036 pds5 Smyd4−4 Cyp6g1 CG13170 CG13164 CG30049 CG30037 CG8321 OSCP1 Cyp6g2 CG43315 CG17739 CG30043 GalT1 CG43191 Hen1 RpS11 CG43316 CG30203 CR33013 RnrS 128up CG8878 CG8858 CG43244 CG30046 CG34231 RpIII128 CG8407 Sr-CII Vha36-2 CR44316 rho−7 Drep1 Prp8 CG13171 Cam CG13163 CG8298 CG30039 CG34232 CR44315 CG8841 CR45319 CG44314 Oda CR43900 garz CR45318 Drep3 SmD3 CG8490 otk CR46194 CG34021 Cct5 Den1 EndoG CR45264 CG8860 SmF Cyp6t3 Tajima's D X 2L 2R 3L 3R 0.5 0.0 -0.5 -1.0 -1.5 -2.0 0 5 10 15 20 0 5 10 15 20 0 5 10 15 20 25 0 5 10 15 20 25 0 5 10 15 20 25 30 B -2.0 -1.5 -1.0 -0.5 0.0 0.5 D 16400000 16500000 16600000 16700000 16800000 CG44836 aos CG13033 l(3)73Ah Baldspot CG9692 CG42852 nudC fax CG43373 Syx8 Galphaf Dab CG13025 TMS1 st CG4229 CR18217 CR46033 CG13026 CG9705 GluRS-m Mipp1 Abl CR45897 CR46214 CG13024 TyrRS CG43206 Prosbeta6 CR45241 CG13031 CG9706 CG33158 mbf1 zetaCOP PGRP−SB2 Int6 dsx-c73A nxf2 CG13032 PGRP−SB1 CG9674 Smn CG43954 Nrt CR43727 CR45179 Rpn12 Dbp73D UQCR−C2 CG9701 CR45895 CG32163 tra spd-2 CG32165 CG32164 CG4098 Su(P) CG4101 CR45865 Mo25 FIG. 2. Candidate signals of selective sweeps in European populations. The central panel shows the distribution of Tajima’s D in 50-kb sliding windows with 40-kb overlap, with red and green dashed lines indicating Tajima’s D ¼ 0 and 1, respectively. The top panel shows a detail of a genomic region on chromosomal arm 2R in the vicinity of Cyp6g1 and Hen1 (highlighted in red), genes reportedly involved in pesticide resistance. This strong sweep signal is characterized by an excess of low-frequency SNP variants and overall negative Tajima’s D in all samples. Colored solid lines depict Tajima’s D for each sample (see supplementary fig. S2, Supplementary Material online, for color codes, Supplementary Material online); the black dashed line shows Tajima’s D averaged across all samples. The bottom panel shows a region on 3L previously identified as a potential target of selection, which shows a similar strong sweep signature. Notably, both regions show strongly reduced genetic variation (supplementary fig. S1, Supplementary Material online). of strong longitudinal differentiation among European pop- (Wright 1949) and using recent estimates of Ne in European ulations. Similarly, model-based spatial clustering also showed populations (Ne  3.1106; Duchen et al. 2013; for pairwise that populations were separated mainly by longitude (fig. 3C; migration rates see supplementary table S3, Supplementary using ConStruct, with K ¼ 3 spatial layers chosen based on Material online). Within the Western European cluster and model selection procedure via cross-validation). We also in- between the clusters, 4 Nem was similar (4 Nem-WE ¼ 43.76, ferred levels of admixture among populations from this anal- 4 Nem-between ¼ 45.97); in Eastern Europe, estimates of ysis, based on the relationship between FST and migration rate 4 Nem indicate significantly higher levels of admixture, despite 2666 Genomic Analysis of European D. melanogaster Populations . doi:10.1093/molbev/msaa120 MBE A 0 10 20 30 40 50 60 ●● ●● 60 60 ● ● Frequency ● ●● ● ●● Latitude Latitude ● ●● ● ●● ●● ● ●● ● 50 ● ● 50 ● ● ● ●●● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● 40 ● ● 40 ● ● ● 30 30 −10 0 10 20 30 40 0.01 0.03 0.05 −10 0 10 20 30 40 Longitude Pairwise FST Longitude B France -0.3 Portugal Dim 2 (5.83 %) France France Spain Ukraine Ukraine -0.1 Germany Ukraine Spain Ukraine Ukraine Switzerland Ukraine Ukraine Germany Ukraine Ukraine Denmark Ukraine Switzerland Austria Ukraine Ukraine UK UK Sweden Ukraine Austria Russia Ukraine Ukraine Ukraine Ukraine 0.1 Denmark UK Germany Ukraine Ukraine Germany Turkey Austria Turkey Finland Finland 0.3 Cyprus Finland -0.2 -0.1 0.0 0.1 0.2 Dim 1 (16.71 %) C Admixture proportion Cross-validation FR_Vil_14_05 UK_Mar_14_12 Predictive accuracy DE_Mun_14_31 UK_Sou_14_10 60 FI_Aka_14_37 -4 ES_Gim_14_34 UK_Lut_14_13 FR_Vil_14_07 DE_Bro_14_14 CH_Cha_14_43 SE_Lun_14_30 FI_Ves_14_38 55 FR_Got_14_08 -5 TR_Yes_14_04 PT_Rec_14_33 CH_Cha_14_42 Latitude FI_Aka_14_36 DE_Mun_14_32 CY_Nic_14_11 TR_Yes_14_03 50 DK_Kar_14_39 AT_See_14_44 -6 ES_Gim_14_35 DE_Bro_14_15 DK_Kar_14_41 UA_Kyi_14_49 UA_Dro_14_27 UA_Yal_14_16 45 UA_Cho_14_28 UA_Ode_14_19 UA_Kha_14_45 -8 UA_Kyi_14_23 UA_Ode_14_21 UA_Yal_14_18 UA_Ode_14_20 UA_Cho_14_48 40 AT_Mau_14_02 UA_Ode_14_22 UA_Kyi_14_24 -10 UA_Cho_14_29 UA_Uma_14_50 UA_Kha_14_46 RU_Val_14_51 35 AT_Mau_14_01 UA_Var_14_25 UA_Cho_14_47 UA_Pyr_14_26 -10 0 10 20 30 40 0.0 0.2 0.4 0.6 0.8 1.0 2 4 6 8 10 Longitude Admixture Values of K FIG. 3. Genetic differentiation among European populations. (A) Average FST among populations at putatively neutral sites. The center plot shows the distribution of FST values for all 1,128 pairwise population comparisons, with the FST values for each comparison obtained from the mean across all 4,034 SNPs used in the analysis. Plots on the left and the right show population pairs in the lower (blue) and upper (red) 5% tails of the FST distribution. (B) PCA analysis of allele frequencies at the same SNPs reveals population substructuring in Europe. Hierarchical model fitting using the first four PCs showed that the populations fell into two clusters (indicated by red and blue), with cluster assignment of each population subsequently estimated by k-means clustering. (C) Admixture proportions for each population inferred by model-based clustering with ConStruct are highlighted as pie charts (left plot) or structure plots (center). The optimal number of three spatial layers (K) was inferred by cross-validation (right plot). the larger geographic range covered by these samples (4 Nem WorldClim database (Hijmans et al. 2005), were associated ¼ 74.17; Mann–Whitney U test; P < 0.001). This result sug- with the genetic structure in our samples. These climatic gests that longitudinal differentiation in Europe might be variables represent interpolated averages across 30 years of partly driven by high levels of genetic exchange in Eastern observation at the geographic coordinates corresponding to Europe, perhaps due to migration and recolonization after our sampling locations. As many of these variables are highly harsh winters in that region. However, these estimates of gene intercorrelated, we analyzed their joint effects on genetic var- flow must be interpreted with caution, as unknown demo- iation, by using PCA to summarize the information they cap- graphic events can confound estimates of migration rates ture. The first three climatic PC axes capture >77% of the from FST (Whitlock and McCauley 1999). variance in the 19 climatic variables (supplementary table S4, In addition to restricted gene flow between geographic Supplementary Material online). PC1 explained 36% of the areas, local adaptation may explain population substructure, variance and was strongly correlated (r > 0.75 or r <0.75) even at neutral sites, if nearby and closely related populations with climatic variables differentiating “hot and dry” from are responding to similar selective pressures. We investigated “cold and wet” climates (e.g., maximum temperature of the whether any of 19 climatic variables, obtained from the warmest month, r ¼ 0.84; mean temperature of warmest 2667 Kapun et al. . doi:10.1093/molbev/msaa120 MBE X 2L 2R 3L 3R ∞ −log10(q−value) 6 PC1 4 In(2L)t In(2R)NS 2 In(3L)P ∞ In(3R)C −log10(q−value) In(3R)K 6 In(3R)Mo PC2 In(3R)P 4 2 0 10 20 30 0 10 20 30 0 10 20 30 0 10 20 30 0 10 20 30 Genomic Position (Mbp) twin Eip93F InR ETHR sgg cry Eip75B fru Eip74EF cpo sick Ilp3 Ubx Machado et al., 2016 Tor Clk Abd−B Samuel Egfr Eip63E foxo Fabian et al., 2012 tim Treh mth timeout FIG. 4. Manhattan plots of SNPs with q values <0.05 in BayeScEnv association tests with PC1 or PC2 of bioclimatic variables. Vertical lines denote the breakpoints of common inversions. The gene names highlight some candidate genes found in our study and which have previously been identified as varying clinally by Fabian et al. (2012) and Machado et al. (2016) along the North American east coast. Note that q values of 0 (which are infinite on a log-scale) are plotted at the top of each figure, above the gray dash-dotted horizontal lines in order to separate them from the other candidates with q values >0. These zero values are unlikely to be spurious as the densities of these infinite values tend to line up with peaks of log10(q) below the dashed line, suggesting that they represent highly significant continuations of these peaks. quarter, r ¼ 0.86; annual mean temperature, r ¼ 0.85; precip- Brodin 2005; Macholan et al. 2008; Knief et al. 2019), we found itation during the warmest quarter, r ¼ 0.87). Conversely, low genome-wide differentiation among eastern and western PC2 (27.3% of variance explained) distinguished climates with populations (average max. FST ~ 0.045), perhaps indicating low and high differences between seasons (e.g., isothermality, that the longitudinal division of European D. melanogaster is r ¼ 0.83; temperature seasonality, r ¼ 0.88; temperature an- not the result of postglacial secondary contact. nual range, r ¼0.78; precipitation in coldest quarter, r ¼ 0.79). PC1 was strongly correlated with latitude (linear Climatic Predictors Identify Genomic Signatures of regression: R2 ¼ 0.48, P < 0.001), whereas PC2 was strongly Local Climate Adaptation associated with longitude (R2 ¼ 0.58, P < 0.001). PC2 was also To further explore climatic patterns, and to identify signa- correlated with latitude (R2 ¼ 0.11, P < 0.05) and altitude (R2 tures of local adaptation caused by climatic differences ¼ 0.12, P < 0.01). among populations independent of neutral demographic We next asked whether any of these climate PCs explained effects, we tested for associations of SNP alleles with climatic any of the genetic structure uncovered above. Pairwise linear PC1 and PC2 using BayeScEnv (de Villemereuil and Gaggiotti regressions of the first three PC axes based on allele frequen- 2015). The total number of SNPs tested and the number of cies of intronic SNPs against the first three climatic PCs “top SNPs” (q value<0.05) are given in supplementary table revealed only one significant correlation after Bonferroni cor- S5A, Supplementary Material online. A large proportion of rection: between climatic PC2 (“seasonality”) versus genetic the top SNPs were intergenic (PC1: 33.5%; PC2: 32.2%) or PC1 (longitude; adjusted a ¼ 0.017; R2 ¼ 0.49, P < 0.001). intronic variants (PC1: 50.1%; PC2: 50.5%). Manhattan plots This suggests that longitudinal differentiation along the of q values for all SNPs are shown in figure 4. These figures European continent might be partly driven by the transition show some distinct “peaks” of highly differentiated SNPs from oceanic to continental climate, possibly leading to local along with some broader regions of moderately differentiated adaptation to gradual changes in temperature seasonality SNPs (fig. 4). For example, the circadian rhythm gene timeout and the severity of winter conditions. and the ecdysone signaling genes Eip74EF and Eip75B all lie Interestingly, the central European division into an eastern near peaks associated with climatic PC1 (“hot/dry” vs. “cold/ and a western clade of D. melanogaster closely resembles wet”; fig. 4, top panels). We note that the corresponding known hybrid zones of organisms which form closely related genes have been identified in previous studies of clinal (lati- pairs of sister taxa. These biogeographic patterns have been tudinal) differentiation in North American D. melanogaster associated with long-term reductions of gene flow between (Fabian et al. 2012; Machado et al. 2016). We found a signif- eastern and western populations during the last glacial max- icant overlap between genes associated with PC1 and PC2 imum, followed by postglacial recolonization of the continent (both of which are correlated with latitude) in our study and from southern refugia (Hewitt 1999). However, in contrast to candidate gene sets from these previous studies of latitudinal many of these taxa, which often exhibit pronounced pre- and clines (SuperExactTest; P < 0.001; Fabian et al. 2012; Machado postzygotic isolation (Szymura and Barton 1986; Haas and et al. 2016). For example, out of 1,974 latitudinally varying loci 2668 Genomic Analysis of European D. melanogaster Populations . doi:10.1093/molbev/msaa120 MBE G1.1 G1.2 G2.1 G2.2 G2.3 A B 1.0 0.8 FPKM (%) 0.6 0.4 0.2 0.0 AT_Mau_14_01 CY_Nic_14_11 UA_Ode_14_21 DE_Mun_14_31 DK_Kar_14_41 RU_Val_14_51 AT_Mau_14_02 UK_Sou_14_10 UK_Mar_14_12 UA_Ode_14_20 UA_Ode_14_22 SE_Lun_14_30 DE_Mun_14_32 CH_Cha_14_42 UA_Uma_14_50 FR_Vil_14_07 TR_Yes_14_03 FR_Got_14_08 UK_She_14_09 UK_Lut_14_13 UA_Yal_14_18 UA_Ode_14_19 UA_Kyi_14_23 UA_Dro_14_27 UA_Cho_14_28 UA_Cho_14_29 PT_Rec_14_33 FI_Aka_14_37 FI_Ves_14_38 DK_Kar_14_39 CH_Cha_14_43 UA_Cho_14_47 UA_Cho_14_48 UA_Kyi_14_49 UA_Yal_14_16 UA_Pyr_14_26 FI_Aka_14_36 UA_Kha_14_46 FR_Vil_14_05 DE_Bro_14_15 UA_Var_14_25 ES_Gim_14_35 UA_Kha_14_45 TR_Yes_14_04 DE_Bro_14_14 UA_Kyi_14_24 ES_Gim_14_34 AT_See_14_44 FIG. 5. Mitochondrial haplotypes. (A) Network showing the relationship of five common mitochondrial haplotypes. (B) Estimated frequency of each mitochondrial haplotype in 48 European samples. along the North American east coast identified by Fabian associated genes (1,205) were also shared between PC1 and et al. (2012), we found 403 (20%) and 505 (26%) of them PC2. Some genes have indeed been previously implicated in to also be associated with PC1 and PC2 in European popula- climatic and clinal adaptation, such as the circadian rhythm tions, respectively (supplementary table S5B and C, genes timeless, timeout, and clock, the sexual differentiation Supplementary Material online). Moreover, the BayeScEnv gene fruitless, and the couch potato locus which underlies the analysis and Pool-hmm analysis together identify four regions latitudinal cline in reproductive dormancy in North America with both climatic associations and evidence for continent- (Tauber et al. 2007; Schmidt et al. 2008; Fabian et al. 2012). wide selective sweeps (supplementary table S5B and C, Notably, these also include the major insulin signaling genes Supplementary Material online). Finally, four other insulin-like receptor (InR) and forkhead box subgroup O (foxo), BayeScEnv candidate genes were previously identified as tar- which have strong genomic and experimental evidence im- gets of selection in African and North American populations plicating these loci in clinal, climatic adaptation along the based on significant McDonald–Kreitman tests (Langley et al. North America east coast (Paaby et al. 2010, 2014; Fabian 2012; see supplementary table S5B and C, Supplementary et al. 2012; Durmaz et al. 2019). Thus, European populations Material online). share multiple potential candidate targets of selection with We next asked whether any insights into the targets of North American populations (cf. Fabian et al. 2012; Machado local selection could be gleaned from examining the functions et al. 2016; also see Bozicevic et al. 2016). We next turned to of genes near the BayeScEnv peaks. We examined annotated examining polymorphisms other than SNPs, that is, mito- features within 2 kb of significantly associated SNPs (PC1: chondrial haplotypes as well as inversion and TE 3,545 SNPs near 2,078 annotated features; PC2: 5,572 SNPs polymorphisms. near 2,717 annotated features; supplementary table S5B and C, Supplementary Material online). First, we performed a GO Mitochondrial Haplotypes Also Exhibit Longitudinal term analysis with GOwinda (Kofler and Schlötterer 2012) to Population Structure ask whether SNPs associated with climatic PCs are enriched Mitochondrial haplotypes also showed evidence of longitudi- for any gene functions. For PC1, we found no GO term en- nal demographic structure in European population. We iden- richment. For PC2, we found enrichment for “cuticle devel- tified two main alternative mitochondrial haplotypes in opment” and “UDP-glucosyltransferase activity.” Next, we Europe, G1 and G2, each with several subhaplotypes (G1.1 performed functional annotation clustering with DAVID and G1.2 and G2.1, G2.2, and G2.3). The two subtypes, G1.2 (v6.8; Huang et al. 2009), and identified 37 and 47 clusters and G2.1, are separated by 41 mutations (fig. 5A). The fre- with an enrichment score > 1.3 for PC1 and PC2, respectively quencies of the alternative G1 and G2 haplotypes varied (supplementary table S5D and E, Supplementary Material among populations between 35.1% and 95.6% and between online, Huang et al. 2009). PC1 was enriched for categories 4.4% and 64.9%, respectively (fig. 5B). Qualitatively, three such as “sex differentiation” and “response to nicotine,” types of European populations could be distinguished based whereas PC2 was enriched for functional categories such as on these haplotypes: 1) central European populations, with a “response to nicotine,” “integral component of membrane,” high frequency (>60%) of G1 haplotypes, 2) Eastern and “sensory perception of chemical stimulus” (supplemen- European populations in summer, with a low frequency tary table S5D and E, Supplementary Material online). (<40%) of G1 haplotypes, and 3) Iberian and Eastern We also asked whether the SNPs identified by BayeScEnv European populations in fall, with a frequency of G1 haplo- show consistent signatures of local adaptation. Many types between 40% and 60% (supplementary fig. S4, 2669 Kapun et al. . doi:10.1093/molbev/msaa120 MBE LTR Satellite SINE LINE others Simple_repeat DNA na rRNA Helitron Genome Proportion 0.20 0.15 0.10 0.05 0.00 AT_Mau_14_01 CY_Nic_14_11 UA_Ode_14_21 DE_Mun_14_31 DK_Kar_14_41 RU_Val_14_51 AT_Mau_14_02 UK_Sou_14_10 UK_Mar_14_12 UA_Ode_14_20 UA_Ode_14_22 SE_Lun_14_30 DE_Mun_14_32 CH_Cha_14_42 UA_Uma_14_50 TR_Yes_14_03 FR_Vil_14_07 FR_Got_14_08 UK_She_14_09 UK_Lut_14_13 UA_Yal_14_18 UA_Ode_14_19 UA_Kyi_14_23 UA_Dro_14_27 UA_Cho_14_28 UA_Cho_14_29 PT_Rec_14_33 FI_Aka_14_37 FI_Ves_14_38 DK_Kar_14_39 CH_Cha_14_43 UA_Cho_14_47 UA_Cho_14_48 UA_Kyi_14_49 UA_Yal_14_16 UA_Pyr_14_26 FI_Aka_14_36 UA_Kha_14_46 FR_Vil_14_05 DE_Bro_14_15 UA_Var_14_25 ES_Gim_14_35 UA_Kha_14_45 TR_Yes_14_04 DE_Bro_14_14 UA_Kyi_14_24 ES_Gim_14_34 AT_See_14_44 Cumulative Frequency 0.0 0.2 0.4 0.6 0.8 1.0 In(2L)t In(3R)Mo In(2R)NS In(3R)Payne In(3L)P In(3R)C FIG. 6. Geographic patterns of structural variants. The upper panel shows stacked bar plots with the relative abundances of transposable elements (TEs) in all 48 population samples. The proportion of each repeat class was estimated from sampled reads with dnaPipeTE (two samples per run, 0.1 coverage per sample). The lower panel shows stacked bar plots depicting absolute frequencies of six cosmopolitan inversions in all 48 population samples. Supplementary Material online). Analyses of mitochondrial (supplementary table S6, Supplementary Material online). haplotypes from a North American population (Cooper et al. LTRs have been previously described as being the most abun- 2015) as well as from world-wide samples (Wolff et al. 2016) dant TEs in the D. melanogaster genome (Kaminker et al. also revealed high levels of haplotype diversity. 2002; Bergman et al. 2006). Correspondingly, variation in Although there was no correlation between the frequency the proportion of LTRs best explained variation in total TE of G1 haplotypes and latitude, G1 haplotypes and longitude content (LINEþLTRþDNA) (Pearson’s r ¼ 0.87, P < 0.01, vs. were weakly but significantly correlated (r2 ¼ 0.10; P < 0.05). DNA r ¼ 0.58, P ¼ 0.0117, and LINE r ¼ 0.36, P < 0.01 and We thus divided the data set into an eastern and a western supplementary fig. S5A, Supplementary Material online). subset along the 20 meridian, corresponding to the division For each of the 1,630 euchromatic TE insertion sites an- of two major climatic zones, temperate (oceanic) versus cold notated in the D. melanogaster reference genome v.6.04, we (continental) (Peel et al. 2007). This split revealed a clear estimated the frequency at which a copy of the TE was pre- correlation (r2¼0.5; P < 0.001) between longitude and the sent at that site using T-lex2 (Fiston-Lavier et al. 2015; see frequency of G1 haplotypes, explaining as much as 50% of supplementary table S7, Supplementary Material online). On the variation in the western group (supplementary fig. S4B, an average, 56% were fixed in all samples. The remaining Supplementary Material online). Similarly, in eastern popula- polymorphic TEs mostly segregated at low frequency in all tions, longitude and the frequency of G1 haplotypes were samples (supplementary fig. S5B, Supplementary Material on- correlated (r2 ¼ 0.2; P < 0.001), explaining 20% of the var- line), potentially due to purifying selection (Gonzalez et al. iance (supplementary fig. S4B, Supplementary Material on- 2008; Petrov et al. 2011; Kofler et al. 2012; Cridland et al. 2013; line). Thus, these mitochondrial haplotypes appear to follow a Blumenstiel et al. 2014). However, 246 were present at inter- similar east–west population structure as observed for the mediate frequencies (>10% and <95%) and located in nuclear SNPs described above. regions of nonzero recombination (Fiston-Lavier et al. 2010; Comeron et al. 2012; see supplementary table S7, The Frequency of Polymorphic TEs Varies with Supplementary Material online). Although some of these Longitude and Altitude insertions might be segregating neutrally at transposition– To examine the population genomics of structural variants, selection balance (Charlesworth et al. 1994; see supplemen- we first focused on TEs. Similar to previous findings, the re- tary fig. S5B, Supplementary Material online), they are likely petitive content of the 48 samples ranged from 16% to 21% of enriched for candidate adaptive mutations (Rech et al. 2019). the nuclear genome size (Quesneville et al. 2005; fig. 6). The In each of the 48 samples, TE frequency and recombination vast majority of detected repeats were TEs, mostly long- rate were negatively correlated genome-wide (Spearman rank terminal repeat elements (LTRs; range 7.55–10.15%) and sum test; P < 0.01), as has also been previously reported for long interspersed nuclear elements (LINEs; range 4.18– D. melanogaster (Bartolome et al. 2002; Petrov et al. 2011; 5.52%), along with a few DNA elements (range 1.16–1.65%) Kofler et al. 2012). This remains true when fixed (population 2670 Genomic Analysis of European D. melanogaster Populations . doi:10.1093/molbev/msaa120 MBE Table 3. Clinality and/or Seasonality of Chromosomal Inversions. Factor Latitude Longitude Altitude Season Moran’s I In(2L)t 2.2 10.09** 43.94*** 0.89 20.92 In(2R)NS 0.25 14.43*** 2.88 2.43 1.25 In(3L)P 21.78*** 2.82 0.62 3.6 21.61 In(3R)C 18.5***(***) 0.75 1.42 0.04 2.79** In(3R)Mo 0.3 0.09 0.35 0.03 20.9 In(3R)Payne 43.47*** 0.66 1.69 1.55 20.89 NOTE.—The values represent F ratios from binomial generalized linear models to account for frequency data. Underlined italic type indicates deviance values that were significant after Bonferroni correction (adjusted a0 ¼0.0071). Asterisks in parentheses indicate significance when accounting for spatial autocorrelation by spatial error models. These models were only calculated when Moran’s I test, as shown in the last column, was significant. **P < 0.01; ***P < 0.001. frequency 95%) TE insertions were excluded from the anal- Inversions Exhibit Latitudinal and Longitudinal Clines ysis, although it was not statistically significant for some chro- in Europe mosomes and populations (supplementary table S8, Polymorphic chromosomal inversions, another class of struc- Supplementary Material online). In both cases, the correlation tural variants besides TEs, are well-known to exhibit pro- was stronger when broad-scale (Fiston-Lavier et al. 2010) nounced spatial (clinal) patterns in North American, rather than fine-scale (Comeron et al. 2012) recombination Australian, and other populations, possibly due to spatially rate estimates were used, indicating that the former may best varying selection (reviewed in Kapun and Flatt 2019; also see capture long-term population recombination patterns (see Mettler et al. 1977; Knibb et al. 1981; Lemeunier and Aulard supplementary materials and methods and table S8, 1992; Hoffmann and Weeks 2007; Fabian et al. 2012; Kapun Supplementary Material online). et al. 2014; Adrion et al. 2015; Rane et al. 2015; Kapun, Fabian, We next tested whether variation in TE frequencies among et al. 2016). However, in contrast to North America and samples was associated with spatially or temporally varying Australia, inversion clines in Europe remain very poorly char- factors. We focused on 111 TE insertions that segregated at acterized (Lemeunier and Aulard 1992; Kapun and Flatt intermediate frequencies, were located in nonzero recombi- 2019). We therefore sought to examine the presence and nation regions, and that showed an interquartile range (IQR) frequency of six cosmopolitan inversions (In(2L)t, In(2R)NS, >10 (see supplementary materials and methods, In(3L)P, In(3R)C, In(3R)Mo, In(3R)Payne) in our European Supplementary Material online). Of these insertions, 57 samples, using a panel of highly diagnostic inversion-specific were significantly associated with an at least one variable of marker SNPs, identified through sequencing of cytologically interest after multiple testing correction (supplementary ta- determined karyotypes by Kapun et al. (2014) (also see ble S9A, Supplementary Material online): 13 were significantly Kapun, Fabian, et al. 2016). All 48 samples were polymorphic associated with longitude, 13 with altitude, five with latitude, for one or more inversions (fig. 6). However, only In(2L)t seg- three with season, and 23 insertions with more than one of regated at substantial frequencies in most populations (aver- these variables (supplementary table S9A, Supplementary age frequency ¼ 20.2%); all other inversions were either Material online). These 57 TEs were mainly located inside absent or rare (average frequencies: In(2R)NS ¼ 6.2%, genes (42 out of 57; Fisher’s exact test, P > 0.05; supplemen- In(3L)P ¼ 4%, In(3R)C ¼ 3.1%, In(3R)Mo ¼2.2%, In(3R)Payne tary table S9A, Supplementary Material online). ¼ 5.7%) (cf. Kapun, Fabian, et al. 2016; Kapun and Flatt 2019). The 57 TEs significantly associated with these environmen- Despite their overall low frequencies, several inversions tal variables were enriched for two TE families: the LTR 297 showed pronounced clinality, in qualitative agreement with family with 11 copies, and the DNA pogo family with five findings from other continents (Lemeunier and Aulard 1992; copies (v2 values after Yate’s correction < 0.05; supplemen- Kapun and Flatt 2019). For the analyses below, we tested for tary table S9B, Supplementary Material online). Interestingly, potentially confounding effects of significant residual spatial 17 of the 57 TEs coincided with previously identified adaptive autocorrelation among samples; all of these tests were negative, candidate TEs, suggesting that our data set might be enriched except for In(3R)C (Moran’s I  0, P > 0.05 for all tests; table 3). for adaptive insertions (SuperExactTest, P < 0.001), several of We observed significant latitudinal clines in Europe for In(3L)P, which exhibit spatial frequency clines that deviate from neu- In(3R)C, and In(3R)Payne (binomial generalized linear model: tral expectation (SuperExactTest, P < 0.001, supplementary Inversion frequency  Latitude þ Longitude þ Altitude þ table S9A, Supplementary Material online; cf.; Rech et al. Season; effect of Latitude: P < 0.001 for all; see table 3). Clines 2019). Moreover, 18 of the 57 TEs also show significant cor- for In(3L)P and In(3R)Payne were similar between Europe and relations with either geographical or temporal variables in North America (with frequencies for both decreasing with lat- North American populations (SuperExactTest, P < 0.001, itude, P < 0.05; see supplementary table S10, Supplementary supplementary table S9A, Supplementary Material online; Material online). However, all inversions differed in their fre- cf. Lerat et al. 2019). quency at the same latitude between North America and Europe (P < 0.001 for the Latitude  Continent interaction; supplementary table S10, Supplementary Material online). 2671 Kapun et al. . doi:10.1093/molbev/msaa120 MBE 1.0 Relative Abundance 0.8 0.6 0.4 0.2 0.0 AT_Mau_14_01 CY_Nic_14_11 UA_Ode_14_21 DE_Mun_14_31 DK_Kar_14_41 RU_Val_14_51 AT_Mau_14_02 UK_Sou_14_10 UK_Mar_14_12 UA_Ode_14_20 UA_Ode_14_22 SE_Lun_14_30 DE_Mun_14_32 CH_Cha_14_42 UA_Uma_14_50 UK_She_14_09 TR_Yes_14_03 FR_Vil_14_07 FR_Got_14_08 UK_Lut_14_13 UA_Yal_14_18 UA_Ode_14_19 UA_Kyi_14_23 UA_Dro_14_27 UA_Cho_14_28 UA_Cho_14_29 PT_Rec_14_33 FI_Aka_14_37 FI_Ves_14_38 DK_Kar_14_39 CH_Cha_14_43 UA_Cho_14_47 UA_Cho_14_48 UA_Kyi_14_49 UA_Yal_14_16 UA_Pyr_14_26 FI_Aka_14_36 UA_Kha_14_46 FR_Vil_14_05 DE_Bro_14_15 UA_Var_14_25 ES_Gim_14_35 UA_Kha_14_45 TR_Yes_14_04 DE_Bro_14_14 UA_Kyi_14_24 ES_Gim_14_34 AT_See_14_44 Wolbachia Leishmania Saccharomyces unclassified Yersinia Gluconobacter Trypanosoma Klebsiella Escherichia Pectobacterium Gluconacetobacter Serratia Lactobacillus Pantoea Acinetobacter Acetobacter Pseudomonas Enterobacter Mycobacterium Salmonella FIG. 7. Microbiota associated with Drosophila. Relative abundance of Drosophila-associated microbes as assessed by MGRAST classified shotgun sequences. Microbes had to reach at least 3% relative abundance in one of the samples to be represented. Latitudinal inversion clines previously observed along the reads in our data set using MGRAST, which identifies and North American and Australian east coasts (supplementary counts short protein motifs (features) within reads (Meyer fig. S6 and table S10, Supplementary Material online; Kapun, et al. 2008). We examined 262 million reads in total. Of these, Fabian, et al. 2016) have been attributed to spatially varying most were assigned to Wolbachia (mean 53.7%; fig. 7 and selection, especially in the case of In(3R)Payne (Anderson et al. supplementary table S1, Supplementary Material online), a 2005; Umina et al. 2005; Kennington et al. 2006; Rako et al. well-known endosymbiont of Drosophila (Werren et al. 2008). 2006; Kapun, Fabian, et al. 2016; Kapun, Schmidt, et al. 2016; The abundance of Wolbachia protein features relative to Durmaz et al. 2018; Kapun and Flatt 2019). Similar to patterns other microbial protein features (relative abundance) varied in North America (Kapun, Fabian, et al. 2016), we observed strongly between samples, ranging from 8.8% in a sample that clinality of the three inversion polymorphisms was mark- from Ukraine to almost 100% in samples from Spain, edly stronger than for putatively neutral SNPs in short introns Portugal, Turkey, and Russia (supplementary table S12, (see supplementary table S11, Supplementary Material on- Supplementary Material online). Similarly, Wolbachia loads line), suggesting that these polymorphisms are maintained varied 100-fold between samples, as estimated from the ratio nonneutrally. Together, these findings suggest that latitudinal of Wolbachia protein features to Drosophila protein features inversion clines in Europe are shaped by spatially varying se- (supplementary table S12, Supplementary Material online). In lection, as they are in North America (Kapun, Fabian, et al. contrast to a previous study (Kriesner et al. 2016), there was 2016; Kapun and Flatt 2019). no evidence for clinality of Wolbachia loads (P ¼ 0.13, longi- We also detected longitudinal clines for In(2L)t and tude; P ¼ 0.41, latitude; Kendall’s rank correlation). However, In(2R)NS, with both polymorphisms decreasing in frequency these authors measured infection frequencies whereas we from east to west (see table 3; P < 0.01; also cf. Kapun and measured Wolbachia loads in pooled samples. Because the Flatt 2019). Longitudinal clines for these two inversions have frequency of infection does not necessarily correlate with also been found in North America (cf. Kapun and Flatt 2019). microbial loads measured in pooled samples, we might not One of these inversions, In(2L)t, also changed in frequency have been able to detect such a signal in our data. with altitude (table 3; P < 0.001). These longitudinal and al- Acetic acid bacteria of the genera Gluconobacter, titudinal inversion clines did, however, not deviate from neu- Gluconacetobacter, and Acetobacter were the second largest tral expectation (supplementary table S11, Supplementary group, with an average relative abundance of 34.4% among Material online). microbial protein features. Furthermore, we found evidence for the presence of several genera of Enterobacteria (Serratia, Yersinia, Klebsiella, Pantoea, Escherichia, Enterobacter, European Drosophila Microbiomes Contain Salmonella, and Pectobacterium). Serratia occurs only at low Entomophthora, Trypanosomatids, and Previously frequencies or is absent from most of our samples, but Unknown DNA Viruses reaches a very high relative abundance among microbial pro- The microbiota can affect life-history traits, immunity, hor- tein features in the Nicosia (Cyprus) summer collection monal physiology, and metabolic homeostasis of their fly (54.5%). This high relative abundance was accompanied by hosts (Martino et al. 2017; Trinder et al. 2017) and might an 80 increase in Serratia bacterial load. thus reveal interesting patterns of local adaptation. We there- We also detected several eukaryotic microorganisms, al- fore examined the bacterial, fungal, protist, and viral micro- though they were less abundant than the bacteria. We found biota sequence content of our samples. To do this, we trypanosomatids, previously reported to be associated with characterized the taxonomic origin of the non-Drosophila Drosophila in other studies (Wilfert et al. 2011; Chandler and 2672 Genomic Analysis of European D. melanogaster Populations . doi:10.1093/molbev/msaa120 MBE James 2013; Hamilton et al. 2015), in 16 of our samples, on an (25/48), followed by Linvill Road virus (7/48) and Viltain virus average representing 15% of all microbial protein features (5/48), with Esparto virus and the entomopox-like virus being identified in these samples. the rarest (2/48 and 1/48, respectively). Because Wolbachia Fungal protein features make up <3% of all but three can protect Drosophila from viruses (Teixeira et al. 2008), we samples (from Finland, Austria, and Turkey; supplementary hypothesized that Wolbachia loads might correlate negatively table S12, Supplementary Material online). This is somewhat with viral loads, but found no evidence of such a correlation surprising because yeasts are commonly found on rotting (P ¼ 0.83 Kallithea virus; P ¼ 0.76 Esparto virus; P ¼ 0.52 fruit, the main food substrate of D. melanogaster, and co- Viltain virus; P ¼ 0.96 Vesanto 1 virus; P ¼ 0.93 Vesanto 2 occur with flies (Barata et al. 2012; Chandler et al. 2012). virus; P ¼ 0.5 Linvill Road virus; Kendall’s rank correlation). This result suggests that, although yeasts can attract flies Perhaps this is because the Kallithea virus, the most prevalent and play a role in food choice (Becher et al. 2012; Buser virus in our data set, is not expected to be affected by et al. 2014), they might not be highly prevalent in or on D. Wolbachia (Palmer et al., 2018). Similarly, Shi et al. (2018) melanogaster bodies. One reason might be that they are ac- found no link between Wolbachia and the prevalence or tively digested and thus not part of the microbiome. We also abundance of RNA viruses in data from individual flies. found the fungal pathogen Entomophthora muscae in 14 The variation in bacterial microbiomes across space and samples, making up 0.18% of the reads (Elya et al. 2018). time reported here is analyzed in more detail in Wang et al. Our data also allowed us to identify DNA viruses. Only one (2020); this study suggests that some of this variation is struc- DNA virus has been previously described for D. melanogaster tured geographically (cf. Walters et al. 2020). Thus, micro- (Kallithea virus; Webster et al. 2015; Palmer et al. 2018) and biome composition might contribute to phenotypic only two additional ones from other Drosophilid species differences and local adaptation among populations (D. innubila Nudivirus, Unckless 2011; Invertebrate (Haselkorn et al. 2009; Richardson et al. 2012; Staubach Iridovirus 31 in D. obscura and D. immigrans, Webster et al. et al. 2013; Kriesner et al. 2016; Wang and Staubach 2018). 2016). In our data set, 2 million reads came from Kallithea nudivirus (Webster et al. 2015), allowing us to assemble the first complete Kallithea genome (>300-fold coverage in the Conclusions Ukrainian sample UA_Kha_14_46; GenBank accession num- Here, we have comprehensively sampled and sequenced ber KX130344). European populations of D. melanogaster for the first time We also found reads from five additional DNA viruses that (fig. 1). We find that European D. melanogaster populations were previously unknown (supplementary table S13, are longitudinally differentiated for putatively neutral SNPs, Supplementary Material online). First, around 1,000 reads mitochondrial haplotypes as well as for inversion and TE in- come from a novel nudivirus closely related to both sertion polymorphisms. Potentially adaptive polymorphisms Kallithea virus and to D. innubila nudivirus (Unckless 2011) also show this pattern, possibly driven by the transition from in sample DK_Kar_14_41 from Karensminde, Denmark (sup- oceanic to continental climate along the longitudinal axis of plementary table S13, Supplementary Material online). As the Europe. We note that this longitudinal differentiation quali- reads from this virus were insufficient to assemble the ge- tatively resembles the one observed for human populations nome, we identified a publicly available data set (SRR3939042: in Europe (Cavalli-Sforza 1966; Xiao et al. 2004; Francalacci 27 male D. melanogaster from Esparto, California; Machado and Sanna 2008; Novembre et al. 2008). Given that et al. 2016) with sufficient reads to complete the genome D. melanogaster is a human commensal (Keller 2007; (provisionally named “Esparto Virus”; KY608910). Second, Arguello et al. 2019), it is thus tempting to speculate that we also identified two novel Densoviruses (Parvoviridae). the demographic history of European populations might have The first is a relative of Culex pipiens densovirus, provisionally been influenced by past human migration. Outside Europe, named “Viltain virus,” found at 94-fold coverage in sample east–west structure has been previously found in sub- FR_Vil_14_07 (Viltain; KX648535). The second is “Linvill Road Saharan Africa populations of D. melanogaster, with the split virus,” a relative of Dendrolimus punctatus densovirus, repre- between eastern and western African populations having oc- sented by only 300 reads here, but with high coverage in data curred 70 ka (Michalakis and Veuille 1996; Aulard et al. set SRR2396966 from a North American sample of 2002; Kapopoulou, Pfeifer, et al. 2018), a period that coincides D. simulans, permitting assembly (KX648536; Machado with a wave of human migration from eastern into western et al. 2016). Third, we detected a novel member of the Africa (Nielsen et al. 2017). However, in contrast to the pro- Bidnaviridae family, “Vesanto virus,” a bidensovirus related nounced pattern observed in Europe, African east–west to Bombyx mori densovirus 3 with 900-fold coverage in structure is relatively weak, explaining only 2.7% of varia- sample FI_Ves_14_38 (Vesanto; KX648533 and KX648534). tion, and is primarily due to an inversion whose frequency Finally, in one sample (UA_Yal_14_16), we detected a sub- varies longitudinally. In contrast, our demographic analyses stantial number of reads from an Entomopox-like virus, which are based on SNPs located in >1 Mb distance from the break- we were unable to fully assemble (supplementary table S13, points of the most common inversions and excluding the Supplementary Material online). inversion bodies, making it unlikely that the longitudinal pat- Using a detection threshold of >0.1% of the Drosophila tern we observe is driven by inversions. genome copy number, the most commonly detected viruses Our extensive sampling was feasible only due to synergistic were Kallithea virus (30/48 of the pools) and Vesanto virus collaboration among many research groups. Our efforts in 2673 Kapun et al. . doi:10.1093/molbev/msaa120 MBE Europe are paralleled in North America by the Dros-RTEC Rwanda and 40 from Siavonga (Zambia). We mapped these consortium (Machado et al. 2019), with whom we are collab- data to the D. melanogaster reference genome using the same orating to compare population genomic data across conti- pipeline as for our own data above, and built consensus nents. Together, we have sampled both continents annually sequences for each haploid sample by only considering alleles since 2014; we aim to continue to sample and sequence with >0.9 allele frequencies. We converted consensus European and North American Drosophila populations with sequences to VCF and used VCFtools (Danecek et al. 2011) increasing spatiotemporal resolution in future years. With for downstream analyses. Additional steps in the mapping these efforts, we hope to provide a rich community resource and variant calling pipeline and further downstream analyses for biologists interested in molecular population genetics and of the data are detailed in supplementary materials and meth- adaptation genomics. ods, Supplementary Material online. Materials and Methods Supplementary Material Supplementary data are available at Molecular Biology and A detailed description of the Materials and Methods is pro- Evolution online. vided in supplementary materials and methods, Supplementary Material online; here, we give a brief overview of the data set and the basic methods used. The 2014 DrosEU Acknowledgments data set represents the most comprehensive spatiotemporal We acknowledge support of the publication fee by the CSIC sampling of European D. melanogaster populations to date Open Access Publication Support Initiative through its Unit (fig. 1 and supplementary table S1, Supplementary Material of Information Resources for Research (URICI). We thank online). It comprises 48 samples of D. melanogaster collected three anonymous reviewers and the editors for their helpful from 32 geographical locations across Europe at different time comments on a previous version of our article. We are grate- points in 2014 through a joint effort of 18 research groups. ful to the members of the DrosEU and Dros-RTEC consortia Collections were mostly performed with baited traps using a and to Dmitri Petrov (Stanford University) for support and standardized protocol (see supplementary materials and discussion. DrosEU is funded by a Special Topic Networks methods, Supplementary Material online). From each collec- (STN) grant from the European Society for Evolutionary tion, we pooled 33–40 wild-caught males. We used males as Biology (ESEB). Computational analyses were partially exe- they are more easily distinguishable morphologically from cuted at the Vital-IT bioinformatics facility of the University similar species than females. Despite our precautions, we of Lausanne (Switzerland), the computing facilities of the CC identified a low level of D. simulans contamination in our LBBE/PRABI in Lyon (France), the bwUniCluster of the state sequences; we computationally filtered these sequences of Baden–Württemberg (bwHPC), and the University of St from the data prior to further analysis (see Supplementary Andrews Bioinformatics Unit which is funded by a Wellcome Material online). To sequence these samples, we extracted Trust ISSF award (Grant No. 105621/Z/14/Z). We are also DNA and barcoded each sample, and sequenced the 40 flies grateful to Simon Boitard and Oscar Gaggiotti for their helpful per sample as a pool (Pool-Seq; Schlötterer et al. 2014), as technical advice on Pool-hmm and BayeScEnv analyses, re- paired-end fragments on a Illumina NextSeq 500 sequencer at spectively. This study was supported by University of Freiburg the Genomics Core Facility of Pompeu Fabra University. Research Innovation Fund 2014, Deutsche Samples were multiplexed in five batches of ten samples, Forschungsgemeinschaft (DFG) (Grant No. STA1154/4- except for one batch of eight samples (supplementary table 1Project 408908608 to F.S.); Academy of Finland (Grant S1, Supplementary Material online). Each multiplexed batch Nos. 268241 and 272927 to M.K.); Russian Foundation of was sequenced on four lanes at 50 raw coverage per Basic Research (Grant No. 15-54-46009 CT_a to E.G.P.); sample. The read length was 151 bp, with a median insert Danish Natural Science Research Council (Grant No. 4002- size of 348 bp (range 209–454 bp). Our genomic data set is 00113 to V.L.); Ministerio de Economia y Competitividad available under NCBI Bioproject accession PRJNA388788. (Grant Nos. CTM2017-88080 [AEI/FEDER, UE] and Sequences were processed and mapped to the CGL2013-42432-P to M.P. and M.P.G.G., D. melanogaster reference genome (v.6.12) and reference respectively); Centre National de la Recherche Scientifique sequences from common commensals and pathogens. Our (Grant No. UMR 9191 to C.M.-M.); Vetenskapsrådet (Grant bioinformatic pipeline is available at https://github.com/ Nos. 2011-05679 and 2015-04680 to J.A.); Emmy Noether capoony/DrosEU_pipeline (last accessed May 22, 2020). To Programme of the DFG (Grant No. PO 1648/3-1 to N.P.); call SNPs, we developed custom software (PoolSNP; see sup- National Institute of Health (Grant Nos. R35GM119686 and plementary materials and methods, Supplementary Material R01GM100366 to A.O.B. and P.S.S., respectively); Scientific online; https://github.com/capoony/PoolSNP, last accessed and Technological Research Council of Turkey (TUBITAK) May 22, 2020), using stringent heuristic parameters. In addi- (Grant No. 214Z238 to B.S.O.); Agence Nationale de la tion, we obtained genome sequences from African flies from Recherche Exhyb (Grant No. 14-CE19-0016 to C.V.); FP6, the Drosophila Genome Nexus (DGN; http://www.johnpool. Network of Excellence LifeSpan (Grant No. FP6 036894 to net/genomes.html, last accessed May 22, 2020; see supple- B.J.Z.); FP7, IDEAL (Grant No. FP7/2007-2011/259679 to mentary table S14, Supplementary Material online, for SRA B.J.Z.); Israel Science Foundation (Grant No. 1737/17 to accession numbers). We used data from 14 individuals from E.T.); Deutsche Forschungsgemeinschaft (Grant No. PA 903/ 2674 Genomic Analysis of European D. melanogaster Populations . doi:10.1093/molbev/msaa120 MBE 8-1 to J.P.); Austrian Science Fund (FWF) (Grant Nos. P32275 transposable element sequences in the Drosophila melanogaster ge- and P27048 to M.K. and A.J.B., respectively); Biotechnology nome. Genome Biol. 7(11):R112. Bilder D, Irvine KD. 2017. Taking Stock of the Drosophila Research and Biological Sciences Research Council (BBSRC) (Grant No. Ecosystem. Genetics 206(3):1227–1236. BB/P00685X/1 to A.J.B.); Swiss National Science Foundation Blumenstiel JP, Chen X, He M, Bergman CM. 2014. An age-of-allele test of (SNSF) (Grant Nos. PP00P3_133641, PP00P3_165836, and neutrality for transposable element insertions. Genetics 31003A_182262 to T.F.); European Commission (Grant No. 196(2):523–538. H2020-ERC-2014CoG-647900 to J.G.); Secretaria Boitard S, Kofler R, Françoise P, Robelin D, Schlötterer C, Futschik A. 2013. Pool-hmm: a Python program for estimating the allele fre- d’Universitats i Recerca, Department Economia i quency spectrum and detecting selective sweeps from next gener- Coneixement, Generalitat de Catalunya (Grant No. GRC ation sequencing of pooled samples. Mol Ecol Resour. 13(2):337–340. 2017 SGR 880 to J.G.); Ministerio de Economia y Boitard S, Schlötterer C, Nolte V, Pandey RV, Futschik A. 2012. Detecting Competitividad/FEDER (Grant No. BFU2014-57779-P to selective sweeps from pooled next-generation sequencing samples. J.G.); and Ministerio de Ciencia e Innovacion/AEI/FEDER, EU Mol Biol Evol. 29(9):2177–2186. Boussy IA, Itoh M, Itoh M, Rand D, Woodruff RC. 1998. Origin and decay (Grant No. BFU2017-82937-P to J.G.) of the P element-associated latitudinal cline in Australian Drosophila melanogaster. Genetica 104(1):45–57. References Bozicevic V, Hutter S, Stephan W, Wollstein A. 2016. Population genetic evidence for cold adaptation in European Drosophila melanogaster Adrian AB, Comeron JM. 2013. The Drosophila early ovarian transcrip- populations. Mol Ecol. 25(5):1175–1191. tome provides insight to the molecular causes of recombination rate Buser CC, Newcomb RD, Gaskett AC, Goddard MR. 2014. Niche con- variation across genomes. BMC Genomics 14(1):794. struction initiates the evolution of mutualistic interactions. Ecol Lett. Adrion JR, Hahn MW, Cooper BS. 2015. Revisiting classic clines in 17(10):1257–1264. Drosophila melanogaster in the age of genomics. Trends Genet. Caracristi G, Schlötterer C. 2003. Genetic differentiation between 31(8):434–444. American and European Drosophila melanogaster populations could Akashi H. 1995. Inferring weak selection from patterns of polymorphism be attributed to admixture of African alleles. Mol Biol Evol. and divergence at “silent” sites in Drosophila DNA. Genetics 20(5):792–799. 139:1067–1076. Cavalli-Sforza LL. 1966. Population structure and human evolution. Proc Anderson AR, Hoffmann AA, McKechnie SW, Umina PA, Weeks AR. R Soc Lond B Biol Sci. 164(995):362–379. 2005. The latitudinal cline in the In(3R)Payne inversion polymor- Chandler JA, Eisen JA, Kopp A. 2012. Yeast communities of diverse phism has shifted in the last 20 years in Australian Drosophila mel- Drosophila species: comparison of two symbiont groups in the anogaster populations. Mol Ecol. 14(3):851–858. same hosts. Appl Environ Microbiol. 78(20):7327–7336. Arguello JR, Laurent S, Clark AG. 2019. Demographic history of the hu- Chandler JA, James PM. 2013. Discovery of trypanosomatid parasites in man commensal Drosophila melanogaster. Genome Biol Evol. globally distributed Drosophila species. PLoS One 8(4):e61937. 11(3):844–854. Charlesworth B, Sniegowski P, Stephan W. 1994. The evolutionary dy- Aulard S, David JR, Lemeunier F. 2002. Chromosomal inversion polymor- namics of repetitive DNA in eukaryotes. Nature 371(6494):215–220. phism in Afrotropical populations of Drosophila melanogaster. Genet Cheng C, White BJ, Kamdem C, Mockaitis K, Costantini C, Hahn MW, Res. 79(1):49–63. Besansky NJ. 2012. Ecological genomics of Anopheles gambiae along Barata A, Santos SC, Malfeito-Ferreira M, Loureiro V. 2012. New insights a latitudinal cline: a population-resequencing approach. Genetics into the ecological interaction between grape berry microorganisms 190(4):1417–1432. and Drosophila flies during the development of sour rot. Microb Ecol. Clemente F, Vogl C. 2012. Unconstrained evolution in short introns? – 64(2):416–430. An analysis of genome-wide polymorphism and divergence data Bartolome C, Maside X, Charlesworth B. 2002. On the abundance and from Drosophila. J Evol Biol. 25(10):1975–1990. distribution of transposable elements in the genome of Drosophila Comeron JM, Ratnappan R, Bailin S. 2012. The many landscapes of melanogaster. Mol Biol Evol. 19(6):926–937. recombination in Drosophila melanogaster. PLoS Genet. Becher PG, Flick G, RozpeR dowska E, Schmidt A, Hagman A, Lebreton S, 8(10):e1002905. Larsson MC, Hansson BS, Piskur J, Witzgall P, et al. 2012. Yeast, not Cooper BS, Burrus CR, Ji C, Hahn MW, Montooth KL. 2015. Similar fruit volatiles mediate Drosophila melanogaster attraction, oviposi- efficacies of selection shape mitochondrial and nuclear genes in tion and development. Funct Ecol. 26(4):822–828. both Drosophila melanogaster and Homo sapiens. G3 (Bethesda) Begun DJ, Holloway AK, Stevens K, Hillier LW, Poh Y-P, Hahn MW, Nista 5:2165–2176. PM, Jones CD, Kern AD, Dewey CN, et al. 2007. Population genomics: Cridland JM, Macdonald SJ, Long AD, Thornton KR. 2013. Abundance whole-genome analysis of polymorphism and divergence in and distribution of transposable elements in two Drosophila QTL Drosophila simulans. PLoS Biol. 5(11):e310. mapping resources. Mol Biol Evol. 30(10):2311–2327. Behrman EL, Howick VM, Kapun M, Staubach F, Bergland AO, Petrov Daborn PJ, Yen JL, Bogwitz MR. 2002. A single p450 allele associated with DA, Lazzaro BP, Schmidt PS. 2018. Rapid seasonal evolution in innate insecticide resistance in Drosophila. Science 297(5590):2253–2256. immunity of wild Drosophila melanogaster. Proc R Soc B. Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, 285(1870):20172599. Handsaker RE, Lunter G, Marth GT, Sherry ST, et al.. 2011. The Beisswanger S, Stephan W, De Lorenzo D. 2006. Evidence for a selective variant call format and VCFtools. Bioinformatics 27(15):2156–2158. sweep in the wapl region of Drosophila melanogaster. Genetics David JR, Capy P. 1988. Genetic variation of Drosophila melanogaster 172(1):265–274. natural populations. Trends Genet. 4(4):106–111. Bergland AO, Behrman EL, O’Brien KR, Schmidt PS, Petrov DA. 2014. de Jong G, Bochdanovits Z. 2003. Latitudinal clines in Drosophila mela- Genomic evidence of rapid and stable adaptive oscillations over nogaster: body size, allozyme frequencies, inversion frequencies, and seasonal time scales in Drosophila. PLoS Genet. 10(11):e1004775. the insulin-signalling pathway. J Genet. 82(3):207–223. Bergland AO, Tobler R, Gonzalez J, Schmidt P, Petrov D. 2016. Secondary Dobzhansky T. 1970. Genetics of the evolutionary process. New York: contact and local adaptation contribute to genome-wide patterns of Columbia University Press. clinal variation in Drosophila melanogaster. Mol Ecol. Duchen P, Zivkovic D, Hutter S, Stephan W, Laurent S. 2013. 25(5):1157–1174. Demographic inference reveals African and European admixture Bergman CM, Quesneville H, Anxolabehere D, Ashburner M. 2006. in the North American Drosophila melanogaster population. Recurrent insertion and duplication generate networks of Genetics 193(1):291–301. 2675 Kapun et al. . doi:10.1093/molbev/msaa120 MBE Durmaz E, Benson C, Kapun M, Schmidt P, Flatt T. 2018. An inversion Huang DW, Sherman BT, Lempicki RA. 2009. Systematic and integrative supergene in Drosophila underpins latitudinal clines in survival traits. analysis of large gene lists using DAVID bioinformatics resources. Nat J Evol Biol. 31(9):1354–1364. Protoc. 4(1):44–57. Durmaz E, Rajpurohit S, Betancourt N, Fabian DK, Kapun M, Schmidt P, Hudson RR, Kreitman M, Aguade M. 1987. A test of neutral molecular Flatt T. 2019. A clinal polymorphism in the insulin signaling tran- evolution based on nucleotide data. Genetics 116(1):153–159. scription factor foxo contributes to life-history adaptation in Hutter S, Li H, Beisswanger S, De Lorenzo D, Stephan W. 2007. Distinctly Drosophila. Evolution 73(9):1774–1792. different sex ratios in African and European populations of Elya C, Lok TC, Spencer QE, McCausland H, Martinez CC, Eisen MB. 2018. Drosophila melanogaster inferred from chromosomewide single nu- Robust manipulation of the behavior of Drosophila melanogaster by cleotide polymorphism data. Genetics 177(1):469–480. a fungal pathogen in the laboratory. eLife 7:e34414. Kaminker JS, Bergman CM, Kronmiller B, Carlson J, Svirskas R, Patel S, Fabian DK, Kapun M, Nolte V, Kofler R, Schmidt PS, Schlötterer C, Flatt Frise E, Wheeler DA, Lewis SE, Rubin GM, et al. 2002. The transpos- T. 2012. Genome-wide patterns of latitudinal differentiation among able elements of the Drosophila melanogaster euchromatin: a geno- populations of Drosophila melanogaster from North America. Mol mics perspective. Genome Biol. 3(12):research0084. Ecol. 21(19):4748–4769. Kao JY, Zubair A, Salomon MP, Nuzhdin SV, Campo D. 2015. Fiston-Lavier A-S, Barron MG, Petrov DA, Gonzalez J. 2015. T-lex2: gen- Population genomic analysis uncovers African and European ad- otyping, frequency estimation and re-annotation of transposable mixture in Drosophila melanogaster populations from the south- elements using single or pooled next-generation sequencing data. eastern United States and Caribbean Islands. Mol Ecol. Nucleic Acids Res. 43(4):e22. 24(7):1499–1509. Fiston-Lavier A-S, Singh ND, Lipatov M, Petrov DA. 2010. Drosophila Kapopoulou A, Kapun M, Pavlidis P, et al. 2018. Early split between melanogaster recombination rate calculator. Gene 463(1–2):18–20. African and European populations of Drosophila melanogaster. Francalacci P, Sanna D. 2008. History and geography of human Y-chro- bioRxiv. doi: https://doi.org/10.1101/340422. mosome in Europe: a SNP perspective. J Anthropol Sci. 86:59–89. Kapopoulou A, Pfeifer S, Jensen J, Laurent S. 2018. The demographic Futschik A, Schlötterer C. 2010. The next generation of molecular history of African Drosophila melanogaster. Genome Biol Evol. markers from massively parallel sequencing of pooled DNA samples. 10(9):2338–2342. Genetics 186(1):207–218. Kapun M, Fabian DK, Goudet J, Flatt T. 2016. Genomic evidence for Gonzalez J, Karasov TL, Messer PW, Petrov DA. 2010. Genome-wide adaptive inversion clines in Drosophila melanogaster. Mol Biol Evol. patterns of adaptation to temperate environments associated with 33(5):1317–1336. transposable elements in Drosophila. PLoS Genet. 6(4):e1000905. Kapun M, Flatt T. 2019. The adaptive significance of chromosomal in- Gonzalez J, Lenkov K, Lipatov M, Macpherson JM, Petrov DA. 2008. High version polymorphisms in Drosophila melanogaster. Mol Ecol. rate of recent transposable element–induced adaptation in 28(6):1263–1282. Drosophila melanogaster. PLoS Biol. 6(10):e251. Kapun M, Schmidt C, Durmaz E, Schmidt PS, Flatt T. 2016. Parallel effects Haas F, Brodin A. 2005. The Crow Corvus corone hybrid zone in southern of the inversion In(3R)Payne on body size across the North Denmark and northern Germany. Ibis. 147(4):649–656. American and Australian clines in Drosophila melanogaster. J Evol Haddrill PR, Charlesworth B, Halligan DL, Andolfatto P. 2005. Patterns of Biol. 29(5):1059–1072. intron sequence evolution in Drosophila are dependent upon length Kapun M, van Schalkwyk H, McAllister B, Flatt T, Schlötterer C. 2014. and GC content. Genome Biol. 6(8):R67. Inference of chromosomal inversion dynamics from Pool-Seq data in Hales KG, Korey CA, Larracuente AM, Roberts DM. 2015. Genetics on natural and laboratory populations of Drosophila melanogaster. Mol the fly: a primer on the Drosophila model system. Genetics Ecol. 23(7):1813–1827. 201(3):815–842. Keller A. 2007. Drosophila melanogaster’s history as a human commen- Halligan DL, Keightley PD. 2006. Ubiquitous selective constraints in the sal. Curr Biol. 17(3):R77–R81. Drosophila genome revealed by a genome-wide interspecies com- Kennington JW, Partridge L, Hoffmann AA. 2006. Patterns of diversity parison. Genome Res. 16(7):875–884. and linkage disequilibrium within the cosmopolitan inversion Hamilton PT, Votypka J, Dostalova A, Yurchenko V, Bird NH, Lukes J, In(3R)Payne in Drosophila melanogaster are indicative of coadapta- Lemaitre B, Perlman SJ. 2015. Infection Dynamics and Immune tion. Genetics 172(3):1655–1663. Response in a Newly Described Drosophila-Trypanosomatid Kimura M. 1984. The neutral theory of molecular evolution. Cambridge, Association. mBio. 6(5):e01356–01315. United Kingdom: Cambridge University Press. Harpur BA, Kent CF, Molodtsova D, Lebon JMD, Alqarni AS, Owayss AA, Knibb WR, Oakeshott JG, Gibson JB. 1981. Chromosome inversion poly- Zayed A. 2014. Population genomics of the honey bee reveals strong morphisms in Drosophila melanogaster. I. Latitudinal clines and signatures of positive selection on worker traits. Proc Natl Acad Sci U associations between inversions in Australasian populations. S A. 111(7):2614–2619. Genetics 98(4):833–847. Haselkorn TS, Markow TA, Moran NA. 2009. Multiple introductions of Knief U, Bossu CM, Saino N, Hansson B, Poelstra J, Vijay N, the Spiroplasma bacterial endosymbiont into Drosophila. Mol Ecol. Weissensteiner M, Wolf J. 2019. Epistatic mutations under divergent 18(6):1294–1305. selection govern phenotypic variation in the crow hybrid zone. Nat Haudry A, Laurent S, Kapun M. 2020. Population genomics on the fly: Ecol Evol. 3(4):570–576. recent advances in Drosophila. In: Dutheil JY, editor. Statistical pop- Kofler R, Betancourt AJ, Schlötterer C. 2012. Sequencing of pooled DNA ulation genomics. Vol. 2090. New York: Springer. p. 357–396. samples (Pool-Seq) uncovers complex dynamics of transposable el- Hewitt GM. 1999. Post-glacial re-colonization of European biota. Biol J ement insertions in Drosophila melanogaster. PLoS Genet. Linn Soc. 68(1–2):87–112. 8(1):e1002487. Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. 2005. Very high Kofler R, Orozco-terWengel P, De Maio N, Pandey RV, Nolte V, Futschik resolution interpolated climate surfaces for global land areas. Int J A, Kosiol C, Schlötterer C. 2011. PoPoolation: a toolbox for popula- Climatol. 25(15):1965–1978. tion genetic analysis of next generation sequencing data from Hoffmann AA, Weeks AR. 2007. Climatic selection on genes and traits pooled individuals. PLoS One 6(1):e15925. after a 100 year-old invasion: a critical look at the temperate-tropical Kofler R, Schlötterer C. 2012. GOwinda: unbiased analysis of gene set clines in Drosophila melanogaster from eastern Australia. Genetica enrichment for genome-wide association studies. Bioinformatics 129(2):133–147. 28(15):2084–2085. Hohenlohe PA, Bassham S, Etter PD, Stiffler N, Johnson EA, Cresko WA. Kolaczkowski B, Kern AD, Holloway AK, Begun DJ. 2011. Genomic dif- 2010. Population genomics of parallel adaptation in threespine stick- ferentiation between temperate and tropical Australian populations leback using sequenced RAD tags. PLoS Genet. 6(2):e1000862. of Drosophila melanogaster. Genetics 187(1):245–260. 2676 Genomic Analysis of European D. melanogaster Populations . doi:10.1093/molbev/msaa120 MBE Kreitman M. 1983. Nucleotide polymorphism at the alcohol Paaby AB, Bergland AO, Behrman EL, Schmidt PS. 2014. A highly pleio- dehydrogenase locus of Drosophila melanogaster. Nature tropic amino acid polymorphism in the Drosophila insulin receptor 304(5925):412–417. contributes to life-history adaptation. Evolution 68(12):3395–3409. Kriesner P, Conner WR, Weeks AR, Turelli M, Hoffmann AA. 2016. Paaby AB, Blacket MJ, Hoffmann AA, Schmidt PS. 2010. Identification of Persistence of a Wolbachia infection frequency cline in Drosophila a candidate adaptive polymorphism for Drosophila life history by melanogaster and the possible role of reproductive dormancy. parallel independent clines on two continents. Mol Ecol. Evolution 70(5):979–997. 19(4):760–774. Lachaise D, Cariou M-L, David JR, et al. 1988. Historical biogeography of Palmer WH, Medd NC, Beard PM, Obbard DJ. 2018. Isolation of a natural the Drosophila melanogaster species subgroup. In: Hecht MK, DNA virus of Drosophila melanogaster, and characterisation of host Wallace B, Prance GT, editors. Evolutionary biology. Boston: resistance and immune responses. PLoS Pathog. 14(6):e1007050. Springer. p. 159–225. Parsch J, Novozhilov S, Saminadin-Peter SS, Wong KM, Andolfatto P. Lack JB, Cardeno CM, Crepeau MW, Taylor W, Corbett-Detig RB, Stevens 2010. On the utility of short intron sequences as a reference for the KA, Langley CH, Pool JE. 2015. The Drosophila Genome Nexus: a detection of positive and negative selection in Drosophila. Mol Biol population genomic resource of 623 Drosophila melanogaster Evol. 27(6):1226–1234. genomes, including 197 from a single ancestral range population. Peel MC, Finlayson BL, McMahon TA. 2007. Updated world map of the Genetics 199(4):1229–1241. Köppen-Geiger climate classification. Hydrol Earth Syst Sci. Lack JB, Lange JD, Tang AD, Corbett-Detig RB, Pool JE. 2016. A thousand 11(5):1633–1644. fly genomes: an expanded Drosophila Genome Nexus. Mol Biol Evol. Petrov DA, Fiston-Lavier AS, Lipatov M, Lenkov K, Gonzalez J. 2011. 33(12):3308–3313. Population genomics of transposable elements in Drosophila mela- Langley CH, Stevens K, Cardeno C, Lee YCG, Schrider DR, Pool JE, Langley nogaster. Mol Biol Evol. 28(5):1633–1644. SA, Suarez C, Corbett-Detig RB, Kolaczkowski B, et al. 2012. Genomic Pool JE, Braun DT, Lack JB. 2016. Parallel evolution of cold tolerance variation in natural populations of Drosophila melanogaster. Genetics within Drosophila melanogaster. Mol Biol Evol. 34:349–360. 192(2):533–598. Pool JE, Corbett-Detig RB, Sugino RP, Stevens KA, Cardeno CM, Crepeau Lawrie DS, Messer PW, Hershberg R, Petrov DA. 2013. Strong purifying MW, Duchen P, Emerson JJ, Saelao P, Begun DJ, et al. 2012. selection at synonymous sites in D. melanogaster. PLoS Genet. Population genomics of sub-Saharan Drosophila melanogaster: 9(5):e1003527. African diversity and non-African admixture. PLoS Genet. Lemeunier F, Aulard S. 1992. Inversion polymorphism in Drosophila 8(12):e1003080. melanogaster. In: Krimbas CB, Powell JR, editors. Drosophila inversion Powell JR. 1997. Progress and prospects in evolutionary biology: the polymorphism. New York: CRC Press. p. 339–405. Drosophila model. New York and Oxford: Oxford University Press. Lerat E, Goubert C, Guirao-Rico S, Merenciano M, Dufour A-B, Vieira C, Quesneville H, Bergman CM, Andrieu O, Autard D, Nouaud D, Gonzalez J. 2019. Population-specific dynamics and selection pat- Ashburner M, Anxolabehere D. 2005. Combined evidence annota- terns of transposable element insertions in European natural pop- tion of transposable elements in genome sequences. PLoS Comp Biol. ulations. Mol Ecol. 28(6):1506–1522. 1(2):e22. Lewontin RC. 1974. The genetic basis of evolutionary change. New York: Rako L, Anderson AR, Sgro CM, Stocker AJ, Hoffmann AA. 2006. The Columbia University Press. association between inversion In(3R)Payne and clinally varying traits Li H, Stephan W. 2006. Inferring the demographic history and rate of in Drosophila melanogaster. Genetica 128(1–3):373–384. adaptive substitution in Drosophila. PLoS Genet. 2(10):e166. Rane RV, Rako L, Kapun M, Lee SF, Hoffmann AA. 2015. Genomic ev- Machado H, Bergland AO, Taylor R. 2019. Broad geographic sampling idence for role of inversion 3RP of Drosophila melanogaster in facil- reveals predictable, pervasive, and strong seasonal adaptation in itating climate change adaptation. Mol Ecol. 24(10):2423–2432. Drosophila. bioRxiv. doi: https://doi.org/10.1101/337543. Rech GE, Bogaerts-Marquez M, Barron MG, Merenciano M, Villanueva- Machado HE, Bergland AO, O’Brien KR, Behrman EL, Schmidt PS, Petrov Ca~nas JL, Horvath V, Fiston-Lavier A-S, Luyten I, Venkataram S, DA. 2016. Comparative population genomics of latitudinal variation Quesneville H, et al. 2019. Stress response, behavior, and develop- in Drosophila simulans and Drosophila melanogaster. Mol Ecol. ment are shaped by transposable element-induced mutations in 25(3):723–740. Drosophila. PLoS Genet. 15(2):e1007900. Macholan M, Baird SJ, Munclinger P, Dufkova P, Bımova B, Pialek J. 2008. Richardson MF, Weinert LA, Welch JJ, Linheiro RS, Magwire MM, Jiggins Genetic conflict outweighs heterogametic incompatibility in the FM, Bergman CM. 2012. Population genomics of the Wolbachia mouse hybrid zone? BMC Evol Biol. 8(1):271. endosymbiont in Drosophila melanogaster. PLoS Genet. Martino ME, Ma D, Leulier F. 2017. Microbial influence on Drosophila 8(12):e1003129. biology. Curr Opin Microbiol. 38:165–170. Rogers RL, Hartl DL. 2012. Chimeric genes as a source of rapid evolution Mateo L, Rech GE, Gonzalez J. 2018. Genome-wide patterns of local in Drosophila melanogaster. Mol Biol Evol. 29(2):517–529. adaptation in Western European Drosophila melanogaster natural Schlötterer C, Tobler R, Kofler R, Nolte V. 2014. Sequencing pools of populations. Sci Rep. 2018;8:16143. individuals – mining genome-wide polymorphism data without big McDonald JH, Kreitman M. 1991. Adaptive protein evolution at the Adh funding. Nat Rev Genet. 15(11):749–763. locus in Drosophila. Nature 351(6328):652–654. Schmidt PS, Paaby AB. 2008. Reproductive diapause and life-history Mettler LE, Voelker RA, Mukai T. 1977. Inversion clines in populations of clines in north American populations of Drosophila melanogaster. Drosophila melanogaster. Genetics 87(1):169–176. Evolution 62(5):1204–1215. Meyer F, Paarmann D, D’Souza M, Olson R, Glass EM, Kubal M, Paczian Schmidt PS, Zhu C-T, Das J, Batavia M, Yang L, Eanes WF. 2008. An T, Rodriguez A, Stevens R, Wilke A, et al. 2008. The metagenomics amino acid polymorphism in the couch potato gene forms the basis RAST server – a public resource for the automatic phylogenetic and for climatic adaptation in Drosophila melanogaster. Proc Natl Acad functional analysis of metagenomes. BMC Bioinformatics 9(1):386. Sci U S A. 105(42):16207–16211. Michalakis Y, Veuille M. 1996. Length variation of CAG/CAA trinucleo- Shi M, White VL, Schlub T, Eden J-S, Hoffmann AA, Holmes EC. 2018. No tide repeats in natural populations of Drosophila melanogaster and detectable effect of Wolbachia w Mel on the prevalence and abun- its relation to the recombination rate. Genetics 143:1713–1725. dance of the RNA virome of Drosophila melanogaster. Proc R Soc B. Nielsen R, Akey JM, Jakobsson M, Pritchard JK, Tishkoff S, Willerslev E. 285(1883):20181165. 2017. Tracing the peopling of the world through genomics. Nature Singh ND, Arndt PF, Clark AG, Aquadro CF. 2009. Strong evidence for 541(7637):302–310. lineage and sequence specificity of substitution rates and patterns in Novembre J, Johnson T, Bryc K, Kutalik Z, Boyko AR, Auton A, Indap A, Drosophila. Mol Biol Evol. 26(7):1591–1605. King KS, Bergmann S, Nelson MR, et al. 2008. Genes mirror geogra- Sprengelmeyer QD, Mansourian S, Lange JD, Matute DR, Cooper BS, Jirle phy within Europe. Nature 456(7218):98–101. EV, Stensmyr MC, Pool JE. 2020. Recurrent collection of Drosophila 2677 Kapun et al. . doi:10.1093/molbev/msaa120 MBE melanogaster from wild African environments and genomic insights Wang Y, Staubach F. 2018. Individual variation of natural into species history. Mol Biol Evol. 37(3):627–638. D. melanogaster-associated bacterial communities. FEMS Microbiol Staubach F, Baines JF, Künzel S, Bik EM, Petrov DA. 2013. Host species Lett. 365(6):fny017. and environmental effects on bacterial communities associated with Webster CL, Longdon B, Lewis SH, Obbard DJ. 2016. Twenty-five new Drosophila in the laboratory and in the natural environment. PLoS viruses associated with the Drosophilidae (Diptera). Evol Bioinform One 8(8):e70749. Online. 12:13–25. Szymura JM, Barton NH. 1986. Genetic analysis of a hybrid zone between Webster CL, Waldron FM, Robertson S, Crowson D, Ferrari G, Quintana the fire-bellied toads, Bombina bombina and B. variegata, near JF, Brouqui J-M, Bayne EH, Longdon B, Buck AH, et al. 2015. The Cracow in Southern Poland. Evolution 40(6):1141–1159. discovery, distribution, and evolution of viruses associated with Tauber E, Zordan M, Sandrelli F, Pegoraro M, Osterwalder N, Breda C, Drosophila melanogaster. PLoS Biol. 13(7):e1002210. Daga A, Selmin A, Monger K, Benna C, et al. 2007. Natural selection Werren JH, Baldo L, Clark ME. 2008. Wolbachia: master manipulators of favors a newly derived timeless allele in Drosophila melanogaster. invertebrate biology. Nat Rev Microbiol. 6(10):741–751. Science 316(5833):1895–1899. Whitlock MC, McCauley DE. 1999. Indirect measures of gene flow and Teixeira L, Ferreira A, Ashburner M. 2008. The Bacterial Symbiont migration: FST¼6 1/(4Nmþ1). Heredity 82(2):117–125. Wolbachia Induces Resistance to RNA Viral Infections in Wilfert L, Longdon B, Ferreira AGA, Bayer F, Jiggins FM. 2011. Drosophila melanogaster. PLoS Biol. 6(12):e2. Trypanosomatids are common and diverse parasites of Drosophila. Trinder M, Daisley BA, Dube JS, Reid G. 2017. Drosophila melanogaster as Parasitology 138(7):858–865. a high-throughput model for host-microbiota interactions. Front Wolff JN, Camus MF, Clancy DJ, Dowling DK. 2016. Complete mi- Microbiol. 8:751. tochondrial genome sequences of thirteen globally sourced Turner TL, Levine MT, Eckert ML, Begun DJ. 2008. Genomic analysis of strains of fruit fly (Drosophila melanogaster) form a powerful adaptive differentiation in Drosophila melanogaster. Genetics model for mitochondrial research. Mitochondrial DNA A. 179(1):455–473. 27(6):4672–4674. Umina PA, Weeks AR, Kearney MR, McKechnie SW, Hoffmann AA. Wright S. 1949. The genetical structure of populations. Ann Eugen. 2005. A rapid shift in a classic clinal pattern in Drosophila reflecting 15(1):323–354. climate change. Science 308(5722):691–693. Xiao F-X, Yotova V, Zietkiewicz E, Lovell A, Gehl D, Bourgeois S, Moreau Unckless RL. 2011. A DNA virus of Drosophila. PLoS One 6(10):e26564. C, Spanaki C, Plaitakis A, Moisan J-P, et al. 2004. Human X-chromo- Villemereuil P, Gaggiotti OE. 2015. A new FST-based method to uncover somal lineages in Europe reveal Middle Eastern and Asiatic contacts. local adaptation using environmental variables. Methods Ecol Evol. Eur J Hum Genet. 12(4):301–311. 6(11):1248–1258. Yukilevich R, True JR. 2008a. Incipient sexual isolation among cosmo- Walters AW, Hughes RC, Call TB, Walker CJ, Wilcox H, Peterson SC, politan Drosophila melanogaster populations. Evolution Rudman SM, Newell PD, Douglas AE, Schmidt PS, et al. 2020. The 62(8):2112–2121. microbiota influences the Drosophila melanogaster life history strat- Yukilevich R, True JR. 2008b. African morphology, behavior and pher- egy. Mol Ecol. 29:639–653. omones underlie incipient sexual isolation between us and Wang Y, Kapun M, Waidele L, Kuenzel S, Bergland AO, Staubach F. 2020. Caribbean Drosophila melanogaster. Evolution 62(11):2807–2828. Common structuring principles of the Drosophila melanogaster Zanini F, Brodin J, Thebo L, Lanz C, Bratt G, Albert J, Neher RA. 2015. microbiome on a continental scale and between host and substrate. Population genomics of intrapatient HIV-1 evolution. eLife Environ Microbiol Rep. 12(2):220–228. 4:e11282. 2678