Skip to main content
Download PDF

Vector-virus interaction affects viral loads and co-occurrence

BMC Biology volume 20, Article number: 284 (2022) Cite this article

Abstract

Background

Vector-borne viral diseases threaten human and wildlife worldwide. Vectors are often viewed as a passive syringe injecting the virus. However, to survive, replicate and spread, viruses must manipulate vector biology. While most vector-borne viral research focuses on vectors transmitting a single virus, in reality, vectors often carry diverse viruses. Yet how viruses affect the vectors remains poorly understood. Here, we focused on the varroa mite (Varroa destructor), an emergent parasite that can carry over 20 honey bee viruses, and has been responsible for colony collapses worldwide, as well as changes in global viral populations. Co-evolution of the varroa and the viral community makes it possible to investigate whether viruses affect vector gene expression and whether these interactions affect viral epidemiology.

Results

Using a large set of available varroa transcriptomes, we identified how abundances of individual viruses affect the vector’s transcriptional network. We found no evidence of competition between viruses, but rather that some virus abundances are positively correlated. Furthermore, viruses that are found together interact with the vector’s gene co-expression modules in similar ways, suggesting that interactions with the vector affect viral epidemiology. We experimentally validated this observation by silencing candidate genes using RNAi and found that the reduction in varroa gene expression was accompanied by a change in viral load.

Conclusions

Combined, the meta-transcriptomic analysis and experimental results shed light on the mechanism by which viruses interact with each other and with their vector to shape the disease course.

Background

Vector-borne viral diseases are transmitted by infected vectors, typically arthropods such as mosquitoes, ticks, and mites [1, 2]. They drive species evolution in both natural and managed ecosystems and threaten a wide range of organisms. Most disease-causing plant viruses depend on vectors for spread and survival [3, 4], but animals also suffer from vector-borne viruses [5]. For humans, these include Dengue, Chikungunya, Chagas disease, Japanese encephalitis, Zika, and yellow fever, leading to hundreds of thousands of deaths worldwide, especially in developing countries [6, 7]. Moreover, ongoing globalization and climate change are expected to increase their outbreak frequency [8, 9]. The most efficient and sustainable measure of coping with emerging vector-borne diseases is to control their vectors [10].

Vectors are often viewed as merely a passive syringe injecting the virus. However, to promote replication and transmission viruses may regulate the vector’s immune system and behavior, which requires interaction with the vector. Indeed, studies showed that viral infection can alter the vector feeding behavior, fecundity, longevity, and survival [11,12,13,14]. The molecular mechanism underlying these interactions was studied in both cell cultures [15,16,17] and in in vivo experiments [18, 19], revealing vector immune-related genes whose expression is regulated following viral infection and shows that a crosstalk between vector genes and the virus is imminent for successful viral replication and transmission [20].

Most of the studies have investigated the interaction of a vector with a single virus [20], and a few considered up to two co-occurring viruses [21, 22]. However, vectors can usually carry diverse viruses that may interact with each other [23,24,25,26]. Studies on hosts infected by several viruses have shown that coinfection affects viral traits such as virulence and transmission, either by direct virus-virus interaction or indirectly by competition or cooperation [27,28,29]. Therefore, multi-infection may have a profound effect on viral evolution, diversity, and pathogenicity [30, 31].

In contrast, the effect of multiple infections on vectors has received far less attention, although is expected to have a considerable effect on the disease epidemiology [25, 26]. Therefore, virus-virus interaction is one of the main challenges in virology today [32]. We can expect that as in the case of a host co-infected by several viruses, multi-infection will also have an effect on the vector-virus interaction. Yet, we cannot simply extrapolate from host-virus studies to vector-viruses, as the two are undergoing different evolutionary processes: while host-pathogen interactions are antagonistic by definition, vector-virus interactions are less definite and can fall anywhere on the continuum between antagonistic and mutualistic [33]. Therefore, the interaction between vector-borne viruses and how they affect their relationship with the vector remains largely unknown.

Varroa destructor (hereafter “varroa”) is a parasitic mite which vectors honeybee pathogenic viruses and is routinely co-infected by multiple viruses [34]. The introduction of varroa to the European honey bee (Apis mellifera) has dramatically changed the honeybee viral landscape, leading to worldwide colony collapses threatening global food security [35]. In varroa-free colonies, the viral diversity is high, while viral loads generally remain low [36]. Normally, 1–2 years following varroa introduction to a naive colony, virus diversity and load shifts greatly, giving rise to high levels of a few viral strains as the colony is dwindling to its final collapse [37]. Despite these observations, not much attention has been given to varroa-virus interaction, and how it may explain the varroa key role as a vector of honeybee viral disease. In addition, varroa carries over 20 diverse viruses commonly co-occurring (Fig. 1), including bee-pathogenic viruses, of which some are highly associated with varroa, while others are of unknown pathogenicity to either varroa or bee [38]. Therefore, the varroa-virus system is a unique opportunity to investigate how different viruses interact with each other and with their vector.

Fig. 1
figure 1

Most of the known bee viruses are present in both honey bee and mite, but only a few replicates in the mite, or shown to be varroa-vectored. Virus presence in bees and varroa, and its ability to replicate and vectored by the mite between bees are marked by: (+) confirmed; (--) not demonstrated; (S) suspected; (?) unknown. A star mark (*) next to the virus abbreviated name indicates that this virus was not detected in any of the varroa libraries in the current study. The viruses' full names: Deformed wing virus type a (DVWa), Deformed wing virus type b (DVWb), Deformed wing virus type c (DVWc), Slow bee paralysis virus (SBPV), Sacbrood virus (SBV), Varroa destructor virus 2 (VDV2), Israel acute paralysis virus of bees (IAPV), Kashmir bee virus (KBV), Acute bee paralysis virus (ABPV), Black queen cell virus (BQCV), Bee Macula-like virus (BMV), Varroa Tymo-like virus (VTLV), Apis flavivirus (AFV), Lake Sinai virus (LSV), Varroa destructor virus 3 (VDV3), Varroa destructor virus 5 (VDV5), Chronic bee paralysis virus (CBPV), Apis mellifera nora virus 1 (ANV), Apis mellifera rhabdovirus-1 (ARV-1), Apis mellifera rhabdovirus-2 (ARV-2), Varroa orthomyxovirus-1 (VOV-1), Varroa destructor virus 4 (VDV4), Varroa mite associated genomovirus 1 isolate VPVL_36 (VPVL_36), Apis mellifera Filamentous virus (AmFV), Varroa mite associated virus 1 isolate VPVL_46 (VPVL_46). More details on each virus sequence and references can be found in Additional file 1.

Full size image

Here we studied the interaction of diverse viruses with varroa’s transcriptional network, hypothesizing that different viruses have distinct effects. We explored a large set of RNAseq data from the vector aiming (1) to determine if viruses compete or cooperate in their vector and (2) to test whether interactions between individual viruses and the vector correlate with viral epidemiological traits, such as viral load in the vector, and co-occurrence with other viruses. We found no evidence of competition between viruses, but positive correlations in the abundance of specific viruses. Furthermore, we found a strong correlation between viral co-occurrence and the manner in which they interact with the vector’s genes, and have experimentally validated some of these predictions. These results show that multi-viral infection can occur not only in the host, but also in the vector, and that these interactions have implications on the virus relationship with its vector and therefore on the disease course.

Results

The viral landscape is heterogeneous across varroa libraries

We found that each of the final 66 varroa RNAseq libraries contains at least three types of viruses (Figs. 1 and 2). The most prevalent viruses belong to the Iflaviridae family, including the bee pathogenic viruses DWVa and DWVb, and the varroa-specific virus, VDV2. On the other hand, three viruses whose presence in varroa was reported before were not detected in any of the libraries (CBPV, VPVL_46, and VPVL_36) (Figs. 1 and 2). Viral load homogeneity across libraries varied greatly between viruses: a few viruses were extremely variable (e.g., DWVa was not detected in a few libraries while reaching < 400,000 TPMs in others), while other viruses (such as ARV-1, ARV-2, and VDV2), were found in somewhat similar loads in most libraries. Interestingly, VDV2 and ARV-1 are the only viruses that were present in all varroa libraries.

Fig. 2
figure 2

