Investigating the viral ecology of global bee communities with high-throughput metagenomics

Bee viral ecology is a fascinating emerging area of research: viruses exert a range of effects on their hosts, exacerbate impacts of other environmental stressors, and, importantly, are readily shared across multiple bee species in a community. However, our understanding of bee viral communities is limited, as it is primarily derived from studies of North American and European Apis mellifera populations. Here, we examined viruses in populations of A. mellifera and 11 other bee species from 9 countries, across 4 continents and Oceania. We developed a novel pipeline to rapidly and inexpensively screen for bee viruses. This pipeline includes purification of encapsulated RNA/DNA viruses, sequence-independent amplification, high throughput sequencing, integrated assembly of contigs, and filtering to identify contigs specifically corresponding to viral sequences. We identified sequences for (+)ssRNA, (−)ssRNA, dsRNA, and ssDNA viruses. Overall, we found 127 contigs corresponding to novel viruses (i.e. previously not observed in bees), with 27 represented by >0.1% of the reads in a given sample, and 7 contained an RdRp or replicase sequence which could be used for robust phylogenetic analysis. This study provides a sequence-independent pipeline for viral metagenomics analysis, and greatly expands our understanding of the diversity of viruses found in bee communities.

in areas with increased density of infected managed colonies [12][13][14][15][16] . Viruses found to be pathogenic in A. mellifera (such as deformed wing virus, Israeli acute paralysis virus, acute bee paralysis virus, and Kashmir bee virus) have also been found to be pathogenic in wild bee species 14,17,18 . Targeted studies have identified subsets of the previously identified honey bee viruses (including deformed wing virus, black queen cell virus, sacbrood virus, Israeli acute paralysis virus, acute bee paralysis virus, Lake Sinai virus, chronic bee paralysis virus, and Kashmir bee virus) in bee populations in Asia 19 , Africa 20 , the Middle East 21 , Australia 22 , and South America 23,24 , suggesting a world-wide distribution of these viruses.
Since the most well-studied viruses were originally identified in A. mellifera, it has often been assumed that A. mellifera are the primary hosts and other bees are only secondarily infected. However, because many viruses seem to be so readily shared, it is not possible to determine which bee species, if any, serves as the primary host. Indeed, based on sequence analyses, the viral species that have been examined thus far seem to be shared within a geographic region rather than within a particular bee species, suggesting that there may not necessarily be bee species-specific viral strains 13,14 . Thus, it is likely that non-Apis populations also harbor viruses which can spill over into A. mellifera.
Similarly, recent metagenomics studies of bees have identified many viruses that were originally described in plants 25,26 . It has been assumed that these are simply contaminants from the pollen and nectar in the bees' diets, but these viruses may in fact be using bees as hosts. Indeed, tobacco ringspot virus (TRSV) was recently shown to infect and replicate in A. mellifera 27 . Further studies are needed to fully examine the intriguing possibility that pollinators are more than simply carriers of plant viruses.
Thus, fully understanding the viral community in managed and wild bee species, and the plants they feed on, requires a holistic approach that examines these communities in diverse geographic regions. Distinct types and species of viruses undoubtedly remain unidentified in populations of A. mellifera and other bee species across the world. Moreover, due to increased globalization, it is likely only a matter of time before these viruses (and other pathogens and parasites) are spread to new regions. For example, recently a colony of an exotic stingless bee species (Plebeia frontalis) was found to be established in California 28 , potentially harboring pathogens that could be spread to other bee species within the region.
With the development and optimization of inexpensive high-throughput sequencing technologies, it is now possible to sequence the genetic material of the entire community of viruses present within a host, without any prior knowledge of the viral genome sequences (viral metagenomics analysis) 29 . This technology allows rapid and efficient identification of previously uncharacterized types and strains of viruses, bypassing labor-intensive steps of purifying and culturing the viruses. Viral metagenomics has been utilized by many studies seeking to identify novel viruses and characterize viral communities in diverse organisms (reviewed in Mokili et al. 30 ).
Using this approach, several recent studies have identified novel viruses in populations of bees (A. mellifera, Bombus terrestris, Bombus pascuorum, Osmia cornuta, and Andrena vaga) collected from sites in the Netherlands, South Africa, Tonga, Israel and Belgium 6,26,31 . However, while these studies have provided valuable information about viral communities within particular species or regions, they are limited in scope, in that they focus on a single bee species in multiple locations (A. mellifera populations in Israel 31 or the Netherlands, South Africa, and Tonga 6 ) or a selection of wild bees in a single region (Belgium 26 ).
To improve our ability to detect viruses, we have developed a protocol for unbiased, sequence-independent detection of novel viruses within diverse bee populations from around the globe. Additionally, we used a viral enrichment protocol to purify the viral particles with the aim of increasing the proportion of sequence reads obtained from viruses versus the host bee species, compared to previous studies that did not use this protocol 6,26 . We used this approach to characterize the viral communities of 12 bee species collected from 9 countries across 4 different continents and Oceania (see Table 1 for bee species and their respective locations). This approach allowed us to search for novel viruses, begin to determine the distribution of these viruses (and viral families) across the world, and begin to evaluate the extent to which A. mellifera and wild bee populations share these viruses. As in other studies 25,26 , we identified several viruses that share sequence homology to plant viruses and thus may be the result of contamination, but given that these were found in multiple populations and locations, it raises the question of whether or not these viruses may use bees as hosts.

