High genetic diversity of DWV in US Varroa-infested honey bee colonies reflects postbottleneck expansion

Honey bee viromes comprehensively characterized by next-generation sequencing (NGS) (S1 Table) had surprisingly high DWV genetic diversity in US Varroa-infested colonies compared with prior population surveys of the Varroa-associated DWV carried out in the UK in 2013 [13] (Fig 2, S2 Table, and S1 Fig). In particular, the levels of DWV genetic diversity in individual US honey bees with high DWV levels (Fig 2A, group 4) were significantly higher than in overtly infected UK bees with high DWV levels [13] (Fig 2A, group 2), contradicting the previously reported nearly clonal nature of DWV in Varroa-infested colonies in the Hawaiian Islands [10] and UK [13]. In fact, US honey bees with both high and low DWV levels (Fig 2A, groups 3 and 4) showed virus diversity levels similar to those observed in covertly infected UK bees with low DWV levels (Fig 2A, group 1). Significantly higher virus diversity in US honey bees compared with UK bees was also observed at the colony level when the pool of 50 bees was analyzed by NGS. Whereas the UK colony-level DWV population from a Warwickshire Varroa-infested colony sourced in July 2013 was nearly clonal (Fig 2A and 2B; “Colony UK 2013,” black arrows), the diversity of DWV in the US Varroa-infested colony sourced in October 2015 in Maryland (Fig 2A and 2B; “Colony USA 2015,” green arrow) was even higher than in the UK bees with low DWV levels, suggesting that DWV populations in the US have more genetic variants.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. The connection between virus levels in the US and UK honey bees.

(A) The columns show the average Shannon's diversity index for the NS region of DWV genomic RNA for individual honey bees from the UK (UK 2013 wild-type) and the USA (USA 2015–2017 wild-type) with low DWV or high DWV levels (less than 1% or more than 1% of DWV reads in NGS library, correspondingly) and honey bees infected with the clone-derived US DWV variants. Error bars indicate standard deviation. The groups with significantly different (P < 0.01) Shannon’s diversity index values are indicated above the bars with different letters. (B) Average Shannon's diversity index for the NS region plotted against the proportion of DWV in the NGS library. (A, B) Arrows at the right side of the graphs indicate Shannon’s diversity index values of the colony-level DWV populations in the UK and the USA were calculated for the NGS libraries for virus preparations from 50 bees. The DWV NS region corresponds to the positions 5,500–9,847 of the US DWV isolate DWV-304 (GenBank accession number MG831200). NGS libraries’ IDs (S2 Table) for the groups are shown in brackets: group 1 (UK -A1, UK-A3, UK-B1, UK-B2, UK-F3, UK-F7), group 2 (UK-E7, UK-E8, UK-E10, UK-E11, UK-INJ4, UK-INJ5, UK-INJ6), group 3 (V113_WT15_0, V31_PBS), group 4 (V52_Pupa15, V56_Pupa15, V57_Pupa15, V62_Pupa15, V114_WT15_3), group 5 (V22_304TR, V25_306TR, V16_422TR, V19_702TR, V21_703TR), colony UK 2013 (UK-VIR2), and colony USA 2015 (V99_VIROCT15). The libraries’ statistics is provided in S2 Table, NGS coverage and Shannon’s diversity profiles for the analyzed region are shown in S1 Fig, and numerical values underlying the summary graphs are provided in S1 Data. DWV, deformed wing virus; ID, identifier; NGS, next-generation sequencing; NS, nonstructural.

https://doi.org/10.1371/journal.pbio.3000502.g002

Indeed, analysis of the divergent position distribution in a virulent DWV population from a US colony that collapsed 3 months after sampling showed that it had 5.6% polymorphic nucleotides across the viral genome, in which the proportion of the alternate nucleotide exceeded 3%, generating 1.8% amino acid substitutions (Fig 1B, “USA 2015, 50 bees”; S1 Table, V99). NGS analysis also showed that this collapsed honey bee colony included exclusively DWV genotypes of type A (DWV-A), most closely related to US isolates DWV-PA [1] and DWV-Ame711 [18], and did not harbor VDV1 type, which has spread rapidly in the US in the last decade [19]. DWV polymorphism levels and distribution of diversity in individual honey bee pupae (Fig 1B, “USA 2015, individual bees”) closely matched this single-colony sample (Fig 1B, “Virus-USA 2015”), with diversity profile Pearson’s correlation coefficients ranging from 0.7453 to 0.9549 (S3 Table). Moreover, bees collected in 2010 in two apiaries in Texas and Maryland had DWV diversity profiles similar to the 2015 Maryland samples (Fig 1B, “USA 2010, average”; 2010 and 2015 Pearson’s correlation coefficients ranged from 0.5287 to 0.7432, S3 Table). Given the phylogenetic relatedness between the DWV consensus sequences from 2010 and 2015 (Fig 3 and S2A Fig, sequences labeled with suffix “-Cons”), a stable divergent DWV population, like the one we report in this study, is the most widespread, or at least very common, in Varroa-infested bees across the mainland US and has existed since at least 2010.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Phylogeny of genomic segments of DWV RNAs.

