Short reads from honey bee (Apis sp.) sequencing projects reflect microbial associate diversity

Screening for microbes in Apis short read sequencing libraries

A graphical overview of our screening process can be found in Fig. 1. We compiled two sets of reference sequences to be used as databases for subsequent screens (Fig. 1). First, for 15 common Apis associated symbionts and pathogens, we used one short signature sequence each instead of complete genomes, in order to reduce the computational expense of all following steps. These were of slowly evolved ribosomal RNA genes to allow a range of diversity to be recovered through sequence similarity to the bait. Second, for 11 common viral associates, we compiled complete genomes. As all of these were very small (∼10,000 bp), this was unlikely to significantly impact computational run times. Fasta files of compiled database sequences are available under https://github.com/gerthmicha/symbiont-sra.

Next, we searched for honey bee sequencing projects in NCBI’s short read archive, using the search term ‘Apis’, and created lists of the hits for DNA and RNA sequencing projects, which were subsequently processed separately (Fig. 1). The great majority of projects were from Apis mellifera samples, but a few RNA sequencing projects from Apis cerana and Apis florea were also included. We will from here on refer to all of those samples as ‘Apis’. From the list of DNA sequencing projects, we excluded metagenome and amplicon projects, as these are specifically aimed at microbiome taxa. At the time of the search (October 2015), 18 projects matched these criteria. We downloaded all short read libraries associated with these projects (993 in total, Table S1) and mapped all reads of each of the libraries to the ‘symbionts & pathogens’ reference database using NextGenMap version 0.4.12 (Sedlazeck, Rescheneder & Von Haeseler, 2013). If at least 1,000 reads of a library were aligned to one or more sequence baits, we extracted the matching reads and assembled them using SPAdes version 3.7 (Bankevich et al., 2012). Contigs resulting from this assembly were then subject to taxonomic annotation via BLAST+ (Camacho et al., 2009) searches against a local copy of the NCBI ‘nt’ database, and the Blobtools package (Kumar et al., 2013).

For the RNA dataset, we removed any library from projects investigating the effects of viruses or Varroa infections, as Varroa treatment usually results in viral infections (Allen & Ball, 1996). Next, the reads of all remaining RNA sequencing libraries (492 libraries from 58 projects in total, Table S2) were mapped against the viral database as described above, and the number of mapped reads was counted using samtools version 1.4 (Li et al., 2009). Mapped reads were then assembled using IDBA-UD (Peng et al., 2012) and MEGAHIT version 1.1.1 (Li et al., 2015), which were designed to handle highly uneven sequencing depths, as expected from viral data (Visser et al., 2016). The ‘better’ assembly from these two methods was chosen by means of viral contig length comparison, i.e., the assembly that created the single longest contig was kept. This was more appropriate than choosing the highest N50 value, as most libraries contained only a single viral associate. If the two assembly methods produced overlapping contigs, which were identical in the overlap, the contigs were merged. A detailed description of all steps used to screen for microbial associates can be found under https://github.com/gerthmicha/symbiont-sra.

Assessing Lactobacillus diversity in Apis DNA sequencing libraries

Since our initial screen of DNA sequencing libraries yielded a high number of hits to various Lactobacillus species, we repeated the entire procedure using a database of 620 16S bait sequences from Lactobacillus only. These sequences were taken from a previously compiled dataset of Lactobacilli associated with Apis, other Hymenoptera, and other Lactobacillus sequences retrieved from public databases (McFrederick et al., 2013). All hits shorter than 250 bp were discarded, and remaining contigs were combined with the reference sequences. We used SSU-ALIGN version 0.1 (Nawrocki, 2009) to align and mask this dataset based on conserved secondary structure. Original and masked alignments are available from https://github.com/gerthmicha/symbiont-sra. A maximum likelihood phylogeny was reconstructed from the complete 16S alignment (740 sequences in total) using IQ-TREE version 1.3.10 (Nguyen et al., 2015) with automated model selection and 1,000 ultrafast bootstraps (Minh, Nguyen & Von Haeseler, 2013) to assess node support. The resulting tree was visualized using the online tool Evolview (He et al., 2016). Furthermore, as an approximate measure for the number of Lactobacillus OTUs recovered with our approach, we used the average neighbor clustering algorithm as implemented in mothur version 1.34.4 (Schloss et al., 2009).

