Genus-Wide Characterization of Bumblebee Genomes Provides Insights into Their Evolution and Variation in Ecological and Behavioral Traits

Abstract Bumblebees are a diverse group of globally important pollinators in natural ecosystems and for agricultural food production. With both eusocial and solitary life-cycle phases, and some social parasite species, they are especially interesting models to understand social evolution, behavior, and ecology. Reports of many species in decline point to pathogen transmission, habitat loss, pesticide usage, and global climate change, as interconnected causes. These threats to bumblebee diversity make our reliance on a handful of well-studied species for agricultural pollination particularly precarious. To broadly sample bumblebee genomic and phenotypic diversity, we de novo sequenced and assembled the genomes of 17 species, representing all 15 subgenera, producing the first genus-wide quantification of genetic and genomic variation potentially underlying key ecological and behavioral traits. The species phylogeny resolves subgenera relationships, whereas incomplete lineage sorting likely drives high levels of gene tree discordance. Five chromosome-level assemblies show a stable 18-chromosome karyotype, with major rearrangements creating 25 chromosomes in social parasites. Differential transposable element activity drives changes in genome sizes, with putative domestications of repetitive sequences influencing gene coding and regulatory potential. Dynamically evolving gene families and signatures of positive selection point to genus-wide variation in processes linked to foraging, diet and metabolism, immunity and detoxification, as well as adaptations for life at high altitudes. Our study reveals how bumblebee genes and genomes have evolved across the Bombus phylogeny and identifies variations potentially linked to key ecological and behavioral traits of these important pollinators.


Introduction
Bumblebees (Hymenoptera: Apidae) are a group of pollinating insects comprising the genus Bombus, which are economically important for crop pollination (Velthuis and Van Doorn 2006;Garibaldi et al. 2013;Martin et al. 2019). Bumblebees are also ecologically important pollinators, serving as the sole or predominant pollinators of many wild plants (Fontaine et al. 2005;Goulson et al. 2008). They are particularly charismatic social insects that exhibit complex behaviors such as learning through observation (Alem et al. 2016) and damaging leaves to stimulate earlier flowering (Pashalidou et al. 2020). Global and local environmental changes have resulted in some species declining in range and abundance and others remaining stable or even increasing (Cameron et al. 2011;Bartomeus et al. 2013;Koch et al. 2015;Cameron and Sadd 2020). Decline in bumblebee abundance and distribution resulting from habitat loss, pathogen transmission, climate change, and agrochemical exposure is threatening pollination services to both wild plants and crops, raising concerns for bumblebees, the plant species they service, food security, and ecosystem stability (Grixti et al. 2009;Williams and Osborne 2009;Potts et al. 2010;Goulson et al. 2015;Cameron and Sadd 2020;Soroye et al. 2020).
Bumblebees comprise $250 extant species classified into 15 subgenera (Williams 1998;Cameron et al. 2007;Williams et al. 2018). The initial diversification of Bombus lineages occurred $25-40 Ma, near the Eocene-Oligocene boundary $34 Ma (Williams 1998;Hines 2008). Bumblebees display considerable interspecific diversity in morphology, color patterning, food preference, and pathogen incidence and exhibit diverse life histories and ecologies (Williams 1994;Sikora and Kelm 2012;Persson et al. 2015;Arbetman et al. 2017;Cameron and Sadd 2020). Members of the subgenus Mendacibombus, the sister group to all other extant bumblebees, are high-elevation specialists with distributions centered on the Qinghai-Tibetan plateau (Williams et al. 2018). Species in the subgenus Psithyrus exhibit social parasitism; they do not have a worker caste, and they feed on food collected by workers of their host species . Bumblebees are distributed across the globe, from Greenland to the Amazon Basin and from sea level to altitudes of 5,640 m in the Himalayas, where they occupy diverse habitats, from alpine meadows to lowland tropical forest (Williams 1985;Williams et al. 2018). Much remains to be learned about bumblebees. For example, little is known about the underlying genetic and genomic variation that gives rise to these diverse phenotypes, including their differential responses to changing environments.
About half of the $250 extant species, representing 14 of the 15 Bombus subgenera, are found in China, making it a hotspot of bumblebee species richness (Williams 1998;Williams et al. 2017). To broadly sample the genomic and phenotypic diversity of bumblebees, we selected representative species from China for whole-genome sequencing based on their phylogeny, ecology, behavior, geography, and specimen availability. To complete subgenus sampling, we additionally selected Bombus polaris from the subgenus Alpinobombus, which is endemic to Arctic/subarctic regions. In total, we performed de novo sequencing and assembly of the genomes of 17 bumblebee species, representing all 15 subgenera within the genus Bombus. Integrating these data sets with two previously published bumblebee genomes, we performed comparative analyses of genome structures, genome contents, and gene evolutionary dynamics across the phylogeny. Our results characterize patterns of molecular and genomic evolution across the Bombus phylogeny and provide the first genus-wide quantification of genetic and genomic variation potentially underlying key ecoethological traits.