Maximum likelihood phylogenetic trees were generated for the sequences coding for the (A) LP (positions 1.1–1.7 kb), (B) structural proteins (positions 1.7–4.3 kb), and (C) major nonstructural proteins (positions 4.3–10.1 kb). Nodes connecting the cloned sequences (“Clone-“), are shown in red; stars indicate DWV genotypes with tested infectivity. Bootstrap values for 1,000 replicates are shown for groups with more than 50% bootstrap support. Scale shows genetic distance (%). DWV, deformed wing virus; LP, leader protein.

https://doi.org/10.1371/journal.pbio.3000502.g003

Phylogenetic analysis of the nearly full-length cloned DWV genomes (9.7-kb sections containing the full protein-coding sequences) produced for this study from a Varroa-infested Maryland colony (S2 Fig, sequences labeled with prefix “Clone-”) revealed that they formed a cluster that also included the source colony’s NGS consensus (S1 Table, NGS library V99_VIROCT15), consensus sequences from individual bees from the same apiary, and the US 2010 bees from Maryland and Texas (S2 Fig, sequences labeled with suffix “-Cons”). This clade, which had 98% bootstrap support, was rooted in a branch with all complete DWV-A genomes from the US and the NGS-derived DWV sequences for 2010 and 2017 US bees from this study, which had 100% bootstrap support (S2A Fig). All cloned viral genomes belonged to DWV-A type (S2B Fig), which was in good agreement with NGS analysis of the source Maryland colony (S1 Table, NGS library V99_VIROCT15). Specific genome sections revealed similar phylogenetic relationships (Fig 3), further corroborating that the US DWV-A population was generated following a strong bottleneck event and subsequent diversification from a single, or closely related, strain(s) from Europe.

Design of the reverse-genetics system for a virulent US DWV population

Stable coexistence of multiple variants within a virulent DWV population prompted questions about interactions between virus genotypes and driving forces maintaining high diversity. We used a molecular approach to investigate these interactions, which involved producing a series of full-length infectious cDNA clones of genomes of distinct isolates coexisting in a DWV population from a single colony (Fig 3 and S2 Fig, sequences labeled with a star). Together, these clones captured a significant proportion of the genetic diversity present in a typical US DWV-A population (Fig 1C, blue line). The colony used for virus preparation was Varroa infested, contained honey bees showing wing deformities consistent with high DWV levels [2,6], and did not survive the following winter season, indicating that it harbored a virulent DWV strain typical of declining US colonies. The selected clones (n = 5) and parental population (“USA-2015, 50 bees”) had highly similar distributions of divergent nucleotides (n = 286) and amino acids (n = 36) (Fig 1C and 1B). Shannon’s diversity profile for equal proportions of each clone (Fig 1C, blue line) was positively correlated (Pearson’s R = 0.804; S3 Table) with that of the parental DWV population (Fig 1B, black line). When injected into honey bee pupae, full-length DWV in vitro RNA transcripts generated from the cDNA clones were infectious. The clone-derived DWV isolates replicated to similarly high levels, about 10¹⁰–10¹¹ genome copies per bee, as wild-type RNA of the same DWV preparation used to design the clones. In addition, typical 30-nm DWV virus particles were observed in filtered tissue extracts from the in vitro transcript-injected pupae (S3 Fig) 4 days postinoculation (dpi). Significantly lower (Wilcoxon rank sum test, P < 0.01) virus levels were observed in pupae injected with PBS or with the mutant transcripts 304Δ and 306Δ, in which essential viral genes required for replication were deleted (Fig 4A). Because DWV is widespread and present at low levels, about 10⁵–10⁷ copies per insect in virtually all honey bees, including Varroa free [5,13], this virus load was detected by quantitative reverse-transcription PCR (qRT-PCR) even in control samples. Therefore, it was important to determine that the virus that replicated to high levels in the DWV transcript-injected honey bees was indeed clone derived. The identity of clone-derived DWV progeny in the injected pupae was confirmed both by using the unique restriction site introduced into the clones (Fig 4B) and by NGS analysis (S1 Table), which showed that the consensus nucleotides of the clone DWV sequences were identical to their respective parental cloned cDNA (Fig 4C), while retaining low diversity (Fig 1C, red lines; Fig 2, group 5—“USA 2015 cDNA clone-derived USA-2015”; S1 Table), thereby proving that the cloned DWV isolates were infectious.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Infectivity of clone-derived DWV isolates.

