Virome compositions indicate that viral spillover is a dead-end between the western honey bee and the common eastern bumblebee

doi:10.21203/rs.3.rs-4802694/v1

Download PDF

Article

Virome compositions indicate that viral spillover is a dead-end between the western honey bee and the common eastern bumblebee

https://doi.org/10.21203/rs.3.rs-4802694/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

The potential of viral spillover from the western honey bee (Apis mellifera) to other insects is well established. New variants should inevitably emerge following a host expansion, yet to our knowledge no study has shown this within this system. To investigate the outcome of viral spillover, we sequenced the meta-transcriptomes of sympatric A. mellifera (n = 389) and common eastern bumblebee Bombus impatiens (n = 117) over three years. Distinct viromes occurred within each bee species throughout the study duration, with honey bee viruses forming a minor fraction of the bumblebee virome. Viruses shared by both bees shared over 98 % nucleotide identity, and no bumblebee-specific strains of honey bee viruses occurred, as expected if spillover led to a true host expansion involving bumblebee-bumblebee transmission. We conclude that the honey bee viruses, namely deformed wing virus, black queen cell virus, and sacbrood virus, were present in the bumblebees due to environmental exposure or dead-end spillover, and not spillover host expansion.

Biological sciences/Microbiology/Virology/Viral transmission

Biological sciences/Ecology/Biodiversity

The interspecific transmission of pathogens occurs when a pathogen is transmitted to a host previously not infected by that pathogen (1–3), one possible outcome of which is spillover, in which the pathogen establishes a new maintenance host. Viral spillover is widely considered a threat to wild insect pollinators, due to the possibility of virus spillover from managed western honey bee (Apis mellifera) colonies. There has been much focus on proposed viral spillover from A. mellifera to bumblebees (genus Bombus), many of which are experiencing population declines with some species facing a high likelihood of extinction (4, 5).

Spillover is a multifaceted phenomenon (6) and a term that is used to refer to a variety of distinct scenarios. We use the terms exposure, dead-end spillover, and spillover as defined in Fig. 1. Briefly, the first and most common scenario is exposure, which is the incidental presence of a virus in or on an organism due to the sharing of a common environment with the true host, and is unrelated to infection (Fig. 1). If a virus exposed to a new host is capable of at least initiating infection, but cannot be sufficiently transmitted between the new hosts, then the result is dead-end spillover (2), in which the virus will inevitably be eliminated from the new host population (Fig. 1).

To establish a new host, a virus must overcome a complex set of barriers, including molecular interactions, as well as host ecology and other biological factors. As a result, exposure and dead-end spillover are vastly more common than spillover, as is seen in other animal systems (7, 8). Even viruses which can overcome important barriers to establishing a new host may still end up at dead-end spillover. For example, even influenza viruses that can complete replication and produce large quantities of virions in a new host can have insufficient virus transmission to establish a new maintenance host (8).

We reserve the term spillover for the scenario in which the virus overcomes these barriers and is able to maintain itself in the new host population through transmission between the new hosts (Fig. 1). As it can maintain itself, the presence of the virus in the new host is not dependent on contact with viral reservoirs in the origin host.

A key feature of spillover is heightened evolution in the pathogen (9) and the rapid emergence of variants that are specific to the new host (Fig. 1), whereas in the exposure and dead-end spillover scenarios, the virus sequence identity remains essentially identical to that of the origin host due to unsuccessful propagation of the virus.

Although there is a large body of research proposing spillover between honey bees and various other arthropods (10), spillover-related scenarios have been poorly distinguished. The scenario of concern discussed in the literature is spillover, as the establishment of a new maintenance host poses the greatest threat to host populations (Fig. 1). However, the majority of spillover evidence is based on the assessment of virus presence, which has limited capability to distinguish between exposure, dead-end spillover, and spillover. As a result, it is likely that cases of exposure or dead-end spillover have been widely mischaracterized as spillover. Correctly identifying the spillover-related scenarios occurring in bees is crucial, as they pose different degrees of risk to pollinators, with the risk increasing from exposure to dead-end spillover, to spillover.

A relatively narrow range of viruses have been studied in the context of bee virus spillover, mainly the A. mellifera viruses deformed wing virus (DWV, Iflavirus aladeformis), black queen cell virus (BQCV, Triatovirus nigereginacellulae) and sacbrood virus (SBV, Iflavirus sacbroodi) (10). A key example of proposed viral spillover is the prominent focus on the spillover of DWV from A. mellifera to B. impatiens, which currently remains unproven. Wild Bombus species are commonly RT-PCR-positive for DWV, and DWV is able to complete its replication cycle in some members of this genus, at least in artificial experiment systems (11), and DWV replication can be detected in wild Bombus species (12). The main hypothesised route of spillover for DWV, and other viruses, is through oral ingestion of viruses deposited on shared flowers, as supported by the RT-PCR detection of DWV on flowers visited by A. mellifera (12). There is, however, also a case against DWV spillover. Firstly, there is little evidence of any negative outcomes of proposed viral spillover in wild pollinators (10, 13). Other studies suggest that successful DWV transmission between bumblebees is implausible in natural contexts (11, 14), as DWV transmission is highly dependent on the Varroa mite which does not parasitise Bombus species (15) and DWV infection is highly ineffective via the oral-faecal route (12, 16, 17). Spillover has been proposed for a variety of other honey bee viruses, albeit with less available evidence than for DWV (10).