Viral load is diverse across the different varroa libraries. Members of the Iflavirus family are the most prevalent, yet while some viruses are homogenous across libraries (VDV2), others are highly diverse (e.g., DWVa and DWVb). Values are log10 transformed of the reads’ TPM (transcript per million). Zero values are marked in grey (i.e., none of the reads in this library mapped to this virus). The viruses’ names are abbreviated as described in the caption of Fig. 1

Full size image

No evidence for competition between the different viruses

Multi-infection by several viruses raises the obvious possibility of interactions between them. Interestingly, all significant correlations between viral loads were positive (Pearson correlation followed by Benjamini-Hochberg FDR-correction, P < 0.1) (Fig. 3a).

Fig. 3
figure 3

Intra-viral interactions can predict virus-vector interactions. a Correlation between viral loads. b Correlation between viral loads and varroa modules (eigengenes). c Correlation model between viral load correlations (a), and the distance-matrix of the module-virus correlations (b). In a and b, viruses and modules are ordered according to hierarchical clustering; P values of the Pearson coefficient in a and b are adjusted according to FDR-correction; correlation significance marked by (*) 0.1 < P < 0.01; (**) P < 0.01. For analysis in c, the Mantel test for correlation between two matrices was conducted using 1000 permutations. All viruses’ loads are transformed by log10 of the TPM, transcripts per million.

Full size image

Varroa modules interact with viral loads

Gene network analysis of the 66 varroa RNAseq libraries clustered the 10,247 varroa genes into 12 modules, each module containing co-expressed genes. The modules were numbered from the largest (module number one which contains 4375 co-expressed genes) to the smallest (module number 12 which contains 37 genes) (Additional file 2). We found significant interactions between specific modules and viruses when correlating module eigengenes to viral loads (transcripts per million, TPM) (Pearson correlation followed by Benjamini-Hochberg FDR-correction, P < 0.1, Fig. 3b). The interaction direction and strength, represented by the correlation coefficient, varied between virus-module pairs. While VDV2, VDV4, ARV-2, and ARV-1 show positive interaction with some modules (modules 2, 3, 6, and 9), ‘bee-viruses’ (DWVc and DWVa) show negative interaction with a different module (modules 1 and 9). Interestingly, the only virus that shows both negative and positive interactions with vector’s modules is the bee-pathogenic virus, DWVa.

GO terms of interacting modules

The significantly interacting modules consist of genes that are enriched for regulatory-related GO terms such as “Regulation of gene expression” (module 7), “Regulation of metabolic processes” (modules 6 and 9), and immune response-related GO terms, such as “regulation of immunoglobulin secretion” GO:0051023 (module 6). In addition, modules 3 and 6 are enriched for viral and symbiont-related GO terms, such as “viral process” GO:0016032 and “modulation by virus of host process” GO:0019054, and “modulation by symbiont of host cellular process,” GO:0044068. For the list of all significant GO terms in the significantly interacting modules, please see Additional file 3. Of the 17 RNAi-gene homologs recently identified in varroa [39], we found that 13 genes belong to modules interacting with viruses: module 1 (negatively interacting with DWVa and DWVc), module 2 (positively interacting with ARV-1), and module 3 (positively interacting with ARV-2 and ARV-1) (Additional file 4).

Viruses that are found together interact with the vector’s gene co-expression modules in similar ways

Among the module-virus interactions (Fig. 3b), we can detect viruses that share a similar interaction with the same module. For example, VDV2 and VDV4 abundance positively correlated with modules 6 and 9. Interestingly, abundances of these two viruses also significantly correlate with each other (Fig. 3a). Similarly, ARV-2 and ARV-1 also have positively correlated abundances and positive interactions with module 3 (Fig. 3a and b). Likewise, DWVa and DWVc are positively correlated in abundance and are negatively correlated to modules 1 and 9. We evaluated if this pattern between virus-virus interaction and virus-module interaction can be generalized across all virus-virus-module interactions in our data set. Indeed, we found a significant positive correlation between the two distance matrices (Mantel test for correlation between two distance-matrices [40], Mantel statistic r = 0.44, p < 0.001) (Fig. 3c). In other words, there is a correlation between viral co-occurrence and the manner in which they interact with the vector.

Validating varroa-virus interaction: gene silencing is accompanied by a change in viral load

To experimentally validate the varroa-virus interaction as predicted by the gene-network analysis, we silenced the mite’s genes using RNAi and tested its viral load, compared to the control (GFP-treated mites). For silencing, we selected candidate genes based on both the gene-network analysis and on the literature. From module 7, which interacts with the bee-pathogenic virus DWVa, and contains GO terms of “Regulation of gene expression,” we selected five genes which are highly connected to other genes in the module. In addition, these genes also possessed a relevant annotation, based on literature survey and/or the presence of conserved domain in the predicted coded protein (e.g. immune response-related domains, and genes that were previously reported to interact with host/vector) (for the list of silenced genes, their module-membership and annotation, see Additional file 5). Four of the five tested genes were successfully silenced, i.e., for these genes, a significant reduction in relative gene expression was measured in mites treated with dsRNA, compared to control mites (P < 0.05, Wilcoxon signed-ranks test, followed by FDR-correction, Additional file 6). We then compared the relative viral loads of DWVa, VDV2, and ARV-2 using quantitative PCR (qPCR). DWVa was predicted to negatively interact with the genes in module 7, while VDV2 and ARV-2 were not predicted to interact with the module (Fig. 3b). Of all silenced genes, only mites that were treated with dsRNA of Cuticle protein type 8-like gene (short name: CuP8, accession: LOC111248360) showed a significant change in viral load, compared to control mites. Mites with decreased expression of CuP8 had lower VDV2 and ARV-2 viral loads, compared to control mites (P = 0.02, Wilcoxon signed-ranks test, followed by FDR-correction), while DWVa viral load has also decreased, but not significantly (P = 0.2, Wilcoxon signed-ranks test, followed by FDR-correction) (Fig. 4a and b).

Fig. 4
figure 4

Validating varroa-virus interaction using gene-silencing (RNAi). a Relative expression of CuP8 gene in control (GFP-dsRNA, n=9) and silenced mites (treated with CuP8-dsRNA, n=13). b Relative viral load of VDV2, ARV-2, and DWVa in control and silenced mites. Significant difference between control and silenced mites was tested using the Wilcoxon signed-rank test, followed by FDR-correction. The boxplot shows the interquartile of the data values, the inner thicker line of the box represents the median value, and the dots are potential outliers

Full size image

The dsRNA soaking treatment did not affect mite survival, compared to the control GFP-dsRNA treated mites (Fisher exact test for goodness of fit, p < 0.05) (Additional file 7).

Discussion

Vector-borne viruses rely on another organism for transmission. As the vector plays a key role in viral fitness, we hypothesized that different viruses have distinct effects on a vector, observable as changes in gene expression. This was indeed true for viruses associated with varroa mites. Namely, loads of individual viral species were associated with specific changes to the vector’s gene expression network (Fig. 3b). Interestingly, co-occurring viruses affected the vector gene expression network in similar ways (Fig. 3c), suggesting a link between viral epidemiological traits (their titer and co-occurrence) and its relationship with the vector. We further propose that vector-virus interactions are evolving rapidly, as we found that closely related viruses have distinct effects on the vector’s gene expression. Interestingly, we found no evidence of viral competition within the vector. On the contrary, abundances of some viruses are positively correlated in co-infected vectors, suggesting a potential for viral cooperation. These complex dynamics underscore the role of vector-virus interactions for viral fitness.

Viruses interact with specific modules in the vector’s transcriptional network

As we hypothesized, viral presence affects vector gene expression. We found that specific modules in the vector’s gene expression network respond to changes in viral titers in a species-specific way (Fig. 3b). This indicates that the host responds to viral presence. Some of these responses may represent antiviral defenses [16, 17, 19, 20] or general stress responses [41] that limit damage to the vector, though it is also possible that viruses trigger additional responses in a way that benefits their spread.