High-Quality Genomic Resources for All 15 Bombus Subgenera
Sequencing and assembly strategies resulted in high-quality genomic resources with 12 scaffold-level and five chromosome-level genome assemblies (table 1). Criteria including phylogenetic position, species traits, and geographic distribution were applied to select species for whole-genome sequencing from across the genus. For the five species for which sufficient samples could be collected, highthroughput chromatin conformation capture (Hi-C) (Belton et al. 2012) was used to produce chromosome-level genome assemblies (table 1). A total of 17 species were selected (supplementary table S1 and fig. S1, Supplementary Material online), which span all 15 subgenera of the genus Bombus (Williams et al. 2008). Among these, two species (Bombus superbus and B. waltoni) are from Mendacibombus, the earliest split in the Bombus phylogeny; four species (B. superbus, B. waltoni, B. skorikovi, and B. difficillimus) inhabit high elevations (>4,000 m above sea level); two species (B. turneri and B. skorikovi) exhibit social parasitism; three species (B. pyrosoma, B. picipes, and B. superbus) are endemic to China; and one species (B. polaris) is endemic to Arctic/subarctic regions (Williams et al. 2019). In addition, species traits including range size, tongue length, parasite incidence, and decline status vary across the selected species (Williams 1994;Arbetman et al. 2017;Cameron and Sadd 2020).
Sequencing and assembly strategies included generating two Illumina sequencing data sets for each species: 1) overlapping paired-end reads (2Â 250 bp) from one small-insert fragment library using one single haploid drone per species (insert size: 400 or 450 bp) and 2) paired-end reads (2Â 150 bp) from four large-insert jump libraries using 3-5 individuals per species (insert sizes: 4, 6, 8, and 10 kb, respectively; supplementary table S2, Supplementary Material online). Whole-genome overlapping paired-end reads from fragment libraries were assembled into continuous sequences (contigs) using the software DISCOVAR de novo (Love et al. 2016), then scaffolded with reads from jump libraries using the software BESST (Sahlin et al. 2014). The resulting assemblies have a mean contig N50 of 325 kb, ranging up to 579 kb for B. breviceps; the mean scaffold N50 is 4.0 Mb, ranging up to 6.9 Mb for B. superbus (table 1). Genome assembly quality in terms of expected gene content was evaluated by Benchmarking Universal Single-Copy Ortholog (BUSCO) analysis (Waterhouse et al. 2018

Genome-Scale Phylogeny of Bumblebees
The species-level molecular phylogeny ( fig. 1A) estimated from maximum-likelihood analysis with IQ-TREE (Minh, Schmidt, et al. 2020) is largely consistent with previously inferred phylogenetic relationships of the 15 subgenera based on five genes (Cameron et al. 2007;Williams et al. 2008), showing only two topological differences. The results support previous conclusions that 1) subgenus Mendacibobus (labeled Md in fig. 1A) is the sister group to all the other subgenera and 2) species of Psithyrus (labeled Ps in purple in fig. 1A and genome assembly size of each sequenced species; fraction of TEs (brown) in each genome. (B) Bar plots show total gene counts for each bumblebee partitioned according to their orthology profiles, from ancient genes found across bumblebees to lineage-restricted and speciesspecific genes. (C, D) The contribution of TE and CDS to genome size variation across bumblebees, respectively. Differences in the total content of TEs (C) and CDS (D) of the 19 genomes relative to that of Bombus superbus (which has the smallest genome assembly size) are plotted against their genome size differences (relative to that of B. superbus).
Genus-Wide Characterization of Bumblebee Genomes . doi:10.1093/molbev/msaa240 et al. 2014) and tomatoes (Pease et al. 2016) and have been attributed to a variety of sources, such as ILS and introgression (Maddison 1997). A lack of informative sites, only 24%, compared with 47% in a similar data set of 25 drosophilids (Da Lage et al. 2019), possibly due to the relatively recent diversification of bumblebees (Hines 2008), may also cause discordance. Gene and site concordance factor (sCF) analysis (Minh, Hahn, et al. 2020) was employed to quantify the amount of discordance between gene trees and the IQ-TREE species tree (node labels in fig. 1A). For each node in the IQ-TREE species tree, gene concordance factors (gCFs) reflect the percentage of gene trees that contain that node as defined by its descendant taxa, and sCFs reflect the percentage of informative sites that support that node via parsimony. On average across the Bombus phylogeny, nodes in the IQ-TREE species tree show a gCF of 38.4%, indicating that on average a node is present in only two-fifths of gene trees.

