A sting in the spit: widespread cross-infection of multiple RNA viruses across wild and managed bees.

Declining populations of bee pollinators are a cause of concern, with major repercussions for biodiversity loss and food security. RNA viruses associated with honeybees represent a potential threat to other insect pollinators, but the extent of this threat is poorly understood. This study aims to attain a detailed understanding of the current and ongoing risk of emerging infectious disease (EID) transmission between managed and wild pollinator species across a wide range of RNA viruses. Within a structured large-scale national survey across 26 independent sites, we quantify the prevalence and pathogen loads of multiple RNA viruses in co-occurring managed honeybee (Apis mellifera) and wild bumblebee (Bombus spp.) populations. We then construct models that compare virus prevalence between wild and managed pollinators. Multiple RNA viruses associated with honeybees are widespread in sympatric wild bumblebee populations. Virus prevalence in honeybees is a significant predictor of virus prevalence in bumblebees, but we remain cautious in speculating over the principle direction of pathogen transmission. We demonstrate species-specific differences in prevalence, indicating significant variation in disease susceptibility or tolerance. Pathogen loads within individual bumblebees may be high and in the case of at least one RNA virus, prevalence is higher in wild bumblebees than in managed honeybee populations. Our findings indicate widespread transmission of RNA viruses between managed and wild bee pollinators, pointing to an interconnected network of potential disease pressures within and among pollinator species. In the context of the biodiversity crisis, our study emphasizes the importance of targeting a wide range of pathogens and defining host associations when considering potential drivers of population decline.

Introduction

The ongoing biodiversity crisis threatens human health and global food security (Cardinale et al. 2012). Emerging infectious diseases (EIDs) have contributed significantly to species declines (Daszak, Cunningham & Hyatt 2000), with lethal chytridiomycosis in amphibians (Fisher, Garner & Walker 2009) and white‐nose syndrome in bats (Blehert et al. 2009) representing prominent recent examples. Infectious diseases may emerge through association with a host species (a ‘reservoir’) in which pathogens have become established, or where disease epidemiology may have recently changed due to perturbation (e.g. through arrival of a novel disease, or disease vector). The switching of pathogens between host species is a major cause of epidemics in humans and other vertebrate hosts (Woolhouse, Haydon & Antia 2005), and EIDs have potentially profound impacts on invertebrates providing important ecosystem services, which secure food production. However, the extent to which EIDs are an issue in invertebrates – and in insect pollinators particularly – is not clear.

Bees provide an essential ecosystem service in the form of crop pollination (Klein et al. 2007), but they are under pressure globally (Biesmeijer et al. 2006; Goulson, Lye & Darvill 2008; Brown & Paxton 2009; Williams & Osborne 2009; Potts et al. 2010; Vanbergen et al. 2013). Bumblebees are major wild pollinators in northern temperate climates (Goulson 2009), but they are declining in both the Old World (Williams 1982; Fitzpatrick et al. 2007; Kosior et al. 2007) and the New World (Bartomeus et al. 2013), with EIDs implicated as a cause of these declines (Cameron et al. 2011; Meeus et al. 2011; Fürst et al. 2014; Schmid‐Hempel et al. 2014). EIDs are known to be a major threat to the most widely used commercial pollinator, the honeybee (Apis mellifera), with the exotic ectoparasitic mite, Varroa destructor, meriting particular attention. The mite has risen to prominence due to its ability to act as a vector of several RNA viruses that previously persisted relatively benignly in honeybee colonies, most notably deformed wing virus (DWV), but also viruses belonging to the acute bee paralysis virus (ABPV) complex (Genersch & Aubert 2010) and slow bee paralysis virus (SBPV) (Carreck, Ball & Martin 2010; Santillán‐Galicia et al. 2014). In the case of DWV, the arrival of V. destructor has been directly linked to increased prevalence and virus loads in honeybees (Martin et al. 2012).

Wild pollinators harbour pathogens previously associated with honeybees (Genersch et al. 2006; Singh et al. 2010; Peng et al. 2011; Evison et al. 2012; Graystock et al. 2013; Levitt et al. 2013; Ravoet et al. 2014), and for at least one emerging RNA virus, disease in managed honeybees and wild bumblebees is linked (Fürst et al. 2014). The association of pathogens with managed honeybees is in part a reflection of study bias, but the trend may also point to an emerging problem of infectious RNA viruses in wild bees – triggered, perhaps, by the arrival of V. destructor mites in the western honeybee some 40 years ago (Rosenkranz, Aumeier & Ziegelmann 2010). For the great majority of RNA viruses in wild bees, detailed knowledge of prevalence and level of infection (pathogen load) is still lacking. This represents a significant gap in understanding, particularly given the prominent role that RNA virus diseases are believed to play in causing managed honeybee colony loss (Schroeder & Martin 2012).

