Honeybee colonies compensate for pesticide-induced effects on royal jelly composition and brood survival with increased brood production

Sublethal doses of pesticides affect individual honeybees, but colony-level effects are less well understood and it is unclear how the two levels integrate. We studied the effect of the neonicotinoid pesticide clothianidin at field realistic concentrations on small colonies. We found that exposure to clothianidin affected worker jelly production of individual workers and created a strong dose-dependent increase in mortality of individual larvae, but strikingly the population size of capped brood remained stable. Thus, hives exhibited short-term resilience. Using a demographic matrix model, we found that the basis of resilience in dosed colonies was a substantive increase in brood initiation rate to compensate for increased brood mortality. However, computer simulation of full size colonies revealed that the increase in brood initiation led to severe reductions in colony reproduction (swarming) and long-term survival. This experiment reveals social regulatory mechanisms on colony-level that enable honeybees to partly compensate for effects on individual level.

3 SFig 1. Sampling scheme. Each week (W1-7), all combs were photographed and samples of worker bees, larvae, and worker jelly were taken. To determine the size of the hypopharyngeal gland, 15 newly-hatched and marked worker bees were introduced into each colony during week 1 and recovered in week 2 (age-defined marked bees, green arrows). In week 7, randomly chosen worker bees of undefined age were taken from brood combs for the same reason (random bees).
The total number of eggs, larvae and pupae were recorded for every week and colony (brood quantification). To determine the survival rate of individually tracked larvae, at least 50 individual cells per colony were followed over 4 weeks in the first half of the experiment (first brood cycle, brood survival phase A) and in the second half of the experiment (second brood cycle, brood survival phase B).

SFig 4. Protein composition of worker jelly.
For protein analysis, samples of three colonies of each treatment group were analyzed from week 4.
To test, whether the concentration of a protein subfraction could be affected in clothianidin exposed bees, 15 µg of protein of each sample were analyzed by SDS-PAGE in 12% polyacrylamide gels followed by staining with Coomassie Brilliant Blue. The concentration and composition of major proteins in worker jelly was unaffected by clothianidin exposure. Although, there was no statistically significant difference in the protein concentration (see above), we observed a decrease in the concentration in major worker jelly proteins for the 100 µg/L treatment group. This could be related to a possible increase of protein degradation in this treatment group. This would lead to an increase of peptides, which appear in the total protein measurement, but are too small to be visualized on this gel. Alternatively, the general decrease of protein levels in the SDS PAGE in the 100 µg/L group could also be due to high molecular weight aggregates that are unable to enter the gel.

Clothianidin uptake
Suppl. We analyzed the sugar syrup used to feed the experimental colonies to assess its exact clothianidin content. For all four experimental groups, the final spiked pesticide levels were close to the target concentrations (Suppl. Table 1). The consumption of the spiked syrup was recorded every week to estimate the total clothianidin uptake over the study period. The provided amount of 400 mL (= 540 g) syrup each week per colony was completely consumed each week with the exception of week 6 in the highest concentration treatment. All five colonies in the 100 µg/L clothianidin treatment group were visibly weakened and did not consume all the sugar syrup, which was provided during week 6 (residual syrup: 76. 7, 114.5, 142.2 and 18.0 g, respectively). To determine the clothianidin levels in bees exposed to the different pesticide levels, 10 randomly chosen worker bees from each hive were analyzed on week 7 (Suppl. Table 1). Clothianidin was detected in the bee sample from only one colony of the control group (colony I, 0.004 ng/bee).

Analytical Method
LC/MS/MS was used for the identification and quantification of the substances in the samples. The system used was a Prominence UFLC XR HPLC (SHIMADZU) coupled to a triple quadrupole mass spectrometer 4000 Q TRAP (AB SCIEX) equipped with an electro spray ionization (ESI) source.
Clothianidin and its metabolites clothianidin-metabolite TZMU and clothianidin-metabolite TZNG were identified by their retention time and three MRM transitions. The residues in the samples were quantified with reference standards in matrix (concentrations: 0. 1, 0.5, 1, 5, 10, 25, 50 and 100 pg/µL). The quantification was carried out by the internal standard method. The values shown for Self-made helper functions are "sourced": > source("HelperFcts.R") 1 1 Raw data The original raw data had been saved in several MS-Excel-sheets and exported as 'comma-separated values' (CSV) files, whose fields have been separated by semicolon (;). The decimal sign is the dot (.). The columns in the CSV-file possess names (in their first rows).

Hyperpharingal gland sizes
From the imported files we extract the data on honeybee hyperpharingal gland sizes and create two data frames named "W2HPGland" and "W7HPGland". We check data consistency by controlling the first few rows of the created data frames as well as their summaries.

