Contrasting impacts of a novel specialist vector on multihost viral pathogen epidemiology in wild and managed bees

Abstract Typically, pathogens infect multiple host species. Such multihost pathogens can show considerable variation in their degree of infection and transmission specificity, which has important implications for potential disease emergence. Transmission of multihost pathogens can be driven by key host species and changes in such transmission networks can lead to disease emergence. We study two viruses that show contrasting patterns of prevalence and specificity in managed honeybees and wild bumblebees, black queen cell virus (BQCV) and slow bee paralysis virus (SBPV), in the context of the novel transmission route provided by the virus‐vectoring Varroa destructor. Our key result is that viral communities and RNA virus genetic variation are structured by location, not host species or V. destructor presence. Interspecific transmission is pervasive with the same viral variants circulating between pollinator hosts in each location; yet, we found virus‐specific host differences in prevalence and viral load. Importantly, V. destructor presence increases the prevalence in honeybees and, indirectly, in wild bumblebees, but in contrast to its impact on deformed wing virus (DWV), BQCV and SBPV viral loads are not increased by Varroa presence, and do not show genetic evidence of recent emergence. Effective control of Varroa in managed honeybee colonies is necessary to mitigate further disease emergence, and alleviate disease pressure on our vital wild bee populations. More generally, our results highlight the over‐riding importance of geographical location to the epidemiological outcome despite the complexity of multihost‐parasite interactions.

disease persistence and transmission to sympatric hosts (Kilpatrick, Daszak, Jones, Marra, & Kramer, 2006). The addition of a new transmission route, such as a vector, into these complex multihost transmission cycles can have important implications to disease dynamics, and the vector's behaviour and population dynamics may need to be taken into account (Dobson, 2004). By parasitising multiple hosts with varying transmission potential, a generalist vector can dilute the transmission potential of a key host, as demonstrated by the effect of increased mammal diversity on Lyme disease incidence: Higher mammalian host diversity reduced the disease transmission role of the white-footed mouse, the most competent reservoir host of Borrelia burgdorferi (the causative agent of Lyme disease) (Kilpatrick et al., 2006), thereby reducing Lyme disease risk (LoGiudice, Ostfeld, Schmidt, & Keesing, 2003). In contrast, a specialist vector can increase the transmission potential of a key host; for example, the American Robin is believed to be responsible for the West Nile virus epidemic in New York due to preferential feeding behaviour of a mosquito vector on this relatively rare but highly competent host (Kilpatrick et al., 2006). Thus, the epidemiology of pathogens across multiple host species, may be critical in driving disease emergence, and ultimately for disease control.
Pollinators are host to a large number of RNA viruses that are known to be pathogenic to honeybees, but have more recently been identified as multihost pathogens, prevalent in wild bee populations (Evison et al., 2012;Fürst, McMahon, Osborne, Paxton, & Brown, 2014;Levitt et al., 2013;Manley, Boots, & Wilfert, 2015;Manley et al., 2019;McMahon et al., 2015;Singh et al., 2010). There are known differences in viral prevalence and viral load across honeybee and bumblebee hosts for many RNA viruses (McMahon et al., 2015). Deformed wing virus (DWV) is undergoing a global epidemic in honeybees and is an emerging disease in bumblebees (Fürst et al., 2014;Manley et al., 2019;Wilfert et al., 2016); its prevalence and viral load is highest in honeybees (Apis mellifera). Black queen cell virus (BQCV) is also closely linked with honeybees (McMahon et al., 2015) but with high prevalence and viral load found across bumblebee species, particularly when apiaries are present in the area (Alger, Burnham, Boncristiani, & Brody, 2019). Slow bee paralysis virus (SBPV) on the other hand shows higher prevalence and viral load in bumblebee species than in honeybees (McMahon et al., 2015). Alongside DWV, SBPV has been assigned to the genus Iflavirus (de Miranda et al., 2010), while BQCV is a member of the Dicistroviridae family. Deformed wing virus is linked to high overwinter mortality of honeybee hives and increased worker mortality in bumblebees (e.g., Berthoud, Imdorf, Haueter, Radloff, & Neumann, 2010;Dainat, Evans, Chen, Gauthier, & Neumann, 2012;Fürst et al., 2014;Genersch et al., 2010;Highfield et al., 2009;Natsopoulou et al., 2017). Slow bee paralysis virus has been shown to infect Bombus terrestris, significantly reducing longevity under nutritional stress (Manley et al., 2015).
RNA viruses differ in their association with V. destructor (hereafter referred to as Varroa for simplicity), which will influence the risk of virus emergence in wild bumblebees. There is currently no clear evidence associating BQCV with transmission by Varroa (Locke, Forsgren, Fries, & De Miranda, 2012;Ribière, Ball, & Aubert, 2008;Tentcheva et al., 2004), although one study found a weak correlation of BQCV titre with Varroa infestation rates in New Zealand (Mondet, De Miranda, Kretzschmar, Le Conte, & Mercer, 2014). Slow bee paralysis virus differs from this pattern: Varroa has been shown experimentally to be capable of transmitting SBPV (Santillán-Galicia, Ball, Clark, & Alderson, 2014), a virus which is more prevalent in Varroa-positive colonies (Carreck, Ball, & Martin, 2010). However, in the wild, SBPV has been found at higher prevalence in certain bumblebee species (specifically in Bombus hortorum) than in A. mellifera (McMahon et al., 2015), suggesting that A. mellifera is not the reservoir host for this virus.
Varroa has invaded the entire European mainland with the exception of several island refuges off the coast of the British Isles and the French coast. Using single-molecule RNA sequencing, we take advantage of this natural experiment to examine how Varroa has impacted on the diversity and composition of RNA viromes across honeybee and bumblebee hosts. Further, we focus on the epidemiology of two multihost RNA viruses, SBPV and BQCV -two viruses with different apparent patterns of host specificity -to examine whether these differences affect the epidemiological history of these viruses and whether this ultimately results in different scenarios for how the acquisition of a specialist virus affects pathogen transmission.