We therefore conducted a comprehensive field analysis of honeybee and wild bumblebee populations across Great Britain and the Isle of Man to (i) understand the contemporary landscape prevalence of common RNA viruses thought to be associated with honeybees, (ii) quantify and compare the individual infection levels of RNA viruses in bee foragers, (iii) assess the extent to which RNA virus spillover is occurring between honeybees and bumblebees (in either direction). We show that multiple RNA viruses are prevalent in wild bee populations and present evidence for recent and widespread circulation of viral diseases between Britain's primary managed and wild bee pollinators.

Materials and methods

Field Sampling and RNA Extraction

Field sampling methodology and RNA extraction follows Fürst et al. (2014). Briefly, we collected free flying honeybees and bumblebees from flowers at 26 sites (A‐Z) across Great Britain and the Isle of Man, each separated by at least 30 km (mean ± SD distance in km = 69·21 ± 26·39). The collection area covered at least 1000 m2 at each location, and where possible, all bees were collected within a single day. Time taken (in minutes) to collect 20 A. mellifera and 20 Bombus spp. individuals was recorded as an estimate of abundance. Honeybees and up to four species of bumblebees from each site were then screened for the presence and quantity of a range of viruses. Honeybee or bumblebee abdomens were bisected longitudinally, one‐half of which was submerged in RLT buffer and disrupted in a Tissue lyser II (Qiagen, Manchester, UK) at 30 Hz for 2 min followed by 20 Hz for 2 min prior to RNA isolation. Total RNA was extracted manually using the RNeasy mini kit (Qiagen, Manchester, UK) following manufacturer's instructions.

Pathogen Detection

We screened for a wide range of known positive‐sense single‐stranded RNA viruses, by employing multiple ligation‐dependent probe amplification (MLPA) using the RT‐MLPA® kit (MRC‐Holland, Amsterdam, Netherlands). We used probes designed for the positive strand of the following six composite positive‐sense single‐stranded RNA virus targets (De Smet et al. 2012): (i) black queen cell virus (BQCV); (ii) deformed wing virus, Varroa destructor virus and kakugo virus (DWV/KV/VDV‐1); (iii) acute bee paralysis virus, Israeli acute paralysis virus and Kashmir bee virus (ABPV/IAPV/KBV); (iv) slow bee paralysis virus (SBPV); (v) sac brood virus (SBV); and (vi) chronic bee paralysis virus (CBPV) and, as a housekeeping (control) gene, β‐actin. Notable viruses such as the Lake Sinai viruses (LSV 1 and 2) have been recorded in North America (Runckel et al. 2011) and Europe (Granberg et al. 2013). While these were not included in the current study, we acknowledge that they may also be transmitted across species. Amplified fragments were resolved by capillary electrophoresis on a QIAxcel (Qiagen, Hilden, Germany), using a positive virus acceptance threshold of 0·1 relative fluorescence units. Samples were excluded from further analysis if the housekeeping gene, β‐actin, fell below this threshold.

For each MLPA positive virus target, samples were analysed by qRT‐PCR to identify the specific virus (in the case of DWV/KV/VDV‐1 and ABPV/IAPV/KBV) and to estimate individual viral load. This allowed us to differentiate between VDV‐1 and DWV/KV and between ABPV, IAPV and KBV. KV and DWV are very closely related and were not differentiated by qRT‐PCR. Total cDNA was synthesized using M‐MLV Revertase (Promega, Mannheim, Germany) following manufacturer's instructions, using 500 ng of sample RNA. For absolute quantification, duplicate qRT‐PCR was performed for each sample with a Bio‐Rad C1000, using SYBRgreen Sensimix (Bioline, Luckenwalde, Germany) in the following program: 5 min at 95 °C, followed by 40 cycles of 10 s at 95 °C, 30 s at 57 °C and 30 s at 72 °C (read). Duplicate β‐actin reactions were also amplified for all samples as an internal reference marker. A negative control containing RNA‐free HPLC water and a virus‐positive sample were included as controls in each reaction run. To account for potential variation in sample quality, an upper cycle threshold (Ct) of 35 was set for β‐actin, above which samples were not included in quantitative analysis. Given the previous positive detection of virus by MLPA, an upper threshold for virus quantification by qRT‐PCR was not applied. We used specific primers for the following viruses: BQCV; DWV; VDV‐1; ABPV; IAPV; KBV; SBPV; SBV (see Table S1, Supporting information). Following PCR, DNA was denatured for 1 min at 95 °C and cooled to 55 °C for 1 min. A melting profile was generated from 55 to 95 °C (0·5 °C per second increments). Quantification was calculated using duplicate DNA standard curves of purified flanking PCR products (DWV, VDV‐1, Table S1, Supporting information for primers) or plasmids (BQCV; ABPV; SBPV), with efficiencies of 98·4% (DWV), 99·9% (VDV‐1), 96·2% (BQCV), 101·3% (ABPV) and 93·1% (SBPV), and correlation coefficients (R 2) from 0·995 to 0·999.

