A Diverse Range of Novel RNA Viruses in Geographically Distinct Honey Bee Populations

Abstract

Understanding the diversity and consequences of viruses present in honey bees is critical for maintaining pollinator health and managing the spread of disease. The viral landscape of honey bees (Apis mellifera) has changed dramatically since the emergence of the parasitic mite Varroa destructor, which increased the spread of virulent variants of viruses such as deformed wing virus. Previous genomic studies have focused on colonies suffering from infections by Varroa and virulent viruses, which could mask other viral species present in honey bees, resulting in a distorted view of viral diversity. To capture the viral diversity within colonies that are exposed to mites but do not suffer the ultimate consequences of the infestation, we examined populations of honey bees that have evolved naturally or have been selected for resistance to Varroa This analysis revealed seven novel viruses isolated from honey bees sampled globally, including the first identification of negative-sense RNA viruses in honey bees. Notably, two rhabdoviruses were present in three geographically diverse locations and were also present in Varroa mites parasitizing the bees. To characterize the antiviral response, we performed deep sequencing of small RNA populations in honey bees and mites. This provided evidence of a Dicer-mediated immune response in honey bees, while the viral small RNA profile in Varroa mites was novel and distinct from the response observed in bees. Overall, we show that viral diversity in honey bee colonies is greater than previously thought, which encourages additional studies of the bee virome on a global scale and which may ultimately improve disease management.IMPORTANCE Honey bee populations have become increasingly susceptible to colony losses due to pathogenic viruses spread by parasitic Varroa mites. To date, 24 viruses have been described in honey bees, with most belonging to the order Picornavirales Collapsing Varroa-infected colonies are often overwhelmed with high levels of picornaviruses. To examine the underlying viral diversity in honey bees, we employed viral metatranscriptomics analyses on three geographically diverse Varroa-resistant populations from Europe, Africa, and the Pacific. We describe seven novel viruses from a range of diverse viral families, including two viruses that are present in all three locations. In honey bees, small RNA sequences indicate that these viruses are processed by Dicer and the RNA interference pathway, whereas Varroa mites produce strikingly novel small RNA patterns. This work increases the number and diversity of known honey bee viruses and will ultimately contribute to improved disease management in our most important agricultural pollinator.

Keywords: RNA interference; arthropod vectors; insect viruses; metagenomics; negative-strand RNA virus; phylogenetic analysis; plus-strand RNA virus.

FIG 1

Lake Sinai virus variants identified in this study. (A) Genome structures of LSV strains identified in The Netherlands (GenBank accession number KY354242), Tonga (accession number KY354241), and South Africa (accession numbers KY354243 to KY354244), compared to previously characterized genomes of LSV-1 and LSV-2. Open reading frames are blue, and conserved functional domains are indicated (NCBI protein sequence accession numbers ARO50053 to ARO50067). (B) Maximum likelihood phylogenetic tree of nucleotide alignment of LSV strains from The Netherlands, Tonga, and South Africa with LSV-1 and -2 and other strains described previously (26).

FIG 2

Genome structures of novel viruses. (A) Genome structures of rhabdo-like viruses, showing the genome size (nucleotides) and rhabdovirus open reading frames (N, P, M, G, and L/RdRp proteins) of ARV-1 (GenBank accession numbers KY354230 to KY354233) and ARV-2 (accession number KY354234) relative to the structure of the previously characterized Farmington virus. (B) Genome structure of bunya-like viruses, showing the identified L segment sizes (nucleotides) and the ORFs of ABV-1 (GenBank accession number KY354236) relative to LepmorLBV1 and of ABV-2 (accession number KY354237) relative to Wuhan mosquito virus 1. (C) Genome structure of a flavi-like virus, AFV (GenBank accession number KY354238), showing the genome size (nucleotides) and ORF relative to GKaV. (D) Genome structure of a dicistro-like virus, ADV (GenBank accession number KY354239), showing the genome size (nucleotides) and two ORFs of ADV relative to Drosophila C virus. (E) Genome structure of a Nora-like virus, ANV (GenBank accession number KY354240), showing the putative 5′-truncated genome size (nucleotides) and four ORFs relative to Drosophila Nora virus. (The NCBI protein sequence accession numbers are ARO50020 to ARO50052.)