Considering insights from other animal virus spillover systems (7, 8), the evidence for viral spillover from A. mellifera to other insects, appears to be highly uncertain. A meta-transcriptomic approach can address some of this uncertainty, specifically whether viral variants adapted to a new host have emerged following a spillover event. By comparing the genomes of the total viral community, variants that have emerged within a new host can be identified, which may indicate whether the presence of A. mellifera viruses in other insects is a result of exposure, dead-end spillover, or spillover. To determine the state of viral spillover between A. mellifera and B. impatiens, we surveyed the meta-transcriptomes of sympatric A. mellifera and B. impatiens from the field study program in Minnesota, USA. Using this approach we tested for evidence of spillover for all viruses identified in both bees by attempting to identify B. impatiens-specific variants of A. mellifera viruses using a genomic approach.

Sequencing libraries were generated for a total of 408 Apis mellifera and 131 Bombus impatiens samples. For A. mellifera and B. impatiens, respectively, 19 and 14 sample libraries were excluded from the analyses because they generated no contigs or generated only contigs below 1000 bp, leaving a total of 389 and 117 samples with sequencing data. For A. mellifera and B. impatiens, respectively, the average number of reads (that passed the Q-score threshold) per library was 80,563 and 101,004, the average read length was 550 bp and 487 bp, and the percentage of reads that mapped to viral contigs was 9.7% and 4.6%.

Across all samples, viral communities were distinct between A. mellifera and B. impatiens (PERMANOVA Df = 1, R² = 0.007, p = 0.001 ***). The only experimental variables which significantly impacted the community compositional similarity between A. mellifera and B. impatiens was the collection year (PERMANOVA Df = 2, R² = 0.01, p = 0.001 ***), as it was the only variable with a significant interaction with the host genus (Table S1; Fig. S1).

Viral communities in B. impatiens were overall more homogenous between the sampling years than in A. mellifera. The key viruses that distinguished the two bee species were (assessed by PCA, Fig. 2, and Pearson correlation coefficients, Fig. S2, and community composition, Fig. 3), in A. mellifera: Iflavirus aladeformis, apis rhabdovirus 1, Cripavirus ropadi, hubei partiti − like virus 34, Cripavirus mortiferum, lake sinai virus 3, lake sinai virus 6, unclassified sinaivirus, Lake Sinai virus 2, Aparavirus apisacutum, and Iflavirus sacbroodi (Fig. S2). B. impatiens was associated with vespa velutina associated permutotetra − like virus 1, hymenopteran phasma − related virus OKIAV231, hubei picorna − like virus 27, cyclosorus interruptus picorna − like virus, ganda orthophasmavirus, unclassified phasmaviridae, mayfield virus 1, hymenopteran phasma − related virus OKIAV234, elf loch virus, bombus − associated virus reo1, agassiz rock virus, bactrocera tryoni iflavirus 1, bombus − associated virus Vir4, andrena haemorrhoa nege − like virus, and allermuir hill virus 1 (Fig. S2). The remaining viruses had a neutral correlation (r = 0.0 or -0.0) with either host genus (Fig. S2) and low relative abundances of reads (Fig. 3, Fig. S5): Aparavirus israelense, Lake Sinai virus 1, vespa velutina partiti − like virus 2, lake sinai virus, Triatovirus nigereginacellulae, chronic bee paralysis virus, and hymenopteran phasma − related virus OKIAV233.

Across collection years, variation in viral community compositions led to minor differences in similarity between A. mellifera and B. impatiens. Communities in A. mellifera and B. impatiens showed slightly closer clustering in 2022 and 2023 than in 2021 (Fig. S3). Overall, the viruses which contributed most to the A. mellifera and B. impatiens communities becoming more similar were Lake Sinai virus 1, Lake Sinai virus 2, unclassified sinaivirus, lake sinai virus, ganda orthophasmavirus, Cripavirus mortiferum, unclassified phasmaviridae, and hymenopteran phasma − related virus OKIAV233 (Fig. S3). Viral community shifts were observed in the transplanted colonies, but overall the communities remained similar to flower-collected bees in the same year (Fig. S6A). One of the transplanted colonies showed a shift to have a greater proportion of B. impatiens-associated viruses, primarily mayfield virus 1 (Fig. S6B).

For all viruses identified in this study, the A. mellifera virus sequences in B. impatiens, and vice versa, did not form distinct clusters or clades indicating a strain resulting from spillover (Fig. 4, Fig. S7). Instead, they are likely the result of exposure or dead-end spillover (as defined in Fig. 1). This is consistent with other reports of RT-PCR amplicons in Bombus, which follow those in A. mellifera - for example, the B genotype of DWV that has come to be dominant in A. mellifera is also found in Bombus (18) and it is unclear if host species has driven the emergence of A. mellifera virus strains specific to other arthropod hosts (19). Our observations of highly similar A. mellifera viruses in B. impatiens is consistent with other spillover studies which report apparent spillover only when Bombus and other arthropods are sympatric with A. mellifera (10, 19–21), which is indicative of exposure or dead-end spillover, since these viruses do not seem to be maintained in Bombus in the absence of A. mellifera.