Larval Survival
Here, we create the two data frames Survival1 and Survival2, and again check data consistency by controlling the first few rows of the created data frames as well as their summaries.
In the present case, we unfortunately have neither large sample sizes nor a balanced design, so we have to resort to an alternative of the approximate χ 2 -tests. In fact, we shall use two (to double-check the results) alternative, reliable methods for analysing fixed effects: a) Kenward and Roger (KR) provide a modification of a Wald test statistic which has under the null hypothesis asymptotically an F -distribution (whose denominator degrees of freedom need and can be estimated) and is said to yield results more reliable than the F -test mentioned above. (This KR approach is for models fitted with restricted maximum likelihood (REML).) b) A parametric bootstrap (PB) allows to determine the distribution (or moments thereof) of the likelihood ratio test statistic under the null hypothesis. (This is for models fitted with maximum likelihood (ML); so models fitted with REML need to be re-fitted with ML before.) The R-package pbkrtest implements both the KR and the PB approach for tests regarding the fixed effects (in linear mixed models with independent errors). See [5] which describes the methods, their implementation and which contains also examples.
Model diagnostics for the model in fit1 are found on page 11.

Posthoc tests in week 2 -comparisons with the Control
We perform Dunnett's multiple comparisons for the one-sided null hypothesis that the Clothianidin treatments do not yield smaller HPG sizes than seen in the Control group.
Model diagnostics for the model in fit1:  Week 7

HiveID
Gland size in µm Control 1 µg/l 10 µg/l 100 µg/l *** Figure 5: Per week and by treatment: Fixed-effects estimates (solid blue circles) with (non-simultaneous!) two-sided 95 %-confidence intervals (blue) of HPG size. (Week 2 on the left, week 7 on the right.) Black error bars indicate the mean HPG size plus/minus one standard error of the mean (ignoring the grouping structure by hives within each treatment). Note that the fixed-effects estimates differ slightly from the treatment specific "simple" mean HPG sizes. They are closer to the overal mean (why they are sometimes also called "shrinkage estimators"). Asteriscs indicate significant differences in one-sided Dunnett contrasts with "Control" of the mixed-effects ANOVA in §2.

Treatment effect on number of eggs
We analyse the number of eggs by a linear mixed-effects model with a fixed effects linear or parabolic time trend along weeks, with fixed treatment main effects, and fixed interaction effects between weeks and treatment as well as between squared weeks and treatment. This is done for weeks centered at 3 -arbitrarily selected -so that the main effect of treatment represents the estimated number of eggs at week 3. Hives are modelled as random shift effects, thus accounting for the within-hive correlation.
For explanations regarding the testing methodology and for comments regarding the R-code see §2.1.2. We again follow, but also extend §10.6 in [6, Faraway (2016)].
Summary: Allowing (only) a linear time trend along weeks, including an interaction between time and treatment, there is no significant influence of treatment on the time trend of number of eggs (p-value = 0.2633 from the Kenward-Rogers method, and p-value = 0.2813 based on parametric bootstrap).
However, using a parabolic time trend along weeks, including an interaction between squared time and treatment (to account for the visible and diffferent curvatures of the time trends of numbers of eggs), there is a significant influence of treatment on the number of eggs (p-value = 2.587e-06 from the Kenward-Rogers method, and a p-value = 2e-04 based on parametric bootstrap). This means in particular, that the time trends of the number of eggs are not the same in the four treatment groups.
Note: It is actually of little interest to analyse the main effects of treatment since they characterize the situation at only a single point in time, namely here at week 3. Instead, we focus on comparing the treatments with the control with respect to the (local) slope of their trend at week 3 and with respect to their (global) curvature. The first is represented by the interaction effect of treatment and week (centered at 3), and the latter by the by the interaction effect of treatment and squared week.
Simultaneous Tests for General Linear Hypotheses Summary: All Clothianidin treatments are significantly different from the control with respect to both their (local) slope at week 3 and their (global) curvature, with the exception of the 100 µg/l group: it does not differ significantly from the control in respect of its curvature.

Treatment effect on number of larvae
In proceeding analogously to §2.2.2 it turned out that the random effects of hives are in both the linear mixed-effects model with a linear time trend and the one with a parabolic time trend not significant. We tested this using exact restricted likelihood ratio tests implemented in package RLRsim. Analysis and results are not shown.
Hence, we decided to analyse the number of larvae by (purely) fixed-effects models with a linear or parabolic time trend along weeks, with treatment main effects, and interaction effects between weeks and treatment as well as between squared weeks and treatment. As in the previous paragraph this is done for weeks centered at 3 so that the main effect of treatment represents the estimated number of larvae at week 3.    Summary: Only the 100 µg/l group is significantly different from the control with respect to both their (local) slope at week 3 and their (global) curvature.