Methods
Sample collection. Sample collection kits and standardized protocols were provided to collaborators to obtain samples of foraging bees from geographically distant locations. For each geographic region, collaborators were instructed to collect western honey bees (Apis mellifera) and other common bee species (including both Apis and non-Apis bee species) as individuals foraged on flowering plants; no attempt was made to preferentially collected individuals that appeared diseased. Individuals from different bee species were collected in separate tubes to ensure that the species were not mixed. At least 20 individuals were collected onto ice-cold 95% ethanol for each species from each region. Details for each site were recorded, including the GPS coordinates. A summary of different locations and sampled species can be found in Table 1. The samples were then shipped to Pennsylvania State University for processing. We also included samples of A. mellifera with deformed wings (a symptom of high levels of DWV) collected from a colony at the Penn State University apiaries as a control for DWV detection (labeled Positive Control in Table 1).
Virus purification and nucleic acid extraction. 10 individuals from each sample group (representing a specific species collected at a specific location) were placed in individual 2.0 ml microcentrifuge tubes with 150 µl nuclease free H 2 O and 3-5 sterile glass beads and homogenized in a FastPrep FP120 cell disruptor (Thermo-Fisher Scientific, Waltham, MA, USA) for 45s. 100 µl of the homogenate from each of the 10 individuals was pooled in a common 1.5 ml microcentrifuge tube. Due to their large size, 10 Xylocopa spp. samples were pooled and homogenized using a mortar and pestle in 5 ml nuclease free H 2 O. The pooled homogenate from SCIENTIFIC REPORtS | (2018) 8:8879 | DOI:10.1038/s41598-018-27164-z each sample was then dialyzed using a 0.2 µm cell filter (Corning, Tewksbury, MA, USA) to purify the virus extracts as in Hunter et al. 32 and Liu et al. 33 . A 125 µl aliquot of each sample was treated with a nuclease cocktail (including 14U Turbo DNase I, 25U Benzonase, 20U RNaseI, and 10X DNase buffer) and incubated at 37 °C for 1.5 hours to remove all nucleic acid that was not protected by a viral capsid as in He et al. 34 Nucleic acids (both RNA and DNA) from the purified encapsulated viruses were extracted using a MagMAX Viral Isolation Kit (Thermo-Fisher Scientific, Waltham, MA, USA), according to the manufacturer's protocol.
Sequence independent nucleic acid amplification and sequencing. A previously established protocol was used to obtain unbiased random amplification of the extracted nucleic acids 35 . Briefly, the first strand cDNA was synthesized from 30 μl of viral nucleic acid using a random primer design (consisting of a 20 base oligonucleotide sequence followed by a randomized octamer sequence: GACCATCTAGCGACCTCCACNNNNNNNN) in a reverse transcriptase (RT) reaction using a High-Capacity cDNA Reverse Transcription Kit (Thermo-Fisher Scientific, Waltham, MA, USA). To synthesize the second strand, the 19 μl of the initial cDNA was denatured at 95 °C for two minutes and cooled to 4 °C prior to the addition of Klenow fragment DNA polymerase (New England Biolabs, Ipswich, MA, USA) for fragment extension at 37 °C for 60 minutes. Finally, PCR amplification was performed for 40 cycles using the above primer without the randomized octamer sequence (GACCATCTAGCGACCTCCAC) and 2 μl of dsDNA template. The PCR reactions were then purified using a MSB Spin PCRapace Purification Kit (Invitek, Berlin, Germany) using the manufacturer's standard protocol. Quality and concentration of the samples was assessed using a 2100 Bioanalyzer (Agilent, Santa Clara, CA). The samples were then submitted to the Genome Core Facility at Penn State University for sequencing on the Illumina MiSeq, resulting in 1,056,165 (+/− 105,794) reads at 150 bases per sample (see Supplementary Table 1 for more information on the number of reads and assembled contigs in each sample). Although the phenomenon known as "index switching" has been observed to misattribute reads to the wrong sample in multiplexed sequencing libraries in the same run on certain Illumina platforms, this effect is demonstrated to be limited on the MiSeq and is not expected to generate increases in external contaminants 36 . Raw data from this Whole Genome Shotgun project have been deposited at DDBJ/ENA/GenBank under the accession SAMN07695258, along with BioProject (PRJNA407112) BioSamples (SAMN07634920-SAMN07634956).
Data processing. The quality of raw sequencing reads was first manually inspected using FastQC 37 . Illumina adaptor sequences and reads with low quality scores (25) were removed using Trimmomatic v0.30 38 . As the randomized octomer sequences correspond to regions of viral genome and not every read contains them, we did not trim the start of each read, although we note that polymorphisms in these regions may possibly complicate the assembly. To identify novel viruses, the trimmed and filtered sequencing reads were assembled de novo using three different programs: Oases v0.2.8/velvet v1.2.10 (N50: 310 nt) 39 , Trinity v2.1.0 (N50: 326 nt) 40 , and SPAdes v3.7.1 (N50: 315 nt) 41 . The N50 (weighted median length of the contigs) for all assemblies were low (the total read length was 150 nt), so the three assemblies for each sample were merged using the software package Minimus2 (AMOS v3.1.0) 42 , resulting in an improved N50 of 699 nt. For comparison, this N50 is greater than for similar previous studies of A. mellifera and Culex pipiens (mosquito) associated viruses 6,43 , yet we do not exclude the possibility that further optimization of Kmer sizes could yield slight increases in overall contig length. For the three assembly methods, we used Kmer sizes of 21 (Oases/velvet) and 25 (Trinity), while SPAdes uses multiple kmer sizes (21, 33, and 55) to identify the optimal assembly. The assembled contigs were then subjected to a series of processing steps. To remove contamination and retain a set of candidate contigs likely to be viral in origin, each assembly was screened for homologous sequences from three potential sources of contamination: the human genome (version hg38), all bacterial genomes in the NCBI database, and the honey bee genome (version amel4.5) using an E-value cutoff of 10 −5 .The total set of contigs was screened for sequence homology to the three possible sources of contamination using BLAST 44 . Contigs with significant hits to any of the contamination sources were then removed from further analyses (see Supplementary  Table 1 for the number of contigs matching each source of contamination). To identify the closest known viral sequence homology, the remaining assembled contigs were compared to all sequences in the NCBI nr database 45 using DIAMOND v0.8.36 46 with an E-value cutoff of 10 −5 . Viral families were assigned using the 2017 report of the International Committee on Taxonomy of Viruses (ICTV) for the classification and nomenclature of viruses 47 (see Supplementary Table 2).
To determine relative abundances of each virus, the trimmed sequence reads from each sample were then mapped to the set of candidate viral contigs identified in that sample using Tophat2 v2.0.9 48 (see Supplementary  Table 2). Contigs with >0.1% of the total reads successfully mapping in that sample were chosen for further analyses (see Supplementary Table 3). For the contigs with >0.1% of the total reads, ORFfinder 49 was used to identify potential open reading frames within the contigs that had significant BLAST hits (E-value < 1 × 10 −5 ) to known viruses to further characterize the viral contigs from each sample. Viral contigs that showed evidence of a putative RNA dependent RNA polymerase (RdRp) or other replicase proteins were selected for further phylogenetic analyses. Note that while we used a sequence independent amplification approach to generate the libraries for sequencing, amplification bias may still have skewed the total numbers of reads for the viral sequences.