Overall, the occurrence of an A. mellifera virus in B. impatiens, and vice versa, coincided with shorter contig length (Fig. S4, Table S2). Though comparable normalised read counts assigned to certain viruses being shared between A. mellifera and B. impatiens (Fig. S5), the discrepancy in contig lengths between A. mellifera and B. impatiens suggested that the contigs resulting from apparent spillover generated more fragmented genomes. Exceptions to this were lake sinai virus, unclassified phasmaviridae, and Triatovirus nigereginacellulae, for which A. mellifera and B. impatiens did not have significantly different contig lengths (Table S2). Across the three study years, the B. impatiens community maintained a consistent composition (Fig. 3), with the A. mellifera viruses in it remaining as minor portions of it, and never expanded to dominate the B. impatiens community. Of the A. mellifera pathogens identified, it was SBV that formed the largest portion of the B. impatiens virome (Fig. 3) and full SBV genomes were frequently recovered from B. impatiens (Fig. S4). Nevertheless, the sequence identity of these contigs remained at over 98% nucleotide identity between A. mellifera and B. impatiens (Fig. 4), suggesting that these high SBV levels in B. impatiens are simply the result of high exposure of SBV from A. mellifera, and not of spillover to B. impatiens (Fig. 1). Together, the community composition and shorter contigs support the exposure or dead-end spillover scenario.

No B. impatiens-specific clades of A. mellifera viruses were observed in RNA dependent RNA polymerase phylogeny (Fig. 4B, Fig. S7, Fig. S8). Clustering based on the average nucleotide identity of whole contigs (Fig. 4A) and percent identity of RNA dependent RNA polymerase gene sequences (Fig. S9) often resulted in clusters composed of sequences derived from both A. mellifera and B. impatiens, with no clusters showing emerging subclusters or clusters specific to either host species. Reference sequences based networks for sequences over 200 bp length from Apis mellifera, Apis cerana, and Bombus spp. only (Fig. S10, Fig. S11) showed various strain clusters within Apis including A and B for DWV, and SBV strains in A. cerana, whilst networks based on sequences over 4000 bp from a variety of insect taxa showed similar strain distinctions (Fig. S12, Fig. S13).

In addition, the networks generated from public reference sequences showed a similar pattern - at 98% identity, DWV, SBV, and BQCV separate into distinct genotype clusters within A. mellifera, but form no Bombus-specific clusters (Fig. S10, Fig. S11). For comparison, SBV separates into clear A. mellifera and Apis cerana strains (22), which can be seen in our reference sequence networks, demonstrating that this network approach can separate virus strains adapted to different hosts (Fig. S10, Fig. S12, Fig. S13). The lack of evolutionary divergence of A. mellifera viruses identified in insects other than A. mellifera, reinforces the argument that the widespread detection of A. mellifera viruses in other insects may often be the result of exposure or dead-end spillover. There is however, a need to expand upon the viral genomic datasets available for wild pollinators, as a clear strain cluster might only be observed after systematic surveillance has been conducted.

To date, the next generation sequencing-based studies for the surveillance of pathogens in bees are still limited. Here we have shown that bumblebees have a distinct viral community, with the presence of A. mellifera viruses attributable to exposure or dead-end spillover from honey bees to bumblebees. The transplantation of B. impatiens colonies from indoors to sites with A. mellifera did not increase the proportion of A. mellifera viruses in B. impatiens (Fig. S6). Instead the variation in exposure to A. mellifera only seemed to produce a shift in the putative B. impatiens viruses already present in the colonies. This concurs with previous studies that has shown that Bombus colonies transplanted to sites free of A. mellifera eliminate DWV and BQCV and reestablish a Bombus-specific virome (23), and that wild bee species including Bombus possess distinct viral communities, even when they are closely related and share a common environment (24–28). Although we attempted to sample B. impatiens at varying levels of exposure to A. mellifera by sampling at multiple distances from apiaries, the lack of effect of distance on virome composition may have been a result of homogenous exposure to A. mellifera regardless of the distances sampled due to the presence of other apiaries in the area.

The factors that influenced the similarity between honey bee and bumblebee viromes were unclear, as the factors that significantly interacted with host genus and year, did not coincide with any meaningful convergence between the viromes of A. mellifera and B. impatiens. None of the viruses involved in these shifts were key honey bee pathogens, and many were minor members of the viral community (ganda orthophasmavirus, Cripavirus mortiferum, unclassified phasmaviridae). Furthermore, it is not clear if poorly classified viruses such as these originated from honey bees, bumblebees, or other sources, such as other insects. The presence of various lake sinai related viruses is interesting, as one classification (lake sinai virus) did increase in abundance in B. impatiens in 2022, though no B. impatiens-specific cluster was observed (Fig. 4). Since the lake sinai viruses are a complex of diverse viral strains (29) with phylogeny and host ranges that are still being resolved (30), more systematic surveys would be needed to distinguish potential spillover from natural host ranges in Bombus. Similarly, more expansive genomic studies would also be needed for many of the poorly characterised viruses (such as those dominating the B. impatiens virome, Fig. 2 and Fig. 3) to determine the natural host ranges of these viruses.