Major Genomic Rearrangements in Social Parasites
The five Hi-C genome assemblies indicate that four of the five subgenera have 18 chromosomes ( fig. 2A and C; supplementary fig. S10A and B, Supplementary Material online), consistent with previous karyotypic analysis that inferred the ancestral chromosome number is 18 (Owen et al. 1995). However, the social parasite bumblebee, B. turneri, subgenus Psithyrus, has 25 chromosomes ( fig. 2B), consistent with previous cytological work (Owen 1983). Despite the higher chromosome number, its genome size is within the range of other bumblebees ( fig. 1A and table 1). Investigating macrosynteny relationships between B. turneri and the other species with chromosomal-level assemblies revealed three major processes that explain how a 25-chromosome karyotype was derived from the ancestral karyotype of 18 chromosomes. First, some chromosomes descended, structurally unchanged, from ancestral chromosomes (e.g., chromosome 5; fig. 2D in blue). Second, some originated by fission of an ancestral chromosome (e.g., 11 and 25 of B. turneri originated by the fission of ancestral chromosome 11; fig. 2D in red). Lastly, some are derived from fusions of two or more ancestral chromosome segments (e.g., B. turneri MBE genes (supplementary table S10, Supplementary Material online). In total, 352 of these genes are universal single-copy orthologs across the 19 bumblebees whose overall dN/dS values are all <1 (supplementary table S11, Supplementary Material online), suggesting long-term functional constraints. One case of a putative ancient and maintained chimeric TEgene fusion involves a gene with single-copy orthologs across the 19 bumblebees where the C-terminus of the proteins match the sequence of a reverse transcriptase of an R1 retrotransposon (supplementary fig. S15, Supplementary Material online). Aligned reads from RNA-sequencing data continue with similar coverage levels into the putatively TE-derived region at the 3 0 -end of the gene, supporting the prediction and the expression of the full chimera. TE activity has therefore contributed to the evolution of the bumblebee proteincoding gene repertoire. In addition, there are thousands of TEs located within 1 kb of a gene in each species (supplementary

Gene Content Evolution Reflects Foraging and Diet Diversity
Orthology delineation results indicate that a majority of genes are found in one or more copies in nearly all species across bumblebees ( fig. 1B). These include 53 orthologous groups specific to the Bombus genus, which are present in all 19 bumblebees but absent in all four honeybees ( fig. 1B; supplementary table S13, Supplementary Material online), and may play roles in lineage-specific traits. Functional annotation suggests that five of these Bombus-specific genes are associated with protein metabolism and transport (supplementary table S13, Supplementary Material online), potentially linked to the higher protein content of pollen collected by bumblebees than honeybees (Leonhardt and Blüthgen 2012) or the importance of proteins for insect diapause, which is a critical step in the bumblebee life cycle (Denlinger 2002;Colgan et al. 2011). Orthologous groups with the broadest species representation are functionally enriched for core biological processes such as protein transport, signal transduction (e.g., Wnt pathway), (de)ubiquitination, and cytoskeleton organization (supplementary table S14, Supplementary Material online). In contrast, those with sparse or lineage-restricted Genus-Wide Characterization of Bumblebee Genomes . doi:10.1093/molbev/msaa240 MBE species representation are enriched for processes including olfactory and gustatory perception, amino acid biosynthesis, and oxidation-reduction (supplementary table S14, Supplementary Material online). On average, 465 speciesspecific genes (those without an ortholog in any other species) were identified in each bumblebee species (range 137-767) (supplementary table S15, Supplementary Material online), which may contribute to species-specific traits but whose functional roles remain to be explored.
Turnover analysis (gains and losses) of gene repertoires across the Bombus phylogeny (15 species, one per subgenus) using CAFE v3.0 (Han et al. 2013) identified expansions and contractions among 13,828 gene families and quantified variations in gene turnover rates across species (supplementary fig. S16, Supplementary Material online). After error correction, the overall rate of gene turnover in Bombus genomes is 0.0036/gene/My, similar to an analysis of 18 anopheline species and 25 drosophilids (supplementary table S16, Supplementary Material online) (Neafsey et al. 2015;Da Lage et al. 2019). However, these genus-specific gene turnover rates are 2-3 times higher than order-wide rates, which average 0.0011 (supplementary table S16, Supplementary Material online) (Thomas et al. 2020), possibly due to the denser sampling in genus-level studies that allow more events to be captured. Gene gain and loss events, along with the number of rapidly evolving gene families, are summarized for each species (supplementary table S17, Supplementary Material online), with a total of 3,797 rapidly changing gene families. The most dynamic gene families are enriched for processes including smell and taste perception, chitin metabolism, microtubule-based movement, and methylation (supplementary table S18, Supplementary Material online). Complementary analysis using three measures of gene copy number variation also identifies these processes as enriched among the most variable gene families, in contrast to the most stable that are involved in processes related to translation, adhesion, and transport (supplementary table S19, Supplementary Material online). In terms of protein domain copy number evolution, the most highly variable genes are those with protein-protein interaction mediating F-box domains, putatively DNA-binding SAP motifs, and phosphate-transferring guanylate kinases (supplementary table S20, Supplementary Material online).

Stable Intron-Exon Structures with Abundant Stop-Codon Readthrough
Protein-coding potential analysis using B. terrestris as the reference species identified 851 candidate readthrough stop codons (supplementary fig. S17 and table S21, Supplementary Material online), that is, where translation likely continues through stop codons to produce extended protein isoforms. Coding potential was assessed using PhyloCSF ) on whole-genome alignments of all 19 bumblebees and four honeybees. The false discovery rate was estimated using enrichment for the TGA-C stopcodon context, which is favored in readthrough genes, to infer that no more than 30% of the 200 highest-scoring candidates are false positives, and that at least 306 of our 851 candidates undergo functional readthrough. Although readthrough is rare beyond Pancrustacea, hundreds of Drosophila and Anopheles genes undergo readthrough Dunn et al. 2013;Jungreis et al. 2016;Rajput et al. 2019) and our whole-genome-alignment-based results support the prediction ) that insect species have abundant stop-codon readthrough. In contrast, intron-exon boundaries within bumblebee genes are relatively stable. Examining evolutionary histories of intron gains and losses revealed few changes, representing only 3-4% of ancestral intron sites, with more gains than losses (supplementary fig.  S18 and table S22, Supplementary Material online), unlike drosophilids and anophelines where losses dominate (Neafsey et al. 2015), suggesting that bumblebee gene structure has remained relatively stable over the 34 My since their last common ancestor.