Although our aim was not to recover all, but only the highly covered symbiont data from honey bee short reads, we wanted to test if our screening approach yields comparable results to more commonly used metagenomic approaches. To this end, we screened the reads of a metagenomic dataset created from the pooled DNA of 150 honeybee worker hindguts (Engel, Martinson & Moran, 2012; ∼43 M 150 bp paired-end reads, SRA accession: SRR5237156) for Lactobacillus in the same way as described above. We found 6 different Lactobacillus 16S sequences, all within the Firm-4 and Firm-5 Lactobacillus groups (Fig. S1). This was in agreement to the results obtained from taxonomic profiling approaches performed by Engel, Martinson & Moran (2012) and thus confirmed the general effectiveness of our approach (Fig. S1).

Reconstruction of symbiont genomes from Apis DNA sequencing libraries

Next, we aimed to validate that whole symbiont genomes can in principle be recovered from Apis sequencing projects. To this end, we chose one Apis mellifera intermissa sequencing library (SRR1046114, ∼85.5 M 100 bp paired-end reads) that contained ‘contamination’ from two Lactobacillus strains (Lactobacillus kunkeei & Fructobacillus sp.). We performed a de novo assembly using all reads with MEGAHIT version 1.1.1 (Li et al., 2015). All resulting contigs of this assembly were taxonomically assigned to either L. kunkeei, Fructobacillus sp. or ‘other’ based on BLAST searches, GC distributions, and read coverage. Reads matching to contigs from either Lactobacillus strain were then separately re-assembled using SPAdes, and all contigs smaller than 500 bp discarded. Completeness and contamination of the novel draft genomes were assessed based on the presence of conserved marker genes using CheckM version 1.0.6 (Parks et al., 2015), and annotation performed with PROKKA version 1.12 (Seemann, 2014). The annotated draft genomes are available under under https://github.com/gerthmicha/symbiont-sra. To evaluate the evolutionary relationships of newly assembled genomes in a broader taxonomic context, we assessed their phylogenetic placement. Whole-genome datasets were compiled for both strains (13 L. kunkeei genomes, 9 Fructobacillus & Leuconostoc genomes altogether, Table S3). For each of the datasets, single copy orthologs were identified using OrthoFinder version 0.2.8 (Emms & Kelly, 2015). Recombining loci were identified by using the pairwise homoplasy index test (Bruen, Philippe & Bryant, 2006), and removed from subsequent analyses (window size = 20 amino acid positions, significance cutoff at 0.05). Using IQ-TREE, we performed maximum likelihood analysis of two final supermatrices (947 loci and 290,774 aa for the L. kunkeii dataset, 435 loci and 145,069 positions for the Fructobacillus/Leuconostoc dataset). Prior to this, best-fitting partitioning schemes and models were selected using the ‘greedy’ scheme implemented in IQ-TREE (Lanfear et al., 2012).

Using the same approach, we assembled and annotated a Spiroplasma melliferum genome from the Apis mellifera library SRR957082, (∼224.5 M 50 bp single end reads). Phylogenetic analysis was performed based on a dataset of 206 concatenated single copy genes (58,950 amino acid positions) shared among 17 Spiroplasma strains (Table S3). Furthermore, to assess synteny, the newly assembled draft genome was ordered against and aligned with other Spiroplasma melliferum genomes (one genome each of strains IPMB4A and KC3) using the progressiveMauve algorithm of Mauve development snapshot version 2015-02-13 (Darling et al. 2010).