Evidently, wild bumblebees are exposed to a wide range of honey bee viruses, but insights from other systems indicate that only a small minority of these viruses and exposure events could lead to spillover (1–3, 6, 7) and therefore pose a significant threat to wild pollinators. It is likely that many cases of viral exposure between insects have occurred over long evolutionary timescales without spillover, due to biological barriers that maintain distinct viral host ranges, and these dynamics must be distinguished from actual spillover threats. Our metagenome survey demonstrates an approach similar to that used in other spillover systems (9) by which spillover can be distinguished by the background noise of widespread viral exposure and dead-end spillover. Even in the case of DWV, where the virus can replicate in a new host that it frequently encounters in the wild, no spillover may ever occur if the virus cannot surmount the biological and ecological barriers to it establishing a new maintenance host. There is a dynamic and complex virus system in wild insect pollinators that has been little explored, and a need to integrate genomic information into current efforts to study insect virus spillover. Further genomic studies may yet find evidence of viral spillover if they identify honey bee virus strains that are specific to bumblebees and other pollinators.

Although this study shows a lack of spillover between A. mellifera and B. impatiens, protection of bumblebees of conservation concern from potential pathogen spillover from honey bees is still a valid conservation action, particularly when protecting species at risk of extinction. Low genetic diversity, which has been demonstrated in bumblebee species of conservation concern, could increase vulnerability to pathogen spillover (31–33). Due to the lethal nature of sampling for this study, bees of conservation concern were not examined. While this examination of B. impatiens is an important step in understanding the interactions between the viromes of honey bees and bumblebees, B. impatiens appears to be expanding in range and is not expected to be representative of all bumblebee species, particularly those of conservation concern (34, 35). This study adds to a small but expanding view of insect viromes and is an attempt to distinguish the largely overlooked outcomes of cross-species viral exposure between insects. Unsurprisingly, the common eastern bumblebee has a distinct viral community from honey bees, and the question must be raised as to what role virus-virus interactions play in these dynamics. With around 250 species of bumblebee globally, there is also much to learn of differences in viral communities among bumblebee species, particularly those of conservation concern. Establishing what the common state and composition of insect pollinator viromes is essential, as it provides a baseline with which to compare proposed cases of spillover.

Field sampling

Between the years 2021 and 2023 in Twin Cities, Minnesota USA, A. mellifera and B. impatiens were collected from three apiary sites between the months May and November.

At each site (Golf, Crosby, Vets) A. mellifera was sampled from colonies, repeated in replicates of 8 colonies per site every 3 weeks, collecting at least 30 worker bees per replicate. Bees, both A. mellifera and B. impatiens, were also collected from flowers during 3 collection periods (early July, early August, and late August) when B. impatiens colonies are at peak colony size. Bees were collected within 100 m of the apiary site, between 500 and 1000 m of the apiary site, and between 1500 and 2000 m of the apiary site. These distances were chosen to form a gradation of flower sharing, with fewer honey bees from the apiary site being present at farther distances, though there was no control for other potential apiary sites in the area that were not part of this study. When possible, four bees were collected from each distance area for each collecting period. Outside of the apiary sites, bees were collected from areas with flowering plants in roadsides or public parks. Lack of bee forage and low bee abundance in some areas resulted in some collection periods with fewer than 4 individuals per species and distance category. In 2021, bees were collected from all flowers. In 2022 and 2023, floral visitation was assessed before each collection period by observing bee visitation to flowers in the area for 20 minutes. Bee collection focused on flowers on which both honey bees and bumblebees were observed foraging. During times of low bee abundance, bumblebees were collected from flowers without observation of shared foraging by honeybees. Replicate samples were collected by year, collection period, distance from apiary, and bee species.

In 2023, two B. impatiens colonies were also sampled before and after transplantation near the Golf apiary site. One B. impatiens colony was initiated from a field caught queen in the spring using standard rearing techniques (36). The other colony was collected from the field and kept inside for 1 week then transplanted to site Golf. Colonies were sampled before transplantation (date 0), then repeatedly at multiple dates with one week between samplings (between July and August) from date 0 to 5–6 dates after transplantation, collecting 4 bees per sampling on average.

Dissociation and RNA extraction

Bees without ectoparasites and signs of disease were dissociated using a gentleMACS dissociator in 10 mL of molecular grade water using setting RNA 02.01. The resulting homogenate was centrifuged at 21,100 g for ten minutes at room temperature. Next, 150 µL of the supernatant fluid was retained for bead-based RNA extraction using a NucleoMag Virus RNA/DNA Extraction kit (Takara, Shiga, Japan), according to the manufacturer’s instructions and using a KingFisher Flex automated extraction machine (ThermoFisher, Walton, MA). The elution volume was 50 µL of molecular grade water.

cDNA synthesis

The template used for cDNA synthesis was 4 µL of extracted RNA per sample, which was annealed with 1 µL of 10 µM custom N6 template switching oligonucleotide (TSO; GCAGTGGTATCAACGCAGAGTACNNNNNNN) and 1 µL of 10 mM dNTPs. The annealing mix was incubated at 70°C for 5 minutes. Next the total volume of the annealed mix (total volume 6 µL per sample) was used as the template for first strand synthesis, to a total volume of 10 µL per sample. For each sample, 2.5 µL of template switching buffer (New England Biolabs, Ipswich, US), 1 µL of template switching enzyme (New England Biolabs, Ipswich, US), and 0.5 µL of a 75 µM TSO (AAGCAGTGGTATCAACGCAGAGTACrGrGrG) was used. This first strand synthesis mix was then cycled for 90 minutes at 42℃, and finally for five minutes at 85℃. Finally, 25 µL of PrimeSTAR GXL Premix (Takara, Shiga, Japan), two µL of 10 µM TSO primer (AAGCAGTGGTATCAACGCAGAGTAC), and 13 µL of molecular grade water were added to each sample to a total volume of 50 µL. The second strand synthesis was incubated at 94°C for one minute, then for 30 cycles at 98°C for one minute, 30 cycles at 60℃ for 15 seconds, and 30 cycles at 68℃ for 2.5 minutes.

