MATERIAL AND METHODS
Study area and design
From May to June 2017, the study was conducted in Central Valais,
Switzerland. The landscape is dominated by vineyards interspersed with
small patches of dry oak stands, steppe and orchards in the plains.
About 70-80% of the vineyards are intensively managed and support
virtually no ground vegetation due to regular herbicide application,
whereas the remaining 20-30% are cultivated by farmers who have adopted
more environmentally friendly management practices, promoting ground
vegetation . Forty B. terrestris terrestris colonies, purchased
from Andermatt Biocontrol, were allocated to vineyard fields of varying
ground-vegetation coverage, landscape composition and configuration in a
semi-experimental factorial approach (Fig. 1; see Maurer et al. 2020).
Bumble bee colonies
Two-week-old colonies of B. terrestris terrestris were weighed as
a proxy of initial colony size and then randomly allocated to
experimental fields, where they remained from 08-May-2017 to
23-June-2017. Weight gain of colonies, number of workers, queens and
total number of pupal cells (as a sum of healthy, dead and hatched pupal
cells, pollen and nectar pots), as well as infections by parasitic
moths, Aphomia sociella , were assessed. For details, please refer
to Maurer et al. (2020). Workers (N = 5) were collected from 35 colonies
at the start and at the end of the experiment, stored at -20°C until the
transfer to the laboratory and then stored at -80°C until further
processing.
Environmental predictors
To assess the role of the composition and configuration of major
land-use types around the experimental fields, we used the Swiss
landcover data SwissTLM3D (copyright@swisstopo 2020, resolution 1-3m).
As not all municipalities were covered by this data, we additionally
used the CORINE landcover data from 2018 to fill the gaps (© European
Union, Copernicus Land Monitoring Service 2020, European Environment
Agency EEA). We then aggregated the landcover classes (Fig. 1, see SI
section 1.1 for details on class aggregation) and calculated the
relative area covered by farmland, forests, residential areas, and
vineyards within 9 buffer zones ranging from 100 m up to 1’500 m radii
in 100 m increments (100–500m) and 250 m increments (500–1’500m),
respectively (see SI Fig. S1). To do so we used the metric percentage of
landscape (PLAND, in the R package landscapemetrics (Hesselbarth et al.
2019). We also calculated the mean distance between patches of the same
landcover class as a measure of configuration (mean of Euclidean
nearest-neighbor distance ENN_MN). Lastly, we obtained an estimate of
landscape heterogeneity for the same 9 spatial scales by calculating a
Shannon diversity index based on all eight reclassified land-cover
classes, taking both the number of land-cover classes and the abundance
of each class into account (metric Shannon diversity SHD). We
additionally calculated habitat area (PLAND) and fragmentation (patch
density, PD) of vegetated vineyard fields specifically within the same 9
spatial scales as it has been shown to be an important landscape feature
for arthropods, including bumble bees .
To investigate field-scale effects of vineyard management on virus
susceptibility of bumble bees, the ground vegetation density and flower
resources were measured during the period of experimentation in a random
subset of 6 vegetated vineyard fields surrounding the experimental field
within a buffer zone of 250 m radius (see . Additionally, we included
the field size of the experimental field as an explanatory variable
(Table 1).
To address the influence of honey bee hives in relation to virus
transmissions in the surroundings of the bumble bee colonies, we
included the nearest-neighbor distance to the next honey bee hive and
the number of honey bee hives within the 9 spatial scales. All these
metrics were calculated in R
Virus analyses
RNA extraction – For the virus identification RNA was extracted from
bumble bees (N = 5 per colony) collected at the beginning (day 0) and at
the end of the experiment (day 45) using a NucleoSpin® RNA II kit
(Macherey-Nagel, Oensingen, Switzerland) following the manufacturer’s
recommendations. Individual bees were crushed in PBS Buffer (0.5 mg
tissue/μL) with a 5-mm metal bead in 2-mL Eppendorf® (Basel,
Switzerland) tubes using a Retsch® (Haan, Germany) MM 300 mixer mill for
1 min at 25 Hz . Fifty μL of the homogenate was used for the extraction
and RNA got stored in 60 μL elution buffer .
Next generation sequencing – RNA from every colony was used for the
identification of bee viruses present in the bumble bee colonies.
Briefly, after QC evaluation, RNA from each individual was pooled
according to the two sampling sessions (day 0 and day 45). Libraries
comprising each sampling session were prepared using the Corall total
RNA-Seq Library kit (Lexogen, Vienna, Austria). NGS were performed using
an Illumina SP flow cell (100 Mio reads/pool, 300 cycles) in paired-end
mode (2 × 150 bp).
Bioinformatics analysis – Reads were quality-trimmed with fastp Version
0.12.5) and mapped to the Bombus terrestris host genome (iyBomTerr1.2,
ncbi bioproject PRJEB45694) using STAR Version 2.7.3a). Quality-trimmed
and unmapped reads were assembled via SPAdes Version 3.14.0). Resulting
scaffolds were then aligned to virus nucleotide and protein and sequence
databases (Genbank and RefSeq viral nucleotide sequences downloaded on
21-01-2021, UniProt viral amino acid sequences downloaded on 21-01-2021)
using BLASTn , Version 2.0.0+, default settings) and DIAMOND , Version
0.9.18, default settings). To exclude false positives, the scaffolds
with a virus hit were aligned to an in-house non-viral database
consisting of archaeal, bacterial, fungal, mammal, and protozoal
sequences. Scaffolds were considered false positive if they had a longer
hit on a sequence of the in-house database compared to the virus
databases or if they had a nucleotide hit of more than 10% of their own
length to any sequence of the non-viral database.
Reverse transcription and quantitative PCR (RT-qPCR) – cDNA synthesis
was performed for each sample (pooled RNA from 5 individuals from each
colony) using a M-MLV RT Kit (Promega, Dübendorf, Switzerland). In
brief, a thermocycler (Biometra, Analytik Jena, Jena, Germany) was used
to incubate 0.75 μL of a random hexamer oligonucleotide (100 μM) and H2O
for 5 Minutes at 70 °C with 0.5 µg of template RNA. Then, 5 μL of 5x
buffer, 1.125 μL dNTPs (10 mM) and 1 μL reverse transcriptase (M-MLV)
were added followed by incubation at 37 °C for 60 min. For the virus
quantification, the group of bee viruses with the highest number of
reads from the above-described libraries were selected (SI Table S1).
qPCR reactions were prepared using 6 μL of 2X reaction buffer
(SensiFAST™ SYBR® No-ROX Kit, Meridian Bioscience, London, UK), 0.24 μL
forward and reverse primer (SI Table S2), 2.52 μL water and 3 μL of
ten-fold diluted cDNA. The qPCR reactions were performed in a CFX96TM
Real-Time PCR Detection Systems (BioRadⓒ, CA, USA) with the following
conditions: 3 min for 95°C, 40 cycles of 3 sec at 95°C and 30 sec at
57°C. The amplification was followed by a melting curve analysis of the
strand dissociation to verify product specificity. The analysis was
performed by reading the fluorescence at 0.5 °C increments from 55 to 95
°C. Each sample was run in duplicate for each of the targeted virus and
the Rps5 reference gene. Furthermore, each plate was run with a ten-fold
serial dilution of purified PCR products that served as standard curves
and two no-template negatives .
Statistical Analyses
We used the viral load of the nine screened viruses to calculate total
viral load (number of virus genome copies). Viral loads were first log
transformed and then summed up per colony and sampling session (day 0
and day 45) to calculate the total viral load. To answer the first
research question regarding how the viral load, richness and viral loads
of the viruses changed from before to after the field exposure, we ran
linear models with total viral load, virus richness, and viral loads of
the six most abundant viruses (Deformed wing virus-B DWV-B, Black queen
cell virus BCQV, Castleton burn virus CBV, Mayfield virus-1 MV 1, Duke
bunyavirus DBV, Bombus cryptarum densovirus BcDV) as response variables
against the session (i.e. day 0 versus day 45).
To test the effects of initial viral loads on colony development
parameters (research question ii) we used initial viral load (day 0),
and viral loads of the same six most abundant viruses in
single-predictor linear models. The following eight colony development
parameters were used: number of workers, number of queens, number of
total cells, number of hatched pupal cells, weight gain, moth
infestation index, and number of parasitic larvae or pupae.
Next, the influence of landscape structure on the virome change was
assessed (research question iii). For this purpose, change in virus
loads was calculated as\(\operatorname{(log(}{{viral\ load)}_{day\ 45}-\log{{(viral\ load)}_{day\ 0})}}\)for total virus load or for each of the six most abundant viruses
separately. Similarly, for the change in virus richness, we calculated
the number of viruses present at day 45 – the number of viruses present
at day 0. For the turnover and appearance of viruses during the field
exposure we used the function ‘turnover’ (package codyn, .
For each of those response variables, we first ran single-predictor
models for each explanatory variable separately (see Table 1), using
simple linear models given the normal distribution of residuals. For all
explanatory variables, those with p<0.1 were combined in a
full model. For predictors that were measured at multiple spatial
scales, one model per scale was run. If (near-)significant effects
(p<0.1) were found, the best scale based on lowest AIC was
selected for the full model. Thus, variables could enter the full model
at different scales – following a multi-scale analytical framework .
Before running the full model, collinearity among predictors was checked
with a threshold of Pearson<0.6. For collinear variables, the
one with higher AIC in the single-predictor model was dropped from the
full model. For the full models, we did a stepwise backwards selection
until only (near-)significant variables were left in the model. Among
those candidate models, the best one was selected based on best model
performance, which was checked using ‘compare_performance’ (package
‘performance’ . Additionally, we correlated colony development
parameters with the field and landscape variables (see detailed methods
and results in SI section 3). Model assumptions such as normality of
residuals, homoscedasticity, and outliers were checked (package
‘performance’). Spatial autocorrelation of model residuals was tested
for each model using the function ‘testSpatialAutocorrelation’ (package
DHARMa, All data and code are available from zenodo online repository
(Bosco et al. 2023).