Assessing viral diversity in Apis RNA sequencing libraries

All potential viral contigs resulting from the assembly steps outlined above were blasted against a local copy of NCBI’s ‘nt’ database. As most contigs were annotated as deformed wing virus (DWV), we reconstructed the phylogenetic relationships between the detected DWV strains. We used only contigs where at least half of the DWV genome was present (5,000 bp minimum length), and aligned all sequences using MAFFT. As a reference, we added complete genomic sequences for each of the three types of DWV (A, B or VDV-1–Varroa destructor virus 1, and C) and one for Kakugo virus (KV), which is a variant of DWV type A (Mordecai et al., 2016). The phylogeny was estimated via maximum likelihood with IQ-TREE as described above.

Exploring viral effects on Apis gene expression from RNA sequencing libraries

We observed that in some cases, libraries from a single project (i.e., NCBI’s BioProject) displayed pronounced differences in viral loads as measured in our screen. We thus hypothesized that these differences should also be detectable through differentially expressed host genes. To test this hypothesis, four projects were investigated (Table 1) in which viral loads were particularly uneven (proportion of viral reads across libraries 0%–37%, 0%–15%, 0%–53%, and 0%–32%, respectively—see Table S4 for details). Several studies have examined differential gene expression of honey bees in response to viral infection, and a set of genes commonly responding to pathogens was recently identified via meta-analyses (Doublet et al., 2017). We used the 20 genes that were consistently differentially expressed across multiple studies (first 20 ranks in the category ‘differentially-regulated’ from the Doublet et al., 2017 analysis). The relative expression levels of these genes were determined for all libraries of the selected projects by (1) quality trimming of all reads with cutadapt version 1.1.3 (Martin, 2011; options: -q 25,25 –minimum-length 50 –pair-filter=both) (2) mapping the reads with NextGenMap (min. identity 99%); (3) counting the reads with samtools; (4) normalizing read counts to reads per kilobase per million mapped reads (RPKM) and log 2 transforming this value, (5) allocating the libraries from each project to either ‘virus’ or ‘virus-free’ groups based on the proportion of viral reads in the libraries; (6) testing for differentially expressed genes between these groups using the qCML method and exact tests in the R package ‘edgeR’ (Robinson & Smyth, 2007; Robinson, McCarthy & Smyth, 2009; R Core Team, 2015). Genes were regarded as differentially expressed for P-values ≤ 0.05.

Apis DNA sequencing libraries can help to reconstruct the honey bee microbiome

We used bait sequences of Apis symbionts and pathogens to determine if microbial data can be retrieved from DNA sequencing projects targeting Apis (honey bees) and found evidence for the presence of these taxa in 11% of 993 Apis short read libraries. This measure of non-target ‘contamination’ can be considered as conservative, since our approach only reports relatively high levels of contamination (at least 1,000 reads per bait sequence). Our approach revealed that all common gut symbionts of honey bees are also present as ‘contamination’ in DNA sequencing libraries (Fig. 2). Furthermore, the relative frequency of each of the gut colonizers in focused study roughly corresponds to the frequencies at which we detected them. For example, the microbiome of healthy honey bees is dominated by Lactobacilli (Kwong & Moran, 2016), and this is also reflected in our results (Fig. 2, Table S5). Our findings demonstrate that a reasonable understanding of honey bee gut microbiome composition could be gained solely from non-target sequences produced as a by-product of honey bee sequencing projects.

