A Comparison of Deformed Wing Virus in Deformed and Asymptomatic Honey Bees

2.1. Honeybee Samples

For the RT-qPCR viral quantification and wing deformity study, honey bees were collected from an apiary maintained by Universidade Federal do Recôncavo da Bahia (UFRB), Cruz das Almas (12.67° S, 39.1019° W), Bahia state, Brazil. We confirmed that the honey bees from this population were infected by the type A variant of DWV (Supplementary Figure S1) as were the bees from Hawaii [2]. For the NGS study, we used previously collected samples of honey bees from Hawaii (collected in 2012). Three pairs of samples were chosen, each consisting of a single deformed bee and a pool of 30 asymptomatic bees from the same colony. Each colony had been exposed to Varroa infestations for different lengths of time (Oahu = 5 years, Big Island East = 4 years, Big Island South = 3 years). The rarity of bees with deformed wings in both populations made it impossible to compare similar numbers of symptomatic and asymptomatic bees. The vast majority of adults with deformed wings contain high DWV loads ([13], this study), whereas DWV titres in individual asymptomatic bees are more variable [33], hence we used a pooled sample for the NGS study to ensure sufficient DWV genomes were present for sequencing.

2.2. Effect of Viral Load on Wing Deformity

In January 2015, in Brazil, a frame of an emerging brood was removed from three study colonies. Each emerging bee, along with the cell that it was emerging from, was checked for the presence or absence of Varroa mites. A total of 45 parasitised and 45 non-parasitised newly emerged worker bees were collected from the three frames. Only bees seen emerging from a cell were used. However, no emerging deformed bees were found despite 500 bees emerging from the frames. As such, a visual search of the three study colonies was conducted that resulted in just three deformed bees being located. As the vast majority of deformed bees develop from parasitised pupae [34], it is likely that these individuals emerged from infested cells, which is supported by the high DWV titres we detected (Figure 1). All bees were killed by freezing at −20 °C before their forewings were removed and mounted on a glass slide for morphometric analysis. Individual bees were then labelled with a unique label and shipped to the UK in a Dry Vapour Shipper at −186 °C for viral analysis. Each forewing (length and width) was measured using a Leica binocular microscope (×10) magnification fitted with a Leica camera. As these were newly emerged bees, no wing wear was present. As directional asymmetry in wing size in honey bees is well established [35,36], we measured both wings of the parasitised and non-parasitised groups and compared the results using a Mann–Whitney U test since not all wing measurement distributions were normally distributed.

For the RT-qPCR analysis, a random subset of ten parasitised and ten non-parasitised newly emerged asymptomatic bees were chosen along with the three deformed bees. Then each of the 23 individual bee samples was ground in liquid Nitrogen to a fine homogeneous powder and 30 mg material used for RNA extraction using the RNeasy mini kit (Qiagen, Venlo, The Netherlands), according to the manufacturer’s instructions. Total RNA samples were quantified using a Nanodrop 8000 (Thermo Scientific, Waltham, MA, USA). One microgram of isolated RNA was treated with DNase I (Promega, Madison, WI, USA), followed by Nanodrop quantification to standardise the amounts of total RNA to 25 ng/µL, before storage at −80 °C.