Divergence and Selective Constraints of Protein-Coding Genes
Bumblebee genes with elevated sequence divergence and/or relaxed constraints include processes related to smell perception, chitin metabolism, RNA processing, DNA repair, and oxidation-reduction ( fig. 3). Measures of evolutionary rate (amino acid sequence divergence measured as the mean of normalized interspecies ortholog protein sequence identities) and selective constraint (dN/dS) showed similar trends among different functional categories of genes. Most genes are strongly constrained, with median estimates of dN/dS much lower than one. Assignment of GO terms and InterPro domains is usually biased toward slower-evolving,

Codon Usage Bias Driven by at Content
Analysis of codon usage bias showed no evidence for selection on optimal codons, in contrast to drosophilids but similar to anophelines (Vicario et al. 2007;Neafsey et al. 2015). Instead, codon usage bias in bumblebees seems to be driven mainly by AT content, consistent with previous reports in Hymenoptera (Behura and Severson 2012). Optimal codons were estimated in each species and correlation coefficients were computed between relative synonymous codon usage and effective number of codons per gene. All species have a similar preference and intensity of preference; for each amino acid, there was a consistently highly preferred codon and often a secondarily preferred

Evolution of Genes Associated with Bumblebee Ecoethology
Many ecological and environmental factors-for example, shortage of food, pathogen emergence, pesticide exposure, and climate change-are contributing to the overall decline of bumblebees worldwide Goulson et al. 2015;Cameron and Sadd 2020). To begin to explore the complement of genes likely to be involved in bumblebee interactions with their environment, we examined the evolution of gene families associated with their ecology and life histories. Sampling across the Bombus genus enabled the first survey of natural gene repertoire diversity of such families that are likely to be important for bumblebee adaptability and success.

