Area Wide Monitoring of Plant and Honey Bee (Apis mellifera) Viruses in Blueberry (Vaccinium corymbosum) Agroecosystems Facilitated by Honey Bee Pollination

Healthy agroecosystems are dependent on a complex web of factors and inter-species interactions. Flowers are hubs for pathogen transmission, including the horizontal or vertical transmission of plant-viruses and the horizontal transmission of bee-viruses. Pollination by the European honey bee (Apis mellifera) is critical for industrial fruit production, but bees can also vector viruses and other pathogens between individuals. Here, we utilized commercial honey bee pollination services in blueberry (Vaccinium corymbosum) farms for a metagenomics-based bee and plant virus monitoring system. Following RNA sequencing, viruses were identified by mapping reads to a reference sequence database through the bioinformatics portal Virtool. In total, 29 unique plant viral species were found at two blueberry farms in British Columbia (BC). Nine viruses were identified at one site in Ontario (ON), five of which were not identified in BC. Ilarviruses blueberry shock virus (BlShV) and prune dwarf virus (PDV) were the most frequently detected viruses in BC but absent in ON, while nepoviruses tomato ringspot virus and tobacco ringspot virus were common in ON but absent in BC. BlShV coat protein (CP) nucleotide sequences were nearly identical in all samples, while PDV CP sequences were more diverse, suggesting multiple strains of PDV circulating at this site. Ten bee-infecting viruses were identified, with black queen cell virus frequently detected in ON and BC. Area-wide bee-mediated pathogen monitoring can provide new insights into the diversity of viruses present in, and the health of, bee-pollination ecosystems. This approach can be limited by a short sampling season, biased towards pollen-transmitted viruses, and the plant material collected by bees can be very diverse. This can obscure the origin of some viruses, but bee-mediated virus monitoring can be an effective preliminary monitoring approach.


Introduction
Viruses are widespread in plant agricultural systems, but approaches towards virus detection are often biased for symptomatic plants and limited to serological or PCR-based detection assays of individual viruses [1,2]. Randomized field sampling can be used for unbiased early detection of disease agents, but is rarely implemented or limited to

Sample Collection and RNA Extraction
Four different sample types were collected from managed bee colonies temporarily located in two independent blueberry (Vaccinium corymbosum) farms in BC (BC1 and BC2), and one in ON, in the summer of 2021, during blueberry bloom. Queen sources for the bee colonies were locally bred, and were mixtures of Italian and Carniolan-derived stocks. Blueberries were the primary crop grown at the ON site, but other small fruit crops were grown nearby including strawberry (Fragaria x ananassa) and raspberry (Rubus idaeus), while both BC sites were more monocultured. The BC sites were located in the Fraser valley, approximately 4.6 km apart, and both consisted of the variety "Hardy Blue" with BC1 having plants 30 years old, while BC2 had plants that were 20 years old. The four sample types consisted of the following: (1) 25 individual returning forager bees collected outside the hive and contained pollen in the corbicula; (2) approximately 10 mL of pollen collected using pollen traps (ApiHex Beekeeping supplies, Guelph, ON, Canada); (3) 25 individual hive bees (adult workers without visible signs of pollen on their bodies) collected from inside the hive; and, (4) approximately 20 mL of bee bread (stored pollen) collected from individual frames with a sterile spatula. Multiple hives were located at each farm site, and replicates were collected from separate colonies at each location. Four replicates of each sample type were collected in BC sites, while three replicates were collected in ON. In addition, two replicates of mixed flower and leaf samples (plant samples) were collected from 10-20 randomly selected blueberry plants at each site, located near the beehives. Sample names were created based on their location, plant species, and sample type and replicate (Ex. BCFV-BC Fraser Valley, site 1 or 2, or ON for Ontario, followed by crop type BB for blueberries, the sample replicate 1-4, followed by sample type). Forager bees were denoted by F, Pollen-P, Bee bread-B, and hive bee-H. Plant tissue samples were denoted by T). Total RNA (totRNA) was extracted from each hive-related sample using the spectrum total plant RNA extraction kit (Sigma Aldrich, ON, Canada), while dsRNA was extracted from composite plant samples following Kesanakurti et al. (2016) [48].

