Genome of the small hive beetle (Aethina tumida, Coleoptera: Nitidulidae), a worldwide parasite of social bee colonies, provides insights into detoxification and herbivory

Abstract

Background: The small hive beetle (Aethina tumida; ATUMI) is an invasive parasite of bee colonies. ATUMI feeds on both fruits and bee nest products, facilitating its spread and increasing its impact on honey bees and other pollinators. We have sequenced and annotated the ATUMI genome, providing the first genomic resources for this species and for the Nitidulidae, a beetle family that is closely related to the extraordinarily species-rich clade of beetles known as the Phytophaga. ATUMI thus provides a contrasting view as a neighbor for one of the most successful known animal groups.

Results: We present a robust genome assembly and a gene set possessing 97.5% of the core proteins known from the holometabolous insects. The ATUMI genome encodes fewer enzymes for plant digestion than the genomes of wood-feeding beetles but nonetheless shows signs of broad metabolic plasticity. Gustatory receptors are few in number compared to other beetles, especially receptors with known sensitivity (in other beetles) to bitter substances. In contrast, several gene families implicated in detoxification of insecticides and adaptation to diverse dietary resources show increased copy numbers. The presence and diversity of homologs involved in detoxification differ substantially from the bee hosts of ATUMI.

Conclusions: Our results provide new insights into the genomic basis for local adaption and invasiveness in ATUMI and a blueprint for control strategies that target this pest without harming their honey bee hosts. A minimal set of gustatory receptors is consistent with the observation that, once a host colony is invaded, food resources are predictable. Unique detoxification pathways and pathway members can help identify which treatments might control this species even in the presence of honey bees, which are notoriously sensitive to pesticides.

Figure 1:

Aethina tumida (A) adult and (B) larva. Photos courtesy of Alex Wild Photography, used with permission.

Figure 2:

Phylogenetic tree and estimated completeness of the genomes of 11 insect species. (a) The phylogenetic tree was constructed on protein sequences of 181 single-copy orthologs shared among all 11 insect species. All nodes have 100% bootstrap support. AMELL and DMELA were used as outgroups. Branch lengths are shown for each node. (b) Completeness of official protein sets of each insect species were assessed by aligning to the Endopterygota sets of benchmarking universal single-copy orthologs (BUSCOs). For ATUMI, 97.5% of complete BUSCOs were found. (c) The pervasiveness of gene loss during endopterygote evolution. From the domain counts of lost BUSCOs, methyltransferase (MT), glycosyltransferase (GT), and leucine-rich repeats (LRR) are among the top 5% of total domains and are commonly lost from multiple species. Beta-transducin repeats (WD) and zinc finger (ZF) red boxes indicate that the gene is lost, while white boxes indicate that the gene is maintained in each species.

Figure 3:

Gene duplication events plotted against the average gene duplication event per gene. The protein sets of the 11 studied beetle species, as well as honey bee and fruit fly, were searched against the Endopterygota BUSCO set using Basic Local Alignment Search Tool (BLAST). Redundant proteins (including recent paralogs and those with known alternative splicing) were used to quantify the average number of duplication events per gene in each species.

Figure 4:

Glycoside hydrolase (GH) family copy numbers identified from beetle genomes. Genes coding for GHs were identified using Pfam domain assignments [24], and genome assemblies and coding gene predictions were obtained from NCBI (GenBank accession numbers: GCA_000390285.1 ALB, GCA_000355655.1 DPOND, GCA_001412225.1 NIC, GCA_001443705.1 ORY, GCA_000002335.3 TCAS) with the exception of CBB, which was downloaded from [25]. Families are color coded from green to red based on their relative abundance (total count/total number of GH genes), with red representing GH families that are highly abundant (≥25% of the total GH genes) and green representing GH families of lesser abundance (≤0.01%). Notably, the GH profiles of ATUMI and TCAST (neither of which feed on living plant material) differ strongly from the GH profiles of the phytophagous beetles, even though they all belong to the same infraorder, suggesting that diet, in part, might be driving the differences in GH family members and copy numbers. ALB = Asian longhorned beetle (A. glabripennis); CBB = Coffee berry borer (H. hampei); DPOND = Mountain pine beetle (D. ponderosae); NIC = burying beetle (N. vespilloides); ORY = scarab beetle (O. borbonicus); SHB = small hive beetle; and TCAS = red flour beetle (TCAST).