(A) Replication of DWV in honey bee pupae injected with 5 μg of in vitro RNA transcripts from full-length infectious DWV cDNA clones (304, 306, 422, 702, and 703), DWV cDNA constructs with a deletion of essential replication genes (Δ304, Δ306), 5 μg of WT DWV RNA (USA 2015 WT DWV RNA), or buffer control (PBS). The injected pupae, n = 6 per each group, were sampled at 3 dpi, and the DWV loads per individual pupa were quantified by qRT-PCR. Black dots indicate DWV load in individual honey bees. The columns show the average DWV copy numbers for each treatment ± SD. Statistically significant (P < 0.01) and ns differences between treatments groups are indicated above the bars. Numerical values underlying the summary graphs are provided in S1 Data. (B) Presence of the introduced PacI restriction site in the progeny of the clone-derived DWV genomes in the in vitro transcript-injected pupae at 3 dpi. The RT-PCR DNA products corresponding to the 5′ terminal part of the DWV RNA, untreated or digested with PacI, were separated by agarose gel electrophoresis. (C) NGS analysis of the cDNA clone-derived DWV progeny from the in vitro transcript-injected and clone-derived extract-injected honey bee pupae. Shown is the percentage of similarity between the DWV consensus sequences produced from NGS libraries (sample ID indicated in brackets) compared with their respective parental DWV cDNA sequences. (D) NGS analysis of the pupae infected with divergent wild-type DWV and clone-derived DWV (clone 304) immediately after injection (0 dpi) and at 3 dpi. NGS library ID is contained within parentheses. dpi, days postinoculation; DWV, deformed wing virus; ID, identifier; IRES, internal ribosome entry site; NGS, next-generation sequencing; ns, nonsignificant; Prot, viral protease; qRT-PCR, quantitative RT-PCR; RdRpol, RNA-dependent RNA polymerase; RT-PCR, reverse-transcription PCR; VP, structural viral protein; Vpg, genome-linked protein; WT, wild-type.

https://doi.org/10.1371/journal.pbio.3000502.g004

We then investigated replication dynamics of the recovered clone-derived virus isolates. After 48 and 72 hours, DWV levels were the same in pupae injected with 10⁷ genome copies of either clone-derived or wild-type DWV (Fig 5A), further demonstrating that individual DWV genomes replicate to the same level as the parental divergent DWV population (Fig 5A, Wilcoxon rank sum test, P > 0.01). Consensus DWV sequences at 3 dpi were identical to their corresponding cDNA clones in all but one case, likely because of a background DWV infection in the injected pupae (Fig 4C). DWV levels at 3 dpi were significantly higher (Fig 5B, Wilcoxon rank sum test, P < 0.0001) than 0-dpi pupae sampled immediately after injection with the virus, which had the same levels as PBS-injected controls (Fig 5B). The recipient honey bee pupae had low-titer background infection, which caused DWV diversity to be high in a DWV-304 0-dpi pupa injected with clone-derived DWV (Fig 4D, black lines—V39), but the clone-derived progeny had low, nearly clonal genetic diversity at 3 dpi (Fig 4D, orange lines—V42), indicating that only the hemolymph-injected clone-derived virus replicated to high levels, rather than the highly divergent background, in stark contrast with pupae injected with wild-type DWV (Fig 4D, red lines—V114). Notably, in the honey bee pupa injected with buffer (PBS) control, divergent DWV remained at a low level (Fig 4D, green lines—V31). It could be argued that the mutation rate was too low to detect accumulation of mutant variants within 3 dpi; therefore, future studies may wish to extend the postinjection analysis time frame to allow for direct observation diversification of the clone-derived genotypes.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Replication dynamics and interaction between the clone-derived DWV isolates.

(A) Average DWV RNA copies per honey bee pupa ± SD. Individual pupae copy numbers indicated by black dots. (B) Proportion of DWV reads in NGS libraries for individual, mixed, and control infections. Significant and ns differences are indicated. (C) Recombination breakpoints generated in mixed clone-derived infections. Black diamonds show the proportion of recombinant reads in individual mixed infections (left y-axis). Blue bars show average proportion of recombinant reads in a genome section (right y-axis). DWV genetic map is shown below; numerical values underlying the summary graphs are provided in S1 Data. dpi, days postinjection; DWV, deformed wing virus; extr, virus extract; LP, leader protein; NGS, next-generation sequencing; ns, nonsignificant; p.i. postinjection; WT, wild-type.

https://doi.org/10.1371/journal.pbio.3000502.g005

Mixed infection of clone-derived DWV genotypes showed a lack of synergistic effects and competitive exclusion