Sequencing

To confirm the identity of viruses, we cloned and sequenced virus fragments from single honeybees and up to two bumblebees that contained high levels of BQCV, ABPV or SBPV. qRT‐PCR products were purified using the Qiaquick PCR Purification Kit (Qiagen, Hilden, Germany) and cloned directly using the pGEM T Easy Vector system (Promega, Mannheim, Germany) following manufacturer's instructions. Plasmid DNA was isolated using a Spin Miniprep Kit (Qiagen, Hilden, Germany). Up to five clones per sample were sequenced in forward and reverse orientation (GATC Biotech, Constance, Germany), and aligned by eye to genome references of BQCV (NC_003784), ABPV (NC_002548) and SBPV (NC_014137). DWV and VDV‐1 sequences have been analysed previously (Fürst et al. 2014).

Statistical Analysis

Analyses were performed in r v 3.0.2 (R Core Team 2013). RNA virus prevalence differences between pollinator genera were compared in a test of proportions (χ2 test), using a Bonferroni correction (α = 0·003; six species; 15 multiple comparisons) for comparisons between species. So that differences among species of different samples sizes could be meaningfully compared, we estimated true prevalence and 95% confidence intervals using the r library ‘epir’ v0.9‐54, with sensitivity and specificity both set at 95%. Disease prevalence was mapped to sites using the ‘mapplots’ package v1.4, or estimated using Gaussian kernel estimators using the package ‘prevr’ as described previously (Fürst et al. 2014). Distributions of infectious loads were compared using Kolmogorov–Smirnov tests.

To explore possible drivers of RNA virus prevalence in managed and wild bees, we performed generalized linear mixed models (GLMM) with binomial error structure using the package ‘lme4’ v.1.0‐6. Prior to any statistical analysis, we used Moran's I and spline correlograms (package ‘ade4’ v1.6‐2: Dray & Dufour 2007; package ‘ncf’ v1.1‐5: Bjørnstad 2013) to test for potential spatial autocorrelation. The geographical distance between all pairs of sites was calculated, and results indicated there was no significant spatial autocorrelation for any of the RNA viruses in A. mellifera or Bombus spp. (P > 0·05). V. destructor mites have caused both an increase in viral load and prevalence of several RNA viruses in western honeybee populations, including DWV (Martin et al. 2012), ABPV (Genersch & Aubert 2010) and SBPV (Carreck, Ball & Martin 2010; Santillán‐Galicia et al. 2014). We hypothesize that the association of these viruses with A. mellifera has resulted in disease spillover into wild Bombus spp. populations. We therefore modelled Bombus virus prevalence as dependent on A. mellifera virus prevalence, A. mellifera abundance, Bombus abundance, latitude, longitude and landcover type, while treating site and species as random effects. However, to account for uncertainty surrounding the true directionality of pathogen spillover, we also conducted models with A. mellifera virus prevalence as the response variable, retaining all other predictors except species as a random effect. We conducted separate GLMMs for BQCV, DWV and ABPV. SBV and SBPV were not modelled due to insufficient positive samples (n = 4 Bombus and n = 5 A. mellifera individuals, respectively). Site G was removed prior to statistical analysis as no A. mellifera foragers were collected at this site. A. mellifera and Bombus spp. abundance were log‐transformed, and all quantitative predictors were standardized to a mean of zero and standard deviation of one prior to analysis. Models were simplified by backward stepwise selection based on AIC (‘drop1’ function). We used variance inflation factors (VIF) to check for colinearity among our explanatory variables, applying a cut‐off value of 3. Variables with a high VIF were removed one at a time until all VIF values were below 3 (Zuur et al. 2009). Both conditional (r 2 c, all factors) and marginal (r 2 , fixed factors only) values are shown.