Clean up of the cDNA was performed by adding 50 µL of CleanNGS DNA & RNA Clean-Up Magnetic Beads to 50 µL of the PCR-amplified cDNA, then incubated at room temperature for 10 minutes. After incubation, samples were centrifuged at 21,000 g force for ten seconds and then placed on a magnetic tube stand. The clear portion of the liquid was removed and discarded after the beads had collected on the side of the tubes, and 200 µL of 80% ethanol was added, left for 30 seconds, and then removed, and then this ethanol was repeated once more. The tubes were centrifuged for 10 seconds at 1,000 rpm or until excess ethanol collected at the bottom and was then removed. The tubes were then left open to dry at room temperature for 1–3 minutes, then resuspended in 20 µL of molecular grade water, incubated for 10 mins at room temperature before placing back on the magnetic stand and retaining the cDNA eluate.

Sequencing

Sequencing was performed using the Oxford Nanopore Technologies (ONT) sequencing technology, with an input of 12 uL of cDNA per sample, R10 flow cells and SQK-NBD114-96 kit according to the manufacturers' instructions (using the kit SQK-RBK110-96 and its option of pooling the barcoded samples together). The sequencing run durations were 24 hrs and used the MinKNOW (v21.02.5) software and ran on a GridION (Oxford Nanopore Technologies, Oxford, UK). Basecalling was performed using Guppy (v4.2.0) with the high-accuracy model and with all default settings applied. Sequencing reads were filtered to a minimum length (≥ 200 bp) and Q-score (≥ 9) and demultiplexed by MinKNOW v24.02.16.

Contig assembly

Adapter sequences were removed using PoreChop v0.2.4 was used to remove the nanopore barcode adapter sequences. To generate metagenome-assembled contigs, the quality filtered reads were de novo assembled by Canu (v2.2) (37) (with parameters maxInputCoverage = 10000 corOutCoverage = all corMinCoverage = 0 minReadLength = 200 minOverlapLength = 25 corMhapSensitivity = high genomeSize = 5000). All contigs below 1000 base pair length were excluded from all analyses. To compare the lengths of viral contigs between A. mellifera and B. impatiens, contig lengths were visualised (Fig. S4) using ggbeeswarm (v0.7.2).

Viral Sequence Identification

To identify and bin the viral contigs, the outputs of several tools were combined for assessment. To taxonomically classify the contigs, Kaiju (v1.9.2) was used (default parameters against both the Reference Viral Database (38) (RVDB; v27) and NCBI non-redundant (nr) datasets built into Kaiju). To assist with contig binning, MMseqs2 (v15.6f452; easy-cluster --min-seq-id 0.9 -c 0.8 --cov-mode 5) was used to assign viral contigs to clusters. All contigs were queried against the Virus Orthologous Groups Database (VOGDB, release 217) and RVDB using diamond BLASTx (v2.1.8.162; --range-culling --max-target-seqs 1 -F 15 --evalue 1e-20 --sensitive). For all contigs, checkV (v1.0.1; end_to_end) was run, and the genome quality assessment was also used to filter and bin the contigs. The combined output of these tools was used to assess the contigs and place them into bins. Binning the contigs down to species level (or as low a taxonomic level as possible) was performed according to the Kaiju classifications, and the final bins were named after these classifications. These classifications represent the nearest best hit in the databases, and do not represent accurate classification at lower taxonomic levels for poorly represented viral classifications. The contig binning assessment was performed as follows: contigs were binned as Non-viral if they had the CheckV warning “no viral genes”, and no viral hit against Kaiju, and no BLASTx hit against RVDB, and no BLASTx hit against VOGDB. Otherwise, contigs were considered viral. Some contigs were removed from the viral bin if they had a combination of CheckV “no viral genes” or CheckV “Low-quality” and no BLASTx hits versus VOGDB and RVDB. Furthermore, contigs with certain VOG hits associated with cellular organisms, such as gag pol retrotransposons were removed from viral bins. Due to taxonomic discrepancies between NCBI and the International Committee on Taxonomy of Viruses (ICTV), various assignments by Kaiju were manually corrected to align with those of ICTV.

Phylogeny

Gene predictions were generated using Prodigal (v2.6.3, parameters: -g 11 -p meta), followed by Pfam annotation of genes using Interproscan (v5.23-62.0, protein query). Protein sequences were grouped by bin and Pfam, and sequences corresponding to RNA dependent RNA polymerase (RdRp) were aligned using MAFFT (v7.525) (parameters: --auto). Before alignment, proteins less than half the average length of the viral cluster were removed. Alignments were trimmed using TrimAl (v1.4.rev15, parameters: -automated1) and phylogeny performed using IQTREE (v2.2.6, parameters: -m MFP) with 1000 bootstraps. Trees were visualised using ggtree (v3.6.2) in R for iflaviruses (Fig. 4B) and other taxa (Fig. S7: apis rhabdovirus, Dicistroviridae, Iflavirus aladeformis, bactrocera tryoni iflaviruses, Iflavirus sacbroodi, nege-like virus; Fig. S8: partiti-like virus, unclassified phasmaviridae, picorna-like, Sinaivirus, virga-like).

