Chronic sublethal stress causes bee colony failure

Current bee population declines and colony failures are well documented yet poorly understood and no single factor has been identified as a leading cause. The evidence is equivocal and puzzling: for instance, many pathogens and parasites can be found in both failing and surviving colonies and field pesticide exposure is typically sublethal. Here, we investigate how these results can be due to sublethal stress impairing colony function. We mathematically modelled stress on individual bees which impairs colony function and found how positive density dependence can cause multiple dynamic outcomes: some colonies fail while others thrive. We then exposed bumblebee colonies to sublethal levels of a neonicotinoid pesticide. The dynamics of colony failure, which we observed, were most accurately described by our model. We argue that our model can explain the enigmatic aspects of bee colony failures, highlighting an important role for sublethal stress in colony declines.


INTRODUCTION
The social bees (honeybees, bumblebees and stingless bees) provide an important ecosystem service, pollinating both wildflowers and agricultural crops (Klein et al. 2007;Winfree et al. 2008;Potts et al. 2010). Reports of social bee colony losses (Oldroyd 2007;Ratnieks & Carreck 2010) and global bee population declines (Biesmeijer et al. 2006;Goulson et al. 2008;Brown & Paxton 2009;Cameron et al. 2011) are therefore of major concern. Several stressors have been implicated in bee declines (Vanbergen & The Insect Pollinators Initiative 2013), including pesticides (Desneux et al. 2007), disease and parasites (Brown & Paxton 2009) and habitat change and loss Potts et al. 2010). While there is still debate over which stressors are most detrimental, no single factor has emerged as an overall primary cause (Ratnieks & Carreck 2010;Vanbergen & The Insect Pollinators Initiative 2013).
Recent evidence has indicated that many environmental stressors can affect bees even when they do not cause direct mortality (so called sublethal impact). For instance, exposure to field level pesticides can affect worker mobility, memory, orientation and foraging performance (Desneux et al. 2007;Gill et al. 2012;Schneider et al. 2012;Williamson & Wright 2013), and parasites can impose energetic stress (Moret & Schmid-Hempel 2000;Mayack & Naug 2009) as well as impair both learning (Iqbal & Mueller 2007) and thermoregulation (Schafer et al. 2011). While these sublethal effects need not kill individual bees, they may have profound effects on the dynamics and functioning of the whole colony. Despite this, we currently have little understanding of how chronic sublethal stress, experienced by individual bees, can cause colonies to fail.
There is an important difference between the effects of lethal and sublethal stress when considering the functioning of social bee colonies. Redundancy in bee colonies allows them to lose a significant proportion of their worker force without any apparent significant impact on colony function and productivity (Schmid-Hempel & Heeb 1991;M€ uller & Schmid-Hempel 1992;Rueppell et al. 2009). However, if bees become impaired rather than die, the impairment may impose a load on the colony and lead to a cumulative effect on normal colony function. Indeed, there is evidence that pesticide induced behavioural impairment can detrimentally feedback to colony eclosion (birth) and death rates (Gill et al. 2012), the production of sexuals (Whitehorn et al. 2012) and increase the prevalence of disease (Brown et al. 2000). This feedback onto birth and death rates generates Allee effects which lead to colony dynamics that are uncertain, such as those with multiple outcomes and breakpoints (May 1977). We show here that effects of sublethal stress on colony function are important to explain colony dynamics. We did this by formulating a mathematical model for colony dynamics that includes colony function. We subsequently fitted this model to novel empirical data from bumblebee colonies that were sublethally exposed to a pesticide and show that the description of this model is superior to models which do not account for impairment to colony function.