Phylogenetic analysis.
Reference sequences corresponding to the RdRp or other replicase proteins were downloaded from NCBI GenBank 50 for viruses that were closely related to the significant BLAST hit to identify the relationship of the novel viruses discovered here to other known viruses. Multiple amino acid sequence alignments were performed using MUSCLE 51 in MEGA7 52 with default settings. Phylogenetic trees for each novel virus identified in Fig. 1 were created using the maximum-likelihood method with 500 bootstrap replications in MEGA7 52 , by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using a JTT model, and then selecting the topology with the superior log likelihood value.

Results
Identification of Viral Sequences. On average, 6.87% of assembled contigs in each sample had evidence of being viral in origin (see Supplementary Table 1). Our enrichment protocol is aimed to maximize the amount of viral nucleic acids in each sample, which are typically far fewer than found in host or other microbes 53 . Indeed, in previous studies of viruses found in bee species, the total proportion of reads determined to originate from viruses range from 0.1 to 2.7% 6,26 . However, we recognize further improvements are needed to increase the abundance of viral genomic material in metagenomics studies of bees, and note that in smaller insect species such as Ixodes scapularis ticks the total proportion of viral reads has been reported to be as high as 35% 54 .
With our approach, we identified 479 different contigs among the 37 samples that shared sequence homology with previously identified viral sequences (see Supplementary Table 2 for information on each of the contigs identified in each sample). Of these, 308 contigs correspond to "novel" viruses (defined as viruses that were not previously found in bees), and 223 contigs correspond to known bee viruses. Many of these contigs were identified in multiple locations and/or species. When duplicate/overlapping contigs are removed, the number of unique novel virus contigs was reduced to 127. Next, we examined the number of contigs that were supported by >0.1% of the total reads for the given sample in which it was found, as these may represent more prevalent or abundant viruses in bee populations. Only 27 contigs matched these criteria (see Supplementary Table 3). Sequence homology analyses of these 27 novel virus contigs revealed 9 different viral families among these bee associated viruses, including 6 viral families that had not been heretofore observed in bees (noted in Supplementary Table 3). Of these, one virus from the family Rhaboviridae (BRV) was recently described in A. mellifera populations in the United States, Israel, South Africa, Tonga and the Netherlands, and B. impatiens populations in the United States 6,55 . Of the 27 novel virus contigs, 7 were found to contain a putative RdRp or replicase that could be used for phylogenetic analyses.
The distribution and relative quantities of the known viruses and newly identified viruses with RdRp domains (that were present at >0.1% of the reads in each sample) across each sample are shown in Fig. 1A,B provides information on the geographic distribution of these viruses and their presence in A. mellifera samples versus non-A. mellifera samples at each site. Details about each of these novel viruses, including their phylogenetic characterization, are presented below. Interestingly, 5 novel viruses were identified in North American samples, 4 in New Zealand, and 2 in each of the other continents; thus, though the North American samples are relatively well-studied compared to the other regions, at least in our analysis, these samples still generated more novel viruses. In regions where both A. mellifera and non-A. mellifera bees were sampled, A. mellifera had higher numbers of viruses in four regions (Kenya, Nicaragua, Switzerland, and Pennsylvania, US), lower numbers in two regions (India and Washington, US) and equivalent numbers in Panama. Only in the United States and India were  Table 2 for more details on the locations and species in which these viruses were found).