Divergent virus populations might have increased fitness compared with individual genotypes because of complementation between genotypes that could reinforce each other to achieve a higher replication rate [20]. The replication dynamics of the individual clone-derived DWV isolates, their mixtures, and divergent wild-type DWV were compared. When honey bees were infected by hemolymph injection with extracts containing a total of 10⁷ genome equivalents of clone-derived virus particles, there were no significant differences (Wilcoxon rank sum test) between DWV levels in pupae infected with either individual clone-derived isolates, all possible pairwise combinations of these isolates, or a mixture composed of all five isolates (Fig 5B). The DWV Shannon’s diversity profile in pupae injected with the mixture of all five clonal isolates was similar to that predicted, assuming equal proportions of each component (Fig 1C, black and blue lines, respectively; S3 Table, Pearson’s correlation coefficient P = 0.9825). Such a tight correlation between theoretical and observed profiles for mixed infections suggests that each of five coinfected genotypes replicated to similar levels in the reassembled population. NGS analysis of infections induced by pairwise clone-derived isolate mixtures (S1 Table) allowed a genome-wide view of strain success via the ratios of alternate nucleotides at the expected divergent positions (S4 Fig). Although some clones showed higher proportions than others, no complete competitive exclusion [21] was observed over 72 hours when resident genotypes of DWV were coinoculated.

Widespread recombination between the DWV isolates in mixed infections

NGS of the pairwise clone-derived isolate mixtures also allowed us to investigate recombination within DWV populations. Analysis of the changes in the proportion of divergent nucleotides along the DWV genome (S4 Fig, lower panels, orange bars) revealed widespread generation of novel DWV variants as a result of recombination events between the clone-derived genotypes. This was confirmed by an analysis for structural variants that identified recombination breakpoints between clones (S4 Fig, upper panels). Recombination sites were clustered mostly in the regions preceding the leader protein (LP) region, the LP region itself (the main 3′-proximal nonstructural block in the helicase region), and the border of the regions coding for the structural viral protein (VP) 1 and VP3 (Fig 5C). Independent evidence for a high recombination frequency between DWV genotypes came from the phylogenetic analysis of the functional sections of the cloned full-length DWV genomes (Fig 3). It showed different topologies in the phylogenetic trees, supported by high bootstrap values, for the genome sections coding for LP, the structural proteins, and the nonstructural proteins. For example, clone 304 and clone 18 have almost identical structural gene sections (Fig 3B), but their LP gene sections are in different clades (Fig 3A). Similarly, the nonstructural genes of clone 119 and clone 304 are very close (Fig 3C), but their structural genes are in separate branches (Fig 3B). Such differences in tree topologies could be explained by reshuffling of genome sections between members of the DWV population. Recombination events are known to be involved in the generation of novel variants of picorna-like viruses [22]. The fact that recombination events are often located at the borders between the functional blocks reflects the modular nature of the DWV genome. The distribution of recombination sites observed in the mixed infection experiment (Fig 5C and S4 Fig) resembles that found in previous studies. This includes positions flanking the LP region and in the 5′-proximal part of the nonstructural gene block in the recombinants between DWV-A and DWV-B [12–14,23,24], as well as the protease and RNA-dependent RNA polymerase (RdRpol) region of the nonstructural block in the recombinants between Asian DWV-A isolates [25].

High genetic heterogenicity of DWV populations as an RNA interference evasion mechanism

We found that individual DWV clones from a virulent population replicated to the same levels as did either mixtures of these clones or the parental divergent population (Fig 5B), indicating the absence of obvious synergetic effects. We suggest that diversification of DWV in honey bees and the maintenance of high genetic diversity in wild-type DWV populations might be driven by host antiviral RNA interference (RNAi) defenses [26] and maintained as a way to evade sequence-specific RNAi similar to West Nile virus (WNV) in mosquitoes and Powassan virus in ticks [27–29]. These studies showed that rare genetic variants of these viruses avoid control by RNAi targeting major variants because of the lack of complementarity between the minor strains and the guiding small interfering RNAs (siRNAs) derived from the major viral strains [27–29]. Analysis of the distribution of polymorphic sites in DWV populations provided several findings in support of this suggestion.

We found that most of the nucleotide changes in sampled DWV populations were silent, with a mean ratio of synonymous (dS) to nonsynonymous (dN) substitutions (dS/dN) in the tested infectious clones reaching 24.7, indicating strong purifying selection (S4 Table). Accumulation of a high number of nucleotide changes throughout DWV genomes, constrained by the need to maintain the same coding capacity, suggests that the virus is forced to explore sequence space to maintain its high diversity. This is consistent with the hypothesis that diversification of DWV RNA is driven mostly by selection of novel sequence variants of DWV RNA capable of evading specific RNAi targeting because of mismatches with siRNAs. Analysis of the NGS libraries confirmed that polymorphic sites are present throughout DWV genomes (Fig 1B), in agreement with previous reports that the entire DWV genome could be targeted by RNAi [13,24] and therefore subjected to diversification. Still, a significantly higher density of polymorphic positions was consistently observed in the LP region and in the 3′-proximal part of DWV coding for the nonstructural block (S5 Fig and S4 Table). Additionally, the cloned infectious DWV isolates showed a subset of codons (n = 31) with dS/dN < 1, which are potentially subjected to diversifying positive and/or relaxed selection. Interestingly, 12 of these putative positively selected codons were located in the LP region, a 10-fold higher density than in the rest of the genome (S5 Fig and S4 Table).