Modelling
Our model accounts for SubLethal Stress (henceforth, the SLS Model) and describes healthy bees (S ) and impaired bees (I ) to represent behavioural impairment by a sublethal stressor. Healthy bees become impaired at rate b, and the level of behavioural impairment is expressed by the parameter c ≤ 1 which reflects the reduced contribution of the impaired bees to colony function, so that the effective operational size of the colony is N = S + cI. The model describes a growing colony, where the eclosion rate of new bees in the colony is bN, with b as a rate constant. To capture individual effects of stress, known to increase the mortality rate of bees (Cartar & Dill 1991;Brown et al. 2000;Dainat et al. 2012), we introduced mortality of impaired bees at a constant rate m to our model. We assume that the per capita death rate is inversely proportional to the effective size of the colony, l/(N + φ), where decreasing φ sets how sharply the death rate increases at low effective colony sizes and l adjusts the overall rate. This reflects that colonies containing more impaired bees are less good at maintaining essential colony functions (e.g. foraging, thermoregulation, defence and hygienic behaviour). This gives the following equations: To visualise the dynamics of the SLS Model, phase plots were generated by numerically integrating the model over time and recording the trajectories of the variables. Many numerical runs were done, iterating over different initial values of S and I while keeping all other parameters the same. Each run was over a fixed time period (20 days).
The feedback of stressed bees onto colony function in the SLS Model distinguishes it from other mathematical models that model disease and increased mortality in bee colonies (Schmid-Hempel 1998;Sumpter & Martin 2004;Khoury et al. 2011Khoury et al. , 2013Becher et al. 2013). To assess the validity of the SLS Model, we compared it with three other alternative models (their equations are listed in Table 1) summarised as follows: (1) A model described by Khoury et al. (2011), that we therefore dub the Khoury model, describes how colony dwelling bees can transition between hive and forager roles. The model also describes the rate at which hive bees eclose and the rate at which forager bees experience mortality, see Khoury et al. (2011) for a full description. It has been used previously to describe colonies treated with pesticides (Henry et al. 2012).
(2) The SLS Variant Model is similar to the SLS Model, but has a linear per capita death rate which is fixed at l rather than l/(N + φ). It will be used as a comparison to test the importance of the nonlinear death rate term which introduces positive density dependence in the SLS Model.
(3) The final model is the Larvae Adult Model (henceforth the LA Model) which considers that the pesticide could have a toxic effect on larvae, increasing its death rate. It models larva (number L) and adult bees (of number A). Larvae hatch at rate g, die at rate k, eclose into adult bees at rate c and adult bees die at rate ξ.