RNA Sequencing
Extracted totRNA was treated with an rRNA depletion step using the RiboMinus TM Plant Kit for RNA-Seq (Invitrogen, Waltham, MA, USA) as per the manufacturer's instructions. Ribo depleted totRNA and dsRNA HTS libraries were generated using the Illumina TruSeq Stranded mRNA Library Prep kit, following the manufacturer's protocol, starting after mRNA selection steps [48]. An Illumina NextSeq500 was used to generate single-ended 75 base read files to each sample. Libraries were dual indexed using the IDT for Illumina TruSeq RNA UD Indexes (Illumina, San Diego, CA, USA) and normalized. Then, 24-32 sample libraries were pooled and sequenced using a NextSeq500 high output kit v2.5, 75 cycles (Illumina) which generated between 10 and 16 million reads, on average for each sample. RNAseq files were uploaded to the Sequence Read archive under the bioprojects number PRJNA967701.

Bioinformatics
HTS sample files were imported into Virtool [48] (www.virtool.ca, accessed on 30 September 2021) for sample management, quality control (QC) and data analysis. Reads passing QC using FASTQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/, accessed on 30 September 2021), were mapped to known virus species in plant and bee virus databases updated December 2021 using the Pathoscope 2 pipeline [49]. Reads were aligned to representative isolates of all known plant viruses pulled from Genbank. Viruses with representative isolates receiving at least one mapped read then had reads mapped against all of their known isolates. Bowtie2 2.3.2 [50] was used in local mode for both rounds of read mapping with minimum score (-min-score) set to "L, 20, 1.0", seed length (-L) to 15, and mismatches per seed (-N) set to 0. In the second round, the maximum number of alignments returned (-k) was set to 100. Reads matching viruses were also mapped to a host reference genome. Except for a 0 value for mismatches per seed (-N), default Bowtie2 parameters were used for mapping. Reads were eliminated from the analysis if they had a greater or equal alignment score to the host versus the virus. Multi-mapping reads were handled using a refactored derivative of the Pathoscope2 identification module, which exactly matches the output of the published module. Pathoscope2 makes the assumption that uniquely mapped reads indicate most likely true source genomes. Read values are fractionally reassigned from least likely source genomes to most likely [51]. Virus identification based on Virtool was used to create sample-specific pathogen profiles. For the purposes of this study, a minimum of 10% genome coverage was required for a virus species to be considered a positive detection from both totRNA and dsRNA samples. Sample profiles were combined to create site-specific profiles which included calculating the frequency of detection for each sample type at each site, and the average frequency of detection across all samples, average genome coverage, and viral reads per million (VRPM) for each virus detected across all samples from BC and ON. VRPM is similar to transcripts per million, and was calculated from the total number reads mapping to each individual virus, dividing by genome length of the virus in kilobasepairs (Kbps), and then normalized for the total number of reads in the sequencing run, per million (Supplemental Files S1 and S2). Due to inconsistencies of Apple hammerhead viroid detections, read counts and VRPM were manually annotated using Geneious prime version 11.0.14.1 (Biomatters Inc., Boston, MA, USA). Genbank accession numbers for reference sequences used for virus detection, metadata and calculations can be found in Supplemental Files S1 and S2.