| MATERIAL S AND ME THODS
We collected foraging bees (355 A. mellifera, 281 Bombus pascuorum, 640 B. terrestris and 38 Bombus lucorum individuals) within a 1 × 1 km area as described in detail in (Manley et al., 2019), from four Varroafree islands; three Varroa-positive islands; and five Varroa-positive mainland sites ( Figure S1, Table S1). We differentiated between B. terrestris and B. lucorum via a mtDNA length polymorphism (Table   S2). We extracted DNA from homogenised gut tissue using Chelex following the manufacturer's instructions, and RNA using Trizol and bromo-chloropropane from individuals (gut homogenate and half the head and thorax (bisected laterally)) as described in detail in Manley et al. (2019). RNA was resuspended in 100 μl (A. mellifera) or 400 μl (Bombus species) of nuclease-free water. cDNA transcription was performed on 2 µl of resuspended RNA using GoScript Reverse Transcriptase (Promega), with random hexamer primers and RNasin to prevent RNA degradation. To determine viral prevalence of BQCV and SBPV we carried out PCR in 20 µl reactions using GoTaq DNA Polymerase (primers and programs in Table S2) and ran the products on 1.5% TAE agarose gel with ethidium bromide nucleic acid staining solution.
We randomly selected 10 BQ CV and SBPV positives per site/species (or total number if there were fewer than 10 positives available) for qRT-PCR analysis (BQCV: n = 148, SBPV: n = 141). We measured RNA concentration (Qbit Fluorometer) and quality (Nanodrop 2000 spectrophotometer) for all samples; all samples had a 260/280 nm ratio between 1.8 and 2.1. We performed cDNA transcription on 2 µl of 400 ng RNA template using GoScript Reverse Transcriptase and diluted the cDNA 1:10 prior to qRT-PCR. Duplicate reactions were run for each sample on a Strategene machine (Mx3005P) using GoTaq qPCR Master mix for dye-based detection (Promega ,   Table S3), alongside two no-template negative controls. We calculated viral copy number using duplicate eight-point standard curves of plasmid DNA, of known quantity, in a 1:10 serial dilution on each plate. We generated BQCV and SBPV plasmids using Promega pGEM-T Easy Vector to clone a 257 bp fragment of ORF 2 of the BQCV genome and 186 bp fragment of the VP2-gene in SBPV from purified PCR products (primer details in Table S3), selecting successful transformants via blue/white screening. Plasmids were extracted using GeneJET Plasmid Miniprep Kit (ThermoFisher).
We used M13 primers (designed to sequence inserts inside pGEM-T Easy Vector: forward 5′-GTTTTCCCAGTCACGAC-3′, reverse 5′-CAGGAAACAGCTATGAC-3′) to confirm the correct product had been cloned. We linearlised the plasmids using the restriction enzyme Apa 1 (New England Biolabs), according to the manufacturer's instructions, and diluted them 1:1,000 with nuclease-free water.
We ran Jmodeltest to compare and select an appropriate evolutionary substitution model for each alignment based on the Bayesian Information Criterion (Alizon & Fraser, 2013) for use in phylogenetic reconstructions (Table S4). We fit discrete trait models with asymmetric substitution models for host species and geographic location, which allows transitions to and from a host or location to occur at different rates (trait rate and indicators operators weight = 1), implemented in Beast v1.8 (Drummond & Bouckaert, 2015). We ran and compared nine models concurrently for each alignment with different demography and molecular clock rates, and used the path sampling maximum likelihood estimator, implemented in Beast 1.8, to determine the best model (Table S5). For BQCV, the preferred model used a HKY+I+G substitution model with an exponential relaxed molecular clock and a constant population size prior. For SBPV, the preferred model used a HKY+G substitution model, with a lognormal relaxed molecular clock and a constant population size prior.
Models were run without prior knowledge of the evolutionary rate.
We excluded B. lucorum samples (n = 38) and the single B. pascuorum from Quiberon from prevalence analyses because of low sample size. We used RStudio (v0.99.896) for all statistical analyses. To account for imperfect detection of viruses by PCR assays, we calculated true prevalence with 95% confidence intervals (using R library epiR v0.9-82 and the function epi.prev, with the confidence level conservatively set at 95% [Reiczigel, Foldi, & Ozsvari, 2010] and confidence intervals calculated based on methods in Blaker (2000)). To test if BQCV and SBPV prevalence was affected by Varroa-presence we ran generalised linear mixed models (GLMMs) lme4 package (v1.1-12) (Bates, Maechler, Bolker, & Walker, 2015) with binomial error distribution and logit link function. In the full models we included three-way interaction between fixed effects; Varroa-presence/absence, species (a factor with three levels: A. mellifera, B. terrestris and B. pascuorum) and island/ mainland location; latitude and sunshine hours duration were included as proxies for favourable disease transmission as additional fixed effects (Fürst et al., 2014;Manley et al., 2019): field site and individual were included as random effects (individual was added to account for over-dispersion in the models [Harrison, 2014]). We removed nonsignificant terms using the 'ANOVA' function to determine the minimum adequate model (MAM). We compared the MAM with the null model (which only included random effects) using ANOVA to test the full effect of our predictors and examined residual plots to assess model fit. We ran GLMMs to examine if viral load was affected by Varroa-presence: BQCV and SBPV were tested in separate models with Gamma error distribution and inverse link function. Viral load data were log transformed before analysis because the data varied across orders of magnitude from 10 3 to 10 10 .
We sampled on Varroa-positive islands, as well as paired Varroapositive mainland sites, to test for a possible island effect on disease prevalence. Further, we ran models on reduced data sets to rule out the effect of island location on disease prevalence (a)  Varroa-presence affects rates of coinfection.