In addition to individual GLMMs, we summed the prevalence of each virus at each site and modelled the resulting total virus prevalence data in a general linear model (GLM) to explore the overall relationship of virus prevalence between A. mellifera and Bombus spp. We modelled Bombus prevalence as dependent on A. mellifera virus prevalence, A. mellifera abundance, Bombus abundance, latitude, longitude and landcover type. As before, abundance was log‐transformed and quantitative predictors were standardized. As before, we also conducted a GLM with A. mellifera virus prevalence as the response variable. Model selection was performed using an automatic approach (package ‘glmulti’, Calcagno & De Mazancourt 2010) using the AICc method. Nagelkerke r values are shown.

Results

Data Summary and Virus Composition

Of 792 sampled bees, the following passed β‐actin quality control (for sample sizes, collection times and species composition by site, see Table S2, Supporting information): 92% A. mellifera (n = 237); 100% B. hortorum (n = 30); 100% B. jonellus (n = 1); 90% B. lapidarius (n = 169); 93% B. lucorum (n = 89); 90% B. pascuorum (n = 55); 100% B. monticola (n = 7); 100% B. pratorum (n = 3); and 64% B. terrestris (n = 92). B. terrestris β‐actin was identified to contain a ligation‐site sequence mismatch, and the proportion of samples passing quality control was lower for this species. In future cross‐species comparisons, MLPA probes should be designed for a wider range of housekeeping genes from which uniform markers across bee species can be selected.

For both A. mellifera and Bombus spp., we detected VDV‐1 and DWV/KV from the DWV/VDV‐1/KV complex, but only ABPV from the ABPV/IAPV/KBV complex. DWV complex strains are closely related at proteolytic sites (de Miranda & Genersch 2010) and naturally recombine (Moore et al. 2011). We therefore refer to the DWV complex as ‘DWV’ from hereon. Although ABPV/IAPV/KBV are thought to be distinct viruses, we also refer to the ABPV complex as ‘ABPV’ from hereon, due to the inability to detect either IAPV or KBV in any sample (n = 54 individuals were positively detected in MLPA, of which 47 were positive for ABPV, but none were positive for IAPV or KBV. Samples that were negative for all three qPCR targets (n = 7) could be attributable to qRT‐PCR primer mismatches preventing amplification). Analysis of nucleotide sequences further confirmed the sequence identity of BQCV, ABPV and SBPV in infected A. mellifera and Bombus spp. foragers (Fig. S1, Supporting information). Unrooted trees for SBPV and BQCV sequences are given in Fig. S2 (Supporting information) (ABPV is not displayed as all sequences were identical). SBPV clones from B. pascuorum and A. mellifera are similar or identical, whereas B. terrestris is represented by two diverging haplotypes. For BQCV, clones from each species were more clearly separated, but this is unsurprising given that each individual host bee originated from a different site. Interestingly, B. terrestris was again represented by two distinct haplotypes.

Both MLPA and qRT‐PCR assays did not specifically amplify the negative strand of RNA viruses, and as such, they did not test for actively replicating virus directly. Nevertheless, our methods provide a reliable indicator as to the presence and potential severity of viral infections in bee foragers by employing a multiplexed presence/absence screen followed by quantification.

Prevalence

In an analysis combining all RNA viruses as a single response, the true prevalence was 51% (95% CI: 44%, 58%) in A. mellifera and 23% (95% CI: 19–27%) in Bombus spp. (χ2 1 = 50·0, P < 0·0001). Most viruses occurred singly, with co‐occurrence of two and three viruses being detected in, respectively, 7% (95% CI: 4–12%) and 1% (95% CI: 0–3%) of A. mellifera individuals, and 3% (95% CI: 2–5%) and 0·2% (95% CI: 0–1%) of Bombus spp. individuals (Fig. S3, Supporting information). The proportion of coinfected individuals did not depart from null expectations (A. mellifera: χ2 3 = 0·5, P = 0·918; Bombus: χ2 3 = 4·8, P = 0·189). The most prevalent virus was DWV in honeybees (36%, 95% CI: 30–43%) and ABPV in bumblebees (11%, 95% CI: 8–14%; Table 1). CBPV was not recorded from any sample.

Table 1

Virus prevalence in per cent for six virus targets, with 95% CI in square brackets. Sample numbers for each pollinator genus are shown in parentheses