Treatment effect on number of capped brood cells
In the analysis of capped brood cells we experienced the same phenomenon as with the number of larvae in §2.2.4: the random effects of hives are in both the linear mixed-effects model with a linear time trend and the one with a parabolic time trend not significant (using again exact restricted likelihood ratio tests implemented in package RLRsim; neither analysis nor results shown).
So, the number of caped brood cells is also analysed by (purely) fixed-effects models with a linear or parabolic time trend along weeks, with treatment main effects, and interaction effects between weeks and treatment as well as between squared weeks and treatment. As in the previous paragraph this is done for weeks centered at 3 so that the main effect of treatment represents the estimated number of larvae at week 3.    Exploratory data analysis using q-q plots revealed that the raw data (displayed in fig. 14) are neither normally distributed nor homoscedastic. The variance of the data seems to increase with decreasing proportions of surviving larvae. To compensate that we take the logarithm of the proportions of dead larvae (= 1 -31 2.3 Question 3: Impact on the larval survival proportion of surviving larvae) to obtain (approximately) normally distributed and homoscedastic values.

Survival−Phase 1
Week Log(Proportion of dead larvae)  fig. 14), we exclude the last week of phase 1 from the following analysis. This is similar in phase 2 where proportions remain constant in the last two weeks in all but only two hives (one in the Control group and one in the 1 µg/l treatment; see bottom part of fig. 14). Since in addition, phase 2 has 100 % dead larvae from week five on in three of the four hives in the 100 µg/l treatment (see bottom part of fig. 14), we exclude this treatment group in phase 2 completely from the following analysis.

Treatment effect on larval survival in phase 1
We analyse the log-proportions of dead larvae using a two-factorial mixed-effects ANOVA with treatment and week as fixed-effects factors including interactions and with hive as random-effects grouping factor. (This model is the same as a two-factorial repeated measures ANOVA with treatment and week as fixed-effects factors and hive as random-effects factor for grouping the repeated measurements.) We restrict the following analysis to weeks 2 and 3 of phase 1, and we proceed as in §2.1.2; for explanations and comments see also there.
Since there is a significant interaction effect between treatment and week (p-value = 0.02929 from Kenward-Rogers method, and p-value = 0.02859 based on parametric bootstrap), the main effect of each factor has to be interpreted as averaged across the levels of the other factor: The averages (across the four treatments) of log-proportions of dead larvae are significantly different between the two weeks, but the averaged logproportions of dead larvae (across the two weeks) are not significantly different between the four treatment groups. A bit more concrete: Interpreting the week main effect averaged across the levels of treatment implies that we compare the two points in time without separating the treatments. Even more precise: in each week we consider the population mean (of log-proportions of dead larvae) after "pooling" over the four treatment groups, i.e., the data from the four treatments are combined in each week as if there were only a single treatment group. In turn, interpreting the treatment main effect averaged across the levels of week implies that we compare the four treatments without separating the weeks. Even more precise: in each treatment group we consider the population mean (of log-proportions of dead larvae) after "pooling" week 2 and week 3, i.e., the data from the two weeks are combined in each treatment group as if there were only a single point in time.

