Insect pollination is essential for the global agricultural industry, and the European honey bee (Apis mellifera) is a vital contributor. In the US honey bee pollination has been valued at over $50 billion annually1,2. However, in recent years colony health has globally declined, raising concerns about the future of food security. The Bee Informed Partnership established a long-term surveillance program of commercial operations in the US and produces annual reports on US beekeeper’s experiences which has previously included the 40.9–46.1% colony loss between 2019 and 20223,4. Many factors can contribute to colony loss, such as bee genetics, pathogen loads, and pesticide levels5,6,7,8.

Thirteen viral families have been ‘associated’ with A. mellifera, including diverse RNA viruses9,10. Prior nationwide surveys of US honey bee viruses have predominantly relied on PCR-only studies11,12,13 which is cost-effective and accessible, but has limited capability to identify emerging viruses or changes in existing variants. RNA viruses have extremely high mutation rates and new variants emerge, either due to genome recombination/reassortment, selection, or the accumulation of point mutations due to the highly error-prone RNA-dependent RNA polymerase (RdRp)14. As the genotypic distribution of the virus shifts as a result of RNA evolution, PCR primers can lose sensitivity, which is reported in human viruses such as SARS-COV-215. A multi-year meta-transcriptomic survey of over 2000 viromes from China during 2016–2019 identified 23 novel viruses from both honey bees and mites16, demonstrating one of the many benefits of conducting meta-virome studies. Should the commercial honey bee industry continue to experience colony health issues, the genomic resolution at a strain level divergence will become vital to monitoring and tracking honey bee diseases across the US, and allow for new PCR primer designs utilizing the most recently available genomes.17

Deformed wing virus (DWV, Iflavirus aladeformis) is of particular concern in the US and is regarded as a global pathogen of honey bees18. DWV is vectored by ectoparasitic mites belonging to the genus Varroa19, with Varroa destructor being the predominant honey bee mite in the US20. Furthermore, the presence of V. destructor has led to changes in the prevalence of specific DWV master variants (type A and B), through unknown mechanisms19,21. As each variant is capable of co-infecting and replicating within the same host cell, this creates an opportunity for recombination events which are known to be involved in the generation of novel variants in the picorna-like viruses22. Type B arrived much later in the US than type A17, and until recently, type A was the most common variant in the US and potentially worldwide12,23, as reported from RdRp RT-PCR amplicon surveys. Then between 2016 and 2019 there was an observed increase in type B and a decrease in A17,24.

In this study, we aimed to examine the total genomic diversity of RNA viruses in both A. mellifera and V. destructor and identify any recombinant DWV genomes within commercial beekeeping operations across major hubs in the US. Our findings support prior work, confirming the ongoing displacement of DWV type A variants by type B, but additionally suggest that genomes displaying some degree of recombination are incredibly common.

A. mellifera and V. destructor viromes

Pooled samples of A. mellifera (n = 30) and, when available, V. destructor (n = 1–10), representing a single colony were collected from 55 apiaries across 11 states in the US (Fig. 1a). We generated 555 cDNA libraries for Oxford Nanopore Technologies sequencing, comprising 383 A. mellifera and 173 V. destructor viromes (Supplementary Data Table 1). Average mite densities in the colonies were low, 4.23 mites per 100 bees (standard deviation (s.d.) 土 5.33, Supplementary Data Table 1), limiting our access to V. destructor from all colonies. We generated 71.8 million sequencing reads, with each virome yielding an average of 128,322 reads (s.d. 土 88,930), with an average length of 415 bp long (s.d 土 124 bp). De novo assembly generated 13,233 contigs, averaging 22 contigs per A. mellifera virome (s.d. 土 12.2) and 27 per V. destructor virome (s.d. 土 13.4, Supplementary Data Table 1). Most reads were host transcripts (82–88%), while viral reads accounted for an average of 0.237% (Fig. 1b). Of the 1747 binned viral contigs, there were 45 complete and 1702 partial genomes (236 high-quality, 489 medium-quality, 789 low-quality, and 188 unassigned) as assessed by CheckV25 (Supplementary Data Table 2). To ensure the inclusion of shorter genome fragments, we clustered viral genomes using 95% sequence identity across 90% of the genome, retaining the longest contig26. This generated 997 viral unique clusters. Only genomes of medium and higher completion levels were considered for most of the downstream analysis unless specified.