Pollinator	BQCV	DWV	ABPV	SBPV	SBV	CBPV
A. mellifera (237)	15 [10, 20]a	36 [30, 43]a	5 [2, 9]	2 [1, 5]	2 [1, 4]	0 [0, 2]
Bombus spp. (453)	6 [4, 8]	3 [2, 5]	11 [8, 14]a	5 [3, 7]	1 [0, 2]	0 [0, 1]

Significantly higher virus prevalence in a test of proportions (A. mellifera vs. Bombus spp.).

Five RNA virus targets were detected in both A. mellifera and Bombus spp. (Table 1). In a test of proportions, BQCV (χ2 1 = 13·2, P < 0·001) and DWV (χ2 1 = 126·4, P < 0·0001) were more prevalent in A. mellifera, whereas ABPV was more prevalent in Bombus spp. (χ2 1 = 6·3, P < 0·05). Although SBPV and SBV were more prevalent in Bombus spp. and A. mellifera, respectively, differences between host genera were not statistically significant (χ2 1 = 2·1, P = 0·15; χ2 1 = 0·32, P = 0·57, for SBPV and SBV, respectively). In a comparison of virus prevalence among the five commonest host species (n > 10 collected individuals), we found that DWV and SBPV were significantly more prevalent in A. mellifera and B. hortorum, respectively (Fig. 1). ABPV also occurred at significantly higher prevalence in B. lapidarius compared with A. mellifera, B. lucorum and B. pascuorum, but not B. hortorum or B. terrestris (Fig. 1).

Figure 1

Prevalence of each virus mapped by individual species, showing mean true estimates and 95% CIs. Bonferroni‐corrected chi‐square test for multiple comparisons: *DWV in Apis mellifera significantly higher prevalence compared with all other species. § SBPV in Bombus hortorum significantly higher prevalence compared with all other species. ABPV: significant multiple comparisons indicated by letters a and b. Note that ‘all Bombus’ is not included in statistical comparison (see Table 1).

We mapped the prevalence of both individual and combined RNA viruses by site and pollinator genus (Fig. 2). These indicated that disease prevalence between managed (A. mellifera) and wild (Bombus spp.) bees were linked. In GLMMs of individual viruses, we found that prevalence of BQCV (final model r 2 = 0·28; r 2 = 0·28), ABPV (final model r 2 = 0·51; r 2 = 0·08) and DWV (final model r 2 = 0·39; r 2 = 0·12) in A. mellifera had a positive effect on BQCV, ABPV and DWV prevalence in Bombus spp., respectively (Table 2, Fig. S4, Supporting information), although this effect was marginally not significant in the final DWV model. Additionally, abundance of Bombus spp. and A. mellifera had a negative and positive effect, respectively, on BQCV prevalence in Bombus spp. When we reconstructed models with A. mellifera virus prevalence as dependent on Bombus spp. pathogen prevalence, both the final models and significant predictors closely matched the original models that treated Bombus spp. prevalence as the response variable (Table S3, Supporting information). In a GLM of combined RNA viruses (where prevalence data were summed for all five positively detected RNA viruses), we found that disease prevalence in A. mellifera also had an overall positive effect on disease prevalence in Bombus spp. (final model Nagelkerke r 2  = 0·62, Table 2, Fig. S4, Supporting information), and that latitude was also a significant predictor. Again, when we reconstructed the GLM with A. mellifera virus prevalence as dependent on Bombus spp., the final models matched the original GLM (Table S3, Supporting information). Sites harbouring highest overall RNA virus prevalence were concentrated in SE England. On the other hand, those harbouring lowest RNA virus prevalence were located in remote western regions, including two V. destructor mite‐free islands (Y: island of Colonsay; Z: Isle of Man) that contained the lowest overall disease prevalence across all sites (Fig. 3).

Figure 2

Prevalence of (a) combined and (b) individual RNA viruses mapped by site and pollinator genus (%). Apis mellifera and Bombus spp. are represented as light red and dark blue circles, respectively.

Table 2

(a) Best model explaining individual virus prevalence in Bombus spp. using GLMMs and AIC for model selection. Note that the sign of the parameter estimates for abundance is opposite to the direction of the relationship between variables due to the way in which abundance was measured (see Materials and methods). (b) Best model explaining total RNA virus prevalence in Bombus spp., using a GLM and AICc for model selection