Figure 5:

Maximum likelihood cladogram for gustatory receptor genes from three coleopteran species. The small hive beetle, Aethina tumida (Atum; green labels/lines), the Asian longhorned beetle, Anoplophora glabripennis (Agla; red labels/lines), and the red flour beetle, Tribolium castaneum (Tcas; blue labels/lines). Individual genes are labeled with species identifier and GenBank accession number. Scale bar for branch lengths represents 0.1 amino acid substitutions per site. Ring around cladogram indicates gene families coded for perceiving bitter (yellow) and sweet (pink) tastants.

Figure 6:

Maximum likelihood phylogenetic tree of ABC proteins from ATUMI (At), TCAST (Tc), and DMELA (Dm). ATUMI genes are marked in blue, TCAST in green, and DMELA in purple. ABC subfamilies are indicated with colored lines to the right of the tree. Names for DMELA proteins were taken from Flybase [52] and include the Flybase number for reference. TCAST names were taken from the two articles in which the genes were identified [50, 51], with the NCBI Refseq accession number provided for reference. ATUMI names were generated for this article by combining the subfamily of the identified sequence with the scaffold on which the encoding gene may be found. If multiple ABC genes of a particular subfamily were found on the same scaffold, the sequences were given an additional letter designation based on their relative location, reading left to right on the scaffold as shown in WebApollo. For reference, the scaffold number and base coordinates for the gene have also been included.

Figure 7:

Maximum likelihood phylogenetic tree of glutathione-S-transferase (GST) proteins. The bootstrap consensus tree inferred from 1,000 replicates is taken to represent the evolutionary history of the taxa A. tumida (ATUMI) in green, A. mellifera (AMELL) in black, D. melanogaster (DMELA) in blue, and T. castaneum (TCAST) in red, identified manually using the Uniprot and Pfam databases. Branches corresponding to partitions reproduced in less than 50% bootstrap replicates are collapsed. Initial tree(s) for the heuristic search were obtained automatically by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using a JTT model and then selecting the topology with superior log likelihood value. All positions with less than 95% site coverage were eliminated. The tree was annotated and visualized with the iToL web tool (itol.embl.de/) [58].

Figure 8:

Maximum likelihood phylogenetic tree of the cytochrome P450 detoxification system. The bootstrap consensus tree inferred from 1,000 replicates is taken to represent the evolutionary history of the taxa A. tumida (ATUMI) in green, A. mellifera (AMELL) in black, D. melanogaster (DMELA) in blue, and T. castaneum (TCAST) in red, identified manually using the Uniprot and Pfam databases. Branches corresponding to partitions reproduced in less than 50% bootstrap replicates are collapsed. Initial tree(s) for the heuristic search was obtained automatically by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using a JTT model and then selecting the topology with superior log likelihood value. All positions with less than 95% site coverage were eliminated. P450s are clustered to CYP2, CYP3, CYP4, and mitochondrial clans. The tree was annotated and visualized with the iToL web tool (itol.embl.de/) [58].

Figure 9:

Maximum likelihood phylogenetic tree of carboxylesterase (COE) genes. The maximum likelihood bootstrap consensus tree (1,000 replicates) showing the relationships among COE genes from the genomes of A. tumida (ATUMI) in green, A. mellifera (AMELL) in black, D. melanogaster (DMELA) in blue, and T. castaneum (TCAST) in red, identified manually using the Uniprot and Pfam databases. Branches corresponding to partitions recovered in less than 50% of bootstrap replicates are collapsed. Starting tree(s) for the heuristic search was obtained automatically using neighbor-joining and BioNJ algorithms applied to a matrix of pairwise distances estimated using a JTT model and then selecting the topology with the superior log likelihood value. All positions with less than 95% site coverage were eliminated. The phylogenetically distinct clusters were named according to established nomenclature for COE genes [12]. The tree was annotated and visualized with the iToL web tool [24].