When we next targeted our screen of DNA sequencing libraries only at Lactobacillus, our protocol detected 25 taxonomically different Lactobacillus strains. Phylogenetic reconstruction of Lactobacillus relationships based on 16S rRNA generally reflected the current understanding of this genus’ taxonomy (Felis & Dellaglio, 2007; Salvetti & Torriani, 2012), and revealed that most Lactobacilli known to be associated with honey bees are also present in Apis short read libraries. This includes Firm-4 and Firm-5 Lactobacilli, both of which are honey bee hindgut colonizers, and L. kunkeei, which is common in nectar and hive material, and sometimes found in honey bee crops (Kwong & Moran, 2016). Furthermore, we found Fructobacillus, which share an ecological niche with L. kunkeei, i.e., they are found in flowers, nectar, and in honey bee guts (Endo, Futagawa-Endo & Dicks, 2009; Endo & Salminen, 2013). Although not classified as such, recent phylogenomic evidence suggests that Fructobacillus (and the closely related Leuconostoc) are part of the Lactobacillus radiation (Sun et al., 2015). Here, we also infer Fructobacillus grouping within, rather than outside of Lactobacillus (Fig. 2A).

In addition to the gut symbionts, three common honey bee pathogens were detected with our approach: Nosema, Crithidia/Lotmaria, and Spiroplasma. Nosema are microsporidian gut parasites of various honey bee species, and while the sampling of our screen is not representative, this finding corroborates the recognition of Nosema as widespread pathogen of honey bee colonies worldwide (Nixon, 1982; Klee et al., 2007). Crithidia/Lotmaria (Trypanosomatidae), another gut pathogen of Apis and related bee species (Schwarz et al., 2015) was detected at an even higher frequency (Fig. 1, Table S5). We further found Spiroplasma melliferum in one of the investigated sequencing libraries. Spiroplasma are common symbiotic bacteria of many invertebrates (Duron et al., 2008) and have been connected to pathogenicity in honey bees (Clark, 1977). The bait sequences of all of these pathogens showed a high coverage in our screen, suggesting that novel genetic variants can be recovered from already available data, or from data that will become available as by-product of future honey bee sequencing projects.

In addition to the targeted microbes, a number of taxa that are likely not part of the natural Apis microbiome were detected. For example, we detected Aspergillus in several sequencing libraries that originated from museum material, which likely represents post mortem saprophytic growth. We also retrieved hits to plant sequences which might originate from co-amplified and sequenced pollen DNA (Fig. 2). This ‘false discovery’ illustrates an important caveat in our approach: the differentiation between host-associated microbes and microbes from other sources may not always be possible, and will be particularly difficult for museum specimens. Though not problematic in the examples we present, the situation is likely more complicated for study of host species with a less well-investigated microbiome, or for symbionts that are very similar to environmental taxa. In these cases, the approach will establish candidates that will then require direct validation.

Finally, as demonstrated previously in other taxa with very similar approaches (Kumar et al., 2013), we show that draft genomes of microbial symbionts can be recovered from Apis short reads. For example, inspecting the non-target components of just a single Apis mellifera intermissa sequencing library produced novel, highly covered, and near complete draft genomes of Lactobacillus kunkeii and a Fructobacillus strain (Figs. 3B–3D). Although the 16S sequence of the Fructobacillus strain best matched F. fructosus, our analysis suggests it belongs to a species so far not represented by genomic sequences in public databases, or even a novel species (Fig. 3C). Conceivably, many additional Lactobacillus variants could be retrieved from the libraries investigated here, potentially providing a more complete picture of the Apis microbiome composition and function. It should be noted that draft genomes reconstructed this way must be regarded as ‘population consensus’ genomes, as opposed to genomes sequenced from cultured bacterial clones. While these genomes cannot be linked to a bacterial clone, they still provide information of metabolic capacities within the Apis microbiome.

In summary, our analyses suggest that ‘contamination’ from honey bee sequencing projects can help to characterize the honey bee microbiome, and inform about its composition, abundance of specific taxa, and through recovery of genomes, about metabolic capacities. Short read repositories thus provide plenty of biological information on known associates, but also on undescribed ones. If, for example, a novel microbial associate of honey bees is discovered, one could potentially use the available short read data (and associated metadata) to learn about its distribution, prevalence, and genome. Or, alternatively, one could confirm the absence of certain taxa with more confidence. Kwong et al. (2017) identified Bombiscardovia sp. and Schmidhempelia sp. as gut microbes of bumble bees (Bombus sp.) that are typically not present honey bees. We also did not find evidence for these taxa in Apis with our screening approach and thus provide further support for their absence based on a large sample size.

