TY - JOUR
AU - Schroeder, Declan C.
AB - Introduction Honey bees are important pollinators of flowering plants and various managed crops [1–3]. Over the past two decades, a large number of reports from different parts of the world show the continued decline in the health of western honey bee, Apis mellifera, populations [4–7]. Multiple biotic and abiotic factors such as pests, parasites and pathogens, pesticides, habitat alteration, poor nutrition and lack of genetic diversity were found to contribute to poor colony health [8–10]. Two specific pests and pathogens of honey bees are the parasitic mite, Varroa destructor, and the virus, Deformed wing virus (DWV) [11, 12]. The majority of viruses found in the honey bee virome have a single-stranded positive-sense RNA genome (+ssRNA) [13]. In particular, DWV master variants A & B (DWV-A & DWV-B), Acute bee paralysis virus (ABPV), and Israeli acute paralysis virus (IAPV) have been shown to have a major impact on colony survival [14–17]. DWV predominates in many honey bee viromes [13, 18–21]. Pests and pathogens move between colonies within and between apiaries [22–24]. The majority of honey bee viruses cause asymptomatic infections in honey bees [13]; consequently, timely testing for viruses is critical to monitor a colony’s ’ health status [17]. Techniques developed for detecting viruses in honey bees are based on various approaches including enzyme-linked immunosorbent assay (ELISA), oligonucleotide microarray, cell lines, quantitative PCR (qPCR and RT-qPCR) and metagenomic next generation sequencing (mNGS) [25–32]. PCR remains one of the most favored and applied techniques for virus detection and disease diagnosis. It is routinely used for virus diagnostics in many research fields and is a standard for evaluating viral loads in biomedical research, agricultural and environmental sectors [33–38]. NGS technology is an additive and sometimes alternative method, capable of identifying unknown viruses in the sample, including variants and quasispecies [39–41], and is a reliable tool for validation of PCR results [32, 39, 40, 41]. To perform virus screening with qPCR, nucleic acid (RNA or DNA) extraction is often the first step in the screening process. Nucleic acid isolation is initiated with mechanical or chemical disruption of cells, following a purification step by precipitation with ethanol/isopropanol, affinity purification columns or magnetic beads-based technology [33–38, 42, 43]. A common approach in bee research is based on using manual methods of extraction utilizing organic solvents or affinity purification columns kits, or a combination of both. These methods are reliable, but time-consuming [21, 43–48]. Acid phenol RNA extraction based on using organic solvents such as phenol and chloroform, or commercially available TRIzol®, involves using reagents which are highly volatile and toxic [25] and require special and, hence, expensive disposal. In addition, an acid-phenol phase separation method for isolating RNA from DNA requires many steps including multiple steps of incubation, vortexing, centrifugation, rescuing aqueous phase, storing samples at -20C overnight, and washing. As a result, nucleic acid extraction is often a labor intensive and lengthy process [43–45]. Although the affinity column purification method has fewer steps compared to the acid phenol RNA extraction, being based on the manual approach, it is limited in its capacity to provide a rapid solution for large-scale extractions for surveys and diagnostic purposes [21, 46–49]. Moreover, in contrast to the automated approach, manual methods bear a risk of cross-contamination and human-related bias. The magnetic beads-based technology has been widely used in biomedical research [33–37], and has become one of the most used methods to extract viral RNA for screening for SARS-CoV-2 [36], as rapid diagnosis of COVID-19 is essential to restrict its spreading. A same extraction method combined with an automated approach, which is rapid, high-throughput, uniform and bears low cross contamination risk [50–52] may be employed for honey bee research, in particular for virus screening surveys to screen multiple colonies in a limited timeframe. The magnetic beads-based technology has previously been used for screening bees for viruses, either for the detection of viruses in individual [53, 54] or pooled samples [53, 55, 56]. However, to our knowledge this method has yet to be compared or validated against the column-based extraction method. Moreover, in all these studies sample preparation remained manual, requiring grinding bees’ samples by hand, and using additional lysis solutions for dissociation of tissues. Here we report on developing a semi-automated high-throughput approach for the detection of honey bee viruses which can be scaled for the simultaneous extraction of 96 pooled bee samples at a time. This method is based on automated dissociation of bee samples in a phosphate buffered saline solution, vacuum manifold-based sterile filtration and a robotic extraction of viral nucleic acids using magnetic beads-based technology. To quantify DWV loads we employ RT-qPCR and confirm our findings by using the Oxford Nanopore Technologies (ONT) GridION sequencing platform. Materials and methods Samples collection and storage Samples of 50–100 adult honey bees, A. mellifera were collected from the colonies located at Agricultural Experimental station in Rosemount, Minnesota, US, during September 2020 and April 2021 (Table 1). Colonies were treated to control Varroa: with formic acid pads (Formic Pro®, Mann Lake Ltd) in late August, and oxalic acid dribble in late October. Bee samples were stored at -80°C immediately after washing with 70% Ethanol. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 1. Forty-eight honey bees pooled samples collected from 20 colonies collected between September 2020 and April 2021. https://doi.org/10.1371/journal.pone.0297623.t001 Sample processing Homogenization. Pooled samples of 30 bees representing a single colony (n = 24) were transferred to a 50 ml sterile gentleMACS™ M Tube (Miltenyi Biotec Inc. Auburn, CA, USA), 15 ml of sterile 1X PBS was added and the samples were dissociated using the gentleMACS Dissociator (V1.02, Miltenyi Biotec Inc. Auburn, CA, USA, RNA_02.01). Following centrifugation at 4,700x g, for 5 minutes, at room temperature (RT), 2 ml of the homogenate were transferred to a sterile 2 ml centrifuge tube, and centrifuged at 21,100x g, for 5 minutes, at RT. A 200 μl aliquot of the homogenate was transferred to a sterile 1.5 ml Eppendorf tube for subsequent analysis. All samples were stored at -80°C prior to proceeding to the next step. Filtration. To exclude particles larger than 0.45 um, a filtration step was carried out. 150 μl of previously obtained homogenate from each of the 24 samples were applied to a 0.45 um sterile filtration plate (MultiScreen, MilliporeSigma, USA) and went through a filtration using a Vacuum Manifold (MultiScreenHTS, MilliporeSigma, USA) for a simultaneous filtration of the 24 samples. Viral RNA extraction. To isolate viral RNA, two different approaches were taken: manual isolation with a NucleoSpin Virus kit (Macherey-Nagel, Düren, Germany); and automated extraction using the NucleoMag Virus kit (Macherey-Nagel, Düren, Germany) and a Magnetic Particle Processor (MPP) (KingFisher Flex, Thermo Fisher Scientific, USA). Manufacturer’s instructions were followed for both protocols. For RNA isolation with NucleoSpin kit, 100 ul of sterile 1X PBS was added to 50 ul bee filtrate for a total volume of 150 ul. For NucleoMag kit, 150 ul 1XPBS was added to 50 ul filtrate for a total volume of 200 ul. RNA purity and concentration were assessed with the NanoDrop Spectrophotometer (NanoDrop™, Thermo Fisher Scientific, USA). RT-qPCR. To screen for DWV-A and DWV-B viral loads, Power-Up™ SYBR® Green RNA-to-Ct 1-Step Kit (Applied Biosystems, Foster City, CA, USA) was used. Each reaction was performed in duplicate employing a Bio-Rad real-time PCR machine (CFX96, Bio-Rad, USA) following the ABC assay protocol as described earlier [57]. To quantify viral loads in the RNA samples 5 μl of SYBR mix, 3.92 μl of RNA diluted in molecular grade water to 50 ng/μl, 0.08 μl reverse transcriptase, 0.5 μl (10 pmol) reverse primer (DWV-A or B) and 0.5 μl (10 pmol) universal forward primer were used. Reverse transcription occurred at 45°C for 10 min and denaturation occurred at 95°C for 10 min, followed by 35 cycles of denaturation at 95°C for 15 s, annealing at 58°C for 15 s, and extension at 72°C for 15 s. A high-resolution melt analysis was performed between 72°C and 90°C, at 0.1°C increments, each with a 5 s hold period. Viral genome copies were calculated as described previously [57] and expressed as log10 viral RNA copies per sample. Library preparation and sequencing. To synthesize double stranded cDNA, 10 ul of viral RNA, 1 ul of N6 Primer II A (24 uM, TakaraBio, USA), 1 ul of SMARTer IIA Oligo (24 uM, TakaraBio, USA), 1 ul of 10x Template Switching Reverse Transcriptase (New England Biolabs, MA, USA) for the synthesis of the 1st strand, and 1 ul of Primer IIA (12 uM, TakaraBio, USA) and 1 ul of PrimeSTAR GXL polymerase (TakaraBio, USA) for the synthesis of the 2nd strand were used following manufacturer user manual instructions (TakaraBio, USA; New England Biolabs, MA, USA). The synthesized cDNA was purified using SPRI AMPure beads according to manufacturer’s instructions (TakaraBio, USA), and its integrity and quality were assessed with Qubit 4 Fluorometer (Qubit™4 Fluorometer, Thermo Fisher Scientific, USA) in accordance with One time dsDNA HS assay kit user manual. To prepare the library for nanopore sequencing, the ONT Rapid Barcoding Sequencing Kits (SQK-RPB004, Oxford Nanopore Technologies, UK) were used as per the manufacturer’s guidelines. Libraries were pooled on FLO-MIN106 flow cells and run on the GridION. Sequencing performance was monitored and was terminated after 24 h. Metagenomic analysis. Sequencing reads were filtered to a minimum length (≥200 bp) and Q-value (≥9) by MinKNOW v4.3.4. Basecalling and demultiplexing was performed using Guppy v6.4 with the high accuracy model. Guppy is only available to NanoPore customers through their community site (https://community.nanoporetech.com). PoreChop v0.2.4 [58] was used to remove the nanopore barcode adapter sequences. To assemble metagenome-assembled contigs, the quality filtered reads were assembled by Canu v2.2 [59] using the following assembly parameters: -nanopore maxInputCoverage = 2000 corOutCoverage = all corMinCoverage = 0 corMhapSensitivity = high minoverlap = 50 minread = 200 genomesize = 5000. Of 48 viromes, 42 were successfully assembled to generate contigs (24 NucleoSpin® and 18 NucleoMag™). Contigs with a minimum length of 1 kbp were binned manually with the anvi’o v7.1 [60, 61] interactive interface. Briefly, anvi’o profiled the contigs using Prodigal v2.6.3 [62], with default parameters, then reads were mapped to the contig database using Minimap2 v2.24 [63], and the read recruitment was stored as a BAM file using samtools. Anvi’o profiles each BAM file, estimating the read coverage and detection statistics of each contig, we then normalized coverage as Reads per kilobase of transcript per Million reads mapped (RPKM) using the python package bioinforkit v2.1.0 [64]. Coverage and RPKM was combined into a merged profile database. The contigs database was populated with additional data, incorporating HMMER results against Virus Orthologous Groups (VOGs; https://vogdb.org/) in addition to the standard anvi’o HMMR profiles, NCBI COGs and KEGG KOfam database [65]. Contig taxonomy was predicted by running Kraken2 v2.1.2 [66, 67] using the non-redundant NCBI database on the gene calls. Finally, merged profiles were clustered with the automatic binning algorithm CONCOCT, and the anvi’o profile was visualized for manual binning. Binning was guided by sequence composition similarity (visualized as a dendrogram in the Anvi’o interface), and the presence of viral HMM hits to the VOG database. Validation of honey bee RNA virus genomes from binned contigs. Binned contigs were size filtered to 5 kbp, then aligned to reference genomes using Minimap2. References for complete genomes included DWV-A (NC_004830.2) and DWV-B/VDV-1 (NC_006494.1). Alignments were used to identify probable viral genomes. Prodigal was utilized to identify coding regions and regions were annotated using BLASTn (Nucleotide Basic Local Alignment Search Tool). Data analysis and visualization. Data was visualized and statistically analyzed using R v4.1.0 in RStudio build 576 and Microsoft Excel software. Statistical tests were performed using base features in R and data visualized with the package ggplot2 v3.4.1. Numeric values of the read length, number of reads and average quality of each read was acquired using SeqKit v2.3.0 [68]. Welch approximation t-test was used to compare data generation outputs from NucleoSpin® versus NucleoMag™ (i.e., RNA yields, number of reads generated, read length, and average read quality, average read mapped to reference genomes). Samples collection and storage Samples of 50–100 adult honey bees, A. mellifera were collected from the colonies located at Agricultural Experimental station in Rosemount, Minnesota, US, during September 2020 and April 2021 (Table 1). Colonies were treated to control Varroa: with formic acid pads (Formic Pro®, Mann Lake Ltd) in late August, and oxalic acid dribble in late October. Bee samples were stored at -80°C immediately after washing with 70% Ethanol. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 1. Forty-eight honey bees pooled samples collected from 20 colonies collected between September 2020 and April 2021. https://doi.org/10.1371/journal.pone.0297623.t001 Sample processing Homogenization. Pooled samples of 30 bees representing a single colony (n = 24) were transferred to a 50 ml sterile gentleMACS™ M Tube (Miltenyi Biotec Inc. Auburn, CA, USA), 15 ml of sterile 1X PBS was added and the samples were dissociated using the gentleMACS Dissociator (V1.02, Miltenyi Biotec Inc. Auburn, CA, USA, RNA_02.01). Following centrifugation at 4,700x g, for 5 minutes, at room temperature (RT), 2 ml of the homogenate were transferred to a sterile 2 ml centrifuge tube, and centrifuged at 21,100x g, for 5 minutes, at RT. A 200 μl aliquot of the homogenate was transferred to a sterile 1.5 ml Eppendorf tube for subsequent analysis. All samples were stored at -80°C prior to proceeding to the next step. Filtration. To exclude particles larger than 0.45 um, a filtration step was carried out. 150 μl of previously obtained homogenate from each of the 24 samples were applied to a 0.45 um sterile filtration plate (MultiScreen, MilliporeSigma, USA) and went through a filtration using a Vacuum Manifold (MultiScreenHTS, MilliporeSigma, USA) for a simultaneous filtration of the 24 samples. Viral RNA extraction. To isolate viral RNA, two different approaches were taken: manual isolation with a NucleoSpin Virus kit (Macherey-Nagel, Düren, Germany); and automated extraction using the NucleoMag Virus kit (Macherey-Nagel, Düren, Germany) and a Magnetic Particle Processor (MPP) (KingFisher Flex, Thermo Fisher Scientific, USA). Manufacturer’s instructions were followed for both protocols. For RNA isolation with NucleoSpin kit, 100 ul of sterile 1X PBS was added to 50 ul bee filtrate for a total volume of 150 ul. For NucleoMag kit, 150 ul 1XPBS was added to 50 ul filtrate for a total volume of 200 ul. RNA purity and concentration were assessed with the NanoDrop Spectrophotometer (NanoDrop™, Thermo Fisher Scientific, USA). RT-qPCR. To screen for DWV-A and DWV-B viral loads, Power-Up™ SYBR® Green RNA-to-Ct 1-Step Kit (Applied Biosystems, Foster City, CA, USA) was used. Each reaction was performed in duplicate employing a Bio-Rad real-time PCR machine (CFX96, Bio-Rad, USA) following the ABC assay protocol as described earlier [57]. To quantify viral loads in the RNA samples 5 μl of SYBR mix, 3.92 μl of RNA diluted in molecular grade water to 50 ng/μl, 0.08 μl reverse transcriptase, 0.5 μl (10 pmol) reverse primer (DWV-A or B) and 0.5 μl (10 pmol) universal forward primer were used. Reverse transcription occurred at 45°C for 10 min and denaturation occurred at 95°C for 10 min, followed by 35 cycles of denaturation at 95°C for 15 s, annealing at 58°C for 15 s, and extension at 72°C for 15 s. A high-resolution melt analysis was performed between 72°C and 90°C, at 0.1°C increments, each with a 5 s hold period. Viral genome copies were calculated as described previously [57] and expressed as log10 viral RNA copies per sample. Library preparation and sequencing. To synthesize double stranded cDNA, 10 ul of viral RNA, 1 ul of N6 Primer II A (24 uM, TakaraBio, USA), 1 ul of SMARTer IIA Oligo (24 uM, TakaraBio, USA), 1 ul of 10x Template Switching Reverse Transcriptase (New England Biolabs, MA, USA) for the synthesis of the 1st strand, and 1 ul of Primer IIA (12 uM, TakaraBio, USA) and 1 ul of PrimeSTAR GXL polymerase (TakaraBio, USA) for the synthesis of the 2nd strand were used following manufacturer user manual instructions (TakaraBio, USA; New England Biolabs, MA, USA). The synthesized cDNA was purified using SPRI AMPure beads according to manufacturer’s instructions (TakaraBio, USA), and its integrity and quality were assessed with Qubit 4 Fluorometer (Qubit™4 Fluorometer, Thermo Fisher Scientific, USA) in accordance with One time dsDNA HS assay kit user manual. To prepare the library for nanopore sequencing, the ONT Rapid Barcoding Sequencing Kits (SQK-RPB004, Oxford Nanopore Technologies, UK) were used as per the manufacturer’s guidelines. Libraries were pooled on FLO-MIN106 flow cells and run on the GridION. Sequencing performance was monitored and was terminated after 24 h. Metagenomic analysis. Sequencing reads were filtered to a minimum length (≥200 bp) and Q-value (≥9) by MinKNOW v4.3.4. Basecalling and demultiplexing was performed using Guppy v6.4 with the high accuracy model. Guppy is only available to NanoPore customers through their community site (https://community.nanoporetech.com). PoreChop v0.2.4 [58] was used to remove the nanopore barcode adapter sequences. To assemble metagenome-assembled contigs, the quality filtered reads were assembled by Canu v2.2 [59] using the following assembly parameters: -nanopore maxInputCoverage = 2000 corOutCoverage = all corMinCoverage = 0 corMhapSensitivity = high minoverlap = 50 minread = 200 genomesize = 5000. Of 48 viromes, 42 were successfully assembled to generate contigs (24 NucleoSpin® and 18 NucleoMag™). Contigs with a minimum length of 1 kbp were binned manually with the anvi’o v7.1 [60, 61] interactive interface. Briefly, anvi’o profiled the contigs using Prodigal v2.6.3 [62], with default parameters, then reads were mapped to the contig database using Minimap2 v2.24 [63], and the read recruitment was stored as a BAM file using samtools. Anvi’o profiles each BAM file, estimating the read coverage and detection statistics of each contig, we then normalized coverage as Reads per kilobase of transcript per Million reads mapped (RPKM) using the python package bioinforkit v2.1.0 [64]. Coverage and RPKM was combined into a merged profile database. The contigs database was populated with additional data, incorporating HMMER results against Virus Orthologous Groups (VOGs; https://vogdb.org/) in addition to the standard anvi’o HMMR profiles, NCBI COGs and KEGG KOfam database [65]. Contig taxonomy was predicted by running Kraken2 v2.1.2 [66, 67] using the non-redundant NCBI database on the gene calls. Finally, merged profiles were clustered with the automatic binning algorithm CONCOCT, and the anvi’o profile was visualized for manual binning. Binning was guided by sequence composition similarity (visualized as a dendrogram in the Anvi’o interface), and the presence of viral HMM hits to the VOG database. Validation of honey bee RNA virus genomes from binned contigs. Binned contigs were size filtered to 5 kbp, then aligned to reference genomes using Minimap2. References for complete genomes included DWV-A (NC_004830.2) and DWV-B/VDV-1 (NC_006494.1). Alignments were used to identify probable viral genomes. Prodigal was utilized to identify coding regions and regions were annotated using BLASTn (Nucleotide Basic Local Alignment Search Tool). Data analysis and visualization. Data was visualized and statistically analyzed using R v4.1.0 in RStudio build 576 and Microsoft Excel software. Statistical tests were performed using base features in R and data visualized with the package ggplot2 v3.4.1. Numeric values of the read length, number of reads and average quality of each read was acquired using SeqKit v2.3.0 [68]. Welch approximation t-test was used to compare data generation outputs from NucleoSpin® versus NucleoMag™ (i.e., RNA yields, number of reads generated, read length, and average read quality, average read mapped to reference genomes). Homogenization. Pooled samples of 30 bees representing a single colony (n = 24) were transferred to a 50 ml sterile gentleMACS™ M Tube (Miltenyi Biotec Inc. Auburn, CA, USA), 15 ml of sterile 1X PBS was added and the samples were dissociated using the gentleMACS Dissociator (V1.02, Miltenyi Biotec Inc. Auburn, CA, USA, RNA_02.01). Following centrifugation at 4,700x g, for 5 minutes, at room temperature (RT), 2 ml of the homogenate were transferred to a sterile 2 ml centrifuge tube, and centrifuged at 21,100x g, for 5 minutes, at RT. A 200 μl aliquot of the homogenate was transferred to a sterile 1.5 ml Eppendorf tube for subsequent analysis. All samples were stored at -80°C prior to proceeding to the next step. Filtration. To exclude particles larger than 0.45 um, a filtration step was carried out. 150 μl of previously obtained homogenate from each of the 24 samples were applied to a 0.45 um sterile filtration plate (MultiScreen, MilliporeSigma, USA) and went through a filtration using a Vacuum Manifold (MultiScreenHTS, MilliporeSigma, USA) for a simultaneous filtration of the 24 samples. Viral RNA extraction. To isolate viral RNA, two different approaches were taken: manual isolation with a NucleoSpin Virus kit (Macherey-Nagel, Düren, Germany); and automated extraction using the NucleoMag Virus kit (Macherey-Nagel, Düren, Germany) and a Magnetic Particle Processor (MPP) (KingFisher Flex, Thermo Fisher Scientific, USA). Manufacturer’s instructions were followed for both protocols. For RNA isolation with NucleoSpin kit, 100 ul of sterile 1X PBS was added to 50 ul bee filtrate for a total volume of 150 ul. For NucleoMag kit, 150 ul 1XPBS was added to 50 ul filtrate for a total volume of 200 ul. RNA purity and concentration were assessed with the NanoDrop Spectrophotometer (NanoDrop™, Thermo Fisher Scientific, USA). RT-qPCR. To screen for DWV-A and DWV-B viral loads, Power-Up™ SYBR® Green RNA-to-Ct 1-Step Kit (Applied Biosystems, Foster City, CA, USA) was used. Each reaction was performed in duplicate employing a Bio-Rad real-time PCR machine (CFX96, Bio-Rad, USA) following the ABC assay protocol as described earlier [57]. To quantify viral loads in the RNA samples 5 μl of SYBR mix, 3.92 μl of RNA diluted in molecular grade water to 50 ng/μl, 0.08 μl reverse transcriptase, 0.5 μl (10 pmol) reverse primer (DWV-A or B) and 0.5 μl (10 pmol) universal forward primer were used. Reverse transcription occurred at 45°C for 10 min and denaturation occurred at 95°C for 10 min, followed by 35 cycles of denaturation at 95°C for 15 s, annealing at 58°C for 15 s, and extension at 72°C for 15 s. A high-resolution melt analysis was performed between 72°C and 90°C, at 0.1°C increments, each with a 5 s hold period. Viral genome copies were calculated as described previously [57] and expressed as log10 viral RNA copies per sample. Library preparation and sequencing. To synthesize double stranded cDNA, 10 ul of viral RNA, 1 ul of N6 Primer II A (24 uM, TakaraBio, USA), 1 ul of SMARTer IIA Oligo (24 uM, TakaraBio, USA), 1 ul of 10x Template Switching Reverse Transcriptase (New England Biolabs, MA, USA) for the synthesis of the 1st strand, and 1 ul of Primer IIA (12 uM, TakaraBio, USA) and 1 ul of PrimeSTAR GXL polymerase (TakaraBio, USA) for the synthesis of the 2nd strand were used following manufacturer user manual instructions (TakaraBio, USA; New England Biolabs, MA, USA). The synthesized cDNA was purified using SPRI AMPure beads according to manufacturer’s instructions (TakaraBio, USA), and its integrity and quality were assessed with Qubit 4 Fluorometer (Qubit™4 Fluorometer, Thermo Fisher Scientific, USA) in accordance with One time dsDNA HS assay kit user manual. To prepare the library for nanopore sequencing, the ONT Rapid Barcoding Sequencing Kits (SQK-RPB004, Oxford Nanopore Technologies, UK) were used as per the manufacturer’s guidelines. Libraries were pooled on FLO-MIN106 flow cells and run on the GridION. Sequencing performance was monitored and was terminated after 24 h. Metagenomic analysis. Sequencing reads were filtered to a minimum length (≥200 bp) and Q-value (≥9) by MinKNOW v4.3.4. Basecalling and demultiplexing was performed using Guppy v6.4 with the high accuracy model. Guppy is only available to NanoPore customers through their community site (https://community.nanoporetech.com). PoreChop v0.2.4 [58] was used to remove the nanopore barcode adapter sequences. To assemble metagenome-assembled contigs, the quality filtered reads were assembled by Canu v2.2 [59] using the following assembly parameters: -nanopore maxInputCoverage = 2000 corOutCoverage = all corMinCoverage = 0 corMhapSensitivity = high minoverlap = 50 minread = 200 genomesize = 5000. Of 48 viromes, 42 were successfully assembled to generate contigs (24 NucleoSpin® and 18 NucleoMag™). Contigs with a minimum length of 1 kbp were binned manually with the anvi’o v7.1 [60, 61] interactive interface. Briefly, anvi’o profiled the contigs using Prodigal v2.6.3 [62], with default parameters, then reads were mapped to the contig database using Minimap2 v2.24 [63], and the read recruitment was stored as a BAM file using samtools. Anvi’o profiles each BAM file, estimating the read coverage and detection statistics of each contig, we then normalized coverage as Reads per kilobase of transcript per Million reads mapped (RPKM) using the python package bioinforkit v2.1.0 [64]. Coverage and RPKM was combined into a merged profile database. The contigs database was populated with additional data, incorporating HMMER results against Virus Orthologous Groups (VOGs; https://vogdb.org/) in addition to the standard anvi’o HMMR profiles, NCBI COGs and KEGG KOfam database [65]. Contig taxonomy was predicted by running Kraken2 v2.1.2 [66, 67] using the non-redundant NCBI database on the gene calls. Finally, merged profiles were clustered with the automatic binning algorithm CONCOCT, and the anvi’o profile was visualized for manual binning. Binning was guided by sequence composition similarity (visualized as a dendrogram in the Anvi’o interface), and the presence of viral HMM hits to the VOG database. Validation of honey bee RNA virus genomes from binned contigs. Binned contigs were size filtered to 5 kbp, then aligned to reference genomes using Minimap2. References for complete genomes included DWV-A (NC_004830.2) and DWV-B/VDV-1 (NC_006494.1). Alignments were used to identify probable viral genomes. Prodigal was utilized to identify coding regions and regions were annotated using BLASTn (Nucleotide Basic Local Alignment Search Tool). Data analysis and visualization. Data was visualized and statistically analyzed using R v4.1.0 in RStudio build 576 and Microsoft Excel software. Statistical tests were performed using base features in R and data visualized with the package ggplot2 v3.4.1. Numeric values of the read length, number of reads and average quality of each read was acquired using SeqKit v2.3.0 [68]. Welch approximation t-test was used to compare data generation outputs from NucleoSpin® versus NucleoMag™ (i.e., RNA yields, number of reads generated, read length, and average read quality, average read mapped to reference genomes). Results RNA extraction and nanopore sequence data generation RNA was isolated twice from the set of 24 honey bee pooled samples each representing a colony, once using a manual column affinity method (NucleoSpin®), and second time employing an automated magnetic beads technology-based method (NucleoMag™). In total RNA was extracted from 48 honey bee pools (Table 1), and the yield was quantified (ng/μl) and assessed for quality (260/280 and 260/230 ratios, S1 Table). The NucleoMag™ protocol did not have comparable yields to NucleoSpin® (t(45) = -4.4263, p-value = 0.0001; Fig 1A), or purity represented as 260/280 ratios (t(45) = -2.329, p-value = 0.0259; Fig 1B), as assessed by Welch t-test (S2 Table). Despite the differences between the means, both methods fell within an acceptable range that ensures the greatest likelihood of successful sequencing. Furthermore, NucleoMag™ benefited from a high consistency in RNA yields. Overall, we observed a similar distribution between the two datasets in regard to read length and quality (Fig 1C and 1D), with NucleoMag™ generating a greater abundance of data. Honey bee cDNA samples were sequenced, generating on average 41,098 and 113,719 reads per sample from NucleoSpin® and NucleoMag™, respectively. A total of 2,024 contigs were assembled from the dataset, with more contigs assembled from samples extracted with NucleoMag™ (1,444) than NucleoSpin® (580) (Fig 1E). To explore this dataset further we conducted multiple Welch t-tests, comparing the means of the two extraction protocols for the number, quality and the length of sequenced reads and contigs. Of the five tests conducted, four revealed significant differences between the means of the two groups (p-value <0.05) (S2 Table). Specifically, we found differences between the means of NucleoMag™ and NucleoSpin® on: the number of reads generated (t(45) = 3.4161, p-value = 0.001); read length (t(100275) = -8.443659, p-value = <2.2x10-16); average read quality (t(100998) = 69.5, p-value = <2.2x10-16) and number of contigs assembled per sample (t(18.427) = 3.5612, p-value = 0.002; Fig 1E). No differences between the means were observed in contig lengths (t(771.86) = -0.01815, p-value = 0.9855). With the exception of read length and contig length, we observe the trend of NucleoMag™ having higher means as compared to column extraction method. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 1. RNA extraction with NucleoMag™ is comparable, and by some measure preferable, to NucleoSpin®. Overview of the RNA quantity and quality as represented by (A) RNA concentration (ng/μl) and (B) quality measured by spectrophotometry (260/280 ratio). Opaque red zone indicates the range of RNA purity (~2–2.2) which is the best practice quality to ensure the greatest likelihood of sequencing success. Histograms with bin size of 100 of (C) read lengths, and (D) average read quality generated from RNA extracted by either NucleoMag™ (dark gray) or NucleoSpin® (magenta), and the (E) number of de novo contigs assembled by CANU. https://doi.org/10.1371/journal.pone.0297623.g001 RT-qPCR. Forty-eight RNA samples were screened for DWV-A and DWV-B variants’ presence and quantity applying RT-qPCR-based ABC assay [57]. Both methods provided similar analytical results with a slight difference in the strength of amplification signal reflected as 1.0–1.1 in log10 viral genome copies per colony more for RNA isolated with NucleoSpin® compared to NucleoMag™. The analysis of correlation between the Ct values obtained from both RNA sets, as well as between viral genome copies per colony showed a high correlation for both DWV strains (R2 = 0.9992 and R2 = 0.9988 for DWV-A; R2 = 0.9976 and R2 = 0.9985 for DWV-B, for Ct values and viral genome copies respectively, Fig 2). Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 2. Correlation between the two methods of extraction by Ct value and log10 viral genome copies per colony. (A) R2 = 0.9988 for viral genome copies, R2 = 0.9992 for Ct values for DWV-A; (B) R2 = 0.9985 for viral genome copies, R2 = 0.9976 for Ct values for DWV-B. https://doi.org/10.1371/journal.pone.0297623.g002 Out of 24 the colonies 11 were confirmed to be DWV-A and DWV-B free or below the limits of detection. All 11 colonies were sampled in April 2021 and included colonies 1, 2, 4, 9, 11–14, 16, 17 and 19, corresponding to RNA samples 8, 12, 22, 9, 10, 26, 23, 24, 27, 37 and 48, respectively (Table 1). Three out of these colonies–colony 1, 2 and 4 (RNA samples 64, 70 and 63 respectively, Table 1)—were sampled in fall 2020 as well. Colonies 1 and 2 showed moderate to high levels of DWV-B, (RNA samples ## 64 and 70, Table 1) showing 8.8 (NucleoMag™) and 9.4 (NucleoSpin®), and 7.9 (NucleoMag™) and 8.8 (NucleoSpin®) log10 DWV-B copies for colony 1 and 2, respectively (Table 1). Colony 4 was positive for both DWV-A and DWV-B (6.9 and 9.1 (NucleoMag™), and 7.8 and 9.8 (NucleoSpin®), respectively, corresponding to RNA sample 63, Table 1). Another colony sampled twice, in fall 2020 and spring 2021, colony 6, was found positive for both DWV master variants A (7.5 (NucleoMag™) and 8.4 (NucleoSpin®) in fall 2020; 7.3 (NucleoMag™) and 8.2 (NucleoSpin®) in spring 2021) and B, increasing drastically its DWV-B loads in spring 2021 (7.3 (NucleoMag™) and 7.8 (NucleoSpin®) in fall 2020; 10 (NucleoMag™) and 10.6 (NucleoSpin®) in spring 2021), (Table 1). Two more colonies sampled in September 2020 only, colony 3 and 5, showed moderate to high levels of DWV-A (7.5 (NucleoMag™) and 8.6 (NucleoSpin®), and 9.8 (NucleoMag™) and 10.3 (NucleoSpin®) log10 genome copies for colony 3 and 5, respectively); and high loads of DWV-B (9.1 (NucleoMag™) and 9.5 (NucleoSpin®), and 9.7 (NucleoMag™) and 9.9 (NucleoSpin®) log10 genome copies for colony 3 and 5 respectively), corresponding to RNA samples 1 and 2 (Table 1). In summary, fall virus levels observed were relatively high as expected. The rest of the colonies were sampled in April 2021 only, showing presence of either both DWV master variants A and B (colonies 7 and 18), or A (colonies 8 and 10), or B (colonies 15 and 20). These spring viral quantities observed were relatively high for the season, ranging from 6.7 to 8.8 (NucleoMag™) and 7.6 to 9.6 (NucleoSpin®), and 6.4 to 8.9 (NucleoMag™) and 7.3 to 9.2 (NucleoSpin®), log10 DWV-A and DWV-B genome copies per colony respectively (Table 1, RNA samples 7, 8, 10, 15, 18, 20). Coverage across DWV genomes. A total of 21 contigs between 2.23 and 10.14 kbp in length, were binned initially as DWV genomes or genome fragments. After removing contigs less than 5 kbp and aligning contigs to DWV reference genomes, ten contigs were putatively classified, after BLASTn annotation, as DWV-B (min = 6.57 kbp, max = 10.14), a single DWV-A contig 7.00 kbp in length, and a two recombinant DWV contigs, one of 6.6 kbp and the other 9.7 kbp in length. Virome reads were mapped to reference DWV genomes (DWV-A, NC_004830.2; DWV-B/VDV-1, NC_006494.1). DWV-B was the most prevalent genome variant within the dataset, as demonstrated by both the abundance of DWV-B contigs assembled, proportion of sample with reads mapping to the reference genome and RT-qPCR data (Fig 3). Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 3. Both extraction protocols can generate high coverage of both DWV-A and DWV-B genomes that are generally in agreement with RT-qPCR data. The number of reads mapping across the reference genomes of (A) DWV-A and (B) DWV-B. Data were filtered to include only samples that mapped with >10 reads on average, removing noisy low-coverage data. The entire (unfiltered) dataset is presented as the number of reads mapped is summarized by (C) RT-qPCR detection and (D) by the total data mapped to the two genomes. In all plots genome copy numbers from the amplification of RdRp from RT-qPCR are represented as the boxplot fill color, while the line color indicates the extraction method. (Pos., Positive; Neg., Negative; *, colony sampled for DNA extraction in the year 2020). https://doi.org/10.1371/journal.pone.0297623.g003 Overall, we observed a similar abundance of reads mapping to both DWV-A (Fig 3A) and DWV-B (Fig 3B) from both extraction protocols (Fig 3D). A Welch t-test was performed to examine the effect of extraction protocol on the average reads per sample mapped to the reference genomes. No significant differences were found between the two extraction protocol for read mapped to DWV-A (t(26.12) = -0.711, p-value = 0.483) or DWV-B (t(27.40) = -1.351, p-value = 0.187; Fig 3D), and the means from NucleoMag™ were larger for both measures. We do observe colonies with absent mapping for one extraction method, such as colony 7 (Fig 3A). This may be explained by uneven sequencing library sizes (colony 7: NucleoSpin 1.99x105 reads and NucleoMag 6.34x105 reads). Regardless of extraction protocols, full genomes of DWV-A and B could be recovered from colonies (Fig 4). This does not always concur with the results of RT-qPCR screening, and individual inspection of the read mapping revealed likely fragmented or deteriorated genomes within the samples, as mapping did not occur around the RdRp region of the genome. Furthermore, both protocols were able to accurately identify a recombinant strain, notable due to reads mapping to regions of both DWV-A and B (Colony Number 8; Fig 3A and 3B), but detection by RT-qPCR to DWV-A (S1 Fig). In general, we observed a congruence between the RT-qPCR results and read mapping depth. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 4. Both extraction protocols can generate high coverage of both DWV-A and DWV-B genomes that is reflected in RT-qPCR detection. Colony #5, Location: Hill, Year: 2020. Coverage histograms (left) represented as the number of reads mapped to two Deformed wing virus’ strains from NCBI, (top) DWV-A (NC_004830.2) and (bottom) DWV-B/VDV-1 (NC_006494.1). Copy number of RNA-dependent RNA polymerase for the two DWV strains as detected by RT-qPCR (right). Error bars in the barplot represent the standard deviation of the technical replicas. https://doi.org/10.1371/journal.pone.0297623.g004 RNA extraction and nanopore sequence data generation RNA was isolated twice from the set of 24 honey bee pooled samples each representing a colony, once using a manual column affinity method (NucleoSpin®), and second time employing an automated magnetic beads technology-based method (NucleoMag™). In total RNA was extracted from 48 honey bee pools (Table 1), and the yield was quantified (ng/μl) and assessed for quality (260/280 and 260/230 ratios, S1 Table). The NucleoMag™ protocol did not have comparable yields to NucleoSpin® (t(45) = -4.4263, p-value = 0.0001; Fig 1A), or purity represented as 260/280 ratios (t(45) = -2.329, p-value = 0.0259; Fig 1B), as assessed by Welch t-test (S2 Table). Despite the differences between the means, both methods fell within an acceptable range that ensures the greatest likelihood of successful sequencing. Furthermore, NucleoMag™ benefited from a high consistency in RNA yields. Overall, we observed a similar distribution between the two datasets in regard to read length and quality (Fig 1C and 1D), with NucleoMag™ generating a greater abundance of data. Honey bee cDNA samples were sequenced, generating on average 41,098 and 113,719 reads per sample from NucleoSpin® and NucleoMag™, respectively. A total of 2,024 contigs were assembled from the dataset, with more contigs assembled from samples extracted with NucleoMag™ (1,444) than NucleoSpin® (580) (Fig 1E). To explore this dataset further we conducted multiple Welch t-tests, comparing the means of the two extraction protocols for the number, quality and the length of sequenced reads and contigs. Of the five tests conducted, four revealed significant differences between the means of the two groups (p-value <0.05) (S2 Table). Specifically, we found differences between the means of NucleoMag™ and NucleoSpin® on: the number of reads generated (t(45) = 3.4161, p-value = 0.001); read length (t(100275) = -8.443659, p-value = <2.2x10-16); average read quality (t(100998) = 69.5, p-value = <2.2x10-16) and number of contigs assembled per sample (t(18.427) = 3.5612, p-value = 0.002; Fig 1E). No differences between the means were observed in contig lengths (t(771.86) = -0.01815, p-value = 0.9855). With the exception of read length and contig length, we observe the trend of NucleoMag™ having higher means as compared to column extraction method. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 1. RNA extraction with NucleoMag™ is comparable, and by some measure preferable, to NucleoSpin®. Overview of the RNA quantity and quality as represented by (A) RNA concentration (ng/μl) and (B) quality measured by spectrophotometry (260/280 ratio). Opaque red zone indicates the range of RNA purity (~2–2.2) which is the best practice quality to ensure the greatest likelihood of sequencing success. Histograms with bin size of 100 of (C) read lengths, and (D) average read quality generated from RNA extracted by either NucleoMag™ (dark gray) or NucleoSpin® (magenta), and the (E) number of de novo contigs assembled by CANU. https://doi.org/10.1371/journal.pone.0297623.g001 RT-qPCR. Forty-eight RNA samples were screened for DWV-A and DWV-B variants’ presence and quantity applying RT-qPCR-based ABC assay [57]. Both methods provided similar analytical results with a slight difference in the strength of amplification signal reflected as 1.0–1.1 in log10 viral genome copies per colony more for RNA isolated with NucleoSpin® compared to NucleoMag™. The analysis of correlation between the Ct values obtained from both RNA sets, as well as between viral genome copies per colony showed a high correlation for both DWV strains (R2 = 0.9992 and R2 = 0.9988 for DWV-A; R2 = 0.9976 and R2 = 0.9985 for DWV-B, for Ct values and viral genome copies respectively, Fig 2). Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 2. Correlation between the two methods of extraction by Ct value and log10 viral genome copies per colony. (A) R2 = 0.9988 for viral genome copies, R2 = 0.9992 for Ct values for DWV-A; (B) R2 = 0.9985 for viral genome copies, R2 = 0.9976 for Ct values for DWV-B. https://doi.org/10.1371/journal.pone.0297623.g002 Out of 24 the colonies 11 were confirmed to be DWV-A and DWV-B free or below the limits of detection. All 11 colonies were sampled in April 2021 and included colonies 1, 2, 4, 9, 11–14, 16, 17 and 19, corresponding to RNA samples 8, 12, 22, 9, 10, 26, 23, 24, 27, 37 and 48, respectively (Table 1). Three out of these colonies–colony 1, 2 and 4 (RNA samples 64, 70 and 63 respectively, Table 1)—were sampled in fall 2020 as well. Colonies 1 and 2 showed moderate to high levels of DWV-B, (RNA samples ## 64 and 70, Table 1) showing 8.8 (NucleoMag™) and 9.4 (NucleoSpin®), and 7.9 (NucleoMag™) and 8.8 (NucleoSpin®) log10 DWV-B copies for colony 1 and 2, respectively (Table 1). Colony 4 was positive for both DWV-A and DWV-B (6.9 and 9.1 (NucleoMag™), and 7.8 and 9.8 (NucleoSpin®), respectively, corresponding to RNA sample 63, Table 1). Another colony sampled twice, in fall 2020 and spring 2021, colony 6, was found positive for both DWV master variants A (7.5 (NucleoMag™) and 8.4 (NucleoSpin®) in fall 2020; 7.3 (NucleoMag™) and 8.2 (NucleoSpin®) in spring 2021) and B, increasing drastically its DWV-B loads in spring 2021 (7.3 (NucleoMag™) and 7.8 (NucleoSpin®) in fall 2020; 10 (NucleoMag™) and 10.6 (NucleoSpin®) in spring 2021), (Table 1). Two more colonies sampled in September 2020 only, colony 3 and 5, showed moderate to high levels of DWV-A (7.5 (NucleoMag™) and 8.6 (NucleoSpin®), and 9.8 (NucleoMag™) and 10.3 (NucleoSpin®) log10 genome copies for colony 3 and 5, respectively); and high loads of DWV-B (9.1 (NucleoMag™) and 9.5 (NucleoSpin®), and 9.7 (NucleoMag™) and 9.9 (NucleoSpin®) log10 genome copies for colony 3 and 5 respectively), corresponding to RNA samples 1 and 2 (Table 1). In summary, fall virus levels observed were relatively high as expected. The rest of the colonies were sampled in April 2021 only, showing presence of either both DWV master variants A and B (colonies 7 and 18), or A (colonies 8 and 10), or B (colonies 15 and 20). These spring viral quantities observed were relatively high for the season, ranging from 6.7 to 8.8 (NucleoMag™) and 7.6 to 9.6 (NucleoSpin®), and 6.4 to 8.9 (NucleoMag™) and 7.3 to 9.2 (NucleoSpin®), log10 DWV-A and DWV-B genome copies per colony respectively (Table 1, RNA samples 7, 8, 10, 15, 18, 20). Coverage across DWV genomes. A total of 21 contigs between 2.23 and 10.14 kbp in length, were binned initially as DWV genomes or genome fragments. After removing contigs less than 5 kbp and aligning contigs to DWV reference genomes, ten contigs were putatively classified, after BLASTn annotation, as DWV-B (min = 6.57 kbp, max = 10.14), a single DWV-A contig 7.00 kbp in length, and a two recombinant DWV contigs, one of 6.6 kbp and the other 9.7 kbp in length. Virome reads were mapped to reference DWV genomes (DWV-A, NC_004830.2; DWV-B/VDV-1, NC_006494.1). DWV-B was the most prevalent genome variant within the dataset, as demonstrated by both the abundance of DWV-B contigs assembled, proportion of sample with reads mapping to the reference genome and RT-qPCR data (Fig 3). Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 3. Both extraction protocols can generate high coverage of both DWV-A and DWV-B genomes that are generally in agreement with RT-qPCR data. The number of reads mapping across the reference genomes of (A) DWV-A and (B) DWV-B. Data were filtered to include only samples that mapped with >10 reads on average, removing noisy low-coverage data. The entire (unfiltered) dataset is presented as the number of reads mapped is summarized by (C) RT-qPCR detection and (D) by the total data mapped to the two genomes. In all plots genome copy numbers from the amplification of RdRp from RT-qPCR are represented as the boxplot fill color, while the line color indicates the extraction method. (Pos., Positive; Neg., Negative; *, colony sampled for DNA extraction in the year 2020). https://doi.org/10.1371/journal.pone.0297623.g003 Overall, we observed a similar abundance of reads mapping to both DWV-A (Fig 3A) and DWV-B (Fig 3B) from both extraction protocols (Fig 3D). A Welch t-test was performed to examine the effect of extraction protocol on the average reads per sample mapped to the reference genomes. No significant differences were found between the two extraction protocol for read mapped to DWV-A (t(26.12) = -0.711, p-value = 0.483) or DWV-B (t(27.40) = -1.351, p-value = 0.187; Fig 3D), and the means from NucleoMag™ were larger for both measures. We do observe colonies with absent mapping for one extraction method, such as colony 7 (Fig 3A). This may be explained by uneven sequencing library sizes (colony 7: NucleoSpin 1.99x105 reads and NucleoMag 6.34x105 reads). Regardless of extraction protocols, full genomes of DWV-A and B could be recovered from colonies (Fig 4). This does not always concur with the results of RT-qPCR screening, and individual inspection of the read mapping revealed likely fragmented or deteriorated genomes within the samples, as mapping did not occur around the RdRp region of the genome. Furthermore, both protocols were able to accurately identify a recombinant strain, notable due to reads mapping to regions of both DWV-A and B (Colony Number 8; Fig 3A and 3B), but detection by RT-qPCR to DWV-A (S1 Fig). In general, we observed a congruence between the RT-qPCR results and read mapping depth. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 4. Both extraction protocols can generate high coverage of both DWV-A and DWV-B genomes that is reflected in RT-qPCR detection. Colony #5, Location: Hill, Year: 2020. Coverage histograms (left) represented as the number of reads mapped to two Deformed wing virus’ strains from NCBI, (top) DWV-A (NC_004830.2) and (bottom) DWV-B/VDV-1 (NC_006494.1). Copy number of RNA-dependent RNA polymerase for the two DWV strains as detected by RT-qPCR (right). Error bars in the barplot represent the standard deviation of the technical replicas. https://doi.org/10.1371/journal.pone.0297623.g004 RT-qPCR. Forty-eight RNA samples were screened for DWV-A and DWV-B variants’ presence and quantity applying RT-qPCR-based ABC assay [57]. Both methods provided similar analytical results with a slight difference in the strength of amplification signal reflected as 1.0–1.1 in log10 viral genome copies per colony more for RNA isolated with NucleoSpin® compared to NucleoMag™. The analysis of correlation between the Ct values obtained from both RNA sets, as well as between viral genome copies per colony showed a high correlation for both DWV strains (R2 = 0.9992 and R2 = 0.9988 for DWV-A; R2 = 0.9976 and R2 = 0.9985 for DWV-B, for Ct values and viral genome copies respectively, Fig 2). Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 2. Correlation between the two methods of extraction by Ct value and log10 viral genome copies per colony. (A) R2 = 0.9988 for viral genome copies, R2 = 0.9992 for Ct values for DWV-A; (B) R2 = 0.9985 for viral genome copies, R2 = 0.9976 for Ct values for DWV-B. https://doi.org/10.1371/journal.pone.0297623.g002 Out of 24 the colonies 11 were confirmed to be DWV-A and DWV-B free or below the limits of detection. All 11 colonies were sampled in April 2021 and included colonies 1, 2, 4, 9, 11–14, 16, 17 and 19, corresponding to RNA samples 8, 12, 22, 9, 10, 26, 23, 24, 27, 37 and 48, respectively (Table 1). Three out of these colonies–colony 1, 2 and 4 (RNA samples 64, 70 and 63 respectively, Table 1)—were sampled in fall 2020 as well. Colonies 1 and 2 showed moderate to high levels of DWV-B, (RNA samples ## 64 and 70, Table 1) showing 8.8 (NucleoMag™) and 9.4 (NucleoSpin®), and 7.9 (NucleoMag™) and 8.8 (NucleoSpin®) log10 DWV-B copies for colony 1 and 2, respectively (Table 1). Colony 4 was positive for both DWV-A and DWV-B (6.9 and 9.1 (NucleoMag™), and 7.8 and 9.8 (NucleoSpin®), respectively, corresponding to RNA sample 63, Table 1). Another colony sampled twice, in fall 2020 and spring 2021, colony 6, was found positive for both DWV master variants A (7.5 (NucleoMag™) and 8.4 (NucleoSpin®) in fall 2020; 7.3 (NucleoMag™) and 8.2 (NucleoSpin®) in spring 2021) and B, increasing drastically its DWV-B loads in spring 2021 (7.3 (NucleoMag™) and 7.8 (NucleoSpin®) in fall 2020; 10 (NucleoMag™) and 10.6 (NucleoSpin®) in spring 2021), (Table 1). Two more colonies sampled in September 2020 only, colony 3 and 5, showed moderate to high levels of DWV-A (7.5 (NucleoMag™) and 8.6 (NucleoSpin®), and 9.8 (NucleoMag™) and 10.3 (NucleoSpin®) log10 genome copies for colony 3 and 5, respectively); and high loads of DWV-B (9.1 (NucleoMag™) and 9.5 (NucleoSpin®), and 9.7 (NucleoMag™) and 9.9 (NucleoSpin®) log10 genome copies for colony 3 and 5 respectively), corresponding to RNA samples 1 and 2 (Table 1). In summary, fall virus levels observed were relatively high as expected. The rest of the colonies were sampled in April 2021 only, showing presence of either both DWV master variants A and B (colonies 7 and 18), or A (colonies 8 and 10), or B (colonies 15 and 20). These spring viral quantities observed were relatively high for the season, ranging from 6.7 to 8.8 (NucleoMag™) and 7.6 to 9.6 (NucleoSpin®), and 6.4 to 8.9 (NucleoMag™) and 7.3 to 9.2 (NucleoSpin®), log10 DWV-A and DWV-B genome copies per colony respectively (Table 1, RNA samples 7, 8, 10, 15, 18, 20). Coverage across DWV genomes. A total of 21 contigs between 2.23 and 10.14 kbp in length, were binned initially as DWV genomes or genome fragments. After removing contigs less than 5 kbp and aligning contigs to DWV reference genomes, ten contigs were putatively classified, after BLASTn annotation, as DWV-B (min = 6.57 kbp, max = 10.14), a single DWV-A contig 7.00 kbp in length, and a two recombinant DWV contigs, one of 6.6 kbp and the other 9.7 kbp in length. Virome reads were mapped to reference DWV genomes (DWV-A, NC_004830.2; DWV-B/VDV-1, NC_006494.1). DWV-B was the most prevalent genome variant within the dataset, as demonstrated by both the abundance of DWV-B contigs assembled, proportion of sample with reads mapping to the reference genome and RT-qPCR data (Fig 3). Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 3. Both extraction protocols can generate high coverage of both DWV-A and DWV-B genomes that are generally in agreement with RT-qPCR data. The number of reads mapping across the reference genomes of (A) DWV-A and (B) DWV-B. Data were filtered to include only samples that mapped with >10 reads on average, removing noisy low-coverage data. The entire (unfiltered) dataset is presented as the number of reads mapped is summarized by (C) RT-qPCR detection and (D) by the total data mapped to the two genomes. In all plots genome copy numbers from the amplification of RdRp from RT-qPCR are represented as the boxplot fill color, while the line color indicates the extraction method. (Pos., Positive; Neg., Negative; *, colony sampled for DNA extraction in the year 2020). https://doi.org/10.1371/journal.pone.0297623.g003 Overall, we observed a similar abundance of reads mapping to both DWV-A (Fig 3A) and DWV-B (Fig 3B) from both extraction protocols (Fig 3D). A Welch t-test was performed to examine the effect of extraction protocol on the average reads per sample mapped to the reference genomes. No significant differences were found between the two extraction protocol for read mapped to DWV-A (t(26.12) = -0.711, p-value = 0.483) or DWV-B (t(27.40) = -1.351, p-value = 0.187; Fig 3D), and the means from NucleoMag™ were larger for both measures. We do observe colonies with absent mapping for one extraction method, such as colony 7 (Fig 3A). This may be explained by uneven sequencing library sizes (colony 7: NucleoSpin 1.99x105 reads and NucleoMag 6.34x105 reads). Regardless of extraction protocols, full genomes of DWV-A and B could be recovered from colonies (Fig 4). This does not always concur with the results of RT-qPCR screening, and individual inspection of the read mapping revealed likely fragmented or deteriorated genomes within the samples, as mapping did not occur around the RdRp region of the genome. Furthermore, both protocols were able to accurately identify a recombinant strain, notable due to reads mapping to regions of both DWV-A and B (Colony Number 8; Fig 3A and 3B), but detection by RT-qPCR to DWV-A (S1 Fig). In general, we observed a congruence between the RT-qPCR results and read mapping depth. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 4. Both extraction protocols can generate high coverage of both DWV-A and DWV-B genomes that is reflected in RT-qPCR detection. Colony #5, Location: Hill, Year: 2020. Coverage histograms (left) represented as the number of reads mapped to two Deformed wing virus’ strains from NCBI, (top) DWV-A (NC_004830.2) and (bottom) DWV-B/VDV-1 (NC_006494.1). Copy number of RNA-dependent RNA polymerase for the two DWV strains as detected by RT-qPCR (right). Error bars in the barplot represent the standard deviation of the technical replicas. https://doi.org/10.1371/journal.pone.0297623.g004 Discussion The NucleoMag™ extraction protocol employed in this study for isolation of viral RNA that was used for subsequent virus screening with RT-qPCR and sequencing was shown to produce as comprehensive and reproducible results as the NucleoSpin® protocol. Despite the differences in RNA yields and purity, both methods fell within an acceptable range required for successful downstream applications for accurate virus detection such as real-time virus quantification assay and cDNA library construction used for sequencing. In this study we demonstrate that automated RNA extraction using NucleoMag™ achieves similar quality and quantity of RNA. However, with respect to sequence data, magnetic beads-based technology could exceed the quantity generated by NucleoSpin®. Spin columns technology for the extraction of total RNA from honey bees is widely used [69–71], while automated magnetic bead-based extraction protocols are not, despite benefiting from reduced manual labor (high-throughput) and improved consistency between yields. Our data demonstrates that the sensitivity of DWV detection by both RT-qPCR and nanopore sequencing is comparable after total RNA extraction by NucleoSpin® Virus and NucleoMag™ Virus isolation kits. Both methods provided similar analytical results obtained from RT-qPCR showing that both DWV strain loads and Ct values from honey bee pooled samples processed by affinity purification columns and magnetic beads-based technology method correlated well, which is in agreement with previous work [37, 38]. Differences observed in the strength of amplification signal could be due to higher RNA yields and purity delivered by NucleoSpin®. In addition, sample viscosity is known to impact magnetic beads performance by impeding migration of the beads (personal communication with a manufacturer), hence impacting both the RNA yields and purity. Yet, samples extracted with NucleoSpin® generated lower sequencing data outcome, possibly due to a decreased ability of the affinity purification columns to recover small fragments of nucleic acids efficiently, as small fragments bind tightly with the silica matrix [42]. Overall, these findings provide insights into the differences RNA extraction protocol can have on the data generated from sequencing honey bee viromes. While differences did occur in the RNA yields between the two strategies, this did not negatively impact sequencing success or achieving sufficient coverage of DWV genomes. This outcome agrees with previous work, such as a study which compared NucleoMag™ with NucleoSpin® Tissue DNA extraction kits on a range of forensic samples (e.g. human tissues), which likewise demonstrated that NucleoMag™ is suitable without compromising RNA yields or quality [72]. Elsewhere, and more relevant to the study of RNA viruses, similar results have been presented for the extraction of Enterovirus RNA, comparing automated and column extraction protocols for the purpose of RT-qPCR [73]. In general, there was a congruence between the RT-qPCR results and read mapping depth. Still, in some cases we found a mismatch between RT-qPCR and sequencing results revealing a limitation of the real-time assay to identify a recombinant DWV strain (S1 Fig, colony 8). Being a sensitive and accurate standard method for virus detection in many research fields, including honey bee research, RT-qPCR technique remains limited due to its specificity to the target genome location, while sequencing allows to target a whole genome. Due to a small sample set we were unable to further investigate the outcome of this study in terms of viral load comparison between fall and spring samples. However, our findings confirm compatibility and accuracy of automated magnetic beads-based technology for RNA extraction and subsequent virus detection both with real-time assay and nanopore sequencing, demonstrating comparable results with manual affinity purification column method, suitable for screening up to 96 samples at a time. Supporting information S1 Table. 48 samples of viral RNA isolated from 24 honey bee pooled samples. Queen label represents a queen line. Viral load expressed as log10 viral RNA copies/ a pooled sample. Total sample volume was 19 ml; Sample input volume was 200ul (50ul of filtrate + 150ul of 1xPBS) and 150ul (50ul of filtrate + 100ul of 1xPBS) for magnetic beads and columns affinity extraction methods respectively. https://doi.org/10.1371/journal.pone.0297623.s001 (DOCX) S2 Table. Summary of results from multiple Welch two-sample t-tests. Significant p-value results (<0.05) are in bold. (M, mean; SD, standard deviation; N, number; df, degrees of freedom). https://doi.org/10.1371/journal.pone.0297623.s002 (DOCX) S1 Fig. Both extraction protocols are capable of identifying chimeric DWV genomes. Colony # 8, Location: Hill, Year: 2021. Coverage histograms (left) represented as the number of reads mapped to two Deformed wing virus’ strains from NCBI, (top) DWV-A (NC_004830.2) and (bottom) DWV-B/VDV-1 (NC_006494.1). Copy number of RNA-dependent RNA polymerase for the two DWV strains as detected by RT-qPCR (right). https://doi.org/10.1371/journal.pone.0297623.s003 (DOCX)
TI - A semi-automated and high-throughput approach for the detection of honey bee viruses in bee samples
JF - PLoS ONE
DO - 10.1371/journal.pone.0297623
DA - 2024-03-14
UR - https://www.deepdyve.com/lp/public-library-of-science-plos-journal/a-semi-automated-and-high-throughput-approach-for-the-detection-of-gkYQgTN1fg
SP - e0297623
VL - 19
IS - 3
DP - DeepDyve
ER -