Colony experiments
We provided bumblebee (Bombus terrestris) colonies a specified amount of sucrose solution (every 2 days and 3 days over the weekend) containing a sublethal concentration of a systemic neonicotinoid pesticide (10 ppb imidacloprid, see Supporting Information for further details). Neonicotinoids are used extensively on horticultural plants and flowering agricultural crops throughout the world. Imidacloprid is applied to crops such as oilseed rape (canola), sunflower, maize, linseed, peas and beans, cucurbits and orchards which are key food sources for bees (Cresswell 1999;Thompson 2001;Holzschuh et al. 2011). The concentration of imidacloprid that we used in this experiment was 10 ppb, which falls near the upper end of the field realistic range reported for nectar and pollen in agricultural crop species (Bonmatin et al. 2003(Bonmatin et al. , 2005Chauzat et al. 2006Chauzat et al. , 2009Krischik et al. 2007;Cresswell 2011;Blacquiere et al. 2012). We also provided colonies with a specified amount of untreated pollen every 2 days (3 days over the weekend). All colonies were provided with increasing levels of pollen and sucrose throughout the experiment to account for increasing demand as the colonies grew (see Supporting Information Section 1.2 for more details on pesticide treatment and feeding regime).
Colonies were housed in wooden nest boxes isolated from one another, and the outside world, in laboratory conditions throughout the experiment. We recorded the growth and development of 16 early-stage queen-right colonies (eight control colonies and eight pesticide-exposed colonies). To eliminate any bias, colonies were sorted by size at the start so that there was no significant difference in the number of workers (mean AE SEM: control = 6.5 AE 0.8; neonicotinoid = 8.25 AE 1.64; t-test: t = À0.96, P = 0.36) or pupae (mean AE SEM: control = 6.4 AE 1.8; neonicotinoid = 6.8 AE 1.53; t-test: t = 0.16, P = 0.88) per colony between control and neonicotinoid assigned colonies (see Supporting Information Section 1.1 for more details on the colony set-up). B. terrestris has relatively small colonies, compared to honeybees, making it possible to track all individual eclosions and deaths. Throughout the 42-day experiment, all colonies were inspected on a day-to-day basis to record the numbers of new eclosions and dead bees (see Supporting Information Section 1.3 for more details on experimental measurements).

Model fitting
Our data collected from the experimental colonies consisted of dayto-day observations of the numbers of bees that had newly eclosed and died in each colony. We assessed four models against these data (see Table 1). To fit a colony's data to solutions of the models, we used discrete stochastic versions of the models with a time interval of one day. The numbers of bees that eclose, die or transition between the two states of each model (the SLS Model and SLS Variant Model each have two variables being Susceptible and Impaired; the Khoury Model has two variables being Hive and Forager; the LA Model has two variables being Larva and Adult) were drawn from a Poisson distribution with means according to the corresponding rates in the deterministic models.
We used the Numerically Integrated State Space (NISS) algorithm to calculate a likelihood for a model and a solution given the data Table 1 The four models we fitted against the empirical data from the experimental bumblebee colonies for eclosion and death rates for each colony algorithm (De Valpine & Hastings 2002). The NISS algorithm calculates the probability of observations for each time step P(O t ), given the observations at the previous time step P(O t-1 ), a model and a solution. This is done by using the model to generate probabilities of all possible states at time step t, given all possible states at the previous time step t − 1.
Given states of the model, probabilities of observations were calculated by summing the Poisson mass functions for the death and eclosion rates from the model for each state. The log likelihood for a set of parameters was calculated by adding the logarithms of the probabilities of observations at every time step and for each colony. The maximum likelihood was then estimated using the adaptive Differential Evolution optimisation algorithm (Brest et al. 2006) with 80 solutions per generation. Optimisation stopped when 25 generations passed with no new best solution. The initial state of the model is given by the initial state of each experimental colony (number of workers) as follows: SLS Model, all initial bees were healthy so S = number of workers, I = 0; Khoury Model, colony split in half with equal numbers of hive and forager bees, H = 0.5 9 number of workers, F = 0.5 9 number of workers; LA Model, A = number of workers, the initial number of larvae per adult was set by a parameter optimised by the Differential Evolution algorithm giving L = A 9 initial_larvae_per_adult.

RESULTS
The SLS Model (Eqns 1 and 2) has fundamentally different dynamics for different values of the parameters. After an initial phase of growth the dynamics will show either: (1) growth for all colonies; (2) multiple outcomes with some colonies growing and others failing dependent on the initial conditions; or (3) all colonies decline and fail (see Fig. 1; Supporting Information Section 2.1). The presence of multiple outcomes relies on the l/(N + φ) term (Supporting Information Section 2.2) which introduces positive density dependence in the model. With a simpler term (see SLS Variant Model) there is still a sharp boundary between growth (Fig. 1a) and failure (Fig. 1c).
The pesticide treated colonies in our empirical experiment followed a similar growth pattern to that observed for SLS Model runs ending in colony failure. By the end of the 42-day experiment, we found a significant difference in colony size between control and treatment colonies (Fig. 2 and data in Tables S1-S6). While all colonies grew at a similar rate during the first 3 weeks, only control colonies continued growing throughout the 42-day study, whereas treatment colonies began to shrink (at 33 days the average colony size of the neonicotinoid treated colonies was and continued to be significantly lower than control colonies; t-test for all comparisons: P < 0.05). Fig. 2 also indicates that there is impaired colony function in the pesticide treated colonies because the birth rates decreased (relative to control colonies) while colonies were still growing, and the death rates also increased during the period in which treated colonies were in decline. This analysis of the data from our experiment thus directly shows that sublethal pesticide exposure decreases colony size after a lagged growth period, and also indicates that this may be due to effects of impairment on colony function rather than direct mortality.
To assess whether colony function is important for explaining colony dynamics, we used a model fitting approach to investigate the patterns of colony sizes, birth rates and death rates found in the treatment colonies. We did this by comparing the fit of the SLS Model with two alternative models. The Khoury Model incorporates lethal stress, but not the impairment and feedback caused by sublethal stress (Khoury et al. 2011). The LA Model incorporates toxic effects from pesticides upon larvae, which could introduce delayed effects on colony size by killing larvae. We fitted the three models to the empirical data using the NISS algorithm (De Valpine & Hastings 2002), which calculates a likelihood value for a model based on all possible trajectories. Parameters were optimised using Differential Evolution (Brest et al. 2006). The best fits found for the SLS Model and the Khoury Model are shown in Fig. 3 (see Fig. S2 for a comparison of the LA Model and SLS Model). Using Akaike weights we selected the best model, and found that the SLS Model describes the data overwhelmingly better with essentially no support for the Khoury Model or the LA Model compared with the SLS Model (Table 2). Of the three models tested, only the SLS model (that incorporates feedback of colony function on birth and death rates) matched the pattern of birth rates decreasing and death rates increasing in the treatment colonies.
This leads us to conclude that colony function is a key element to explain the dynamics of our treatment colonies, and suggests an explanation for the mechanism by which sublethal effects lead to colony failure. Social bee colonies depend on the efficient cooperative performance of multiple individual workers so that essential tasks such as foraging, thermoregulation and brood care, sustain and enhance overall colony function. They have many workers and are able to buffer some effects of stress. However, if too many bees become behaviourally impaired, irrespective of the reason, the colony reaches a tipping point and is set on a path to failure through moderate, but chronic, levels of stress. If the stress level is below this tipping point, the colony can continue growing. There is, con- There is a breaking point (dashed line) between two basins of attraction depending on the numbers of bees present in the colony and the level of impaired bees in the colony; depending on initial conditions colonies will either grow or fail. (c) All colonies can grow at first but will eventually fail. sequently, a critical stress level, where a small change in the amount of stress can mean the difference between colony growth or failure. Next, we investigated whether there is evidence for positive density dependence in the dynamics of our colonies. To do so, we also fitted the data to a variant of the SLS Model in which we replaced the positive density dependence term l/(N + φ) with a linear death rate l. The variant model does not allow for multiple outcomes (see Supporting Information Section 2.2), but there is still a critical stress level for colony failure when impairment is too high (Fig. 1a and c). Because the best fit of the variant model was much worse than the full SLS Model (Fig. S3), with Akaike weights showing essentially no support for the variant model (see Table 2), we can decisively infer that social bee colony dynamics are subject to positive density dependence and thus capable of showing multiple outcomes.
The SLS Model predicts that the level of pesticide is a crucial factor in setting the dynamics of colonies. In our empirical experiment bees received an imidacloprid concentration, in their isolated colonies, near the upper range of that typically found in field realistic conditions (Cresswell 2011;Blacquiere et al. 2012). In a previous experiment (Gill et al. 2012), bumblebees were treated with the same concentration of imidacloprid, but we would expect them to have received a lower level of exposure as they were free to forage for both pollen and nectar in the field. The imidacloprid treated colonies continued to grow throughout the Gill et al. (2012) study, a dynamic consistent with lower exposure in the SLS Model (Fig. 1a). In this study, treated colonies showed a dynamic consistent with higher pesticide exposure in the SLS Model (Fig. 1c), as all treatment colonies failed. Our model predicts that multiple outcomes are possible at a mid-range of imidacloprid exposure. To simulate lower exposure in our model we can reduce parameter b from that found in the fit against the treatment colonies, keeping all other parameters the same, and we observe multiple outcomes in the model's dynamics (see Fig. 4).

DISCUSSION
We have shown here that bumblebee colonies fail when exposed to sustained sublethal levels of pesticide, and that this can be explained (a) Both treatment and control (plot shows mean numbers of workers AE 95% CI) have similar growth at first, but trajectories begin to diverge after approximately 3 weeks. (b, c) Worker birth and death rates (mean number of births and deaths per colony are shown) for treatment and control colonies illustrate the impact of the pesticide on colony health. In the treatment colonies, the birth rate decreases while colony size increases during the first 21 days, and after 21 days the death rate increases while colony size decreases.
by a decrease in colony function. By testing a suite of models against data collected from failing colonies, we are able to make several inferences. We infer that social bee colonies have positive density dependence, they are subject to an Allee effect, and also that there is a critical stress level for the success of a colony such that a small increase in the level of stress can make the difference between failure or success.
Our results provide two explanations as to why it has been so difficult to explain colony losses. First, although we imposed a specific (pesticide) stress in our experiment, the argument about reduced colony function applies to any stressor that reduces the contribution made by individual bees to colony function. This suggests that multifactorial stress can cause colony failure (Vanbergen & The Insect Pollinators Initiative 2013), but that failure is the result of a critical stress level from the accumulation of multiple sublethal factors (e.g. disease, weather and anthropogenic influences) without the need for synergy between these effects. This can explain why finding the link between colony failures and a single specific stress factor has so far proved elusive.
Second, it appears that a number of potential causes for the failure of bee colonies have been dismissed because the presence of the stressor is not a good predictor of colony failure (Cox-Foster  Our study found that the best fitting model to our empirical data was capable of showing multiple outcomes where the fate of the colony can depend on its initial size; however, our experimental results did not directly display these multiple outcomes. It is known to be highly challenging to empirically demonstrate such multiple outcomes (Ives et al. 2008;Mumby et al. 2013). In our case, the difficulty could arise from there being a very precise level of stress at which multiple outcomes exist and since bee colony populations are subject to random occurrences, the level of colony replication required would make such experiments difficult to perform. Even in the event that such experiments could be successfully conducted, the results would still require model fitting (such as that done in this paper) to verify these colonies did not fail due to random chance (process error). The results of such an experiment would in essence simply reconfirm what we have already shown here. The evidence for multiple outcomes that we find here comes from the fact that the SLS Model incorporates feedback from the population density onto the death rate (the Allee effect). By showing that this nonlinear response is important, we expect that any other models that have as strong a fit to our data as the SLS Model, will also show critical stress levels and multiple outcomes.
This study demonstrates two key aspects of how stress on individual bees can disrupt colony function and lead to colony failure. First, a stressor must have a chronic impact (over a period of several weeks) before we see any noticeable effect: meaning that risk assessments of a stressor's impacts must be over a similar time scale. Second, we show how a stressor that impairs colony function can cause an Allee effect which makes colonies especially susceptible to failure from stress at earlier points in their life cycles. This has important implications for the number of colonies that need to be tested to assess the impacts of stressors to adequately measure the proportion of colonies that fail.
The dominance of social bees as crucial pollinators stems primarily from their social organisation: large colony sizes are supported by the efficient coordination of tasks across group members, such that colony performance is better than a collection of uncoordinated individuals. It is intriguing that the social organisation that leads to the success of social bees may also be a key factor in their declines and colony failures. Phase plane demonstrating multiple outcomes in the dynamics of the SLS Model between colony growth and failure. Trajectories denote 20 days of the model starting from different initial conditions (see Methods for more information) and are depicted by arrows coloured according to whether they lead to continued growth (black arrows) or colony failure (red arrows). Parameters were the same as those in Fig. 3 with b = 0.12.