Chemosensory Receptor Diversity
Chemosensation plays a critical role in locating food and nests, communicating with nestmates, and identifying other environmental cues (Ayasse and Jarau 2014   . For example, all species have two genes encoding Gram-negative bacteria-binding proteins, whereas peptidoglycan-recognition proteins are more variable with between four and six gene copies. Comparing dN/dS ratios between immune genes and all single-copy orthologous genes in bumblebees showed that immune genes exhibit slightly higher dN/dS ratios (P ¼ 0.04, Wilcoxon rank sum test), and among immune genes, recognition and signaling genes have higher dN/dS ratios than effector genes ( fig. 4B). In addition, despite widespread evidence for purifying selection, a total of 52 immune genes showed evidence of positive selection in a subset of bumblebee species (supplementary table  S30, Supplementary Material online).
Genes Involved in High-Elevation Adaptation Bombus superbus, B. waltoni, B. difficillimus, and B. skorikovi are four species collected at elevations >4,000 m that represent three subgenera ( fig. 1). No genes show signatures of positive selection in all high-elevation species but none of the low-elevation species. However, six genes show evidence of positive selection in species representing two of the three high-elevation subgenera, but none of the low-elevation species (supplementary table S31, Supplementary Material online). One encodes CPAMD8, which is involved in eye development (Cheong et al. 2016). As bumblebees detect flowers visually (Meyer-Rochow 2019), signatures of selection might be related to fine tuning eye development for optimal foraging in high-altitude light conditions. Three genes encode a histone deacetylase, synaptotagmin-12, and a heterogeneous nuclear ribonucleoprotein, which are involved in maintaining muscle integrity and keeping "flight state," which is critical for undertaking long-distance food searching (Liu et

Sex Determination
Evolutionary analysis of sex-determination genes in bumblebees and related species indicated that all bumblebee genomes share a duplicated copy of feminizer (fem), named fem 1 (fig. 4C). Compared with fem, fem 1 shows a higher level of divergence among bumblebees (fem Bombus dN/dS ¼ 0.24; fem 1 Bombus dN/dS ¼ 0.77; fig. 4C). These ratios are close to the range observed for Apis, in which fem has evolved under purifying selection and the paralogous gene complementary sex determiner (csd) has evolved by neofunctionalization ( fig. 4C) (Hasselmann et al. 2008). A hypothesis branch-site testing framework (RELAX; Wertheim et al. 2015) identifies evidence for relaxation of selection in fem 1 Bombus compared with fem Bombus (P < 0.001, LR ¼ 36.34). Moreover, the spurious action of diversifying selection on branches was predominantly found in fem 1 Bombus (fig. 4C). A mixed effect model of evolution (Murrell et al. 2012) was applied to identify individual sites that were subject to episodic diversifying selection, and at least 15 sites (P < 0.05) were found to be under positive selection, with some being located in known motifs (supplementary fig. S26, Supplementary Material online). The results of these selection analyses suggest that both fem and fem 1 contribute to the bumblebee sexdetermination pathway. For the transformer 2 (tra-2) gene, consistent amino acid changes between Bombus and Apis were found within the RNA recognition domain (supplementary fig. S27, Supplementary Material online), supporting a previous hypothesis of a regulatory modification between honeybees and bumblebees (Biewer et al. 2015). Boxplots showing dN/dS ratios for different categories of immune genes and all single-copy genes in bumblebee (All genes). Elevated dN/dS ratios among immune-related genes are driven by higher ratios for genes involved in recognition and signaling processes. Notched boxes show medians of orthologous group values with the limits of the upper and lower quartiles. (C). The evolutionary history of fem genes of bees including their paralogs fem1 in bumblebees (Bombus) and csd in honeybees (Apis). Global nonsynonymous to synonymous rate ratio (x) was calculated for fem Bombus (reference, blue) and fem1 Bombus (test, red), including a branch-site testing framework with model fitting and Likelihood Ratio Tests, showing evidence for relaxation of selection in fem1 Bombus (P < 0.001, LR ¼ 36.34). Spurious actions of diversifying selection on branches predominantly found in fem1 Bombus are marked in red. For comparison, x for fem and csd in Apis is given, known as striking example of neofunctionalization.

Discussion
Comparative analysis of multiple genomes in a phylogenetic framework substantially improves the precision and sensitivity of evolutionary inference and provides robust results identifying stable and dynamic features. In this study, we performed comparative analyses of genome structures and contents, as well as global and family-targeted gene evolutionary dynamics across the phylogeny of the genus Bombus, using 17 annotated de novo assemblies and two previously published genomes.
Many attributes of bumblebee genomes are highly conserved across species. For example, overall genome size and genome structure, the number of protein-coding genes and noncoding RNAs, gene intron-exon structures, and the pattern of codon usage are all very similar across these 19 genomes. However, other aspects of genome biology are dynamically evolving. TEs are a major contributor to genome size variation ( fig. 1) as well as a potential source of coding and regulatory sequences (supplementary tables S10-S12, Supplementary Material online). Differential gene gain and loss also contribute to gene content variation across bumblebees and lead to lineage-specific gene repertoires ( fig. 4A; supplementary fig. S25 and table S17, Supplementary Material online). Finally, for genes shared by all species, the action of positive selection is different across species (supplementary tables S26, S28, S30, and S31, Supplementary Material online), which can lead to gene functional divergence possibly reflecting key ecoethological differences.
An exception to the otherwise overall conserved genome structure is the set of species in the subgenus Psithyrus. These bumblebees exhibit social parasitism; they do not have a worker caste, and it is not necessary for them to forage for nectar and pollen to provision developing larvae . Originally, this subgenus was argued to be a separate genus due to distinct behavior and higher chromosome number; however, subsequent phylogenetic analysis placed Psithyrus within the genus Bombus (Cameron et al. 2007;Williams et al. 2008). Here, based on a much larger genomic data set, we confirm that species in the subgenus Psithyrus (labeled Ps in purple in fig. 1A) form a monophyletic group within the genus Bombus. In addition, we show that, although Psithyrus species have an increased chromosome number, their genome sizes are within the range of those of the other bumblebees ( fig. 1A), and their 25 chromosomes reflect a mix of fission, fusion, and retention of the 18 ancestral bumblebee chromosomes ( fig. 2; supplementary fig. S10, Supplementary Material online). Chromosome rearrangements (e.g., fissions, fusions, and inversions) have been posited to play roles in speciation (Ayala and Coluzzi 2005) and thus may explain the diversification and social parasitic behavior of Psithyrus. In addition to genome structure variation, we identified a net loss of 11 OR genes in the common ancestor of Psithyrus species (fig. 4), which could be a cause or consequence of their socially parasitic behavior.
Bumblebee species exhibit different food preferences (Goulson and Darvill 2004;Sikora and Kelm 2012;Somme et al. 2015), but the genetic basis underlying such variation is unknown. Like in other insects, smell and taste are used to distinguish different food sources (Kunze and Gumbert 2001;Ruedenauer et al. 2015). In this study, we found out that genes involved in olfactory and gustatory perception are among the fastest-evolving gene categories, both in copy number variation and in sequence divergence ( fig. 3; supplementary fig. S20 and tables S18 and S19, Supplementary Material online). Therefore, the dynamic evolution of genes involved in olfactory and gustatory perception likely relates to different food preferences, improved understanding of which could inform the use of new species in agricultural settings. The importance of chemical perception for social Hymenoptera in mating (Ayasse et al. 2001), determination of reproductive status (Monnin 2006), and recognition of kin (Page and Breed 1987) could also contribute to rapid evolution of genes in this category.
Bumblebees exhibit rich morphology differences across species including mandible, labrum, tibia, and basitarsus structures, as well as patterns of keels, ridges, and grooves formed by the cuticle (Williams 1994), and they show species-specific responses to insecticides (Baron et al. 2017). Chitin is a major component of the insect cuticle and peritrophic matrix, and chitin metabolic processes are related to morphogenesis, resistance to insecticides, and the tolerance of toxins in food (Barbehenn 2001;Merzendorfer and Zimoch 2003;Zhu et al. 2016;Erlandson et al. 2019). Genes related to chitin metabolism are also among the fastest-evolving functional categories in bumblebees, both in copy number variation and in sequence divergence ( fig. 3; supplementary fig. S20 and tables S18 and S19, Supplementary Material online). These variable patterns of chitin-related gene evolution potentially underlie observed differences in exoskeleton morphological characters and insecticide resistance, but not color pattern variation, which is determined primarily by hair pigmentation (Tian et al. 2019). Across bumblebee genomes some of the fastest-evolving genes are also related to processes including protein glycosylation, methylation, proteolysis, and tRNA aminoacylation for protein translation ( fig. 3; supplementary fig. S20 and tables S18 and S19, Supplementary Material online). Protein glycosylation is involved in multiple physiological processes including growth, development, circadian rhythms, immunity, and fertility (Walski et al. 2017). tRNA aminoacylation for protein translation process is involved in response to the changing environment (Pan 2013).
Some genes that are not among the fastest-evolving categories-for example, immune and detoxification genes, which are involved in the interaction of bumblebees with external environments-show differential patterns of positive selection in subsets of species, which can lead to gene functional divergence. The positive diversifying selection found in some detoxification genes of a subset of species could also contribute to species differences in response to insecticide exposure (Baron et al. 2017 MBE these species, where selective pressures encountered in different ecological niches may vary among species. These differences could affect susceptibility to emerging and re-emerging infectious diseases and explain observed species-specific differences in contemporary pathogen prevalence (Cameron et al. 2016;Cameron and Sadd 2020). Taken together, identification of the fastest-evolving genes and those showing patterns of differential positive selection reveals substantial genetic variation across bumblebee genera. Future experimental investigations will be required to determine how the identified genetic variation is linked to specific differences in traits such as food preference, morphogenesis, insecticide and pathogen resistance, and the response to changing environments.
In addition to our discoveries regarding protein-coding genes, we found that TE-related sequences likely contribute to the variation of coding and regulatory repertoires ( fig. 1; supplementary tables S10-S12, Supplementary Material online). Compared with non-Mendacibombus bumblebees, Mendacibombus species have smaller genomes ( fig. 1) and relatively narrow geographical distributions (Williams et al. 2016). Considering TEs are the major determinant of genome size difference, with evidence that they were potentially domesticated in bumblebee genomes (supplementary fig. S15 and tables S10-S12, Supplementary Material online), TEs may be implicated in the dispersal of non-Mendacibombus species across the globe, as they have been in other taxa (Casacuberta and Gonz alez 2013;Baduel et al. 2019;Schrader and Schmitz 2019).
More recent range expansions or contractions are driven, at least in part, by global climate change. In response to a warming climate, there is evidence that the ranges of some bumblebees are being pushed northward or to higher elevations (Ploquin et al. 2013;Kerr et al. 2015;Biella et al. 2017;Soroye et al. 2020). The sequenced genomes of species collected at high-elevation sites (>4,000 m) and others collected at low elevations (<2,000 m) ( fig. 1) represent high-quality genomic resources for investigating genes involved in highelevation adaptation. For example, population-genomics studies facilitated by the A. cerana reference assembly identified candidate high-altitude adaptation loci in that species (Montero-Mendieta et al. 2019) and the B. impatiens genome was used to identify climate adaptation across latitude and altitude in two montane North American bumblebee species (Jackson et al. 2020). We identified genes showing signs of positive selection in at least two subgenera of high-elevation species but not in any of the low-elevation species (supplementary table S31, Supplementary Material online). These include genes putatively involved in eye development, flight muscle integrity, and metabolism, highlighting the importance of successful food searching in high-elevation habitats where food is scarce. Signatures of positive selection in neuromuscular genes mirror findings from the population genomic study on the two North American montane species (Jackson et al. 2020). Beyond specific genes, comparing high-and lowelevation species showed a consistent pattern of faster genome-wide evolutionary rates in those occupying lower elevations (Lin et al. 2019). Exploring these trends and genes further and identifying additional genomic features linked to life at high altitudes will help to understand differential successes of bumblebee species in a changing world.

Conclusions
We have produced highly complete and accurate genome assemblies of 17 bumblebee species, including representatives from all 15 Bombus subgenera. Our genus-wide comparative analysis of bumblebee genomes revealed how genome structures, genome contents, and gene evolutionary dynamics vary across bumblebees and identified genetic variations that may underlie species trait differences in foraging, diet and metabolism, morphology and insecticide resistance, immunity and detoxification, as well as adaptations for life at high altitudes. Our work provides genomic resources that capture genetic and phenotypic variation, which should advance our understanding of bumblebee success and help identify potential threats. These resources form a foundation for future research, including resequencing and population-genomics studies for functional gene positioning and cloning, which will inform the use of bumblebees in agriculture, as well as the design of strategies to prevent the decline of these important pollinators.

Materials and Methods
For detailed methods, please see supplementary Materials and Methods, Supplementary Material online. Genomic DNA purified from one single haploid drone of each species was used to generate one "fragment" library to produce overlapping paired-end shotgun reads (2Â 250 bp). Genomic DNA purified from three to five individuals of each species was used to generate four "jumping" libraries (insert sizes: 4, 6, 8, and 10 kb, respectively). For each species, the 250-bp overlapping paired-end shotgun reads from the fragment library were assembled using the software DISCOVAR de novo (Love et al. 2016) to produce contiguous sequences (contigs). Shotgun reads from jumping libraries were used to scaffold the contigs using the software BESST (Sahlin et al. 2014). Hi-C sequencing libraries were generated as described previously (Belton et al. 2012;Lin et al. 2018), and the 3D-DNA pipeline (Dudchenko et al. 2017) was applied to assemble the scaffold sequences to chromosome level. The contiguity of the genome assemblies was compared with other genomic features (supplementary fig. S28, Supplementary Material online). The completeness of the genome assemblies was evaluated using BUSCO v3 (Waterhouse et al. 2018). TEs were de novo identified by LTRharvest (Ellinghaus et al. 2008), MGEScan-non-LTR (Rho and Tang 2009), and RepeatScout (Price et al. 2005). Protein-coding genes were annotated with the MAKER pipeline (Cantarel et al. 2008), based on ab initio gene predictions, transcript evidence (supplementary table S33, Supplementary Material online), and homologous protein evidence, and assessed with BUSCO v3 (Waterhouse et al. 2018). Orthologous groups were delineated by using the OrthoDB software (Kriventseva et al. 2015), which was also used to compute evolutionary rates. IQ-TREE (Minh, Schmidt, et al. 2020) and ASTRAL  were employed for Genus-Wide Characterization of Bumblebee Genomes . doi:10.1093/molbev/msaa240 MBE phylogeny analysis. Testing for introgression employed the D statistic (Huson et al. 2005), which follows the same logic as the ABBA-BABA site patterns used to calculate D-statistics, but uses tree topologies instead of alignment sites. MCScanX (Wang et al. 2012) was used to perform gene synteny analysis. Global gene family evolution was estimated with CAFE v3.0 (Han et al. 2013). PAML 4 (Yang 2007) was used to detect positive selection and calculate dN/dS ratios across alignments of single-copy orthologs. This was followed by aBSREL  analysis for gene families of special interest (i.e., chemosensory genes, detoxification genes, immune genes, sex-determination genes, and piRNA pathway genes) and for genes identified in the initial PAML scans as potentially showing signatures of adaptation to high elevation. Sex-determination genes were additionally analyzed using RELAX ) and mixed effect model of evolution (Murrell et al. 2012). Intron history analyses were performed using Malin (Cs} uros 2008). Whole-genome alignments were produced using Cactus (Paten et al. 2011). PhyloCSF ) was used to study stop-codon readthrough. Codon usage bias was analyzed as described previously (Vicario et al. 2007).

Supplementary Material
Supplementary data are available at Molecular Biology and Evolution online.