Newly Identified Positive Sense RNA Viruses. Dicistro-like virus.
We identified one Dicistro-like contig that shared sequence homology with known viruses from the family Dicistroviridae, a viral family that includes several known bee-infecting viruses, such as IAPV, BQCV, and Aphid lethal paralysis virus (ALPV). This virus was found in Apis florea samples from India, with a contig length of 3,818 nt, which included a putative RdRp region, as well as in A. mellifera samples from New Zealand with a length of 1,382 nt. A phylogenetic analysis of this region revealed that this particular contig was most closely related to aphid lethal paralysis virus (ALPV), and forms a clade with ALPV and other dicistroviruses identified in a previous study using other invertebrates, including odonates, heteropterans, myriapods, and orthopterans ( Fig. 2A) 56 .

Seco-like virus.
We identified a Seco-like virus that shares sequence homology to Secoviridae, which is in the same viral order as Dicistroviridae (Picornavirales). Overlapping viral contigs were found in four bee species across 4 continents and Oceania (A. mellifera in Georgia: 1,988 nt, Kenya: 971 nt, Panama: 1,295 nt, Nicaragua: 1,955 nt, Switzerland: 1,105 nt, and New Zealand: 1,799, Megachile rotundata in the US (Washington): 1,500 nt, Noda-like virus. We identified a novel (+)ssRNA virus within the family Nodaviridae in two different species, A. mellifera (1,216 nt) and B. impatiens (1,831 nt). Interestingly, these were from two distant geographic locations (Hamilton, New Zealand and Pennsylvania, United States). Using a BLASTx 44 approach on the open reading frames identified within the assembled contigs, we identified a putative RdRp that was used for phylogenetic analysis. The phylogenetic analysis revealed that this virus was most closely related to other viruses within the Nodaviridae family of viruses, forming a clade with members of the viral genus Alphanodavirus (Fig. 2C), which includes several insect viruses, including flock house virus, black beetle virus, and nodamura virus. viruses been identified in bees 6,26 . One of these, Apis mellifera Rhabdovirus 1 (ARV-1), from the family Rhabdoviridae 6 , was detected in A. mellifera populations in South Africa, Tonga, and the Netherlands. It was also recently identified in bee populations in Israel and the United States 55 . Our data sets related to this virus were published elsewhere 55 and thus we will only briefly summarize this information here.

Tymo-like virus. A Tymo-like virus with overlapping contigs was identified in
We identified contigs corresponding to this Rhabdovirus in two different species (A. mellifera and Bombus impatiens) that were collected in Pennsylvania, United States, suggesting this virus can easily switch hosts. Thus, Levin et al. 55 proposed a different nomenclature, Bee Rhabdovirus (BRV). The two assembled contigs for each species were 13,844 nt and 14,589 nt in A. mellifera and B. impatiens respectively, and were >99% similar to ARV-1. These assembled contigs contained the 5 genes that are typically carried by Rhabdoviruses, including the large protein (L), a glycoprotein (G), nucleoprotein (N), phosphoprotein (P), and matrix protein (M). Further information about this virus can be found in Levin et al. 55 .

Newly Identified Double Stranded RNA Viruses. Partiti-like virus.
Overlapping contigs with sequence homology to other viruses within the viral family Partitiviridae were identified in A. mellifera samples obtained from California, United States (1,139 nt), Switzerland (409 nt), and Georgia (385 nt). A phylogenetic analysis of a putative RdRp region from the assembled contigs, determined by a BLASTx 44 screen, revealed sequence homology that is similar to other identified partitiviruses, but most closely related to the newly identified arthropod partitiviruses 56 (Fig. 2E).

Newly Identified Single Stranded DNA viruses. Circo-like virus.
Sequences corresponding to a circo-like virus in A. mellifera from three different locations (Pennsylvania, United States: 521 nt, Kenya: 688 nt, and Nicaragua: 1,531 nt) with sequence homology to viruses from the Circoviridae family. Interestingly, a phylogenetic analysis of the replicase protein revealed two different viruses, with the sequences from the US and Nicaraguan A. mellifera populations being closely related to circular ssDNA viruses found in bats and dragonflies, respectively (Fig. 2F). The Circo-like virus contig that was discovered in Kenya shared sequence homology with bat circovirus (similar to the virus in the Pennsylvania samples), but did not overlap with the replicase protein. Figure 3. Genome of the newly discovered Tymo-like virus. Genome structure of the Tymo-like virus discovered in Bombus impatiens and Apis mellifera from the United States and New Zealand. This viral contig contains a long polyprotein encoding a methyltransferase (MTR), proteinase (PRO), helicase (HEL), and an RNA dependent RNA polymerase (RdRP) that is common among viruses from the Tymovirus genus. Additionally, a separate overlapping protein encoding a movement protein at the 5′ end and coat proteins was detected at the 3′ end of the viral contig. Log read counts were calculated using a sliding window of 100 bases. Confirmation of viruses using PCR. We used PCR to confirm the presence of the expected viral contigs in the samples for BRV (see Levin et al. 55 ), Circo-like 2, Noda-like, Tymo-like, Seco-like, and Dicistro-like viruses (see Supplementary Materials for primers, amplicon size, bee sample set).

Discussion
In this study, we developed a viral metagenomics pipeline to comprehensively examine the viromes of specimens of 12 different bee species collected from a broad array of locations (9 countries, across 4 continents and Oceania). This dataset dramatically enhances our current knowledge of bee associated viruses, by identifying 27 novel (heretofore unidentified in bees) viral contigs from 9 viral families that were present at >0.1% of the reads within a sample. Six of these contigs contained a putative RdRp or replicate protein region (identified by BLASTx 44 and ORFfinder 49 ) that could be used for phylogenetic analyses to determine relationships between previously identified viruses (see Fig. 2). Moreover, using our pipeline, we identified a wide variety of new types of viruses that had not been previously described (or only recently described) in bees, including (−)ssRNA, dsRNA, and ssDNA viruses.
This study represents the largest effort to identify novel pathogens in global bee samples to date. However, it is important to note that this study still had relatively small samples sizes, given the diversity of bee species and their global distributions. Sampling 10 bees per population (as in our study) is sufficient to detect more prevalent viruses within a given area, however, this sampling effort may be limited in the detection of low-prevalence viruses 57 . Indeed, by focusing on viruses that were represented in >0.1% of the reads in a sample, we may have further removed data from low abundance viruses. We were also unable to assemble full length genomes for most of the viruses detected, which may be due to a number of factors, including low sequencing coverage, uneven sequence coverage across viral genomes, and low sequence homology to sequences in current databases. Moreover, while our approach to enrich our samples for encapsulated viruses may have improved our ability to detect viruses by reducing the reads coming from host material (indeed, in 30/37 samples, more contigs matched to the viral databases than the honey bee genome, see Supplementary Table 1), this method involves additional processing steps which may have degraded some types of viruses. Furthermore, our conservative bioinformatics filtering procedure (E-value < 1 × 10 −5 ) retained only a set of contigs with highly significant matches to viruses described in the NCBI nr database, thereby potentially removing sequences with homology to uncharacterized viral families (in only 3/37 samples, more contigs matched the viral databases than all potential sources of contamination and unmapped sequences combined). To limit the amount of genomic material generated from microbial contaminants, it may also be possible to avoid using midgut tissue which likely harbors a greater number of bacteria than other tissues. Thus, additional screens using larger sample sizes, both pooled and individual samples, alternative extraction techniques, and greater sequencing depth may help improve the recovery of larger contigs, potentially identifying additional viruses and likely demonstrating a broader geographic and species range for the viruses we identified in this study.
Additionally, the goal of this study was to identify viruses that may be circulating among co-foraging bee populations. Thus, we exclusively sampled foraging bees in the landscapes in the summer months. Viral titers can vary greatly among bees in different developmental stages, behavioral states, or seasons 8 . Moreover, highly diseased bees exhibit reduced foraging 58 . Additional studies using a broader diversity of bee samples -including samples of bees exhibiting clear symptoms of disease -may result in a broader array of identified viruses.
It is important to note that we did not assess replication of these viruses in our samples, and this will be necessary to confirm that bees are indeed hosts of these viruses. Viruses may be present in bee samples due to contamination from the environment or their diet; indeed, viruses can infect or contaminate pollen grains and be transmitted to individuals foraging on or consuming this pollen 12,59 . We selected viruses that were present at >0.1% of the reads within a sample (an average of ~1,000 reads per sample), to reduce the likelihood that we were examining a trace contaminant. Additionally, all of the viruses were identified in multiple species and multiple sites in independent collections and extractions, and thus it seems somewhat unlikely that these viruses were simply contaminants. Moreover, recent screens of samples collected as part of the US National Honey Bee Disease Survey 60 have identified several of these viruses in samples from across the country (Ray et al., in preparation).
Positive sense RNA ((+)ssRNA) viruses are the most common type of virus that have been identified in bees thus far, with a vast majority of the viruses possessing this genome architecture. Most of these viruses belong to the viral order Picornavirales (including DWV, BQCV, IAPV, and SBV), but recently additionally studies have added to the list of viral orders that infect bees, including Tymovirales 6,26 . Here we have identified novel virus sequences that share sequence homology with three (+)ssRNA virus families, Dicistroviridae, Secoviridae, and Tymoviridae.
Dicistroviruses are commonly found in honey bee colonies, with the most common viruses being IAPV and BQCV 61 . Unsurprisingly, we found BQCV to be the most widely distributed virus in our dataset, while IAPV was also fairly common, following only behind BQCV and DWV. Another recent study also identified a novel dicistrovirus in A. mellifera from Amsterdam 6 , however this virus is phylogenetically different than the virus we described here. We found this viral contig (containing a putative RdRp) in A. florea from India, however, we also detected another viral contig from A. mellifera from New Zealand that shared sequence homology with a different region of the same virus (Hubei picorna-like virus 15; accession number: NC_032757.1). Unfortunately, there was no overlap between these two contigs, so we cannot say for certain whether these two contigs are a part of the same virus or if they represent two different Dicistroviruses.
We identified a contig corresponding to the Secovirus family in four bee species collected from 4 continents and Oceania. This virus shares its closest homology with viruses detected primarily in plant species 62 . While additional studies are needed to determine if this virus can replicate in bees, its broad distribution suggests that it is closely associated with bee species. Other recent studies have identified viruses with sequence similarity to plant viruses in bee samples 26,27 . In one case, for tobacco ringspot virus, replication was confirmed in bee samples 27  remains to be determined if these viruses are truly plant-associated viruses that simply contaminate bees or if in fact bees are an important host. Regardless, the fact that bees carry plant-associated viruses so frequently suggests that they may be important vectors of plant diseases in the field. Indeed, recent studies found that A. mellifera are able to distribute plant viruses through pollination services 63 .
A study by de Miranda et al. characterized the first virus from the viral family Tymoviridae in A. mellifera honey bees, and also found a virus from this family replicating in Varroa destructor mites (Bee macula-like virus and Varroa Tymo-like virus, respectively) 5 . Our phylogenetic analysis suggested that the virus we observed is distinct from the previously described bee virus in this family, and is more closely related to viruses within the Tymovirus genus (as opposed to the Maculavirus genus for the other previously identified bee virus in this family). We found evidence of this Tymo-like virus in two distant geographic locations (New Zealand and United States), suggesting that this virus is widespread. Additionally, we found this virus in two different bee species (Apis mellifera and Bombus impatiens) indicating that this virus can infect both Apis and non-Apis bees.
Until recently, none of the viruses found to infect bees were (−)ssRNA viruses. However, Remnant et al. found four novel (−)ssRNA viruses, including a novel Rhabdovirus in A. mellifera honey bees from South Africa, Tonga, and the Netherlands 6 . This group proposed the name Apis mellifera Rhabdovirus (ARV). We also detected this same Rhabdovirus (>99% similarity) in Apis mellifera and Bombus impatiens samples from the US (Pennsylvania). Data on these sequences and detection of the virus in Israel honey bee populations has recently been published 55 . Given that the virus is found in two bee species, we recommend referring to it as Bee Rhabdovirus (BRV).
Partitiviruses have a dsRNA genome, and have been observed mainly in plants and fungi 64 . A recent study examining the invertebrate RNA virosphere 56 identified a number of Partitiviruses within animal hosts from a wide range of arthropod species, including insects. A phylogenetic analysis of our novel Partitivirus revealed that it is most closely related to this newly identified group of Arthropod Partitiviruses. The distribution of this virus, as detected here, suggests that it is also widespread, being found in Europe (Switzerland and Georgia) and the US (California). To the best of our knowledge, this is the first Partitivirus identified in bees, and further, this is the first evidence of a dsRNA virus in bees.
Viruses that utilize a DNA genome are very rare in bees, with very few being observed to this point 26,65 . Here, we identified a viral contig that is closely related to other members of the Circoviridae family of viruses, a viral family that utilizes a circular single stranded DNA (ssDNA) genome. These viruses have been found primarily within birds and mammals, but recent studies have begun to identify circular viruses within insects as well 66 , including 24 diverse circular viruses found to be infecting odonates (dragonflies and damselflies) 67 . It is possible that bees harbor a similarly diverse and robust set of circular viruses. The identified Circovirus was phylogenetically related to Dragonfly circularisvirus and to the best of our knowledge, represents the first evidence of a circular ssDNA virus observed in bees.
Overall, we found 127 unique contigs corresponding to novel viruses (i.e. viruses previously not observed in bees), with 27 contigs represented by >0.  Table 3), which thus represent important viruses for future analyses. Of these contigs, 7 contained RdRp or replicase regions that could be used for phylogenetic analyses (see Fig. 1). Our study suggests that, despite the focus on North American and European A. mellifera viral populations, these populations still remain to be fully characterized: indeed, the greatest number of newly identified viruses were found in North American samples. Additionally, when considering all the sample sets, it appears that all the newly identified viruses were found in multiple bee species (note that non-overlapping contigs for a novel dicistrovirus found in non-A. mellifera samples in India and A. mellifera samples in New Zealand -whether these represent distinct viruses or the same virus requires further sequence information). Finally, the viruses were largely evenly shared across the A. mellifera and non-A. mellifera bee species in a region, though overall, more sites (4 out of 7) had higher numbers of viruses identified in the A. mellifera samples than in the non-A. mellifera samples.
With our increasing awareness of the importance of viruses to bee health, and the complexity of bee viral ecology -where viruses are shared across insect and potentially plant species -there has been an extensive effort to fully characterize the bee viral community. Prior to 2010, more than 24 viruses were identified in bee populations [2][3][4][5][6][7] . Since then, several additional viruses have been identified using metagenomics approaches 6,26,31 . Together with our results, more than 30 viruses have been found in bees, however, this number likely represents a small subset of the viral diversity present in bee communities. Viral metagenomics for virus discovery in bees has thus far been fruitful, and our protocol has provided a foundation for future studies to continue identifying novel pathogens infecting global bee populations using an inexpensive, unbiased method for the sequence-independent detection of novel viruses. Additional studies are needed to confirm that these viruses replicate within bees, obtain the full viral genome sequences, examine any negative (or positive) impacts infection may have on bees, and examine how these viruses move within bee and plant communities. Our results suggest, however, that viral communities are generally shared globally across species, with species from geographically distant locations hosting similar suites of viruses and viral families.