Metagenomic analyses

Prior to mapping, all assembled contigs were de-replicated using CD-HIT (v4.8.1, parameters: -c 0.9 -s 0.8). Reads were mapped using Minimap2 (v2.26, parameters: -ax map-ont --secondary = no), mapping the reads of each sample independently against the concatenated contigs from all samples. Secondary read mappings were removed from the mapping results using SAMtools (v1.18; -h -F 2308). To normalise the read counts by library size, the R package metagenomeSeq (v1.40.0) was used to perform cumulative sum scaling (CSS), without transformation of the CSS counts. To examine the yield of reads for viral taxa, the CSS counts of reads mapped against the classified contigs were visualised as a heatmap (Fig. S5). The CSS counts of mapped reads were converted to a phyloseq object in R (v4.2.2) using the package phyloseq (v1.42.0), with the incorporation of the Kaiju classifications as taxa, and the sample metadata. Using this phyloseq object as a basis, the R package microViz (v0.11.0) was used to perform principal component analysis (PCA), redundancy analysis (RDA), correlation heatmaps, and composition plots. All plots were visualised using R’s ggplot2 (v3.5.0) or base R. The PCA and RDA plots were generated using different variables to identify which viral taxa or experimental variables were contributing to the composition of viral communities. They were used to identify how viral community similarity was affected by viral taxa across all samples (Fig. 2 - shown split by collection year in Fig. S3), for the transplanted B. impatiens colonies (Fig. S6A), and for the apiary, distance from apiary, and flower genus from which the samples were collected (Fig. S1). To further examine the association of specific viral taxa with host genus and collection year, correlation heatmaps were also generated, with the Pearson correlation coefficients indicating the strength of association with either bee species (Fig. S2). The composition of viral communities were visualised as barplots for all samples (Fig. 3) and the transplanted B. impatiens colonies (Fig. S6B).

Network analyses

To generate percent similarity scores, fastANI (v1.33, parameters: --fragLen 200) was run on the whole contigs, all versus all. For RdRp, the protein sequences were queried against themselves using diamond BLASTp. The similarity scores were imported into Cytoscape (v3.10.2) and visualised with the organic layout, and hiding all edges below 98% identity. Networks were generated for this study’s sequences with some reference sequences for average nucleotide identity (Fig. 4A) and protein percent identity (Fig. S9), and also for reference virus sequences downloaded from NCBI, for nucleotide sequences over 200 bp length and recovered from Bombus spp., Apis mellifera, and Apis cerana (Fig. S10) and RdRp protein sequences predicted from these sequences (Fig. S11). Nucleotide and RdRp networks were also generated for reference sequences over 4000 bp length and derived from various insect taxa (Fig. S12, Fig. S13). Reference sequences were accessed in June 2024 from NCBI Virus, selecting nucleotide sequences under the taxids for Iflavirus sacbroodi: 3048309, 89463, 886669, Iflavirus aladeformis: 198112, 3047792, 2951794, 2498590, and Triatovirus nigereginacellulae: 3047684.

Statistical Analysis

Permanova was performed using adonis2 (method="bray") in R, using an input matrix of the CSS read counts for each contig classified by viral taxa. For permanova, interactions between the variables host genus, collection year, distance from apiary, apiary, and flower genus were included (Table S1). These variables were included to assess whether viral community composition varied across the years, between sites (apiaries) between A. mellifera and B. impatiens. Distance was considered because the viral community was expected to be impacted by decreasing A. mellifera presence with increased distance from the apiaries. Flower genus was expected to alter viral community composition due to differences in the physical properties of the flowers and the levels of A. mellifera and B. impatiens visiting different flower types. In addition, other flower properties were included in the permanova, including flower depth (shallow = < 3 mm, medium = 4–6 mm, deep = > 8 mm), flower shape (classified as composite, cup, pea, tube), and the ratio of A. mellifera and B. impatiens observed visiting that flower species during the study. Contig lengths were tested for differences in length using a Wilcox test (wilcox.test) in base R (Table S2).

Funding: Funding provided by the Environment and Natural Resources Trust Fund as recommended by the Legislative-Citizen Commission on Minnesota Resources. The UMN Bee Squad Program received private donations from individuals and the Minnesota Saint Paul Airport (MSP) to support beekeeping activities taken at the Crosby, Golf, and Vet sites.

Author contributions: Conceptualization: DCS, MS. Methodology: DCS, MS, EE, JH, JW, RZ, BM, RM, AB, MN. Software: DAM, PHB. Validation: DAM. Formal analysis: DAM, DCS Investigation: DAM, EE, JH, JW, RZ, BM, RM. Resources: DCS, MS. Data Curation: DAM, PHB. Writing - Original Draft: DAM, DCS. Writing - Review & Editing: All. Visualisation: DAM. Supervision: DCS. Project administration: DCS. Funding acquisition: DCS.

Competing interests: Authors declare that they have no competing interests.

Data and materials availability: All data are available in the main text or the supplementary materials. Sequencing reads were submitted to NCBI SRA under BioProject PRJNA1119970. All contig assemblies are available in FigShare at https://figshare.com/s/a632154b4c9284e20b4a. All scripts used in analyses can be found at a GitHub repository at https://github.com/dmckeow/bioinf, specifically the following scripts in the bin/ folder: bioinf-assembly-canu.sh, bioinf-binning.sh, bioinf-phylogeny.py, spillover-CovDepth.R, spillover-annotations.R, spillover-annotations.sh, spillover-main.R, spillover-main.sh, spillover-network.R, spillover-network.sh, spillover-phylogeny.R.