In accordance with this work, we found that the interacting modules involve both specific antiviral and non-specific stress responses. The modules included genes within the RNAi pathway (see Additional file 4), the main arthropod antiviral response [42]. These modules were enriched for infection-specific GO terms such as immune response and viral replication. At the same time, the modules were also enriched for regular cell maintenance GO terms such as cell metabolism, and gene expression regulation. It should be noted that most of the modules responding to viral infection are not obviously related either to stress or antiviral gene expression, and their function vis a vis vector or viral fitness is unclear.

Following the network analysis, we have empirically validated the module-virus interactions using gene-silencing experiments of genes with high module connectivity, these are “hub-genes” [43]. While we found that experimentally altering host gene expression does affect viral titer, the interactions were not in the direction the network analysis predicted. Namely, we found that the knockdown of CuP8, a cuticular protein that has been found to assist in plant viral transmission by binding to viruses [44], was accompanied by a significant reduction in viral load of VDV2 and ARV-2, and a non-significant reduction in DWVa (Fig. 4b). Yet, CuP8 was a hub-gene in a module negatively correlated with DWVa (Fig. 3b), suggesting that it’s expression should be negatively correlated with that of the virus. The experimental result verifies on one hand the importance of the network hub-genes in the vector-virus interaction, while on the other hand, it illustrates the high complexity of the molecular mechanism underlying the vector response, which may involve other factors such as gene-to-gene transcriptional regulation, and interaction with the environment other microorganisms [25].

A link between viral epidemiological traits and their relationship with the vector

Interestingly, viruses with similar effects on the vector’s transcriptional network tended to co-occur (Fig. 3c). This observation can be explained by two non-mutually exclusive explanations. On one hand, a high load of viral RNA and proteins can trigger defense or stress responses by the vector [41, 45]. In addition, viruses could trigger these responses as a way to manipulate the vector’s ability to spread [46, 47]. While the former mechanism is borne out by our data (see discussion of stress and antiviral pathways above), we can make several predictions about the possible manipulation of vectors by the viruses. First, we would expect a more viral species-specific response by the vector, rather than a generalized response caused by increased virus titer. Second, we predict that the effect of the virus on the vector will evolve rapidly, with closely related viruses showing different effects on the vector. We explore these predictions below.

Species-specificity of vector-virus interactions

Viruses indeed show different patterns of interaction with varroa that seem linked to their ecology. For example, some of the more pathogenic bee viruses (DWVa and DWVc) interact with the same module in the host gene expression network (modules 1 and 9, Fig. 3b). These viruses are known to be associated with varroa infestation [34, 48,49,50], and have a positive correlation with each other (Fig. 3a). A different set of modules is affected by viruses that show a high level of expression in varroa, though are detected only at low levels in honey bees (VDV2, ARV-2, and ARV-1 (Fig. 3a)) [51,52,53], suggesting that these viruses may be infecting the mites, rather than using them as a vector. Interestingly, DWVa, a major driver of bee population declines, showed strong interactions with the varroa gene expression network, but modules associated with stress and antiviral responses were largely unaffected, suggesting that this virus may avoid them (Fig. 3b). This could be because the virus does not replicate in varroa [54, 55], though it does trigger a few other gene expression responses. These data suggest that viruses interact with varroa in diverse ways and that varroa gene expression is not likely solely driven by their titer.

Our study does not shed light on the mechanisms affecting viral titer in varroa. It could be a direct function of interactions with the varroa gene regulatory network for some viruses, as suggested by our RNAi manipulation (Fig. 4). It is also possible that some viruses may passively accumulate in the mites as a result of feeding on the honey bees. However, mite-feeding behavior may be a trait that viruses could affect making the coevolutionary interactions between bees, mites, and viruses harder to disentangle.

Vector-virus interactions modified by relatively small changes in viral sequence

Virus interaction with its host can drive viral evolution and speciation [56], as host shift is a main force in viral diversification [57, 58], and can change the pathogen virulence [59]. Our findings show that this also can occur in vector-virus interactions, as we found that closely related viruses have distinctly different effects on the vector’s gene expression. The two most well-known variants of the DWV population, DWVa and DWVb are associated with varroa mites and colony collapses [60]. Varroa infection drives the DWV evolution and the relative prevalence of the two variants [37]. Despite some sequence similarity, DWVb is less associated with the mite’s infestation level in the colony, compared to DWVa [61], while on the other hand, DWVb is the only strain with concrete evidence for replicating in the mite [55]. We therefore expected that the two variants, although similar in sequence, will have a different interaction with the mite’s gene regulatory networks. Indeed, we found that the two variants show contrasting interactions with the vector modules: while DWVa shows both positive and negative interactions with the mite’s modules, DWVb shows no significant interactions with varroa modules (Fig. 3b), which is surprising since DWVb was shown to replicate in the mite [55] and should therefore interact with the vector genes. However, we should note that the analysis can be biased by viruses with high titers, such as DWVa, which extreme high viral loads in many samples may mask other virus-module and virus-virus interactions. Our findings suggest that a relatively small modification in the virus’s sequence can lead to a great change in the virus interaction with its vector, and that the vector-virus interactions are continuously and rapidly changing, resulting in a diverse viral community. Experimental infections of varroa by individual viruses can help further refine how individual viruses interact with the vector.

Co-infecting viruses do not compete in the vector

The “competitive exclusion principle” states that when competing in the same niche, one species will always suppress the other [62]. This was demonstrated also to happen in vectors, at least in mosquito cell lines co-infected with two viruses [63, 64]. In the light of the vector being a limited nourishing source, we could have expected the multi-infecting viruses to compete with each other in the varroa mite. However, we found exactly the opposite, and the only significant intra-viral interactions are positive ones (Fig. 3a). This result, along with the high viral diversity in the varroa mite, supports the model by Leeks et al. [65], which suggests that beneficial multi-viral interactions help to maintain high viral diversity. Still, our findings do not exclude the “short-sighted evolution” model of virulence, which argues that in diverse infections, faster-growing (more virulent) strains are favored because they compete for limited resources [66]. As the virulence of these viruses to varroa mite were not tested so far, we cannot conclude at this point the evolutionary mechanism which led to the observed multi-viral dynamic in the mite.

The positive correlation between some of the viruses may even imply mutualistic interactions, a phenomenon observed before for different strains of West Nile virus co-infecting mosquitoes [67]. A few mechanisms for mutualistic virus-virus relationships were suggested before such as cross-immunity, in multi-viral infection of influenza and other respiratory viral diseases [28, 29], and structural protein complementation in measles virus [68, 69], and in mutants of Dengue virus infecting mosquito cells [70]. However, the latter is the only example for vector-borne viruses. As viral colonization is the bottleneck for transmission, cooperative interaction between viruses in the vector has direct implications for the viral community dynamics, as they can favor specific viral strains that are not necessarily more virulent in a single, or even double-infected vector. This further emphasizes the need to study multi-viral infections and their molecular mechanism.

Conclusions

Our study contributes to the ongoing investigation of the way viruses interact with their vector, and how this affects the disease course, specifically multi-viral infection, a current major gap in vector-borne viral research [25, 26]. Our results imply a link between the virus epidemiological traits and its relationship with the vector. In addition, not only that co-occurring viruses interact with each other, but their abundance may predict the way these viruses will regulate their vector, and potentially, its ability to successfully infect a new host. However, experiments on these viruses’ biology and effect on the vector are needed to have the full ecological context of these findings. Hopefully, the gene network pipeline established here can be adopted to other vector-borne diseases, opening the way to study the vector and its associated pathogenic and mutualistic symbionts [71, 72].

Methods

In this study, we investigated virus-virus and vector-virus interactions by two approaches: [1] meta-transcriptomic and [2] gene-silencing experiments. In the meta-transcriptomic analysis, we looked at the overall viral landscape in the different libraries, detected vector modules (co-expressing genes using gene-network analysis), and correlated the viruses’ abundance to the vector modules. Last, we tested if the virus abundance matrix can predict the virus-vector interaction. In the second step of the study, we experimentally validated the vector-virus interaction on specific virus-gene combinations, selected based on the meta-transcriptomic analysis, by RNAi-silencing varroa hub genes and measuring the viral load. All analyses were carried out in the R statistical environment [73]. All meta-transcriptomic analyses are available and reproducible directly from the online supplementary data [74].

Vector-virus interaction, meta-transcriptomic analysis

Sequence Read Archive (SRA) data collection