Implications of viral reads as common contaminations in Apis RNA sequencing libraries

We found viral reads in about half of the 498 investigated RNA libraries, and DWV type A as the most common viral associate, followed by DWV type B (=Varroa destructor virus-1), and all other viruses we screened for (Fig. 5). DWV has a global distribution (Genersch & Aubert, 2010) and its transmission is facilitated by Varroa mites, which are also globally distributed (Martin et al., 2012; Wilfert et al., 2016). Our findings are thus in line with these previous observations. The detection of DWV despite our deliberate exclusion of RNA sequencing experiments which investigated Varroa or virus treatments can be explained by the fact that, in many cases, DWV infections are asymptomatic (Lanzi et al., 2006). The low proportion of DWV reads found in many samples (as a proxy for viral load) also argues for asymptomatic infections in the bees from which the RNA was extracted (Fig. 5). However, we found a number of RNA sequencing libraries in which DWV read numbers made up a very large proportion of the reads (10–54%, Fig. 5). Very high viral titers are typically found in symptomatic bees (Chen, Higgins & Feldlaufer, 2005; Lanzi et al., 2006; Brettell et al., 2017), and inspection of the metadata associated with RNA libraries showing high proportions of viral reads revealed that many of these were created from heads or brains (out of 58 libraries with >10% viral reads, 25 from heads, 10 from other tissues, 23 undisclosed, Table S4). This is notable because high DWV titres in honey bee heads have been found to correlate with deformed wing phenotypes (Yue, 2005; Zioni, Soroker & Chejanovsky, 2011). This suggests that for these specific experiments, RNA may have been extracted from DWV infected bees despite obvious symptoms (i.e, crippled wings).

Varroa mites and DWV are considered one of the main factors driving the decline in managed honey bee populations (VanEngelsdorp et al., 2009; Conte, Ellis & Ritter, 2010), and molecular data have been used to trace the global spread of these pathogens (Martin et al., 2012; Wilfert et al., 2016). We here showed that almost complete DWV draft genomes can be extracted from short read data, revealing considerable genetic diversity between the strains (Fig. 6). Archived RNA sequencing reads of honey bees are thus a valuable, largely untapped resource for genomic data of viruses. As described above for symbionts in DNA sequencing libraries, these data inform about host association, distribution, and genetic diversity of viruses.

Furthermore, the detection of large amounts of viral reads in honey bee RNA sequencing libraries has other, more direct implications. Differential gene expression analysis was previously used to study e.g., development, learning, caste differences, and other behavioral traits in honey bees (Whitfield, 2003; Liu et al., 2011; Liang et al., 2012; Cameron, Duncan & Dearden, 2013; Naeger & Robinson, 2016). It has further been established that viral infections lead to transcriptomic responses in honey bees, some of which are conserved and universal for all pathogens, while others are specific with respect to the pathogen present (summarized in Doublet et al., 2017). Following from this, the presence of viruses alters transcriptomic profiles of honey bees and therefore incorporated into analyses of differential gene expression, e.g., by adding viral load as one factor in the interpretation. In most experiments analysed here, the presence of viruses was likely not a problem, as viral reads were uniformly low across all samples. However, we did identify some examples in which uneven viral loads very likely had an effect on transcriptomic profiles (Tables 1 and 2, Fig. 7 & Figs. S2–S5). In these examples, the presence of viruses, if unaccounted for, potentially leads to erroneous interpretations of the RNAseq data. Consequently, viruses should be screened for and identified in any RNA sequencing experiment investigating differential gene expression.