Total RNA was analysed in duplicate for each sample using the one-step SensiFAST SYBR No ROX One-step kit (Bioline, London, UK). RT-qPCR reactions contained 50 ng RNA, 1× SYBR one-step Sensimix, 2.5 mM MgCl₂, 5 units of RNase inhibitor, and 7.5 pmol of each primer: DWVQ-F1 and R1 for DWV (primers bind within the RdRp gene) with Actin F1 and R1 as the reference gene [11] (Supplementary Table S1). Reactions were run on a Rotor-Gene Q Thermocycler (Qiagen) with an initial reverse transcription stage at 49 °C for 30 min and a denaturation step of 95 °C for 10 min, followed by 40 cycles of denaturation for 15 s at 95 °C, annealing for 30 s at 54 °C for DWV, and 58 °C for Actin, and extension for 20 s at 72 °C. The SYBR green signal was measured on the green channel after each extension step. A final dissociation melt curve was performed between 65 °C and 95 °C, at 0.5 °C increments, each with a 10 s hold and acquisition to the green channel. The melt curve was used to ensure that a single targeted product was amplified, and that no contamination was present in the reverse transcription negative controls or in the no-template controls. The threshold cycle (C_t) value was determined for each sample using the QIAGEN Rotor—Gene Q Series Analysis software. All samples were run in duplicate and the average was taken. Those samples which had a standard deviation of ≥1 C_t were re-run to obtain duplicates with standard deviation <1 C_t. Each sample was normalised against Actin, and then presented relative to the asymptomatic non-parasitised bees as ΔΔC_t values. Statistical differences were calculated using a pair wise Mann–Whitney U test when the data were not normally distributed.

2.3. Next Generation Sequencing, Assembly and Data Normalisation

RNA was again extracted from 30 mg of material from each of the six Hawaiian samples (three colonies that each contained a single deformed bee and 30 asymptomatic bees) using the RNeasy mini kit (see above). Total RNA was used for a cDNA amplification step using oligo dT priming followed by sequencing. Illumina sequencing (Hi-Seq 100 bp paired end reads) was carried out by The Earlham Institute, Norwich. A Bioinformatics pipeline designed to accommodate the large amount of variation found within DWV, first described in Mordecai et al. [30] was applied. This involved using reads which mapped to a custom BLAST database of DWV master variants type A (NC_004830.2 and Kakugo virus NC_005876.1), B (AY251269.2) and C (CEND01000001.1) database using an e value of 10e−05 to assemble DWV-like contigs using VICUNA which was specifically developed to deal with highly variable data. Read data were uploaded to the NCBI Sequence Read Archive under study number SRP095247.

Viral contigs were imported into Geneious (Version 7.04, created by Biomatters, Aukland, New Zealand) and the “Map to Reference tool” was used to align the contigs against the DWV type A (NC_004830.2), B (AY251269.2) and C (ERS657949) reference genomes. These contigs were used to assess the breadth of genome coverage (Figure 2) as well as to analyse the dominant variants in each sample (Figure S3). The phylogenetic trees (Figure 3) were created within Geneious (Version 7.04, created by Biomatters) using a Tamura-Nei Genetic Distance model and a neighbor joining tree building method. In order to ensure that the contigs produced truly represented the viral populations, geneious competitive alignments were performed in which the raw sequencing read 1 files in fasta format were competitively aligned against DWV types A, B and C reference genomes (allowing for 5% mismatches and no gaps with reads with multiple best matches being discarded) to produce coverage graphs for reads corresponding to each type.

To quantify the number of DWV reads in each sample, the number of reads unambiguously mapping to each master variant were counted. These were then normalised against the number of reads that matched to part of the Actin gene used for RT-qPCR by Highfield et al. [11] using the geneious “map to reference” tool.

3.1. Viral Quantification and Wing Deformity in Honey Bees

All 23 individual bees from Brazil tested positive for DWV using RT-qPCR. HRM analysis and Sanger sequencing indicated that all bees were dominated by the DWV type A (Supplementary Figures S1 and S2). The highest loads were consistently detected in the three bees with wing deformities (Figure 1b). This was followed by asymptomatic bees that had been parasitised by Varroa mites as pupae. The lowest DWV loads were detected in asymptomatic bees that had developed free from Varroa (non-parasitised bees). Due to the low viral load in the non-parasitised bees, primer dimer was also amplified along with the DWV RdRp diagnostic fragment and this would have led to an overestimate in viral load, so the actual amount may be lower than shown. Despite these significant differences in DWV load between the three groups, their wing morphology did not follow the same trend. The bees with deformed wings had the highest viral load, as expected. However, no significant differences in wing length (Left wing, U = 668, Z = 1.21, p = 0.22: Right wing, U = 709, Z = −1.53, p = 0.12) or wing width (Left wing, U = 694, Z = 0.91, p = 0.36: Right wing, U = 871, Z = 0.28, p = 0.78) were found (Figure 1c,d) between the asymptomatic, non-parasitised and parasitised groups of bees. In both groups, directional asymmetry was detected in wing length but not wing width. That is, the bees’ left wing was significantly longer than their right wing in both the non-parasitised (U = 493, Z = 3.71, p = 0.0002) and Varroa-parasitised (U = 402, Z = −3.52, p = 0.0004) groups (Figure 1c).