Response (Model)	Virus	Parameters	Estimate	SE	z‐value	P‐value
(a)	BQCV	Intercept	−3·212	0·305	−10·542
		Apis BQCV	0·542	0·186	2·917	0·004a
		Apis abundance	−0·686	0·335	−2·046	0·041a
Bombus virus prevalence		Bombus abundance	0·813	0·345	2·356	0·018a
(GLMM)	DWV	Intercept	−4·185	0·506	−8·275	0·055
		Apis DWV	0·818	0·426	1·918
	ABPV	Intercept	−3·725	0·586	−6·356	0·008a
		Apis ABPV	0·727	0·274	2·654
(b)	ALL	Intercept	0·050	0·063	0·792
Bombus virus prevalence		Apis all viruses	0·367	0·097	3·777	0·001a
(GLM)		Longitude	0·091	0·039	2·350	0·028a

Significant variables.

Figure 3

(a) Combined Apis mellifera and Bombus spp. virus prevalence by site. Mean and 95% CIs calculated from presence/absence individual data (single and multiple infections treated equally). Sites with ≤5% and >50% RNA virus prevalence are emphasized (red and blue, respectively). *island of Colonsay (site Y) and §Isle of Man (site Z). (b) RNA virus prevalence mapped by Gaussian kernel estimation, with site locations overlaid. Bar graphs summarize the prevalence of individual viruses at blue sites, with proportions derived from Bombus spp (white) and A. mellifera (black) indicated. bq = BQCV; dw = DWV; ab = ABPV; sp = SBPV; sb = SBV.

Virus Load

We quantified BQCV, DWV, ABPV, SBPV and SBV from the positively detected A. mellifera and Bombus spp. foragers (Fig. 4). For BQCV, virus loads between A. mellifera and Bombus spp. were not significantly different (two‐sided Kolmogorov–Smirnov, D 9,17  = 0·38, P = 0·63), with putatively low‐level infections (104–106 virus particles) predominating in bees from both genera. For DWV, viral loads in A. mellifera foragers were greater than in Bombus spp. (one‐sided Kolmogorov–Smirnov, D 8,45  = 0·63, P < 0·05), consistent with the presence of high‐level infections in A. mellifera (1010–1011 virus particles) vs. low‐level infections in Bombus spp. (104–106 virus particles). We detected a wide range of ABPV and SBPV virus loads in Bombus spp. (104–1011 ABPV and 105–1011 SBPV particles), but sample sizes were inadequate in A. mellifera (ABPV and SBPV: n = 2 each) to be able to compare distributions between pollinator genera. SBV was not detected in any of the samples by qRT‐PCR (n = 8 individuals positively detected by MLPA).

Figure 4

Comparison of relative frequencies (%) of inferred absolute virus loads in Apis mellifera and Bombus spp. individual foragers.

Discussion

In a comprehensive field survey of managed and wild bee pollinators, we found that common RNA viruses previously associated with honeybees are widespread in bumblebee populations, and that viruses vary substantially in terms of individual pathogen load and population‐level prevalence. Significantly, we demonstrate a positive association in disease prevalence between managed and wild bees, indicating that disease spillover may be an important general aspect of RNA virus epidemiology in bee pollinators.

RNA Viruses are Widespread in Wild Bees

RNA viruses are prevalent in wild bee populations and occur broadly in the landscape. More than one in every five bumblebee foragers sampled contained at least one of the RNA viruses that were screened for. We note that the total impact of RNA viruses on wild bees is likely to be higher than suggested from our prevalence data, as severely affected individuals may be less likely to fly and forage.

Of the targeted viruses, we found that BQCV, DWV, ABPV and SBPV occurred commonly in bumblebee foragers. Previous studies have indicated that wild bees could harbour RNA viruses typically associated with honeybees (Genersch et al. 2006; Singh et al. 2010; Evison et al. 2012), and in the case of DWV, that spillover is likely occurring between managed and wild bee populations (Fürst et al. 2014). By combining a structured survey of multiple RNA viruses with a quantitative analysis of pathogen load, we show that not only DWV but also BQCV is widespread in wild bee populations, and that bumblebee foragers largely harbour low levels of these viruses. In honeybee foragers, we find that BQCV occurs at similarly low levels, but that levels of DWV are significantly higher, as would be expected given the prominent role of V. destructor mites in vectoring this virus. Surprisingly, we find that ABPV and SBPV are more prevalent in bumblebee than honeybee foragers (although this difference is not significant for SBPV), and that bumblebees harbour a wide range of pathogen loads, including a substantial proportion of putatively high infections (>109 virus particles per individual).

