Whole Genome Sequencing and Assembly of the Asian Honey Bee Apis dorsata.

The Asian honey bee (Apis dorsata) is distinct from its more widely distributed cousin Apis mellifera by a few key characteristics. Most prominently, A. dorsata, nest in the open by forming a colony clustered around the honeycomb, whereas A. mellifera nest in concealed cavities. Additionally, the worker and reproductive castes are all of the same size in A. dorsata. In order to investigate these differences, we performed whole genome sequencing of A. dorsata using a hybrid Oxford Nanopore and Illumina approach. The 223 Mb genome has an N50 of 35 kb with the largest scaffold of 302 kb. We have found that there are many genes in the dorsata genome that are distinct from other hymenoptera and also large amounts of transposable elements, and we suggest some candidate genes for A. dorsata's exceptional level of defensive aggression.

Introduction

The Asian honey bee (Apis dorsata) is an important pollinator and source of honey throughout its natural range in Southeast and East Asia (Sheppard and Berlocher 1989; Thapa 2006; Chantawannakul et al. 2018). This species is distinct from the European honey bee (Apis mellifera) both morphologically and behaviorally. Unlike A. mellifera, which is now found throughout the world, A. dorsata is found only in Asia (Chantawannakul et al. 2018). Apis dorsata is much larger than A. mellifera, and A. dorsata workers are almost twice as long as A. mellifera workers. In addition, A. dorsata is distinguished by a lack of body size variation between castes as in other honey bees.

The flavor of A. dorsata honey differs from that of A. mellifera. Biochemical assays of honey have demonstrated variation in antibacterial activity, protein content, and glucose content (Iftikhar et al. 2011; Pattamayutanon et al. 2015). While much of this variation can be explained by the nectar source that bees use to produce the honey, honey characteristics also differ based on honeybee species.

From an anthropocentric point of view, the most striking difference between A. dorsata and other honey bees is in their nesting behavior. Apis dorsata do not nest in cavities but instead build exposed nests that hang from tree branches or cliffs. The exposed comb is covered at all times by a blanket of up to 100,000 worker bees. Unlike species that nest in enclosed spaces like tree trunks and rock crevices, A. dorsata colonies are unwilling to live in the traditional Langstroth hives used by commercial beekeepers (Koeniger et al. 2010). This distinct lifestyle, together with an exceptional level of defensive aggression (Breed et al. 2004), has limited the ability of humans to domesticate the species and to commercially produce its honey (Koeniger et al. 2010).

To fully understand the genetic basis of important Apis traits requires comparative data from a variety of Apis species. In this study, we combine short-read and long-read sequencing technologies to generate a high-quality assembly of the A. dorsata.

Materials and Methods

Sample Collection

Multiple A. dorsata drones were collected by on the campus of Chiang Mai University (Chiang Mai, Thailand). The drones were collected in the early evening when they visited artificial light sources near the University’s gymnasium. Specimens were sacrificed immediately. One exemplar was kept intact; each of the remaining specimens was cut into three sections and placed in RNAlater Stabilization Solution (Invitrogen). Tissues were held in RNAlater at 4C for at least 24 h to allow the solution to saturate the tissue, then shipped to the American Museum of Natural History (AMNH). The exemplar specimen was accessioned into the Ambrose Monell Cryo Collection at AMNH. The remaining specimens were sent to the Cold Spring Harbor Laboratories for DNA extraction and sequencing.

DNA Extraction and Sequencing

Individual heads and thoraces were chopped into small pieces and placed in 1.7 ml centrifuge tubes (∼25 mg of tissue per tube). DNA was isolated using the MagAttract HMW DNA kit (Qiagen), following the “Manual Purification of High-Molecular-Weight Genomic DNA from Fresh or Frozen Tissue” protocol. Several tubes were pooled after extraction.

Extracted DNA was sheared to 20 kb using a Covaris G-tube following manufacturer’s recommendations. The sheared DNA was prepared for Nanopore sequencing using the Ligation Sequencing 1D kit (SQK-LSK108). Briefly, 1–1.5 µg of fragmented DNA was end-repaired with the NEB FFPE repair kit, followed by end repair and A-tailing with the NEB Ultra II end-prep kit. Following an Ampure clean-up step, prepared fragments were ligated to ONT-specific adapters with the NEB blunt/TA master mix kit. After a final Ampure clean-up, the library was loaded on to a MinION R9.0 SpotON flowcell as per manufacturer’s instructions. The flowcell was sequenced with standard parameters for 2 days. The resulting sequence data were processed with the Albacore pipeline (Oxford Nanopore).