To obtain varroa RNAseq data, we searched for “varroa” term in the SRA databases (NCBI, January 2020) with the following filtering criteria: “RNA” (Source filter), “RNA-seq” (Library Strategy filter) using Illumina technology, and the terms “TRANSCRIPTOMIC” and “VIRAL RNA” (Library Source filter). In total, the filtration yielded 71 libraries of varroa transcriptomes from 11 different studies. The libraries vary significantly in mite species, library preparation method and total library size, number of mites per sample, collected geographical region, sex, physiological stage, and the host species and developmental stage from which the mite was collected (for libraries details please see Additional file 8).

Reads mapping and transcripts quantification

The reads were mapped to both available varroa genome (Vdes_3.0, accession number: GCF_002443255 [75]), and to the genomes of 25 selected viruses (Fig. 1). The alignment and estimation of transcript and virus abundances in transcripts per million (TPM), was carried out using Kallisto [76] (version 0.46.1 with default options). Kallisto pseudo-aligned reads to 35, 659 identified varroa isoforms from 10,247 genes.

Mapping virus reads

The viruses were selected following a survey in the literature and NCBI genome database for viruses related to honey bee and/or varroa mite. Among the hundreds of virus sequences, we included viruses that were previously detected in varroa, in addition to honey bee viruses not detected in varroa before, as a “negative control” (Fig. 1). The final 25 viruses are mostly positive ssRNA viruses, five are negative ssRNA viruses (ARV-1 and ARV-2 (Rhabdoviridae), VOV-1 (Orthomyxoviridae), VDV4 and VDV5 (unclassified)), and another three DNA viruses (one circular dsDNA filamentous, and two ssDNA, found only in varroa mite and not in honey bees [77]). Although many DNA viruses were also found in bees [78], their abundance and importance to varroa/bee health is unknown. Therefore, these sequences were not included in the current study.

Filtering data set

Given the diversity of sources, we wanted to make sure that the input data were as homogeneous as possible. A Principal Component Analysis (PCA) of the different libraries based on varroa gene expression showed that five libraries are obvious outliers (Additional file 9). These were excluded from further analysis: library SRR8867385 from Brattel et al. [79]; libraries SRR5109825 and SRR5109827 from Remnant et al. [80], were deep sequenced for small RNA; and library SRR3927496 from a study by Levin et al. [53], in which a specific virome extraction procedure was implemented prior to library preparation. Last, library SRR533974 by Cornman et al. [81], was made of a pull of 1000 mites, which is exceptionally higher than the number of mites used in most of the studies (a pool of one - five mites). All of these exceptional sample preparation procedures may account for these libraries’ deviation from the majority of the data sets in the PCA plot. The remaining 66 libraries are distributed somewhat homogeneously in a subsequent PCA (Additional file 9), and their reads were used for further analyses.

Viral abundance analysis

Viral abundance for all 25 viruses revealed that five viruses were not detected in any of the libraries (CBPV, AFV, ANV, VPVl_46, and VPVL_36) (Fig. 2), and so were not included in further analyses. To check for virus-virus interactions within the mite, we conducted correlation matrix abundance using the Pearson correlation method for all viruses’ pairs on the log10 transformed TPM (Fig. 3a).

Weighted gene network co-expression analysis (WGCNA)

We used a network analysis approach to identify groups of genes that share a similar expression pattern across a large set of available varroa transcriptomic data (RNAseq). To construct the gene network, a weighted gene co-expression analysis was carried out using the WGCNA package in R, following the authors’ tutorial [43, 82]. The WGCNA included 4 main steps: (1) network construction and module detection; (2) correlating modules to external information, the varroa viral load; (3) identifying important genes; and [4] GO term enrichment analysis for varroa modules which interact with the viral load.

  1. (1)

    Network construction and module detection

Based on the analysis of network topology, for the construction of the network, we set our threshold for merging of modules to 0.25, minimum number of 30 genes per module, and the power β of 12. This power is the lowest for which the scale-free topology fit index curve flattens out upon reaching a high value, in this case, when Rsq reaches 0.886 (Additional file 2). We then performed hierarchical clustering of the genes based on topological overlap (sharing of network neighborhood) to identify groups of genes that co-expressed across libraries, these are the network modules (Additional file 2).

  1. (2)

    Correlating modules to viral load

To test if the varroa modules interact with the different viruses it carries, we correlated the module eigengenes to the viruses’ load (log10TPM). We used the Pearson correlation method and adjusted the p-values for multiple comparisons using the Benjamini–Hochberg method to control the false discovery rate [83] (Fig. 3b).

  1. (3)

    Identifying hub genes in network modules

Hub genes have high connectivity within a module (Module Membership) [43], and their annotation (based on sequence similarity to homologous genes). The Module Membership is calculated by Pearson correlation of the module eigengene and the gene expression. Genes with high Module Membership and relevant annotation are likely to play a role in the vector-virus interaction and are good candidates for later experimental validation.

  1. (4)

    GO term enrichment analysis for varroa modules

GO terms enrichment analyses for the genes in the significantly interacting modules, were conducted with R package GOstats using the hypergeometric test for the association of categories and genes [84]. The test parameters for each species and each ontology (biological process (BP)) using gene ID from NCBI were as follows: p-value cutoff < 0.05, not conditional, and with detection of over-represented GO terms (testDirection = over).

Correlating virus interaction and virus-varroa interactions

To test if we can predict the virus-varroa interaction given the virus abundance, we used Mantel test for correlation between two distance-matrices [40].

Vector-virus interaction, mites’ gene-silencing experiment

In the second step of the study, we experimentally validated specific gene-virus interactions, as predicted by the meta-transcriptomic model. To check if indeed the gene expression is correlated to the varroa viral load, we silenced hub genes with high Module Membership and relevant annotation (see the “Identifying hub genes in network modules” section in the meta-transcriptomic method part). We then checked for the viral load in the silenced mites. For the list of the selected genes for silencing, their accession, and annotation please see Additional file 5.

Mites and honey bee collection

Mites and bees (A. mellifera liguistica) were collected from the same colonies, at the apiary of the Okinawa Institute of Science and Technology (OIST). The hives were not treated against mites and were supplemented with sugar solution and 70% pollen cakes as necessary. Mites were collected from drone and worker pupa of different stages and were kept on bees until soaking, up to five hours from collection.

RNAi silencing of varroa genes

DsRNA preparation

For the dsRNA preparation, we first synthesized a T7-promotor attached dsDNA of each of the targeted genes by PCR amplification of cDNA prepared from a pool of 5–10 mites, with specific primers with T7-promotor attached to the 5′ end (see Additional file 10 for the primers information, and section “varroa genes primer design” for details on primer design). For PCR amplification we used Phusion™ High-Fidelity DNA Polymerase (Thermo Fisher Scientific) in a PTC-200 Peltier Thermal Cycler, MJ Research (BioRad, Toronto, Ontario) with the following steps: an initial denaturation at 98 °C for 30 s followed by 30 cycles of denaturation at 98 °C for 10 s, annealing at 60 °C for 10 s, an extension at 72 °C for 30 s, and a final extension of 72 °C for additional 5 min. We checked that the size of the amplicons matches the expected length, by running 3 μl of the PCR product in 1% agarose gel (135V, 20 min), then verified the sequence by purifying the product and Sanger sequencing on ThermoFisher SeqStudio Genetic Analyzer, using the original reverse primer and BigDye® Direct Cycle Sequencing kit (Thermo Fisher Scientific), following the manufacturer’s instructions. Prior to dsRNA synthesis, we purified the PCR products using a MinElute PCR purification kit (QIAGEN), measured their concentration by Qubit™4 Fluorometer (Life Technologies) with dsDNA HS Assay Kit (Invitrogen), and checked the product size by 4200 Tapestation (Agilent, Tokyo, Japan).