Phylogenetics and Sequence Analysis
Using host genome-subtracted de novo assembled contigs for each sample, sequences were aligned to the BlShV, PDV, and BQCV reference sequence (Supplemental Table S1) using Geneious Prime. Samples with full coat protein sequence coverage were used for pairwise and phylogenetic analysis. Pairwise nucleotide distance comparisons were constructed using Geneious prime. Maximum likelihood phylogenetic trees were constructed using MEGA 11 with 1000 bootstrap replications [52]. Viral coat protein sequences identified in this study were uploaded to Genbank (Supplemental Table S1).
The plant virus profile from the ON site was quite different from BC, with fewer virus species detected. Nine viruses were identified in 12 samples. Two viruses in the genus Nepovirus, tomato ringspot virus (ToRSV) and tobacco ringspot virus (TRSV), were the most frequently detected viruses at 42 and 33%, respectively (Figure 1b; Table 3). Four viruses were detected in all three sites across both provinces: CVA (BC: 44%, ON: 17%), BCCV1 (BC:  Figure 1c). The plant tissue virome was also quite different from BC samples, and even bee-collected samples from ON (Figure 1b; Table 4); four viruses were detected in ON blueberry leaf and flower tissue including blueberry green mosaic-associated virus (BGMaV; genus Vitivirus), BlMaV, BLV, and blueberry virus A (BVA; unassigned Closteroviridae family), which were not observed in bee-collected samples (Figure 1b). BlMaV and BLV were not detected in ON bee samples, whereas in BC both viruses were detected in plant and bee samples.   Bee-infecting viruses were also identified from the same samples (Table 5). In BC samples, 10 viruses were identified with BQCV being detected in every sample tested. Varroa destructor virus 1 (VDV1; genus Iflavirus) and israeli acute paralysis virus (IAPV; genus Aparavirus) were common in BC samples (94% and 63%, respectively), but were absent in the ON samples tested. In ON, only four viruses were identified, with BQCV being detected in 92% of the samples tested (Table 6)  in 92% of the samples tested (  (Figure 2c). Very few bee viruses were detected in plant leaf and flower tissue: Four BC and two ON plant samples tested positive for LSV, and BQCV was also detected in ON plant tissue samples (Table 5; Figure 2a; Supplemental Tables S3 and  S4).

Comparison of Viral Profiles from Bee Pollen, Forager Bees, Hive Bees and Bee Bread
Efficiency of virus detection was evaluated by calculating the average number of viruses identified per sample at all three sites. Each site was treated as an independent replicate, and samples were only compared within each site for statistical tests. Typically, plant virus diversity was significantly higher in bee bread samples relative to other sample types at both BC sites sample (Figure 3a; two-way ANOVA, p < 0.05). While the trend was consistent in samples collected at the ON site, bee bread plant virus species richness was not significantly different from that of pollen or hive bees at this site (Figure 3a). In BC1, bee bread tissues had an average of 6.5 plant viruses detected per sample, while other samples ranged from 2.75-3.25 (Figure 1a; ANOVA p < 0.05). In BC2, bee bread had an average of 9.25 viruses detected per sample, with other sample types ranging from 4.5-6.5 (Figure 1a). Bee bread consistently had the most diverse virome, and significantly more plant viruses were detected in this sample type, although this trend was not as pronounced in ON. When examining the diversity of bee-infecting viruses detected, forager bees and hive bees typically had significantly more viruses detected than other sample types. However, in BC2, a similar number of bee viruses were detected in the pollen sample (Figure 1b).  Error bars indicate standard error. Letters correspond to categories of significant differences between samples based on a two-way ANOVA with Tukey's post hock test (p < 0.05). Each site was analyzed independently, and statistical analysis was performed within the dataset for each site independently.

Viral Diversity in Bee-Collected Samples
A number of plant viruses were detected in relatively high frequency, and with good genome coverage and sequencing depth. Sequences corresponding to the coat protein regions of BlShV, PDV, and BQCV were further examined to better understand their diversity in bee-collected samples. BlShV is only known to infect Vaccinium hosts, while PDV is not known to infect Vaccinium spp. [43]. BlShV was prominent in both BC sites, detected in a total of 9 samples from site 1 and in all 16 samples from site 2, with an average of 57.5% genome coverage (Table 1). Ten samples in total had complete coverage of the BlShV coat protein region, and was used to create a consensus sequence for pair-wise sequence analysis and phylogenetic analysis ( Figure 4A,B). Very little sequence variation was de- Error bars indicate standard error. Letters correspond to categories of significant differences between samples based on a two-way ANOVA with Tukey's post hock test (p < 0.05). Each site was analyzed independently, and statistical analysis was performed within the dataset for each site independently.

Viral Diversity in Bee-Collected Samples
A number of plant viruses were detected in relatively high frequency, and with good genome coverage and sequencing depth. Sequences corresponding to the coat protein regions of BlShV, PDV, and BQCV were further examined to better understand their diversity in bee-collected samples. BlShV is only known to infect Vaccinium hosts, while PDV is not known to infect Vaccinium spp. [43]. BlShV was prominent in both BC sites, detected in a total of 9 samples from site 1 and in all 16 samples from site 2, with an average of 57.5% genome coverage (Table 1). Ten samples in total had complete coverage of the BlShV coat protein region, and was used to create a consensus sequence for pair-wise sequence analysis and phylogenetic analysis ( Figure 4A,B). Very little sequence variation was detected between different samples, with all BlShV CP sequences being 99.9-100% identical to each other in all bee collected samples, and~98.5% identical to the reference sequence ( Figure 4B). All BlShV CP sequences recovered from bee samples clustered together according to phylogenetic analysis, while the cranberry (Vaccinium macrocarpon) isolates formed a separate clade ( Figure 4A). PDV was identified in both BC sites with an average genome coverage of 52.3%, but was not detected in ON samples (Table 1). Complete PDV CP sequences were recovered from 10 different samples, and were more divergent than the BlShV CP sequences ( Figure  5A,B). Complete coverage was mostly obtained from bee bread or hive bee samples, with only one forager bee sample having complete CP coverage (Supplemental Table S2). Following alignment of these sequences with other major PDV CP sequences from three separate phylogroups [53], pairwise distances were calculated and a phylogenetic tree was constructed. Eight samples clustered tightly together and were 97.4-99.2% identical to each other, and also clustered closely with a PDV isolate from a BC cherry tree that is located in phylogroup 2 (AF208741, 98.2-99.5% identity Figure 5A, B). Three other samples also branched within phylogroup 2, with 96-99.2% identity to the BC cherry isolate. BCFV2-BB-B2 sample clustered more closely with two PDV isolates from sweet cherry in the USA (AF208740, GU066792). A final isolate, BCFV1-BB-3B, clustered in phylogroup 1, and was 91.2-92.8% identical to other isolates from this site and 92.3% identical to the BC cherry isolate (Figure 5A,B) [54], suggesting multiple PDV isolates circulating in this system. PDV was identified in both BC sites with an average genome coverage of 52.3%, but was not detected in ON samples (Table 1). Complete PDV CP sequences were recovered from 10 different samples, and were more divergent than the BlShV CP sequences ( Figure 5A,B). Complete coverage was mostly obtained from bee bread or hive bee samples, with only one forager bee sample having complete CP coverage (Supplemental Table S2). Following alignment of these sequences with other major PDV CP sequences from three separate phylogroups [53], pairwise distances were calculated and a phylogenetic tree was constructed. Eight samples clustered tightly together and were 97.4-99.2% identical to each other, and also clustered closely with a PDV isolate from a BC cherry tree that is located in phylogroup 2 (AF208741, 98.2-99.5% identity Figure 5A,B). Three other samples also branched within phylogroup 2, with 96-99.2% identity to the BC cherry isolate. BCFV2-BB-B2 sample clustered more closely with two PDV isolates from sweet cherry in the USA (AF208740, GU066792). A final isolate, BCFV1-BB-3B, clustered in phylogroup 1, and was 91.2-92.8% identical to other isolates from this site and 92.3% identical to the BC cherry isolate ( Figure 5A,B) [54], suggesting multiple PDV isolates circulating in this system.
Of the known bee viruses, BQCV was detected in all bee-collected samples in BC, and in most samples in ON. This virus was also prominent in plant samples, suggesting a ubiquitous presence in these ecosystems (Tables 5 and 6). Complete capsid protein 4 sequences were recovered from 14 samples, predominantly from forager and hive bee samples. BQCV capsid protein 4 sequences were generally highly similar to each other, ranging from 98.4-99.9% identity within bee-collected samples ( Figure 6A). The lone ON isolate branched closely with BC isolates, suggesting low diversity of this virus across Canada ( Figure 6B). Of the known bee viruses, BQCV was detected in all bee-collected samples in BC, and in most samples in ON. This virus was also prominent in plant samples, suggesting a ubiquitous presence in these ecosystems (Tables 5 and 6). Complete capsid protein 4 sequences were recovered from 14 samples, predominantly from forager and hive bee samples. BQCV capsid protein 4 sequences were generally highly similar to each other, ranging from 98.4-99.9% identity within bee-collected samples ( Figure 6A). The lone ON isolate branched closely with BC isolates, suggesting low diversity of this virus across Canada ( Figure 6B).

Discussion
The pollen virome is a largely untapped resource for virus identification and discovery. Recent reports using bee-mediated monitoring to survey for invasive virus species [7,9], and a pollen virome study of wild flowering plants, identified many similar viruses to this study including Ilarviruses ApMV, BlChRSV PNRSV, SNSV, TSV, and nepoviruses ToRSV and TRSV [41]. Bee viruses were also identified in the wild plant pollen virome including DWV [41]. The transmission of viruses through pollen in blueberries, and how this relates to crop and bee health has long been a major research interest [42,43,[55][56][57][58]. Here, we report the first blueberry pollen and bee metavirome, with unique viral profiles based on sample type, individual farms, and geographic regions, demonstrating a powerful area-wide monitoring approach for agricultural pathogens.

Plant Virome
Unsurprisingly many pollen-associated small berry-infecting viruses were detected

Discussion
The pollen virome is a largely untapped resource for virus identification and discovery. Recent reports using bee-mediated monitoring to survey for invasive virus species [7,9], and a pollen virome study of wild flowering plants, identified many similar viruses to this study including Ilarviruses ApMV, BlChRSV PNRSV, SNSV, TSV, and nepoviruses ToRSV and TRSV [41]. Bee viruses were also identified in the wild plant pollen virome including DWV [41]. The transmission of viruses through pollen in blueberries, and how this relates to crop and bee health has long been a major research interest [42,43,[55][56][57][58]. Here, we report the first blueberry pollen and bee metavirome, with unique viral profiles based on sample type, individual farms, and geographic regions, demonstrating a powerful area-wide monitoring approach for agricultural pathogens.

Plant Virome
Unsurprisingly many pollen-associated small berry-infecting viruses were detected in this study. BlShV was the most prominent virus identified in BC, consistent with previous reports of this virus being prevalent in the Pacific Northwest and readily pollen transmitted [44]. Although not known to be pollen transmitted, BlScV is a serious pathogen in blueberry production systems and is primarily vectored by aphids [56]. BlScV was only detected twice in forager bee-collected samples, but was prominent in plant samples. Association of this virus with forager bees could be due to direct flower tissue contact as opposed to associations with pollen [56]. BlShV, BlScV, BlMaV and BLV were all detected in both bee and plant tissues, suggesting a good overlap between bee-related identification of plant viruses and viruses infecting the host of interest.
While there was good overlap between viruses detected in bee and plant samples at both BC sites, there was no such overlap between the viruses detected from samples collected from ON. ToRSV and TRSV were common in bee samples but absent in plants, while BLV, BVA, BlMaV, and BGMaV were identified in plant tissues, but not bee-collected samples. ToRSV and TRSV could have been missed in our random sampling of composite leaf and flower tissues, or they could have originated from other host plants [45,59]. ToRSV can infect a wide host range, and cause serious economic losses [60]. BGMaV is a member of the Vitivirus genus, while BVA is a member of the Closterovirus genus, and neither genus is commonly associated with pollen transmission [61,62]. BLV and BlMaV were both detected in bee-collected and plant samples in BC but were not detected in ON bee samples, only in plant samples. The number of reads mapping to BlMaV was much lower in ON plant tissues compared to BC, which could reflect differences in the levels of infection at each site. Further research is required to better understand the specificity and sensitivity of plant virus detection through bees relative to the proportion of plants infected by these viruses in the immediate area.
Ilarviruses, in particular, were very prominent and widely detected. Many identified ilarviruses are commonly associated with small berry production such as SNSV or BCRV, and are likely infecting hosts aside from blueberries, such as Rubus spp., in the immediate area. BCRV and GaIV have not previously been reported in Canada, further demonstrating the advantages of bee monitoring in detecting new and emerging pathogens [63][64][65]. However, BCRV and GalV were identified with low read counts and genome coverage, further research should confirm and validate their presence in these systems. Other identified ilarviruses are more commonly associated with tree fruits such as PNRSV and PDV, and have not been reported to directly infect Vaccinium plants, to the best of our knowledge. It is unlikely these sequences were derived from blueberry plants, and are more likely to be present in other plant species nearby. Curiously, there was very little sequence variation for recovered BlShV CP sequences, which could reflect a founder effect due to clonal propagation practices in blueberry production. Other ilarviruses have been reported to have low CP sequence diversity including American plum line pattern virus [53,66]. Very few BlShV sequences are available in GenBank, and little is known regarding the population diversity of BlShV. Other ilarviruses such as PDV did show more diverse sequences. RNA3 of PDV separates into three phylogroups [53], we identified PDV CP sequences corresponding to all three phylogroups, suggesting diverse origins of PDV and potentially multiple strains of this virus circulating in the environment.
We expected the greatest diversity of plant viruses to be detected in plant tissues (pollen, bee bread, and plant tissue), and more bee viruses to be identified in bee samples (forager bee, hive bee), and for the most part this was consistent. Bee bread displayed the highest diversity of plant viruses, and also the highest sequence diversity of individual viruses. This is likely due to bee bread being collected over the seasonal flowering period of different plant species. Many of the viruses identified in bee bread are often associated with apple or tree fruit production systems such as ASGV, citrus concave gum-associated virus (genus Coguvirus), cherry virus F (genus Fabavirus) and citrus virus A (genus Coguvirus) that may have originated from plants flowering concurrently, before or after blueberries. Bee bread represents a highly diverse collection of plant pollen samples, and could be useful for further study of plant virus diversity. Surprisingly, hive bees also showed an unexpected richness in plant virus species compared to the other sample types which in theory do not directly interact with plants, suggesting the potential for dissemination of plant viruses within a bee colony. BlScV and BlMaV were only identified in forager or hive bee samples, and could have been acquired through direct flower contacts, as opposed to pollen transmission.

Bee Virome
Bee health and monitoring the viral load of bee colonies are other important aspects of ecosystem health monitoring [9,67,68]. BQCV and LSV were common in BC and ON sites, while VDV1 and IAPV were frequently detected at both BC sites and absent in ON. BQCV and IAPV have previously been detected in pollen and pollen pellets, and can be foodborne transmitted [69,70]. LSV, DWV, and SBV have been found in pollen and pollen pellets [71][72][73]. Peculiarly, LSV was widely prevalent in bee samples but absent in bee bread, and yet was also identified in plant leaf and flower tissue. Fecal-oral routes have been suggested to play an important role in LSV dissemination, and its relation to plant food source is important [74]. One dsDNA virus, apis mellifera filamentous virus, was detected in both BC and ON. Surprisingly, multiple samples had over 60% genome coverage for this virus, demonstrating the ability to identify DNA viruses from totRNA extracted samples.

Metagenomics-based Monitoring of Plant Viruses
A large number of plant and bee viruses were detected widely in different sample types, sites, and regions. An arbitrary genome coverage threshold of 10% was established to remove low-confidence identifications. Other genomics-based studies have argued for genome coverage cut off levels between 10 and 15% for genomics, and 20% for metagenomics-based surveys [2,41,51]. Genome coverage threshold cut offs can be arbitrary; yet, with metagenomics studies, it is important to cast a wide net to understand the diversity of viruses within the local population, and could help with the identification of emerging viruses. The trade-off to using a lower threshold is a greater chance of false positives inherent to the method [75]. Other studies have employed similar biovigilance approaches through RT-PCR-based assays to detect major viruses of concern in Prunus production systems [76]. Metagenomics-based approaches can complement more focused assays through identification of all pollen-associated viruses, and provide valuable viral sequence diversity information.
As evidenced from the wide diversity of plant viruses identified, pathogen monitoring through bee-mediated pollination activities is an excellent approach to understand the diversity of viruses present in an ecosystem. Previous studies of blueberry viruses have identified~12 high priority viruses in the US and BC [44]. Many of the viruses listed in this report were identified through bee/pollen screening including BlShV, BlScV, BLV, BlMaV, TRSV and ToRSV. Other viruses identified in plant tissue but not bee samples include BGMaV and BVA, which have not been previously reported in Canada. While many viruses are detectable in this system, the exact risks of transmission for each specific virus can be unclear. Understanding pollen transmission risks for each virus must still be individually determined. Bee-mediated virus monitoring can perform as efficiently or better than randomly collected samples while also providing useful data regarding sequence diversity, demonstrating the utility of a bee-mediated metagenomics-based pathogen monitoring approach in agricultural systems. Despite these advantages, it is important to understand the limitations of bee-mediated pathogen monitoring: only pollen-transmitted viruses or those contaminating bees acquired through interactions with flowers can be detected.

Supplementary Materials:
The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/v15051209/s1. Table S1: Blueberry shock virus, prune dwarf virus, and black queen cell virus coat and capsid protein nucleotide sequence genbank accession numbers and number of mapped reads; Table S2: Sequence read archive accession numbers for RNAseq files; Table S3: Bee viruses detected in leaf and flower samples from two BC blueberry farms; Table S4: Bee viruses detected in leaf and flower samples from on ON blueberry farm; Supplemental File S1; Supplemental File S2.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.