To estimate the potential of US DWV populations to evade specific RNAi established to a single given DWV strain, we analyzed DWV NGS libraries from a typical colony, a single honey bee sampled in 2015, and a pool of eight honey bees sampled in 2010 (S1 Table; Libraries B2_PBS, V57, and V99) for sequence variants capable of avoiding specific RNAi targeting by recording the number of alternate nucleotides (occurring above 1% and 10% levels of the read coverage) in a sliding 22-nt window (Fig 6)—i.e., the size of a potential siRNA target in the honey bee [13,24]. Within the colony-level DWV population, most of the genome, 81.80%, showed at least one sequence variant present at a frequency of 1% or higher, and 38.81% of the genome showed variants present at 10% or higher. Notably, the polymorphic 22-nt windows covered a significant proportion of the viral genome in the DWV populations from the individual bee sampled in 2015, with 35.53% and 41.31% of genomic windows showing sequence variants present at 1% or higher and 10% or higher, respectively. Similarly, the DWV population from a pool of 2010 US honey bees showed 57.42% and 33.01% divergent coverage. Clearly, divergent DWV populations harbor a pool of readily available genetic variants, some of which may possess reduced vulnerability to specific RNAi. Therefore, we propose that widespread recombination events between divergent coexisting DWV genotypes in a population contribute to the generation of variants less effectively targeted by established RNAi, thereby helping DWV to evade this sequence-specific antivirus response.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Potential contribution of DWV population diversity to evasion of sequence-specific RNAi targeting.

Distribution of the number of alternate positions in 22-nt-long potential siRNA target windows for the polymorphic sites with more than 1% (light blue graphs) and 10% (dark blue graphs) of NGS coverage. A DWV genetic map is shown above the graphs. Summaries of the DWV genome coverage (%) with different number of mismatches in a 22-nt sliding window for each of the NGS libraries is shown at the right of the graph. DWV, deformed wing virus; IRES, internal ribosome entry site; LP, leader protein; NGS, next-generation sequencing; Prot, viral protease; RdRpol, RNA-dependent RNA polymerase; RNAi, RNA interference; small interfering RNA; VP, structural viral protein.

https://doi.org/10.1371/journal.pbio.3000502.g006

Natural differences between DWV genotypes might significantly affect efficiency of antiviral RNAi targeting

We sought to experimentally test whether natural variation between DWV genotypes coexisting in a single population was sufficient to significantly impact the effectiveness of antiviral RNAi and, indeed, could act as an RNAi evasion mechanism. It was shown previously that antiviral RNAi in honey bees could be induced by introduction of double-stranded RNA (dsRNA) molecules orally [30]. To do this, we produced dsRNA specific to two coexisting isolates, US DWV-304 and DWV-422. The targeted 283-nt region of DWV genomic RNA (positions 1,242–1,524 of the GenBank accession MG831200) had 5.3% divergent nucleotides between isolates DWV-304 and DWV-422 (S1 Text), which mostly resulted in one mismatch in any given 22-nt sliding window (Fig 7A). In two independent experiments (design summarized in Fig 7D), newly emerged adult honey bees from colonies with low Varroa infestation levels were fed with 1 μg of the dsRNA (ds304 or ds422) in 5 μL of 50% sucrose either with or without the virus extract containing 10⁸ genome equivalents of the clone-derived DWV-304 or DWV-422. In order to control for the impact of non-sequence-specific triggering of antiviral response by dsRNA, a phenomenon reported in honey bees [31,32], the two dsRNAs used were 94.7% identical at the nucleotide level and therefore likely to have similar nonspecific antiviral triggering capacity. The DWV-304- and DWV-422-derived dsRNAs therefore acted as controls against each other, making it possible to investigate the impact of the 5.3% nt variation between them on DWV infection dynamics. The bees were maintained for 7 days prior to sampling and quantification of DWV loads in individual insects. The identity of clone-derived DWV variants was confirmed by sequencing the RT-PCR fragments corresponding to the LP region.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Impact of specificity of RNAi targeting on outcome of DWV infections.

