Genome Sequencing of Museum Specimens Reveals Rapid Changes in the Genetic Composition of Honey Bees in California

Genome Sequencing of Museum Specimens Reveals Rapid Changes in the Genetic Composition of Honey... The western honey bee, Apis mellifera, is an enormously influential pollinator in both natural and managed ecosystems. In North America, this species has been introduced numerous times from a variety of different source populations in Europe and Africa. Since then, feral populations have expanded into many different environments across their broad introduced range. Here, we used whole genome sequencing of historical museum specimens and newly collected modern populations from California (USA) to analyze the impact of demography and selection on introduced populations during the past 105 years. We find that populations from both northern and southern California exhibit pronounced genetic changes, but have changed in different ways. In northern populations, honey bees underwent a substantial shift from western European to eastern European ancestry since the 1960s, whereas southern populations are dominated by the introgression of Africanized genomes during the past two decades. Additionally, we identify an isolated island population that has experienced comparatively little change over a large time span. Fine-scale comparison of different populations and time points also revealed SNPs that differ in frequency, highlighting a number of genes that may be important for recent adaptations in these introduced populations. Key words: Apis mellifera, population genomics, demography. Introduction Apis mellifera was first introduced into North America in The western honey bee, Apis mellifera, is the world’s most 1622 (Crane 1999) when populations were brought over by important managed pollinator (Klein et al. 2007) and an iconic French and English colonists. These populations probably home garden visitor. It has been used throughout its native belonged to the M lineage that occurs throughout western range in Europe, Africa, and western Asia as a source of Europe (Ruttner’s designations, Ruttner 1988). In 1853, honey and wax since the Paleolithic era, and the practice of A. mellifera was first introduced to California, when hives keeping bees in hives likely emerged in Egypt between 3000 from the eastern United States arrived in the city of San and 5000 BCE (Crane 1999). Although beekeeping in the Jose (Watkins 1968). Subsequently, in 1859, Italian A. melli- honey bee’s native range has primarily relied on the use of fera bees belonging to the C lineage (Ruttner 1988)were local populations of A. mellifera, beekeepers outside the na- imported to North America and were quickly transported to tive range rely on various introduced, and genetically distinct, California (Crane 1999). Despite occurring in close geographic populations. These populations are now the most important proximity in Europe, the M and C lineages represent two pollinators of agricultural crops (Klein et al. 2007) and under- separate and genetically distinct radiations of A. mellifera standing temporal and spatial changes in genetic variation out of Africa into Europe (Ruttner 1988; Whitfield et al. within and between these populations is vital to the agricul- 2006). However, these two major lineages breed freely and tural industry. there is evidence of admixture in their modern-day contact The Author(s) 2018. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution. This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com 458 Genome Biol. Evol. 10(2):458–472. doi:10.1093/gbe/evy007 Advance Access publication January 15, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/458/4810442 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Genome Sequencing of Museum Specimens GBE zone across central Europe (Cridland et al. 2017). These two the USDA reported a 15% loss of colonies in California be- lineages are the primary contributors to the modern managed tween January and March (USDA National Agriculture honey bee stocks present in California. Statistics Service, 2016). These losses are believed to be influ- In addition to large populations of managed A. mellifera enced by a number of factors that currently threaten man- there are also substantial feral populations of bees in aged honey bee populations including pesticides (Brandt et al. California. Feral bees are ultimately derived from a variety of 2016), diseases and parasites (vanEngelsdorp and Meixner managed stocks and are continuously replenished as man- 2010; Sanchez-Bayo et al. 2016). However, there are indica- aged bees routinely escape and establish new feral colonies. tions that genetically diverse hives may be more resistant to Although managed and feral honey bees share genetic histo- disease (Tarpy and Seeley 2006). ries, they experience very different sets of selective pressures. Several previous studies have attempted to quantify the Managed honey bee populations experience significant selec- contributions of different A. mellifera lineages to the popula- tive pressures for desirable traits (Winston et al. 1983), includ- tions in California and the United States generally (Schiff et al. ing honey production (Guzman-Novoa and Page 1999), 1994; Magnus and Szalanski 2010; Harpur et al. 2012; Kono hygenic behavior, and pathogen resistance (Buchler et al. and Kohn 2015; Rangel et al. 2016), but previous studies used 2010). At the same time, managed honey bees are shielded either mitochondrial markers or small sets of sequenced genes from other selective pressures as beekeepers provide supple- to assess population structure. Here, we present a whole- mental food, shelter, and medication against common patho- genome based study of the contributions of European and gens and parasites. Feral bees, however, experience different African honey bee populations to introduced (feral and man- selective pressures relative to managed honey bees, and thus aged) Californian bees over a 105 year period (1910–2015). may provide insight into the processes of adaptation to local Our analyses of how ancestry and genes change through conditions. These populations may even serve as a genetic time, across the landscape, and among populations provide reservoir for future managed stocks (Sheppard and Huttel insights into adaptive loci and genetic divergence, and are the 1988; Schiff et al. 1994; De la Rua et al. 2009). Feral bees first steps towards understanding ecologically relevant traits also contribute to the pollination of agricultural crops, though and local adaptation (Luikartetal. 2003; Stinchcombe and the amount of that contribution is unclear (Losey and Vaughn Hoekstra 2008; Allendorf et al. 2010). To conduct our analy- 2006). sis, we used paired historical and modern samples from The introduction of tropical African A. mellifera scutellata throughout the state, spanning 65% of the time that honey (lineage A, Ruttner 1988) to Brazil in 1956 (Kerr 1967; Crane bees have existed in California. We hypothesize substantial 1999), and its subsequent expansion into California provided changes in the contributions of various European and African an additional, and dissimilar, source of genetic variation to ancestral lineages to the genetic profile of Californian popu- feral Californian A. mellifera populations. These populations lations over both spatial and temporal dimensions. In partic- have spread rapidly and have quickly altered the genetic land- ular, we expect to see an increase in the contribution of scape of feral population in the southern United States African lineages to southern Californian populations over though they have not entirely replaced European-derived feral time, including the identification of admixed individuals, par- populations (Pinto et al. 2004). Africanized bees were first alleling the arrival and spread of Africanized bees. We further documented in southern California in 1994. Since then these hypothesized a reduction in genetic diversity in northern populations have introgressed northward, and a 2015 study Californian populations in years following Varroa mite out- found evidence of African derived alleles only 40 km south of breaks as feral colonies suffer substantial losses and the Sacramento (Kono and Kohn 2015). expected proportion of feral colonies that are recently derived California is the seasonal home to nearly half of the man- from managed populations increases. Finally, we expect to aged bee colonies in the continental United States (1,140,000 identify candidate genes that may be involved in local adap- in January of 2016 [USDA National Agriculture Statistics tation to the state of California’s diverse ecological regions. Service, 2016]). This makes information about changes in the genetic composition of bee populations in California, Materials and Methods both feral and managed, of particular interest to bee breeders We acquired 29 museum samples of A. mellifera collected in and the agricultural industry. Apis mellifera populations in California between 1910 and 2011 (fig. 1 and table 1). In California have experienced a number of stresses in the last addition, we included a set of A. mellifera individuals from several decades. First, there are concerns about the spread of two publicly available data sets (Harpur et al. 2014; Wallberg Africanized honey bees and the potential of these bees to et al. 2014). We used a SNP data set that we previously gen- interbreed with managed colonies, producing offspring with erated from these two resources (Cridland et al. 2017)to infer unfavorable traits (i.e., aggressive behavior, lower honey pro- patterns of ancestry in introduced A. mellifera populations duction) (Rinderer et al. 1985). Second, the mortality rate of from California. This set included samples from Africa (47 managed colonies has increased in the United States between individuals), western Europe (39 individuals), eastern Europe 1944 and 2008 (vanEngelsdorp and Meixner 2010). In 2016, Genome Biol. Evol. 10(2):458–472 doi:10.1093/gbe/evy007 Advance Access publication January 15, 2018 459 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/458/4810442 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Cridland et al. GBE FIG.1.—Sampling Locations in California. Sky Valley and Idyllwild are collectively referred to as Riverside. (29 individuals), the Middle East (10 individuals), and ten for modern specimens and a reduced elution volume of A. cerana individuals from Japan. 130 ml for museum specimens. A vacuum centrifuge concen- DNA extraction and library preparation for museum and trator was used to increase DNA concentration for samples modern samples were similar. Benchtop, pipettes and extrac- with a Qubit reading less than 2.5 ng/ml. tion forceps were cleaned with bleach and rinsed with distilled Whole-genome library preparation was completed using water prior to use. Filter tips were used of all steps. Museum Nextera DNA Library Preparation Kit (Illumina) following samples were handled individually and separately from mod- standard protocol for low plexity index pooling. Library suc- ern samples. DNA extractions were completed using DNeasy cess was confirmed by Bioanalyzer High Sensitivity DNA kit Blood and Tissue kit (Qiagen) following the standard protocol (Agilent Technologies). DNA sequencing of 100 bp single-end 460 Genome Biol. Evol. 10(2):458–472 doi:10.1093/gbe/evy007 Advance Access publication January 15, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/458/4810442 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Genome Sequencing of Museum Specimens GBE Table 1 Samples Location Longitude Latitude Date N Avalon, Catalina Island, Los Angeles county 118.32 33.34 June 1910, September 2014 2, 5 UC Berkeley Campus, Alameda county 122.26 37.87 July 2011 1 Arcata, Humboldt county 124.07 40.87 January 2015 6 Blue Lake, Humboldt county 123.95 40.88 February 1966 6 Placerita Canyon Nature Area, Los Angeles county 118.46 34.37 August 1999, September 2014 5, 6 Sky Valley, Riverside county 116.27 33.84 November 1983, September 2014 2, 4 Idyllwild, Riverside county 116.72 33.74 May 1999, September 2014 2, 4 San Pedro, Los Angeles 118.29 33.73 January 2002 2 La Grange, Stanislaus county 120.46 37.66 September 1976, September 2014 2, 6 Stebbins Cold Canyon Reserve, Solano county 122.08 38.3 March 1996, September 2014 5, 5 UC Davis Campus, Yolo county 121.75 38.54 1968 and January 2015 2, 6 Noble Apiaries, Dixon CA, Yolo county 121.86 38.43 March 2015 3 C.F.Koehnen Apiary, Knight’s Landing, Yolo county 121.72 38.8 March 2015 2 Wooten’s Golden Queens Apiary, Palo Cedro, Shasta county 122.23 40.61 March 2015 1 reads on Illumina HiSeq2000 was performed at Vincent J. those positions in all samples using the same set of coverage Coates Genomics Sequencing Laboratory at UC Berkeley. requirements as above. A total of 3,890,276 SNPs were iden- We sequenced each of these individuals on between one- tified used for downstream analyses. sixth and one-third of a lane and generated between 0.38 and 32.7 mean coverage of the genomes. Sequences were Ancestry Analyses aligned to the A. mellifera reference genome version 4.5, avail- We ran ADMIXTURE (Alexander et al. 2009) for the complete able from beebase.org, using Bowtie2 with the very-sensitive- set of European, Middle Eastern, African, and Californian local alignment parameters Langmead and Salzberg (2012). individuals for K population values between 2 and 6 to exam- The variation in coverage is largely due to the difficulties in ine population splits at different assumed values. We initially extracting high quality DNA for sequencing from preserved ran this analysis with O group individuals included, but these museum specimens and our oldest sequences, especially those individuals were never separated by ADMIXTURE from C from 1910, have the lowest mean coverage (supplementary group individuals. Running the analysis with the O group re- table 1, Supplementary Material online). Older samples were moved produced the same result with respect to the also more likely to have missing data, though all samples but Californian individuals. We therefore conclude that either one had calls for at least half of the SNPs included. the O group does not contribute to the Californian popula- tions or we are not able to distinguish contributions from the SNP Sets O group from contributions from the C group. We found that We generated two sets of SNPs for analyzing the Californian individuals from the two locations in Riverside County (Sky samples. The first SNPs were a previously generated set of Valley and Idyllwild) were very similar in each analysis and SNPs identified in native range African and European bees we combined these two sites into a single population for (Cridland et al. 2017) and was used for all ancestry analyses. downstream analyses. We used samtools/vcftools to generate SNP calls for each in- We calculated F for pairs of populations collected in ST dividual, requiring a quality score of 30 for reads to be in- California in 2014 and the populations from Africa, western cluded. A minimum coverage of 7 was required to make Europe and eastern Europe (Harpur et al. 2014, Wallberg genotype calls for an individual and we required 2 coverage et al. 2014)using Dadi (Gutenkunst et al. 2009). To examine of an alternate base to make heterozygote calls. The second the patterns of relatedness between individuals and between set was generated from the samples of bees collected in major groups we ran a nonparametric multidimensional scal- California in 2014–2015 and was used to examine differences ing analysis (nMDS) using the R package ecodist version 1.2.2 between Californian populations. Because many of the his- and calculated the stress value and R value for each analysis. torical samples had lower levels of coverage and thus fewer An nMDS analysis represents similarity between individuals by sites where there is genotype information we used the 2014 positioning them in multidimensional space. We ran the ado- samples (41 feralþ 6 domestic) to identify a set of SNP set for nis function from the vegan version 2.3–4 package in R to downstream analyses. To include a SNP in the data set we identify the effect of major groupings on the distance matrix. required that we were able to make a genotype call in 40/47 We performed formal tests of admixture using samples from 2014/2015. We then made genotype calls at ADMIXTOOLS (Patterson et al. 2012). We calculated Genome Biol. Evol. 10(2):458–472 doi:10.1093/gbe/evy007 Advance Access publication January 15, 2018 461 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/458/4810442 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Cridland et al. GBE f statistics for all pairs of source populations for each poten- Benjamini corrected P value (Benjamini and Hochberg 1995) tial target population to identify populations with evidence of to identify significantly enriched terms within clusters. admixture, where admixture is signified by a negative f sta- tistic. For tests where we identified a negative f statistic, we Results calculated the upper and lower bound on the admixture pro- We called SNPs in a set of museum specimens and modern portion. We then kept all results for which we found evidence populations of A. mellifera collected throughout the state of of admixture with an f statistic less than 0.01 and a z-score California between 1910 and 2015 (fig. 1). We acquired mu- at least four standard deviations from the mean. F ratio es- seum specimens for a number of sites in northern and south- timation was performed on the same set of f statistics that ern California as well as from Avalon, Catalina Island off the passed our filters to estimate the proportion of ancestry in the southern coast of California. Sampling bees from the same admixed populations. sites in 2014–2015 generated a modern collection. In addi- tion, we included domestic drones from three different bee Geographic Differentiation breeders in northern California. To identify patterns of varia- tion in Californian populations with respect to changes in We calculated F per site for pairs of populations between ST contributions from ancestral populations we used a previously African, western European, and eastern European groups and curated set of SNPs generated from native range, African and Californian 2014 populations. For these tests we combined European, bees for inferring the demographic history of the southern Californian populations of Placerita and Riverside A. mellifera (Cridland et al. 2017). into a single population and the northern Californian popula- tions of Davis, Humboldt, Stanislaus, and Stebbins into a sin- gle population. We required genotype information for at least Californian Populations Derive Ancestry from Multiple 80% of individuals for each population in a given population Native Range Honey Bee Lineages comparison to include that site in the analysis. We then cal- We tested the hypothesis that Californian A. mellifera derive culated the 95th percentile of F . ST ancestry from some combination of A. mellifera populations We identified all sites for which an alternate SNP was only from eastern Europe (C lineage), western Europe (M lineage), ever seen in one of the ancestral populations. We then com- and Africa (A lineage). We calculated mean p (nucleotide di- pared the SNPs found in only one ancestral population at versity) along 10 kb windows for each or our modern popu- intermediate to high frequency, 10% or greater, to the fre- lations (table 2). The African population had the highest level quencies of those SNPs in the 2014 northern and southern of diversity followed by the southern Californian populations. Californian populations and used a Fisher’s Exact Test to ex- The domesticated drone grouping showed the lowest nucle- amine the proportions of SNPs between each ancestral pop- otide diversity in California. However, it exhibited similar levels ulation and the southern versus northern Californian of variation to that observed in the eastern and western populations. We identified for each site in each population European populations. We find that northern Californian pop- comparison if the site was in an exonic region and, if so, if the ulations have diversity levels intermediate to eastern and west- SNP difference produces a synonymous or nonsynonymous ern European populations. amino acid. For each site identified in an exon we also identified which Native Range Populations Are Differentially Represented amino acid it coded for and the alternate amino acid encoded across California for by the alternate base. For each pair of populations for the modern Californian Modern populations of bees from northern California, south- populations as well as for the Humboldt 1966 versus 2015 ern California, and Avalon were genetically distinct from each populations we calculated a mid-P value based on a 2 by 2 other, and each population was most similar to a different contingency table for the number of reference and alternate native range population. To quantify genetic differentiation allele calls for the two populations (epitools R-package). We across populations from distinct geographic regions across then performed a q-value correction for the set of p-values California, we calculated F between pairs of modern pop- ST using the R package qvalue. ulations using Dadi (Gutenkunst et al. 2009). The southern Californian populations (Riverside County and Placerita Canyon) were more similar to the African population than Gene Ontology Analysis they were to either the central or the western European pop- We used the DAVID version 6.8 Functional Annotation tool ulations (supplementary table 1, Supplementary Material on- (Huang et al. 2009a, 2009b) to identify enriched gene ontol- line). In contrast, the northern Californian populations and the ogy terms. We used the medium stringency for gene classifi- domesticated individuals were most similar to the eastern cation in DAVID and used an enrichment score of 2.5 as the European population. Avalon was most similar to the western minimum score for identifying enriched clusters and the European population. Within California, we found that 462 Genome Biol. Evol. 10(2):458–472 doi:10.1093/gbe/evy007 Advance Access publication January 15, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/458/4810442 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Genome Sequencing of Museum Specimens GBE Table 2 populations, we find evidence of admixture between the Mean p in Californian Populations Western European (M) group the domestic individuals and Population Location Mean p the eastern European (C) group as indicated by negative f scores and large negative Z-scores. In the modern southern Avalon Island 7.14E-05 Davis Northern California 7.62E-05 Californian populations, we observed evidence of admixture Humboldt Northern California 6.43E-05 between the African (A) group and the domesticated popu- Placerita Canyon Southern California 1.05E-04 lations and western European populations (table 3). Riverside Southern California 1.14E-04 For source population pairs where we found evidence of Stanislaus Northern California 7.54E-05 admixture, we calculated the F ratio (Patterson et al. 2012)to Stebbins Northern California 7.46E-05 estimate ancestry proportions for the source populations. We Domestic Northern California 6.38E-05 observed that the proportion of the European M group in A Africa 1.63E-04 modern northern Californian populations ranges from around C Central Europe 5.17E-05 53.5%6 8.3% to 71.4%6 9.2% depending upon the other M Western Europe 8.38E-05 source population considered (table 4). In modern southern Californian populations we found that the contribution of African populations was between 37.6%6 7.3% and modern populations from the Central Valley (Stanislaus, 73.3%6 3.9% depending on the other source population Stebbins, and Davis) were most similar both to each other considered. The variation in percentages is likely a reflection of and to the domesticated populations, also from the Central the populations in question having more than two ancestral Valley (supplementary table 2, Supplementary Material on- populations contributing to the target population. line). The southern Californian populations from Riverside County andPlacerita Canyon were most similar toeachother. Patterns of Differentiation between Modern Californian We tested the hypotheses that Californian bees exhibit and Native Range Populations admixture with between 2 and 5 source populations using To investigate the patterns of differentiation between modern the program ADMIXTURE (Alexander et al. 2009). In the fol- California and native range populations, we conducted an F ST lowing analysis, we also included two outgroup populations: analysis across all sites between the modern Californian pop- A. cerana, a sister species, and Y lineage A. mellifera,which ulation and the African, eastern European, and western are not thought to have contributed to the genetic makeup of European populations. For these population comparisons Californian bees (Magnus and Szalanski 2010). The samples we grouped the modern Humboldt, Stanislaus, Stebbins, from the four native range populations (African [A], western and Davis populations into a single, northern Californian pop- European [M], eastern European [C], and Middle Eastern [Y]) ulation. Similarly, we grouped the modern Placerita Canyon plus A. cerana clustered into five distinct lineages, as expected and the Riverside populations into a single mainland southern and as previously reported (Harpur et al. 2014; Wallberg et al. California population. We found that the southern Californian 2014; Cridland et al. 2017)(fig. 2 and supplementary fig. 1, population was the least differentiated from the African line- Supplementary Material online). Comparing the native range age population, 95th percentile (supplementary table 2, bee population to introduced Californian populations showed Supplementary Material online). For northern California the that the latter derive ancestry from the African, eastern, and F differences were lowest for the eastern European lineage. ST western European populations, which is in line with our Avalon was most differentiated from the eastern European expectations based on historical data documenting the initial lineage. introduction of honey bees to California. We examined the relatedness between the modern indi- Shared Patterns of High Frequency Alternate SNPs in viduals in California and the African, eastern European, and Modern Populations western European population representatives by creating a distance matrix from the genotype data and performing a We examined nonreference SNPs that were at high frequency nonmetric Multidimensional Scaling Analysis (nMDS, (10%) in one or more native range populations and also at stress¼ 0.071, R ¼ 0.989). All Californian individuals were high frequency in one or more of the Californian populations placed intermediate to the native range individuals (fig. 4A) (table 5). Alleles that are observed at high frequency in both a and southern Californian individuals were placed closer to Californian and a native range population can be interpreted African individuals than northern Californian or Avalon as alleles in California that are most likely to have originated individuals. from the corresponding native range population(s). In gen- We performed formal tests of admixture using eral, we observed a higher proportion of SNPs at high fre- ADMIXTOOLS (Patterson et al. 2012) to estimate contribu- quency in both the African population and the southern tions of native range populations to the Californian individuals Californian population than in the African population and (table 3). In the modern Davis, Humboldt, and Stanislaus the northern Californian population. We also observed a Genome Biol. Evol. 10(2):458–472 doi:10.1093/gbe/evy007 Advance Access publication January 15, 2018 463 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/458/4810442 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Cridland et al. GBE Domestic Avalon Northern California Southern California 1.0 0.8 0.6 Ancestral Contributions A lineage M lineage 0.4 C lineage 0.2 0.0 FIG.2.—ADMIXTURE results for K¼ 5for Californianbees. Table 3 Admixture in Modern Californian Populations Source 1 Source 2 Target f_3 std.err Z SNPs a a L U CM Davis 0.052211 0.001976 26.424 31070 0.769 0.841 M Domestic Davis 0.014021 0.003084 4.546 25749 0.037 0.379 C M Humboldt 0.025731 0.002279 11.289 30161 0.804 0.867 A Domestic Placerita 0.038386 0.00093 41.287 42501 0.435 0.493 A M Placerita 0.018291 0.001474 12.405 51054 0.284 0.708 AMRiverside 0.024526 0.001092 22.468 59331 0.481 0.712 A Domestic Riverside 0.024418 0.000992 24.612 51768 0.591 0.603 CM Riverside 0.010761 0.001425 7.554 56143 0.445 0.883 CM Stanislaus 0.065307 0.001788 36.535 31404 0.689 0.789 A Domestic Stanislaus 0.013122 0.001643 7.985 25962 0.092 0.237 Table 4 Estimated Admixture Proportions in Modern Californian Populations Outgroup Ougroup 2 Source 1 Source 2 Target alpha std.err Z Score Cerana Y MC Davis 0.645738 0.088711 7.279 Cerana Y Domestic M Davis 0.59859 0.122403 4.89 Cerana Y MC Humboldt 0.534668 0.083113 6.433 Cerana Y ADomestic Placerita 0.577514 0.039873 14.484 Cerana Y AM Placerita 0.376642 0.073008 5.159 Cerana Y AM Riverside 0.603464 0.060694 9.943 Cerana Y ADomestic Riverside 0.732793 0.03965 18.482 Cerana Y MC Stanislaus 0.714388 0.091494 7.808 Cerana Y Domestic A Stanislaus 0.826711 0.038753 21.333 464 Genome Biol. Evol. 10(2):458–472 doi:10.1093/gbe/evy007 Advance Access publication January 15, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/458/4810442 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Domestic_2014 Modern Domestic Domestic_2014 Modern Domestic Domestic_2014 Modern Domestic Domestic_2014 Modern Domestic Domestic_2014 Modern Domestic Domestic_2014 Modern Domestic Avalon 1910 Avalon_1910 Avalon 1910 Avalon_1910 Modern A Avalon_2014 valon Modern A Avalon_2014 valon Modern A Avalon_2014 valon Modern A Avalon_2014 valon Modern A Avalon_2014 valon Humboldt 1966 Humboldt_1966 Humboldt_1966 Humboldt 1966 Humboldt_1966 Humboldt 1966 Humboldt_1966 Humboldt 1966 Humboldt_1966 Humboldt 1966 Humboldt_1966 Humboldt 1966 Davis 1968 Davis_1968 Davis 1968 Davis_1968 Stanislaus 1976 Stanislaus_1976 Stanislaus 1976 Stanislaus_1976 Stebbins_1996 Stebbins 1996 Stebbins_1996 Stebbins 1996 Stebbins_1996 Stebbins 1996 Stebbins_1996 Stebbins 1996 Stebbins 1996 Stebbins_1996 Berkeley_2011 Berkeley 2011 Davis_2014 Modern Davis Davis_2014 Modern Davis Davis_2014 Modern Davis Davis_2014 Modern Davis Modern Davis Davis_2014 Modern Humboldt Humboldt_2014 Modern Humboldt Humboldt_2014 Modern Humboldt Humboldt_2014 Modern Humboldt Humboldt_2014 Modern Humboldt Humboldt_2014 Modern Humboldt Humboldt_2014 Modern Stanislaus Stanislaus_2014 Modern Stanislaus Stanislaus_2014 Modern Stanislaus Stanislaus_2014 Modern Stanislaus Stanislaus_2014 Stanislaus_2014 Modern Stanislaus Stanislaus_2014 Modern Stanislaus Modern Stebbins Stebbins_2014 Modern Stebbins Stebbins_2014 Modern Stebbins Stebbins_2014 Modern Stebbins Stebbins_2014 Modern Stebbins Stebbins_2014 Riverside 1983 Riverside_1983 Riverside 1983 Riverside_1983 Placerita Canyon 1999 Placerita_1999 Placerita_1999 Placerita Canyon 1999 Placerita_1999 Placerita Canyon 1999 Placerita_1999 Placerita Canyon 1999 Riverside 1999 Riverside_1999 Riverside_1999 Riverside 1999 San Pedro 2002 SanPedro_2002 San Pedro 2002 SanPedro_2002 Modern Placerita Canyon Placerita_2014 Modern Placerita Canyon Placerita_2014 Modern Placerita Canyon Placerita_2014 Modern Placerita Canyon Placerita_2014 Modern Placerita Canyon Placerita_2014 Modern Placerita Canyon Placerita_2014 Riverside_2014 Modern Riverside Riverside_2014 Modern Riverside Riverside_2014 Modern Riverside Riverside_2014 Modern Riverside Riverside_2014 Modern Riverside Riverside_2014 Modern Riverside Modern Riverside Riverside_2014 Modern Riverside Riverside_2014 Genome Sequencing of Museum Specimens GBE Table 5 was largely supplanted by contributions from the eastern High FrequencySNPs European (C) lineage. Moreover, in contrast to southern Ancestral Population Northern Southern Avalon Californian populations (below), modern populations in California (%) California (%) (%) northern California show very little introgression of African All SharedHighFrequency SNPs (A) alleles. Two individuals, one from Stebbins and one from A 2.43 4.82 3.04 Stanislaus County show very small contributions from Africa, C 10.76 7.33 8.73 but this is not seen in any other samples from northern M 10.50 9.55 13.78 California, either historical or modern (fig. 2). Unique Shared High Frequency SNPs We grouped individual bees based on their location into A 0.05 1.29 0.48 the categories of northern California (Davis, Humboldt, C 0.95 0.47 0.33 Stanislaus, and Stebbins), southern California (Placerita, M 0.38 0.39 0.85 Riverside, and San Pedro) and Avalon. We then generated boxplots for each regional population for each time period for which we have data (fig. 3). We find that there is a decline higher proportion of high frequency SNPs in both the in mean western European ancestry in northern California European (C) population and southern California than in over time and at the same time an increase in the mean European (C) population and southern California. eastern European ancestry in these individuals. In southern We then restricted our high frequency native range California, we find an increase in the mean African ancestry SNPs to the set of SNPs we only observe in one of the over time along with a corresponding increase in mean west- native range lineages. We found a greater proportion of ern European ancestry. The Avalon population appears fairly SNPs that are both unique and at high frequency in the consistent over time with modest changes in the mean con- African ancestral population and in southern California tributions of western and eastern European populations. than in Africa and northern California (Fisher’s Exact These changes are also reflected in the nMDS analysis. We Test [FET],¼0) (table 5). Similarly, the southern find that the modern southern Californian populations are Californian population also had a greater proportion of most similar to the African individuals (fig. 4A) than the his- high frequency African SNPs than the Avalon population torical southern Californian populations (fig. 4B). Additionally, (FET, P¼0). We detected a higher proportion of eastern we find that modern Avalon individuals are positioned closer European SNPs in northern California than in southern to the western European individuals than to the eastern California (FET, P¼ 7.62809e-26). However, we did not European individuals (fig. 4A). observe any difference in the proportions of western The more recent arrival and spread of Africanized bees European SNPs between northern and southern (descended from A lineage A. mellifera scutellata) in southern California, but we observed a higher proportion of west- California left a clear signature in our genomic data. The ear- ernEuropean SNPsinAvalonthanwe found in either liest specimens in our data set from this region, the Riverside northern California (FET, 9.803935e-43) or southern individuals from 1983, were collected prior to the arrival of California (FET, 6.480419e-45). Africanized bees and, accordingly, showed no evidence of To determine whether enriched clusters of genes exhibit African ancestry (fig. 2). Individuals collected 5 years after SNPs at high frequency in both native range populations and the arrival of Africanized bees (Riverside 1999), however, be- within a particular Californian population, we performed a gin to show evidence of African introgression with 33% of GO analysis using DAVID version 6.8 (Huang et al. 2009a, individuals deriving some of their ancestry from the African 2009b). We found enriched clusters of high frequency SNPs lineage. Then, southern Californian specimens collected more only in the Africa-southern California comparison (supple- recently (in Riverside and Placerita Canyon in 2014) showed a mentary table 3, Supplementary Material online). These in- substantially greater contribution from Africa in all individuals cluded a number of clusters that are associated with examined with 100% of individuals collected after 2002 de- developmental processes such as growth factor binding pro- riving some of their ancestry from the African lineage (fig. 2 tein and calcium-binding. and table 3). Although also located in southern California, the individu- als from Avalon exhibited lower introgression from African- Temporal Patterns of Native Range Contributions derived alleles, likely due to their isolation on Santa Catalina We observed two main patterns of temporal variation in our Island, 22 miles from the mainland. We observed that the ADMIXTURE results (fig. 2). The oldest honey bee samples modern Avalon population received a relative higher contri- from northern California (Humbolt 1966; Stebbins 1996) ap- bution from western Europe than did Riverside and Placerita pear to have extensive western European (M) ancestry. Over Canyon in 2014. Moreover, only one individual from Avalon, time, however, these populations appear to have undergone out of the five individuals we sampled, shows any evidence of a dramatic genetic shift, as this western European ancestry African ancestry. Genome Biol. Evol. 10(2):458–472 doi:10.1093/gbe/evy007 Advance Access publication January 15, 2018 465 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/458/4810442 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Cridland et al. GBE FIG.3—Boxplots of ADMIXTURE estimates for regional populations over time. Genetic Diversity within California as belonging to intronic or exonic regions within this sub- set. In all of our population comparisons, we found that In addition to the set of SNPs used for our ancestry analysis, there were significantly fewer SNPs in exonic regions com- we identified an additional set of 3,890,279 SNPs in California pared with the number expected if SNPs were distributed that we used to examine both temporal and spatial popula- randomly throughout the genome (supplementary table tion changes within California. The southern populations ap- 4 4, Supplementary Material online). In contrast to exonic pear to be more genetically diverse (mean p is 1.09 e vs. 5 SNPs, we find more SNPs in intronic regions than expected 7.02 e in northern California), likely due to the recent con- in most of our comparisons (supplementary table 4, tribution from African-derived alleles. Supplementary Material online). There are more differen- ces between the Avalon population and the mainland Geographic Patterns of Variation populations than between the all mainland populations, We identified all sites that exhibited differentiation between further supporting the idea that the Avalon population all modern Californian population pairs. For each comparison, has remained genetically isolated from the mainland for we specified a false discovery rate of 1% and identified SNPs a relatively long period of time. 466 Genome Biol. Evol. 10(2):458–472 doi:10.1093/gbe/evy007 Advance Access publication January 15, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/458/4810442 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Genome Sequencing of Museum Specimens GBE FIG.4.—Non-metric Multidimensional Scaling Analyses for ancestral population representatives, modern, and historical Californian individuals. (A) Modern populations. (B) museum collected. We examined the genomic location of differences and populations within California and are good candidates for found that, between northern and southern California, further study. there are peaks of SNPs in genic regions that pass the FDR cutoff. We found a total of 208 genes with nonsy- Temporal Patterns nonymous SNPs within this subset of the data. Vitellogenin, for example, is found in one of these peaks We compared population pairs for each of the seven locations on chromosome 4 with nine synonymous SNPs and one for which we had both historical and a modern samples. The nonsynonymous SNP that pass the FDR cutoff (fig. 5). temporal difference between sampling times ranged from 15 Although most genes have only one nonsynonymous to 104 years. We found that genetic differentiation increased SNP, there were four genes with five or more. In the with increasing time between sampling dates at a given loca- gene zormin, which is involved in actin binding, we found tion (regression R ¼ 0.5649, P¼ 0.01184) (supplementary 16 nonsynonymous SNPs at the 1% FDR level. We found table 5, Supplementary Material online). However, there are ten SNPs in cubilin, which is involved in nephrocyte filtra- more differences between modern northern and modern tion in Drosophila,nine SNPs LOC100578512; hydroceph- southern Californian populations (F ¼ 0.0193) than be- ST alus-inducing protein-like,and five in LOC725833,which tween northern Californian populations from the 1960s and is an uncharacterized gene. 1970s (combined) and the modern populations (mean We also examined the subset of genes with nonsynony- F ¼ 0.0084). The mean F between southern Californian ST ST mous SNPs in the comparisons between modern Avalon and populations from 1999 and 2014 was 0.0148. the mainland populations. Within this differentiated set, We compared all sites where there was a difference be- 2,064 genes exhibited differentiation between Avalon and tween the population from Humboldt County in 1966 and in northern California and 906 genes showed differentiation be- 2014, the only historical to modern comparison where we tween Avalon and southern California. A total of 807 genes have even sample sizes and high quality sequence for all indi- showed differentiation in both comparisons. Odorant recep- viduals. We found that there are fewer SNPs at a FDR of 1% in tor (OR) genes occur in the list of differentiated genes with six exonic regions and more in intronic regions than expected. OR genes differentiated between northern California and We found 306 genes with nonsynonymous SNPs in the Avalon and two OR genes differentiated between southern Humboldt comparison, nine of which have five or more non- California and Avalon. A further four OR genes are differen- synonymous SNPs. Most of these genes were uncharacter- tiated in both comparisons. Vitellogenin is differentiated in ized, but hemolectin, which is involved in clotting, and both population comparisons (fig. 4) and its receptor, yolkless, LOC551016, which is described as midasin-like, are in this is differentiated between Avalon and northern California. set. The overall set of genes with nonsynonymous SNPs is These genes may represent local adaptation between very different between the Humboldt and the modern Genome Biol. Evol. 10(2):458–472 doi:10.1093/gbe/evy007 Advance Access publication January 15, 2018 467 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/458/4810442 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Cridland et al. GBE FIG.5.—Differentiated SNPs in Vitellogenin in modern population comparisons. Vg is a minus strand gene and there are no SNPs observed in the 500 bp upstream of the 5 end. Californian comparison with only 29 genes in both compar- of feral A. mellifera and potential local adaptation of feral isons, and 456 that are in only one comparison. Of the genes populations. that do occur in both comparisons, one is yolkless,a vitello- genin receptor. Modern Populations Populations examined in this study were collected away from Discussion agricultural areas, but cannot be confirmed to be feral bees. Here, we document patterns of genetic diversity of Apis mel- Examining the patterns of genetic differences between pop- lifera populations over temporal and spatial scales in the state ulations can indicate the relatedness of these individuals to of California. We document substantial genetic diversity in the managed stocks. Regardless, these bees serve as indicators of A. mellifera populations that are likely the result of 1) changes the genetic diversity of A. mellifera found in regions through- in beekeeping practices over time, 2) the introduction of out California. We identified three main geographically and Africanized bees to California in 1994, 3) the introduction genetically distinct groups of honey bee populations in of the Varroa mite to California in 1987, and 4) the life history California. First, there are populations of northern 468 Genome Biol. Evol. 10(2):458–472 doi:10.1093/gbe/evy007 Advance Access publication January 15, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/458/4810442 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Genome Sequencing of Museum Specimens GBE Californian bees that closely resemble managed stocks. Prior 1994 provide a source of genetic variation that is currently to the introduction of Varroa mites, feral colonies appear to lacking in northern California, though Africanized bees ap- be generally stable over time, but the demographics of feral pear to be spreading northwards (Kono and Kohn 2015) bees can still change relatively rapidly. The density of feral and our analyses support this conclusion. colonies has been measured in Ithaca, NY was found to We find that African populations contribute between change slowly as established colonies are likely to survive 37.6% and 73.3% of the ancestry of modern southern through the winter (79% survival; Seeley 1978), which is Californian populations. In addition, because of concerns of similar to the survival of feral populations in Arizona (80%; Africanization, modern beekeepers in southern California Taber 1979), where winters are milder. Baum et al. (2005), have altered breeding practices to reduce or eliminate the who examined feral A. mellifera colonies in Texas between introduction of Africanized alleles into domesticated lineages. 1993 and 2000, found that, although colony turnover was 5– This makes it very likely that the southern Californian bees are 30% per year, over that time Africanized colonies increased genuine feral bees. Compared with northern Californian pop- from an initial colony discovered in 1993 to 80.3% of colonies ulations, modern southern Californian populations also ap- Africanized in 2000 (Baum et al. 2005). This baseline rate of pear to have more substantial genetic contributions from turnover gives us a good idea of what we should expect in western European (M) lineages, a pattern which has been terms of the genetic relationships between feral populations observed previously (Whitfield et al. 2006; Rangel et al. and nearby managed bees if we assume that managed pop- 2016). This may be due to a decrease in gene flow between ulations continue to generate new feral colonies and feral managed and feral populations in southern Californian, selec- colonies are not experiencing any additional pressures such tive pressures favoring western European alleles, or hitchhik- as Varroa mites or African introgression. ing of M alleles with African spread. Further studies are The introduction of the Varroa mite has substantially al- required to determine what processes may be at work here. tered the population dynamics of feral colonies. After being This produces a population with high diversity and contribu- first discovered in California in 1987, Varroa reduced the feral tions from both western and eastern European bees as well as colony population of California to less than one fifth its pre- Africanized strains and is similar to the population dynamics of vious size by 1994 (Kraus and Page 1995). Weexpectthatthis feral bees in the Yucatan in between 1993 and 1998 as substantial reduction in the feral bee population would result Africanization of the local populations occurred (Clarke in reduced genetic diversity of modern feral bees both by et al. 2002). Our results support the hypothesis that popula- wiping out most feral colonies that existed prior to Varroa, tions from southern California are experiencing rapid ongoing and also increasing the turnover of feral colonies by reducing changes, possibly due to local adaptation and admixture from the mean feral colony survival time, as was shown in Kraus multiple lineages. and Page (1995). These predictions match well with the pat- The increased genetic diversity we observed in mainland terns we documented in northern California. Historical pop- southern Californian populations bears directly on some of ulations exhibited a shift over time, concurrent with changing the current issues facing A. mellifera populations in beekeeping preferences from M lineage to C lineage popula- California. Genetic diversity has been shown to be an impor- tions (Crane 1999), and resembling more closely domesti- tant factor in colony health, with genetically diverse colonies cated populations, derived mostly from European (C) exhibiting higher disease resistance (Tarpy and Seeley 2006), lineages with some European (M) ancestry as well. Many of increased thermoregulation capacity (Jones et al. 2004), in- these populations were taken from the area at the northern creased diversity of worker responses to environmental stimuli end of California’s Central Valley, an area with high mite in- (Page et al. 1995), and higher productivity and fitness (Mattila festation rates (Kraus and Page 1995). and Seeley 2007). Managed bee colonies have recently been Our analysis based on southern Californian bees revealed found to have higher levels of genetic variation than honey that these populations derive substantial genetic ancestry bee populations in the eastern or western Europe (Harpur from African lineages and this influx of new alleles result in et al. 2012). However, those studies were based on colonies populations that are much more genetically diverse than that exhibited much lower genetic diversity than the northern Californian bees. This diversity is consistent with Africanized populations we documented here. We find that observations in Texas of hybrid swarms that were found in within California, Africanized southern Californian popula- the contact zone of African and European derived bees, post tions exhibit more genetic diversity than managed or modern Africanization (Pinto et al. 2005). The southern Californian northern Californian populations, but the managed and mod- mainland populations appear to be primarily influenced by ern northern Californian populations show similar levels of different factors than northern Californian populations, genetic diversity both with respect to each other and to west- though the timing of the Varroa mite outbreak and the intro- ern and eastern European bees. duction of Africanized bees into the United States may have Africanized bees have been found to be more resistant to played a role in favoring Africanized bee spread (Pinto et al. Varroa mites, likely due to a higher grooming rate and shorter 2004). The introduction of Africanized bees to California in period between egg laying and capping of the larvae Genome Biol. Evol. 10(2):458–472 doi:10.1093/gbe/evy007 Advance Access publication January 15, 2018 469 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/458/4810442 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Cridland et al. GBE (Camazine 1986; Moritz and Mautz 1990; Emsen et al. 2012), western European (M) like bees to primarily eastern European Moreover, it has been hypothesized that Africanization of (C) like bees, but we suspect that it was driven by the reduc- Californian honey bees may assist feral A. mellifera popula- tion of the feral population in the mid 1990s and changes in tions in surviving Varroa by introducing these beneficial breeding practices over that time period. behaviors into southern Californian populations (Kraus and Page 1995). Environmental variables may also play an impor- Differentiated Genes tant role in facilitating mite infestations. Medina-Flores et al. In addition to the broad-scale changes in the genetic compo- (2014) observed an association between drier conditions and sition of modern honey bee populations of California we iden- lower infestation rates in Mexico. Our study is congruent with tified a substantial number of genes that are differentiated this observation, and may help explain the observed patterns between populations. Some genes, like vitellogenin, show of higher genetic diversity in southern California where envi- differentiation between multiple population comparisons. ronmental conditions are similar. Further studies should inves- These genes constitute a candidate sets for both local adap- tigate the relative role of genetic background and tation to the various environmental conditions along environmental conditions with respect to honey bee resis- California, and future studies may investigate their potential tance Varroa and other health threats. role in resistance to Varroa or other diseases faced by honey Our analysis uncovered the existence of a third modern bee populations. We also find some interesting candidates in population in Catalina Island, which is located off the coast the comparison of the 1966 and modern Humboldt popula- of California near Los Angeles. This population is genetically tion, though these populations are composed of only a few distinct from the mainland populations and bears more re- individuals each and the differentiation observed may be due semblance to the western European (M) lineages and the to sampling. The gene spatzle 2,a toll receptor, is highly dif- Humboldt County 1966 population. This population is likely ferentiated between the Humboldt 1966 and Humboldt 2015 to have been on the island for some time, possibly since 1890, populations. It was found to be downregulated in high Varroa when A. mellifera was introduced on nearby Santa Cruz Island treatment groups in a previous study of the effects of de- (Wenner and Thorp 1994) and there are no managed bees on formed wing virus and Varroa on honey bees (Ryabov et al. the island. Additionally, this population is likely to be free of 2014). Other functions such as neuronal development, neu- Varroa mites. The Varroa mite was never found on Santa Cruz ronal sensitivity,and olfaction have been implicated in studies Island until it was introduced in an effort to biologically control of the genetic basis of tolerance to Varroa (Zakar et al. 2014). honey bees (Wenner and Thorp 1994) and is therefore also unlikely to have made it to Santa Catalina Island. Moreover, the genetic differences between the mainland population and Conclusions Santa Catalina Island indicate that there is little to no gene Our study provides an in-depth examination of temporal and flow between the populations, and thus little opportunity for spatial changes in honey bee populations in California. We the transfer of mites. This population provides an interesting identified substantial genetic differences across these two look into past honey bee populations and has likely under- dimensions in populations that are likely the result of both gone adaptation to the local environment. biotic and abiotic factors. Our analysis uncovers the potential genetic signatures that beekeeping practices and the intro- duction of Varroa have had on the changes experienced by Temporal Changes bee populations over time and suggests that southern The temporal changes we observed were dramatically differ- Californian populations may harbor higher genetic variation ent between southern California and northern California. The that couldbe harnessedtoimprove the healthof A. mellifera first Africanized bees to reach California in 1994 have, over and its important role in the agriculture of California. the past couple of decades, contributed substantially to the population structure. In contrast, historical populations of A. Supplementary Material mellifera prior to the 1970s in northern California differ from more modern northern Californian populations primarily due Supplementary data are available at Genome Biology and to a reduction in the contribution from western European (M) Evolution online. lineages and a concomitant increase in the contribution from eastern European (C) lineages. Modern northern Californian Acknowledgments populations in the Central Valley (Davis, Stanislaus, and Stebbins), are also very similar to modern domesticated We would like to thank Dr. Michael Mesler for help in collect- bees, F ¼ 0.059–0.065, which likely reflects frequent dis- ing bees. We would also like to thank C.F. Koehnen and Sons, ST persal of domesticated individuals into feral populations as Phil Hofland, and Glenda Wooten for providing drones. This escaped queens and drones establish new colonies. It is diffi- work was supported by the Berkeley Initiative for Global cult to assess the precise timing of this change from more Change Biology; the Office of the Vice Chancellor for 470 Genome Biol. Evol. 10(2):458–472 doi:10.1093/gbe/evy007 Advance Access publication January 15, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/458/4810442 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Genome Sequencing of Museum Specimens GBE Kerr EW. 1967. The history of the introduction of African bees in Brazil. Research at UC Berkeley; the Gordon and Betty Moore South Afric Bee J. 39:33–35. Foundation [grant number GBMF2983], the Vincent Coates Klein AM, et al. 2007. Importance of pollinators in changing landscapes for Genome Sequencing Facility; National Institutes of Health, S10 world crops. Proc Biol Sci. 274(1608):303–313. Instrumentation Grants [grant numbers S10RR029668 and Kono Y, Kohn JR. 2015. Range and frequency of Africanized honey bees in S10RR027303]; the USDA National Institute of Food and Calfornia (USA). PLoS One 10(9):e0137407. Kraus B, Page RE Jr. 1995. Effect of Varroa jacobsoni (Mesostigmata: Agriculture, Hatch Project [grant number CA-B-INS-0087-H]. Varroidae) on feral Apis mellifera (Hymenopter: Apidae) in California. Environ Entomol. 24(6):1473–1480. Literature Cited Langmead B, Salzberg S. 2012. Fast gapped-read alignment with Bowtie Allendorf FW, Hohenlohe PA, Luikart G. 2010. Genomics and the future of 2. Nat Methods 9(4):357–359. Li H, 1000 Genome Project Data Processing Subgroup, et al. 2009. The conservation genetics. Nat Rev Genet. 11(10):697–709. Alexander DH, Novembre J, Lange K. 2009. Fast model-based estimation sequence alignement/map (SAM) format and SAMtools. Bioinformatics 25(16):2078–2079. of ancestry in unrelated individuals. Genome Res. 19(9):1655–1664. Losey JE, Vaughn M. 2006. The economic value of ecological services Baum KA, Rubink WL, Pinto MA, Coulson RN. 2005. Spatial and temporal distribution and nest site characteristics of feral honey bee provided by insects. BioScience 56:311–323. Luikart G, England PR, Tallmon D, Jordan S, Taberlet P. 2003. The power (Hymenoptera: Apidae) colonies in a coastal prairie landscape. Environ Entomol. 34(3):610–618. and promise of population genomics: from genotyping to genome typing. Nat Rev Genet. 4(12):981–994. Benjamini Y, Hochberg Y. 1995. Controlling the false discovery rate: a Magnus R, Szalanski AL. 2010. Genetic evidence for honey bees (Apis practical and powerful approach to multiple testing. J Roy Stat Soc B 1:289–300. mellifera L.) of middle eastern lineage in the United States. Sociobiology 55:285–296. Brandt A, Gorenflo A, Siede R, Meixner M, Bu ¨ chler R. 2016. The neon- icotinoids thiacloprid, imidacloprid, and clothianidin affect the immu- Mattila HR, Seeley TD. 2007. Genetic diversity in honey bee colonies nocompetence of honey bees (Apis mellifera L.). J Ins Phys 86:40–47. enhances productivity and fitness. Science 317(5836):362–364. Moritz RFA, Mautz D. 1990. Development of Varroa jacobsoni in colonies Buchler R, Berg S, Le Conte Y. 2010. Breeding for resistance to Varroa destructor in Europe. Apidologie 41(3):393–408. of Apis mellifera capensis and Apis mellifera carnica. Apidologie 21(1):53–58. Camazine S. 1986. Differential reproduction of the mite Varroa jacobsoni (Mesostigmata: Varroidae), on Africanized and European bees Medina-Flores CA, Guzman-Novoa E, Hamiduzzaman MM, (Hymenoptera: Apidae). Ann Entomol Soc Am. 79(5):801–803. Are ´ chiga-Flores CF, Lopez-Carlos MA. 2014. Africanized honey bees (Apis mellifera) have low infestation levels of the mite Clarke KE, Rinderer TE, Franck P, Quezada-Euan JG, Oldroyd BP. 2002. The Africanization of honeybees (Apis mellifera L.) of the Yucatan: a study Varroa destructor in different ecological regions in Mexico. Genet Mol Res. 13(3):7282–7293. of a massive hybridization event across time. Evolution 56:1462–1474. Crane E. (1999) The world history of beekeeping and honey hunting. UK: Page RE, Robinson GE, Fondrk MK, Nasr ME. 1995. Effects of worker Tayor and Francis. p. 1–682. genotypic diversity on honey bee colony development and behavior (Apis mellifera L.). Behav Ecol Sociobiol. 36(6):387. Cridland JM, Tsutsui ND, Ramirez SR. 2017. The complex demographic history and evolutionary origin of the western honey bee, Apis melli- Patterson N, et al. 2012. Ancient admixture in human history. Genetics. doi: 10.1534/genetics.112.145037. fera. Genome Biol Evol. doi: 10.1093/gbe/evx009. ~ Pinto AM, Rubink WL, Coulson RN, Patton JC, Johnston JS. 2004. Temporal De la Rua P, Jaffe ´ R, Dall’Olio R, Munoz I, Serrano J. 2009. Biodiversity, conservation and current threats to European honeybees. Apidologie pattern of Africanization in a feral honeybee population from Texas inferred from mitochondrial DNA. Evolution 58(5):1047–1055. 40(3):263–284. Emsen B, Petukhova T, Guzman-Novoa E. 2012. Factors Limiting the Pinto AM, Rubink WL, Patton JC, Coulson RN, Johnston JS. 2005. Growth of Varroa destructor Populations in Selected Honey Bee Africanization in the United States. Genetics 170(4):1653–1665. Rangel J, et al. 2016. Africanization of a feral honey bee (Apis mellifera) (Apis mellifera L.) Colonies. J Anim Vet Adv. 11:4519–4525. Gutenkunst RN, Hernandez RD, Williamson SH, Bustamante CD. 2009. population in south Texas: does a decade make a difference?. Ecol Evol. 6(7):2158–2169. Inferring the joint demographic history of multiple populations from multidimensional SNP data. PLoS Genet. 5(10):e1000695. Rinderer TE, Collins AM, Tucker KW. 1985. Honey production and under- Guzman-Novoa E, Page RE. 1999. Selective breeding of honey bees lying nectar harvesting activities of Africanized and European honey- bees. J Apicult Res. 24(3):161–167. (Hymenoptera: Apidae) in Africanized areas. J Econ Entomol. 92(3):521–525. Ruttner F. (1988) The western honeybee Apis mellifera L. classification and natural history. Biography and taxonomy of honeybees. Berlin Harpur BA, et al. 2014. Population genomics of the honey bee reveals strong signatures of positive selection on worker traits. Proc Natl Acad Heidelberg GmbH: Springer-Verlag. P. 163–175. Sci U S A. 111(7):2614–2619. Ryabov EV, et al. 2014. A virulent strain of deformed wing virus (DWV) of honeybees (Apis mellifera) prevails after Varroa destructor-mediated, or Harpur BA, MinaeI S, Kent CF, Zayed A. 2012. Management increases genetic diversity of honey bees via admixture. Mol Ecol. in vitro, transmission. PLoS Pathogn. doi: 10.1371/journal.ppat.1004230. Sanchez-Bayo F, et al. 2016. Are bee diseases linked to pesticides?—A 21(18):4414–4421. brief review. Environ Inter. 89–90:7–11. Huang DW, Sherman BT, Lempicki RA. 2009a. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Schiff NM, Sheppard WS, Loper GM, Shimanuki H. 1994. Genetic diversity of feral honey bee (Hymenoptera: Apidae) populations in the southern Protoc. 4:44–57. Huang DW, Sherman BT, Lempicki RA. 2009b. Bioinformatics enrichment United States. Ann Entomol Soc Am. 87(6):842–848. Seeley TD. 1978. Life history strategy of the honey bee, Apis mellifera. tools: paths towards the comprehensive functional analysis of large Oecologia 32(1):109–118. gene lists. Nucleic Acids Res. 37:1–13. Jones JC, Myerscough MR, Graham S, Oldroyd BP. 2004. Honey bee nest Sheppard WS, Huettel MD. 1988. Biochemical genetic markers, intraspe- cific variation, and population genetics of the honey bee, Apis melli- thermoregulation: diversity promotes stability. Science 305(5682):402–404. fera. Environ Entomol. 22(1):183–189. Genome Biol. Evol. 10(2):458–472 doi:10.1093/gbe/evy007 Advance Access publication January 15, 2018 471 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/458/4810442 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Cridland et al. GBE Stinchcombe JR, Hoekstra HE. 2008. Combining population genomics and Watkins LH. 1968. California’s first honey bees. Am Bee J. quantitative genetics: finding the genes underlying ecologically impor- 108:190–191. tant traits. Heredity 100(2):158–170. Wenner AM, Thorp RW. 1994. Removal of feral honey bee (Apis mellifera) Taber S III. 1979. A population of feral honey bee colonies. Am Bee J. colonies from Santa Cruz Island. In: Halvorson WL, Meander GL, edi- 119:842–847. tors. The Fourth California Islands Symposium: update on the status of Tarpy DR, et al. 2010. Mating frequencies of Africanized honey bees in the resources. Santa Barbara (CA): Santa Barbara Museum of Natural south western USA. J Apicult Res. 49(4):302–310. History. p. 513–522. Tarpy DR, Seeley TD. 2006. Lower disease infections in honeybee (Apis Winston ML, Taylor OR, Otis GW. 1983. Some differences between tem- mellifera) colonies headed by polyandrous vs monoandrous queens. perate European and tropical African and South American honeybees. Naturwissenschaften 93(4):195. Bee World 64(1):12–21. US Department of Agriculture, National Agricultural Statistics Service. Whitfield CW, et al. 2006. Thrice Out of Africa: Ancient and Recent 2016. Honey Bee Colonies (US Department of Agriculture, Expansions of the Honey Bee, Apis mellifera.Science Washington, DC). Available from: http://usda.mannlib.cornell.edu/ 314(5799):642–645. MannUsda/viewDocumentInfo.do? documentID¼1943; last accessed Zakar E, Javor A, Kusza S. 2014. Genetic basis of tolerance to Varroa February 2017. destructor in honey bees (Apis mellifera L.). Insect Soc. vanEngelsdorp D, Meixner MD. 2010. A historical review of managed 61(3):207–215. honey bee populations in Europe and the United States and the factors Associate editor: Brandon Gaut that may affect them. J Invertebr Pathol. 103:S80–S95. Wallberg, A, et al. 2014. A worldwide survey of genome sequence vari- ation provides insight into the evolutionary history of the honeybee Apis mellifera. Nat. Genet. 46:1081–1088. 472 Genome Biol. Evol. 10(2):458–472 doi:10.1093/gbe/evy007 Advance Access publication January 15, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/458/4810442 by Ed 'DeepDyve' Gillespie user on 16 March 2018 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Genome Biology and Evolution Oxford University Press