FIG 3

Evolutionary relationships of novel viruses. Shown are maximum likelihood phylogenies of the novel rhabdoviruses ARV-1 and ARV-2 (A), the novel bunyaviruses ABV-1 and ABV-2 (B), the novel flavivirus AFV (C), the novel dicistrovirus ADV (D), and the novel Nora virus ANV (E). (See Fig. S1 to S3 in the supplemental material for detailed trees for panels A to C.)

FIG 4

Apis mellifera rhabdovirus 1 variants. Shown is a maximum likelihood phylogenetic tree of the nucleotide alignment of ARV-1 variant genomes isolated from The Netherlands (GenBank accession number KY354230), South Africa (accession number KY354231), and Tonga (accession number KY354232).

FIG 5

Small RNA analysis of ARV-1 and ARV-2 in honey bees. (A and B) Size distributions (15 to 37 nt) and 5′-nucleotide compositions of small RNAs arising from ARV-1 (A) and ARV-2 (B). The sample from which the small RNA library was produced is shown above each graph. Bars plotted above the x axis represent reads that map to the positive strand, and those plotted below represent those that map to the negative strand. Bars are colored according to the proportions of reads starting with A, C, G, and T. (C and D) Mapping of 20- to 23-nt-long viral RNAs (viRNAs) to the genomes of ARV-1 (C) and ARV-2 (D). The cartoon shows the domains of the viral genomes as shown in Fig. 2. (E and F) Phasing score over 8 phasing cycles for each position within a 22-nt phasing window. The top and bottom graphs show the phasing scores for the sense and antisense reads, respectively. A high score indicates that many small RNAs fall into that phase position (indicated with arrowheads). This analysis was performed by using the 21- to 23-nt-long reads from panels A and B. (G and H) Observed 5′ nucleotide (Obs) compared with that expected (Exp) from the base compositions of the viral genomes for ARV-1 (G) and ARV-2 (H). Sense (S) and antisense (AS) reads were compared by using a chi-squared test. **, P value of <0.01.

FIG 6

Small RNA analysis of ARV-1 and ARV-2 in Varroa. Left panels show ARV-1, and right panels show ARV-2. (A and B) Size distributions (15 to 37 nt) and 5′-nucleotide compositions of small RNAs in mite 1 (M1) arising from ARV-1 (A) and ARV-2 (B). Bars plotted above the x axis represent reads that map to the positive strand, and those plotted below represent those that map to the negative strand. Bars are colored according to the proportions of reads starting with A, C, G, and T. (C and D) Mapping of 23- to 25-nt-long viRNAs to the genomes of ARV-1 (C) and ARV-2 (D). The cartoon shows the domains of the viral genomes as shown in Fig. 2. (E and F) Size distributions and 5′-nucleotide compositions of the sense small RNAs from panels A and B, respectively, normalized to the number of sense reads present. Mapping of the 19- to 24-nt sense reads to the viral genomes is also shown. (G and H) Phasing scores over 8 phasing cycles for each position within a 24-nt phasing window for ARV-1 (G) and ARV-2 (H). The top and bottom graphs show the phasing scores for the sense and antisense reads, respectively. This analysis was performed by using the 24-nt-long reads only. (I and J) Observed 5′ nucleotides (Obs) compared with those expected (Exp) from the base compositions of the viral genomes of ARV-1 (I) and ARV-2 (J). Sense (S) and antisense (AS) reads were compared by using a chi-squared test. *, P value of <0.05; ns, not significant. (K and L) Distances between the 5′ ends of overlapping reads on opposite strands (left) and the base compositions of each nucleotide position (right) for the 23- to 25-nt-long reads of ARV-1 (K) and ARV-2 (L).