(A) Design of the 283-nt dsRNA specific to the DWV isolates 304 and 422 (ds304 and ds422) corresponding to the regions 1,242–1,524 of DWV genome (GenBank accession MG831200). (B, C) DWV loads in the individual adult honey bees 7 days postfeeding with 10⁸ clone-derived virus particles and 1 μg of the DWV-304- or DWV-422-specific dsRNA. Columns show average copy number per treatment groups ± SD; black dots indicate DWV copy numbers in individual honey bees. The groups with significantly different (P < 0.01, ANOVA) DWV copy numbers are indicated above the bars with different letters. The red lines indicate the threshold level of 10⁹ DWV genome copies per honey bee, which was used to categorize DWV levels as low or high. Numerical values underlying the summary graphs are provided in S1 Data. (D) Analysis of the DWV infection outcomes in two oral inoculation experiments with the strain-specific dsRNAs. Circles indicate experimental groups; numbers of the bees in each group with high (>10⁹ copies) or low DWV (below 10⁹ DWV copies) levels of DWV are shown on the top or bottom of the circles. A 2 × 2 contingency table analysis was used to test significance of the differences in proportions of the bees with low or high DWV levels. dsRNA, double-stranded RNA; DWV, deformed wing virus; IRES, internal ribosome entry site; LP, leader protein; NS, not significant; Prot, viral protease; RdRpol, RNA-dependent RNA polymerase; RNAi, RNA interference; VP, structural viral protein; Vpg, genome-linked protein.

https://doi.org/10.1371/journal.pbio.3000502.g007

The outcome of DWV-304 infection was significantly influenced by the specificity of the dsRNA. There were significantly higher levels of DWV in the DWV-304-infected bees, which received ds422 (Fig 7B and 7C, group DWV304/ds422), compared with the group that received ds304 (Fig 7B and 7C, group DWV-304/ds304), P < 0.05, ANOVA. This was observed in two independent experiments using honey bees from different colonies with different wild-type background DWV levels (Fig 7C and 7D, mock infected). Similar variation in “background” wild-type DWV loads was reported previously [33]. Virus levels across individual bees showed a bimodal distribution, which allowed us to categorize DWV loads as low or high, with a 10⁹-copy threshold, and carry out a 2 × 2 contingency table analysis. We found that, in the case of bees infected with DWV-304, a significantly reduced number of bees developed a high-level DWV infection in the group “DWV-304/ds304” (which received ds304 perfectly matching virus isolate) compared with the group “DWV-304/ds422” (which received ds422 mismatching the inoculated virus), P = 0.000185 (Fig 7D). Notably, the proportion of bees with high DWV levels in the group “DWV-304/ds304” were not significantly different from that in control group “Mock-infected/ds304” (Fig 7D). Similarly, ds422 suppressed development of high DWV levels more efficiently in the case of inoculation with DWV-422 compared with DWV-304, P = 0.0115 (Fig 7D, groups “DWV-304/ds422” and “DWV-422/ds422”), although it cannot be excluded that, in this case, differences of infectivity of these virus isolates might play a role in development of the observed pattern of infection. These results suggest that natural variation between DWV genotypes within a single population (e.g., DWV-304 and DWV-422) might significantly reduce the efficiency of RNAi targeting.

Honey bees

Varroa-infested honey bee colonies maintained by the USDA-ARS Bee Research Laboratory (Beltsville, Maryland, USA), which had high Varroa infestation (15%–20% of pupae infested with mites), and bees showing wing deformities consistent with high DWV levels were used as the source of mite-exposed honey bee pupae for NGS analysis and isolation of DWV particles in October 2015. Unparasitized honey bee pupae from colonies with low Varroa infestation (with below 1% of mite-parasitized pupae) were used in the injection experiments. Newly emerged pupae from the Maryland colonies with low Varroa infestation were used in the oral infection experiments. A colony from Texas sourced in January 2017 was used in the in vitro RNA transcript and wild-type RNA injection experiment. A colony from Florida sourced in February 2017 was used in the experiment involving injection of virus extracts, both clone-derived and wild-type.

Sequence analysis