Illumina libraries were prepared with the Kapa hyper prep kit, following manufactures guidelines. Prepared libraries pooled and sequenced on an Illumina NextSeq 500. The sequencing was done with PE150 High Output format.

Genome Assembly and Annotation

The Nanopore and Illumina reads were assembled together using the SPAdes software (Antipov et al. 2016). Structural annotation of the assembled genome was performed with MAKER as described in Campbell et al. (2014). Briefly, in the first round, transcript sequences from a previous A. dorsata study (NCBI BioBroject PRJNA174631) were provided to MAKER as EST evidence. Protein evidence was obtained using protein sequences from closely related hymenopteran species from NCBI, including A.dorsata (PRJNA174631), A. mellifera (PRJNA477511), Apiscerana (PRJNA324433), Apisflorea (PRJNA45871), Polistes dominula (PRJNA307991), and Papiliocanadensis PRJNA301748). In addition, the full Swiss-Prot database was included in the analysis (Boutet et al. 2016). Transcript sequences from the non-dorsata species were provided to MAKER as alternative EST evidence. Repetitive elements were identified using (Chen 2004) (A.F.A. Smit).

Coding sequences were predicted in the MAKER pipeline with Exonerate (Slater and Birney 2005) and BLAST (Camacho et al. 2009). The output of the first MAKER run was used to train the gene predictors SNAP (Korf 2004) and AUGUSTUS (Keller et al. 2011). The trained models were then provided to MAKER and a new set of training models generated from the result used for another round of running MAKER. Another sets of SNAP and AUGUSTUS models were trained with the new outputs. The assembled contig sequences were used to train GeneMark-ES (Lomsadze et al. 2005). The SNAP, AUGUSTUS, and GeneMark-ES models were used in a final MAKER run, and the resulting predicted coding sequences used in all subsequent analyses. The completeness of the assembly and annotation were analyzed with BUSCO. Simão et al. (2015) using the Insecta data set (1,658 single-copy conserved genes) as a reference. We used the assembly-stats tool (https://github.com/sanger-pathogens/assembly-stats) to determine assembly statistics for the dorsata and other Apis genomes.

Functional Annotation

Annotation of the predicted genes was undertaken with a variety of tools. First, we conducted BLASTp (BLAST+ 2.8.1) searches against all Hymenoptera entries in the NCBI nr database (Wheeler et al. 2007) (default parameters except that we used the -taxidlist option to restrict the search to species within NCBI Taxonomy ID 7399 [Hymenoptera]; nr database accessed on 09/09/2019). For sequences that had no hit, we did a second BLASTp search against all sequences in the NCBI nr database. We used InterProScan 5 (Jones et al. 2014) to classify genes into protein families and identify functional domains(Mitchell et al. 2015). The results were imported into Blast2GO (Conesa et al. 2005), where we mapped them to Gene Ontology (GO) terms (Harris et al. 2004).

We also compared the content of the assembled A. dorsata genome with that of A. mellifera.

Ortholog Identification

We used OrthoPipe, a stand-alone pipeline version of OrthoDB 2.3.1 (Kriventseva et al. 2015), to classify the protein-coding genes from A. dorsata, A. mellifera, A. cerana, and A. florea into orthologous clusters. To evaluate the resulting clusters, we used the reference topology presented by Han et al. (2012).

Identification of Species-Specific A. dorsata Genes

Species-specific genes were narrowly defined as sequences that lacked orthologs in any other Apis species. A more stringent definition requires a lack of orthologs in any published data set (Tautz and Domazet-Lošo 2011), and we used BLASTp searches against NCBI’s nr database to test whether our narrowly defined species-specific genes were legitimate orphans.

Sequence Comparison of A. dorsata Genes to Others

To examine how genes from A. dorsata have diverged from those of the other examined Apis species, we generated multiple sequence alignments with MAFFT (v7.312) (Katoh and Standley 2013). All alignments were generated using “high accuracy” mode (–maxiterate 1000 –localpair). The resulting alignments were used to generate cladograms to visualize the relationship between different genes. We used FastTree (v 2.1.10) (Price et al. 2009) to approximately-maximum-likelihood phylogenetic trees, with the options –nome (to compute minimum-evolution bootstrap support values) and -slow (to conduct exhaustive searches).

Functional Enrichment Tests

Because raw frequency-based comparisons are difficult to interpret in evolutionary terms, we used enrichment tests to examine whether any gene functions were overrepresented in the genome of A. dorsata as compared with A. mellifera. We conducted Fisher’s exact tests (q-value < 0.01) to identify GO terms or InterPro signatures enriched in a test set of annotated genes relative to a reference set. Since GO annotations, which are assigned transitively via BLAST hits, can sometimes lead to incorrect functional annotation (Yon Rhee et al. 2008), we used InterPro signatures as the primary source of functional annotation and comparison.

Results and Discussion

Genome Assembly

Our assembly results were consistent with previous Apis genome sequencing projects (fig. 1): The total assembly size was 224 Mb (vs 225–230 Mb in other apid bees, Consortium 2006; Wang et al. 2013; Park et al. 2015), with a contig N50 of 28 kb (largest contig: 302 kb). This assembly puts the A. dorsata on the same scale as the A. florea and A. cerana genomes, which have contig N50s of 20–30 kb and largest contigs of over 300 kb. In contrast to the previous A. dorsata assembly from 2013 (PRJNA174631), we have an improved assembly from a lower level of sequencing coverage. None of these genomes compared with A. mellifera, an intensely studied model organism with a very high-quality contiguous genome assembly. BUSCO analysis of our assembly showed that the assembly was 98.9% complete. We identified a wide variety of repetitive elements (supplementary table 1, Supplementary Material online), and these comprised ∼7% of the genome. Using MAKEr, we identified 13,517 protein coding genes (vs 12,145–12,940 genes in other apids, Consortium 2006; Wang et al. 2013; Park et al. 2015).

. 1.

—(A) Phylogenetic relationship between the four Apis species used for ortholog identification, relative sizes of drones are shown. Inset bar charts indicate the fraction of genes from each species that were orthologous at different levels: Universal orthologs (present in all four species), species-specific genes (those found only in one species), and others. (B) Assembly statistics for Apis species.

More than 80% (11,264) of the predicted A. dorsata genes had matches to other Hymenopteran sequences. The remaining 2,253 genes had no hits to other hymenoptera (or only had hits to A. dorsata sequences). When we blasted these genes against the entire nr database, 182 of them had hits but 2,071 still had no hit. The results of InterProScan analyses were similar: 1,956 genes had no InterPro signature. The overlap between these annotations left us with 1,251 genes that lacked functional annotation.

We identified 3,451 single copy universal ortholog clusters that contained one gene from each of the four species included in ortholog analysis. Consistent with the idea that single copy universal orthologs represent highly conserved genes, every A. dorsata single copy universal ortholog had a hit to other Hymenoptera in the nr database. Seventy-three percent of the A. dorsata genes (9,900) had orthologs in at least one other species, versus 93% of the A. mellifera genes. In the specific comparison between A. mellifera and A. dorsata, 22% (3,034) of the A. dorsata genes had no ortholog in A. mellifera, and 2,177 of these had no ortholog in any of the included Apis species and no hit to any other hymenopteran in our BLAST searches (fig. 1).

Species-Specific A. dorsata Genes

Some 3,617 A. dorsata genes did not have orthologs in any of the included Apis species. We used BLAST searches against all hymenopteran sequences in the nr database to determine whether these were species-specific genes. We found that 1,424 of these genes had a BLAST hit to another hymenopteran species, and an additional 75 had a hit to a nonhymenopteran species. This left 2,118 genes that appear to be either pseudogenes or species-specific A. dorsata genes because they had no hit to any sequence in the nr database.

Divergence in 5-Hydroxytryptamine Receptor 2A (5-HT2a)

We were particularly interested in identifying genes that might contribute to A. dorsata’s behavioral divergence from other Apis species. We focused on receptors for the neurotransmitter serotonin (5-hydroxytryptamine, 5-HT) because of serotonin’s role as a neuromodulator (Turlejski 1996). The effects of serotonin depend on receptor-specific binding, and in Drosophila and other insects 5-HT receptors have been associated with variation in locomotion, feeding behavior, learning and memory, and aggression (Majeed et al. 2016; Huser et al. 2017).

We identified four 5-HT2A receptors in A. dorsata (fig. 2). Although three of these had orthologs in other Apis species, one was sufficiently divergent that it did not cluster with any other Apis genes. Because serotonin has a stimulatory effect on aggression in some insects (Bubak et al. 2019), this gene is an attractive candidate for future studies on the genetic basis of aggression in A. dorsata.

. 2.

—(A) Cladogram with bootstrap support values and (B) sequence alignment of putative 5-hydroxytryptamine (serotonin) receptor 2A (5-HT2A) genes from Apis dorsata and other Apis species. (C) Cladogram with bootstrap support values and (D) sequence alignment of A. dorsata species-specific genes with InterPro signature IPR001127 (Lactobacillus PTS carbohydrate transport). Lactobacillus apis genes are from the gut of Apis mellifera. Both cladograms (A and C) are based on full-length gene alignments; in the interests of space, only a representative section of the full alignment is shown (B and D).

Functional Enrichment

Single Copy Universal Orthologs versus All Other Genes

These highly conserved genes were enriched for functions related to development, including a growth factor responsible for cell proliferation and differentiation (IPR001111), and juvenile hormone binding proteins (JHBP; IPR010562), which regulates embryogenesis, larval development, and reproductive maturation in the adult forms (supplementary table 2, Supplementary Material online).

Species-Specific A. mellifera Genes versus All Other Genes

The A. mellifera genes that lacked A. dorsata homologs were enriched for signal transduction functions. Several enriched signatures related to small GTPases, which are critical components in cellular signal transduction pathways,

Species-Specific A. dorsata Genes versus All Other Genes

The InterPro terms that were enriched in the A. dorsata genes that lacked A. mellifera homologs fell into two main functional categories. First, there were 28 genes with Ty1/copia-like retrotransposon signatures. The “domestication” of transposable elements can lead to the emergence of species-specific genes by providing a source of biochemically active elements such as transcription factor-binding sites, and by generating genomic rearrangements (Toll-Riera et al. 2008; Kaessmann 2010; Lisch 2013; Göke and Ng 2016), and the overrepresentation of transposable element-related signatures in the genome of A. dorsata, suggests that the evolution of “new” genes may have played a role in the behavioral diversification of A. dorsata from other honeybees.

Second, there were 23 sequences that appear to be part of A. dorsata’s microbiome assemblage. These all had high-scoring hits to bacterial species in the nr database, and when aligned with bacterial sequences known to occur in Hymenopteran microbiomes showed high levels of similarity. The enriched InterPro signature IPR001127, a domain that is associated with the carbohydrate transport system in Lactobacillus bacteria, was found in six of the species-specific A. dorsata genes. We compared the A. dorsata sequences to Lactobacillus apis genes from the microbiome of A. mellifera (fig. 2).

In honeybees, the microbiome of worker bees is of especial importance. All worker gut microbiomes studied to date contain abundant carbohydrate-processing genes, but the exact components of the microbiome community show variation between honeybee strains and species (Engel et al. 2012). This variation suggests that some of the dietary differences between honeybees in different localities (in their use of pollen/nectar from different plant species) might be controlled by microbiome components (Disayathanoowat et al. 2012; Saraithong et al. 2015, 2017). Previous research has shown that Lactobacillus species, which were the most common top hits for our A. dorsata-only sequences, play a major role in carbohydrate digestion (Lee et al. 2015; Moran 2015). Thus, the specific microbiome components present in different colonies or species might determine their nutritional ecology and/or their ability to cope with dietary toxins and could explain the geographic distribution of A. dorsata as a case of coevolution between bees, bacteria, and flowering plants.

The microbiome of A. dorsata may also be involved in the distinctive flavor of the honey they produce. Lactobacillus species are thought to prevent the fermentation of stored honey by retarding the growth of fermenting yeasts (Olofsson and Vásquez 2008). Variation in the activity of Lactobacillus strains in honeybee microbiomes can affect the flavor of honey since flavor is partially dependent on lactic acid bacterial metabolites (Olofsson et al. 2016). Thus, the “terroir” of different honeys (Marchese and Flottum 2013) may result as much from the honeybee microbiome (Tajabadi et al. 2011) as from the particular locale where the honeybees forage.

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.

Click here for additional data file.