Drosophila Evolution over Space and Time (DEST): A New Population Genomics Resource Martin Kapun ,*,†,‡,1,2 Joaquin C.B. Nunez,†,3 Marıa Bogaerts-Marquez,†,‡,4 Jes us Murga-Moreno,†,‡,5,6 †,‡,7 Margot Paris, Joseph Outten, Marta Coronado-Zamora, Courtney Tern, Omar Rota-Stabelli,‡,8 3 ‡,4 3 Maria P. Garcıa Guerreiro,‡,5 S onia Casillas ,‡,5,6 Dorcas J. Orengo,‡,9,10 Eva Puerma,‡,9,10 Maaria Kankare, Lino Ometto ,‡,12 Volker Loeschcke,‡,13 Banu S. Onder ,‡,14 Jessica K. Abbott,‡,15 ‡,11 Stephen W. Schaeffer ,§,16 Subhash Rajpurohit ,§,17,18 Emily L. Behrman,§,17,19 Mads F. Schou ,‡,13,15 Thomas J.S. Merritt,§,20 Brian P. Lazzaro,§,21 Amanda Glaser-Schmitt ,‡,22 Eliza Argyridou,‡,22 Fabian Staubach ,‡,23 Yun Wang,‡,23 Eran Tauber,‡,24 Svitlana V. Serga ,‡,25,26 Daniel K. Fabian,§,27 Kelly A. Dyer,§,28 Christopher W. Wheat,‡,29 John Parsch,‡,22 Sonja Grath ,‡,22 Marija Savic Veselinovic,‡,30 Marina Stamenkovic-Radak,‡,30 Mihailo Jelic,‡,30 Antonio J. Buendıa-Ruız,‡,31 Maria Josefa G omez-Julian,‡,31 Maria Luisa Espinosa-Jimenez, Francisco D. Gallardo-Jimenez, Aleksandra Patenkovic ,‡,33 ‡,31 ‡,32 Katarina Eric,‡,33 Marija Tanaskovic,‡,33 Anna Ullastres,‡,4 Lain Guio,‡,4 Miriam Merenciano,‡,4 Sara Guirao-Rico,‡,4 Vivien Horvath,‡,4 Darren J. Obbard ,‡,34 Elena Pasyukova,‡,35 Vladimir E. Alatortsev,‡,35 Cristina P. Vieira,‡,36,37 Jorge Vieira,‡,36,37 Jorge Roberto Torres,‡,38 Iryna Kozeretska,‡,25,26 Oleksandr M. Maistrenko,‡,25,39 Catherine Montchamp-Moreau,‡,40 Dmitry V. Mukha,‡,41 Heather E. Machado,§,42,43 Keric Lamb,3 T^ania Paulo,‡,44 Leeban Yusuf,‡,45 Antonio Barbadilla ,‡,5,6 Dmitri Petrov,*,‡,42 Paul Schmidt,*,‡,§,16 Josefa Gonzalez,*,‡,§,4 Thomas Flatt ,*,‡,§,7 and Alan O. Bergland,*,‡,§,3 1 Department of Evolutionary Biology and Environmental Studies, University of Zürich, Switzerland 2 Department of Cell & Developmental Biology, Center of Anatomy and Cell Biology, Medical University of Vienna, Vienna, Austria 3 Department of Biology, University of Virginia, Charlottesville, VA, USA 4 Institute of Evolutionary Biology, CSIC-Universitat Pompeu Fabra, Barcelona, Spain 5 Department of Genetics and Microbiology, Universitat Autonoma de Barcelona, Barcelona, Spain 6 Institute of Biotechnology and Biomedicine, Universitat Autonoma de Barcelona, Barcelona, Spain 7 Department of Biology, University of Fribourg, Fribourg, Switzerland 8 Center Agriculture Food Environment, University of Trento, San Michele all’ Adige, Italy 9 Departament de Genètica, Microbiologia i Estadıstica, Facultat de Biologia, Universitat de Barcelona, Barcelona, Spain 10 Article Institut de Recerca de la Biodiversitat (IRBio), Universitat de Barcelona, Barcelona, Spain 11 Department of Biological and Environmental Science, University of Jyv€askyl€a, Jyv€askyl€a, Finland 12 Department of Biology and Biotechnology, University of Pavia, Pavia, Italy 13 Department of Biology, Aarhus University, Aarhus, Denmark 14 Department of Biology, Hacettepe University, Ankara, Turkey 15 Department of Biology, Lund University, Lund, Sweden 16 Department of Biology, The Pennsylvania State University, University Park, PA, USA 17 Department of Biology, University of Pennsylvania, Philadelphia, PA, USA 18 Division of Biological and Life Sciences, School of Arts and Sciences, Ahmedabad University, Ahmedabad, India 19 Janelia Research Campus, Ashburn, VA, USA 20 Department of Chemistry & Biochemistry, Laurentian University, Sudbury, ON, Canada 21 Department of Entomology, Cornell University, Ithaca, NY, USA 22 Division of Evolutionary Biology, Faculty of Biology, Ludwig-Maximilians-Universit€at, Munich, Germany 23 Department of Evolution and Ecology, University of Freiburg, Freiburg, Germany 24 Department of Evolutionary and Environmental Biology, Institute of Evolution, University of Haifa, Haifa, Israel 25 Department of General and Medical Genetics, Taras Shevchenko National University of Kyiv, Kyiv, Ukraine 26 State Institution National Antarctic Scientific Center, Ministry of Education and Science of Ukraine, Kyiv, Ukraine 27 Department of Genetics, University of Cambridge, Cambridge, United Kingdom ß The Author(s) 2021. 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-NonCommercial License (https://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 5782 Mol. Biol. Evol. 38(12):5782–5805 doi:10.1093/molbev/msab259 Advance Access publication September 1, 2021 DEST Population Genomics . doi:10.1093/molbev/msab259 MBE 28 Department of Genetics, University of Georgia, Athens, GA, USA 29 Department of Zoology, Stockholm University, Stockholm, Sweden 30 Faculty of Biology, University of Belgrade, Belgrade, Serbia 31 IES Eladio Caba~ nero, Tomelloso, Spain 32 IES Jose de Mora, Baza, Spain 33 Institute for Biological Research “Sinisa Stankovic”, National Institute of Republic of Serbia, University of Belgrade, Belgrade, Serbia 34 Institute of Evolutionary Biology, University of Edinburgh, Edinburgh, United Kingdom 35 Institute of Molecular Genetics of the National Research Centre “Kurchatov Institute”, Moscow, Russia 36 Instituto de Biologia Molecular e Celular (IBMC), Porto, Portugal 37 Instituto de Investigaç~ao e Inovaç~ao em Saude, Universidade do Porto, Porto, Portugal 38 La ciència al teu mon, Barcelona, Spain 39 Structural and Computational Biology Unit, European Molecular Biology Laboratory, Heidelberg, Germany 40 UMR Evolution, Genomes, Comportement et Ecologie, Universite Paris-Saclay, CNRS, Gif-sur-Yvette, France 41 Vavilov Institute of General Genetics, Russian Academy of Sciences, Moscow, Russia 42 Department of Biology, Stanford University, Stanford, CA, USA 43 Wellcome Trust Sanger Institute, Hinxton, United Kingdom 44 Departamento de Biologia Animal, Instituto Gulbenkian de Ci^encia, Oeiras, Portugal 45 Center for Biological Diversity, University of St. Andrews, St Andrews, United Kingdom † These authors contributed equally to this work. ‡ The European Drosophila Population Genomics Consortium (DrosEU). § The Drosophila Real-Time Evolution Consortium (DrosRTEC). *Corresponding authors: E-mails: martin.kapun@uzh.ch; aob2x@virginia.edu; thomas.flatt@unifr.ch; josefa.gonzalez@csic.es; dpetrov@stanford.edu; schmidtp@upenn.edu. Associate editor: Rasmus Nielsen Abstract Drosophila melanogaster is a leading model in population genetics and genomics, and a growing number of whole- genome data sets from natural populations of this species have been published over the last years. A major challenge is the integration of disparate data sets, often generated using different sequencing technologies and bioinformatic pipe- lines, which hampers our ability to address questions about the evolution of this species. Here we address these issues by developing a bioinformatics pipeline that maps pooled sequencing (Pool-Seq) reads from D. melanogaster to a hologe- nome consisting of fly and symbiont genomes and estimates allele frequencies using either a heuristic (PoolSNP) or a probabilistic variant caller (SNAPE-pooled). We use this pipeline to generate the largest data repository of genomic data available for D. melanogaster to date, encompassing 271 previously published and unpublished population samples from over 100 locations in >20 countries on four continents. Several of these locations have been sampled at different seasons across multiple years. This data set, which we call Drosophila Evolution over Space and Time (DEST), is coupled with sampling and environmental metadata. A web-based genome browser and web portal provide easy access to the SNP data set. We further provide guidelines on how to use Pool-Seq data for model-based demographic inference. Our aim is to provide this scalable platform as a community resource which can be easily extended via future efforts for an even more extensive cosmopolitan data set. Our resource will enable population geneticists to analyze spatiotemporal genetic patterns and evolutionary dynamics of D. melanogaster populations in unprecedented detail. Key words: Drosophila melanogaster, population genomics, SNPs, evolution, adaptation, demography. Introduction colonization history and demography (Caracristi and The vinegar fly Drosophila melanogaster is one of the oldest Schlötterer 2003; Li and Stephan 2006; Duchen et al. 2013; and most important genetic model systems and has played a Grenier et al. 2015; Bergland et al. 2016; Arguello et al. 2019; key role in the development of theoretical and empirical pop- Kapopoulou et al. 2020), and spatiotemporal diversification ulation genetics (e.g., Schneider 2000; Hales et al. 2015; Haudry patterns of this species and its close relatives (Kolaczkowski et al. 2020). Through decades of work, we now have a basic et al. 2011; Fabian et al. 2012; Bergland et al. 2014; Kapun et al. picture of the evolutionary origin (David and Capy 1988; 2016, 2020; Lack et al. 2016; Machado et al. 2016, 2021). Lachaise et al. 1988; Keller 2007; Sprengelmeyer et al. 2020), The availability of high-quality reference genomes (Adams 5783 Kapun et al. . doi:10.1093/molbev/msab259 MBE 2000; Celniker and Rubin 2003; dos Santos et al. 2015) based on PoolSNP (Kapun et al. 2020) and SNAPE-pooled and genetic tools (Schneider 2000; Duffy 2002; Jennings (Raineri et al. 2012), respectively. We also provide tools for 2011; Hales et al. 2015; Haudry et al. 2020) facilitates placing performing in silico pooling of existing inbred (haploid) lines evolutionary studies of flies in a mechanistic context, allow- that exist as part of other Drosophila population genomic ing for the functional characterization of ecologically rele- resources (Langley et al. 2012; Pool et al. 2012; Grenier et al. vant polymorphisms (e.g., de Jong and Bochdanovits 2003; 2015; Kao et al. 2015; Lack et al. 2015, 2016). This pipeline is Paaby et al. 2010, 2014; Mateo et al. 2014; Kapun et al. 2016; also designed to be flexible, facilitating the streamlined addi- Durmaz et al. 2018, 2019; Ramaekers et al. 2019). tion of new population samples as they arise. Recently, work on the evolutionary biology of Drosophila Using this pipeline, we generated a unified data set of has been fueled by a growing number of population genomic pooled allele frequency estimates of D. melanogaster sampled data sets from field collections across a large portion of D. across a large portion of its world-wide distribution, including melanogaster’s range (Kapun et al. 2020; Grenier et al. 2015; Europe, North America, Africa, Australia, and Asia. This data Arguello et al. 2019; Guirao-Rico and Gonzalez 2019; Machado set is the result of the collaborative efforts of the European et al. 2021). These genomic data consist either of re- DrosEU (Kapun et al. 2020) and DrosRTEC (Machado et al. sequenced inbred (or haploid) individuals (e.g., Langley et al. 2021) consortia and combines both novel and previously 2012; Mackay et al. 2012; Grenier et al. 2015; Lack et al. 2015, published population genomic data. Our data set combines 2016; Mateo et al. 2018; Kapopoulou et al. 2020) or pooled samples from 100 localities, 55 of which were sampled at two sequencing (Pool-Seq) of outbred population samples (Pool- or more time points across the reproductive season (10–15 Seq; e.g., Kolaczkowski et al. 2011; Fabian et al. 2012; Bastide et generations/year) for one or more years. Collectively, these al. 2013; Campo et al. 2013; Bergland et al. 2014; Machado et samples represent >13,000 individuals, cumulatively se- al. 2016, 2021; Kapun et al. 2016, 2020). Pooled resequencing quenced to >16,000 coverage or 1 per fly. The cost provides accurate and precise estimates of allele frequencies effectiveness of Pool-Seq has enabled us to estimate across most of the allele frequency spectrum (Zhu et al. 2012; genome-wide allele frequencies over geographic space (con- Lynch et al. 2014; Schlötterer et al. 2014) at a fraction of the tinental and subcontinental) and time (seasonal, annual, and cost of individual-based sequencing. Although Pool-Seq decadal) scales, thus making our data a unique resource for retains limited information about linkage disequilibrium advancing our understanding of fundamental adaptive and (LD) relative to individual sequencing (Feder et al. 2012), neutral evolutionary processes. We provide data in two file Pool-Seq data can be used to infer complex demographic formats (VCF and GDS: Danecek et al. 2011; Zheng et al. histories (e.g., Cheng et al. 2012; Bergland et al. 2016; Deitz 2017), thus allowing researchers to utilize a variety of tools et al. 2016; Corbett-Detig and Nielsen 2017; Gould et al. 2017; for computational analyses. Our data set also contains sam- Giesen et al. 2020), characterize levels of diversity (Kofler, pling and environmental metadata to enable various down- Orozco-terWengel, et al. 2011; Kofler, Pandey, et al. 2011), stream analyses of biological interest. We further employed and infer genomic loci involved in recent adaptation demographic modeling to investigate the evolutionary history in nature (Flatt 2016; Kapun et al. 2016, 2020; Gould et al. of two distinct genetic clusters in Europe using the Drosophila 2017; Bogaerts-Marquez et al. 2021; Machado et al. 2021) Evolution over Space and Time (DEST) Pool-Seq data set and and during experimental evolution (e.g., Turner et al. 2011; developed guidelines for using Pool-Seq data for model-based Burke 2012; Orozco-terWengel et al. 2012; Kofler and demographic inference using the python package moments. Schlötterer 2014). However, the rapidly increasing number of genomic data sets processed with different bioinformatic Results pipelines makes it difficult to compare results across studies and to jointly analyze multiple data sets. Differences among Integrating a Worldwide Collection of D. melanogaster bioinformatic pipelines include filtering methods for the raw Population Genomics Resources reads, mapping algorithms, the choice of the reference ge- We developed a modular and standardized pipeline for nome, or SNP calling approaches, potentially generating biases generating allele frequency estimates from pooled resequenc- when combining processed data sets from different sources ing of D. melanogaster genomes (supplementary fig. S1, for joint analyses (e.g., Gautier et al. 2013; Hoban et al. 2016). Supplementary Material online). Using this pipeline, we assem- To address these issues, we have developed a modular bled a data set of allele frequencies from 271 D. melanogaster bioinformatics pipeline to map Pool-Seq reads to a hologe- populations sampled around the world (fig. 1A and supple- nome consisting of fly and microbial genomes, to remove mentary table S1, Supplementary Material online). Many of reads from potential Drosophila simulans contaminants, these samples were collected at the same location, at different and to estimate allele frequencies using two complementary seasons and over multiple years (fig. 1B). The nature of the SNP callers. Our pipeline is available as a Docker image (avail- genomic data for each population varies as a consequence of able from https://dest.bio, last accessed September 6, 2021) to biological origin (e.g., inbred lines or Pool-Seq), library prepa- standardize versions of software used for filtering and map- ration method, and sequencing platform. ping, to make the pipeline available independently of the To assess whether these features affect basic attributes of operating system used, and to facilitate future updates and the data set, we calculated six basic quality metrics focusing modification of the pipeline. In addition, our pipeline allows on the Pool-Seq samples (fig. 1C and supplementary table S2, using either heuristic or probabilistic methods for SNP calling, Supplementary Material online). On average, median read 5784 DEST Population Genomics . doi:10.1093/molbev/msab259 MBE A C DrosEU DrosRTEC Read Depth 200 150 100 50 50 DGN Latitude 0.20 % Missing DrosEU 0 0.15 DrosRTEC 0.10 0.05 0 −5 150 Num. Flies 100 00 0 0 0 10 20 −1 50 Longitude 90 B Eff. RD 2009 2010 2011 2012 2013 2014 2015 2016 60 30 60 0.8 PCR Dup. % 0.6 Latitude 50 0.4 Europe 0.2 North America 0.0 40 0.6 Prop. Sim. 0.4 30 0.2 0.0 Ol Ja t n Ol Ja t n Ol Ja t n Ol Ja t n Ol Ja t n Ol Ja t n Ol Ja t n Ol Ja t n Population Ju Ju Ju Ju Ju Ju Ju Ju c c c c c c c c Collection Date Europe North America FIG. 1. Sampling location, dates, and quality metrics. (A) Map showing the 271 sampling localities forming the DEST data set. Colors denote the data sets of origin (DGN, DrosEU, or DrosRTEC). (B) Collection dates for localities sampled more than once. (C) General sample features of the DEST data set. The x-axis represents the population sample, ordered by the average read depth. depth across samples is 62x (range: 10–217x). The per- generally absent (<1% D. simulans specific reads; supplemen- nucleotide missing allele frequency rate was less than 7% tary table S1, Supplementary Material online). for most (95%) of the samples. Excluding populations with high missing data rate (>7%), the proportion of sites with Identification and Quality Control of SNP missing data was positively correlated with read depth Polymorphisms (P ¼ 1.2  109, R2 ¼ 0.4). The positive correlation between In order to determine appropriate SNP calling and filtering read depth and missing data rate is primarily due to an in- parameters, and to identify potentially problematic popula- creased sensitivity to identify indels. The number of flies per tion samples, we first calculated the ratio of the number of sample varied from 33 to 205, with considerable heterogene- nonsynonymous polymorphisms to the number of synony- ity among the DrosRTEC samples [standard deviation (SD) ¼ mous polymorphisms (pN/pS) for each population sample 30], but not among DrosEU samples (SD ¼ 0.04). Variation in across the whole genome. Because nonsynonymous changes the number of flies and in sequencing depth is reflected in the are expected to be under strong purifying selection (Kreitman effective coverage (NEff) of each pool, an estimate of the num- 1983), the pN/pS metric can reflect the presence of sequencing ber of independent reads after accounting for double bino- errors that would disproportionately inflate pN relative to pS. mial sampling that occurs during Pool-Seq (Kolaczkowski et Our primary goal was not to provide novel estimates of pN/pS al. 2011; Feder et al. 2012; fig. 1C). There was considerable but rather to ensure that all population samples have esti- variation in PCR duplicate rate among samples, with notable mates that are consistent with estimates generated from in- differences between batches of DrosEU samples (6% in dependent Drosophila data sets (Mackay et al. 2012). 2014 vs. 18% in 2015/16; t-test, P ¼ 1.8  1019) and For the PoolSNP data set, we varied the global minor allele DrosRTEC samples (3% in samples collected as part of count (MAC) and global minor allele frequency (MAF) and Bergland et al. 2014 vs. 14% in samples collected as part then calculated pN/pS. MAC thresholds <50 resulted in large of Machado et al. 2021; P ¼ 6.37  103). Curiously, the variances of pN/pS caused by 20 outlier populations charac- 2015/2016 DrosEU samples were made with a PCR-free kit, terized by unusually high pN/pS ratios and numbers of private suggesting that the observed PCR duplicates were optical SNPs (supplementary table S3, Supplementary Material on- duplicates and not amplification artifacts. Contamination of line, and fig. 2A and B) indicating that there may be elevated samples by D. simulans varied among populations but was numbers of sequencing errors in some samples. Some 5785 Kapun et al. . doi:10.1093/molbev/msab259 MBE A B PoolSNP 2.50 6 2.00 pN pS 1.50 1.00 Log10(No. of private SNPs) 5 0.50 0.00 0 25 50 75 100 4 Minor allele count (MAC) C SNAPE 2.5 2.0 3 pN pS 1.5 1.0 0.5 2 0.0 0.5 1.0 1.5 2.0 0.00 0.05 0.10 0.15 0.20 pN pS Minor allele frequency (MAF) Status Exclude Keep FIG. 2. Quality control of SNPs called with SNAPE-pooled and PoolSNP. Panel (A) shows genome-wide pN/pS ratios and the log10-scaled number of private SNPs for all Pool-Seq samples based on SNP calling with SNAPE-pooled. We highlight 20 outlier samples in red, which are characterized by exceptionally high values of both metrics. The dashed black lines indicate the 95% confidence limits (average þ 1.96 SD) for both statistics. The vertical green dashed line highlights the empirical estimate of pN/pS calculated from individual sequencing data of the DGRP freeze2 data set (Mackay et al. 2012). The green diamond shows the corresponding value of the DGRP population, which was pool-sequenced as part of the DrosRTEC data set (NC_ra_03_n; Zhu et al. 2012). Panels (B) and (C) show the effects of heuristic MAC and MAF thresholds on pN/pS ratios in SNP data based on PoolSNP and SNAPE-pooled, respectively. Blue lines in both panels show average genome-wide pN/pS ratios across 271 and 246 populations, respectively. The blue ribbons depict the corresponding standard deviations. The 20 outlier samples, which are characterized in panel (A), are highlighted red. In addition, pN/pS ratios of the DGRP Pool-Seq sample (NC_ra_03_n) are shown at different cut-offs as green diamonds and the empirical values from the DGRP freeze2 data set are indicated as dashed green lines. (n ¼ 17) of these samples had previously been found to show ratios. Accordingly, we excluded these 20 samples from fur- positive values of Tajima’s D across the whole genome ther analyses of low-frequency variants and private SNPs and (Kapun et al. 2020). We observed that, as expected, pN/pS applied a conservative local MAF filter of 5% for the remain- was negatively correlated with MAC (linear regression; der of the SNAPE-pooled analysis to avoid misclassification of P < 0.001; fig. 2B) and that applying a MAC threshold of 50 sequencing errors as low-frequency variants. Our SNAPE- reduced the elevated pN/pS ratios of the 20 aforementioned pooled results identified 8,541,651 polymorphisms segregat- outlier samples to values similar to the rest of the data set, ing among the remaining 226 samples. Below, we discuss the suggesting that potential sequencing errors had been largely geographic distribution and global frequency of SNPs identi- removed. To minimize false-positive variant calling, we chose fied using these two methods in order to provide insight into MAC ¼ 50 and MAF ¼ 0.001 as conservative threshold the marked discrepancy in the number of SNPs that they parameters for SNP calling with PoolSNP. Using these param- identify. eters, PoolSNP identified 4,381,144 polymorphisms segregat- ing among the 271 D. melanogaster samples (Pool-Seq plus Similarity of SNP Polymorphisms Detected with DGN), and 4,042,456 polymorphisms segregating among the PoolSNP and SNAPE-Pooled 246 Pool-Seq samples (excluding DGN). We calculated three metrics related to the amount of poly- In contrast to PoolSNP, SNAPE-pooled calls variants in morphism discovered by our pipelines: the abundance of each sample separately using a probabilistic approach which polymorphisms segregating in n populations across each integrates allelic information across all populations for heu- chromosome (fig. 3A), the difference of discovered polymor- ristic SNP calling. To quantify the number of putative se- phisms between SNAPE-pooled and PoolSNP (defined as the quencing errors among low frequency variants we varied absolute value of PoolSNP minus SNAPE-pooled; fig. 3B), and the local MAF threshold per sample and calculated pN/pS the amount of polymorphism discovered per MAF bin (fig. for each sample in the SNAPE-pooled data set. Similar to 3C). We evaluated these three metrics across a 2  2 filtering PoolSNP, we found that elevated pN/pS was negatively corre- scheme: two MAF filters (0.001, 0.05) and two sample sets lated with a local MAF threshold (linear regression; P < 0.001; (the whole data set of 246 samples; and the 226 samples that fig. 2C) and that the 20 aforementioned problematic samples passed the sequencing error filter in SNAPE-pooled; see also had a strong effect on the variance and mean of pN/pS Identification and Quality Control). Notably, PoolSNP was 5786 DEST Population Genomics . doi:10.1093/molbev/msab259 MBE FIG. 3. Polymorphism data in the PoolSNP and SNAPE data sets. (A) Number of polymorphic sites discovered across populations. The x-axis shows the number of populations that share a polymorphic site. The y-axis corresponds to the number of polymorphic sites shared by any number of populations, on a log10 scale. The colored lines represent different chromosomes and are stacked on top of each other. (B) The difference of discovered polymorphisms between SNAPE-pooled and PoolSNP. (C) Number of polymorphic sites as a function of allele frequency and the number of populations in which the polymorphisms are present. The color gradient represents the number of variant alleles from low to high (black to green). The x-axis is the same as in (A), and the y-axis is the MAF. The 2  2 filtering scheme is shown on the right side of the figure. biased toward identification of common SNPs present in Mutation-Class Frequencies multiple samples, whereas SNAPE-pooled was more sensitive We estimated the percentage of mutation classes (e.g., A ! to the identification of polymorphisms that appeared in few C, A ! G, A ! T, etc.) accepted as polymorphisms in both populations only (fig. 3B). For example, at a MAF filter of our SNP calling pipelines and classified these loci as being 0.001, SNAPE-pooled discovered more polymorphisms that either “rare” (i.e., AF <5% and shared in less than 50 popu- were shared in less than 25 populations (relative to PoolSNP), lations) or “common” (AF >5% and shared in more than 150 and these accounted for 79% of all polymorphisms discov- populations). For this analysis, we classified the minor allele as ered by the pipeline. Likewise, at a MAF filter of 0.05, SNAPE- the derived allele. Figure 4A shows the percentage of each pooled discovered more polymorphisms that were shared in mutation class for the 226 populations which passed filters in less than 97 populations; these accounted for 71% of all both SNAPE-pooled and PoolSNP. In addition, we overlaid, as discovered polymorphisms. SNAPE-pooled identifies fewer a horizontal line, the expected mutation frequencies for rare polymorphic sites that are shared among a large number of (blue; Assaf et al. 2017) and common (red; Mackay et al. 2012) populations than PoolSNP does because SNAPE-pooled does mutations. In general, our SNP discovery pipelines produced not integrate information across multiple populations. mutation-class relative frequencies of rare and common Consequently, SNAPE-pooled can fail to identify SNPs that mutations that are consistent with empirical expectations, are at low overall frequencies and get called as monomorphic however, there were some exceptions to this pattern. For or missing in a subset of populations given the posterior example, the frequencies of the C/G rare mutation-class were consistently underestimated by both callers, a phenom- probability thresholds that we employed (see Materials and enon that might be related to the known GC bias of modern Methods). sequencing machines (Benjamini and Speed 2012). The cor- We also compared AF estimates between the two callers relation between SNP calling pipelines was high across both using the data set of 226 populations applying a local MAF common and rare mutation classes, with marginal discrep- filter of 0.05 in the SNAPE-pooled data set (see supplementary ancies observed for rare variants (fig. 4B). table S2, Supplementary Material online). Among the posi- tions identified as polymorphic by both calling methods, our Inversion Frequencies frequency estimates were identical for the great majority of Using a set of inversion-specific marker SNPs (Kapun et al. SNPs (92–99.67%) in all samples analyzed. Between 0.1% 2014), we estimated the frequencies of seven cosmopolitan and 7.1% of the polymorphic SNPs differed by less than inversion polymorphisms (In(2L)t, In(2R)NS, In(3L)P, In(3R)C, 5% frequency between the two methods, 0.003–2.1% of poly- In(3R)K, In(3R)Mo, and In(3R)Payne). We found that most of morphic SNPs differed by 5–10% frequency and only up the 271 populations were polymorphic for at least one or to 0.3% varied >10% frequency (supplementary table S4, more chromosomal inversions (supplementary table S1, Supplementary Material online). Finally, on average 13.32% Supplementary Material online). Although most inversions of the positions analyzed were called as polymorphic by were either absent or rare (average frequencies: In(2R)NS ¼ PoolSNP whereas there were monomorphic or no data 5.2% [6 4.7% SD], In(3L)P ¼ 3.1% [6 4.3% SD], according to SNAPE-pooled, consistent with the use of a In(3R)C ¼ 2.5% [6 2.3% SD], In(3R)K ¼ 1.8% [6 7.4% SD], hard threshold of the posterior-probability in the SNAPE call- In(3R)Mo ¼ 2.2% [6 3.6% SD] and In(3R)Payne ¼ 5.7% [6 ing step (supplementary table S4, Supplementary Material 7.1% SD]), only In(2L)t segregated at substantial frequencies in online). most populations (average frequency ¼ 18.3% [6 11% SD]). 5787 Kapun et al. . doi:10.1093/molbev/msab259 MBE A B AC AG AT CA CG CT 0.16 Common Corr. = 0.998, P << 0.01 0.20 Rare Corr. = 0.951, P << 0.01 GA 0.15 CT AG SNAPE Mutation Frequency CT 0.10 TC 0.12 Mutation Frequency 0.05 CA TA 0.00 GT TA 0.08 AT AT GA GC GT TA TC TG 0.20 CA AG AC TC 0.15 GC GT GC CG TG 0.10 0.04 AC CG 0.05 TG 0.00 0.04 0.08 0.12 0.16 rare PoolSNP common PoolSNP PoolSNP Mutation Frequency rare SNAPE common SNAPE Rare Common FIG. 4. Frequencies of observed nucleotide polymorphism in the DEST data set (226 populations common to PoolSNP and SNAPE-pooled). (A) Each panel represents a mutation type. The red color indicates common mutations (AF >0.05, and common in more than 150 populations) whereas the blue color indicates rare mutations (AF <0.05, and shared in less than 50 populations). The dark colors correspond to the PoolSNP pipeline and the soft colors correspond to the SNAPE-pooled pipeline. The hovering red and blue horizontal lines represent the estimated mutation rates for common and rare mutations, respectively. (B) Correlation between the observed mutation frequencies seen in SNAPE-pooled and PoolSNP. The one-to-one correspondence line is shown as a black-dashed diagonal. Correlation estimates (Pearson’s correlation) and P values for common and rare mutations are shown. We found that our novel inversion frequency estimates of the that has passed quality control, and NEff (Kofler, Orozco- DrosEU data from 2014 were highly consistent with previous terWengel, et al. 2011; Kolaczkowski et al. 2011; Feder et al. estimates from Kapun et al. (2020) as coefficients of determi- 2012; Schlötterer et al. 2014). Similar to AF estimates, the nation (R2) ranged from 91% to 99%. Pearson correlation coefficients for both coverage and effec- tive coverage were large (0.92, 0.95, 0.90 for Machado et al. Comparison to Previously Published Data Sets [2021], Kapun et al. [2020], and Bergland et al. [2014], respec- We compared the allele frequency and read depth estimates tively; see supplementary figs. S7–S12, Supplementary from the DEST data set (based on PoolSNP) to previously Material online), indicating that sample identity was pre- published estimates by Bergland et al. (2014), and Kapun et al. served appropriately. However, the concordance correlation (2020), Machado et al. (2021). For these data sets, we coefficients were substantially lower between the data sets employed two types of correlations: the nominal correlation (0.24, 0.88, 0.79, respectively), indicating systematic differences (i.e., Pearson’s correlation; CO) and the concordance correla- in read depth between the DEST data set and previously tion coefficient (CCC; Lin 1989; Liao and Lewis 2000). The published data. Indeed, read depth estimates were on average CCC determines how much the observed data deviate from 12%, 14%, and 20% lower in the DEST data set as com- the line of perfect concordance (i.e., the 45 degree-line on a pared with the previously published data in Machado et al. square scatter plot). (2021), Kapun et al. (2020), and Bergland et al. (2014), respec- Estimates of allele frequency were strongly correlated and tively. The lower read depth and effective read depth esti- consistent with previously published data. The strongest cor- mates in the DEST data set reflect our more stringent quality relation of DEST AF and previously published AF was ob- control and filtering. served with the data of Kapun et al. (2020) (average CO and CCC >0.99; fig. 5, top row and supplementary fig. S4, Genetic Diversity Supplementary Material online). AF correlations with We estimated nucleotide diversity (p), Watterson’s h, and Machado et al. (2021) are also generally high (average CO Tajima’s D for both the PoolSNP and SNAPE-pooled data and CCC >0.98; fig. 5, top row and supplementary fig. S5, sets (supplementary table S5, Supplementary Material on- Supplementary Material online). AF correlations with the line). Results for the African, European, and North data from Bergland et al. (2014) were lower (0.94; supplemen- American population samples are presented in figure 6 tary fig. S6, Supplementary Material online), likely reflecting (also see supplementary fig. S13, Supplementary Material on- differences in data processing and quality control. line for estimates by chromosome arm). All estimates were We also examined two aspects of read depth, that is, nom- positively correlated between PoolSNP and SNAPE-pooled inal coverage (COV), the number of reads mapping to a site (P < 0.001), with Pearson’s correlation coefficients of 0.90, 5788 DEST Population Genomics . doi:10.1093/molbev/msab259 MBE Machado et al. 2019 Kapun et al. 2020 Bergland et al. 2014 200 20 75 150 15 50 100 10 AF 25 50 5 0 0 0 0.7 0.8 0.9 1.0 0.80 0.85 0.90 0.95 1.00 0.6 0.7 0.8 0.9 15 6 Density 15 10 4 10 5 COV 5 2 0 0 0 0.0 0.5 1.0 0.00 0.25 0.50 0.75 1.00 0.0 0.5 1.0 15 20 5 15 4 10 3 10 2 NEFF 5 5 1 0 0 0 0.0 0.5 1.0 0.00 0.25 0.50 0.75 1.00 0.0 0.5 1.0 Value CCC CO type Actual Random FIG. 5. Correlations between DEST data set and previously published data sets. Correlations between allele frequencies (AF), Nominal Coverage (COV), and Effective Coverage (NEFF) between the DEST data set (using the PoolSNP method) and the three previous Drosophila data sets: Machado et al. (2021), Kapun et al. (2020), and Bergland et al. (2014). For each data set, we show the distribution of two types of correlation coefficients: the nominal (Pearson’s) correlation (CO; dashed lines) and the concordant correlation (CCC; solid lines). In addition to the actual correlations between the data sets (red distributions), we show the distributions of correlations estimated with random population pairs (green distributions). A B PoolSNP SNAPE PoolSNP SNAPE Watterson's θ 0.007 0.006 0.005 Nucleotide diversity (π) 0.004 0.006 0.003 Africa North Europe North Europe America America C 0.005 PoolSNP SNAPE 1.500 Tajima's D 1.000 0.500 0.004 0.000 Africa North Europe North Europe Africa North Europe North Europe America America America America FIG. 6. Population genetic estimates for African, European, and North American populations. Shown are genome-wide estimates of (A) nucleotide diversity (p), (B) Watterson’s h and (C) Tajima’s D for African populations using the PoolSNP data set, and for European and North American populations using both the PoolSNP and SNAPE-pooled (SNAPE) data sets. As can be seen from the figure, estimates based on PoolSNP versus SNAPE-pooled (SNAPE) are highly correlated (see main text). Genetic variability is seen to be highest for African populations, followed by North American and then European populations, as previously observed (e.g., see Lack et al. [2016] and Kapun et al. [2020]). 0.83, and 0.70 for p, Watterson’s h, and Tajima’s D, respec- in North American populations (linear models, all P > 0.05), tively. Higher values of genetic diversity were obtained for the suggesting that data from populations with heterogeneous SNAPE-pooled data set, probably due to its higher sensitivity pool sizes can be safely merged for accurate population ge- for detecting rare variants (see Patterns of Polymorphism nomic analysis. between PoolSNP and SNAPE-Pooled). Pool size had no sig- The highest levels of genetic diversity were observed for nificant effect on the four summary statistics in European or ancestral African populations (mean p ¼ 0.0060, mean 5789 Kapun et al. . doi:10.1093/molbev/msab259 MBE h ¼ 0.0059); North American populations exhibited higher significantly correlated with data type (PC1: P ¼ 0.352, PC2: genetic variability (mean p ¼ 0.0054, mean h ¼ 0.0054) than P ¼ 0.223, PC3: P ¼ 0.998). European populations (mean p ¼ 0.0049, mean h ¼ 0.0048). These results are consistent with previous observations based Geographic Proximity Analysis on individual genome sequencing (e.g., see Lack et al. [2016] The geographic distribution of our samples allows leveraging and Kapun et al. [2020]). Our observations are also consistent basic principles of phylogeography and population genetics with previous estimates based on pooled data from three to assess the biological significance of rare SNPs (Wright 1943; North American populations (mean p ¼ 0.00577, mean Battey et al. 2020). We expect to observe young neutral alleles h ¼ 0.00597; Fabian et al. 2012) and 48 European populations at low frequencies among geographically close populations, (mean p ¼ 0.0051, mean h ¼ 0.0052; Kapun et al. 2020). reflecting isolation by distance. We tested this hypothesis by Estimates of Tajima’s D were positive when using PoolSNP, estimating the average geographic distance among pairs of and slightly negative using SNAPE. These results are expected populations that share SNPs only occurring in these two given biases in the detection of rare alleles between these two populations (doubletons), among three populations that SNP calling methods. In addition, our estimates for p, share tripletons, and so forth. Without imposing a MAF filter, Watterson’s h and Tajima’s D were positively correlated with both SNAPE-pooled and PoolSNP pipelines produced pat- previous estimates for the 48 European populations analyzed terns concordant with the expectation. That is, populations by Kapun et al. (2020) (all P < 0.01). Notably, slightly lower in close proximity were more likely to share rare mutations levels of Tajima’s D in North America as compared with relative to random chance pairings (fig. 8A). Notably, SNPs both Africa and Europe (fig. 6C) may be indicative for admix- identified in less than 25 populations tend to be geographi- ture (Stajich and Hahn 2005), which has been identified previ- cally closer in PoolSNP, relative to SNAPE-pooled. The pri- ously along the North American east coast (Caracristi and mary source of this discrepancy between callers occurs when Schlötterer 2003; Kao et al. 2015; Bergland et al. 2016). evaluating SNPs shared by just two populations (fig. 8B). In the case of PoolSNP, only 0.0006% of all SNPs are private to Phylogeographic Clusters in D. melanogaster just two populations and the mean geographical distance is We performed PCA on the PoolSNP variants using samples 702 km. In the case of SNAPE-pooled, 9.3% of all SNPs are from the North American (DrosRTEC), European (DrosEU), private to two populations and the mean distance is and African (DGN) data sets (excluding all Asian and 2,000 km. Aside from the case of n ¼ 2, the difference in Oceanian samples). Prior to analysis, we filtered the joint proximity estimates between the callers is minimal. These data sets to include only high-quality biallelic SNPs. Because findings suggest that some of the SNAPE-pooled SNPs which LD decays rapidly in Drosophila (Comeron et al. 2012), we only only segregate in two populations or less might be false pos- considered SNPs at least 500 bp away from each other. PCA on itives. To further evaluate these geographical patterns, we the resulting 100,000 SNPs revealed evidence for discrete phy- estimated the probability that any given population pair logeographic clusters that correspond to geographic regions belongs to a particular phylogeographic cluster (supplemen- (supplementary fig. S14B, Supplementary Material online). tary fig. S16, Supplementary Material online) as a function of PC1 (24% variance explained [VE]) partitions samples between their shared variants. Our results indicate that rare variants, Africa and the other continents (fig. 7A). PC2 (9% VE) separates private to geographically proximate populations, are strong European from North American populations, and both PC2 predictors of phylogeographic provenance (see fig. 8C). and PC3 (4% VE) divide Europe into two population clusters (fig. 7B). As expected, North American samples are intermedi- Geographically Informative Markers ate to European and African samples, presumably due to re- An inherent strength of our broad biogeographic sampling is cent secondary contact (Kao et al. 2015; Pool 2015; Bergland et the potential to generate a panel of core demography SNPs to al. 2016). Notably, these spatial relationships become evident investigate the provenance of current and future samples. We when PCA projections from each sample are plotted onto a created a panel of geographically informative markers (GIMs) world map (fig. 7C). Interestingly, the emergent clusters in by conducting a discriminant analysis of principal compo- Europe are not strictly defined by geography. For example, nents (DAPC) to discover which loci drive the phylogeo- the western cluster (diamonds in fig. 7D) includes Western graphic signal in the data set. We trained two separate Europe as well as Finland, Turkey, Cyprus, and Egypt. The east- DAPC models: the first utilized the four phylogeographic ern cluster, on the other hand, consists of several populations clusters identified by principal components (PCs; fig. 6A collected in previous Soviet republics as well as Poland, and B and supplementary fig. S16 and table S1, Hungary, Serbia and Austria. Below, we use demographic Supplementary Material online); the second utilized the geo- modeling to resolve the split time between these clusters. graphic localities where the samples were collected (i.e., coun- A unique feature of this data set is that it contains a mix- tries in Europe and the U.S. states). This optimization ture of Pool-Seq and inbred (or haploid) genome data. For indicated that the information contained in the first 40 PCs some geographic regions, the DEST data set contains both maximizes the probability of successful assignment (fig. 9A). data types. Inbred and Pool-Seq samples from nearby geo- This resulted in the inclusion of 30,000 GIMs, most of which graphic regions clustered in the same regions of PC space were strongly associated with PCs 1–3 (fig. 9B inset). (supplementary fig. S15, Supplementary Material online). Moreover, the correlations were larger among the first 3 Excluding the DGN-derived African samples, no PC was PCs and decayed monotonically for the additional PCs (fig. 5790 DEST Population Genomics . doi:10.1093/molbev/msab259 MBE A South Africa B P Zimbabwe Ukraine Rwanda U 20 Russia 10 PC 3 (3.9% VE) erland PC 2 (9% VE) Guinea Ken It pt South Africa Sou Nigeria gypt T Zimbabwe Kenya 0 0 es B H ada Guadeloupe g France Guadeloupe Guinea Canada S −20 United King gdom Unit om −10 0 50 100 0 50 100 PC 1 (24% VE) PC 1 (24% VE) C PC 1 projections onto map D PC 3 projections onto map 60°N 60°N 40°N 55°N Latitude 20°N 50°N 0° 45°N 40°N 20°S 35°N 40°S 100°W 50°W 0° 10°W 0° 10°E 20°E 30°E 40°E Projections on PC1 Projections on PC3 Longitude Longitude 0 30 60 90 -10 -5 0 5 10 15 FIG. 7. Demographic signatures of the DrosEU, DrosRTEC, and DGN data (using the PoolSNP pipeline). (A) PCA dimensions 1 and 2. The mean centroid of a country’s assignment is labeled. (B) PCA dimensions 1 and 3. (C) Projections of PC1 onto a World map. PC1 projections define the existence of continental level clusters of population structure (indicated by the shapes circles: Africa; triangles: North America; diamonds and squares: Europe). (D) Projections of PC3 onto Europe. These projections show the existence of a demographic divide within Europe: the diamond shapes indicate a western cluster, whereas the squares represent an eastern cluster. For panels (C) and (D), the intensity of the color is proportional to the PC projection. The black dashed line shows the two-cluster divide. A 8000 B Number of C PoolSNP populations 0.75 Mean Distance in SNAPE (Km) 4000 200 6000 0.50 P(Phylogenetic Cluster) 150 100 0.25 Distance (Km) 3500 50 24 22 0.00 4000 20 18 16 23 SNAPE 3000 11 15 21 19 2000 12 17 0.75 10 Expected 9 14 0.50 2500 8 7 13 SNAPE 0.25 0 56 PoolSNP 3 4 0.00 2 2000 0 50 100 150 200 1000 2000 3000 4000 0 20 40 60 Mean Distance in PoolSNP (Km) Number of populations Number of populations FIG. 8. Geographic proximity analysis. (A) Average (local regression; LOESS) geographic distance between populations that share a polymorphism at any given site for PoolSNP and SNAPE-pooled. The x-axis represents the number of populations considered; the y-axis is the mean geographic distance among samples. The yellow line represents the random expectation calculated as random pairings of the data. The band around the lines is the standard deviation of the estimator. (B) Correlation graph showing the different mean distance estimate for both callers as a function of the number of populations (the groups from n ¼ 2 to n ¼ 25 are labeled in the graph). A 1-to-1 line is also shown. (C) Probability that all populations containing a polymorphic site come from the same phylogeographic cluster (as defined by PC space, fig. 7 and supplementary fig. S14, Supplementary Material online). The y-axis is the probability of “x” populations belonging to the same phylogeographic cluster. The axis only shows up to 60 populations since, after 40 populations, the probabilities approach 0. The colors are consistent across panels. 5791 Kapun et al. . doi:10.1093/molbev/msab259 MBE FIG. 9. Geographically informative markers. (A) Number of retained PCs which maximize the DAPC model’s capacity to assign group membership. Model trained on the phylogeographic clusters (dashed lines) or the country/state labels (solid line). (B) Absolute correlation for the 33,000 individual SNPs with highest weights onto the first 40 components of the PCA. Inset: Number of SNPs per PC. (C) Location of the 33,000 most informative demographic SNPs across the chromosomes. (D) LOOCV of the DAPC model trained on the phylogeographic clusters. (E) LOOCV of the DAPC model trained on the phylogeographic state/country labels. For panels (D) and (E), the y-axis shows the highest posterior produced by the prediction model and the x-axis is the posterior assigned to the actual label classification of the sample. Also, for (D) and (E), marginal histograms are shown. 9B). Lastly, our GIMs were uniformly distributed across the fly time between pairs of populations sampled throughout genome (fig. 9C). Europe. This is motivated by the observation that the two We assessed the accuracy of our GIM panel using a leave-one- European clusters have different levels of genetic variation. out cross-validation approach (LOOCV). We trained the DAPC The eastern cluster (E) is largely self-contained to Eastern model using all but one sample and then classified the excluded Europe and harbors the lowest levels of hp(0.0049, 95% CI sample. We performed LOOCV separately for the phylogeo- ¼ 0.0047–0.0050). The western cluster (W), on the other graphic cluster groups, as well as for the state/country labels. hand, contains populations from Western Europe as well as The phylogeographic model used all DrosRTEC, DrosEU, and Finland, Turkey, Cyprus, and Egypt, thus making it geograph- DGN samples (excluding Asia and Oceania with too few indi- ically heterogeneous. The western cluster harbors higher lev- viduals per sample); the state/country model used only samples els of hp relative to its eastern counterpart (0.0052, 95% CI ¼ for which each label had at least three or more samples. Our 0.0050–0.0054). Consequently, both clusters harbor statisti- results showed that the model is 100% accurate in terms of cally different levels of genetic variation (t-test; t-value ¼ resolving samples at the phylogeographic cluster level (fig. 9D) 5.22, degrees of freedom [df] ¼ 332.96, P ¼ 3.1010-7), and 89% at the state/country level (fig. 9E). We anticipate that thus suggesting potentially different demographic histories. this set of GIMs will be useful to validate the geographic origin of We tested whether the split time between eastern and west- samples in future sequencing efforts (i.e., identify sample swaps; ern D. melanogaster populations was older than within clus- Nunez et al. 2021) and to study patterns of migration. We note ters, and whether split time was positively correlated with that although Drosophila populations evolve over short time- geographic distance. Prior to addressing this hypothesis, we scales in temperate orchards, samples collected over multiple first evaluated the behavior of the PoolSNP and SNAPE- years were predicted with 89% accuracy in our LOOCV analysis, pooled data sets in demographic inference and also evaluated suggesting that these markers will be valuable for future samples. different methods for converting Pool-Seq data for use with We provide a tutorial on the usage of the GIMs in supplemen- site frequency spectrum-based analysis. tary methods, Supplementary Material online. Prior to estimating divergence times between and among the European clusters, we assessed the behavior of our mo- Estimating the Divergence Time between European ment implementations using the summary statistic h across Genetic Clusters models. We chose h because it has a well-estimated value in The DEST data set can be used to test comparative hypoth- D. melanogaster (h¼4Nel¼0.005; Lack et al. 2016) and thus eses of demographic history. We examined the divergence can serve as a biologically informed calibration parameter. We 5792 DEST Population Genomics . doi:10.1093/molbev/msab259 MBE conducted these preliminary assessments in our simplest example for nui is shown in fig. 10B). Thus, for the remainder model, SþSyM. Our results reveal conspicuous differences of our analysis, we only report the binomial method. between SNAPE-pooled and PoolSNP. PoolSNP produces pre- We used AIC to test which of the four demographic mod- cise estimates of h around the biological expectation, yet els best fit the data: population divergence with symmetric SNAPE-pooled estimates are imprecise and often converge migration (SþSyM), population divergence with asymmetric to the bounds of the estimator (fig. 10A). This behavior is migration (SþAsyM), population divergence followed by a consistent for both AF discretization methods (binomial and bottleneck and growth with symmetric migration counts). PoolSNP results also vary as a function of the AF (SþBGþSyM), or population divergence followed by a bot- discretization method. Based on these results, we chose to tleneck and growth with asymmetric migration use only PoolSNP data for the implementation of our demo- (SþBGþAsyM). We find that the SþAsyM was the best graphic inference. model 71.5% of the time, followed by SþSyM 26.6% of the To further evaluate the behavior of PoolSNP’s estimates as time. Our more complex models (SþBGþSyM and a function of the AF discretization method, we explored SþBGþAsyM) were not generally favored by AIC (fig. 10C). values of the raw parameter outputs by moments. We ex- We also evaluated dAIC, the difference in AIC between the plored the values of the nui parameter (the ancestral popu- best and all other models. We found that SþAsyM and lation size; see Materials and Methods). In general, the counts SþSyM are generally the best models, whereas SþBGþSyM method produced nui estimates which are sparser and less and SþBGþAsyM underperform by at least four orders of stable, by an order of magnitude, relative to binomial draws magnitude in terms of AIC (fig. 10D). We further evaluated (nui sdbinom¼0.938, nui sdcounts¼3.72). In addition, nui AIC performance as a function of the number of completed generated from the counts method produce highly skewed runs. As described in the Materials and Methods, these de- distributions, particularly for jSFS estimated for population mographic inferences are computationally expensive and not pairs in eastern Europe (fig. 10B). Similar to SNAPE, estimates all models ran 50 times in the allotted time. This is of partic- from the counts method also showed the problematic ten- ular concern because all SþAsyM/SyM models ran 50 itera- dency to converge toward the parameter bounds (an tions, whereas SþBGþSyM/AsyM ran, on average, 44.7 and FIG. 10. Optimizing demographic models. (A) Estimates of h from moments as a function of input data: PoolSNP (positive distribution) or SNAPE (negative distribution). We also show the AF discretization method (binomial, “binom,” top; counts, bottom). (B) Distribution of the parameter nui produced by moments as a function of AF discretization strategy. The three colors represent pairwise comparisons done within and across demographic clusters identified via PCA above. Specifically, pink: within eastern clusters (EE), blue: between clusters (EW), and green: within western clusters (WW). (C) Proportion of times a given model was determined to be the best according to AIC. (D) Distribution of d(AICbest), the difference between the best model’s AIC, and all other evaluated models. The y-axis shows the proportion of time a given model appeared in a given d(AICbest) bin. Because the models were Log10transformed, all values were shifted by þ1 (to avoid Log10(0)¼Undefined). Colors corre- spond to model type as labeled in the plot. 5793 Kapun et al. . doi:10.1093/molbev/msab259 MBE 35.2 times, respectively. As such, there is an inherent risk that the effective number of migrants per generation was higher the more complex models (SþBGþSyM/AsyM) did not find for west-into-east migration (Mwest!east ¼ 0.209 flies/gen; the best possible solution. We explored this possibility by 95% CI ¼ 0.169–0.250) as compared with the opposite direc- partitioning dAIC as a function of the number of runs com- tion (Meast!west ¼ 0.178 flies/gen; 95% CI ¼ 0.161–0.196). pleted (supplementary fig. S17, Supplementary Material on- line). Our results indicate that the dAIC of the SþBGþSyM/ Discussion AsyM does not improve among population pairs which run Here we have presented a new, modular, and unified bioin- 40þ or the full 50 iteration cycles. This suggests that our AIC formatics pipeline for processing, integrating and analyzing behavior is not a byproduct of the computational limit on SNP variants segregating in population samples of D. mela- iteration times. We also evaluated the residuals for the four nogaster. We have used this pipeline to assemble the largest demographic models, averaged across all population pairs worldwide data repository of genome-wide SNPs in D. mela- that we contrasted (supplementary fig. S18, Supplementary nogaster to date, based both on previously published data Material online). These results show that, in general, all mod- (DGN: Africa; Lack et al. 2015, 2016) as well as on new data els slightly underestimate rare variants (<10%) and slightly collected by our two collaborating consortia (DrosRTEC: overestimate variants between 10% and 35%. For the remain- mostly North America; Machado et al. 2021; DrosEU: mostly der of analysis, we used the SþAsyM model to estimate di- Europe; Kapun et al. 2020). We assembled this data set using vergence times among populations. two SNP calling strategies that differ in their ability to identify Our analyses suggest that the eastern and western demo- rare polymorphisms, thereby enabling future work studying graphic clusters diverged, on average, 1,013 years ago (95% CI ¼ the evolutionary history of this species. We are dubbing this 887–1,139 years; median ¼ 715 years; fig. 11A). Consistent with data repository and the supporting bioinformatics tools DEST. biological expectation, divergence estimates within population The DEST data repository was built using two different clusters were lower than between clusters. For example, the SNP calling pipelines, SNAPE-pooled (Raineri et al. 2012) eastern cluster is estimated to have a mean divergence within and PoolSNP (Kapun et al. 2020). These two methods differ populations of 294 years (95% CI ¼ 225–362 years; median ¼ fundamentally in their approach to SNP identification, yield 231 years). The western cluster has a mean divergence within data sets amenable to different types of analyses and each populations of 648 years (95% CI ¼ 627–668 years; median ¼ approach has its own specific limitations. The fundamental 626 years). We evaluated the relationship between spatial dis- difference between the data sets produced by these methods tance and divergence time. Similar to our proximity analysis is the number of rare and endemic SNPs identified. This dif- (fig. 8), the biological expectation is that populations in close ference will result in biased estimates of parameters from site proximity are likely to display low divergence estimates. Our frequency spectrum-based demographic models. As a conse- results fit with this expectation, with neighboring populations quence, some care should be taken when interpreting differ- within clusters displaying low divergence estimates (fig. 11B). ent analyses based on these data sets. Lastly, we estimated other population genetic parameters of SNAPE-pooled treats each Pool-Seq sample separately and these population clusters such as effective population size (NE) calculates the posterior probability that a site is polymorphic and migration rates (M). Our estimates of NE suggest that the based on read depth, alternate allele count, and a prior esti- western cluster has larger NE (NE j west¼ 84,921; 95% CI ¼ mate of nucleotide diversity; this approach was designed to 83,373–86,468) relative to the eastern cluster (NE j east¼ identify rare polymorphisms and has been validated using 62,287; 95% CI ¼ 60,207–64,368). In terms of asymmetrical both simulations and empirical approaches (Guirao-Rico migration rates between clusters, our findings show that and Gonzalez 2021). Here, we also provide evidence that A B Log10(Divergence Time) Demographic Cluster 648 y. 4 WW 3 1013 y. EW 2 294 y. 1 CorEE = 0.58, P = 2.2x10-16 EE CorWW = 0.40, P = 2.2x10-16 0 CorEW = 0.22, P = 1.0x10-10 0 1 2 3 4 0 10 20 30 40 Log10(Divergence Time) Distance (km x 100) FIG. 11. Demographic inference of European clusters. (A) Estimates of divergence time between and within the European clusters, pink: within eastern clusters (EE), blue: between clusters (EW), and green: within western clusters (WW). (B) Divergence time as a function of the geographic distance between population pairs. Color palette is consistent with panel (A). Correlation values are shown in the figure. 5794 DEST Population Genomics . doi:10.1093/molbev/msab259 MBE rare and private SNPs identified by SNAPE-pooled are Pool-Seq samples (supplementary fig. S15, Supplementary enriched for true positives (fig. 8) after applying rigorous fil- Material online and discussed in the Results). Thus, there tering and excluding 20 population samples likely affected by does not appear to be significant bias when combining these problems during library preparation which may have resulted data sets, at least when integrating information across the in elevated error rates. genome. Nonetheless, some care should be taken when inter- The dataset based on SNAPE-pooled could therefore be preting allele frequency differences based on data sets gener- useful for studies that rely on rare SNPs, such as those inves- ated by different means. However, any real-time monitoring tigating recent demographic events (Keinan and Clark 2012). activity will likely suffer from the rapidly changing landscape SNAPE-pooled has several limitations though. First, it is only of sequencing technologies. capable of handling Pool-Seq data. Second, because of the One of the biggest challenges in the present “omics” era is hard filtering that we are imposing with our posterior prob- the rapidly growing number of complex large-scale data sets ability cut-off, some true SNPs are being called as missing which require technically elaborate bioinformatics know-how data (see Materials and Methods). This problem is apparent to become accessible and utilizable. This hurdle often prohib- when comparing the number of polymorphisms identified its the exploitation of already available genomics data sets by by SNAPE-pooled and PoolSNP (fig. 3). Third, any demo- scientists without a strong bioinformatics or computational graphic inference done with SNAPE must be limited to cases background. To remedy this situation for the Drosophila evo- where a SNP is discovered in at least three populations or lution community, our bioinformatics pipeline is provided as more, because the caller appears to produce too many false a Docker image (to standardize across software versions, as positives when only two populations are considered (see fig. well as make the pipeline independent of specific operating 8B and our demographic inference with moments, which systems) and a new genome browser makes our SNP data set uses a pairwise, two-population, model). In addition, studies available through an easy-to-use web interface (see supple- that rely on the SNAPE-pooled data set should exclude the mentary figs. S2 and S3, Supplementary Material online; avail- 20 samples we flagged here (fig. 2A and supplementary table able at https://dest.bio, last accessed September 6, 2021). S1, Supplementary Material online). The DEST data repository and platform will enable the PoolSNP, on the other hand, is useful for analysis of com- population genomics community to address a variety of mon variants and allows studying aspects of population longstanding, fundamental questions in ecological and evo- structure and local adaptation based on shared polymor- lutionary genetics. The current data set might for instance phism. Such analyses could include the inference of migra- be valuable for providing a more accurate picture of the tion out of Africa (Kapopoulou et al. 2020), admixture demographic history of D. melanogaster populations, in (Bergland et al. 2016), and back migration to Africa (Pool particular in Europe and North America, and with respect and Aquadro 2006). PoolSNP is an extension of the ap- to multiple bouts of out-of-Africa migration and recent proach developed elsewhere (Kofler, Orozco-terWengel, patterns of admixture. Such analyses can be strongly af- et al. 2011; Kofler, Pandey, et al. 2011). PoolSNP necessarily fected by chromosomal inversions that are known to im- has a limited capacity to identify rare and private SNPs be- pact LD and haplotype variation (Kapun and Flatt 2019; cause it imposes global MAC and allele frequency filters. Durmaz et al. 2020). We have therefore provided frequency Therefore, the more populations that are used for SNP call- estimates for the seven most common cosmopolitan inver- ing by PoolSNP, the less likely PoolSNP is to identify private sions (In(2L)t, In(2R)NS, In(3L)P, In(3R)C, In(3R)K, In(3R)Mo, polymorphisms. Because PoolSNP filters out rare and private and In(3R)Payne; Lemeunier and Aulard 1992), which polymorphisms, it is less sensitive to sequencing or library allows accounting for the effects of inversions in popula- preparation errors. Notably, the 20 flagged populations do tion genetic inference (e.g., Kapopoulou et al. 2020). not have elevated pN/pS with MAC > 50. Additionally, The DEST data set will likewise be useful for an improved Kapun et al. (2020) demonstrated that these problematic understanding of the genomic signatures underlying both samples did not affect population genetic inference based global and local adaptation, including a more fine-grained on common SNPs. The problematic samples derived from view of selective sweeps, their evolutionary origin and distri- the DrosRTEC studies likely do not have a major impact on bution (e.g., see Glinka et al. 2003; Beisswanger et al. 2006; their results either as both Bergland et al. (2014) and Ometto et al. 2005; Stephan 2016; Kapun et al. 2020). In terms Machado et al. (2021) imposed stringent MAF filters. of local adaptation, the broad spatial sampling across latitu- PoolSNP has the added advantage that it can incorporate dinal and longitudinal gradients on the North American and in-silico pooled data sets wherein haplotype or genotype in- European continents, encompassing a broad range of climate formation are collapsed into allele frequencies (see Materials zones and areas of varying degrees of seasonality, will allow and Methods). We took this approach by incorporating the examining the parallel nature of local (clinal) adaptation in Drosophila Genome Nexus data set (DGN; Lack et al. 2016), a response to similar environmental factors in greater depth data set that amalgamates whole-genome sequencing of in- than possible before (e.g., Turner et al. 2008; Kolaczkowski et bred line data and haploid embryos from samples collected al. 2011; Fabian et al. 2012; Bergland et al. 2014, 2016; around the world. Although the DGN data was originally Reinhardt et al. 2014; Kapun et al. 2016, 2020; Waldvogel et generated by multiple labs and run through a different map- al. 2020; Bogaerts-Marquez et al. 2021; Machado et al. 2021). ping pipeline than what we used for the Pool-Seq data, these Another major opportunity provided by the DEST data set samples appear to cluster tightly with geographically close lies in studying the temporal dynamics of evolutionary 5795 Kapun et al. . doi:10.1093/molbev/msab259 MBE change. Sampling at dozens of localities across the growing across comparisons. This behavior is due to the inherent noise season and over multiple years will help to advance our un- of directly converting Pool-Seq AFs (which are heavily af- derstanding of the short-term population and evolutionary fected by coverage) to counts. Based on these observations, dynamics of flies living in diverse environments, thereby pro- we recommend the use of the binomial method of AF dis- viding novel insights into the nature of temporally varying cretization for Pool-Seq analysis (Thia and Riginos 2019). selection (Bergland et al. 2014; Wittmann et al. 2017; Second, we also observed a difference in the estimator’s be- Machado et al. 2021) and evolutionary responses to climate havior based on whether the PoolSNP or SNAPE-pooled data change (e.g., Umina 2005; Rodrıguez-Trelles et al. 2013; were used to build the SFS. In general, PoolSNP generated h Waldvogel et al. 2020). estimates which converge toward 0.005, the biological expec- Moreover, by integrating these worldwide estimates of al- tation for Drosophila. SNAPE-pooled estimates, on the other lele frequencies, those from lab- and field-based “evolve and hand, produced h distributions with high variance as well as a resequence” experiments (E&R; Turner et al. 2011; reviewed in tendency to converge toward the edge of the prior. Kofler and Schlötterer 2014; Schlötterer et al. 2014; Flatt 2020) Interestingly, this type of run-to-the-edge pathological behav- and those from mesocosm experiments (e.g., Rudman et al. ior has been previously characterized (Rosen et al. 2018) and is 2019; Erickson et al. 2020), we might be able to gain deeper generally caused by two possible reasons: over-specified mod- insights into the genetic basis and evolutionary history of els, or, alternatively, noisy input SFS data. Given the relative variation in fitness components (e.g., Flatt 2020). simplicity of the model used for optimization (SþSyM; diver- In addition to analyses of selection, the DEST data set can gence with symmetrical-migration), it is likely that SNAPE’s also be used for preliminary demographic inference. Although SNP calling approach is producing a high number of false Pool-Seq data sets lack important haplotype information, positives which affect model convergence (see also they have been successfully used in the past to generate de- Geographic Proximity Analysis and fig. 8). We therefore rec- mographic and biogeographic insights into both model and ommend PoolSNP over SNAPE-pooled for the purposes of non-model species (e.g., Gautier et al. 2021; Machado et al. exploring or testing demographic hypotheses in cases where 2021; Nunez et al. 2021; see fig. 7). Our analyses suggest that only two populations are considered. Pool-Seq data can be used for demographic model inference. Although our analyses of the DEST sequencing data al- A major caveat in this endeavor is that, to the best of our ready led to novel insights into the evolutionary history of knowledge, Pool-Seq has not been exhaustively benchmarked Drosophila, we believe that the real value of the DEST data set for demographic inference. As such, and until proper valida- lies in the future: its long-term utility will grow as natural and tion has been completed, we present our results as tools for experimental populations are continually being sampled, re- hypothesis generation and exploration. sequenced and added to the repository by the community of Our results from moments are in full agreement with Drosophila evolutionary geneticists. The pipeline that we have basic biological expectations. For example, our estimates of established will make future updates to the data repository h are concordant with previously reported values (0.005; Lack et al. 2016). Moreover, our estimate of mean divergence straightforward. Furthermore, because it is not easily feasible time between the eastern and western European clusters of for any single research group to sample flies densely through D. melanogaster is 1,013 years. This estimate is subject to time and across a broad geographic range, the growing value caveats, given the nature of Pool-Seq data and that future of the DEST data set will depend upon the synergistic collab- validation may need to be done using different types of data. oration among research groups across the globe, as exempli- Nevertheless, we note that this value is plausible as it is fied by the DrosRTEC and DrosEU consortia. Importantly, in well within the newer estimates for Drosophila’s expansion an era of rapidly decreasing sequencing costs, comprehensive into Europe from Africa (4,139 years; Kapopoulou et al. population genomic analyses are no longer limited by genetic 2020). Although previous studies estimated D. melanogast- marker density but by the availability of biological samples er’s European expansion to have occurred around from standardized, collaborative long-term collection efforts 13,000 years ago (e.g., Li and Stephan 2006; Hutter et al. through space and time (e.g., Kapun et al. 2020; Machado et 2007; Laurent et al. 2011), Kapopoulou et al. (2020) showed al. 2021). In this vein, the collaborative framework presented that accounting for the role of asymmetric migration and here might allow us, as a global community, to fill some im- admixture reduces the estimated divergence time between portant gaps in the current data repository: for example, continents. Moreover, our mean estimates of NE for each many areas of the world (notably Asia and South America) cluster (NE j east¼ 62,287, NE j west¼ 84,921) are also remain largely uncharted territory in Drosophila population within Kapopoulou et al.’s (2020) confidence interval for genomics, and the addition of phased sequencing data (e.g., modern European D. melanogaster NE (67,444633,186). providing information on haplotypes, LD, linked selection) Our analyses also revealed two notable behaviors that are will be crucially important for future analyses of demography, relevant to demographic analysis of Pool-Seq data. First, we selection, and their interplay. observed a remarkable difference between the method used We are convinced that the DEST platform will become a to discretize AFs from Pool-Seq, prior to SFS estimation. valuable and widely used resource for scientists interested in Discretizing the data based on direct counts results in noisier Drosophila evolution and genetics, and we actively encourage demographic estimates. Discretizing based on binomial prob- the community to join the collaborative effort we are seeking abilities, on the other hand, produced consistent results to build. 5796 DEST Population Genomics . doi:10.1093/molbev/msab259 MBE Materials and Methods substrate or using baited traps (see Behrman et al. 2018; Machado et al. 2021 for details). Samples were either field- Data Sources caught flies (n ¼ 227), from F1 offspring of wild-caught The genomic data set presented here has been assembled females (n ¼ 7), from a mixture of F1 and wild-caught flies from a combination of Pool-Seq libraries and in silico pooled (n ¼ 7), or from flies kept as isofemale lines in the laboratory haplotypes. We combined 246 Pool-Seq libraries of popula- for five generations or less (n ¼ 4); see supplementary table tion samples from Europe, North America, and the Caribbean S1, Supplementary Material online for more information. To that were sampled through space and time by two collabo- minimize cross-contamination with the closely related sym- rating consortia in North America (DrosRTEC: https://web. patric sister species D. simulans, we only sequenced male D. sas.upenn.edu/paul-schmidt-lab/dros-rtec/, last accessed melanogaster specimens, allowing for higher confidence dis- September 6, 2021) and Europe (DrosEU: http://droseu.net, crimination between the two species based on the morphol- last accessed September 6, 2021) between 2003 and 2016. Of ogy of male genitalia (Capy and Gibert 2004; Markow and these 246 Pool-Seq samples, 121 samples represent previously O’Grady 2006). Samples were stored in 95% ethanol at 20 C unpublished samples generated by DrosEU, 48 DrosEU sam- before DNA extraction. ples previously reported in Kapun et al. (2020), and 77 sam- ples previously reported in Machado et al. (2021). In addition, DNA Extraction and Sequencing we integrated genomic data from >900 inbred or haploid The DrosEU and DrosRTEC consortia centralized extractions genomes from 25 populations in Africa, Europe, Australia, and from pools of flies. DNA was extracted either using North America available from the Drosophila Genome Nexus chloroform/phenol-based (DrosEU: Kapun et al. 2020) or lith- data set (DGN v1.1; Pool et al. 2012; Langley et al. 2012; ium chloride/potassium acetate extraction protocols Grenier et al. 2015; Kao et al. 2015; Lack et al. 2015, 2016) (DrosRTEC: Bergland et al. 2014; Machado et al. 2021) after We further included the D. simulans haplotype (w501; Hu et al. homogenization with bead beating or a motorized pestle. 2013), built as part of the DGN data set, as an outgroup, DrosEU samples from the 2014 collection were sequenced making this repository of 272 (246 Pool-Seq þ 25 DGN þ on an Illumina NextSeq 500 sequencer at the Genomics 1 D. simulans) whole-genome sequenced samples the largest Core Facility of the Pompeu Fabra University in Barcelona, data set of genome-wide SNP polymorphisms available for D. Spain. Libraries of the previously unpublished DrosEU samples melanogaster to date. from 2015 and 2016 were constructed using the Illumina TruSeq PCR Free library preparation kit following the manu- Metadata facturer’s instructions and sequenced on the Illumina HiSeq X We assembled uniform metadata for all samples (supplemen- platform as paired-end fragments with 2  150 bp length at tary table S1, Supplementary Material online). This informa- NGX Bio (San Francisco, California, USA). The previously tion includes collection coordinates, collection date, and the published samples of the DrosRTEC consortium were pre- number of flies per sample. Samples are also linked to biocli- pared and sequenced on GAIIX, HiSeq2000, or HiSeq3000 matic variables from the nearest WorldClim (Hijmans et al. platforms, as described in Bergland et al. (2014) and 2005) raster cell at a resolution of 2.5 and to weather stations Machado et al. (2021). For information on DNA extraction from the Global Historical Climatology Network (GHCND; and sequencing methods of the various DGN samples, see ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/daily/) to allow for fu- Lack et al. (2016) and others (Langley et al. 2012; Pool et al. ture analyses of the environmental drivers that might underlie 2012; Grenier et al. 2015; Kao et al. 2015). genetic change. We also provide summaries of basic attributes of each sample derived from the sequencing data including Mapping Pipeline average read depth, PCR duplicate rate, D. simulans contam- The joint analysis of genomic data from different sources ination rate, relative abundances of non-synonymous versus requires the application of uniform quality criteria and a com- synonymous polymorphisms (pN/pS), the number of private mon bioinformatics pipeline. To accomplish this, we devel- polymorphisms, diversity statistics (Watterson’s h, p, and oped a standardized pipeline that performs filtering, quality Tajima’s D), and estimates of inversion frequencies. control and mapping of any given Pool-Seq sample (see sup- plementary fig. S1, Supplementary Material online). This pipe- Sample Collection line performs quality filtering of raw reads, maps reads to a Most population samples contributed by the DrosEU and the hologenome (see below), performs realignment and filtering DrosRTEC consortia were collected in a coordinated fashion around indels, and filters for mapping quality. The output of to generate a consistent data set with minimized sampling this pipeline includes quality control metrics, bam files, pileup bias. In brief, fly collections were performed exclusively in files, and allele frequency estimates for every site in the ge- natural or seminatural habitats, such as orchards, vineyards, nome (gSYNC, see below). Our pipeline is provided as a and compost piles. For most European collections, flies were Docker image and will facilitate the integration of future collected using mashed banana, or apples with live yeast as samples to extend the worldwide D. melanogaster SNP data bait in traps placed at sampling sites for multiple days to set presented here. attract flies, or by sweep netting (see Kapun et al. 2020 for The mapping pipeline includes the following major more details). For North American collections, flies were col- steps. Prior to mapping, we removed sequencing adapters lected by sweep-net, aspiration, or baiting over natural and trimmed the 30 ends of all reads using cutadapt 5797 Kapun et al. . doi:10.1093/molbev/msab259 MBE (Martin 2011). We enforced a minimum base quality score computed FASTA files (“Consensus Sequence Files” from 18 (-q flag in cutadapt) and assessed the quality of https://www.johnpool.net/genomes.html, last accessed raw and trimmed reads with FASTQC. Trimmed reads September 6, 2021) and calculated allele frequencies at every with minimum length <75 bp were discarded and only intact site, for each population, using custom bash scripts. We cal- read pairs were considered for further analyses. Overlapping culated allele frequencies for each population by summing paired-end reads were merged using bbmerge (v. 35.50; reference and alternative allele counts across all individuals Bushnell et al. 2017). Trimmed reads were mapped against using the precomputed haplotype FASTA files. Because esti- a compound reference genome (“hologenome”) consisting mates of allele frequencies and total allele counts for the DGN of the genomes of D. melanogaster (v.6.12) and D. simulans samples only consider unambiguous IUPAC codes, heterozy- (Hu et al. 2013) as well as genomes of common commensals gous sites or sites masked as N’s in the original FASTA files and pathogens, including Saccharomyces cerevisiae (GCF_ were converted to missing data. We used liftover (Kuhn et al. 000146045.2), Wolbachia pipientis (NC_002978.6), 2013) to translate genome coordinates to Drosophila refer- Pseudomonas entomophila (NC_008027.1), Commensalibacter ence genome release 6 (dos Santos et al. 2015) and formatted intestine (NZ_AGFR00000000.1), Acetobacter pomorum (NZ_ them to match the gSYNC format (described below). Scripts AEUP00000000.1), Gluconobacter morbifer (NZ_AGQV000 for reformatting the DGN data can be found in the GitHub 00000.1), Providencia burhodogranariea (NZ_AKKL0000 repository for this project (https://github.com/DEST-bio/ 0000.1), Providencia alcalifaciens (NZ_AKKM01000049.1), DEST_freeze1, last accessed September 6, 2021). Providencia rettgeri (NZ_AJSB00000000.1), Enterococcus faecalis (NC_004668.1), Lactobacillus brevis (NC_008497.1), and SNP Calling Strategies Lactobacillus plantarum (NC_004567.2), using bwa mem (v. We used two complementary approaches to perform SNP 0.7.15; Li 2013) with default parameters. We retained reads calling. The first was PoolSNP (Kapun et al. 2020), a heuristic with mapping quality greater than 20 as well as those with tool which identifies polymorphisms based on the combined no secondary alignment using samtools (Li et al. 2009). PCR evidence from multiple samples. This approach is similar to duplicate reads were removed using Picard MarkDuplicates other common Pool-Seq variant calling tools (Koboldt et al. (v.1.109; http://broadinstitute.github.io/picard/, last accessed 2009, 2012; Kofler, Orozco-terWengel, et al. 2011; Kofler, September 6, 2021). Sequences were realigned in the proximity Pandey, et al. 2011). PoolSNP integrates allele counts across of insertions–deletions (indels) with GATK (v3.4-46; McKenna multiple independent samples and applies stringent MAC et al. 2010). We identified and removed any reads that mapped and MAF thresholds for variant detection. PoolSNP is to the D. simulans genome using a custom python script, fol- expected to be good at detecting variants present in multiple lowing methods outlined previously (Kapun et al. 2020; populations but is not very sensitive to rare private alleles. The Machado et al. 2021; for a more in-depth analysis of D. simulans second approach was SNAPE-pooled (Raineri et al. 2012), a contamination, see Wallace et al. 2021). Although this method tool that identifies polymorphic sites based on Bayesian in- of decontamination by D. simulans accurately estimates con- ference for each population independently using pairwise tamination rate and removes the vast majority of D. simulans nucleotide diversity estimates as a prior. SNAPE-pooled is reads (Machado et al. 2021), care should be taken when ana- expected to be more sensitive to rare private polymorphisms lyzing samples with higher contamination rates at sites that are (Raineri et al. 2012; Guirao-Rico and Gonzalez 2021). The SNP shared polymorphisms between the two species. calling step is built using the snakemake (Mölder et al. 2021) pipeline and the parameters to run the two callers can be Incorporation of the DGN Data Set found at https://github.com/DEST-bio/DEST_freeze1 (last We incorporated population allele frequency estimates de- accessed September 6, 2021). rived from inbred line and haploid embryo sequencing data from populations sampled throughout the world using an in gSYNC Generation and Filtering silico pooling approach. These samples have been previously Our pipeline utilizes a common data format to encode allele collected and sequenced by several groups (Langley et al. counts for each population sample (SYNC; Kofler, Pandey, 2012; Mackay et al. 2012; Pool et al. 2012; Grenier et al. et al. 2011). A “genome-wide SYNC” (gSYNC) file records 2015; Kao et al. 2015; Lack et al. 2015, 2016) and together the number of A, T, C, and G for every site of the reference form the Drosophila Genome Nexus data set (DGN; Lack et al. genome. Because gSYNC files for all populations have the 2015, 2016). We included 25 DGN populations with 5 indi- same dimension, they can be quickly combined and passed viduals per population, plus the D. simulans haplotype w501 to a SNP calling tool. They can be filtered and are also rela- built as part of the DGN data set. The DGN populations that tively small for a given sample (500 Mb), enabling efficient we used are primarily from Africa (n ¼ 18) but also include data sharing and access. The gSYNC file is analogous to the populations from Europe (n ¼ 2), North America (n ¼ 3), gVCF file format as part of the GATK HaplotypeCaller ap- Australia (n ¼ 1), and Asia (n ¼ 1). The complete list of proach (McKenna et al. 2010) but is specifically tailored to DGN populations, and samples, used in this data set can be Pool-Seq samples. found in supplementary table S1, Supplementary Material We generated gSYNC files for both PoolSNP and SNAPE. online. To generate a PoolSNP gSYNC file, we first converted BAM To incorporate the DGN populations into the DrosEU files to the MPILEUP format with samtools mpileup using the - and DrosRTEC Pool-Seq data sets, we used the pre- B parameter to suppress recalculations of per-base alignment 5798 DEST Population Genomics . doi:10.1093/molbev/msab259 MBE qualities and filtered for a minimum mapping quality with the used this value to tune our choice of MAF and MAC and parameter -q 25. Next, we converted the MPILEUP file con- to identify egregious outlier samples. We found that a global taining mapped and filtered reads to the gSYNC format using MAC ¼ 50 provided qualitatively identical estimates of pN/pS custom python scripts. To generate a SNAPE-pooled gSYNC across all populations (fig. 2B) and that the results were in- file, we ran the SNAPE-pooled version specific to Pool-Seq sensitive to MAF (results not shown). We therefore used data for each sample in MPILEUP format with the following these parameters for genome-wide variant calling (see parameters: h ¼ 0.005, D ¼ 0.01, prior¼‘informative’, Identification and Quality Control of SNP Polymorphisms). fold¼‘unfolded’, and nchr¼number of flies (x2 for autosomes We kept a third heuristic parameter, the missing data rate, and x1 for the X and Y chromosomes) following Guirao-Rico constant at a minimum of 50%. and Gonzalez (2021). We converted the SNAPE-pooled out- To generate the SNAPE-pooled VCF files, we combined the put file to a gSYNC file containing the counts of each allele 246 masked SNAPE-pooled gSYNC files into a two- per position and the posterior probability of polymorphism as dimensional matrix, as described above, and generated a defined by SNAPE-pooled using custom python scripts. We VCF formatted output based on allele counts for any site only considered positions with a posterior probability 0.9 as found to be polymorphic in one or more populations. We being polymorphic and with a posterior probability 0.1 as evaluated pN/pS across a range of local MAF thresholds (fig. being monomorphic. In all other cases, positions were marked 2C) and found that pN/pS is largely insensitive to local MAF, as missing data. once accounting for some problematic samples (see below). We masked gSYNC files for PoolSNP and SNAPE-pooled Final VCF files with annotations from SNPeff (version 4.3; using a common set of filters. Sites were filtered from gSYNC Cingolani et al. 2012) were stored in VCF and BCF (Danecek et files if they had: 1) minimum read depth <10; 2) maximum al. 2011) file formats alongside an index file in TABIX format read depth >the 95% coverage percentile of a given chromo- (Li 2011). Besides VCF files, we also stored SNP data in the somal arm and sample; 3) located within repetitive elements GDS file format using the R package SeqArray (Zheng et al. as defined by RepeatMasker; 4) within 5-bp distance up- and 2017). downstream of indel polymorphisms identified by the GATK IndelRealigner. Filtered sites were converted to missing data in Inversion Frequency Estimates the gSYNC file. The location of masked positions for every We estimated the frequencies of seven cosmopolitan inver- sample was recorded as a BED file. sion polymorphisms (In(2L)t, In(2R)NS, In(3L)P, In(3R)C, In(3R)K, In(3R)Mo, In(3R)Payne) based on a previously pub- VCF Generation lished panel of diagnostic SNP markers that are in tight LD We generated three versions of the variant files, which differ with the corresponding inversions (Kapun et al. 2014). As in their inclusion of the DGN samples and the SNP calling previously described (Kapun et al. 2016), we isolated the strategy. For PoolSNP variant calling, we generated two var- positions in the VCF file of all marker SNPs and estimated iant tables: the first version incorporates all 272 samples of the the frequency of each inversion as the mean frequency of Pool-Seq (DrosRTEC, DrosEU) and in silico Pool-Seq popula- inversion-specific alleles at all marker SNPs. tions (DGN). The second version only considers the 246 Pool- Seq samples excluding the DGN samples (used for compar- Population Genetic Analyses ison to the SNAPE-pooled version). The third file is based on We estimated allele frequencies for each site across popula- SNAPE-pooled and contains 246 Pool-Seq samples only. tions as the ratio of the alternate allele count to the total site To generate the PoolSNP versions, we combined the coverage. We also calculated per-site averages for nucleotide masked PoolSNP-gSYNC files into a two-dimensional matrix, diversity (p, Nei 1987), Watterson’s h (Watterson 1975) and where rows correspond to each position in the reference Tajima’s D (Tajima 1989) across all sites or in nonoverlapping genome and columns describe chromosome, position, and windows of 100, 50, and 10 kb length. To estimate these sum- reference allele, followed by allele counts in SYNC format for mary statistics, we converted masked gSYNC files (with posi- every sample in the data set. This combined matrix was then tions filtered for repetitive elements, low and high read depth, subjected to variant calling using PoolSNP, resulting in a VCF- and proximity to indels; see gSYNC Generation and Filtering) formatted file. We performed SNP calling only for the major back to the MPILEUP format using custom-made scripts. The chromosomal arms (X, 2L, 2R, 3L, 3R) and the 4th (dot) chro- MPILEUP files were processed using npstat v.1 (Ferretti et al. mosome. Data for heterochromatic arms of the autosomes, 2013) with parameters -maxcov 10000 and -nolowfreq m ¼ 0 the Y chromosome, and the mitochondrial genome can be in order to include all filtered positions for analysis. We only extracted from the MPILEUP files provided at https://dest.bio considered sites identified as being polymorphic by PoolSNP (last accessed September 6, 2021). or SNAPE-pooled for analysis, using the -snpfile option of We evaluated the choice of two heuristic parameters ap- npstat. For the DGN populations, chromosome-wide sum- plied to PoolSNP: global MAC and global MAF. Using all 272 mary statistics were estimated only for samples with less samples, we varied MAF (0.001, 0.01, 0.05) and MAC (5-100) than 50% missing data per chromosome. Due to small sample and called SNPs at a randomly selected 10% subset of the sizes, Tajima’s D was not estimated for seven African DGN genome. Based on SNP annotations with SNPeff (version 4.3; populations that consisted of only five haploid embryos. To Cingolani et al. 2012), we calculated pN/pS, which is the ratio compare population genetic estimates between the PoolSNP of nonsynonymous to synonymous polymorphisms, and versus SNAPE-pooled data sets, we performed Pearson’s 5799 Kapun et al. . doi:10.1093/molbev/msab259 MBE correlation on 226 populations present in both data sets (see four basic demographic models: 1) population divergence Identification and Quality Control of SNP Polymorphisms) with symmetric migration (SþSyM), 2) population diver- using the stats package of R v.3.6.3. The effects of pool size gence with asymmetric migration (SþAsyM), 3) population (number of individuals sampled per population) on genome- divergence followed by a bottleneck and growth with sym- wide estimates of p, Watterson’s h and Tajima’s DS estimates metric migration (SþBGþSyM), and 4) population diver- were examined for European and North American popula- gence followed by a bottleneck and growth with tions using the PoolSNP data set and a linear model in R asymmetric migration (SþBGþAsyM). We fit these models v.3.6.3. Finally, for 48 European populations we estimated using the python package moments (Jouganous et al. 2017). Pearson’s correlations between p, Watterson’s h and We converted our data to the moments input format using Tajima’s D as estimated from the PoolSNP data set versus the genomalicious (Thia and Riginos 2019) function dadi_in- previous estimates by Kapun et al. (2020) using the stats puts_pools(), using either the “counts” or the “probs” (here- package of R v3.6.3. after “binomial”) methods. These methods are used to Next, we examined patterns of between-population differ- convert Pool-Seq allele frequency data, which has a variable entiation by calculating window-wise estimates of pairwise denominator (read depth), to the integer-based count of the FST, based on the method from Hivert et al. (2018) imple- site frequency spectrum (SFS) used by moments and other mented in the computePairwiseFSTmatrix() function of the R SFS analyses (Gutenkunst et al. 2009). The “counts” method package poolfstat (v1.1.1). This analysis was performed for the rounds the allele counts to the nearest integer based on the data set composed of 271 samples (all samples excluding the number of chromosomes sampled. The “binomial” method D. simulans reference strain) processed with PoolSNP, focus- generates allele counts based on a binomial draw given the ing on SNPs shared across the whole data set. Finally, we observed allele frequency and the number of chromosomes. averaged pairwise FST within and among phylogeographic For all analyses, we used the mean effective coverage (Feder et clusters identified in our analyses: Africa (17 samples), al. 2012) per population as the number of chromosomes North America (76 samples), Eastern Europe (83 samples), sampled. We only focused on autosomal SNPs and only and Western Europe (93 samples). Samples from China and used populations that passed quality control (fig. 2). Australia were not included due to limited sampling. These Our model estimates a different number of parameters FST tracks at windows sizes of 100, 50, and 10 kb are available depending on its type. For instance, the SþSyM model esti- at https://dest.bio (last accessed September 6, 2021; supple- mates three core parameters: the divergence time between mentary figs. S2 and S3, Supplementary Material online). populations (Ts), the migration rate between populations To assess population structure in the worldwide data set, (mi$j) and the ancestral population sizes (nui). The nui, we applied principal components analysis (PCA), population Ts and mi$j parameters are initially drawn from uniform clustering, and population assignment based on a DAPC priors with user-defined upper boundaries of 10, 5, and 50, (Jombart et al. 2010) to all 271 PoolSNP-processed samples. and lower boundaries of 1.010-5, 1.010-5, and 0, respec- For these analyses, we subsampled a set of 100,000 SNPs tively. The SþAsyM model includes all above parameters, but spaced apart from each other by at least 500 bp. We opti- has explicit asymmetric migration parameters (i.e., mi!j and mized our models using cross-validation by iteratively divid- mj!i) which are also parametrized as uniform distributions ing the data as 90% for training and 10% for learning. We with 0–50 parameter bounds. Models SþBGþSyM and extracted the first 40 PCs from the PCA and ran Pearson’s SþBGþAsyM are similar to their SþAsyM and SþSyM coun- correlations between each PC and all loci. We subsequently terparts, with the addition of the initial (nuiB) and final (nuiF) extracted the top 33,000 SNPs with large and significant cor- sizes of each population. These are also parametrized as uni- relations to PCs 1–40. We chose the 33,000 number as a form distribution bounded between 1.010-5 and 10. compromise between panel size and differentiation power. Overall, we explored the behavior of the estimators for two For example, depending on the number of individuals sur- allele frequency (AF) discretization strategies (counts and bi- veyed, these 33,000 loci can discern genetic differentiation (s) nomial) and two SNP callers (SNAPE and PoolSNP). between two populations with parametric FST of 0.001– Our pipeline estimates a joint SFS (jSFS) from the discre- 0.0001 for sample sizes (n) of 10–1,000. These estimates tized AF data for a given population pair. These are always come from the phase change formula: s  FST ¼ 1/(nm)1/2 folded jSFS to account for unknown ancestral states. For (Patterson et al. 2006). Here, the two populations were sam- computational purposes, we did not evaluate every possible pled for n/2 individuals and genotyped at m ¼ 33,000 pairwise combination in the DEST data set. Instead, we ran- markers. Furthermore, we included SNPs as a function of domly sampled 1,200 population pairs drawn from European the percent variance explained by each PC. PCAs, clustering, populations that passed quality filtering (supplementary table and assignment based DAPC analyses were carried out using S1, Supplementary Material online). The moments simula- the R packages FactoMiner (v. 2.3), factoextra (v. 1.0.7) and tions were run with a maximum of 50 iterations. It is impor- adegenet (v. 2.1.3), respectively. tant to note that running these demographic models is computationally expensive and some individual runs fail to Demographic Inference with Moments converge across the 50 iterations, and thus some models did To evaluate the efficacy of PoolSNP and SNAPE-pooled in not run all 50 times. Nevertheless, we explicitly explored the inferring reasonable demographic parameters, we ran pair- consequences of the total number of completed runs in the wise comparisons of European Drosophila populations under performance of the model selection. 5800 DEST Population Genomics . doi:10.1093/molbev/msab259 MBE Model selection was performed using maximum log- Research Council (ERC) under the European Union’s likelihood and Akaike’s information criterion (AIC) for each Horizon 2020 research and innovation programme (H2020- completed simulation run. For each implementation per pop- ERC-2014-CoG-647900) and by the Spanish Ministry of ulation pair, the simulation with lowest AIC was retained as Science and Innovation (BFU-2011-24397); T.F. by the Swiss the “best fit” for later comparison. The model fit was observed National Science Foundation (SNSF grants PP00P3_133641, in a subset of models run via residuals as well. Raw model PP00P3_165836, and 31003A_182262) and a Mercator parameter outputs were converted to interpretable units in Fellowship from the German Research Foundation (DFG), accordance with the moments manual. To this end, we used held as a EvoPAD Visiting Professor at the Institute for known biological constants for Drosophila, namely l, L, and Evolution and Biodiversity, University of Münster; AOB by generations per year (g). The mutation rate, l, was set to the National Institutes of Health (R35 GM119686); M.K. by 2.810-9 (Keightley et al. 2014). L is the sum of the autosomal Academy of Finland grant 322980; V.L. by Danish Natural chromosome arms minus the median of the number of Science Research Council (FNU) (grant no. 4002-00113B); masked sites across all of the European (DrosEU) samples. FS Deutsche Forschungsgemeinschaft (DFG) (grant no. In moments, outputs are scaled in units of 2Nref, where Nref STA1154/4-1), Project 408908608; J.P. by the Deutsche is the ancestral population size (Nref ¼h/4lL). Divergence Forschungsgemeinschaft Projects 274388701 and time (2Nref Ts) was converted to chronological time assum- 347368302; A.U. by FPI fellowship (BES-2012-052999); ET ing 15 generations per year (Pool 2015). Israel Science Foundation (ISF) (grant no. 1737/17); M.S.V., M.S.R. and M.J. by a grant from the Ministry of Education, Web-Based Genome Browser Science and Technological Development of the Republic of Our HTML-based DEST browser (supplementary fig. S2, Serbia (451-03-68/2020-14/200178); A.P., K.E. and M.T. by a Supplementary Material online) is built on a JBrowse grant from the Ministry of Education, Science and Docker container (Buels et al. 2016), which runs under Technological Development of the Republic of Serbia (451- Apache on a CentOS 7.2 Linux x64 server with 16 Intel 03-68/2020-14/200007); and TM NSERC grant RGPIN-2018- Xeon 2.4 GHz processors and 32 GB RAM. It implements a 05551. The authors acknowledge Research Computing at The hierarchical data selector that facilitates the visualization and University of Virginia for providing computational resources selection of multiple population genetic metrics or statistics and technical support that have contributed to the results for all 271 samples based on the PoolSNP-processed data set, reported within this publication (https://rc.virginia.edu, last taking into account sampling location and date. Importantly, accessed September 6, 2021). our genome browser provides a portal for downloading allelic information and precomputed population genetics statistics Author Contributions in multiple formats (supplementary figs. S2A and C and S3, Martin Kapun: Conceptualization, data curation, formal anal- Supplementary Material online), a usage tutorial (supplemen- ysis, funding acquisition, investigation, methodology, project tary fig. S2B, Supplementary Material online) and versatile administration, resources, software, supervision, visualization, track information (supplementary fig. S2D, Supplementary writing—original draft, writing—review and editing. Joaquin Material online). Bulk downloads of full variation tracks are Nunez: formal analysis, methodology, software, visualization, available in BigWig format (Kent et al. 2010) and Pool-Seq files writing—original draft, writing—review and editing. Marıa (in VCF format) are downloadable by population and/or sam- Bogaerts-Marquez: formal analysis, methodology, software, pling date using custom options from the Tools menu (sup- visualization, writing—original draft, writing—review and plementary fig. S2C, Supplementary Material online). All data, editing. Jes us Murga-Moreno: formal analysis, methodology, tools, and supporting resources for the DEST data set, as well software, visualization, writing—original draft, writing—re- as reference tracks downloaded from FlyBase (v.6.12) (dos view and editing. Margot Paris: formal analysis, methodology, Santos et al. 2015), are freely available at https://dest.bio software, visualization, writing—original draft, writing—re- (last accessed September 6, 2021). view and editing. Joseph Outten: software, writing—review and editing. Marta Coronado-Zamora: formal analysis, meth- Supplementary Material odology, software, visualization, writing—original draft, writ- Supplementary data are available at Molecular Biology and ing—review and editing. Aleksandra Patenkovic: resources. Evolution online. Amanda Glaser-Schmitt: resources. Anna Ullastres: resources. Antonio J. Buendıa-Ruız: resources. Banu S. Onder: resources. Acknowledgments Brian P. Lazzaro: resources, writing—review and editing. We thank four reviewers and the handling editor for helpful Catherine Montchamp-Moreau: resources. Christopher W. comments on previous versions of our manuscript. We are Wheat: resources, writing—review and editing. Cristina P. grateful to the members of the DrosEU and DrosRTEC con- Vieira: resources, writing—review and editing. Daniel K. sortia for their long-standing support, collaboration, and for Fabian: resources. Darren J. Obbard: resources. Dmitry V. discussion. DrosEU was funded by a Special Topic Networks Mukha: resources. Dorcas J. Orengo: resources, writing—re- (STN) grant from the European Society for Evolutionary view and editing. Elena Pasyukova: resources. Eliza Argyridou: Biology (ESEB). M.K. was supported by the Austrian Science resources. Emily L. Behrman: resources, writing—review and Foundation (grant no. FWF P32275); J.G. by the European editing. Eran Tauber: resources. Eva Puerma: resources, 5801 Kapun et al. . doi:10.1093/molbev/msab259 MBE writing—review and editing. Fabian Staubach: resources, writ- References ing—review and editing. Francisco D. Gallardo-Jimenez: Adams MD. 2000. The genome sequence of Drosophila melanogaster. resources. Iryna Kozeretska: resources. J. Roberto Torres: Science 287:2185–2195. resources. Jessica K. Abbott: resources. John Parsch: funding Arguello JR, Laurent S, Clark AG. 2019. Demographic history of the hu- acquisition, resources, writing—review and editing. Jorge man commensal Drosophila melanogaster. Genome Biol Evol. 11(3):844–854. Vieira: resources, writing—review and editing. M. Josefa Assaf ZJ, Tilk S, Park J, Siegal ML, Petrov DA. 2017. Deep sequencing of Gomez-Julian: resources. Katarina Eric: resources. Kelly A. natural and experimental populations of Drosophila melanogaster Dyer: resources. Lain Guio: resources. Lino Ometto: writ- reveals biases in the spectrum of new mutations. Genome Res. ing—review and editing. M. Luisa Espinosa-Jimenez: resour- 27(12):1988–2000. ces. Maaria Kankare: resources, writing—review and editing. Bastide H, Betancourt A, Nolte V, Tobler R, Stöbe P, Futschik A, Schlötterer C. 2013. A genome-wide, fine-scale map of natural pig- Mads F. Schou: resources, writing—review and editing. Maria mentation variation in Drosophila melanogaster. PLoS Genet. P. Garcıa Guerreiro: resources, writing—review and editing. 9(6):e1003534. Marija Savic Veselinovic: resources. Marija Tanaskovic: resour- Battey CJ, Ralph PL, Kern AD. 2020. Space is the place: effects of contin- ces. Marina Stamenkovic-Radak: funding acquisition, resour- uous spatial structure on analysis of population genetic data. ces. Mihailo Jelic: resources. Miriam Merenciano: resources. Genetics 215(1):193–214. Behrman EL, Howick VM, Kapun M, Staubach F, Bergland AO, Petrov Oleksandr M. Maistrenko: writing—review and editing. DA, Lazzaro BP, Schmidt PS. 2018. Rapid seasonal evolution in innate Omar Rota-Stabelli: resources. Sara Guirao-Rico: resources, immunity of wild Drosophila melanogaster. Proc R Soc B Biol Sci. writing—review and editing. Sonia Casillas: resources, writ- 285(1870):20172599. ing—review and editing. Sonja Grath: resources. Stephen W. Beisswanger S, Stephan W, Lorenzo D. 2006. Evidence for a selective Schaeffer: resources. Subhash Rajpurohit: resources. Svitlana sweep in the wapl region of Drosophila melanogaster. Genetics 172(1):265–274. V. Serga: resources. Thomas J.S. Merritt: resources. Vivien Benjamini Y, Speed TP. 2012. Summarizing and correcting the GC con- Horvath: resources. Vladimir E. Alatortsev: resources. Volker tent bias in high-throughput sequencing. Nucleic Acids Res. Loeschcke: resources. Yun Wang: resources. Heather E. 40(10):e72. Machado: resources. Keric Lamb: analysis, writing—original Bergland AO, Behrman EL, O’Brien KR, Schmidt PS, Petrov DA. 2014. draft, writing—review and editing. T^ania Paulo: resources. Genomic evidence of rapid and stable adaptive oscillations over seasonal time scales in Drosophila. PLoS Genet. 10(11):e1004775. Leeban Yusuf: analysis. Antonio Barbadilla: software, writ- Bergland AO, Tobler R, Gonzalez J, Schmidt P, Petrov D. 2016. ing—review and editing. Dmitri Petrov: conceptualization, Secondary contact and local adaptation contribute to genome- funding acquisition, project administration, resources, writ- wide patterns of clinal variation in Drosophila melanogaster. Mol ing—review and editing. Paul Schmidt: conceptualization, Ecol. 25(5):1157–1174. Bogaerts-Marquez M, Guirao-Rico S, Gautier M, Gonzalez J. 2021. funding acquisition, project administration, resources, writ- Temperature, rainfall and wind variables underlie environmental ing—review and editing. Josefa Gonzalez: conceptualization, adaptation in natural populations of Drosophila melanogaster. Mol funding acquisition, project administration, resources, super- Ecol. 30(4):938–954. vision, writing—original draft, writing—review and editing. Buels R, Yao E, Diesh CM, Hayes RD, Munoz-Torres M, Helt G, Goodstein Thomas Flatt: conceptualization, funding acquisition, project DM, Elsik CG, Lewis SE, Stein L, et al. 2016. JBrowse: a dynamic web platform for genome visualization and analysis. Genome Biol. 17:66. administration, resources, supervision, writing—original draft, Burke MK. 2012. How does adaptation sweep through the genome? writing—review and editing. Alan Bergland: conceptualiza- Insights from long-term selection experiments. Proc R Soc B Biol tion, data curation, formal analysis, funding acquisition, inves- Sci. 279(1749):5029–5038. tigation, methodology, project administration, resources, Bushnell B, Rood J, Singer E. 2017. BBMerge – accurate paired shotgun software, supervision, visualization, writing—original draft, read merging via overlap. PLoS One. 12(10):e0185056. Campo D, Lehmann K, Fjeldsted C, Souaiaia T, Kao J, Nuzhdin SV. 2013. writing—review and editing. Whole-genome sequencing of two North American Drosophila mel- anogaster populations reveals genetic differentiation and positive Data Availability selection. Mol Ecol. 22(20):5084–5097. All scripts to make figures and perform analyses associated Capy P, Gibert P. 2004. Drosophila melanogaster, Drosophila simulans: so with this manuscript are available at: https://github.com/ similar yet so different. Genetica 120(1-3):5–16. Caracristi G, Schlötterer C. 2003. Genetic differentiation between DEST-bio/data-paper (last accessed September 6, 2021). American and European Drosophila melanogaster populations could All scripts to build the data set, including the mapping be attributed to admixture of African alleles. Mol Biol Evol. pipeline, SNP calling scripts, and metadata are available at: 20(5):792–799. https://github.com/DEST-bio/DEST_freeze1 (last accessed Celniker SE, Rubin GM. 2003. The Drosophila melanogaster genome. September 6, 2021). All output from the DEST pipeline, in- Annu Rev Genomics Hum Genet. 4:89–117. Cheng C, White BJ, Kamdem C, Mockaitis K, Costantini C, Hahn MW, cluding intermediate output files, metadata, etc. can be found Besansky NJ. 2012. Ecological genomics of Anopheles gambiae along at: https://dest.bio (last accessed September 6, 2021). Datafiles a latitudinal cline: a population-resequencing approach. Genetics available via the website can also be downloaded through the 190(4):1417–1432. command-line interface. The genome browser associated Cingolani P, Platts A, Wang LL, Coon M, Nguyen T, Wang L, Land SJ, Lu X, with the DEST data set can be found at: http://dest-bio.uab.- Ruden DM. 2012. A program for annotating and predicting the effects of single nucleotide polymorphisms. Fly 6(2):80–92. cat, last accessed September 6, 2021. The dockerized mapping Comeron JM, Ratnappan R, Bailin S. 2012. The many landscapes of pipeline can be found at https://hub.docker.com/r/destbio/ recombination in Drosophila melanogaster. PLoS Genet. destbiodocker, last accessed September 6, 2021. 8(10):e1002905. 5802 DEST Population Genomics . doi:10.1093/molbev/msab259 MBE Corbett-Detig R, Nielsen R. 2017. A hidden Markov model approach for in Drosophila melanogaster: a multi-locus approach. Genetics simultaneously estimating local ancestry and admixture time using 165(3):1269–1278. next generation sequence data in samples of arbitrary ploidy. PLoS Gould BA, Chen Y, Lowry DB. 2017. Pooled ecotype sequencing reveals Genet. 13(1):e1006529. candidate genetic mechanisms for adaptive differentiation and re- Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, productive isolation. Mol Ecol. 26(1):163–177. Handsaker RE, Lunter G, Marth GT, Sherry ST, et al. 2011. The variant Grenier JK, Arguello JR, Moreira MC, Gottipati S, Mohammed J, Hackett call format and VCFtools. Bioinformatics 27(15):2156–2158. SR, Boughton R, Greenberg AJ, Clark AG. 2015. Global diversity lines– David JR, Capy P. 1988. Genetic variation of Drosophila melanogaster a five-continent reference panel of sequenced Drosophila mela- natural populations. Trends Genet. 4(4):106–111. nogaster strains. G3 5(4):593–603. de Jong G, Bochdanovits Z. 2003. Latitudinal clines in Drosophila mela- Guirao-Rico S, Gonzalez J. 2019. Evolutionary insights from large scale nogaster: body size, allozyme frequencies, inversion frequencies, and resequencing datasets in Drosophila melanogaster. Curr Opin Insect the insulin-signalling pathway. J Genet. 82(3):207–223. Sci. 31:70–76. Deitz KC, Athrey GA, Jawara M, Overgaard HJ, Matias A, Slotman MA. Guirao-Rico S, Gonzalez J. 2021. Benchmarking the performance of Pool- 2016. Genome-wide divergence in the West-African malaria vector seq SNP callers using simulated and real sequencing data. Mol Ecol Anopheles melas. G3 6(9):2867–2879. Resour. 21(4):1216–1229. dos Santos G, Schroeder AJ, Goodman JL, Strelets VB, Crosby MA, Gutenkunst RN, Hernandez RD, Williamson SH, Bustamante CD. 2009. Thurmond J, Emmert DB, Gelbart WM, FlyBase Consortium 2015. Inferring the joint demographic history of multiple populations from FlyBase: introduction of the Drosophila melanogaster release 6 ref- multidimensional SNP frequency data. PLOS Genet. 5(10):e1000695. erence genome assembly and large-scale migration of genome anno- Hales KG, Korey CA, Larracuente AM, Roberts DM. 2015. Genetics on tations. Nucleic Acids Res. 43(Database issue):D690–D697. the fly: a primer on the Drosophila model system. Genetics Duchen P, Z ivkovic D, Hutter S, Stephan W, Laurent S. 2013. 201(3):815–842. Demographic Inference reveals African and European admixture Haudry A, Laurent S, Kapun M. 2020. Population genomics on the fly: in the North American Drosophila melanogaster population. recent advances in Drosophila. In: Dutheil, JY, editor. Statistical pop- Genetics 193(1):291–301. ulation genomics. Methods in molecular biology. Vol. 2090. New Duffy JB. 2002. GAL4 system in Drosophila: a fly geneticist’s Swiss army York: Humana. p. 357–396. knife. Genesis 34(1–2):1–15. Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. 2005. Very high Durmaz E, Benson C, Kapun M, Schmidt P, Flatt T. 2018. An inversion resolution interpolated climate surfaces for global land areas. Int J supergene in Drosophila underpins latitudinal clines in survival traits. Climatol. 25(15):1965–1978. J Evol Biol. 31(9):1354–1364. Hivert V, Leblois R, Petit EJ, Gautier M, Vitalis R. 2018. Measuring genetic Durmaz E, Kerdaffrec E, Katsianis G, Kapun M, Flatt T. 2020. How selec- differentiation from pool-seq data. Genetics 210(1):315–330. Hoban S, Kelley JL, Lotterhos KE, Antolin MF, Bradburd G, Lowry DB, tion acts on chromosomal inversions. eLS. 1(2):307–315. Available Poss ML, Reed LK, Storfer A, Whitlock MC. 2016. Finding the geno- from: 10.1002/9780470015902.a0028745 mic basis of local adaptation: pitfalls, practical solutions, and future Durmaz E, Rajpurohit S, Betancourt N, Fabian DK, Kapun M, Schmidt P, directions. Am Nat. 188(4):379–397. Flatt T. 2019. A clinal polymorphism in the insulin signaling tran- Hu TT, Eisen MB, Thornton KR, Andolfatto P. 2013. A second- scription factor foxo contributes to life-history adaptation in generation assembly of the Drosophila simulans genome provides Drosophila. Evolution 73(9):1774–1792. new insights into patterns of lineage-specific divergence. Genome Erickson PA, Weller CA, Song DY, Bangerter AS, Schmidt P, Bergland AO. Res. 23(1):89–98. 2020. Unique genetic signatures of local adaptation over space and Hutter S, Li H, Beisswanger S, Lorenzo DD, Stephan W. 2007. Distinctly time for diapause, an ecologically relevant complex trait, in different sex ratios in African and European populations of Drosophila melanogaster. PLoS Genet. 16(11):e1009110. Drosophila melanogaster inferred from chromosomewide single nu- Fabian DK, Kapun M, Nolte V, Kofler R, Schmidt PS, Schlötterer C, Flatt cleotide polymorphism data. Genetics 177(1):469–480. T. 2012. Genome-wide patterns of latitudinal differentiation among Jennings BH. 2011. Drosophila – a versatile model in biology & medicine. populations of Drosophila melanogaster from North America. Mol Mater Today. 14(5):190–195. Ecol. 21(19):4748–4769. Jombart T, Devillard S, Balloux F. 2010. Discriminant analysis of principal Feder AF, Petrov DA, Bergland AO. 2012. LDx: estimation of linkage components: a new method for the analysis of genetically structured disequilibrium from high-throughput pooled resequencing data. populations. BMC Genet. 11:94. PLoS One. 7(11):e48588. Jouganous J, Long W, Ragsdale AP, Gravel S. 2017. Inferring the joint Ferretti L, Ramos-Onsins SE, Perez-Enciso M. 2013. Population genomics demographic history of multiple populations: beyond the diffusion from pool sequencing. Mol Ecol. 22(22):5561–5576. approximation. Genetics 206(3):1549–1567. Flatt T. 2016. Genomics of clinal variation in Drosophila: disentangling Kao JY, Zubair A, Salomon MP, Nuzhdin SV, Campo D. 2015. Population the interactions of selection and demography. Mol Ecol. genomic analysis uncovers African and European admixture in 25(5):1023–1026. Drosophila melanogaster populations from the south-eastern Flatt T. 2020. Life-history evolution and the genetics of fitness compo- United States and Caribbean Islands. Mol Ecol. 24(7):1499–1509. nents in Drosophila melanogaster. Genetics 214(1):3–48. Kapopoulou A, Kapun M, Pieper B, Pavlidis P, Wilches R, Duchen P, Gautier M, Vitalis R, Flori I, Estoup A. 2021. f-statistics estimation and Stephan W, Laurent S. 2020. Demographic analyses of a new sample admixture graph construction with Pool-Seq or allele count data of haploid genomes from a Swedish population of Drosophila mel- using the R package poolfstat. biorxiv. 10.1101/2021.05.28.445945v1. anogaster. Sci Rep. 10(1):22415. Gautier M, Foucaud J, Gharbi K, Cezard T, Galan M, Loiseau A, Thomson Kapun M, Barron MG, Staubach F, Obbard DJ, Wiberg RAW, Vieira J, M, Pudlo P, Kerdelhue C, Estoup A. 2013. Estimation of population Goubert C, Rota-Stabelli O, Kankare M, Bogaerts-Marquez M, et al. allele frequencies from next-generation sequencing data: pool-versus 2020. Genomic analysis of European Drosophila melanogaster pop- individual-based genotyping. Mol Ecol. 22(14):3766–3779. ulations reveals longitudinal structure, continent-wide selection, and Giesen A, Blanckenhorn WU, Sch€afer MA, Shimizu KK, Shimizu-Inatsugi previously unknown DNA viruses. Mol Biol Evol. 37(9):2661–2678. R, Misof B, Niehuis O, Podsiadlowski L, Lischer HEL, Aeschbacher S, et Kapun M, Fabian DK, Goudet J, Flatt T. 2016. Genomic evidence for al. 2020. Genomic signals of admixture and reinforcement between adaptive inversion clines in Drosophila melanogaster. Mol Biol Evol. two closely related species of European sepsid flies. bioRxiv. 33(5):1317–1336. 03.11.985903. Available from: 10.1101/2020.03.11.985903 Kapun M, Flatt T. 2019. The adaptive significance of chromosomal in- Glinka S, Ometto L, Mousset S, Stephan W, Lorenzo D. 2003. version polymorphisms in Drosophila melanogaster. Mol Ecol. Demography and natural selection have shaped genetic variation 28(6):1263–1282. 5803 Kapun et al. . doi:10.1093/molbev/msab259 MBE Kapun M, Schalkwyk H, van McAllister B, Flatt T, Schlötterer C. 2014. Subgroup 2009. The sequence alignment/map format and Inference of chromosomal inversion dynamics from Pool-Seq data in SAMtools. Bioinformatics 25(16):2078–2079. natural and laboratory populations of Drosophila melanogaster. Mol Li H, Stephan W. 2006. Inferring the demographic history and rate of Ecol. 23(7):1813–1827. adaptive substitution in Drosophila. PLoS Genet. 2(10):e166. Keinan A, Clark AG. 2012. Recent explosive human population growth Liao JJZ, Lewis JW. 2000. A note on concordance correlation coefficient. has resulted in an excess of rare genetic variants. Science PDA J Pharm Sci Technol. 54(1):23–26. 336(6082):740–743. Lin LI-K. 1989. a concordance correlation coefficient to evaluate repro- Keightley PD, Ness RW, Halligan DL, Haddrill PR. 2014. Estimation of the ducibility. Biometrics 45(1):255–268. spontaneous mutation rate per nucleotide site in a Drosophila mel- Lynch M, Bost D, Wilson S, Maruki T, Harrison S. 2014. Population- anogaster full-sib family. Genetics 196(1):313–320. genetic inference from pooled-sequencing data. Genome Biol Evol. Keller A. 2007. Drosophila melanogaster’s history as a human commen- 6(5):1210–1218. sal. Curr Biol. 17(3):R77–R81. Machado HE, Bergland AO, O’Brien KR, Behrman EL, Schmidt PS, Petrov Kent WJ, Zweig AS, Barber G, Hinrichs AS, Karolchik D. 2010. BigWig and DA. 2016. Comparative population genomics of latitudinal variation BigBed: enabling browsing of large distributed datasets. in Drosophila simulans and Drosophila melanogaster. Mol Ecol. Bioinformatics 26(17):2204–2207. 25(3):723–740. Koboldt DC, Chen K, Wylie T, Larson DE, McLellan MD, Mardis ER, Machado HE, Bergland AO, Taylor R, Tilk S, Behrman E, Dyer K, Fabian Weinstock GM, Wilson RK, Ding L. 2009. VarScan: variant detection DK, Flatt T, Gonzalez J, Karasov TL, et al. 2021. Broad geographic in massively parallel sequencing of individual and pooled samples. sampling reveals the shared basis and environmental correlates of Bioinformatics 25(17):2283–2285. seasonal adaptation in Drosophila. eLife 10:e67577. Koboldt DC, Zhang Q, Larson DE, Shen D, McLellan MD, Lin L, Miller CA, Mackay TFC, Richards S, Stone EA, Barbadilla A, Ayroles JF, Zhu D, Mardis ER, Ding L, Wilson RK. 2012. VarScan 2: somatic mutation Casillas S, Han Y, Magwire MM, Cridland JM, et al. 2012. and copy number alteration discovery in cancer by exome sequenc- The Drosophila melanogaster genetic reference panel. Nature ing. Genome Res. 22(3):568–576. 482(7384):173–178. Kofler R, Orozco-terWengel P, De Maio N, Pandey RV, Nolte V, Futschik Markow TA, O’Grady PM. 2006. Drosophila: a guide to species identifi- A, Kosiol C, Schlötterer C. 2011. PoPoolation: a toolbox for popula- cation and use. Amsterdam, Boston: Elsevier/Academic Press. tion genetic analysis of next generation sequencing data from Martin M. 2011. Cutadapt removes adapter sequences from high- pooled individuals. PLoS One 6(1):e15925. throughput sequencing reads. EMBNET J. 17(1):10–12. Kofler R, Pandey RV, Schlotterer C. 2011. PoPoolation2: identifying dif- Mateo L, Rech GE, Gonzalez J. 2018. Genome-wide patterns of local ferentiation between populations using sequencing of pooled DNA adaptation in Western European Drosophila melanogaster natural samples (Pool-Seq). Bioinformatics 27(24):3435–3436. populations. Sci Rep. 8(1):16143. Kofler R, Schlötterer C. 2014. A guide for the design of evolve and Mateo L, Ullastres A, Gonzalez J. 2014. A transposable element resequencing studies. Mol Biol Evol. 31(2):474–483. insertion confers xenobiotic resistance in Drosophila. PLoS Genet. Kolaczkowski B, Kern AD, Holloway AK, Begun DJ. 2011. Genomic dif- 10(8):e1004560. ferentiation between temperate and tropical Australian populations McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, of Drosophila melanogaster. Genetics 187(1):245–260. Garimella K, Altshuler D, Gabriel S, Daly M, et al. 2010. The Genome Kreitman M. 1983. Nucleotide polymorphism at the alcohol dehydro- Analysis Toolkit: a MapReduce framework for analyzing next- genase locus of Drosophila melanogaster. Nature 304(5925):412–417. generation DNA sequencing data. Genome Res. 20(9):1297–1303. Kuhn RM, Haussler D, Kent WJ. 2013. The UCSC genome browser and Mölder F, Jablonski KP, Letcher B, Hall MB, Tomkins-Tinch CH, Sochat V, associated tools. Brief Bioinform. 14(2):144–161. Forster J, Lee S, Twardziok SO, Kanitz A, et al. 2021. Sustainable data Lachaise D, Cariou M-L, David JR, Lemeunier F, Tsacas L, Ashburner M. analysis with Snakemake. F1000Res 10:33. 1988. Historical biogeography of the Drosophila melanogaster species Nunez JCB, Paris M, Machado H, Bogaerts M, Gonzalez J, Flatt T, subgroup. In: Hecht MK, Wallace B, Prance GT, editors. Evolutionary Coronado M, Kapun M, Schmidt P, Petrov D, et al. 2021. Note: biology. Boston: Springer US. p. 159–225. updating the metadata of four misidentified samples in the Lack JB, Cardeno CM, Crepeau MW, Taylor W, Corbett-Detig RB, Stevens DrosRTEC dataset. bioRxiv 2021.01.26.428249. KA, Langley CH, Pool JE. 2015. The Drosophila genome nexus: a Nunez JCB, Rong S, Damian-Serrano A, Burley JT, Elyanow RG, Ferranti population genomic resource of 623 Drosophila melanogaster DA, Neil KB, Glenner H, Rosenblad MA, Blomberg A, et al. 2021. genomes, including 197 from a single ancestral range population. Ecological load and balancing selection in circumboreal barnacles. Genetics 199(4):1229–1241. Mol Biol Evol. 38(2):676–685. Lack JB, Lange JD, Tang AD, Russell C-DB, Pool JE. 2016. A thousand fly Ometto L, Glinka S, De Lorenzo D, Stephan W. 2005. Inferring the effects genomes: an expanded Drosophila genome nexus. Mol Biol Evol. of demography and selection on Drosophila melanogaster popula- 33(12):3308–3313. tions from a chromosome-wide scan of DNA variation. Mol Biol Evol. Langley CH, Stevens K, Cardeno C, Lee YCG, Schrider DR, Pool JE, Langley 22(10):2119–2130. SA, Suarez C, Corbett-Detig RB, Kolaczkowski B, et al. 2012. Genomic Orozco-terWengel P, Kapun M, Nolte V, Kofler R, Flatt T, Schlötterer C. variation in natural populations of Drosophila melanogaster. Genetics 2012. Adaptation of Drosophila to a novel laboratory environment 192(2):533–598. reveals temporally heterogeneous trajectories of selected alleles. Mol Laurent SJY, Werzner A, Excoffier L, Stephan W. 2011. Approximate Ecol. 21(20):4931–4941. Bayesian analysis of Drosophila melanogaster polymorphism data Paaby AB, Bergland AO, Behrman EL, Schmidt PS. 2014. A highly pleio- reveals a recent colonization of Southeast Asia. Mol Biol Evol. tropic amino acid polymorphism in the Drosophila insulin receptor 28(7):2041–2051. contributes to life-history adaptation. Evolution 68(12):3395–3409. Lemeunier F, Aulard S. 1992. Inversion polymorphism in Drosophila Paaby AB, Blacket MJ, Hoffmann AA, Schmidt PS. 2010. Identification of melanogaster. In: Krimbas CB, Powell JR, editors. Drosophila inversion a candidate adaptive polymorphism for Drosophila life history by polymorphism. Boka Raton: CRC Press. p. 576. parallel independent clines on two continents. Mol Ecol. Li H. 2011. Tabix: fast retrieval of sequence features from generic TAB- 19(4):760–774. delimited files. Bioinformatics 27(5):718–719. Patterson N, Price AL, Reich D. 2006. Population structure and eigena- Li H. 2013. Aligning sequence reads, clone sequences and assembly nalysis. PLoS Genet. 2(12):e190. contigs with BWA-MEM. arXiv. 1303:3997. Pool JE. 2015. The mosaic ancestry of the Drosophila genetic reference Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, panel and the D. melanogaster reference genome reveals a network Abecasis G, Durbin R, 1000 Genome Project Data Processing of epistatic fitness interactions. Mol Biol Evol. 32:3236–3251. 5804 DEST Population Genomics . doi:10.1093/molbev/msab259 MBE Pool JE, Aquadro C. 2006. History and structure of sub-Saharan popula- Stajich JE, Hahn MW. 2005. Disentangling the effects of demography and tions of Drosophila melanogaster. Genetics 174(2):915–929. selection in human history. Mol Biol Evol. 22(1):63–73. Pool JE, Corbett-Detig RB, Sugino RP, Stevens KA, Cardeno CM, Crepeau Stephan W. 2016. Signatures of positive selection: from selective sweeps MW, Duchen P, Emerson JJ, Saelao P, Begun DJ, et al. 2012. at individual loci to subtle allele frequency changes in polygenic Population genomics of sub-Saharan Drosophila melanogaster: adaptation. Mol Ecol. 25(1):79–88. African diversity and non-African admixture. PLoS Genet. Thia JA, Riginos C. 2019. genomalicious: serving up a smorgasbord of R 8(12):e1003080. functions for population genomic analyses. bioRxiv. 667337. Raineri E, Ferretti L, Esteve-Codina A, Nevado B, Heath S, Perez-Enciso M. Turner TL, Levine MT, Eckert ML, Begun DJ. 2008. Genomic analysis of 2012. SNP calling by sequencing pooled samples. BMC adaptive differentiation in Drosophila melanogaster. Genetics Bioinformatics. 13:239–238. 179(1):455–473. Ramaekers A, Claeys A, Kapun M, Mouchel-Vielh E, Potier D, Weinberger Turner TL, Stewart AD, Fields AT, Rice WR, Tarone AM. 2011. S, Grillenzoni N, Dardalhon-Cumenal D, Yan J, Wolf R, et al. 2019. Population-based resequencing of experimentally evolved popula- Altering the temporal regulation of one transcription factor drives tions reveals the genetic basis of body size variation in Drosophila evolutionary trade-offs between head sensory organs. Dev Cell. melanogaster. PLoS Genet. 7(3):e1001336. 50(6):780–792. Umina PA, Weeks AR, Kearney MR, McKechnie SW, Hoffmann AA. Reinhardt J, Kolaczkowski B, Jones C, Begun D, Kern A. 2014. Parallel 2005. A rapid shift in a classic clinal pattern in Drosophila reflecting geographic variation in Drosophila melanogaster. Genetics climate change. Science 308(5722):691–693. 197(1):361–373. Waldvogel A-M, Feldmeyer B, Rolshausen G, Exposito-Alonso M, Rodrıguez-Trelles F, Tarrıo R, Santos M. 2013. Genome-wide evolution- Rellstab C, Kofler R, Mock T, Schmid K, Schmitt I, Bataillon T, et ary response to a heat wave in Drosophila. Biol Lett. 9(4):20130228. al. 2020. Evolutionary genomics can improve prediction of species’ Rosen Z, Bhaskar A, Roch S, Song YS. 2018. Geometry of the sample responses to climate change. Evol Lett. 4(1):4–18. frequency spectrum and the perils of demographic inference. Wallace MA, Coffman KA, Gilbert C, Ravindran S, Albery GF, Abbott J, Genetics 210(2):665–682. Argyridou E, Bellosta P, Betancourt AJ, Colinet H, et al. 2021. The Rudman SM, Greenblum S, Hughes RC, Rajpurohit S, Kiratli O, Lowder DB, discovery, distribution and diversity of DNA viruses associated with Lemmon SG, Petrov DA, Chaston JM, Schmidt P. 2019. Microbiome Drosophila melanogaster in Europe. Virus Evol. 7(1):veab031. composition shapes rapid genomic adaptation of Drosophila mela- Wittmann MJ, Bergland AO, Feldman MW, Schmidt PS, Petrov DA. nogaster. Proc Natl Acad Sci U S A. 116(40):20025–20032. 2017. Seasonally fluctuating selection can maintain polymorphism Schlötterer C, Tobler R, Kofler R, Nolte V. 2014. Sequencing pools of at many loci via segregation lift. Proc Natl Acad Sci U S A. individuals — mining genome-wide polymorphism data without big 114(46):E9932–E9941. funding. Nat Rev Genet. 15(11):749–763. Wright S. 1943. Isolation by distance. Genetics 28(2):114–138. Schneider D. 2000. Using Drosophila as a model insect. Nat Rev Genet. Zheng X, Gogarten SM, Lawrence M, Stilp A, Conomos MP, Weir BS, Laurie 1(3):218–226. C, Levine D. 2017. SeqArray—a storage-efficient high-performance Sprengelmeyer QD, Mansourian S, Lange JD, Matute DR, Cooper BS, Jirle data format for WGS variant calls. Bioinformatics 33(15):2251–2257. EV, Stensmyr MC, Pool JE. 2020. Recurrent collection of Drosophila Zhu Y, Bergland AO, Gonzalez J, Petrov DA. 2012. Empirical validation of melanogaster from wild African environments and genomic insights pooled whole genome population re-sequencing in Drosophila mel- into species history. Mol Biol Evol. 37(3):627–638. anogaster. PLoS One 7(7):e41901. 5805