The NGS libraries for this study were produced by Illumina HiSeq2500, with each containing 12–14 million paired-end reads (S1 Table; NCBI BioProject PRJNA431793, Sequence Read Archive study SRP135682). These libraries, along with the previously published UK NGS libraries [13] (S2 Table 2), were analyzed as described previously [19]. Briefly, Illumina reads were trimmed to remove adapters, contaminant sequences, and low-quality bases. Cleaned data from each sample library were individually aligned to reference sequences of three distinct strains of DWV-like viruses, DWV-A (GenBank accession GU109335), VDV1 (GenBank accession AY251269), and DWV-C [3], using Bowtie2. Nucleotide counts, coverage, and Shannon’s diversity index estimates for each nucleotide position were calculated from SAMtools mpileup output. Although Illumina HiSeq technology has a very low per-nucleotide substitution error rate (0.08% and 0.12% for the first and second reads, correspondingly [60]), to reduce sequencing chemistry errors, the Shannon’s diversity index estimate was further corrected for potential errors using the approach described previously [61]. In brief, the error model was calibrated using data from an Illumina library from a bee injected with a cloned virus of known sequence. Data from positions 5,230 to 6,600 of the DWV genomic RNA, representing a low-diversity region of the virus, were used to calculate an alphahat of 9.88 × 10⁻⁵. To quantitatively assess the degree of similarity between different DWV genetic diversity profiles, we calculated the Pearson’s correlation between the Shannon’s diversity index profiles averaged for 100-nt windows. For detection of recombination breakpoints in the honey bee pupae injected with the mixtures of the clone-derived DWV genotypes, cleaned data from each sample NGS library were individually aligned to the VDV1 and DWV references using SpeedSeq (version 0.1.2). Structural variants were detected using the lumpy smooth script from the LUMPY package (version 0.2.13), and additional genotype and coverage metrics for each variant were calculated using SVTyper, a script within the SpeedSeq package. Breakpoint positions representing a recombination between two different viral isolate sequences were considered for further analysis if they had evidence from more than 10 supporting events (split or discordant reads) and exceeded 0.02% of NGS coverage (the level was chosen to ensure that the reads corresponded to viable—i.e., replicating—recombinant genomes, rather than aberrant variants). Phylogenetic analysis was carried out by aligning the viral nucleotide sequences using CLUSTAL and producing maximum likelihood and neighbor-joining trees, which were bootstrapped using 1,000 replicates with the RAcML and PHYLIP packages, respectively. dS and dN variations in viral genomes were identified using all known full-length North and South American DWV-A sequences (S2 Fig) using tools provided by the Los Alamos National Laboratory’s HIV Databases (https://www.hiv.lanl.gov [62]); putative positive diversifying selection was assumed when dS/dN < 1 (S4 Table).

Analysis of DWV diversity in Hawaiian (Oahu Island) apiary-level samples collected in 2015 as a part of the USDA-APHIS honey bee survey [35] included qRT-PCR amplification of the DWV cDNA fragments corresponding to the LP and the RdRp regions, using the oligonucleotide primers “DWV-LP-For” and “DWV-LP-Rev”, “DWV-POL-For” and “DWV-POL-Rev” (S5 Table), respectively, and direct sequencing of the RT-PCR products to identify polymorphic sites in the electropherograms. The RT-PCR fragments corresponding to the LP region were also cloned into a pGem-T easy vector (Promega) and Sanger sequenced.

DWV cDNA clones and clone-derived viruses

Viral particles were purified by caesium chloride gradient centrifugation from a pool of 50 pupae out of a single Varroa-infested colony at the Bee Research Lab in Beltsville, Maryland, USA, in October 2015 as described previously [12]. Viral RNA was extracted from the virus preparation using the RNeasy kit (Qiagen). The first cDNA strand corresponding to the entire 10-kb-long genomic RNA was produced using Superscript III Reverse Transcriptase (Invitrogen) and the oligonucleotide primer “DWV-PmeI-A27Rev” complementary to the 3′ terminus of the DWV genomic RNA, which was preceded by the 27-nt oligo dT and the PmeI restriction site (S5 Table). The cDNA was used as a template to amplify a nearly full-length, 9.7-kb DWV RT-PCR fragment containing the 686-nt section of 5′ untranslated region adjacent to the viral ORF, the full 8,682-nt protein-coding sequence, and complete 341-nt 3′ untranslated region using the oligonucleotide primers “DWV-PmeI-A27Rev” and “DWV-USA-PacI462F,” identical to positions 472–502, preceded by the PacI restriction site (S5 Table), and the proof-reading thermostable DNA polymerase Phusion (New England Biolabs). The library was cloned into the pCR4Blunt-TOPO vector (Invitrogen) according to the manufacturer’s instructions. Ten clones containing 9.7-kb viral cDNA inserts were Sanger sequenced (GenBank accession numbers MG831200–MG831204 and MH069503–MH069507). A set of five divergent DWV cDNA constructs (GenBank accession numbers MG831200–MG831204) was selected to construct a set of full-length infectious DWV cDNA clones capturing a significant proportion of colony-level diversity. This was accomplished by inserting at the position corresponding to the 5′ part of the DWV cDNA the synthetic fragment “T7-Ribo-5′-DWV” (S5 Table) containing the T7 RNA polymerase promoter sequence, the 5′ ribozyme sequence [63], and the 462-nt-long 5′ terminal part of the DWV genome, using the NotI and PacI restriction sites to create a series of five plasmids with different full-length cDNA inserts. The 462-nt-long 5′ terminal sequence was designed according to the DWV USA NGS data obtained in this study and included an extended 5′ terminal sequence (as previously [64] and in the GenBank accession number KT215904). To produce control constructs p304Δ and p306Δ containing noninfectious DWV cDNA, the regions 7,128–7,998 and 7,128–8,347 coding for the viral protease region were deleted from clones p304 and p306 (GenBank accession numbers MG831200 and MG831201), respectively, creating a frameshift in the viral ORF that abolished translation of the viral RdRpol. Plasmid constructs were linearized by using the PmeI restriction site located at the 3′ end of DWV cDNA, downstream of the A27 sequence, to produce the templates for a full-length 10.2-kb transcript. The in vitro transcripts were produced using HiScribe T7 High Yield RNA Synthesis Kit (New England Biolabs) according to the manufacturer’s instructions. After 3 hours of incubation at +37 °C, the DNA templates were digested using TURBO DNase (Life Technologies). The RNA transcripts were purified using the RNeasy mini kit (Qiagen), eluted with RNAse-free sterile water, and stored at −80 °C prior to use.

RNA transcripts were injected into honey bee pupae without Varroa mite feeding exposure. These pupae were extracted from brood frames with low Varroa infestation (approximately one mite-infested cell per 500) from a Texas apiary in December 2016. The pupae at pink eye stage were injected intra-abdominally into the hemolymph with 5 μg (13.2 × 10¹² copies) of the in vitro transcripts or wild-type DWV RNA in 10 μL of PBS or 10 μL of PBS control using syringes with a 0.3-mm needle G31 (BD Micro-Fine), essentially as described previously [13], and were incubated at +33 °C, relative humidity 80%, for 3 days to allow replication of the virus prior to RNA extraction and sequencing. It was necessary to inject high copy numbers of naked viral RNA to reliably initiate clone-derived infection because, as expected, infectivity of naked RNA was very low compared with an encapsidated RNA. Notably, the injected RNA, either in vitro transcript or extracted from virus particles, was almost entirely degraded in the hemolymph after 3 days; qRT-PCR quantification of DWV in the pupae injected with 9.3 × 10¹² copies of the control nonreplicating DWV in vitro transcripts p304Δ and p306Δ had between 1.3 × 10⁵ and 1.7 × 10⁷ genome copies detected (Fig 4A). For comprehensive characterization of the transcript-derived progeny, NGS libraries were generated from individual pupae each containing 12–17 million paired-end 150-nt reads, with the share of DWV reads ranging between 43% and 63% (S1 Table 1).

For preparation of the clonal DWV extracts containing infectious DWV virus particles, individual transcript-injected pupae were homogenized with 2 mL of PBS at 3 dpi. For each individual pupal extract, 1 mL was used for preparation of the DWV inocula, which included clarification by centrifugation at 3,000g for 5 minutes and filtration through a 0.22 μm nylon filter (Millipore). The DWV concentration in the extracts was quantified by qRT-PCR, and the extracts were stored at −80 °C prior to use. To confirm the identity of the clonal DWV in the extracts with their respective parental DWV cDNA clones, the DWV RNA region, positions 30–1,266 nt, was amplified by RT-PCR, sequenced, and also digested with PacI restriction enzyme.

Pink eye–stage honey bee pupae, with no exposure to Varroa mite feeding, were extracted from a Florida brood frame with low Varroa infestation (1 infested cell per 500) in January 2017. Pupae were intra-abdominally injected into the hemolymph to introduce either the filtered DWV extracts (10⁷ genome copies in 10 μL of PBS), which included the clone-derived recovered DWV variants that originated from transcript-injected pupae, the wild-type DWV, or buffer control (PBS). We sampled time series from “time 0” (immediately after injection) through to 3 dpi, and the DWV RNA copy numbers were quantified in individual bees by qRT-PCR as described previously [38].

Assessment of the effect of strain-specific dsRNA on DWV replication

dsRNAs specific to the DWV clones were produced in vitro with T7 RNA polymerase (HiScribe T7 RNA polymerase, New England Biolabs) using as the templates PCR fragments corresponding to the positions 1,242–1,524 of the clones DWV-304 and DWV-422 (GenBank accession numbers MG831200 and MG831202, respectively) with T7 promoter sequences at both the 5′ and 3′ termini. These PCR fragments were amplified using the US DWV cDNA clones p304 and p422 as the templates and the primers shown in S5 Table. Oral infection included feeding individual newly emerged bees with 5 μL of 50% sucrose containing 1 μg of purified dsRNA (ds304 or ds422) and 10⁸ clone-derived virus particles, DWV-304 or DWV-422, or no virus in the case of control groups (Fig 7D). The honey bees were sourced from Maryland colonies with low Varroa infestation levels (less than 1% of pupae being infested with the Varroa mites) and were starved for an additional 3 hours before feeding. The bees were kept for 7 days essentially as described in [65] in cages at +33 °C, 75% relative humidity, with ad libitum sources of water and 1:1 sucrose/water syrup (w/v), prior to RNA extraction.

Data deposition

Raw sequence data and consensus DWV sequences have been deposited at DDBJ/EMBL/GenBank within BioProject PRJNA431793 under the Sequence Read Archive (SRA) study SRP135682 (accessions SRR6833910–SRR6833954) and Transcriptome Shotgun Assembly (TSA) study GGSE00000000 (accessions GGSE01000001–GGSE01000045), respectively. Cloned viral cDNA sequences are available in GenBank with accession numbers MG831200–MG831204, MH069503–MH069507, and MH594118–MH594121.