| Virome diversity
PacBio single molecule RNAseq data suggest that geographic location, rather than host species or Varroa presence/absence, determines the overall viral composition of pollinator populations as well as the viral diversity. A principal component analysis of viral F I G U R E 1 Principle component analysis on the viral composition of eight bee populations, based on single molecule real-time (SMRT) sequences mapped by BWA-mem to all previously sequenced "honeybee" viruses and fragments from 23 newly discovered bumblebee viruses (Table S6). Apis mellifera populations are shown with white circles and Bombus terrestris populations shown with black circles. 38.8% of the variance is explained by PC1, and the cumulative % variance explained by PC1 and PC2 is 62.7% PacBio reads (Figure 1) shows that populations from Liverpool and the Isle of Man (located in the North of England and the Irish Sea, respectively; see Figure S1) are clearly separated from  (Gauthier et al., 2015), is also present in all populations across host and site (6.8% of reads), although highest in A. mellifera populations from Ushant and Isle of Man and relatively rare in B. terrestris populations (Table S7). Slow bee paralysis virus (Harpenden and Rothamsted strains) (12.9% of reads), Sacbrood virus (4.4% of reads) and BQCV (0.5% of reads) are notably common in the Liverpool/Isle of Man populations but mostly absent from the Le Conquet/Ushant populations (for map of sites see Figure S1); reads from all other viruses (Table S6) are rare (including four newly discovered bumblebee viruses [Pascall et al., 2018]), or absent across all populations. We found no evidence for further novel viruses within these pools from analysis of contigs assembled from reads that did not map to either host or viral genomes.
In GLMMs, host species was a significant predictor of BQCV prevalence, but not for SBPV (Table 2), and was thus removed from the model by model selection (SBPV χ 2 = 2.42 (df = 2), p = .30). Varroa presence was a significant predictor for BQCV and SBPV prevalence (Table 2, Figure 3). Varroa presence predicts a 10-fold increase in BQCV prevalence in A. mellifera and a 20-fold increase in both B. terrestris and B. pascuorum; and an increase in SBPV prevalence across host species by over 100-fold (Table S8) To further confirm that Varroa presence, rather than an island effect, explained viral prevalence, we ran GLMMs on a reduced data set excluding Varroa-free sites (i.e., comparing three Varroa-present islands with five Varroa-present mainland sites) and with island as an explanatory variable: there was no significant difference in BQCV or SBPV prevalence on Varroa-present islands compared to Varroapresent mainland (estimate ± SE of the fixed factor "island/mainland" in the model = 1.17 ± 1.10, p = .28; 2.3 ± 2.1, p = .26, respectively). In addition, we ran the same models on a data set excluding mainland sites (i.e., comparing four Varroa-free islands to three Varroa-present islands), here Varroa presence remained a significant explanatory factor of BQCV and SBPV prevalence, with higher prevalence on Varroapresent islands compared to Varroa-free islands (estimate ± SE of the fixed factor Varroa presence in the model = 3.76 ± 1.49, p = .018; 7.19 ± 3.87, p = .063, respectively). It is notable that SBPV was found at extremely high prevalence on two Varroa-present islands, Guernsey and Jersey, but is absent on Belle-Ile (the third Varroa-present island).
However, prevalence is also high and absent on the closest mainland sites, Cherbourg and Quiberon, respectively, indicating a location difference rather than an island effect.
Host species is a significant factor predicting SBPV viral load (GLMM results reported in Table 3). Bombus pascuorum has significantly higher SBPV loads than A. mellifera (Figure 4, KS-test: D = 0.71, p < .001), with the GLMM predicting SBPV loads in B. pascuorum to be one order of magnitude higher than in A. mellifera and B. terrestris. There were no significant differences in viral load between host species for BQCV (KS-test: D = 0.18, p = .37, GLMM results in Table 3). Although Varroa influences prevalence of both BQCV and SBPV (Table 2, Figure 3), it is not a significant predictor of viral load (Table 3, Figure 4). Latitude, duration of sunshine hours and island/mainland had no influence on viral load and were removed from the models by model selection.   [Hudson et al., 1992;Hudson, 2000] (Figure 6a,b).  (Table   S9a). For SBPV, there is strong support for transmission between all species, but strongest support for transmission between A. mellifera and B. pascuorum, and vice-versa (BF = 30.14, BF = 28.78, respectively, Table S10a). There is no clear pattern of transmission routes between locations; for example, locations nearest each other did not have the strongest support for state transition rates (Tables S9b and S10b, Figure S3; MCC trees for BQCV and SBPV see Figure S4a,b).  However, prevalence of viruses significantly varies across species and viral pathogen. A. mellifera has previously been implicated as the reservoir host for DWV-A and DWV-B (Fürst et al., 2014;Manley et al., 2019), and results of the present study suggest  (Dobson, 2004;Keesing, Holt, & Ostfeld, 2006). While it is known that SBPV and DWV affect worker mortality in B. terrestris, such virulence data is lacking for B. pascuorum as well as for BQCV. To assess the potential impact of these pathogens on pollinator health, it will be essential to experimentally measure virulence in a representative range of insect pollinators. It is important to note that Crithidia bombi, which shows relatively small condition-dependent effects on worker mortality (Brown, Loosli, & Schmid-Hempel, 2000), reduces the fitness of eusocial bumblebees by up to 40% when examining this trait at a colony-level (Brown, Schmid-Hempel, & Schmid-Hempel, 2003;Yourth, Brown, & Schmid-Hempel, 2008 (Mondet et al., 2014). Further, a previous study found no effect of Varroa on SBPV prevalence (Martin et al., 2012), despite associations with Varroa transmission (Carreck et al., 2010;Santillán-Galicia et al., 2014).

| D ISCUSS I ON
Our results suggest that the role of Varroa in the epidemiology of these viruses needs to be reassessed.

Crucially, BQCV and SBPV viral load is not increased by Varroa
presence, suggesting that Varroa does not increase the initial titre by replication or bioaccumulation; rather these results suggest that Varroa can passively vector the virus between individual A. mellifera.
One possible explanation is that the natural titre is high enough to cause maximum effect without amplification by Varroa, but this is unlikely given the broad range of viral titres found across naturally infected individuals (10 3 -10 11 viral particles). It is possible that the indirect effects of Varroa infection such as immunosuppression (Nazzi et al., 2012)  In contrast to the low genetic variation and exponential expansion All evidence combined suggests that SBPV and BQCV are relatively stable infections compared to DWV-B, which has recently emerged and is undergoing epidemic growth (Manley et al., 2019). There is no clear pattern of transmission routes between locations, for example locations nearest each other did not have the strongest support for routes of transmission, which gives additional support to the hypothesis that these are long term stable virus populations with strong geographic structure, obscuring any pattern of geographic movement.
In combination, our results suggest that Varroa presence increases intraspecific BQCV and SBPV prevalence in A. mellifera by direct but passive transmission, and possibly opportunistic infection of individuals weakened by Varroa or Varroa-DWV infection: This increases spillover to wild bumblebees, resulting in increased prevalence across hosts, but has not caused epidemic outbreaks of these two viruses. This is in stark contrast to DWV-B, where Varroa dramatically increases DWV prevalence and viral loads in A. mellifera, leading to increased spillover of high viral loads to competent bumblebee hosts, resulting in an ongoing epidemic across hosts (Manley et al., 2019). It is clear that Varroa plays a complex role in facilitating disease emergence in wild bumblebees. Controlling the Varroa mite infection in managed honeybees is vital to prevent further impact of viral disease in wild bees, already under pressure from habitat loss and pesticide use. These results demonstrate how a specialist vector can increase transmission potential for different multihost-pathogens; the ultimate outcome -a global epidemic and disease emergence for DWV (Fürst et al., 2014;Wilfert et al., 2016) versus increased viral prevalence without increased pathogen load or epidemic growthmay depend on the specifics of host-pathogen interactions.

AUTH O R CO NTR I B UTI O N S
The study was designed by L.W., and R.M., with input from M.B.
R.M. performed field and laboratory research. B.T. contributed bioinformatics analysis of SMRT data. R.M. analysed the data with input from L.W. The paper was written jointly by R.M., and L.W., with input from all authors.

DATA AVA I L A B I L I T Y S TAT E M E N T
BQCV and SBPV sequences have been uploaded to Genbank under embargo, with accession numbers MG265504-MG265572 and MG265573-MG2656050, respectively, and will be released on pub-