Genome Sequencing of Museum Specimens Reveals Rapid Changes in the Genetic Composition of Honey Bees in California

Free
15 pages

Loading next page...
 
/lp/ou_press/genome-sequencing-of-museum-specimens-reveals-rapid-changes-in-the-raerg0MjLK
Publisher
Oxford University Press
Copyright
© The Author(s) 2018. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution.
ISSN
1759-6653
eISSN
1759-6653
D.O.I.
10.1093/gbe/evy007
Publisher site
See Article on Publisher Site

Abstract

The western honey bee, Apis mellifera, is an enormously influential pollinator in both natural and managed ecosystems. In North America, this species has been introduced numerous times from a variety of different source populations in Europe and Africa. Since then, feral populations have expanded into many different environments across their broad introduced range. Here, we used whole genome sequencing of historical museum specimens and newly collected modern populations from California (USA) to analyze the impact of demography and selection on introduced populations during the past 105 years. We find that populations from both northern and southern California exhibit pronounced genetic changes, but have changed in different ways. In northern populations, honey bees underwent a substantial shift from western European to eastern European ancestry since the 1960s, whereas southern populations are dominated by the introgression of Africanized genomes during the past two decades. Additionally, we identify an isolated island population that has experienced comparatively little change over a large time span. Fine-scale comparison of different populations and time points also revealed SNPs that differ in frequency, highlighting a number of genes that may be important for recent adaptations in these introduced populations. Key words: Apis mellifera, population genomics, demography. Introduction Apis mellifera was first introduced into North America in The western honey bee, Apis mellifera, is the world’s most 1622 (Crane 1999) when populations were brought over by important managed pollinator (Klein et al. 2007) and an iconic French and English colonists. These populations probably home garden visitor. It has been used throughout its native belonged to the M lineage that occurs throughout western range in Europe, Africa, and western Asia as a source of Europe (Ruttner’s designations, Ruttner 1988). In 1853, honey and wax since the Paleolithic era, and the practice of A. mellifera was first introduced to California, when hives keeping bees in hives likely emerged in Egypt between 3000 from the eastern United States arrived in the city of San and 5000 BCE (Crane 1999). Although beekeeping in the Jose (Watkins 1968). Subsequently, in 1859, Italian A. melli- honey bee’s native range has primarily relied on the use of fera bees belonging to the C lineage (Ruttner 1988)were local populations of A. mellifera, beekeepers outside the na- imported to North America and were quickly transported to tive range rely on various introduced, and genetically distinct, California (Crane 1999). Despite occurring in close geographic populations. These populations are now the most important proximity in Europe, the M and C lineages represent two pollinators of agricultural crops (Klein et al. 2007) and under- separate and genetically distinct radiations of A. mellifera standing temporal and spatial changes in genetic variation out of Africa into Europe (Ruttner 1988; Whitfield et al. within and between these populations is vital to the agricul- 2006). However, these two major lineages breed freely and tural industry. there is evidence of admixture in their modern-day contact The Author(s) 2018. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution. This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com 458 Genome Biol. Evol. 10(2):458–472. doi:10.1093/gbe/evy007 Advance Access publication January 15, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/458/4810442 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Genome Sequencing of Museum Specimens GBE zone across central Europe (Cridland et al. 2017). These two the USDA reported a 15% loss of colonies in California be- lineages are the primary contributors to the modern managed tween January and March (USDA National Agriculture honey bee stocks present in California. Statistics Service, 2016). These losses are believed to be influ- In addition to large populations of managed A. mellifera enced by a number of factors that currently threaten man- there are also substantial feral populations of bees in aged honey bee populations including pesticides (Brandt et al. California. Feral bees are ultimately derived from a variety of 2016), diseases and parasites (vanEngelsdorp and Meixner managed stocks and are continuously replenished as man- 2010; Sanchez-Bayo et al. 2016). However, there are indica- aged bees routinely escape and establish new feral colonies. tions that genetically diverse hives may be more resistant to Although managed and feral honey bees share genetic histo- disease (Tarpy and Seeley 2006). ries, they experience very different sets of selective pressures. Several previous studies have attempted to quantify the Managed honey bee populations experience significant selec- contributions of different A. mellifera lineages to the popula- tive pressures for desirable traits (Winston et al. 1983), includ- tions in California and the United States generally (Schiff et al. ing honey production (Guzman-Novoa and Page 1999), 1994; Magnus and Szalanski 2010; Harpur et al. 2012; Kono hygenic behavior, and pathogen resistance (Buchler et al. and Kohn 2015; Rangel et al. 2016), but previous studies used 2010). At the same time, managed honey bees are shielded either mitochondrial markers or small sets of sequenced genes from other selective pressures as beekeepers provide supple- to assess population structure. Here, we present a whole- mental food, shelter, and medication against common patho- genome based study of the contributions of European and gens and parasites. Feral bees, however, experience different African honey bee populations to introduced (feral and man- selective pressures relative to managed honey bees, and thus aged) Californian bees over a 105 year period (1910–2015). may provide insight into the processes of adaptation to local Our analyses of how ancestry and genes change through conditions. These populations may even serve as a genetic time, across the landscape, and among populations provide reservoir for future managed stocks (Sheppard and Huttel insights into adaptive loci and genetic divergence, and are the 1988; Schiff et al. 1994; De la Rua et al. 2009). Feral bees first steps towards understanding ecologically relevant traits also contribute to the pollination of agricultural crops, though and local adaptation (Luikartetal. 2003; Stinchcombe and the amount of that contribution is unclear (Losey and Vaughn Hoekstra 2008; Allendorf et al. 2010). To conduct our analy- 2006). sis, we used paired historical and modern samples from The introduction of tropical African A. mellifera scutellata throughout the state, spanning 65% of the time that honey (lineage A, Ruttner 1988) to Brazil in 1956 (Kerr 1967; Crane bees have existed in California. We hypothesize substantial 1999), and its subsequent expansion into California provided changes in the contributions of various European and African an additional, and dissimilar, source of genetic variation to ancestral lineages to the genetic profile of Californian popu- feral Californian A. mellifera populations. These populations lations over both spatial and temporal dimensions. In partic- have spread rapidly and have quickly altered the genetic land- ular, we expect to see an increase in the contribution of scape of feral population in the southern United States African lineages to southern Californian populations over though they have not entirely replaced European-derived feral time, including the identification of admixed individuals, par- populations (Pinto et al. 2004). Africanized bees were first alleling the arrival and spread of Africanized bees. We further documented in southern California in 1994. Since then these hypothesized a reduction in genetic diversity in northern populations have introgressed northward, and a 2015 study Californian populations in years following Varroa mite out- found evidence of African derived alleles only 40 km south of breaks as feral colonies suffer substantial losses and the Sacramento (Kono and Kohn 2015). expected proportion of feral colonies that are recently derived California is the seasonal home to nearly half of the man- from managed populations increases. Finally, we expect to aged bee colonies in the continental United States (1,140,000 identify candidate genes that may be involved in local adap- in January of 2016 [USDA National Agriculture Statistics tation to the state of California’s diverse ecological regions. Service, 2016]). This makes information about changes in the genetic composition of bee populations in California, Materials and Methods both feral and managed, of particular interest to bee breeders We acquired 29 museum samples of A. mellifera collected in and the agricultural industry. Apis mellifera populations in California between 1910 and 2011 (fig. 1 and table 1). In California have experienced a number of stresses in the last addition, we included a set of A. mellifera individuals from several decades. First, there are concerns about the spread of two publicly available data sets (Harpur et al. 2014; Wallberg Africanized honey bees and the potential of these bees to et al. 2014). We used a SNP data set that we previously gen- interbreed with managed colonies, producing offspring with erated from these two resources (Cridland et al. 2017)to infer unfavorable traits (i.e., aggressive behavior, lower honey pro- patterns of ancestry in introduced A. mellifera populations duction) (Rinderer et al. 1985). Second, the mortality rate of from California. This set included samples from Africa (47 managed colonies has increased in the United States between individuals), western Europe (39 individuals), eastern Europe 1944 and 2008 (vanEngelsdorp and Meixner 2010). In 2016, Genome Biol. Evol. 10(2):458–472 doi:10.1093/gbe/evy007 Advance Access publication January 15, 2018 459 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/458/4810442 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Cridland et al. GBE FIG.1.—Sampling Locations in California. Sky Valley and Idyllwild are collectively referred to as Riverside. (29 individuals), the Middle East (10 individuals), and ten for modern specimens and a reduced elution volume of A. cerana individuals from Japan. 130 ml for museum specimens. A vacuum centrifuge concen- DNA extraction and library preparation for museum and trator was used to increase DNA concentration for samples modern samples were similar. Benchtop, pipettes and extrac- with a Qubit reading less than 2.5 ng/ml. tion forceps were cleaned with bleach and rinsed with distilled Whole-genome library preparation was completed using water prior to use. Filter tips were used of all steps. Museum Nextera DNA Library Preparation Kit (Illumina) following samples were handled individually and separately from mod- standard protocol for low plexity index pooling. Library suc- ern samples. DNA extractions were completed using DNeasy cess was confirmed by Bioanalyzer High Sensitivity DNA kit Blood and Tissue kit (Qiagen) following the standard protocol (Agilent Technologies). DNA sequencing of 100 bp single-end 460 Genome Biol. Evol. 10(2):458–472 doi:10.1093/gbe/evy007 Advance Access publication January 15, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/458/4810442 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Genome Sequencing of Museum Specimens GBE Table 1 Samples Location Longitude Latitude Date N Avalon, Catalina Island, Los Angeles county 118.32 33.34 June 1910, September 2014 2, 5 UC Berkeley Campus, Alameda county 122.26 37.87 July 2011 1 Arcata, Humboldt county 124.07 40.87 January 2015 6 Blue Lake, Humboldt county 123.95 40.88 February 1966 6 Placerita Canyon Nature Area, Los Angeles county 118.46 34.37 August 1999, September 2014 5, 6 Sky Valley, Riverside county 116.27 33.84 November 1983, September 2014 2, 4 Idyllwild, Riverside county 116.72 33.74 May 1999, September 2014 2, 4 San Pedro, Los Angeles 118.29 33.73 January 2002 2 La Grange, Stanislaus county 120.46 37.66 September 1976, September 2014 2, 6 Stebbins Cold Canyon Reserve, Solano county 122.08 38.3 March 1996, September 2014 5, 5 UC Davis Campus, Yolo county 121.75 38.54 1968 and January 2015 2, 6 Noble Apiaries, Dixon CA, Yolo county 121.86 38.43 March 2015 3 C.F.Koehnen Apiary, Knight’s Landing, Yolo county 121.72 38.8 March 2015 2 Wooten’s Golden Queens Apiary, Palo Cedro, Shasta county 122.23 40.61 March 2015 1 reads on Illumina HiSeq2000 was performed at Vincent J. those positions in all samples using the same set of coverage Coates Genomics Sequencing Laboratory at UC Berkeley. requirements as above. A total of 3,890,276 SNPs were iden- We sequenced each of these individuals on between one- tified used for downstream analyses. sixth and one-third of a lane and generated between 0.38 and 32.7 mean coverage of the genomes. Sequences were Ancestry Analyses aligned to the A. mellifera reference genome version 4.5, avail- We ran ADMIXTURE (Alexander et al. 2009) for the complete able from beebase.org, using Bowtie2 with the very-sensitive- set of European, Middle Eastern, African, and Californian local alignment parameters Langmead and Salzberg (2012). individuals for K population values between 2 and 6 to exam- The variation in coverage is largely due to the difficulties in ine population splits at different assumed values. We initially extracting high quality DNA for sequencing from preserved ran this analysis with O group individuals included, but these museum specimens and our oldest sequences, especially those individuals were never separated by ADMIXTURE from C from 1910, have the lowest mean coverage (supplementary group individuals. Running the analysis with the O group re- table 1, Supplementary Material online). Older samples were moved produced the same result with respect to the also more likely to have missing data, though all samples but Californian individuals. We therefore conclude that either one had calls for at least half of the SNPs included. the O group does not contribute to the Californian popula- tions or we are not able to distinguish contributions from the SNP Sets O group from contributions from the C group. We found that We generated two sets of SNPs for analyzing the Californian individuals from the two locations in Riverside County (Sky samples. The first SNPs were a previously generated set of Valley and Idyllwild) were very similar in each analysis and SNPs identified in native range African and European bees we combined these two sites into a single population for (Cridland et al. 2017) and was used for all ancestry analyses. downstream analyses. We used samtools/vcftools to generate SNP calls for each in- We calculated F for pairs of populations collected in ST dividual, requiring a quality score of 30 for reads to be in- California in 2014 and the populations from Africa, western cluded. A minimum coverage of 7 was required to make Europe and eastern Europe (Harpur et al. 2014, Wallberg genotype calls for an individual and we required 2 coverage et al. 2014)using Dadi (Gutenkunst et al. 2009). To examine of an alternate base to make heterozygote calls. The second the patterns of relatedness between individuals and between set was generated from the samples of bees collected in major groups we ran a nonparametric multidimensional scal- California in 2014–2015 and was used to examine differences ing analysis (nMDS) using the R package ecodist version 1.2.2 between Californian populations. Because many of the his- and calculated the stress value and R value for each analysis. torical samples had lower levels of coverage and thus fewer An nMDS analysis represents similarity between individuals by sites where there is genotype information we used the 2014 positioning them in multidimensional space. We ran the ado- samples (41 feralþ 6 domestic) to identify a set of SNP set for nis function from the vegan version 2.3–4 package in R to downstream analyses. To include a SNP in the data set we identify the effect of major groupings on the distance matrix. required that we were able to make a genotype call in 40/47 We performed formal tests of admixture using samples from 2014/2015. We then made genotype calls at ADMIXTOOLS (Patterson et al. 2012). We calculated Genome Biol. Evol. 10(2):458–472 doi:10.1093/gbe/evy007 Advance Access publication January 15, 2018 461 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/458/4810442 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Cridland et al. GBE f statistics for all pairs of source populations for each poten- Benjamini corrected P value (Benjamini and Hochberg 1995) tial target population to identify populations with evidence of to identify significantly enriched terms within clusters. admixture, where admixture is signified by a negative f sta- tistic. For tests where we identified a negative f statistic, we Results calculated the upper and lower bound on the admixture pro- We called SNPs in a set of museum specimens and modern portion. We then kept all results for which we found evidence populations of A. mellifera collected throughout the state of of admixture with an f statistic less than 0.01 and a z-score California between 1910 and 2015 (fig. 1). We acquired mu- at least four standard deviations from the mean. F ratio es- seum specimens for a number of sites in northern and south- timation was performed on the same set of f statistics that ern California as well as from Avalon, Catalina Island off the passed our filters to estimate the proportion of ancestry in the southern coast of California. Sampling bees from the same admixed populations. sites in 2014–2015 generated a modern collection. In addi- tion, we included domestic drones from three different bee Geographic Differentiation breeders in northern California. To identify patterns of varia- tion in Californian populations with respect to changes in We calculated F per site for pairs of populations between ST contributions from ancestral populations we used a previously African, western European, and eastern European groups and curated set of SNPs generated from native range, African and Californian 2014 populations. For these tests we combined European, bees for inferring the demographic history of the southern Californian populations of Placerita and Riverside A. mellifera (Cridland et al. 2017). into a single population and the northern Californian popula- tions of Davis, Humboldt, Stanislaus, and Stebbins into a sin- gle population. We required genotype information for at least Californian Populations Derive Ancestry from Multiple 80% of individuals for each population in a given population Native Range Honey Bee Lineages comparison to include that site in the analysis. We then cal- We tested the hypothesis that Californian A. mellifera derive culated the 95th percentile of F . ST ancestry from some combination of A. mellifera populations We identified all sites for which an alternate SNP was only from eastern Europe (C lineage), western Europe (M lineage), ever seen in one of the ancestral populations. We then com- and Africa (A lineage). We calculated mean p (nucleotide di- pared the SNPs found in only one ancestral population at versity) along 10 kb windows for each or our modern popu- intermediate to high frequency, 10% or greater, to the fre- lations (table 2). The African population had the highest level quencies of those SNPs in the 2014 northern and southern of diversity followed by the southern Californian populations. Californian populations and used a Fisher’s Exact Test to ex- The domesticated drone grouping showed the lowest nucle- amine the proportions of SNPs between each ancestral pop- otide diversity in California. However, it exhibited similar levels ulation and the southern versus northern Californian of variation to that observed in the eastern and western populations. We identified for each site in each population European populations. We find that northern Californian pop- comparison if the site was in an exonic region and, if so, if the ulations have diversity levels intermediate to eastern and west- SNP difference produces a synonymous or nonsynonymous ern European populations. amino acid. For each site identified in an exon we also identified which Native Range Populations Are Differentially Represented amino acid it coded for and the alternate amino acid encoded across California for by the alternate base. For each pair of populations for the modern Californian Modern populations of bees from northern California, south- populations as well as for the Humboldt 1966 versus 2015 ern California, and Avalon were genetically distinct from each populations we calculated a mid-P value based on a 2 by 2 other, and each population was most similar to a different contingency table for the number of reference and alternate native range population. To quantify genetic differentiation allele calls for the two populations (epitools R-package). We across populations from distinct geographic regions across then performed a q-value correction for the set of p-values California, we calculated F between pairs of modern pop- ST using the R package qvalue. ulations using Dadi (Gutenkunst et al. 2009). The southern Californian populations (Riverside County and Placerita Canyon) were more similar to the African population than Gene Ontology Analysis they were to either the central or the western European pop- We used the DAVID version 6.8 Functional Annotation tool ulations (supplementary table 1, Supplementary Material on- (Huang et al. 2009a, 2009b) to identify enriched gene ontol- line). In contrast, the northern Californian populations and the ogy terms. We used the medium stringency for gene classifi- domesticated individuals were most similar to the eastern cation in DAVID and used an enrichment score of 2.5 as the European population. Avalon was most similar to the western minimum score for identifying enriched clusters and the European population. Within California, we found that 462 Genome Biol. Evol. 10(2):458–472 doi:10.1093/gbe/evy007 Advance Access publication January 15, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/458/4810442 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Genome Sequencing of Museum Specimens GBE Table 2 populations, we find evidence of admixture between the Mean p in Californian Populations Western European (M) group the domestic individuals and Population Location Mean p the eastern European (C) group as indicated by negative f scores and large negative Z-scores. In the modern southern Avalon Island 7.14E-05 Davis Northern California 7.62E-05 Californian populations, we observed evidence of admixture Humboldt Northern California 6.43E-05 between the African (A) group and the domesticated popu- Placerita Canyon Southern California 1.05E-04 lations and western European populations (table 3). Riverside Southern California 1.14E-04 For source population pairs where we found evidence of Stanislaus Northern California 7.54E-05 admixture, we calculated the F ratio (Patterson et al. 2012)to Stebbins Northern California 7.46E-05 estimate ancestry proportions for the source populations. We Domestic Northern California 6.38E-05 observed that the proportion of the European M group in A Africa 1.63E-04 modern northern Californian populations ranges from around C Central Europe 5.17E-05 53.5%6 8.3% to 71.4%6 9.2% depending upon the other M Western Europe 8.38E-05 source population considered (table 4). In modern southern Californian populations we found that the contribution of African populations was between 37.6%6 7.3% and modern populations from the Central Valley (Stanislaus, 73.3%6 3.9% depending on the other source population Stebbins, and Davis) were most similar both to each other considered. The variation in percentages is likely a reflection of and to the domesticated populations, also from the Central the populations in question having more than two ancestral Valley (supplementary table 2, Supplementary Material on- populations contributing to the target population. line). The southern Californian populations from Riverside County andPlacerita Canyon were most similar toeachother. Patterns of Differentiation between Modern Californian We tested the hypotheses that Californian bees exhibit and Native Range Populations admixture with between 2 and 5 source populations using To investigate the patterns of differentiation between modern the program ADMIXTURE (Alexander et al. 2009). In the fol- California and native range populations, we conducted an F ST lowing analysis, we also included two outgroup populations: analysis across all sites between the modern Californian pop- A. cerana, a sister species, and Y lineage A. mellifera,which ulation and the African, eastern European, and western are not thought to have contributed to the genetic makeup of European populations. For these population comparisons Californian bees (Magnus and Szalanski 2010). The samples we grouped the modern Humboldt, Stanislaus, Stebbins, from the four native range populations (African [A], western and Davis populations into a single, northern Californian pop- European [M], eastern European [C], and Middle Eastern [Y]) ulation. Similarly, we grouped the modern Placerita Canyon plus A. cerana clustered into five distinct lineages, as expected and the Riverside populations into a single mainland southern and as previously reported (Harpur et al. 2014; Wallberg et al. California population. We found that the southern Californian 2014; Cridland et al. 2017)(fig. 2 and supplementary fig. 1, population was the least differentiated from the African line- Supplementary Material online). Comparing the native range age population, 95th percentile (supplementary table 2, bee population to introduced Californian populations showed Supplementary Material online). For northern California the that the latter derive ancestry from the African, eastern, and F differences were lowest for the eastern European lineage. ST western European populations, which is in line with our Avalon was most differentiated from the eastern European expectations based on historical data documenting the initial lineage. introduction of honey bees to California. We examined the relatedness between the modern indi- Shared Patterns of High Frequency Alternate SNPs in viduals in California and the African, eastern European, and Modern Populations western European population representatives by creating a distance matrix from the genotype data and performing a We examined nonreference SNPs that were at high frequency nonmetric Multidimensional Scaling Analysis (nMDS, (10%) in one or more native range populations and also at stress¼ 0.071, R ¼ 0.989). All Californian individuals were high frequency in one or more of the Californian populations placed intermediate to the native range individuals (fig. 4A) (table 5). Alleles that are observed at high frequency in both a and southern Californian individuals were placed closer to Californian and a native range population can be interpreted African individuals than northern Californian or Avalon as alleles in California that are most likely to have originated individuals. from the corresponding native range population(s). In gen- We performed formal tests of admixture using eral, we observed a higher proportion of SNPs at high fre- ADMIXTOOLS (Patterson et al. 2012) to estimate contribu- quency in both the African population and the southern tions of native range populations to the Californian individuals Californian population than in the African population and (table 3). In the modern Davis, Humboldt, and Stanislaus the northern Californian population. We also observed a Genome Biol. Evol. 10(2):458–472 doi:10.1093/gbe/evy007 Advance Access publication January 15, 2018 463 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/458/4810442 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Cridland et al. GBE Domestic Avalon Northern California Southern California 1.0 0.8 0.6 Ancestral Contributions A lineage M lineage 0.4 C lineage 0.2 0.0 FIG.2.—ADMIXTURE results for K¼ 5for Californianbees. Table 3 Admixture in Modern Californian Populations Source 1 Source 2 Target f_3 std.err Z SNPs a a L U CM Davis 0.052211 0.001976 26.424 31070 0.769 0.841 M Domestic Davis 0.014021 0.003084 4.546 25749 0.037 0.379 C M Humboldt 0.025731 0.002279 11.289 30161 0.804 0.867 A Domestic Placerita 0.038386 0.00093 41.287 42501 0.435 0.493 A M Placerita 0.018291 0.001474 12.405 51054 0.284 0.708 AMRiverside 0.024526 0.001092 22.468 59331 0.481 0.712 A Domestic Riverside 0.024418 0.000992 24.612 51768 0.591 0.603 CM Riverside 0.010761 0.001425 7.554 56143 0.445 0.883 CM Stanislaus 0.065307 0.001788 36.535 31404 0.689 0.789 A Domestic Stanislaus 0.013122 0.001643 7.985 25962 0.092 0.237 Table 4 Estimated Admixture Proportions in Modern Californian Populations Outgroup Ougroup 2 Source 1 Source 2 Target alpha std.err Z Score Cerana Y MC Davis 0.645738 0.088711 7.279 Cerana Y Domestic M Davis 0.59859 0.122403 4.89 Cerana Y MC Humboldt 0.534668 0.083113 6.433 Cerana Y ADomestic Placerita 0.577514 0.039873 14.484 Cerana Y AM Placerita 0.376642 0.073008 5.159 Cerana Y AM Riverside 0.603464 0.060694 9.943 Cerana Y ADomestic Riverside 0.732793 0.03965 18.482 Cerana Y MC Stanislaus 0.714388 0.091494 7.808 Cerana Y Domestic A Stanislaus 0.826711 0.038753 21.333 464 Genome Biol. Evol. 10(2):458–472 doi:10.1093/gbe/evy007 Advance Access publication January 15, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/458/4810442 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Domestic_2014 Modern Domestic Domestic_2014 Modern Domestic Domestic_2014 Modern Domestic Domestic_2014 Modern Domestic Domestic_2014 Modern Domestic Domestic_2014 Modern Domestic Avalon 1910 Avalon_1910 Avalon 1910 Avalon_1910 Modern A Avalon_2014 valon Modern A Avalon_2014 valon Modern A Avalon_2014 valon Modern A Avalon_2014 valon Modern A Avalon_2014 valon Humboldt 1966 Humboldt_1966 Humboldt_1966 Humboldt 1966 Humboldt_1966 Humboldt 1966 Humboldt_1966 Humboldt 1966 Humboldt_1966 Humboldt 1966 Humboldt_1966 Humboldt 1966 Davis 1968 Davis_1968 Davis 1968 Davis_1968 Stanislaus 1976 Stanislaus_1976 Stanislaus 1976 Stanislaus_1976 Stebbins_1996 Stebbins 1996 Stebbins_1996 Stebbins 1996 Stebbins_1996 Stebbins 1996 Stebbins_1996 Stebbins 1996 Stebbins 1996 Stebbins_1996 Berkeley_2011 Berkeley 2011 Davis_2014 Modern Davis Davis_2014 Modern Davis Davis_2014 Modern Davis Davis_2014 Modern Davis Modern Davis Davis_2014 Modern Humboldt Humboldt_2014 Modern Humboldt Humboldt_2014 Modern Humboldt Humboldt_2014 Modern Humboldt Humboldt_2014 Modern Humboldt Humboldt_2014 Modern Humboldt Humboldt_2014 Modern Stanislaus Stanislaus_2014 Modern Stanislaus Stanislaus_2014 Modern Stanislaus Stanislaus_2014 Modern Stanislaus Stanislaus_2014 Stanislaus_2014 Modern Stanislaus Stanislaus_2014 Modern Stanislaus Modern Stebbins Stebbins_2014 Modern Stebbins Stebbins_2014 Modern Stebbins Stebbins_2014 Modern Stebbins Stebbins_2014 Modern Stebbins Stebbins_2014 Riverside 1983 Riverside_1983 Riverside 1983 Riverside_1983 Placerita Canyon 1999 Placerita_1999 Placerita_1999 Placerita Canyon 1999 Placerita_1999 Placerita Canyon 1999 Placerita_1999 Placerita Canyon 1999 Riverside 1999 Riverside_1999 Riverside_1999 Riverside 1999 San Pedro 2002 SanPedro_2002 San Pedro 2002 SanPedro_2002 Modern Placerita Canyon Placerita_2014 Modern Placerita Canyon Placerita_2014 Modern Placerita Canyon Placerita_2014 Modern Placerita Canyon Placerita_2014 Modern Placerita Canyon Placerita_2014 Modern Placerita Canyon Placerita_2014 Riverside_2014 Modern Riverside Riverside_2014 Modern Riverside Riverside_2014 Modern Riverside Riverside_2014 Modern Riverside Riverside_2014 Modern Riverside Riverside_2014 Modern Riverside Modern Riverside Riverside_2014 Modern Riverside Riverside_2014 Genome Sequencing of Museum Specimens GBE Table 5 was largely supplanted by contributions from the eastern High FrequencySNPs European (C) lineage. Moreover, in contrast to southern Ancestral Population Northern Southern Avalon Californian populations (below), modern populations in California (%) California (%) (%) northern California show very little introgression of African All SharedHighFrequency SNPs (A) alleles. Two individuals, one from Stebbins and one from A 2.43 4.82 3.04 Stanislaus County show very small contributions from Africa, C 10.76 7.33 8.73 but this is not seen in any other samples from northern M 10.50 9.55 13.78 California, either historical or modern (fig. 2). Unique Shared High Frequency SNPs We grouped individual bees based on their location into A 0.05 1.29 0.48 the categories of northern California (Davis, Humboldt, C 0.95 0.47 0.33 Stanislaus, and Stebbins), southern California (Placerita, M 0.38 0.39 0.85 Riverside, and San Pedro) and Avalon. We then generated boxplots for each regional population for each time period for which we have data (fig. 3). We find that there is a decline higher proportion of high frequency SNPs in both the in mean western European ancestry in northern California European (C) population and southern California than in over time and at the same time an increase in the mean European (C) population and southern California. eastern European ancestry in these individuals. In southern We then restricted our high frequency native range California, we find an increase in the mean African ancestry SNPs to the set of SNPs we only observe in one of the over time along with a corresponding increase in mean west- native range lineages. We found a greater proportion of ern European ancestry. The Avalon population appears fairly SNPs that are both unique and at high frequency in the consistent over time with modest changes in the mean con- African ancestral population and in southern California tributions of western and eastern European populations. than in Africa and northern California (Fisher’s Exact These changes are also reflected in the nMDS analysis. We Test [FET],¼0) (table 5). Similarly, the southern find that the modern southern Californian populations are Californian population also had a greater proportion of most similar to the African individuals (fig. 4A) than the his- high frequency African SNPs than the Avalon population torical southern Californian populations (fig. 4B). Additionally, (FET, P¼0). We detected a higher proportion of eastern we find that modern Avalon individuals are positioned closer European SNPs in northern California than in southern to the western European individuals than to the eastern California (FET, P¼ 7.62809e-26). However, we did not European individuals (fig. 4A). observe any difference in the proportions of western The more recent arrival and spread of Africanized bees European SNPs between northern and southern (descended from A lineage A. mellifera scutellata) in southern California, but we observed a higher proportion of west- California left a clear signature in our genomic data. The ear- ernEuropean SNPsinAvalonthanwe found in either liest specimens in our data set from this region, the Riverside northern California (FET, 9.803935e-43) or southern individuals from 1983, were collected prior to the arrival of California (FET, 6.480419e-45). Africanized bees and, accordingly, showed no evidence of To determine whether enriched clusters of genes exhibit African ancestry (fig. 2). Individuals collected 5 years after SNPs at high frequency in both native range populations and the arrival of Africanized bees (Riverside 1999), however, be- within a particular Californian population, we performed a gin to show evidence of African introgression with 33% of GO analysis using DAVID version 6.8 (Huang et al. 2009a, individuals deriving some of their ancestry from the African 2009b). We found enriched clusters of high frequency SNPs lineage. Then, southern Californian specimens collected more only in the Africa-southern California comparison (supple- recently (in Riverside and Placerita Canyon in 2014) showed a mentary table 3, Supplementary Material online). These in- substantially greater contribution from Africa in all individuals cluded a number of clusters that are associated with examined with 100% of individuals collected after 2002 de- developmental processes such as growth factor binding pro- riving some of their ancestry from the African lineage (fig. 2 tein and calcium-binding. and table 3). Although also located in southern California, the individu- als from Avalon exhibited lower introgression from African- Temporal Patterns of Native Range Contributions derived alleles, likely due to their isolation on Santa Catalina We observed two main patterns of temporal variation in our Island, 22 miles from the mainland. We observed that the ADMIXTURE results (fig. 2). The oldest honey bee samples modern Avalon population received a relative higher contri- from northern California (Humbolt 1966; Stebbins 1996) ap- bution from western Europe than did Riverside and Placerita pear to have extensive western European (M) ancestry. Over Canyon in 2014. Moreover, only one individual from Avalon, time, however, these populations appear to have undergone out of the five individuals we sampled, shows any evidence of a dramatic genetic shift, as this western European ancestry African ancestry. Genome Biol. Evol. 10(2):458–472 doi:10.1093/gbe/evy007 Advance Access publication January 15, 2018 465 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/458/4810442 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Cridland et al. GBE FIG.3—Boxplots of ADMIXTURE estimates for regional populations over time. Genetic Diversity within California as belonging to intronic or exonic regions within this sub- set. In all of our population comparisons, we found that In addition to the set of SNPs used for our ancestry analysis, there were significantly fewer SNPs in exonic regions com- we identified an additional set of 3,890,279 SNPs in California pared with the number expected if SNPs were distributed that we used to examine both temporal and spatial popula- randomly throughout the genome (supplementary table tion changes within California. The southern populations ap- 4 4, Supplementary Material online). In contrast to exonic pear to be more genetically diverse (mean p is 1.09 e vs. 5 SNPs, we find more SNPs in intronic regions than expected 7.02 e in northern California), likely due to the recent con- in most of our comparisons (supplementary table 4, tribution from African-derived alleles. Supplementary Material online). There are more differen- ces between the Avalon population and the mainland Geographic Patterns of Variation populations than between the all mainland populations, We identified all sites that exhibited differentiation between further supporting the idea that the Avalon population all modern Californian population pairs. For each comparison, has remained genetically isolated from the mainland for we specified a false discovery rate of 1% and identified SNPs a relatively long period of time. 466 Genome Biol. Evol. 10(2):458–472 doi:10.1093/gbe/evy007 Advance Access publication January 15, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/458/4810442 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Genome Sequencing of Museum Specimens GBE FIG.4.—Non-metric Multidimensional Scaling Analyses for ancestral population representatives, modern, and historical Californian individuals. (A) Modern populations. (B) museum collected. We examined the genomic location of differences and populations within California and are good candidates for found that, between northern and southern California, further study. there are peaks of SNPs in genic regions that pass the FDR cutoff. We found a total of 208 genes with nonsy- Temporal Patterns nonymous SNPs within this subset of the data. Vitellogenin, for example, is found in one of these peaks We compared population pairs for each of the seven locations on chromosome 4 with nine synonymous SNPs and one for which we had both historical and a modern samples. The nonsynonymous SNP that pass the FDR cutoff (fig. 5). temporal difference between sampling times ranged from 15 Although most genes have only one nonsynonymous to 104 years. We found that genetic differentiation increased SNP, there were four genes with five or more. In the with increasing time between sampling dates at a given loca- gene zormin, which is involved in actin binding, we found tion (regression R ¼ 0.5649, P¼ 0.01184) (supplementary 16 nonsynonymous SNPs at the 1% FDR level. We found table 5, Supplementary Material online). However, there are ten SNPs in cubilin, which is involved in nephrocyte filtra- more differences between modern northern and modern tion in Drosophila,nine SNPs LOC100578512; hydroceph- southern Californian populations (F ¼ 0.0193) than be- ST alus-inducing protein-like,and five in LOC725833,which tween northern Californian populations from the 1960s and is an uncharacterized gene. 1970s (combined) and the modern populations (mean We also examined the subset of genes with nonsynony- F ¼ 0.0084). The mean F between southern Californian ST ST mous SNPs in the comparisons between modern Avalon and populations from 1999 and 2014 was 0.0148. the mainland populations. Within this differentiated set, We compared all sites where there was a difference be- 2,064 genes exhibited differentiation between Avalon and tween the population from Humboldt County in 1966 and in northern California and 906 genes showed differentiation be- 2014, the only historical to modern comparison where we tween Avalon and southern California. A total of 807 genes have even sample sizes and high quality sequence for all indi- showed differentiation in both comparisons. Odorant recep- viduals. We found that there are fewer SNPs at a FDR of 1% in tor (OR) genes occur in the list of differentiated genes with six exonic regions and more in intronic regions than expected. OR genes differentiated between northern California and We found 306 genes with nonsynonymous SNPs in the Avalon and two OR genes differentiated between southern Humboldt comparison, nine of which have five or more non- California and Avalon. A further four OR genes are differen- synonymous SNPs. Most of these genes were uncharacter- tiated in both comparisons. Vitellogenin is differentiated in ized, but hemolectin, which is involved in clotting, and both population comparisons (fig. 4) and its receptor, yolkless, LOC551016, which is described as midasin-like, are in this is differentiated between Avalon and northern California. set. The overall set of genes with nonsynonymous SNPs is These genes may represent local adaptation between very different between the Humboldt and the modern Genome Biol. Evol. 10(2):458–472 doi:10.1093/gbe/evy007 Advance Access publication January 15, 2018 467 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/458/4810442 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Cridland et al. GBE FIG.5.—Differentiated SNPs in Vitellogenin in modern population comparisons. Vg is a minus strand gene and there are no SNPs observed in the 500 bp upstream of the 5 end. Californian comparison with only 29 genes in both compar- of feral A. mellifera and potential local adaptation of feral isons, and 456 that are in only one comparison. Of the genes populations. that do occur in both comparisons, one is yolkless,a vitello- genin receptor. Modern Populations Populations examined in this study were collected away from Discussion agricultural areas, but cannot be confirmed to be feral bees. Here, we document patterns of genetic diversity of Apis mel- Examining the patterns of genetic differences between pop- lifera populations over temporal and spatial scales in the state ulations can indicate the relatedness of these individuals to of California. We document substantial genetic diversity in the managed stocks. Regardless, these bees serve as indicators of A. mellifera populations that are likely the result of 1) changes the genetic diversity of A. mellifera found in regions through- in beekeeping practices over time, 2) the introduction of out California. We identified three main geographically and Africanized bees to California in 1994, 3) the introduction genetically distinct groups of honey bee populations in of the Varroa mite to California in 1987, and 4) the life history California. First, there are populations of northern 468 Genome Biol. Evol. 10(2):458–472 doi:10.1093/gbe/evy007 Advance Access publication January 15, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/458/4810442 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Genome Sequencing of Museum Specimens GBE Californian bees that closely resemble managed stocks. Prior 1994 provide a source of genetic variation that is currently to the introduction of Varroa mites, feral colonies appear to lacking in northern California, though Africanized bees ap- be generally stable over time, but the demographics of feral pear to be spreading northwards (Kono and Kohn 2015) bees can still change relatively rapidly. The density of feral and our analyses support this conclusion. colonies has been measured in Ithaca, NY was found to We find that African populations contribute between change slowly as established colonies are likely to survive 37.6% and 73.3% of the ancestry of modern southern through the winter (79% survival; Seeley 1978), which is Californian populations. In addition, because of concerns of similar to the survival of feral populations in Arizona (80%; Africanization, modern beekeepers in southern California Taber 1979), where winters are milder. Baum et al. (2005), have altered breeding practices to reduce or eliminate the who examined feral A. mellifera colonies in Texas between introduction of Africanized alleles into domesticated lineages. 1993 and 2000, found that, although colony turnover was 5– This makes it very likely that the southern Californian bees are 30% per year, over that time Africanized colonies increased genuine feral bees. Compared with northern Californian pop- from an initial colony discovered in 1993 to 80.3% of colonies ulations, modern southern Californian populations also ap- Africanized in 2000 (Baum et al. 2005). This baseline rate of pear to have more substantial genetic contributions from turnover gives us a good idea of what we should expect in western European (M) lineages, a pattern which has been terms of the genetic relationships between feral populations observed previously (Whitfield et al. 2006; Rangel et al. and nearby managed bees if we assume that managed pop- 2016). This may be due to a decrease in gene flow between ulations continue to generate new feral colonies and feral managed and feral populations in southern Californian, selec- colonies are not experiencing any additional pressures such tive pressures favoring western European alleles, or hitchhik- as Varroa mites or African introgression. ing of M alleles with African spread. Further studies are The introduction of the Varroa mite has substantially al- required to determine what processes may be at work here. tered the population dynamics of feral colonies. After being This produces a population with high diversity and contribu- first discovered in California in 1987, Varroa reduced the feral tions from both western and eastern European bees as well as colony population of California to less than one fifth its pre- Africanized strains and is similar to the population dynamics of vious size by 1994 (Kraus and Page 1995). Weexpectthatthis feral bees in the Yucatan in between 1993 and 1998 as substantial reduction in the feral bee population would result Africanization of the local populations occurred (Clarke in reduced genetic diversity of modern feral bees both by et al. 2002). Our results support the hypothesis that popula- wiping out most feral colonies that existed prior to Varroa, tions from southern California are experiencing rapid ongoing and also increasing the turnover of feral colonies by reducing changes, possibly due to local adaptation and admixture from the mean feral colony survival time, as was shown in Kraus multiple lineages. and Page (1995). These predictions match well with the pat- The increased genetic diversity we observed in mainland terns we documented in northern California. Historical pop- southern Californian populations bears directly on some of ulations exhibited a shift over time, concurrent with changing the current issues facing A. mellifera populations in beekeeping preferences from M lineage to C lineage popula- California. Genetic diversity has been shown to be an impor- tions (Crane 1999), and resembling more closely domesti- tant factor in colony health, with genetically diverse colonies cated populations, derived mostly from European (C) exhibiting higher disease resistance (Tarpy and Seeley 2006), lineages with some European (M) ancestry as well. Many of increased thermoregulation capacity (Jones et al. 2004), in- these populations were taken from the area at the northern creased diversity of worker responses to environmental stimuli end of California’s Central Valley, an area with high mite in- (Page et al. 1995), and higher productivity and fitness (Mattila festation rates (Kraus and Page 1995). and Seeley 2007). Managed bee colonies have recently been Our analysis based on southern Californian bees revealed found to have higher levels of genetic variation than honey that these populations derive substantial genetic ancestry bee populations in the eastern or western Europe (Harpur from African lineages and this influx of new alleles result in et al. 2012). However, those studies were based on colonies populations that are much more genetically diverse than that exhibited much lower genetic diversity than the northern Californian bees. This diversity is consistent with Africanized populations we documented here. We find that observations in Texas of hybrid swarms that were found in within California, Africanized southern Californian popula- the contact zone of African and European derived bees, post tions exhibit more genetic diversity than managed or modern Africanization (Pinto et al. 2005). The southern Californian northern Californian populations, but the managed and mod- mainland populations appear to be primarily influenced by ern northern Californian populations show similar levels of different factors than northern Californian populations, genetic diversity both with respect to each other and to west- though the timing of the Varroa mite outbreak and the intro- ern and eastern European bees. duction of Africanized bees into the United States may have Africanized bees have been found to be more resistant to played a role in favoring Africanized bee spread (Pinto et al. Varroa mites, likely due to a higher grooming rate and shorter 2004). The introduction of Africanized bees to California in period between egg laying and capping of the larvae Genome Biol. Evol. 10(2):458–472 doi:10.1093/gbe/evy007 Advance Access publication January 15, 2018 469 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/458/4810442 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Cridland et al. GBE (Camazine 1986; Moritz and Mautz 1990; Emsen et al. 2012), western European (M) like bees to primarily eastern European Moreover, it has been hypothesized that Africanization of (C) like bees, but we suspect that it was driven by the reduc- Californian honey bees may assist feral A. mellifera popula- tion of the feral population in the mid 1990s and changes in tions in surviving Varroa by introducing these beneficial breeding practices over that time period. behaviors into southern Californian populations (Kraus and Page 1995). Environmental variables may also play an impor- Differentiated Genes tant role in facilitating mite infestations. Medina-Flores et al. In addition to the broad-scale changes in the genetic compo- (2014) observed an association between drier conditions and sition of modern honey bee populations of California we iden- lower infestation rates in Mexico. Our study is congruent with tified a substantial number of genes that are differentiated this observation, and may help explain the observed patterns between populations. Some genes, like vitellogenin, show of higher genetic diversity in southern California where envi- differentiation between multiple population comparisons. ronmental conditions are similar. Further studies should inves- These genes constitute a candidate sets for both local adap- tigate the relative role of genetic background and tation to the various environmental conditions along environmental conditions with respect to honey bee resis- California, and future studies may investigate their potential tance Varroa and other health threats. role in resistance to Varroa or other diseases faced by honey Our analysis uncovered the existence of a third modern bee populations. We also find some interesting candidates in population in Catalina Island, which is located off the coast the comparison of the 1966 and modern Humboldt popula- of California near Los Angeles. This population is genetically tion, though these populations are composed of only a few distinct from the mainland populations and bears more re- individuals each and the differentiation observed may be due semblance to the western European (M) lineages and the to sampling. The gene spatzle 2,a toll receptor, is highly dif- Humboldt County 1966 population. This population is likely ferentiated between the Humboldt 1966 and Humboldt 2015 to have been on the island for some time, possibly since 1890, populations. It was found to be downregulated in high Varroa when A. mellifera was introduced on nearby Santa Cruz Island treatment groups in a previous study of the effects of de- (Wenner and Thorp 1994) and there are no managed bees on formed wing virus and Varroa on honey bees (Ryabov et al. the island. Additionally, this population is likely to be free of 2014). Other functions such as neuronal development, neu- Varroa mites. The Varroa mite was never found on Santa Cruz ronal sensitivity,and olfaction have been implicated in studies Island until it was introduced in an effort to biologically control of the genetic basis of tolerance to Varroa (Zakar et al. 2014). honey bees (Wenner and Thorp 1994) and is therefore also unlikely to have made it to Santa Catalina Island. Moreover, the genetic differences between the mainland population and Conclusions Santa Catalina Island indicate that there is little to no gene Our study provides an in-depth examination of temporal and flow between the populations, and thus little opportunity for spatial changes in honey bee populations in California. We the transfer of mites. This population provides an interesting identified substantial genetic differences across these two look into past honey bee populations and has likely under- dimensions in populations that are likely the result of both gone adaptation to the local environment. biotic and abiotic factors. Our analysis uncovers the potential genetic signatures that beekeeping practices and the intro- duction of Varroa have had on the changes experienced by Temporal Changes bee populations over time and suggests that southern The temporal changes we observed were dramatically differ- Californian populations may harbor higher genetic variation ent between southern California and northern California. The that couldbe harnessedtoimprove the healthof A. mellifera first Africanized bees to reach California in 1994 have, over and its important role in the agriculture of California. the past couple of decades, contributed substantially to the population structure. In contrast, historical populations of A. Supplementary Material mellifera prior to the 1970s in northern California differ from more modern northern Californian populations primarily due Supplementary data are available at Genome Biology and to a reduction in the contribution from western European (M) Evolution online. lineages and a concomitant increase in the contribution from eastern European (C) lineages. Modern northern Californian Acknowledgments populations in the Central Valley (Davis, Stanislaus, and Stebbins), are also very similar to modern domesticated We would like to thank Dr. Michael Mesler for help in collect- bees, F ¼ 0.059–0.065, which likely reflects frequent dis- ing bees. We would also like to thank C.F. Koehnen and Sons, ST persal of domesticated individuals into feral populations as Phil Hofland, and Glenda Wooten for providing drones. This escaped queens and drones establish new colonies. It is diffi- work was supported by the Berkeley Initiative for Global cult to assess the precise timing of this change from more Change Biology; the Office of the Vice Chancellor for 470 Genome Biol. Evol. 10(2):458–472 doi:10.1093/gbe/evy007 Advance Access publication January 15, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/458/4810442 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Genome Sequencing of Museum Specimens GBE Kerr EW. 1967. The history of the introduction of African bees in Brazil. Research at UC Berkeley; the Gordon and Betty Moore South Afric Bee J. 39:33–35. Foundation [grant number GBMF2983], the Vincent Coates Klein AM, et al. 2007. Importance of pollinators in changing landscapes for Genome Sequencing Facility; National Institutes of Health, S10 world crops. Proc Biol Sci. 274(1608):303–313. Instrumentation Grants [grant numbers S10RR029668 and Kono Y, Kohn JR. 2015. Range and frequency of Africanized honey bees in S10RR027303]; the USDA National Institute of Food and Calfornia (USA). PLoS One 10(9):e0137407. Kraus B, Page RE Jr. 1995. Effect of Varroa jacobsoni (Mesostigmata: Agriculture, Hatch Project [grant number CA-B-INS-0087-H]. Varroidae) on feral Apis mellifera (Hymenopter: Apidae) in California. Environ Entomol. 24(6):1473–1480. Literature Cited Langmead B, Salzberg S. 2012. Fast gapped-read alignment with Bowtie Allendorf FW, Hohenlohe PA, Luikart G. 2010. Genomics and the future of 2. Nat Methods 9(4):357–359. Li H, 1000 Genome Project Data Processing Subgroup, et al. 2009. The conservation genetics. Nat Rev Genet. 11(10):697–709. Alexander DH, Novembre J, Lange K. 2009. Fast model-based estimation sequence alignement/map (SAM) format and SAMtools. Bioinformatics 25(16):2078–2079. of ancestry in unrelated individuals. Genome Res. 19(9):1655–1664. Losey JE, Vaughn M. 2006. The economic value of ecological services Baum KA, Rubink WL, Pinto MA, Coulson RN. 2005. Spatial and temporal distribution and nest site characteristics of feral honey bee provided by insects. BioScience 56:311–323. Luikart G, England PR, Tallmon D, Jordan S, Taberlet P. 2003. The power (Hymenoptera: Apidae) colonies in a coastal prairie landscape. Environ Entomol. 34(3):610–618. and promise of population genomics: from genotyping to genome typing. Nat Rev Genet. 4(12):981–994. Benjamini Y, Hochberg Y. 1995. Controlling the false discovery rate: a Magnus R, Szalanski AL. 2010. Genetic evidence for honey bees (Apis practical and powerful approach to multiple testing. J Roy Stat Soc B 1:289–300. mellifera L.) of middle eastern lineage in the United States. Sociobiology 55:285–296. Brandt A, Gorenflo A, Siede R, Meixner M, Bu ¨ chler R. 2016. The neon- icotinoids thiacloprid, imidacloprid, and clothianidin affect the immu- Mattila HR, Seeley TD. 2007. Genetic diversity in honey bee colonies nocompetence of honey bees (Apis mellifera L.). J Ins Phys 86:40–47. enhances productivity and fitness. Science 317(5836):362–364. Moritz RFA, Mautz D. 1990. Development of Varroa jacobsoni in colonies Buchler R, Berg S, Le Conte Y. 2010. Breeding for resistance to Varroa destructor in Europe. Apidologie 41(3):393–408. of Apis mellifera capensis and Apis mellifera carnica. Apidologie 21(1):53–58. Camazine S. 1986. Differential reproduction of the mite Varroa jacobsoni (Mesostigmata: Varroidae), on Africanized and European bees Medina-Flores CA, Guzman-Novoa E, Hamiduzzaman MM, (Hymenoptera: Apidae). Ann Entomol Soc Am. 79(5):801–803. Are ´ chiga-Flores CF, Lopez-Carlos MA. 2014. Africanized honey bees (Apis mellifera) have low infestation levels of the mite Clarke KE, Rinderer TE, Franck P, Quezada-Euan JG, Oldroyd BP. 2002. The Africanization of honeybees (Apis mellifera L.) of the Yucatan: a study Varroa destructor in different ecological regions in Mexico. Genet Mol Res. 13(3):7282–7293. of a massive hybridization event across time. Evolution 56:1462–1474. Crane E. (1999) The world history of beekeeping and honey hunting. UK: Page RE, Robinson GE, Fondrk MK, Nasr ME. 1995. Effects of worker Tayor and Francis. p. 1–682. genotypic diversity on honey bee colony development and behavior (Apis mellifera L.). Behav Ecol Sociobiol. 36(6):387. Cridland JM, Tsutsui ND, Ramirez SR. 2017. The complex demographic history and evolutionary origin of the western honey bee, Apis melli- Patterson N, et al. 2012. Ancient admixture in human history. Genetics. doi: 10.1534/genetics.112.145037. fera. Genome Biol Evol. doi: 10.1093/gbe/evx009. ~ Pinto AM, Rubink WL, Coulson RN, Patton JC, Johnston JS. 2004. Temporal De la Rua P, Jaffe ´ R, Dall’Olio R, Munoz I, Serrano J. 2009. Biodiversity, conservation and current threats to European honeybees. Apidologie pattern of Africanization in a feral honeybee population from Texas inferred from mitochondrial DNA. Evolution 58(5):1047–1055. 40(3):263–284. Emsen B, Petukhova T, Guzman-Novoa E. 2012. Factors Limiting the Pinto AM, Rubink WL, Patton JC, Coulson RN, Johnston JS. 2005. Growth of Varroa destructor Populations in Selected Honey Bee Africanization in the United States. Genetics 170(4):1653–1665. Rangel J, et al. 2016. Africanization of a feral honey bee (Apis mellifera) (Apis mellifera L.) Colonies. J Anim Vet Adv. 11:4519–4525. Gutenkunst RN, Hernandez RD, Williamson SH, Bustamante CD. 2009. population in south Texas: does a decade make a difference?. Ecol Evol. 6(7):2158–2169. Inferring the joint demographic history of multiple populations from multidimensional SNP data. PLoS Genet. 5(10):e1000695. Rinderer TE, Collins AM, Tucker KW. 1985. Honey production and under- Guzman-Novoa E, Page RE. 1999. Selective breeding of honey bees lying nectar harvesting activities of Africanized and European honey- bees. J Apicult Res. 24(3):161–167. (Hymenoptera: Apidae) in Africanized areas. J Econ Entomol. 92(3):521–525. Ruttner F. (1988) The western honeybee Apis mellifera L. classification and natural history. Biography and taxonomy of honeybees. Berlin Harpur BA, et al. 2014. Population genomics of the honey bee reveals strong signatures of positive selection on worker traits. Proc Natl Acad Heidelberg GmbH: Springer-Verlag. P. 163–175. Sci U S A. 111(7):2614–2619. Ryabov EV, et al. 2014. A virulent strain of deformed wing virus (DWV) of honeybees (Apis mellifera) prevails after Varroa destructor-mediated, or Harpur BA, MinaeI S, Kent CF, Zayed A. 2012. Management increases genetic diversity of honey bees via admixture. Mol Ecol. in vitro, transmission. PLoS Pathogn. doi: 10.1371/journal.ppat.1004230. Sanchez-Bayo F, et al. 2016. Are bee diseases linked to pesticides?—A 21(18):4414–4421. brief review. Environ Inter. 89–90:7–11. Huang DW, Sherman BT, Lempicki RA. 2009a. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Schiff NM, Sheppard WS, Loper GM, Shimanuki H. 1994. Genetic diversity of feral honey bee (Hymenoptera: Apidae) populations in the southern Protoc. 4:44–57. Huang DW, Sherman BT, Lempicki RA. 2009b. Bioinformatics enrichment United States. Ann Entomol Soc Am. 87(6):842–848. Seeley TD. 1978. Life history strategy of the honey bee, Apis mellifera. tools: paths towards the comprehensive functional analysis of large Oecologia 32(1):109–118. gene lists. Nucleic Acids Res. 37:1–13. Jones JC, Myerscough MR, Graham S, Oldroyd BP. 2004. Honey bee nest Sheppard WS, Huettel MD. 1988. Biochemical genetic markers, intraspe- cific variation, and population genetics of the honey bee, Apis melli- thermoregulation: diversity promotes stability. Science 305(5682):402–404. fera. Environ Entomol. 22(1):183–189. Genome Biol. Evol. 10(2):458–472 doi:10.1093/gbe/evy007 Advance Access publication January 15, 2018 471 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/458/4810442 by Ed 'DeepDyve' Gillespie user on 16 March 2018 Cridland et al. GBE Stinchcombe JR, Hoekstra HE. 2008. Combining population genomics and Watkins LH. 1968. California’s first honey bees. Am Bee J. quantitative genetics: finding the genes underlying ecologically impor- 108:190–191. tant traits. Heredity 100(2):158–170. Wenner AM, Thorp RW. 1994. Removal of feral honey bee (Apis mellifera) Taber S III. 1979. A population of feral honey bee colonies. Am Bee J. colonies from Santa Cruz Island. In: Halvorson WL, Meander GL, edi- 119:842–847. tors. The Fourth California Islands Symposium: update on the status of Tarpy DR, et al. 2010. Mating frequencies of Africanized honey bees in the resources. Santa Barbara (CA): Santa Barbara Museum of Natural south western USA. J Apicult Res. 49(4):302–310. History. p. 513–522. Tarpy DR, Seeley TD. 2006. Lower disease infections in honeybee (Apis Winston ML, Taylor OR, Otis GW. 1983. Some differences between tem- mellifera) colonies headed by polyandrous vs monoandrous queens. perate European and tropical African and South American honeybees. Naturwissenschaften 93(4):195. Bee World 64(1):12–21. US Department of Agriculture, National Agricultural Statistics Service. Whitfield CW, et al. 2006. Thrice Out of Africa: Ancient and Recent 2016. Honey Bee Colonies (US Department of Agriculture, Expansions of the Honey Bee, Apis mellifera.Science Washington, DC). Available from: http://usda.mannlib.cornell.edu/ 314(5799):642–645. MannUsda/viewDocumentInfo.do? documentID¼1943; last accessed Zakar E, Javor A, Kusza S. 2014. Genetic basis of tolerance to Varroa February 2017. destructor in honey bees (Apis mellifera L.). Insect Soc. vanEngelsdorp D, Meixner MD. 2010. A historical review of managed 61(3):207–215. honey bee populations in Europe and the United States and the factors Associate editor: Brandon Gaut that may affect them. J Invertebr Pathol. 103:S80–S95. Wallberg, A, et al. 2014. A worldwide survey of genome sequence vari- ation provides insight into the evolutionary history of the honeybee Apis mellifera. Nat. Genet. 46:1081–1088. 472 Genome Biol. Evol. 10(2):458–472 doi:10.1093/gbe/evy007 Advance Access publication January 15, 2018 Downloaded from https://academic.oup.com/gbe/article-abstract/10/2/458/4810442 by Ed 'DeepDyve' Gillespie user on 16 March 2018

Journal

Genome Biology and EvolutionOxford University Press

Published: Feb 1, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Unlimited reading

Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere.

Stay up to date

Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates.

Organize your research

It’s easy to organize your research with our built-in tools.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve Freelancer

DeepDyve Pro

Price
FREE
$49/month

$360/year
Save searches from Google Scholar, PubMed
Create lists to organize your research
Export lists, citations
Access to DeepDyve database
Abstract access only
Unlimited access to over
18 million full-text articles
Print
20 pages/month
PDF Discount
20% off