Daszak, P., Cunningham, A.A., Hyatt, A.D.: Emerging Infectious Diseases of Wildlife– Threats to Biodiversity and Human Health. Science. 287, 443–449 (2000)
Ellwanger, J.H., Chies, J.A.B.: Zoonotic spillover: Understanding basic aspects for better prevention. Genet. Mol. Biol. 44, e20200355
Power, A.G., Mitchell, C.E.: Pathogen Spillover in Disease Epidemics. Am. Nat. 164, S79–S89 (2004)
Colla, S.R., Gadallah, F., Richardson, L., Wagner, D., Gall, L.: Assessing declines of North American bumble bees (Bombus spp.) using museum specimens. Biodivers. Conserv. 21, 3585–3595 (2012)
Goulson, D., Lye, G.C., Darvill, B.: Decline and Conservation of Bumble Bees. Annu. Rev. Entomol. 53, 191–208 (2008)
Fenton, A., Pedersen, A.B.: Community Epidemiology Framework for Classifying Disease Threats. Emerg. Infect. Dis. 11, 1815–1821 (2005)
Letko, M., Seifert, S.N., Olival, K.J., Plowright, R.K., Munster, V.J.: Bat-borne virus diversity, spillover and emergence. Nat. Rev. Microbiol. 18, 461–471 (2020)
Wasik, B.R., de Wit, E., Munster, V., Lloyd-Smith, J.O., Martinez-Sobrido, L., Parrish, C.R.: Onward transmission of viruses: how do viruses emerge to cause epidemics after spillover? Philos. Trans. R Soc. B Biol. Sci. 374, 20190017 (2019)
Tan, C.C.S., van Dorp, L., Balloux, F.: The evolutionary drivers and correlates of viral host jumps. Nat. Ecol. Evol. 8, 960–971 (2024)
Nanetti, A., Bortolotti, L., Cilia, G.: Pathogens Spillover from Honey Bees to Other Arthropods. Pathogens. 10, 1044 (2021)
Tehel, A., Streicher, T., Tragust, S., Paxton, R.J.: Experimental cross species transmission of a major viral pathogen in bees is predominantly from honeybees to bumblebees. Proc. R. Soc. B Biol. Sci. 289, 20212255 (2022)
Yañez, O., Piot, N., Dalmon, A., de Miranda, J.R., Chantawannakul, P., Panziera, D., Amiri, E., Smagghe, G., Schroeder, D., Chejanovsky, N.: Bee Viruses: Routes of Infection in Hymenoptera. Front. Microbiol. 11, 943 (2020)
Cameron, S.A., Sadd, B.M.: Global Trends in Bumble Bee Health. Annu. Rev. Entomol. 65, 209–232 (2020)
Tehel, A., Streicher, T., Tragust, S., Paxton, R.J.: Experimental infection of bumblebees with honeybee-associated viruses: no direct fitness costs but potential future threats to novel wild bee hosts. R Soc. Open. Sci. 7, 200480 (2020)
Ryabov, E.V., Wood, G.R., Fannon, J.M., Moore, J.D., Bull, J.C., Chandler, D., Mead, A., Burroughs, N., Evans, D.J.: A Virulent Strain of Deformed Wing Virus (DWV) of Honeybees (Apis mellifera) Prevails after Varroa destructor-Mediated, or In Vitro, Transmission. PLOS Pathog. 10, e1004230 (2014)
Gusachenko, O.N., Woodford, L., Balbirnie-Cumming, K., Ryabov, E.V., Evans, D.J.: Evidence for and against deformed wing virus spillover from honey bees to bumble bees: a reverse genetic analysis. Sci. Rep. 10, 16847 (2020)
Streicher, T., Tehel, A., Tragust, S., Paxton, R.J.: Experimental viral spillover can harm Bombus terrestris workers under field conditions. Ecol. Entomol. 48, 81–89 (2023)
Manley, R., Temperton, B., Doyle, T., Gates, D., Hedges, S., Boots, M., Wilfert, L.: Knock-on community impacts of a novel vector: spillover of emerging DWV-B from Varroa-infested honeybees to wild bumblebees. Ecol. Lett. 22, 1306–1315 (2019)
Dobelmann, J., Felden, A., Lester, P.J.: Genetic Strain Diversity of Multi-Host RNA Viruses that Infect a Wide Range of Pollinators and Associates is Shaped by Geographic Origins. Viruses. 12, 358 (2020)
Alger, S.A., Burnham, P.A., Boncristiani, H.F., Brody, A.K.: RNA virus spillover from managed honeybees (Apis mellifera) to wild bumblebees (Bombus spp). PLOS ONE. 14, e0217822 (2019)
Fürst, M.A., McMahon, D.P., Osborne, J.L., Paxton, R.J., Brown, M.J.F.: Disease associations between honeybees and bumblebees as a threat to wild pollinators. Nature. 506, 364–366 (2014)
Chang, J.-C., Chang, Z.-T., Ko, C.-Y., Chen, Y.-W., Nai, Y.-S.: Genomic Sequencing and Comparison of Sacbrood Viruses from Apis cerana and Apis mellifera in Taiwan. Pathogens. 10, 14 (2020)
Bosco, L., Yañez, O., Schauer, A., Maurer, C., Cushman, S., Arlettaz, R., Jacot, A., Seuberlich, T., Neumann, P., Schläppi, D.: Honey bee associated viruses are unlikely to impact bumble bee colonies while habitat heterogeneity supports their resiliencepreprint, Preprints, ; (2023). https://doi.org/10.22541/au.169384326.69650793/v1
Brettell, L.E., Schroeder, D.C., Martin, S.J.: RNAseq of Deformed Wing Virus and Other Honey Bee-Associated Viruses in Eight Insect Taxa with or without Varroa Infestation. Viruses. 12, 1229 (2020)
Pascall, D.J., Tinsley, M.C., Obbard, D.J., Wilfert, L.: Host evolutionary history predicts virus prevalence across bumblebee species. bioRxiv [Preprint]. (2019). https://doi.org/10.1101/498717
Robinson, C.R.P., Dolezal, A.G., Newton, I.L.G.: Host species and geography impact bee-associated RNA virus communities with evidence for isolation by distance in viral populations. ISME Commun. 4, ycad003 (2024)
Schoonvaere, K., Smet, L.D., Smagghe, G., Vierstraete, A., Braeckman, B.P., de Graaf, D.C.: Unbiased RNA Shotgun Metagenomics in Social and Solitary Wild Bees Detects Associations with Eukaryote Parasites and New Viruses. PLOS ONE. 11, e0168456 (2016)
Schoonvaere, K., Smagghe, G., Francis, F., de Graaf, D.C.: Study of the Metatranscriptome of Eight Social and Solitary Wild Bee Species Reveals Novel Viruses and Bee Parasites. Front. Microbiol. 9 (2018)
Hou, C., Liang, H., Chen, C., Zhao, H., Zhao, P., Deng, S., Li, B., Yang, D., Yang, S., Wilfert, L.: Lake Sinai virus is a diverse, globally distributed but not emerging multi-strain honeybee virus. Mol. Ecol. 32, 3859–3871 (2023)
Bigot, D., Dalmon, A., Roy, B., Hou, C., Germain, M., Romary, M., Deng, S., Diao, Q., Weinert, L.A., Cook, J.M., Herniou, E.A., Gayral, P.: The discovery of Halictivirus resolves the Sinaivirus phylogeny. J. Gen. Virol. 98, 2864–2875 (2017)
Kent, C.F., Dey, A., Patel, H., Tsvetkov, N., Tiwari, T., MacPhail, V.J., Gobeil, Y., Harpur, B.A., Gurtowski, J., Schatz, M.C., Colla, S.R., Zayed, A.: Conservation Genomics of the Declining North American Bumblebee Bombus terricola Reveals Inbreeding and Selection on Immune Genes. Front. Genet. 9 (2018)
Mola, J.M., Pearse, I.S., Boone, M.L., Evans, E., Hepner, M.J., Jean, R.P., Kochanski, J.M., Nordmeyer, C., Runquist, E., Smith, T.A., Strange, J.P., Watson, J., Koch, J.B.U.: Range-wide genetic analysis of an endangered bumble bee (Bombus affinis, Hymenoptera: Apidae) reveals population structure, isolation by distance, and low colony abundance. J. Insect Sci. 24, 19 (2024)
Lively, C.M.: The Effect of Host Genetic Diversity on Disease Spread. Am. Nat. 175, E149–E152 (2010)
Evans, E., Ascher, J.S., Cariveau, D.P., Spivak, M.: A century of change for bees and their floral associations. Elements (2022)
Looney, C., Strange, J.P., Freeman, M., Jennings, D.: The expanding Pacific Northwest range of Bombus impatiens Cresson and its establishment in Washington State. Biol. Invasions. 21, 1879–1885 (2019)
Evans, E., Strange, J., Sadd, B., Tripodi, A., Figueroa, L., Adams, L., Colla, S., Duennes, M., Lehmann, D., Moylett, H., Richardson, L., Smith, J.W., Smith, T., Spevak, E., Inouye, D.W.: Parasites, parasitoids, and hive products that are potentially deleterious to wild and commercially raised bumble bees (Bombus spp.) in North America. J. Pollinat. Ecol. 33, 37–53 (2023)
Koren, S., Walenz, B.P., Berlin, K., Miller, J.R., Bergman, N.H., Phillippy, A.M.: Canu: scalable and accurate long-read assembly via adaptive k-mer weighting and repeat separation. Genome Res. 27, 722–736 (2017)
Goodacre, N., Aljanahi, A., Nandakumar, S., Mikailov, M., Khan, A.S.: A Reference Viral Database (RVDB) To Enhance Bioinformatics Analysis of High-Throughput Sequencing for Novel Virus Detection. mSphere 3, 1–18 (2018)

There is NO Competing Interest.

Supplementarymaterial.docx

Download PDF

Version 1

posted

You are reading this latest preprint version

Virome compositions indicate that viral spillover is a dead-end between the western honey bee and the common eastern bumblebee

Status:

Version 1

Abstract

Figures

Introduction

Results & Discussion

Materials and Methods

Field sampling

Dissociation and RNA extraction

cDNA synthesis

Sequencing

Contig assembly

Viral Sequence Identification

Phylogeny

Metagenomic analyses

Network analyses

Statistical Analysis

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1