Next, 1200ng of the purified dsDNA with T7-promotor attached was used as a template for the dsRNA synthesis, using MEGAscript™ T7 Transcription Kit (Thermo Fisher Scientific). We followed the manufacturer’s protocol, with slight adjustments. The mix in a total volume of 100μl was incubated overnight at 37 °C, following the addition of TURBO DNase buffer to the reaction, and incubation for another 15 mins at 37 °C. We purified and concentrated the RNA mix using MEGAclear Transcription Clean-Up kit (Thermo Fisher Scientific), and finally measured the RNA concentration by Nanodrop, and checked the product size by 4200 Tapestation (Agilent, Tokyo, Japan). To make sure that the dsRNA effect is specific, we also prepared a negative control dsRNA of a non-target gene, green fluorescent protein (GFP), using pET6Xhn-GFPuv vector (Clontech, Takara) as a template.

Soaking experiment

For applying the dsRNA into the mite body, we followed a protocol first developed by Campbell et al. [85] and successfully done by us and others [86,87,88]. We added three mites in a 0.5ml tube containing 20μl dsRNA (2.5μg/μl) in 0.9% NaCl solution. The tubes were kept in 4 °C for 10–15 mins then we checked that all of the mites were soaked, and re-dipped them into the solution, if needed. The mites were kept in 4 °C overnight (~16 h), dried on a filter paper, then put on a bee pupa (all same age), three mites per bee in a gelatin capsule with perforated top (#1, 0.49ml, HF capsules, Matuya, Japan). Following former studies [88, 89], showing optimal silencing effectiveness 48 h post dsRNA treatment, the mites were incubated for 48 h in a controlled, dark environment at 34.5 °C, 60–75% RH, and the pupa was replaced after 24 h. After incubation, each moving-viable mite was separated in a 1.5-ml tube, snap-frozen in liquid nitrogen, and kept in −80 °C until RNA extraction. The experiment was replicated in seven experimental batches, between October and November 2020, and each batch included control group mites soaked in GFP-dsRNA of the same concentration and kept under the same conditions as described above. We checked the effect of the dsRNA treatment on mite viability, by comparing the numbers of live and dead dsRNA-treated mites to that of the control mites in each experimental batch (Fisher’s exact test for goodness of fit, p < 0.05) (Additional file 7).

RNA extraction and cDNA preparation

Each individual mite was processed following the protocol developed in our lab and described previously [90]. Briefly, each individual mite was crushed in a 1.5 ml tube dipped in liquid nitrogen, then RNA was extracted using a slightly modified TRIzol manufacturer’s protocol, with 50% volume of reagents. Total RNA quality and quantity were evaluated using Nanodrop spectrophotometer. Three hundred nanograms of purified RNA was used to synthesize a first-strand cDNA using SuperScript II (Invitrogen) and 1:2 ratio of random hexamer and oligo dT primers following the manufacturer’s protocol.

Measuring varroa gene expression and viral load

For both viruses and genes, sets of primers were designed with the NCBI primer design tool (utilizing Primer3 and BLAST), with default parameters and product size set to 100–400bp. Primer’s sequence, product length size, and gene IDs or viruses’ accession numbers are provided in Additional files 10 and 11 for varroa genes and viruses respectively. The product size of all amplicons was checked by running in 1% agarose gel (135V, 20 min).

Varroa genes primer design

For each of the varroa-selected genes (Additional file 5), we designed two sets of primers, using the gene mRNA sequence as a template. A first set for dsRNA preparation (as described in the “DsRNA preparation” section), and a second, not overlapping primers-set for gene quantification using qPCR.

Identification of local viruses and RdRp primer design

We targeted the conserved gene of RNA-dependent RNA polymerase (RdRp) commonly used for the detection and measurement of different RNA viruses in honey bees and varroa mites [91, 92]. As the genome sequence of bee viruses is slightly different across geographical regions [55, 81, 92], we first wanted to obtain the specific sequence of the strain present in our local mites. For that, for each of the three viruses (VDV2, ARV-2, and DWVa), we amplified and sequenced a wide region of the RdRp gene (amplicon size ~800bp) (see the “DsRNA preparation” section for PCR and sequencing details). To verify the sequence, we nBlasted the product to the nr database (NCBI) [93]. The reverse complementary sequence of viral RdRp was then used as a template to design the qPCR primer sets, as described above. For viruses’ amplicon sequences, nBlast results, and primer sequence and position, see Additional file 12.

qPCR

To evaluate the relative gene expression and viral load we performed a qPCR using a StepOnePlus Real-Time PCR system (Applied Biosystems Japan, Tokyo, Japan) with TB Green® Premix Ex Taq™ II (Tli RNaseH Plus, Takara). The cycling conditions were as follows: 95 °C for 30 s, followed by 40 cycles at 95 °C for 5 s, and 60 °C for 30 s. Data normalization and quantitating was done using StepOnePlus Real-Time system software (Applied Biosystems, Japan), with an automatic threshold set. Small subunit of the ribosomal RNA gene (18S rRNA) of varroa was used as a normalizing gene [94], and a control mite (treated with GFP-dsRNA) from the same experimental batch, used as the normalizing sample. For all qPCR assays a no-template control was included (data not shown). To test for differences in gene expression and viral load of dsRNA-treated mites and control mites, we used a-parametric Wilcoxon signed-ranks test. For all statistical tests, the p-values were adjusted for multiple comparisons using the Benjamini–Hochberg method to control the false discovery rate [83].

Availability of data and materials

Supplementary data material for the reported results can be found in the “Supplementary information.” All meta-transcriptomic analyses are available and reproducible directly from the online supplementary data: https://github.com/nurit-eliash/varroa-virus-networks [74].

Abbreviations

PCA:

Principal Component Analysis

OIST:

Okinawa Institute of Science and Technology

RdRp:

RNA-dependent RNA polymerase

References

  1. Brault AC, Savage HM, Duggal NK, Eisen RJ, Staples JE. Heartland virus epidemiology, vector association, and disease potential. Viruses. 2018;10(9):498.

    Article  Google Scholar 

  2. de la Fuente J, Antunes S, Bonnet S, Cabezas-Cruz A, Domingos AG, Estrada-Peña A, et al. Tick-Pathogen Interactions and Vector Competence: Identification of Molecular Drivers for Tick-Borne Diseases. Front Cell Infect Microbiol. 2017;7:114.

    Google Scholar 

  3. Ng JCK, Falk BW. Virus-vector interactions mediating nonpersistent and semipersistent transmission of plant viruses. Annu Rev Phytopathol. 2006;44:183–212.

    Article  CAS  Google Scholar 

  4. Whitfield AE, Falk BW, Rotenberg D. Insect vector-mediated transmission of plant viruses. Virology. 2015;479–480:278–89.

    Article  Google Scholar 

  5. Olival KJ, Hosseini PR, Zambrana-Torrelio C, Ross N, Bogich TL, Daszak P. Host and viral traits predict zoonotic spillover from mammals. Nature. 2017;546(7660):646–50.

    Article  CAS  Google Scholar 

  6. Weaver SC, Charlier C, Vasilakis N, Lecuit M. Zika, Chikungunya, and Other Emerging Vector-Borne Viral Diseases. Annu Rev Med. 2018;69:395–408.

    Article  CAS  Google Scholar 

  7. World Health Organization. Global vector control response: an integrated approach for the control of vector-borne diseases. 2017. Available from: https://apps.who.int/gb/ebwha/pdf_files/WHA70/A70_R16-en.pdf?ua=1

    Google Scholar 

  8. Rocklöv J, Dubrow R. Climate change: an enduring challenge for vector-borne disease prevention and control. Nat Immunol. 2020;21(5):479–83.

    Article  Google Scholar 

  9. Sutherst RW. Global change and human vulnerability to vector-borne diseases. Clin Microbiol Rev. 2004;17(1):136–73.

    Article  Google Scholar 

  10. Wilson AL, Courtenay O, Kelly-Hope LA, Scott TW, Takken W, Torr SJ, et al. The importance of vector control for the control and elimination of vector-borne diseases. PLoS Negl Trop Dis. 2020;14(1):e0007831.

    Article  CAS  Google Scholar 

  11. Jackson BT, Brewster CC, Paulson SL. La Crosse virus infection alters blood feeding behavior in Aedes triseriatus and Aedes albopictus (Diptera: Culicidae). J Med Entomol. 2012;49(6):1424–9.

    Article  Google Scholar 

  12. Maciel-de-Freitas R, Koella JC, Lourenço-de-Oliveira R. Lower survival rate, longevity and fecundity of Aedes aegypti (Diptera: Culicidae) females orally challenged with dengue virus serotype 2. Trans R Soc Trop Med Hyg. 2011;105(8):452–8.

    Article  CAS  Google Scholar 

  13. Moncayo AC, Edman JD, Turell MJ. Effect of eastern equine encephalomyelitis virus on the survival of Aedes albopictus, Anopheles quadrimaculatus, and Coquillettidia perturbans (Diptera: Culicidae). J Med Entomol. 2000;37(5):701–6.

    Article  CAS  Google Scholar 

  14. Neelakanta G, Sultana H, Fish D, Anderson JF, Fikrig E. Anaplasma phagocytophilum induces Ixodes scapularis ticks to express an antifreeze glycoprotein gene that enhances their survival in the cold. J Clin Invest. 2010;120(9):3179–90.

    Article  CAS  Google Scholar 

  15. Cime-Castillo J, Delannoy P, Mendoza-Hernández G, Monroy-Martínez V, Harduin-Lepers A, Lanz-Mendoza H, et al. Sialic acid expression in the mosquito Aedes aegypti and its possible role in dengue virus-vector interactions. Biomed Res Int. 2015;2015:504187.

    Article  Google Scholar 

  16. Göertz GP, van Bree JWM, Hiralal A, Fernhout BM, Steffens C, Boeren S, et al. Subgenomic flavivirus RNA binds the mosquito DEAD/H-box helicase ME31B and determines Zika virus transmission by Aedes aegypti. Proc Natl Acad Sci U S A. 2019;116(38):19136–44.

    Article  Google Scholar 

  17. Schnettler E, Tykalová H, Watson M, Sharma M, Sterken MG, Obbard DJ, et al. Induction and suppression of tick cell antiviral RNAi responses by tick-borne flaviviruses. Nucleic Acids Res. 2014;42(14):9436–46.

    Article  CAS  Google Scholar 

  18. Luplertlop N, Surasombatpattana P, Patramool S, Dumas E, Wasinpiyamongkol L, Saune L, et al. Induction of a peptide with activity against a broad spectrum of pathogens in the Aedes aegypti salivary gland, following Infection with Dengue Virus. PLoS Pathog. 2011;7(1):e1001252.

    Article  CAS  Google Scholar 

  19. Zink SD, Van Slyke GA, Palumbo MJ, Kramer LD, Ciota AT. Exposure to West Nile Virus Increases Bacterial Diversity and Immune Gene Expression in Culex pipiens. Viruses. 2015;7(10):5619–31.

    Article  CAS  Google Scholar 

  20. Huang Y-JS, Higgs S, Vanlandingham DL. Arbovirus-Mosquito Vector-Host Interactions and the Impact on Transmission and Disease Pathogenesis of Arboviruses. Front Microbiol. 2019;10:22.

    Article  Google Scholar 

  21. Goenaga S, Kenney JL, Duggal NK, Delorey M, Ebel GD, Zhang B, et al. Potential for Co-Infection of a Mosquito-Specific Flavivirus, Nhumirim Virus, to Block West Nile Virus Transmission in Mosquitoes. Viruses. 2015;7(11):5801–12.

    Article  CAS  Google Scholar 

  22. Göertz GP, Vogels CBF, Geertsema C, Koenraadt CJM, Pijlman GP. Mosquito co-infection with Zika and chikungunya virus allows simultaneous transmission without affecting vector competence of Aedes aegypti. PLoS Negl Trop Dis. 2017;11(6):e0005654.

    Article  Google Scholar 

  23. Batovska J, Mee PT, Lynch SE, Sawbridge TI, Rodoni BC. Sensitivity and specificity of metatranscriptomics as an arbovirus surveillance tool. Sci Rep. 2019;9(1):19398.

    Article  CAS  Google Scholar 

  24. Batson J, Dudas G, Haas-Stapleton E, Kistler AL, Li LM, Logan P, et al. Single mosquito metatranscriptomics identifies vectors, emerging pathogens and reservoirs in one assay. Elife. 2021;10. https://doi.org/10.7554/eLife.68353.

  25. Ciota AT. The role of co-infection and swarm dynamics in arbovirus transmission. Virus Res. 2019;265:88–93.

    Article  CAS  Google Scholar 

  26. Vogels CBF, Rückert C, Cavany SM, Perkins TA, Ebel GD, Grubaugh ND. Arbovirus coinfection and co-transmission: A neglected public health concern? PLoS Biol. 2019;17(1):e3000130.

    Article  Google Scholar 

  27. Erez Z, Steinberger-Levy I, Shamir M, Doron S, Stokar-Avihail A, Peleg Y, et al. Communication between viruses guides lysis-lysogeny decisions. Nature. 2017;541(7638):488–93.

    Article  CAS  Google Scholar 

  28. Ferguson NM, Galvani AP, Bush RM. Ecological and immunological determinants of influenza evolution. Nature. 2003;422(6930):428–33.

    Article  CAS  Google Scholar 

  29. Nickbakhsh S, Mair C, Matthews L, Reeve R, Johnson PCD, Thorburn F, et al. Virus-virus interactions impact the population dynamics of influenza and the common cold. Proc Natl Acad Sci U S A. 2019. https://doi.org/10.1073/pnas.1911083116.

  30. Alcaide C, Rabadán MP, Moreno-Pérez MG, Gómez P. Implications of mixed viral infections on plant disease ecology and evolution. Adv Virus Res. 2020;106:145–69.

    Article  CAS  Google Scholar 

  31. Díaz-Muñoz SL. Viral coinfection is shaped by host ecology and virus-virus interactions across diverse microbial taxa and environments. Virus Evol. 2017;3(1):vex011.

    Article  Google Scholar 

  32. Sanjuán R, Illingworth CJR, Geoghegan JL, Iranzo J, Zwart MP, Ciota AT, et al. Five challenges in the field of viral diversity and evolution. 2021; Available from: https://aspace.repository.cam.ac.uk/bitstream/handle/1810/322603/fviro-01-684949.pdf?sequence=2&isAllowed=y

    Book  Google Scholar 

  33. Wilson AJ, Morgan ER, Booth M, Norman R, Perkins SE, Hauffe HC, et al. What is a vector? Philos Trans R Soc Lond Ser B Biol Sci. 2017;372(1719). https://doi.org/10.1098/rstb.2016.0085.

  34. McMenamin AJ, Genersch E. Honey bee colony losses and associated viruses. Curr Opini Insect Sci. 2015;8:121–9.

    Article  Google Scholar 

  35. Steinhauer N, Kulhanek K, Antúnez K, Human H, Chantawannakul P, Chauzat MP, et al. Drivers of colony losses. Curr Opin Insect Sci. 2018;26:142–8.

    Article  Google Scholar 

  36. Carreck NL, Ball BV, Martin SJ. Honey bee colony collapse and changes in viral prevalence associated with Varroa destructor. J Apic Res. 2010;49(1):93–4.

    Article  Google Scholar 

  37. Martin SJ, Highfield AC, Brettell L, Villalobos EM, Budge GE, Powell M, et al. Global honey bee viral landscape altered by a parasitic mite. Science. 2012;336(6086):1304–6.

    Article  CAS  Google Scholar 

  38. Yañez O, Chávez-Galarza J, Tellgren-Roth C, Pinto MA, Neumann P, de Miranda JR. The honeybee (Apis mellifera) developmental state shapes the genetic composition of the deformed wing virus-A quasispecies during serial transmission. Sci Rep. 2020;10(1):5956.

    Article  Google Scholar 

  39. Nganso BT, Sela N, Soroker V. A genome-wide screening for RNAi pathway proteins in Acari. BMC Genomics. 2020;21:791.

    Article  CAS  Google Scholar 

  40. Mantel N. The detection of disease clustering and a generalized regression approach. Cancer Res. 1967;27(2):209–20.

    CAS  Google Scholar 

  41. Rosche KL, Sidak-Loftis LC, Hurtado J, Fisk EA, Shaw DK. Arthropods Under Pressure: Stress Responses and Immunity at the Pathogen-Vector Interface. Front Immunol. 2021;11:629777.

    Article  Google Scholar 

  42. Blair CD, Olson KE. The role of RNA interference (RNAi) in arbovirus-vector interactions. Viruses. 2015;7(2):820–43.

    Article  CAS  Google Scholar 

  43. Langfelder P, Horvath S. WGCNA: An R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559 Available from: http://www.biomedcentral.com/1471-2105/9/559.

    Article  Google Scholar 

  44. Deshoux M, Monsion B, Uzest M. Insect cuticular proteins and their role in transmission of phytoviruses. Curr Opin Virol. 2018;33:137–43.

    Article  CAS  Google Scholar 

  45. Baxter RHG, Contet A, Krueger K. Arthropod Innate Immune Systems and Vector-Borne Diseases. Biochemistry. 2017;56(7):907–18.

    Article  CAS  Google Scholar 

  46. Hurd H. Manipulation of medically important insect vectors by their parasites. Annu Rev Entomol. 2003;48:141–61.

    Article  CAS  Google Scholar 

  47. Targett GAT. Parasites, arthropod vectors, and immune responses. Parasite Immunol. 2006;28(4):117–9.

    Article  CAS  Google Scholar 

  48. Benjeddou M, Leat N, Allsopp M, Davison S. Detection of acute bee paralysis virus and black queen cell virus from honeybees by reverse transcriptase pcr. Appl Environ Microbiol. 2001;67(5):2384–7.

    Article  CAS  Google Scholar 

  49. Daughenbaugh KF, Martin M, Brutscher LM, Cavigli I, Garcia E, Lavin M, et al. Honey bee infecting Lake Sinai viruses. Viruses. 2015;7(6):3285–309.

    Article  CAS  Google Scholar 

  50. Kevill JL, Highfield A, Mordecai GJ, Martin SJ, Schroeder DC. ABC Assay: Method Development and Application to Quantify the Role of Three DWV Master Variants in Overwinter Colony Losses of European Honey Bees. Viruses. 2017;9(11). https://doi.org/10.3390/v9110314.

  51. Chen G, Wang S, Jia S, Feng Y, Hu F, Chen Y, et al. A New Strain of Virus Discovered in China Specific to the Parasitic Mite Varroa destructor Poses a Potential Threat to Honey Bees. Viruses. 2021;13(4):679.

    Article  CAS  Google Scholar 

  52. Levin S, Galbraith D, Sela N, Erez T, Grozinger CM, Chejanovsky N. Presence of Apis rhabdovirus-1 in populations of pollinators and their parasites from two continents. Front Microbiol. 2017;8:2482.

    Article  Google Scholar 

  53. Levin S, Sela N, Chejanovsky N. Two novel viruses associated with the Apis mellifera pathogenic mite Varroa destructor. Sci Rep. 2016;6:37710.

    Article  CAS  Google Scholar 

  54. Posada-Florez F, Childers AK, Heerman MC, Egekwu NI, Cook SC, Chen Y, et al. Deformed wing virus type A, a major honey bee pathogen, is vectored by the mite Varroa destructor in a non-propagative manner. Sci Rep. 2019;9(1):12445.

    Article  Google Scholar 

  55. Gisder S, Genersch E. Direct Evidence for Infection of Varroa destructor Mites with the Bee-Pathogenic Deformed Wing Virus Variant B - but Not Variant A - via Fluorescence-in situ-Hybridization Analysis. J Virol. 2020. https://doi.org/10.1128/JVI.01786-20.

  56. Koskella B, Brockhurst MA. Bacteria-phage coevolution as a driver of ecological and evolutionary processes in microbial communities. FEMS Microbiol Rev. 2014;38(5):916–31.

    Article  CAS  Google Scholar 

  57. Geoghegan JL, Duchêne S, Holmes EC. Comparative analysis estimates the relative frequencies of co-divergence and cross-species transmission within viral families. PLoS Pathog. 2017;13(2):e1006215.

    Article  Google Scholar 

  58. Ricklefs RE, Outlaw DC, Svensson-Coelho M, Medeiros MCI, Ellis VA, Latta S. Species formation by host shifting in avian malaria parasites. Proc Natl Acad Sci U S A. 2014;111(41):14816–21.

    Article  CAS  Google Scholar 

  59. Williams PD, Kamel SJ. The evolution of pathogen virulence: Effects of transitions between host types. J Theor Biol. 2018;438:1–8.

    Article  Google Scholar 

  60. Martin SJ, Brettell LE. Deformed wing virus in honeybees and other insects. Annu Rev Virol. 2019;6(1). https://doi.org/10.1146/annurev-virology-092818-.

  61. Norton AM, Remnant EJ, Tom J, Buchmann G, Blacquiere T, Beekman M. Adaptation to vector-based transmission in a honey bee virus. J Anim Ecol. 2021. https://doi.org/10.1111/1365-2656.13493.

  62. Domingo E. Chapter 6 - Virus Population Dynamics Examined with Experimental Model Systems. In: Domingo E, editor. Virus as Populations. Boston: Academic Press; 2016. p. 197–225.

    Chapter  Google Scholar 

  63. Abrao EP, da Fonseca BAL. Infection of Mosquito Cells (C6/36) by Dengue-2 Virus Interferes with Subsequent Infection by Yellow Fever Virus. Vector Borne Zoonotic Dis. 2016;16(2):124–30.

    Article  Google Scholar 

  64. Pepin KM, Lambeth K, Hanley KA. Asymmetric competitive suppression between strains of dengue virus. BMC Microbiol. 2008;8:28.

    Article  Google Scholar 

  65. Leeks A, Segredo-Otero EA, Sanjuán R, West SA. Beneficial coinfection can promote within-host viral diversity. Virus Evol. 2018;4(2):vey028.

    Article  Google Scholar 

  66. Levin BR, Bull JJ. Short-sighted evolution and the virulence of pathogenic microorganisms. Trends Microbiol. 1994;2(3):76–81.

    Article  CAS  Google Scholar 

  67. Ciota AT, Ehrbar DJ, Van Slyke GA, Willsey GG, Kramer LD. Cooperative interactions in the West Nile virus mutant swarm. BMC Evol Biol. 2012;12:58.

    Article  CAS  Google Scholar 

  68. Shirogane Y, Watanabe S, Yanagi Y. Cooperation between different RNA virus genomes produces a new phenotype. Nat Commun. 2012;3:1235.

    Article  Google Scholar 

  69. Shirogane Y, Watanabe S, Yanagi Y. Cooperative Interaction Within RNA Virus Mutant Spectra. In: Domingo E, Schuster P, editors. Quasispecies: From Theory to Experimental Systems. Cham: Springer International Publishing; 2016. p. 219–29.

    Google Scholar 

  70. Aaskov J, Buzacott K, Thu HM, Lowry K, Holmes EC. Long-term transmission of defective RNA viruses in humans and Aedes mosquitoes. Science. 2006;311(5758):236–8.

    Article  CAS  Google Scholar 

  71. Hegde S, Rasgon JL, Hughes GL. The microbiome modulates arbovirus transmission in mosquitoes. Curr Opin Virol. 2015;15:97–102.

    Article  CAS  Google Scholar 

  72. Rainey SM, Shah P, Kohl A, Dietrich I. Understanding the Wolbachia-mediated inhibition of arboviruses in mosquitoes: progress and challenges. J Gen Virol. 2014;95(Pt 3):517–30.

    Article  CAS  Google Scholar 

  73. Team RC. R: A language and environment for statistical computing. 2013; Available from: https://www.R-project.org/

    Google Scholar 

  74. Eliash N. varroa-virus-networks. 2021. Available from: https://github.com/nurit-eliash/varroa-virus-networks

    Google Scholar 

  75. Techer MA, Rane RV, Grau ML, Roberts JMK, Sullivan ST, Liachko I, et al. Divergent evolutionary trajectories following speciation in two ectoparasitic honey bee mites. Commun Biol. 2019;2(1):357.

    Article  Google Scholar 

  76. Bray NL, Pimentel H, Melsted P, Pachter L. Near-optimal probabilistic RNA-seq quantification. Nat Biotechnol. 2016;34(5):525–7.

    Article  CAS  Google Scholar 

  77. Kraberger S, Visnovsky GA, van Toor RF, Male MF, Waits K, Fontenele RS, et al. Genome Sequences of Two Single-Stranded DNA Viruses Identified in Varroa destructor. Genome Announc. 2018;6(9). https://doi.org/10.1128/genomeA.00107-18.

  78. Kraberger S, Cook CN, Schmidlin K, Fontenele RS, Bautista J, Smith B, et al. Diverse single-stranded DNA viruses associated with honey bees (Apis mellifera). Infect Genet Evol. 2019;71:179–88.

    Article  CAS  Google Scholar 

  79. Brettell LE, Schroeder DC, Martin SJ. RNAseq analysis reveals virus diversity within hawaiian apiary insect communities. Viruses. 2019;11(5). https://doi.org/10.3390/v11050397.

  80. Remnant EJ, Shi M, Buchmann G, Blacquière T, Holmes EC, Beekman M, et al. A diverse range of novel RNA viruses in geographically distinct honey bee populations. J Virol. 2017;91(16):1–19.

    Article  Google Scholar 

  81. Cornman RS, Boncristiani H, Dainat B, Chen Y, VanEngelsdorp D, Weaver D, et al. Population-genomic variation within RNA viruses of the Western honey bee, Apis mellifera, inferred from deep sequencing. BMC Genomics. 2013;14(1):154.

    Article  CAS  Google Scholar 

  82. Langfelder P, Horvath S. Tutorials for WGCNA R package. 2016. Available from: https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/ [cited 20 Jul 2021]

    Google Scholar 

  83. Benjamini Y, Hochberg Y. Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. J R Stat Soc Ser B Stat Methodol. 1995;57(1):289–300.

    Google Scholar 

  84. Falcon S, Gentleman R. Using GOstats to test gene lists for GO term association. Bioinformatics. 2007;23(2):257–8.

    Article  CAS  Google Scholar 

  85. Campbell EM, Budge GE, Bowman AS. Gene-knockdown in the honey bee mite Varroa destructor by a non-invasive approach: studies on a glutathione S-transferase. Parasit Vectors. 2010;3:73.

    Article  Google Scholar 

  86. Garbian Y, Maori E, Kalev H, Shafir S, Sela I. Bidirectional transfer of RNAi between honey bee and Varroa destructor: Varroa gene silencing reduces Varroa population. PLoS Pathog. 2012;8(12):e1003035.

    Article  CAS  Google Scholar 

  87. Nganso BT, Mani K, Eliash N, Rafaeli A, Soroker V. Towards disrupting Varroa -honey bee chemosensing: A focus on a Niemann-Pick type C2 transcript. Insect Mol Biol. 2021. https://doi.org/10.1111/imb.12722.

  88. Singh NK, Eliash N, Stein I, Kamer Y, Ilia Z, Rafaeli A, et al. Identification and gene-silencing of a putative odorant receptor transcription factor in Varroa destructor: possible role in olfaction. Insect Mol Biol. 2016;25(2):181–90.

    Article  CAS  Google Scholar 

  89. Campbell EM, Budge GE, Watkins M, Bowman AS. Transcriptome analysis of the synganglion from the honey bee mite, Varroa destructor and RNAi knockdown of neural peptide targets. Insect Biochem Mol Biol. 2016;70:116–26.

    Article  CAS  Google Scholar 

  90. Hasegawa N, Techer M, Mikheyev AS. A toolkit for studying Varroa genomics and transcriptomics: preservation, extraction, and sequencing library preparation. BMC Genomics. 2021;22(1):54.

    Article  CAS  Google Scholar 

  91. de Miranda JR, Bailey L, Ball BV, Blanchard P, Budge GE, Chejanovsky N, et al. Standard methods for virus research in Apis mellifera. J Apic Res. 2013;52(4):1–56.

    Article  Google Scholar 

  92. Gisder S, Möckel N, Eisenhardt D, Genersch E. In vivo evolution of viral virulence: switching of deformed wing virus between hosts results in virulence changes and sequence shifts. Environ Microbiol. 2018;20(12):4612–28.

    Article  CAS  Google Scholar 

  93. Altschul SF, Gish W, Miller WT, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215(3):403–10.

    Article  CAS  Google Scholar 

  94. Campbell EM, McIntosh CH, Bowman AS. A Toolbox for Quantitative Gene Expression in Varroa destructor: RNA Degradation in Field Samples and Systematic Analysis of Reference Gene Stability. PLoS One. 2016;11(5):e0155640.

    Article  Google Scholar 

Download references

Acknowledgements

We wish to thank Dr. Olesya Gusachenko (University of St Andrews, UK) for consulting regarding the virus identification and viruses’ primer design. We also wish to thank Dr. Maeva Techer for her help in the tedious task of mites’ collection for the silencing experiment, and for all Mikheyev unit members at OIST for stimulating discussion and support throughout.

Funding

Not applicable

Author information

Authors and Affiliations

  1. Shamir Research Institute, University of Haifa, Katzrin, Israel

    Nurit Eliash

  2. Okinawa Institute of Science and Technology, 1919-1 Tancha Onna-son, Okinawa, 904-0495, Japan

    Nurit Eliash, Miyuki Suenaga & Alexander S. Mikheyev

  3. Australian National University, Canberra, ACT, 2600, Australia

    Alexander S. Mikheyev

Authors
  1. Nurit Eliash
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Miyuki Suenaga
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Alexander S. Mikheyev
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

NE and ASM conceptualized the study. ASM mapped the varroa and virus reads and was a major contributor in writing the manuscript. NE analyzed and interpreted the WGCNA and the bioinformatic analysis, supervised by ASM. NE collected bees and mites and designed the silencing experiments and analyzed and interpreted the gene expression and viral load results. MS performed the RNA extraction, qPCR and dsRNA preparation. NE and MS preformed the soaking experiment. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Nurit Eliash or Alexander S. Mikheyev.

Ethics declarations

Ethics approval and consent to participate

Not applicable

Consent for publication

Not applicable

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1.

Viruses mapped to the varroa RNAseq libraries in the meta transcriptomic analysis.

Additional file 2.

Network construction using 10,247 genes of 66 SRA varroa libraries. a. Picking soft threshold. b. Hierarchical clustering dendrogram using merge cut height of 25%, revealing 12 co-expressed genes modules. Each branch of the dendrogram represents a single gene, and the colored bar below denotes its corresponding module, as annotated in the legend to the right. The dendrogram height is the distance between genes.

Additional file 3.

The significant GO-terms based on GO-term enrichment analysis (adjusted p-value<0.05), of eight modules. These modules were significantly interacting with at least one of the viruses (based on the analysis in Fig. 3b).

Additional file 4.

Varroa RNAi pathway homolog genes, and the module it belongs to in the current network analysis.

Additional file 5.

Varroa hub-genes targeted for RNAi silencing.

Additional file 6.

The probability that the difference between RNAi-treated and control mites (treated in GFP-dsRNA solution) have occurred by chance, for relative gene expression and viral load (Wilcoxon signed-ranks test, followed by FDR-correction).

Additional file 7.

Fisher exact test of goodness of fit, testing the difference in live and dead mite distribution between RNAi-treated and control mites (treated in GFP-dsRNA solution), in each experimental batch, ’Date’.

Additional file 8.

Varroa SRA libraries used in the meta-transcriptomic analysis.

Additional file 9.

PCA of varroa SRA libraires based on their genes TPM. a. All initial 71 libraries. The outlier libraries are circled in red. These 5 libraries were excluded from further analysis. b. the final 66 libraries used for the analysis.

Additional file 10.

Primers' sequences for varroa genes' primers for gene quantification using qPCR, and for dsRNA preparation.

Additional file 11.

Primers' sequences for viruses RdRp sequencing and for viral load quantification using qPCR.

Additional file 12.

Viruses’ RdRp amplicon sequences, nBlast results and primer positions, for (a) DWVa, (b) VDV2 and (c) ARV-2.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Eliash, N., Suenaga, M. & Mikheyev, A.S. Vector-virus interaction affects viral loads and co-occurrence. BMC Biol 20, 284 (2022). https://doi.org/10.1186/s12915-022-01463-4

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12915-022-01463-4

Keywords

  • Vector-virus interaction
  • Gene network analysis
  • Varroa
  • RNAi silencing

BMC Biology

ISSN: 1741-7007

Contact us