Surveillance of Viruses in Varroa destructor Samples Collected from Honey Bee Colonies in Ontario, Canada, between 2015 and 2019

Varroa destructor is the main cause of honey bee colony mortality in the northern hemisphere (Guzman-Novoa et al., 2010; Le Conte et al., 2010). It acts as a mechanical and biological vector for viruses, including Deformed wing virus (DWV) and Israeli acute bee paralysis virus (IAPV) (Di Prisco et al., 2011; Gisder & Genersch, 2020). Viral infections have been linked to honey bee colony losses (McMenamin & Genersch, 2015); some cause severe damage to individual bees, like deformities in developing bees infected with DWV-A (de Miranda & Generesch, 2010) and learning impairment in adult bees (Iqbal & Mueller, 2007).

The surveillance of viruses in V. destructor samples is essential to assess risks of emerging virulent variants and evaluate their impact on honey bee health. The Ontario Ministry of Agriculture, Food, and Rural Affairs (OMAFRA) conducted a monitoring project from 2015 to 2019, that included more than 1,000 honey bee colonies located in five different regions of the province; the colonies were monitored for the recording of their strength and V. destructor levels in the spring, summer and fall (OMAFRA, 2016). Overwinter colony mortality was also assessed, and samples of live bees were collected to identify viruses with real time PCR (qPCR), including DWV, Black queen cell virus (BQCV), IAPV, Acute bee paralysis virus (ABPV), Chronic bee paralysis virus (CBPV), Kashmir bee virus (KBV) and Sacbrood virus (SBV). The aim of this study was to determine the presence of pathogens in V. destructor samples with the use of a metagenomic technique. We determined the load and species of viruses in V. destructor mites retrieved from honey bee samples during the OMAFRA monitoring project (OMAFRA, 2016). Our study showed that metagenomic tools assist in active pathogen surveillance in V. destructor populations, which will help understand the effect of varroosis on bee health and prevent extreme colony mortality.

Ten V. destructor mites were randomly selected from each of five honey bee samples collected during the 2015–2019 OMAFRA monitoring project (OMAFRA, 2016) and kept at −80°C. Thus, a total of fifty mites were used for molecular analyses. The molecular analyses of viral pathogens were performed at the Animal Health Laboratory, University of Guelph (Ontario, Canada). The RNA from the fifty individual V. destructor mites was extracted using a SPEX Mini-G homogenizer and MagMAX™ 96 (Applied Biosystems; Mississauga, ON, Canada), eluting 100 μL of nucleic acid template per mite. cDNA was generated with the SuperScript™ III First-Strand Synthesis SuperMix (Invitrogen; Mississauga, ON, Canada), with the use of proprietary random hexamers as primers. For the RNAseq, 1 ng of cDNA was used to construct cDNA libraries with Nextera XT Library Preparation Kit (Illumina; San Diego, CA, USA). RNAseq was performed as 300 bp paired end reads with the use of Illumina Miseq (San Diego, CA, USA). The quality of the raw sequencing reads was assessed with FastQC (Wingett & Andrews, 2010). Taxonomic labels were assigned to metagenomic cDNA sequences with Kraken v2 (Wood & Salzberg, 2014) with a confidence score of 0.65. FastVirome pipeline (Tithi et al., 2018) was used to identify and quantify viral transcript abundance in the samples, with the use of a precomputed Kallisto index (Bray et al., 2016). Edge R in Degust was used to estimate transcript level abundance among known negative samples (control) and twenty-seven mite samples with the presence of at least one virus (Robinson et al., 2010; Powell, 2019; p<0.05).

The RNA of honey bees (pool of ten whole bodies, from the same samples used to retrieve V. destructor mites) and V. destructor mites was extracted and used to calculate the number of viral genome copies (gc) with the use of a LC480 thermocycler (Roche; Mississauga, ON, Canada). Each reaction consisted of 2 μl template, 1 μl forward and reverse primers (500 nM), 10 μl LC480 SYBR mastermix (1×) and 7 μL nuclease-free H₂O. The target gene copy number was calculated as per Bustin et al. (2009) with specific primers for ABPV, BQCV and DWV-A (Van Engelsdorp et al., 2009); CBPV and KBV (Boncristiani et al., 2012); IAPV (this study; F5′CAGTTACAGCAGTCGTATGG and R5′TTGT-GTACGTTGGGAGTATTG) and SBV (this study; F’ACCGATTTGTTTAATGGTTGGG and R5′AT-TCCAGATTCTTCGTCCACTC). The DWV-A and ABPV genome copy numbers in mites and the corresponding bee samples were subjected to a Spearman's correlation test, as the data was not normal based on Shapiro Wilk test (α of 0.05).

Kraken was used to identify eight viruses relevant to honey bee health in thirty-five out of the fifty V. destructor samples, including DWV-A, VDV-1, and Bee macula virus (Tab. 1 and Suppl. Tab. 1). Variability was observed in the number of viruses present in the samples and the number of reads of each virus, which indicated a complex distribution of viruses in the tested samples. Additionally, RNA sequences of DWV-A (NC_004830.2), VDV-1 (NC_006494.1), VDV-2 (NC_040601.1), VDV-3 (NC_040313.1), ABPV (NC_002548.1), Apis filamentous virus (NC_027925.1), KBV (NC_004807.1), Bee macula virus (NC_027631.1), and transcripts of Apis rhabdovirus were identified with the use Kallisto. Only DWV-A (NC_004830.2) and VDV-1 (NC_006494.1) were higher compared to samples with zero viral abundance (control) (7.5 and 5.72 log2 fold, respectively; p<0.05). These results in line with previous studies show that DWV-A and VDV-1 are among the most abundant viruses in honey bee populations (Ryabov et al., 2017). In addition, a positive correlation between the genome copy numbers of DWV-A in mites and the corresponding bee samples (r=0.60, p=0.05) but not for ABPV (r=−0.059, p=0.95) confirms a strong relationship between DWV-A, V. destructor and their host (Bowen-Walker & Gunn, 1999). Only one virus identified in honey bee samples with qPCR was not present in mites (CBPV; Tab. 2). Additionally, the presence of Bee macula virus and Apis rhabdovirus in V. destructor samples indicates that further studies on their transmission and pathogenicity are needed. As previously proposed, high throughput techniques assist in the identification of emerging viruses in honey bee populations (McMenamin & Flenniken, 2018; Abou-Shaara, 2019). Our study shows that RNAseq combined with Kraken and Kallisto, can be used in disease monitoring programs to identify viral pathogens in V. destructor populations. This study identified virulent variants of DWV (VDV-1, VDV-2 and VDV-3; Levin et al., 2016; Natsopoulou et al., 2017), but no positive samples of Slow bee paralysis virus or Lake Sinai virus were found, nor bacterial or fungal diseases relevant to bee health were identified in the metagenomic analysis. The diagnostic tools used in this study can aid in the active surveillance of novel viral variants in V. destructor populations and help prevent extreme honey bee colony losses linked to viral diseases.