3.2. Next Generation Illumina Sequencing

Each of the six samples from Hawaii contained sufficient DWV reads to assemble into contigs that together covered the entire length of the genome (Figure 2). Five of the six samples were dominated by type A reads that yielded full genome coverage. However, a single asymptomatic sample from Big Island East was found to be dominated by type B reads that yielded full genome coverage of the type B genome (Figure 2). The breadth of type B coverage is less in the samples from Oahu where Varroa has been present for a long time where the type A variant dominates. The read depth coverage graphs for all samples showed a strong 3’ bias which can be attributed to the inherent 3’ bias of reverse transcription produced using oligo dT priming in the library preparation [37]. The plots confirmed that the asymptomatic sample from Big Island East was the only one that contained type B read coverage across the entire genome and that all others contained type B coverage at the 3’ end only (Figure 2). The samples from Oahu did contain low amounts of type B reads, however there was insufficient read depth for the assembler to produce contigs. The asymptomatic sample from Oahu contained a small number of DWV type C reads (n = 359). Although type C was not the dominant variant, coverage was sufficiently high to assemble contigs spanning the majority of the genome. Very low numbers of type C reads were found in other samples by counting reads unambiguously mapping to type C and were used when normalising DWV variants to actin (Figure 2), however the read depth was insufficient to generate type C contigs in any other sample. Although it is impossible to rule out that the DWV type C identified in this study is a result of contamination via barcode shifting originating from samples dominated by type C run on the same flow cell lane see [27], as evidenced by the divergent 3’ bias in the type C read density for HB_S67 (Figure 2), we suggest the small number of sequence differences in the HB_S67 assembled contigs are significantly distinct compared to those in [27] (Table S2). However, the presence of small amounts of type C in one sample is interesting but does not influence the findings of this study.

Using the NGS data, both types A and B DWV Vicuna contigs were aligned across the entire genome to look for differences that correlated with deformity. Despite this unprecedented level of detail, there were no regions of the genome where all three deformed or asymptomatic samples grouped together. Neighbor-joining trees were created to examine the phylogeny of DWV variants sequenced and assembled using three regions of the DWV genome: a 4360 bp region spanning the majority of the non-structural block including Helicase and 3C protease, and the majority of the RdRp gene (Figure 3a); and two further regions, both 145 bp in length, that represent a portion of the RdRp gene (Figure 3b and Figure S3a) and the Capsid region (Figure 3c and Figure S3b). The DWV variants are split into the three master sequences; types A, B and C. The phylogenies show that within the type A clade, deformed and asymptomatic samples from the same site never share the same dominant variant of DWV. The low amount of genetic diversity within the type B clade can be attributed to the low viral load.

Although, within each location, the sequences differed for deformed and asymptomatic samples, the variation between colonies is always greater than that within each colony i.e., deformed vs. asymptomatic. Interestingly, the asymptomatic Big Island (East) sample produced an RdRp contig which contained elements of both A and B variants (Supplementary Figure S3a, HB_S21 contig 1) indicating a possible recombination between variants. This was removed prior to creating the phylogeny (Figure 3b). Another possible A–B recombinant was also observed in the deformed Oahu sample within the Helicase gene (Figure 2). However, we could not confirm the precise recombination junction site due to the lack of specific reads covering this region.