Digging into the Genomic Past of Swiss Honey Bees by Whole-Genome Sequencing Museum Specimens

Digging into the Genomic Past of Swiss Honey Bees by Whole-Genome Sequencing Museum Specimens Melanie Parejo 1,2, *, David Wragg 3 , Dora Henriques 4 , Jean-Daniel Charrie`re 1 , and Andone Estonba 2 1 Agroscope, Swiss Bee Research Center, Bern, Switzerland 2 Lab. Genetics, Department of Genetics, Faculty of Science and Technology, University of the Basque Country (UPV/EHU), Leioa, Spain 3 The Roslin Institute, University of Edinburgh, Edinburgh, United Kingdom 4 Instituto Polit ecnico de Braganc¸a, Centro de Investigac¸~ ao de Montanha (CIMO), Braganc¸a, Portugal *Corresponding author: E-mails: [email protected], melanie.par[email protected]. Accepted: 24 August 2020 Data deposition: This project has been deposited at the European Nucleotide Archive (http://www.ebi.ac.uk/ena, last accessed September 5, 2020) under the accession PRJEB34590. Abstract Historical specimens in museum collections provide opportunities to gain insights into the genomic past. For the Western honey bee, Apis mellifera L., this is particularly important because its populations are currently under threat worldwide and have experienced many changes in management and environment over the last century. Using Swiss Apis mellifera mellifera as a case study, our research provides important insights into the genetic diversity of native honey bees prior to the industrial-scale introductions and trade of non-native stocks during the 20th century—the onset of intensive commercial breeding and the decline of wild honey bees following the arrival of Varroa destructor. We sequenced whole-genomes of 22 honey bees from the Natural History Museum in Bern collected in Switzerland, including the oldest A. mellifera sample ever sequenced. We identify both, a historic and a recent migrant, natural or human-mediated, which corroborates with the population history of honey bees in Switzerland. Contrary to what we expected, we find no evidence for a significant genetic bottleneck in Swiss honey bees, and find that genetic diversity is not only maintained, but even slightly increased, most probably due to modern apicultural practices. Finally, we identify signals of selection between historic and modern honey bee populations associated with genes enriched in functions linked to xenobiotics, suggesting a possible selective pressure from the increasing use and diversity of chemicals used in agriculture and apiculture over the last century. Key words: Apis mellifera mellifera, museum genomics, genetic diversity, selection signatures, haplotype phasing, biodiversity. Significance Little is known about native honey bees’ genetic diversity and structure preagricultural and apicultural revolutions during the 20th century—the beginning of commercial bee breeding and decline of wild honey bees following the arrival an invasive ectoparasite. We find no reduction in genetic diversity of a historic honey bee population compared with its contemporary conspecifics. We further identify genes enriched in functions linked to immunity, and the detoxification of possible agrochemicals. The results do not only reveal novel insights into the honey bee genomic past, but also provide valuable baseline genomic data of native populations aiming at making improved conservation management decisions. In addition, our approach to sequence honey bee museum samples serves as a case study for the sequencing of other precious museum specimens. ßThe Author(s) 2020. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution. This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact [email protected] Genome Biol. Evol. 12(12):2535–2551. doi:10.1093/gbe/evaa188 Advance Access publication 2 September 2020 2535 GBE Downloaded from https://academic.oup.com/gbe/article/12/12/2535/5900668 by UNIVERSIDAD DEL PAIS VASCO user on 02 March 2021 Introduction For the most important pollinator of wild and cultured plants, the Western honey bee (Apis mellifera)(Klein et al. 2007; Gallai et al. 2009;IPBES 2016), much has changed in the last decades. Today, it is under pressure globally from new invasive parasites and emergent pathogens, increased use of pesticides, habitat loss, and climate change (Neumann and Carreck 2010;Potts et al. 2010;Vanengelsdorp and Meixner 2010), culminating in major losses of managed honey bee colonies worldwide (Liu et al. 2016;Maggi et al. 2016;Gray et al. 2019;Morawetz et al. 2019). One of the primary factors driving colony losses is the ectoparasitic mite, Varroa destructor, and its associated viruses (Guzm an-Novoa et al. 2010;Dainat et al. 2012). In the late 1970s, this parasite, native to Asia, spread throughout Western Europe and North America decimating wild A. mellifera colonies. Such is the scale of the threat that today the majority of honey bees cannot survive without human intervention (Rosenkranz et al. 2010). It is thus generally accepted that feral honey bee colonies nearly became extinct after the arrival of the mite (Moritz et al. 2007;De la R ua et al. 2009), although there are reports of wild honey bee populations persisting in large woodlands (Seeley 2017;Kohl and Rutschmann 2018). The widespread colony losses and associated population size decline might have potentially resulted in a genetic bottleneck for the remaining wild European honey bees. Such population collapses can lead to loss of genetic diversity and thereby threaten long-term adaptive potential to future environmental changes (Frankham et al. 2002;Allendorfetal. 2013). Honey bees, due to their haplodiploid mating system and social organization, are particularly sensitive to inbreeding depression (Zayed 2009). Studies have demonstrated that high intracolony diversity decreases pathogen load (Desai and Currie 2015) and enhances productivity (Oldroyd et al. 1992;Mattila and Seeley 2007), survivorship (Tarpy et al. 2013), thermoregulation (Jones et al. 2004), and homeostasis (Oldroyd and Fewell 2007). Moreover, it has also been shown that locally adapted honey bees have higher survival (Bu¨chler et al. 2014;Burnham et al. 2019) and lower pathogen levels (Francis et al. 2014), from which follows that there is a need to conserve the underlying genotypic variation (Frankham et al. 2002). In the second half of the 20th century, Europe underwent large-scale agricultural intensification associated with drastic land-use changes, triggering a significant decline in insect biodiversity (Robinson and Sutherland 2002;van Lexmond et al. 2015). One of the major drivers of this decline is the increasing use of pesticides (Le F eon et al. 2010;Goulson et al. 2015). The late 20th century also witnessed an intensification of apiculture with beekeepers beginning to apply chemicals inside the colony to control pests and pathogens (Johnson 2015). Chemicals applied within the hive, such as miticides and antibiotics, as well as agrochemicals acquired externally can persist for many years in beeswax and affect honey bee colonies in the long-term (Mullin et al. 2010).Inthesametime frame, apiculture has further experienced rapid professionalization including migratory beekeeping, increased breeding efforts, and importations of non-native subspecies and selected stock. For instance, throughout large parts of its distributional range the native European M-lineage honey bee subspecies, Apis mellifera mellifera, has been replaced by Clineage bees, mainly Apis mellifera carnica,Apis mellifera ligustica and Buckfast bees preferred by beekeepers (Pinto et al. 2014;Parejo et al. 2016). Past and contemporary populations differ by natural and human-mediated factors and circumstances, such as beekeeping practices, prevailing pathogens, or the pesticide regime on crops. Modern honey bee populations rely heavily on human management, and their genetic composition is therefore influenced by commercial trade (Vanengelsdorp and Meixner 2010) and artificial selection (Wragg et al. 2016; Parejo et al. 2017). Moreover, they are much more exposed to the drastic land-use changes and prevailing agricultural practices of recent times. Throughout much of Europe until the 1950s, honey bees were much less intensively managed, more closely reflecting natural conditions with little or no human selection, and mostly kept by swarm beekeeping, and thereby in constant gene flow with the wild population. Gaining a greater understanding of the genetic diversity in the past can inform our understanding of the impact of the agricultural and apicultural revolutions on honey bee populations. One powerful way to investigate the changes between past and modern populations is by analyzing samples that predate the drastic environmental and human-induced transformations. Museum specimens, therefore, offer an excellent opportunity by providing a window into the past (Lister 2011). Comparing historic and contemporary allelic frequencies is the most direct and powerful way to detect microevolutionary change (Mikheyev et al. 2015). The main caveat of museum samples, however, is the difficulty of obtaining high-quality DNA for molecular genetic analysis (Staats et al. 2013), although improvements in DNA extraction protocols continue to be developed (Tin et al. 2014;Sproul and Maddison 2017). Until recently the majority of studies using museum specimens have been based on PCR-amplification of specific genes or mitochondrial DNA (Habel et al. 2009;La Haye et al. 2012). However, due to DNA degradation, fragments which are shorter than the PCR target region cannot be amplified (Tin et al. 2014). The recent advances in high-throughput sequencing enable us now to overcome the challenges of extracting genomic information from museum specimens, as most methods are designed for short fragmented DNA (Staats et al. 2013;Burrell et al. 2015). In the field of human evolution, protocols for high-throughput sequencing applied to ancient DNA from archaeological sites are well established Parejo et al. GBE 2536 Genome Biol. Evol. 12(12):2535–2551 doi:10.1093/gbe/evaa188 Advance Access publication 2 September 2020 Downloaded from https://academic.oup.com/gbe/article/12/12/2535/5900668 by UNIVERSIDAD DEL PAIS VASCO user on 02 March 2021 (reviewed by Slatkin and Racimo 2016). However, fewer efforts have been made in the application to historical museum collections from animals and plants (e.g., sequencing of the mitogenome as in Miller et al. 2008;Hung et al. 2013; Besnard et al. 2016), and only most recently high-density genome-wide analyses of museum specimens have been reported (Mikheyev et al. 2015;Linck et al. 2017;S anchez Barreiro et al. 2017;Cridland et al. 2018). To date, for honey bees, there have been two studies published using whole-genome sequence data of museum specimen, both of which concerned the introduced range of A. mellifera in North America. Mikheyev et al.(2015) investigated the genomic changes over a 33-year period of a wild honey bee population in Ithaca following the introduction of V. destructor. The authors found evidence of a mitochondrial bottleneck but with little loss of nuclear genetic diversity and population size. Cridland et al. (2018) documented the temporal genetic changes in Californian populations with northern populations experiencing a shift in genetic ancestry from Mto C-lineage since the 1960s and southern populations undergoing Africanization. Here, we hypothesize that the drastic changes regarding population decline, agricultural and apicultural intensification during the last decades have had profound effects on the genetic diversity and ancestry of native honey bee populations and left signatures of selection on their genomes. Using the Swiss dark honey bee population as a case study, we sequenced museum specimens dated between 1879 and 1959 that were provided by the Natural History Museum in Bern, Switzerland, to investigate the genomic past of the native A. m. mellifera bees. Native Swiss honey bees have suffered from a severe population size decline in recent decades due, in part, to Varroa mites, as well as introductions and replacements of non-native stocks. Together with wholegenome sequence data from a previous study of contemporary Swiss bees (Parejo et al. 2016), we investigated genetic diversity, mitochondrial DNA haplotypes, admixture, and selection signatures. To our knowledge, this is the first study to whole-genome sequence historic honey bee specimen from their native range. Materials and Methods Samples Museum Samples Twenty-two Swiss A. m. mellifera museum specimens dated between 1879 and 1959 (61–141 years old) were obtained from the Natural History Museum in Bern, Switzerland (fig. 1). Specimen consisted of dried and pinned worker bees (diploid) stored by the museum but originating from several private collections of Swiss entomologists. Samples have been assigned with a QR-code deposited in the Natural History Museum in Bern (supplementary table S1,Supplementary Material online). Modern Samples Whole-genome sequence data were available from a previous study (Parejo et al. 2016) of which we selected 40 pure Swiss A. m. mellifera drones (haploid). These samples cover a slightly larger geographical range than the available museum bees, but due to breeding activities they are very similar to each other representing the current population. In addition, to investigate overall genetic structure and admixture proportions, 36 honey bees from a different evolutionary lineage widely employed by beekeepers in Switzerland were included in some analyses. These included 24 A. m. carnica and 12 A. m. ligustica drones (supplementary table S1,Supplementary Material online) from Parejo et al. (2016) and Henriques et al. (2018). DNA Extraction and Sequencing Genomic DNA was extracted from the hind legs of museum specimen (fig. 1) carefully rinsed with Ringer solution, using a phenol–chloroform–isoamyl alcohol (25:24:1) method (Ausubel 1988). Pair-end (2 125 bp) libraries (kit) were prepared following manufacturers protocol using the NEBNext Ultra II kit (New England Biolabs, Inc) and sequenced on the Illumina HiSeq3000 platform with 20 samples per lane. Mapping, Variant Calling, and Single Nucleotide Polymorphism Sets Mapping Raw sequence data from modern and historic samples were processed using Cutadapt v1.8 (Martin 2011)toremove Illumina universal adaptors and keep only reads with minimum lengths of 20 bp and a minimum base quality score of 20. Trimmed reads were then mapped against the dark honey bee reference genome INRA_AMelMel_1.0 (www.ncbi. nlm.nih.gov/assembly/GCA_003314205.1, last accessed September 5, 2020) using bwa mem 0.7.10 (Li and Durbin 2009). PCR duplicates were marked using PICARD 2.18.23 (http://broadinstitute.github.io/picard/, last accessed September 5, 2020). Mapping statistics including depth of coverage and percentage of mapped reads were calculated using samtools 1.7 (Li et al. 2009) and GATK v4.1.0.0 (Mckenna et al. 2010;Van Der Auwera et al. 2013). DamageProfiler (https://damageprofiler.readthedocs.io/en/ latest/index.html, last accessed September 5, 2020) was used on museum samples to generate damage profiles of the mapped DNA reads caused by deamination of cytosine over time which leads to misincorporations of G!Aatthe5 0 and C!Tatthe3 0ends (Briggs 2010;Sawyer et al. 2012). Whole-Genome Sequencing Museum Specimens GBE Genome Biol. Evol. 12(12):2535–2551 doi:10.1093/gbe/evaa188 Advance Access publication 2 September 2020 2537 Downloaded from https://academic.oup.com/gbe/article/12/12/2535/5900668 by UNIVERSIDAD DEL PAIS VASCO user on 02 March 2021 Variant Calling To increase variant confidence, single nucleotide polymorphism (SNP) calling for the A. m. mellifera (N MUSEUM ¼22 workers, N MODERN ¼40 drones) was performed using two different software tools (GATK’s Haplotypecaller and SAMtools mpileup). First, following GATK’s best practices, individual GVCFs were produced using Haplotypecaller with parameters: minimum mapping quality ¼20, max alternate alleles ¼2, minimum quality score ¼20, and sampleploidy ¼2 for museum (diploid workers) and sampleploidy ¼1 for modern samples (haploids drones). Subsequently, GVCFs were combined and genotyped to produce a VCF-file containing raw variants for all A. m. mellifera samples. Variants were filtered according GATK’s hardfiltering recommendation (MQ <40.0, FS >60.0, QUAL < 30.0, MQRankSum <12.5, ReadPosRankSum <8.0) except for quality by depth (QD), where a stricter filter was applied (QD >5) after QD distribution was investigated before and after filtering (supplementary figs. S1 and S2, Supplementary Material online). Second, multisample SNP calling was also performed using SAMtools/BCFtools mpileup 1.7 (Li et al. 2009) with parameters mapQuality (q>30), baseQuality (Q>20), and filtering low-quality variants (QUAL <30). Both call sets (GATK and SAMtools) were filtered on depth (minimum 5, maximum 3*average DP) and to include only biallelic SNPs on chromosomes 1–16. Finally, the variants from both call sets were merged using BCFtools isec to keep only SNPs identified in both sets. Variant calling statistics were calculated with BCFtools stats. Annotation No published annotation is available for the dark honey bee reference genome (A. m. mellifera; INRA_AMelMel_1.0). Thus, annotation files (gtf, annotation release 104) from the latest honey bee genome Amel_HAv3.1 (Wallberg et al. 2019) were remapped onto the INRA_AMelMel_1.0 genome using NCBI’s remapping service (www.ncbi.nlm.nih.gov/genome/tools/remap, last accessed September 5, 2020). Finally, a custom database for SnpEff4.3t (Cingolani et al. 2012) as per software instructions was generated to annotate the variants and predict their potential effects excluding intergenic, upand downstream annotations. Haplotype Phasing Phasing genotypes into haplotypes is a fundamental requirement of some analyses, such as that of extended haplotype homozygosity (see below), which seek to exploit linkage disequilibrium (LD) between markers. Statistical phasing can be performed with or without a reference panel, and generally the use of an external reference panel has been shown to increase phasing accuracy (Delaneau et al. 2013). However, using haplotypes from the modern drone (haploid) data set to phase the historic worker (diploid) data set risks potentially Fig. 1.—Sampling sites of the 22 A. m. mellifera museum specimens. Most samples originate from the region around Bern dating between 1941 and 1959, but some are from mountain areas. The oldest sample is from Luzern (1879), Central Switzerland. The second oldest sample (1884) is from Zermatt, Valais, in the Southern Alps. Map created with Datawrapper (www.datawrapper.de, accessed February 2020). Parejo et al. GBE 2538 Genome Biol. Evol. 12(12):2535–2551 doi:10.1093/gbe/evaa188 Advance Access publication 2 September 2020 Downloaded from https://academic.oup.com/gbe/article/12/12/2535/5900668 by UNIVERSIDAD DEL PAIS VASCO user on 02 March 2021 creating artifactual museum haplotypes that are modern recombinants of ancestral variation. We therefore tested the impact of phasing both with and without the use of the modern bees as a reference panel on the museum bees. First, SNPs for all samples were filtered on call rate >0.95 using VCFtools (Danecek et al. 2011) leaving 2,651,904 SNPs and missing SNPs were imputed within SHAPEIT4 (Delaneau et al. 2019). Preliminary analyses masking 5% of the SNPs with 100% call rate revealed that imputation using this approach is highly accurate (>90% accuracy, supplementary fig. S6,Supplementary Material online) while keeping a larger number of SNPs. Subsequently, the data were phased without a reference panel in museum and modern data sets independently, and using the self-imputed modern bees as references we phased the unphased museum data once more. Phasing was performed with SHAPEIT4 (Delaneau et al. 2019) with the sequencing flag, a minimum window size of 0.1 Mb (the minimum permitted by SHAPEIT4) and an effective population (Ne) size of 150,000 which approximates the Ne calculation of Wallberg et al. (2014) for the Northern A. m. mellifera samples they studied. Genetic maps for phasing were generated by lifting over the Amel4.5 reference genome physical positions in the crossover data generated by Liu et al. (2015) to those of INRA_AMelMel_1.0, and used the crossover information to estimate genetic positions for each SNP. The combined phased data sets including either reference or self-phased museum worker bees were filtered on minor allele frequency (MAF >0.05), each leaving 1,828,439 SNPs for signatures of selection analyses and linkage decay estimation. Mitochondrial DNA Variants in the mitochondrial genome of the A. m. mellifera bees were called using SAMtools mpileup 1.7 with the— ploidy 1 option and keeping only SNP variants with high quality (QUAL >30). The resulting 254 SNPs were genotyped in the C-lineage samples using GATK’s Haplotypecaller with mode genotype-given-alleles. C-lineage and A. m. mellifera samples were combined and SNPs with >20% missing calls were removed. This left 205 SNPs to perform Median-Joining network analysis (Bandelt et al. 1999) in PopART (Leigh and Bryant 2015). Population Structure Analyses Combining Haploid Drones to Diploid Individuals Population structure analyses are sensitive when haploid and diploid data sets are analyzed together (Dufresne et al. 2014; Wragg et al. 2016). To get the most informative and unbiased results for population structure analyses, we therefore randomly combined the haploid genotypes of two modern A. m. mellifera drones using a custom script to generate in silico diploids following the procedure in Wragg et al. (2016). For population structure analyses, the data set with the diploid museum samples and “diploidized” modern samples of A. m. mellifera was filtered and pruned on LD using PLINK1.9 (Chang et al. 2015). We applied –indep-pairwise 50 10 0.1 to filter out variants with a correlation of >0.1 in each window of 50 variants as recommended by Alexander et al. (2009) and kept only SNPs with a 100% call rate. This left 59 K independent SNPs, which were then genotyped in the C-lineage drones using GATK’s Haplotypecaller with mode genotype-given-alleles. Finally, similar to the modern A. m. mellifera samples, the haploid C-lineage drones were randomly combined into diploid individuals using a custom script. The combined data set for the population structure analyses comprised 59 K SNPs genotyped in 60 samples (22 museum A. m. mellifera, 20 “diploidized” modern A. m. mellifera, 12 “diploidized” A. m. carnica, and 6 “diploidized” A. m. ligustica). This data set was used to estimate the average genome-wide divergence, model-based ancestry, and in principal component analysis (PCA). Ancestry, PCA, and Population Differentiation To infer the genetic ancestry of each individual, we performed model-based clustering as implemented in ADMIXTURE (Alexander et al. 2009). We ran the analysis unsupervised with 10,000 iterations for 1–5 hypothetical ancestral (K) clusters. Cross-validation error was estimated for each cluster and used to determine the optimal number of K clusters. We also performed PCA to assess the population structure in the absence of a model (Price et al. 2006). PCA was applied to the pairwise genetic relationships between all individuals (N¼60) according to their identity-by-state values computed in PLINK 1.9 (Chang et al. 2015). Admixture and PCA results were processed and plotted in R (R Development Core Team 2013). Based on these results, a single admixed museum bee was identified, which was subsequently excluded from downstream analyses (population differentiation, genetic diversity, linkage, and selection signatures). Population differentiation was estimated as mean pairwise F ST (Weir and Cockerham 1984) per site as implemented in VCFtools (Danecek et al. 2011). The mean and confidence intervals were calculated from 10 randomly selected bootstrap samples of 10 modern A. m. mellifera,10historicA. m. mellifera, and 10 C-lineage bees. Genetic Diversity To infer the adaptive potential within the modern and historic A. m. mellifera populations two genetic diversity measures were employed: 1) Expected heterozygosity (H Exp ) for each individual was calculated using the data set of 59 K unlinked SNPs; and 2) nucleotide diversity (p;Nei 1982), which was calculated from the whole-genome data for each population Whole-Genome Sequencing Museum Specimens GBE Genome Biol. Evol. 12(12):2535–2551 doi:10.1093/gbe/evaa188 Advance Access publication 2 September 2020 2539 Downloaded from https://academic.oup.com/gbe/article/12/12/2535/5900668 by UNIVERSIDAD DEL PAIS VASCO user on 02 March 2021 in window sizes of 10 kb with 5 kb overlap using VCFtools (Danecek et al. 2011). The mean and confidence intervals for H Exp and pwere calculated with 10 randomly selected samplesforeachof10modernand10historicA. m. mellifera. LD Estimation LD between each pair of SNPs in the modern and historic populations was calculated using the selfand referencephased museum bees (N¼21 samples !N¼42 haplotypes) and the haploid modern drones (N¼40) by estimating the Pearson’s squared correlation coefficient (r 2 ) in Plink v1.09 (Chang et al. 2015). Pairwise LD (as r 2 values) were calculated between all the pairs of SNPs within each chromosome based on the exact solution of the Hill equation (Hill 1974;Gaunt et al. 2007) and applying a MAF filter of 0.05. The extent of LD decay was estimated based on physical distance between each SNP pair. LD decay curves were calculated as the average r 2 within bins of 200 base pairs, up to a distance of 10 kbp and plotted in R. Selection Signatures Analyses Genome Scan Several measures are available to infer signatures of selection, employing a range of test statistics each with its own limits and merits (reviewed in Vitti et al. 2013). Here, we employed cross-population extended haplotype homozygosity (XP-EHH) described by Sabeti et al. (2007). The XP-EHH test statistic is basedonextendedhaplotypelengthwhichhasbeenshown to be most suitable to infer recent selection and is especially useful for identifying hard and soft selective sweeps with causative alleles that have not reached fixation (Sabeti et al. 2007;Vitti et al. 2013). The XP-EHH analysis is based on haplotype length and is, thus, sensitive to accurate phasing. Hence, we performed the analysis twice—once comparing the self-phased museum to modern bees, and secondly comparing museum bees phased with the modern bees as a reference panel to the modern bees. The analyses were performed using selscan v1.2.1 (Szpiech and Hernandez 2014) with default MAF and EHH truncation values of 0.05. The obtained XP-EHH scores for each SNP were normalized (Z-transformed) by subtracting genome-wide mean XP-EHH and dividing by the standard deviation (supplementary fig. S4, Supplementary Material online). SNPs with absolute Z-transformed XP-EHH values in the 99th percentile were considered significant. A similar approach has been employed by other studies investigating selection signatures in the honey bee genome (Wallberg et al. 2016;Montero-Mendieta et al. 2019). Finally, genes associated with SNPs annotated with SNPeff as intronic, exonic, 30UTR, 50UTR, and splice variants, that were identified in both genome scans (self-phased and reference-phased worker bees) are considered as putative candidates under selection. Gene Ontology Enrichment Analysis Gene Ontology (GO) annotations provide a convenient means of grouping genes by their known functions and predicted biological roles, enabling enrichment analyses to be conducted. The online resource DAVID v.6.8 (Database for Annotation, Visualization, and Integrated Discovery; Huang et al. 2007) was accessed in June 2020 to test if the candidate genes identified demonstrated enrichment for any particular function (Huang et al. 2009). We used as a background gene set all honey bee genes associated with at least one SNP in our analyses. Functionally related genes were clustered using the gene functional classification tool set to highest stringency. Enrichment for GO category terms was performed with the functional annotation analysis tool using the GO categories of Biological Process, Cellular Component, Molecular Function, and KEGG pathway. The functional annotation clustering tool was subsequently used to cluster similar GO terms. Results Mapping and SNP Calling In total, 686,376,114 sequencing reads were generated from the 22 museum samples. A summary of alignment statistics of these sequence reads in addition to those of the 76 modern samplesisprovidedinsupplementary table S1, and figures S1 and S2, Supplementary Material online. The average mapping rate across the museum samples was 93.7%, with a mean depth of coverage of 13.9, ranging from 5.55to 27.39, in comparison with the modern A. m. mellifera drones with mean 10.3and range 7.3–21.2. Only 34% of reads for sample KirBE_1941 mapped to the reference genome, possibly indicating contamination, nevertheless, the depth of coverage was 7.8and the sample retained in downstream analyses. On average for museum and modern samples, respectively, 85% and 92.2% of the genome per sample was callable, having a depth of coverage >4, whereas the lowest breadth of coverage was observed in BerBE_1947-2 (57.7%). The oldest sample, LuzLU_1879-2, returned a depth and breadth of coverage of 9.02and 69.6%, respectively. Analysis of DNA degradation by DamageProfiler (supplementary table S2 and fig. S5,Supplementary Material online) indicated only minor 50(2.1% 60.6SD)and3 0(2.4% 60.5 SD) misincorporations compared with studies of ancient DNA (e.g., 8%, Peltzeretal.2018). The misincorporations of two the oldest museum samples dating from 1884 to 1879 were estimated at 3% and 4.5% for G!Aatthe5 0end, whereas for the C!Tatthe3 0end it was 2.5% and 4.3%, respectively. However, because mismatches at the ends of the sequence reads are soft-clipped by bwa mem during alignment, no additional read filtering or clipping was performed. Overall, with except for the two oldest samples and BerBE_1947-2 with the lowest breath of coverage, the sequence quality of the museum samples can be considered Parejo et al. GBE 2540 Genome Biol. Evol. 12(12):2535–2551 doi:10.1093/gbe/evaa188 Advance Access publication 2 September 2020 Downloaded from https://academic.oup.com/gbe/article/12/12/2535/5900668 by UNIVERSIDAD DEL PAIS VASCO user on 02 March 2021 comparable with those of the modern samples (supplementary tables S1 and S2,Supplementary Material online). Variant calling was performed using GATK’s Haplotypecaller and SAMtools mpileup, resulting in 4,611,541 and 4,145,606 SNPs, respectively (supplementary table S3,Supplementary Material online). Applying filters to remove low-quality variants and intersecting both variant sets, left 3,252,197 genome-wide SNPs. The average genotyping call rate was 0.94 in the museum samples and 0.97 in the modern bees, with all samples except for three exceeding 90% call rate (supplementary table S3A,figs. S6 and S7, Supplementary Material online). The museum sample with the lowest genotyping rate of 0.78 was BerBE_1947-2, which was also the sample with the lowest mapping rate potentially due to contamination. Different additional filters were applied for subsequent analyses (supplementary table S3,Supplementary Material online): For estimating nucleotide diversity and for phasing the data set was additionally filtered on call rate 0.95 leaving 2,651,904 SNPs, for linkage and XP-EHH a MAF filter of >0.05 was applied leaving 1,828,439 SNPs, and for PCA, model-based admixture, F ST ,andH Exp the data set was filtered for unlinked SNPs with a call rate of 100% leaving 59,320 SNPs. Mitochondrial Haplotype Networks Mitochondrial network analysis revealed evidence of mitonuclear discordance. The native dark honey bee of Switzerland belongs to the M-lineage, but we found that one individual sampled in Liebefeld in 1959 (LieBE_1959-1) and another sampled in the Swiss Alps in 1958 (LoeVS_1958) possess mitochondrial haplotypes that cluster with C-lineage bees (fig. 2). A phylogenetic analysis of SNPs identified on the complete mitochondrial genome showed that all but two modern and historic A. m. mellifera samples cluster in two M-lineage clades (supplementary fig. S8,Supplementary Material online). The two historic samples that are placed in a distant branch in the phylogenetic analysis are the same samples that cluster within the C-linage of the haplotype network. Population Structure Population structure inferred by the model-based ADMIXTURE and PCA each revealed four (sub-) populations representing historic A. m. mellifera, modern A. m. mellifera,A. m. carnica and A. m. ligustica bees (fig. 3). The optimal number of clusters identified from the lowest cross-validation error is K¼2(supplementary fig. S9,Supplementary Material online), which correspond to the two major evolutionary lineages: M- (A. m. mellifera) and C-lineage (A. m. carnica and A. m. ligustica). One of the museum bees, LieBE_1959-1 shows evidence of genetic admixture with the C-lineage (60.5% M-lineage and 39.5% C-lineage ancestry), as indicated by the dual-color genetic background in the ADMIXTURE plot (fig. 3A)andits intermediate placement along PC1 in the PCA plot (fig. 3B). LieBE_1959-1 is one of two bees identified from the mitochondrial analyses to possess C-lineage mtDNA (fig. 2). When estimating population differentiation (F ST )usingthe 59 K SNP data set, the museum bee identified as being admixed was excluded from the analysis. The F ST analysis revealed a high divergence between C-lineage bees and the historic and modern A. m. mellifera populations (F ST >0.4), and, whereas the divergence between the modern and museum samples was expectedly very low, it was significantly different from 0 based on random subsampling (F ST ¼ 0.007, 95% CI 0.004–0.009). Genetic Diversity To investigate differences in genetic diversity and population histories between modern and museum A. m. mellifera,we calculated expected heterozygosity (H Exp ) and nucleotide diversity (p). Heterozygosity differed significantly between both populations (H E(Museum) ¼0.243, 95% CI 0.241–0.245, H E(Modern) ¼0.258, 95% CI 0. 257–0.259). The formula of H Exp is based on allele frequencies and is therefore not influenced by ploidy, nor the generation of in silico diploids (as we tested in preliminary analyses). Moreover, we also estimated nucleotide diversity which is calculated on haploid sequences (phased genotypes) and, thus, insensitive to ploidy. Similarly to H Exp also nucleotide diversity was slightly higher in the modern A. m. mellifera population (p¼0.00241, 95% CI 0.00240–0.00243) compared with the historic (p¼0.00227, 95% CI 0.00224–0.00230). Both measures of genetic diversity being larger in the modern population suggest enhanced adaptive potential. Observed heterozygosity in the identified admixed museum bee (LieBE_1959-1) was considerably higher (H Obs ¼0.42) than the expected heterozygosity (H Exp ¼0.26), and also higher than for all other A. m. mellifera samples, suggesting the admixture to be recent. LD Decay LD decay between SNPs in the modern and historic A. m. mellifera populations as measured using r 2 over increasing distances between pairwise SNPs is shown in figure 4.The maximum average LD for SNPs less than 200 bps apart was 1.5 times as high in the modern population (r 2 ¼0.43) compared with the historic population (r 2 0.28). LD decays quickly for both populations, but long-range LD was found to be considerably lower for the historic (r 2 0.02) than the modern population (r 2 0.12) potentially reflecting the recent population history of a small, inbred or admixed population. There is a slight, but insignificant tendency for LD to be higher in museum bees phased with the drones as references (fig. 4). Whole-Genome Sequencing Museum Specimens GBE Genome Biol. Evol. 12(12):2535–2551 doi:10.1093/gbe/evaa188 Advance Access publication 2 September 2020 2541 Downloaded from https://academic.oup.com/gbe/article/12/12/2535/5900668 by UNIVERSIDAD DEL PAIS VASCO user on 02 March 2021 Selection Signatures between Historic and Modern A. m. mellifera Populations We investigated the presence of signatures of selection between the historic and modern A. m. mellifera populations using the XP-EHH method (Sabeti et al. 2007) and with two different data sets including either the self-phased or reference-phased museum worker haplotypes. There was an 84% overlap of the SNPs (15,354 SNPs) falling within the 99th percentile in both scans (supplementary fig. S10, Fig. 2.—Median-joining network inferred from 205 mtDNA SNPs and 60 samples (N¼22 museum A. m. mellifera,N¼20 modern A. m. mellifera, N¼12 A. m. carnica,andN¼6A. m. ligustica). Hypothetical (unsampled or extinct) haplotypes are denoted as filled black circles. The values in brackets indicate base pair differences between haplotypes. M-lineage samples including modern and historic A. m. mellifera are grouped into two clades, with the exception of two museum samples (LieBE_1959-1 and LoeVS_1958) which cluster with C-lineage bees (denoted by the arrows). Fig. 3.—Population structure inferred from the LD-pruned 59 K SNPs and 60 samples (N¼22 museum A. m. mellifera,N¼20 “diploidized” modern A. m. mellifera,N¼12 “diploidized” A. m. carnica,andN¼6 “diploidized” A. m. ligustica). (A) Genetic ancestry as calculated with ADMIXTURE for K¼2to4 hypothetical ancestral populations. Each color represents one of Kclusters. Each individual is represented by a horizontal bar and colored according to the proportion of the genome that was derived from each cluster. The optimal number of clusters identified by cross-validation is K¼2. (B) PCA of genetic distance between individuals. The first principal component (PC1) explains 97% of the variation indicating strong divergence between Mand C-lineage honey bees, whereas PC2 accounts only for 0.2% of the variance. Parejo et al. GBE 2542 Genome Biol. Evol. 12(12):2535–2551 doi:10.1093/gbe/evaa188 Advance Access publication 2 September 2020 Downloaded from https://academic.oup.com/gbe/article/12/12/2535/5900668 by UNIVERSIDAD DEL PAIS VASCO user on 02 March 2021 Supplementary Material online) which were considered as being associated with putative signatures of selection (fig. 5 and supplementary table S4,Supplementary Material online). In a first screening, we identified an extreme peak on chromosome 11 (supplementary fig. S11,Supplementary Material online), which on further investigation was identified to most likely be a duplication not captured in the reference genome due to 1) the average depth of coverage in the region being 1.6 times the chromosome average, and 2) the presence of reads with both alleles (reference and alternate) in the haploid drones. Furthermore, this peak lies within a lowcomplexity centromere and no gene is annotated within 615 kb. We have thus excluded it from further analyses. The remaining SNPs are associated with 644 candidate genes (supplementary table S5,Supplementary Material online). Many of the genes identified are uncharacterized loci (29.8%) as the annotation of the honey bee genome is far from complete. Of the 15,354 SNPs, 275 and 13 SNPs were predicted by SNPeff to have MODERATE and HIGH impact, respectively (supplementary table S6,Supplementary Material online). These included mostly nonsynonymous base pair changes, and also splice region variants as well as stop lost/ gained variants. GO analysis is a strategy to identify the most important biological processes of candidate regions identified in whole-genome selection scans. However, the power of these analyses depends on the number and quality of annotated genes available for the focal organism (Yon Rhee et al. 2008). We conducted a GO analysis using the online resource Fig. 4.—LD between SNPs as measured by r 2 (yaxis) for increasing distance between SNPs (xaxis) for A. m. mellifera modern drones (N¼40) and A. m. mellifera museum bee haplotypes (N¼42 haplotypes) selfphased and phased using the drones as a reference panel. Fig. 5.—Signatures of selection between historic and modern A. m. mellifera from Switzerland. XP-EHH was performed using 42 haplotypes derived from 21 museum samples (diploid) (A) self-phased and (B) reference-phased, and 40 haplotypes derived from modern drones (haploid). XP-EHH scores are plotted along the 16 honey bee chromosomes with negative values indicating selection in the modern population. The dashed lines denote SNPs in the 99th percentile of the absolute XP-EHH scores. This figure excludes the false positive peak on chromosome 11 (4945317-4945798), which can be seen in figure S11,Supplementary Material online. The five highest peaks of each analysis are labeled with their putative genes under selection. Whole-Genome Sequencing Museum Specimens GBE Genome Biol. Evol. 12(12):2535–2551 doi:10.1093/gbe/evaa188 Advance Access publication 2 September 2020 2543 Downloaded from https://academic.oup.com/gbe/article/12/12/2535/5900668 by UNIVERSIDAD DEL PAIS VASCO user on 02 March 2021 Bonn, Germany: Secretariat of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services. Johnson RM. 2015. Honey bee toxicology. Annu Rev Entomol. 60(1):415–434. Jones JC, Myerscough MR, Graham S, Oldroyd BP. 2004. Honey bee nest thermoregulation: diversity promotes stability. Science 305(5682):402–404. Klein A-M, et al. 2007. Importance of pollinators in changing landscapes for world crops. Proc R Soc B 274(1608):303–313. Vol. p Kohl PL, Rutschmann B. 2018. The neglected bee trees: European beech forests as a home for feral honey bee colonies. PeerJ 6:e4602. La Haye MJJ, Neumann K, Koelewijn HP. 2012. Strong decline of gene diversity in local populations of the highly endangered Common hamster (Cricetus cricetus) in the western part of its European range. Conserv Genet. 13(2):311–322. Le F eon V, et al. 2010. Intensification of agriculture, landscape composition and wild bee communities: a large scale study in four European countries. Agric Ecosyst Environ. 137(1–2):143–150 Leigh JW, Bryant D. 2015. PopART: full-feature software for haplotype network construction. Methods Ecol Evol. 6(9):1110–1116. Lemmon MA, Schlessinger J. 2010. Cell signaling by receptor tyrosine kinases. Cell 141(7):1117–1134. Li H, et al. 2009. The sequence alignment/map format and SAMtools. Bioinformatics 25(16):2078–2079. Li H, Durbin R. 2009. Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics 25(14):1754–1760. Linck EB, Hanna ZR, Sellas A, Dumbacher JP. 2017. Evaluating hybridization capture with RAD probes as a tool for museum genomics with historical bird specimens. Ecol Evol. 7(13):4755–4767. Lister AM. 2011. Natural history collections as sources of long-term datasets. Trends Ecol Evol. 26(4):153–154. Liu H, et al. 2015. Causes and consequences of crossing-over evidenced via a high-resolution recombinational landscape of the honey bee. Genome Biol. 16(1):15. Liu Z, et al. 2016. Survey results of honey bee (Apis mellifera) colony losses in China (2010–2013). J Apic Res. 55(1):29–37. Maggi M, et al. 2016. Honeybee health in South America. Apidologie 47(6):835–854. Martin M. 2011. Cutadapt removes adapter sequences from highthroughput sequencing reads. EMBnet J. 17(1):10–12. Mattila HR, Seeley TD. 2007. Genetic diversity in honey bee colonies enhances productivity and fitness. Science 317(5836):362–364. Mckenna A, et al. 2010. The genome analysis toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 20(9):1297–1303. Miguel I, et al. 2011. Both geometric morphometric and microsatellite data consistently support the differentiation of the Apis mellifera M evolutionary branch. Apidologie 42(2):150–161. Mikheyev AS, Tin MM, Arora J, Seeley TD. 2015. Museum samples reveal rapid evolution by wild honey bees exposed to a novel parasite. Nat Commun. 6:7991. Miller W, et al. 2008. The mitochondrial genome sequence of the Tasmanian tiger (Thylacinus cynocephalus). Genome Res. 19(2):213–220. Montero-Mendieta S, et al. 2019. The genomic basis of adaptation to high-altitude habitats in the eastern honey bee (Apis cerana). Mol Ecol. 28(4):746–760. Morawetz L, et al. 2019. Health status of honey bee colonies (Apis mellifera) and disease-related risk factors for colony losses in Austria. PLoS One 14(7):e0219293. Moreau-Fauvarque C, Taillebourg E, Pr eat T, Dura J-M. 2002. Mutation of linotte causes behavioral defects independently of pigeon in Drosophila. NeuroReport 13(17):2309–2312. Moritz RF, Kraus FB, Kryger P, Crewe RM. 2007. The size of wild honeybee populations (Apis mellifera) and its implications for the conservation of honeybees. J Insect Conserv. 11(4):391–397. Mullin CA, et al. 2010. High levels of miticides and agrochemicals in North American apiaries: implications for honey bee health. PLoS One 5(3):e9754. Nei M. 1982. Evolution of human races at the gene level. Hum Genet A 167–181. Neumann P, Carreck NL. 2010. Honey bee colony losses. J Apic Res. 49(1):1–6. Ng TH, Kurtz J. 2020. Dscam in immunity: a question of diversity in insects and crustaceans. Dev Comp Immunol. 105:103539. Oldroyd BP, Fewell J. 2007. Genetic diversity promotes homeostasis in insect colonies. Trends Ecol Evol. 22(8):408–413. Oldroyd BP, Rinderer T, Harbo J, Buco S. 1992. Effects of intracolonial genetic diversity on honey bee (Hymenoptera: Apidae) colony performance. Ann Entomol Soc Am. 85(3):335–343. Oxley PR, Spivak M, Oldroyd BP. 2010. Six quantitative trait loci influence task thresholds for hygienic behaviour in honeybees (Apis mellifera). Mol Ecol. 19(7):1452–1461. Parejo M, et al. 2016. Using whole-genome sequence information to foster conservation efforts for the European Dark Honey Bee, Apis mellifera mellifera. Front Ecol Evol. 4:140. doi:10.3389/ fevo.2016.00140. Parejo M, Wragg D, Henriques D, Vignal A, Neuditschko M. 2017. Genome-wide scans between two honeybee populations reveal putative signatures of human-mediated selection. Anim Genet. 48(6):704–707. Peltzer A, et al. 2018. Inferring genetic origins and phenotypic traits of George B€ ahr, the architect of the Dresden Frauenkirche. Sci Rep. 8(1). doi:10.1038/s41598-018-20180-z Pinto MA, et al. 2014. Genetic integrity of the Dark European honey bee (Apis mellifera mellifera) from protected populations: a genome-wide assessment using SNPs and mtDNA sequence data. J Apic Res. 53(2):269–278. Potts SG, et al. 2010. Global pollinator declines: trends, impacts and drivers. Trends Ecol Evol. 25(6):345–353. Price AL, et al. 2006. Principal components analysis corrects for stratification in genome-wide association studies. Nat Genet. 38(8):904–909. Priestley CM, Williamson EM, Wafford KA, Sattelle DB. 2003. Thymol, a constituent of thyme essential oil, is a positive allosteric modulator of human GABAA receptors and a homo-oligomeric GABA receptor from Drosophila melanogaster. Br J Pharmacol. 140(8):1363–1372. Prowell DP, McMichael M, Silvain J-F. 2004. Multilocus genetic analysis of host use, introgression, and speciation in host strains of fall armyworm (Lepidoptera: Noctuidae). Ann Entomol Soc Am. 97(5):1034–1044. R Development Core Team. 2013. R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. Robinson RA, Sutherland WJ. 2002. Post-war changes in arable farming and biodiversity in Great Britain. J Appl Ecol. 39(1):157–176. Roetschi A, Berthoud H, Kuhn R, Imdorf A. 2008. Infection rate based on quantitative real-time PCR of Melissococcus plutonius,thecausal agent of European foulbrood, in honeybee colonies before and after apiary sanitation. Apidologie 39(3):362–371. Rosenkranz P, Aumeier P, Ziegelmann B. 2010. Biology and control of Varroa destructor. J Invertebr Pathol. 103:S96–119. Sabeti PC, et al. 2007. Genome-wide detection and characterization of positive selection in human populations. Nature 449(7164):913–918. S anchez Barreiro F, et al. 2017. Characterizing restriction enzymeassociated loci in historic ragweed (Ambrosia artemisiifolia) voucher specimens using custom-designed RNA probes. Mol Ecol Resour. 17(2):209–220. Parejo et al. GBE 2550 Genome Biol. Evol. 12(12):2535–2551 doi:10.1093/gbe/evaa188 Advance Access publication 2 September 2020 Downloaded from https://academic.oup.com/gbe/article/12/12/2535/5900668 by UNIVERSIDAD DEL PAIS VASCO user on 02 March 2021 Sawyer S, Krause J, Guschanski K, Savolainen V, P€ a€ abo S. 2012. Temporal patterns of nucleotide misincorporations and DNA fragmentation in ancient DNA. PLoS One 7(3):e34131. Seeley TD. 2017. Life-history traits of wild honey bee colonies living in forests around Ithaca, NY, USA. Apidologie 48(6):743–754. Sgolastra F, et al. 2020. Bees and pesticide regulation: lessons from the neonicotinoid experience. Biol Conserv. 241:108356. Shi JF, Zhu TT, Guo WC, Li GQ. 2016. Receptor tyrosine kinase genes respond transcriptionally to sublethal doses of five insecticides by a mode-of-action independent way in Leptinotarsa decemlineata (Say). J Asia-Pac Entomol. 19(4):1103–1110. Simon-Delso N, et al. 2015. Systemic insecticides (neonicotinoids and fipronil): trends, uses, mode of action and metabolites. Environ Sci Pollut Res. 22(1):5–34. Slatkin M. 2008. Linkage disequilibrium—understanding the evolutionary past and mapping the medical future. Nat Rev Genet. 9(6):477–485. Slatkin M, Racimo F. 2016. Ancient DNA and human history. Proc Natl Acad Sci USA. 113(23):6380–6387. Smith NM, et al. 2019. Strikingly high levels of heterozygosity despite 20 years of inbreeding in a clonal honey bee. J Evol Biol. 32(2):144–152. Sproul JS, Maddison DR. 2017. Sequencing historical specimens: successful preparation of small specimens with low amounts of degraded DNA. Mol Ecol Resour. 17(6):1183–1201. Staats M, et al. 2013. Genomic treasure troves: complete genome sequencing of herbarium and insect museum specimens. PLoS One 8(7):e69189. Szpiech ZA, Hernandez RD. 2014. selscan: an efficient multithreaded program to perform EHH-based scans for positive selection. Mol Biol Evol. 31(10):2824–2827. Tarpy DR, vanEngelsdorp D, Pettis JS. 2013. Genetic diversity affects colony survivorship in commercial honey bee colonies. Naturwissenschaften 100(8):723–728. Tewhey R, Bansal V, Torkamani A, Topol EJ, Schork NJ. 2011. The importance of phase information for human genomics. Nat Rev Genet. 12(3):215–223. Thompson HM. 2003. Behavioural effects of pesticides in bees—their potential for use in risk assessment. Ecotoxicology 12(1/4):317–330. Tin MM-Y, Economo EP, Mikheyev AS. 2014. Sequencing degraded DNA from non-destructively sampled museum specimens for RAD-tagging and low-coverage shotgun phylogenetics. PLoS One 9(5):e96793. Toews DPL, Brelsford A. 2012. The biogeography of mitochondrial and nuclear discordance in animals. Mol Ecol. 21(16):3907–3930. Tong F, Coats JR. 2010. Effects of monoterpenoid insecticides on [3H]- TBOB binding in house fly GABA receptor and 36Cluptake in American cockroach ventral nerve cord. Pestic Biochem Physiol. 98(3):317–324. Tsuruda JM, Harris JW, Bourgeois L, Danka RG, Hunt GJ. 2012. High-resolution linkage analyses to identify genes that influence Varroa sensitive hygiene behavior in honey bees. PLoS One 7(11):e48276. Vanengelsdorp D, Meixner MD. 2010. A historical review of managed honey bee populations in Europe and the United States and the factors that may affect them. J Invertebr Pathol. 103:80–95. van Lexmond MB, Bonmatin J-M, Goulson D, Noome DA. 2015. Worldwide integrated assessment on systemic pesticides. Environ Sci Pollut Res. 22(1):1–4. Vitti JJ, Grossman SR, Sabeti PC. 2013. Detecting natural selection in genomic data. Annu Rev Genet. 47(1):97–120. Vogel C, Teichmann SA, Chothia C. 2003. The immunoglobulin superfamily in Drosophila melanogaster and Caenorhabditis elegans and the evolution of complexity. Development 130(25):6317–6328. Waliwitiya R, Belton P, Nicholson RA, Lowenberger CA. 2010. Effects of the essential oil constituent thymol and other neuroactive chemicals on flight motor activity and wing beat frequency in the blowfly Phaenicia sericata. Pest Manag Sci. 66(3):277–289. Wallberg A, et al. 2014. A worldwide survey of genome sequence variation provides insight into the evolutionary history of the honeybee Apis mellifera. Nat Genet. 46(10):1081–1088. Wallberg A, et al. 2019. A hybrid de novo genome assembly of the honeybee, Apis mellifera, with chromosome-length scaffolds. BMC Genomics 20(1):275. Wallberg A, Pirk CW, Allsopp MH, Webster MT. 2016. Identification of multiple loci associated with social parasitism in honeybees. PLoS Genet. 12(6):e1006097. Walton C, et al. 2008. Genetic population structure and introgression in Anopheles dirus mosquitoes in South-east Asia. Mol Ecol. 10(3):569–580. Watson FL. 2005. Extensive diversity of Ig-superfamily proteins in the immune system of insects. Science 309(5742):1874–1878. Weir BS, Cockerham CC. 1984. Estimating F-statistics for the analysis of population structure. Evolution 38(6):1358–1370. Wilson EB. 1905. The chromosomes in relation to the determination of sex in insects. Science 22(564):500–502. Winston ML, Dropkin JA, Taylor OR. 1981. Demography and life history characteristics of two honey bee races (Apis mellifera). Oecologia 48(3):407–413. Wragg D, et al. 2016. Whole-genome resequencing of honeybee drones to detect genomic selection in a population managed for royal jelly. Sci Rep. 6:1–13. Wragg D, et al. 2018. Autosomal and mitochondrial adaptation following admixture: a case study on the honeybees of Reunion Island. Genome Biol Evol. 10(1):220–238. Yon Rhee S, Wood V, Dolinski K, Draghici S. 2008. Use and misuse of the gene ontology annotations. Nat Rev Genet. 9(7):509–515. Zayed A. 2009. Bee genetics and conservation. Apidologie 40(3):237–262. Associate editor: Andrea Betancourt Whole-Genome Sequencing Museum Specimens GBE Genome Biol. Evol. 12(12):2535–2551 doi:10.1093/gbe/evaa188 Advance Access publication 2 September 2020 2551 Downloaded from https://academic.oup.com/gbe/article/12/12/2535/5900668 by UNIVERSIDAD DEL PAIS VASCO user on 02 March 2021