ABPV and SBPV were significantly more prevalent in B. lapidarius and B. hortorum, respectively, as compared with most other bee species (including A. mellifera), suggesting that differences in host susceptibility or quality may exist (Ruiz‐González et al. 2012). B. lapidarius is common in England and Wales, while B. hortorum is widespread across Great Britain (Goulson 2009), so there is no clear pattern of pathogen prevalence and bumblebee rarity among our samples. Alternatively, this could indicate that species such as B. lapidarius and B. hortorum are simply more tolerant to infection than others. In addition to factors relating to host immunity, life history parameters linked to phenology, such as the relative abundance of bumblebee foragers and/or reproductives relative to honeybees (the viral titres of which are known to vary temporally, Runckel et al. 2011) may play an important role in between‐species differences in disease prevalence and transmission. Additionally, although V. destructor mites are able to vector both ABPV and SBPV in honeybees, these viruses were found at lower prevalence in honeybee foragers. In the case of ABPV, while there is evidence that the arrival of V. destructor mites has increased the prevalence of ABPV in A. mellifera (reviewed in Genersch & Aubert 2010), the reduced survival of infected pupae could explain why ABPV is detected less frequently than DWV in honeybee foragers (Sumpter & Martin 2004; Schroeder & Martin 2012).

The situation for SBPV is less clear, although both field (Carreck, Ball & Martin 2010) and laboratory experiments (Santillán‐Galicia et al. 2014) suggest that it may be transmitted between honeybees via V. destructor mites and that it may be more virulent than DWV. With respect to wild bees, and bumblebees in particular, virtually nothing is known of the epidemiology of these RNA virus. Furthermore, controlled infection experiments are required to improve basic knowledge of the impacts of these and other RNA virus in non‐Apis bees (e.g. Meeus et al. 2014), and to test competing hypotheses for host species differences in disease prevalence.

Circulation of Viruses between Managed and Wild Bees

We detected a significant association between prevalence of viruses in honeybees and bumblebees, both in a combined RNA virus analysis, and for viruses analysed separately, notably BQCV and ABPV. Interestingly, BQCV prevalence in bumblebees was associated negatively with Bombus abundance but positively with A. mellifera abundance. This might indicate that lower bumblebee abundance is caused by higher BQCV prevalence, itself the result of higher honeybee abundance. However, directions of causality remain equivocal, and given current understanding, we advocate restraint in the extent to which our models are interpreted. Our data also indicated an association between pollinators for DWV, but the relationship was not as strong as the effect detected in a previous study (Fürst et al. 2014). Several factors could be responsible for this variation. First, Fürst et al. (2014) based their analysis on a data set comprised of a random subsample of 10 individuals of the two commonest species, whereas in this study, we included every successfully amplified sample at each site regardless of species. Secondly, many individuals were differentially excluded based on separate quality control measures, which resulted in substantially reduced representation of B. terrestris in our study. Finally, the sensitivity of MLPA, a multiplexed approach based on competitive PCR, is lower than a single RT‐PCR approach, resulting in a higher likelihood of false negatives (de Miranda et al. 2013). Given the last consideration, it is probable that the impact of RNA viruses on bee populations is greater than we report.

Several outstanding questions emerge from our findings that merit further attention, given the pressures faced by bumblebees globally (Fitzpatrick et al. 2007; Bartomeus et al. 2013) and the potential role of pathogens in these declines (Cameron et al. 2011; Meeus et al. 2011). First, the direction of disease spillover between managed and wild bees represents a major unanswered question. By extension, whether honeybee or bumblebee populations are more important as natural reservoirs of RNA virus infections also remains unclear. Based on prevalence alone, BQCV and DWV appear more closely linked with honeybees whereas ABPV and SBPV are more common in bumblebees. Unfortunately, prevalence levels taken at a single time point are not informative with regard to understanding the primary reservoir host(s), or by extension, the principle direction of disease transmission. Spillover may be bidirectional, with both pollinator genera functioning as suitable long‐term reservoirs (although the perennial life cycle and presence of V. destructor mites may favour the honeybee as a more likely long‐term disease source). We also cannot exclude the possibility that one or more unknown species from the wider invertebrate community is the primary disease reservoir. Equally, the main reservoir may be a low‐prevalence host, whereupon entering a second host the virus spreads rapidly, resulting in an epidemic and higher observed prevalence. For example, it is plausible that increased ABPV prevalence in bumblebees is due indirectly to the increased exposure to infectious virus particles emerging from V. destructor‐infested honeybee colonies that contain higher than normal loads of ABPV (Genersch & Aubert 2010). Interestingly, we found that the sites least affected by disease in this study are also those where V. destructor has not yet become established in honeybees, although these are also the most remote island locations.

