Study design

In 2013, a total of 16 fields (8.9 ± 5.4 ha; mean ± s.d.) in southern Sweden, intended for spring-sown oilseed rape (Brassica napus L.) were paired according to geographical proximity (but separated by >4 km) and land use (Fig. 1, see above). The surrounding landscape was inspected for the absence of flowering crops. However in 2013, two fields remained in the study even though another oilseed rape field was present nearby, so as to retain as many farm-pairs as possible34. In each farm pair, one field was randomly assigned to be sown with clothianidin-treated oilseed rape seeds, while the other field was sown with seeds not treated with clothianidin (treated: 8; control: 8). The same paired farms were used in 2014 but with the treatments reversed, i.e., locations with treated fields in 2013 had control fields in 2014 and vice versa (treated: 6; control: 4). Due to crop rotation, different fields within each farm were used in 2013 and 2014 (Fig. 1, see above). In order to create Fig. 1, we downloaded a map from the World Borders Database (downloadable here: http://thematicmapping.org/downloads/world_borders.php).

Information on surrounding landscape variables for the different farms in 2014 is presented in Supplementary Table 7. In 2014, half of the focal fields had additional spring-sown oilseed rape (1–13 ha) within a 2 km radius. The clothianidin-treated seeds for the focal fields were coated with Elado, a trademarked blend of two active ingredients: clothianidin (400 g l−1) and β-cyfluthrin (+80 g l−1), chosen for this study because it was the predominant seed insecticide treatment in oilseed rape in Sweden and in other parts of Europe63. Clothianidin is taken up by the plant and systemically distributed to all its parts for protection against insects64. β-Cyfluthrin is not considered to be systemic and no residues were detected in samples collected in this study34. Both clothianidin-treated and control seeds were coated with the fungicide thiram in 2013.

The participating farmers were instructed not to use other neonicotinoids in the fields during the study, although other insecticide foliar sprays; primarily Plenum (pymetrozine), Avaunt/Steward (indoxacarb) and Mavrik (tau-fluvalinate; also used as varroacide in beekeeping) were used for pest control (Supplementary Table 8). However, Biscaya, tradename for a spray formulation containing the neonicotinoid thiacloprid, was applied to one control field in 2013, followed by a Mavrik spray 1 week later, and to one treated field in 2014, on both occasions at 0.3 L ha−1. Thiacloprid has a considerably lower acute toxicity for bees than clothianidin65 and only trace amounts of thiacloprid were detected in the pollen, nectar and bee samples in 2013 and none in 2014. While Rundlöf et al.34 did not observe any change in results when fields where Biscaya was applied were excluded from the analyses, we detected some qualitative changes (Supplementary Table 9). These changes could be due to the higher thiacloprid residues detected in 2014, but may just as well not relate to Biscaya but rather to the difference between the Biscaya/Mavrik and the alternative insecticide spray combinations used34 or be due to reduced statistical power.

Honeybee colonies

One hundred and sixteen honeybee colonies were prepared at the end of May 2013 by a professional beekeeper in single full-size Langstroth hives containing two combs with mainly sealed brood (with bees), two full honeycombs (with bees), one drawn out empty comb, five combs with wax foundation, bees shaken from two combs and either a 1-year-old (84 experimental colonies) or 2-year-old (12 experimental colonies plus 20 reserve colonies) queen of known descent to produce relatively small, equally sized (3418 ± 123 adult bees; mean ± s.e.m.; n = 96 colonies) colonies with plenty of room for growth that could become strong enough to survive the coming winter, but not outgrow their space during the summer. Six experimental colonies were placed along the field edge in each of the 16 oilseed rape fields (96 colonies in total) between 14 and 28 June 2013 at the onset of oilseed rape flowering (Supplementary Fig. 1). Queen lineage and age was matched between farm-pairs, but colonies were otherwise distributed randomly. Colonies were kept at a 60 ha organically managed winter oilseed rape field in full bloom before placement at the 16 experimental fields (Supplementary Fig. 1) to ensure pre-experiment colony growth was based as much as possible on pesticide-free foraging.

When the oilseed rape bloom in the experimental field had ceased, the colonies were moved between 2 and 31 July (Supplementary Fig. 1) to a common apiary to overwinter. On 10 August, the colonies were given a formic acid vapour treatment against Varroa mites, consisting of 20 ml 60% formic acid soaked into a flat household sponge placed under the inner cover on top of the frames. The colonies were fed a total of 20 kg of sugar per colony in the form of a 55–60% v/v sucrose solution, provided in a feeding box across three occasions during August–September 2013. An additional light Varroa treatment was carried out on 4 December by sprinkling 30 ml 2.6% oxalic acid in 60% sucrose in between the frames, directly onto the bee cluster. In spring 2014, colonies were moved to an organically managed oilseed rape field before placement at the 10 spring-sown oilseed rape fields (Supplementary Fig. 1). Colonies were considered for inclusion in the 2014 part of the study if they had a 2-year-old, egg-laying queen in April 2014 (excluding colonies that died, re-queened or had 3-year-old queens in April 2014) and had not swarmed by the beginning of June 2014. These restrictions, in addition to the requirement that colonies would be exposed/unexposed to clothianidin for both years of the experiment, meant that in 2014 only four colonies could be allocated to each field. Colonies placed by treated fields in 2013 were again placed by treated fields in 2014, so as to assess the cumulative effects of multi-year clothianidin exposure, but were otherwise re-randomized prior to placement to minimize unintended biases. Even so, two control colonies from 2013 had to be placed by a clothianidin-treated field in 2014, due to insufficient qualifying exposed colonies for the six clothianidin-treated fields. Enough colonies were available for the four control fields. The strength of colonies was equalized as described for 2013, but only within each treatment group. The colonies were reduced and equalized a second time (8 June 2014), after some of them grew too large and attempted to swarm (Supplementary Fig. 1). Each reduced colony included 1 full honey comb (with bees), 3 combs with mainly sealed brood (with bees), and the original queen from 2013 and 6 combs with wax foundation. Colonies were moved to the spring oilseed rape fields between 16 and 25 June 2014 and brought back to the common overwintering site between 14 and 22 July 2014 (Supplementary Fig. 1).

Residue analyses

To confirm clothianidin exposure, 24 adult honeybees per field caught at the hive entrance, pollen pellets collected from 5 honeybees foraging in the oilseed rape fields and nectar removed from the stomachs of 5 nectar-foraging honeybees in the oilseed rape fields were analysed from each site for clothianidin residues. Pollen (>25 ml) was collected using pollen traps, which were installed for 1 day on three colonies per site and analysed for plant species origin. Samples were handled and analysed as in Rundlöf et al.34 and collected during the peak bloom assessments in both years (Supplementary Fig. 1), with the concentrations of clothianidin and four other neonicotinoids used in Sweden (Supplementary Table 2) quantified using liquid chromatography coupled with tandem mass spectrometry (LC-MS/MS) and pollen identified to oilseed rape type using light microscopy and a pollen reference library (see Supplementary Table 2 for limits of detection and quantification). For further analyses of the variation in neonicotinoid exposure of honeybees in different sites, we collected 12 honeybees per site from the hive entrance. This sampling was done at three clothianidin-treated sites in 2013. Nectar was extracted from the honey stomach of the collected bees. The concentrations of clothianidin was thereafter quantified in both nectar and bee tissue for each bee individual. More details on the sample treatment for different matrices, LC-MS/MS method and quality controls are given in Supplementary Methods.

Colony development, re-queening and honey production

Honeybee colony development was assessed by the same trained observer and one assistant. The presence of a laying queen was established, as well as the presence of queen cells. If a re-queening event was accompanied by a large loss of adult bees, it was deemed to have swarmed. If no loss of adult bees was observed, the colony was deemed to have re-queened through supersedure. Colony honey production and development was determined by weighing the colonies and by assessing colony strength using the Liebefeld method66, as the total number of adult bees and the area of capped brood over all frames. The number of adult bees was estimated by counting honeybees on both sides of the 10 frames. The number of capped brood cells (amount of brood) was determined by multiplying the proportion of closed brood coverage by 2700, which is the number of cells on one side of the frames used. The colonies were weighed during pre-exposure and post-exposure assessments (using a Mettler Toledo bench scale able to weigh up to 32 kg with 1 g precision), to estimate honey production. Full honey frames were replaced by empty frames during the oilseed rape bloom, to allow the colony to grow and reduce swarming. Both the full and the empty frames were weighed, for inclusion in the calculating of honey production. During post-exposure assessment, as many honey frames as possible were removed (max 10% of the area covered with covered brood) to simulate the beekeepers honey harvest. Pre-exposure assessments were done at the organically managed winter-sown oilseed rape field on 6–17 June 2013 and 9–11 June 2014, and post-exposure assessments at the common overwintering apiary on 29 July to 9 August 2013 and 28–31 July 2014 (Supplementary Fig. 1). Furthermore, a spring colony strength assessment was performed in April 2014 by estimating the total number of adult bees and the number of capped brood (Supplementary Fig. 1). The colony assessor and assistants were blinded during data collection with respect to the treatment regimen of the fields.

Pathogen and parasite sample collection and processing

Samples of around 100 adult honeybees were taken from each colony during pre- and post-exposure assessments in the clothianidin-treated and control experimental oilseed rape fields in both 2013 and 2014 (Supplementary Fig. 1). Bees were taken from the outer comb of each colony and consisted therefore of a mixture of house bees and forager bees67. All bee samples were stored at −20 °C until the laboratory work was performed. The V. destructor infestation rates for each colony were determined by washing the adult bee samples with soapy water to dislodge and count the mites68. The abdomens of 60 adult honeybees per colony (for individual colony analyses) or per apiary (for the 2013 pooled colony analyses, 10 bees per colony) were removed and placed in a polyethylene bag with an inner mesh (BioReba). The abdomens were ground in the bag using a pestle and 30 ml of nuclease-free (Milli-Q) water (0.5 ml per bee) was mixed thoroughly with the sample to create a homogenous suspension. Several 1 ml aliquots of this suspension were removed and frozen immediately at −80 °C for DNA and RNA extraction and as future reference material.

Parasites, pathogens, symbiotic microbes and immune genes

The collected bee samples were assessed for a variety of pathogenic and non-pathogenic parasites and microbes in order to study the impact placement of colonies at clothianidin-treated fields on their prevalence and abundance. The organisms included the ubiquitous ectoparasite Varroa destructor, 13 viruses: acute bee paralysis virus (ABPV), aphid lethal paralysis virus (ALPV), Big Sioux River virus (BSRV), black queen cell virus (BQCV), chronic bee paralysis virus (CBPV), deformed wing virus type-A (DWV-A), deformed wing virus type-B (DWV-B), Israeli acute paralysis virus (IAPV), Kashmir bee virus (KBV), Lake Sinai virus strain 1 (LSV-1) and strain 2 (LSV-2), Sacbrood virus (SBV), slow bee paralysis virus (SBPV); two common honeybee microsporidian gut parasites (Nosema apis and Nosema ceranae) and two symbiotic gut bacteria (Gammaproteobacterium: Gilliamella apicola and Betaproteobacterium: Snodgrassella alvi). For the 2013 samples we also analysed at apiary level the mRNA levels of eight honeybee genes (Amel/LRR, Apidaecin, cSP33, Dorsal-1A, Eater-like, NimC2, PGRP-S2 and SPH51) whose expression had previously been linked to pesticide, pathogen and/or parasite exposure19,44 and (social) immunity in honeybees45.

Nucleic acid extraction

DNA was extracted from the bee homogenates using the protocol for extracting DNA from Nosema spores69, which is sufficiently robust to also extract DNA from bacteria and other microorganisms. A total of 500 µl primary bee homogenate was centrifuged for 5 min in a microfuge at 13,000 rpm. The pellet was repeatedly frozen-thawed with liquid nitrogen and ground with a sterile Teflon micropestle until pulverized. The pulverized pellet was resuspended in 400 µl Qiagen Plant tissues DNeasy AP1 lysis buffer containing 4 µl RNAse-A (10 mg ml−1) and incubated and shaken for 10 min at 65 oC, after which 130 µl P3 neutralization buffer (3.0 M potassium acetate pH 5.5) was added, followed by 5 min of incubation on ice and centrifugation for 5 min at 14,000 rpm to remove the lysis debris. DNA was purified from 500 µl of the supernatant by the Qiagen automated Qiacube extraction robot, following the plant DNeasy protocol and eluting the DNA into 100 µl nuclease-free water. RNA was extracted by the Qiacube robot directly from 100 µl primary honeybee homogenate using the Qiagen Plant RNeasy protocol (including the Qia-shredder for additional homogenization70) and the RNA was eluted into 50 µl nuclease-free water. The approximate nucleic acid concentration was determined by NanoDrop, after which the samples were diluted with nuclease-free water to a uniform 10 ng µl−1 (DNA and LSV-1(RNA)) or 20 ng µl−1 (for all other RNA samples) and stored at −80 °C.

RT-qPCR and qPCR

The various microorganisms and host mRNA targets were detected and quantified by either One-Step reverse transcription-quantitative PCR (RT-qPCR) for pathogens with a RNA genome and the immune and internal reference gene mRNA targets, or by quantitative PCR (qPCR) for organisms with a DNA genome. Details of the assays are shown in Supplementary Table 10, Supplementary Table 11 and Supplementary Table 12. The reverse primer for Amel/LRR was slightly re-designed from Di Prisco et al.19 because the extremely high complementarity between the original forward and reverse primers resulted in high levels of PCR artefacts dominating the quantitative signal. The reactions were conducted in duplicate, in 20 µl (DNA) or 10 µl (RNA) reaction volumes containing 2 µl (DNA) or 1.5 µl (RNA) template, 0.4 µM (DNA) or 0.2 µM (RNA) of forward and reverse primer and either the Bio-Rad Eva Green qPCR mix (DNA) or the Bio-Rad One-Step iTaq RT-qPCR mix with SYBR Green detection chemistry (RNA). The reactions were incubated in 96-well optical qPCR plates in the Bio-Rad CFX connect thermocycler, using the following amplification cycling profiles: 10 min at 50 °C for complementary DNA (cDNA) synthesis (RT-qPCR only): 5 min at 95 °C (to inactivate the reverse transcriptase and activate the Taq polymerase) followed by 40 cycles of 10 s at 95 °C for denaturation and 30 s at 58 °C for primer annealing, extension and data collection. For DNA assays the following amplification cycle profiles were used: 2 min at 98 °C for the initial denaturation followed by 40 cycles of 5 s at 98 °C for denaturation and 10 s at 60 °C for primer annealing, extension and data collection. The amplification cycles were followed by a melting curve analysis to determine the specificity of the amplification by reading the fluorescence at 0.5 °C increments from 65 °C to 95 °C. Included on each reaction plate were positive and negative (non-template) assay controls. For each type of assay (Supplementary Table 10, Supplementary Table 11 and Supplementary Table 12) a calibration curve was prepared through a 10-fold dilution series of a positive control of known concentration covering 6 orders of magnitude, for quantitative data conversion, establishing the reference melting curve profile of the amplicon and estimating the reaction performance statistics.

Data conversion and normalization

The melting curves of individual reactions were evaluated visually in order to separate out non-specific amplifications, which differ in melting temperature profiles from true target cDNA/DNA amplicons. Non-specific amplifications were deleted from the dataset. All assays were run in duplicate, with the mean value of these two duplicates used in further calculations. Both duplicates had to yield a positive quantitative value and pass the melting curve analysis for the data to be included in the dataset. The raw RT-qPCR data of all confirmed amplifications were subsequently converted to estimated copy numbers of each target RNA, using the corresponding calibration curve for the assay. These data were multiplied by the various dilution factors throughout the procedure to calculate the estimated copies of each target per bee69. Since RNA is easily degraded, there is a risk that differences between individual samples in RNA quality (i.e., degradation) can affect the results70. To correct for this, RT-qPCR assays for the mRNA of a common honeybee internal reference gene (RP49) were run on all samples. The data for the RNA targets of interest were then normalized to the average value for RP49 mRNA, thus correcting the data for sample-specific differences in RNA quality with respect to RT-qPCR performance70.

Statistical analyses

The proportions of honeybee-collected pollen that originated from oilseed rape-type plants were compared between treatments (clothianidin seed treatment/untreated) using a generalized linear model assuming binominal distribution and correcting for overdispersion. The clothianidin concentrations in nectar and pollen collected by honeybees and in bee tissue were compared between treatments using Wilcoxon–Mann–Whitney tests. To compare the concentrations of clothianidin in the bee tissue and nectar in individual bees between fields, we used analyses of variance (ANOVA), with field identity as predictor. Furthermore, the concentrations of clothianidin in the tissue and nectar stomach content of individual bees were related using a multiple linear regression with field identity and concentration of clothianidin in nectar as explanatory variables and concentration of clothianidin in bee tissue as response variable.

The study followed in general a Before-After-Control-Impact (BACI) design, with a paired field structure, repeated for two consecutive years at colony level for data on colony development as well as the prevalence and abundance of parasites, pathogens and gut bacteria. The years, 2013 and 2014, were analysed together in one full model, with seed treatment, bloom, year and their interactions as fixed factors. The effect of the clothianidin treatment was assessed by the interaction between bloom and seed treatment, as this term reflects the difference in change between treatments over the oilseed rape bloom(s). If the three-way interaction (bloom × seed treatment × year) was significant (i.e., if the variable responded differently to the clothianidin treatment from one year to the next), the dataset was split by year and year was dropped as a fixed factor. Furthermore, the dataset was analysed for only 1 year if the data consisted of a sample size ≤10 in one year both for microbiota prevalence and abundance (Supplementary Table 1). Colonies that swarmed (8 at control fields and 10 at clothianidin-treated fields in 2013; one at a control field and two at clothianidin-treated fields in 2014) were excluded from the analysis, since swarming has a large effect on colony development. Also excluded is the single colony that lost its queen during transport before field placement in 2013 (treated field). Excluding colonies that swarmed from the analysis qualitatively altered some results (see Supplementary Table 13). Changes in significance level might be due to reduced statistical power, random chance or biological effects.

Linear mixed-effects models (LMM) were used to test the effect of the clothianidin treatment on colony development measured as the number of capped brood cells (amount of brood) and the number of adult bees. Seed treatment (clothianidin or control), bloom (before or after oilseed rape bloom), year (2013 or 2014) and their interactions were fixed factors. Farm pair identity, farm identity and colony identity were included as random factors. Honey production was compared between treatments using a LMM with farm pair identity and farm identity as random factors. Generalized linear mixed models (GLMMs) were used to test the influence of the clothianidin treatment on re-queening and mortality of the colonies with farm identity as a random factor.

The influence of the clothianidin treatment on spring colony development, measured as the number of adult bees and amount of brood, was tested using a LMM and a GLMM, respectively, with seed treatment as fixed factor and farm pair identity and farm identity as random factors. For the number of capped brood cells we used a negative binomial error distribution and logarithmic link function.

The microbiome and Varroa mite data were analysed both on their binomial (presence/absence) and quantitative (abundance) character, using GLMMs (with binomial error distributions and a logit link function) and LMMs (with normal error distributions), respectively, with seed treatment, bloom, year and their interactions as fixed factors. GLMMs on microorganism or Varroa mite prevalence included only colony identity as random factor, as the effective sample size (i.e., the less frequent outcome of the presence/absence data) did not allow for the inclusion of more random factors. Only organisms and years with an (effective) sample size >10 were analysed for both the prevalence and the abundance data. In addition, colonies that did not at least once test positive for a particular microorganism were excluded from the analysis of abundance. Bee pathogen and bacterial abundance were logarithmically (log 10 ) transformed, as they are generally exponentially distributed. LMMs on target organism abundance contained farm pair identity, farm identity and colony identity as random factors. Varroa mite numbers per 100 bees and colony weight were square-root transformed to avoid non-normally distributed residuals. Confidence intervals were calculated based on profile likelihood. For the square-root transformed data, estimates were back-transformed to the original scale for graphical illustrations.

The immune gene transcripts were only available for 2013 and at apiary level but also followed the BACI design. LMMs on gene expression contained seed treatment and bloom as fixed factors and farm identity as random factors.

Statistical data analyses were performed using R except for analyses addressing the verification of clothianidin exposure and land use, for which SAS 9.4 for Windows (SAS Institute Inc.) was used. LMMs were fit using the lmer function of the lme4 package and GLMMs were fit using the glmmTMB function of the glmmTMB package in R. P values from GLMMs were calculated by likelihood ratio tests. P values from LMMs were calculated using the Anova function of the car package, whereby type-III F-tests were used for models containing interactions and type-II F-tests for models without interactions (i.e., spring assessment and honey production). Effects of fixed factors were estimated using sum-to-zero contrasts in all models except those on neonicotinoid residues. Sum-to-zero contrasts allow for the determination of main effects/interactions (i.e., estimation independently of other independent variables) and show the effects of factors as deviations from the grand mean (intercept). For factors with two levels the magnitude of the deviation of each level from the grand mean is the same but the direction differs. We represent effects of factors (seed treatment, bloom, year), as the deviance of the second level (clothianidin, after, 2014) from the grand mean. This was also the case for interactions, so that for example the seed treatment × bloom interaction indicates to what extent clothianidin-exposed colonies differed in change over the oilseed rape bloom from the mean change of both treatments.

Power analysis

We performed power analyses for number of adult bees and the number of capped brood cells and honey production to investigate the effect size we could potentially detect given our design, replication and model choice. Power was determined for a range of effect sizes at a nominal confidence level of α = 0.05 by 1000 Monte Carlo simulations per effect size using the powerSim function of the simr package. Power was calculated for a range of effect sizes, expressed as the change in number of adult bees, number of capped brood cells or honey production. By dividing the effect size with the mean number of bees, the mean number of brood cells or honey production of all control colonies, we obtained effect size expressed as the percentage change of those matrices (Supplementary Fig. 3). This power analysis made it possible to compare our effect size with the effect size presented by Rundlöf et al.34. Using the full model, we could detect an effect size for the number of adult bees of below 5% with a power of 80% compared to the effect size below 20% presented by Rundlöf et al.34. This is even lower than the requirements of an effect size <7% set by EFSA32. As a result of a significant interaction of seed treatment, bloom and year, the dataset for the number of capped brood cells were analysed separately for each year. Therefore, we present here the power analysis for each year. The effect size at which 80% was reached increased from below 10% in 2013 to slightly below 11% in 2014 (Supplementary Fig. 3), likely due to the reduced replication in 2014. We also performed a power analysis of honey production (amount of honey per colony in kg) using the dataset of both years, showing that an effect size of below 20% could be detected with a power of 80% (Supplementary Fig. 3).

Reporting summary

Further information on experimental design is available in the Nature Research Reporting Summary linked to this article.