Simultaneous Tests for General Linear Hypotheses
Fit: lmer(formula = log ( Summary: If we combine the two test families from above into one family all statements from above hold still true on a significance level slightly higher than 5 %: a level of 6 % would yield the same conclusions.

Treatment effect on larval survival in phase 2
Recall the note on p. 32: Proportions of dead larvae remain constant in the last two weeks in almost all hives and the 100 µg/l treatment group has 100 % dead larvae from week five on in three of the four hives (see bottom part of fig. 14). Therefore, we exclude the 100 µg/l treatment group completely from the following analysis and restrict the analysis to weeks 5 and 6 of phase 2. We proceed as in §2.3.2.
Since there is a significant interaction effect between treatment and week (p-value = 0.03776 from Kenward-Rogers method, and p-value = 0.03559 based on parametric bootstrap), the main effect of each factor has to be interpreted as averaged across the levels of the other factor: The averages (across the four treatments) of log-proportions of dead larvae are significantly different between the two weeks, and the averaged logproportions of dead larvae (across the two weeks) are significantly different between the four treatment groups. For help in interpreting the above see the "more concrete" explanations in the summary on page 35.
Recall: the significant interaction effect between treatment and week means that changes (!) in (average) log-proportions of dead larvae along weeks are different between treatments. We proceed completely analogous to §2.3.3.

Question 4: Compensation rate
The main idea is to calculate and show the compensation effort for each treatment. The bees have to compensate for the higher larval mortality and therefore expend more energy.

Larvae-to-eggs ratio: Longitudinal EDA
To estimate/assess the development of the hives along time we calculate larvae-to-eggs ratios within each hive for each week. For this ratio to be sensible, the number of larvae in week w has to be compared with the number of eggs of the previous week w − 1 (as the larvae of week w have been "produced" by the eggs of week w − 1). The time courses of these "shifted" larvae-to-eggs ratios are then compared among treatments.
One point to consider is, that of all eggs seen in week w − 1 approximately two thirds have already turned into capped brood at inspection in the subsequent week w because they were already one or two days old at inspection in week w − 1. In turn, the actually observed number of larvae in week w can approximately be only one third of the number of eggs seen in the previous week. However, we will not take this into consideration since it would complicate things unnecessarily since the deviation is very likely just a constant factor. We therefore take the raw weekly numbers of eggs and larvae as good estimates.
Week Larvae(w)/Eggs(w−1)   The following two figures display just the mean and median larvae-to-eggs ratio profiles, respectively, i.e., without the raw values.

'Shifted' larvae−to−eggs ratios
Week Mean of larvae(w)/eggs(w−1) We again proceed as in §2.1.2, i.e., we follow, but also extend §10.6 in [6, Faraway (2016)]. For explanations regarding the testing methodology and for comments regarding the R-code see §2.1.2.
Note: We exclude the extreme outlier above 15 in treatment 10 µg/l from the following analysis! > EL <-subset(EL0, subset = L2E < 15) We analyse the larvae-to-eggs ratio by a linear mixed-effects model with a fixed effects linear time trend along weeks, with fixed treatment main effects, and fixed interaction effects between weeks and treatment. This is done for weeks centered at 2 -arbitrarily selected -so that the main effect of treatment represents the estimated larvae-to-eggs ratio at week 2. Hives are modelled as random shift effects, thus accounting for the within-hive correlation.  Summary: Allowing a linear time trend along weeks with interaction between time and treatment there is a significant interaction between time and treatment (p-value = 0.007611 from the Kenward-Rogers method, and p-value = 0.008998 based on parametric bootstrap). This means in particular, that the time trends of the larvae-to-eggs ratio are significantly different between the four treatment groups.
However, there is no significant treatment main effect in the model with interaction (p-value = 0.07313 from the Kenward-Rogers method, and p-value = 0.07479 based on parametric bootstrap). This means, the estimated average larvae-to-eggs ratio at week 2 (!) are not significantly different between the four treatment groups.    (201712) # For reproducibility. > # Relevant: p-value in row "PBtest" > # (PB1 <-summary(PBmodcomp(largeM = linfit2, smallM = linfit1, nsim = nsim.PBmodcomp, > # cl = cl))) > (PB2 <-summary(PBmodcomp(largeM = linfit3, smallM = linfit2, nsim = nsim.PBmodcomp, + cl = cl))) Summary: Allowing a linear time trend along weeks with interaction between time and treatment there IS NO significant interaction between time and treatment (p-value = 0.06655 from the Kenward-Rogers method, and p-value = 0.07199 based on parametric bootstrap). This means in particular, that the time trends of the larvae-to-eggs ratio ARE NOT significantly different between the four treatment groups.
However, there IS A significant treatment main effect in the model with interaction (p-value = 0.02906 from the Kenward-Rogers method, and p-value = 0.02619 based on parametric bootstrap). This means, the estimated average larvae-to-eggs ratio at week 2 (!) ARE significantly different between the four treatment groups.
Week Brood(w)/Larvae(w−1)   The following two figures display just the mean and median capped-brood-to-larvae ratio profiles (shifted by lag 1), respectively, i.e., without the raw values.

'Shifted' capped−brood−to−larvae ratios
Week Mean of brood(w)/larvae(w−1) Control 1 µg/l 10 µg/l 100 µg/l Figure 25: Medians (across hives) of "shifted" capped-brood-to-larvae ratios along weeks by treatment. As eggs of week w − 1 are the larvae of week w the numbers of larvae are shifted by one week. Note that the median values in the control group are quite consistently the largest values and that all "low-dose"-groups (incl. Control) stay above 1.5. (File name: Cloth-Q4 MedianB2L along Weeks by Treat.pdf ) 56 2.4.7 Treatment effect on capped-brood-to-larvae ratio We proceed as before in §2.4.2; for details see there.
And there is no significant treatment main effect in the model with interaction (p-value = 0.6764 from the Kenward-Rogers method, and p-value = 0.7097 based on parametric bootstrap). This means, the estimated average capped-brood-to-larvae ratios at week 2 (!) are not significantly different between the four treatment groups.

HiveID
Standard normal quantiles