Our study significantly extends previous findings that suggested horizontal transmission of pathogens between bee pollinators. However, demonstration of the mechanistic basis of host switching in the field is still lacking. Potential transmission routes include direct contacts between bees (bumblebees entering managed honeybee colonies are not uncommon, for example, Genersch et al. 2006), or more likely, via indirect interactions such as through shared use of floral resources (McArt et al. 2014). Unfortunately, the field is largely devoid of observational or experimental data that tackle the issue of pathogen transmission at flowers directly (although see Durrer & Schmid‐Hempel 1994).

Singh et al. (2010) report that pollen pellets carried by honeybee foragers (in addition to stored honey) contain viable RNA virus, despite absence of virus in the forager itself. This suggests that infective inocula reside at flowers and may be collected by flower visitors. However, the probabilities of infectious material being deposited and subsequently acquired by a new host while remaining viable are unknown. As outlined previously (McArt et al. 2014), we expect traits such as flower complexity (Durrer & Schmid‐Hempel 1994), pollinator flower‐handling time and floral secondary compounds [e.g. antimicrobial compounds, host immune modulators (Mao, Schuler & Berenbaum 2013)] to influence the likelihood of infection. Pathogen transmission at flowers remains poorly understood, but bridging this gap in understanding should be a priority for pollinator research.

Emerging Environmental Pressures on Wild Bees

EIDs represent one of several stressors that have been implicated in bee pollinator declines. Other major drivers are thought to include habitat change and loss (reviewed in Potts et al. 2010), and more recently, pesticides (Desneux, Decourtye & Delpuech 2007; Gill, Ramos‐Rodriguez & Raine 2012; Whitehorn et al. 2012; Williamson & Wright 2013). However, an explicit causal link between any single factor and bee declines has not emerged (Vanbergen et al. 2013). It is plausible that several factors acting in synergy serve to amplify pressures on pollinators (González‐Varo et al. 2013), or that a range of different factors may produce similar levels of stress at the colony level (Bryden et al. 2013).

Among the bumblebees in this study, we found a significant proportion of the active foraging workforce (>20%) to harbour RNA viruses, often at putatively low levels. In combination with other stressors, such pervasive disease pressures could have a general and profound impact on the long‐term health of bee populations. Recent studies have reported interactions between sublethal doses of neonicotinoid pesticides and pathogens, with significant impacts on virus replication and host immunity (Di Prisco et al. 2013) and bee mortality (Doublet et al. 2014). Alongside the lethal impacts of environmental stressors (either acting in isolation or in combination), the role of sublethal chronic stress has also attracted attention. Both pathogens (Mayack & Naug 2009) and pesticides (Gill, Ramos‐Rodriguez & Raine 2012) can act as chronic stressors, with negative impacts on social bee colony function, leading eventually to colony failure (Bryden et al. 2013).

Our findings reveal the widespread prevalence in wild bee populations of multiple RNA viruses previously associated with honeybees. We present evidence of ongoing or recent transmission of viral diseases between managed and wild bee populations, but we remain cautious in speculating on the main direction of spread between pollinator genera, or in making predictions about which bee species act as the principle reservoir for infectious disease. The arrival of V. destructor mites heralded a major shift in the epidemiology of several RNA viruses in the western honeybee, with potentially wide implications for disease spillover among wild pollinators. While we show that RNA viruses are widespread in wild bees, it is unclear to what extent viral challenge impacts bumblebees at the population level, either in isolation or in combination with other stressors.

Table S1. List of qRT‐PCR primers used in this study

Click here for additional data file.

Table S2. Detailed sample information by collection site

Click here for additional data file.

Table S3. GLMM and GLM final models with A. mellifera virus prevalence as the response variable

Click here for additional data file.

Fig. S1. BQCV, ABPV and SBPV cloned sequence alignments.

Click here for additional data file.

Fig. S2. Phylogenetic trees of cloned BQCV and SBPV sequences.

Click here for additional data file.

Fig. S3. Summary prevalence (%) of RNA viruses in A. mellifera and Bombus spp.

Click here for additional data file.

Fig. S4. Linear regressions of raw A. mellifera and Bombus spp. RNA virus prevalence data.

Click here for additional data file.