bee colony monitoring at landscape scale in West European intensive agrosystems. J. Apic. Res. 53, 57-66.

Odoux, J.-F., Feuillet, D., Aupinel, P., Loublier, Y., Tasei, J.-N., Mateescu, C., 2012. Territorial biodiversity and consequences on physico-chemical characteristics of pollen collected by honey bee colonies. *Apidologie* 43, 561-575. https://doi.org/10.1007/s13592-012-0125-1

OECD, 2021a. Guidance document on honey bee (Apis mellifera L.) homing flight test, using single oral exposure to sub-lethal doses of test chemical (No. 332), OECD Test Guidelines, Series on Testing and Assessment. Paris, France.

OECD, 2021b. Reports of ring-tests on the Honeybee (Apis mellifera) homing flight test following single exposure to sublethal doses of test chemical, conducted between 2015 and 2019 (No. 333), OECD Test Guidelines, Series on Testing and Assessment. Paris, France.

Oldroyd, B.P., 2007. What's killing American honey bees? PLoS Biol 5, e168. https://doi.org/10.1371/journal.pbio.0050168

Palmer, M.J., Moffat, C., Saranzewa, N., Harvey, J., Wright, G.A., Connolly, C.N., 2013. Cholinergic pesticides cause mushroom body neuronal inactivation in honeybees. Nat. Commun. 4, 1634. https://doi.org/10.1038/ncomms2648

Park, M.G., Blitzer, E.J., Gibbs, J., Losey, J.E., Danforth, B.N., 2015. Negative effects of pesticides on wild bee communities can be buffered by landscape context. Proc. R. Soc. Lond. B Biol. Sci. 282, 20150299. https://doi.org/10.1098/rspb.2015.0299

Parreño, M.A., Alaux, C., Brunet, J.-L., Buydens, L., Filipiak, M., **Henry, M.**, Keller, A., Klein, A.-M., Kuhlmann, M., Leroy, C., Meeus, I., Palmer-Young, E., Piot, N., Requier, F., Ruedenauer, F., Smagghe, G., Stevenson, P.C., Leonhardt, S.D., 2021. Critical links between biodiversity and health in wild bee conservation. Trends Ecol. Evol. https://doi.org/10.1016/j.tree.2021.11.013

Perry, C.J., Søvik, E., Myerscough, M.R., Barron, A.B., 2015. Rapid behavioral maturation accelerates failure of stressed honey bee colonies. Proc. Natl. Acad. Sci. 201422089. https://doi.org/10.1073/pnas.1422089112

Pilling, E., Campbell, P., Coulson, M., Ruddle, N., Tornier, I., 2013. A four-year field program investigating long-term effects of repeated exposure of Honey Bee colonies to flowering crops treated with thiamethoxam. PLoS One 8, e77193. https://doi.org/10.1371/journal.pone.0077193

Poquet, Y., Vidau, C., Alaux, C., 2016. Modulation of pesticide response in honeybees. *Apidologie* 47, 412-426.

Potts, S.G., Biesmeijer, J.C., Kremen, C., Neumann, P., Schweiger, O., Kunin, W.E., 2010. Global pollinator declines: trends, impacts and drivers. Trends Ecol. Evol. 25, 345-353. https://doi.org/10.1016/j.tree.2010.01.007

Powney, G.D., Carvell, C., Edwards, M., Morris, R.K.A., Roy, H.E., Woodcock, B.A., Isaac, N.J.B., 2019. Widespread losses of pollinating insects in Britain. Nat. Commun. 10, 1018. https://doi.org/10.1038/s41467-019-08974-9

Prisco, G.D., Cavaliere, V., Annoscia, D., Varricchio, P., Caprio, E., Nazzi, F., Gargiulo, G., Pennacchio, F., 2013. Neonicotinoid clothianidin adversely affects insect immunity and promotes replication of a viral pathogen in honey bees. Proc. Natl. Acad. Sci. 110, 18466-18471. https://doi.org/10.1073/pnas.1314923110

Requier, F., Odoux, J.-F., **Henry, M.**, Bretagnolle, V., 2017. The carry-over effects of pollen shortage decrease the survival of honeybee colonies in farmlands. J. Appl. Ecol. 54, 1161-1170. https://doi.org/10.1111/1365-2664.12836

Requier, F., Odoux, J.-F., Tamic, T., Moreau, N., **Henry, M.**, Decourtye, A., Bretagnolle, V., 2015. Honey bee diet in intensive farmland habitats reveals an unexpectedly high flower richness and a major role of weeds. Ecol. Appl. 25, 881-890. https://doi.org/10.1890/14-1011.1

Requier, F., Rome, Q., Chiron, G., Decante, D., Marion, S., Menard, M., Muller, F., Villemant, C., **Henry, M.**, 2018. Predation of the invasive Asian hornet affects foraging activity and survival probability of honey bees in Western Europe. J. Pest Sci. https://doi.org/10.1007/s10340-018-1063-0

Requier, F., Rome, Q., Villemant, C., **Henry, M.**, 2020. A biodiversity-friendly method to mitigate the invasive Asian hornet's impact on European honey bees. J. Pest Sci. 93, 1-9. https://doi.org/10.1007/s10340-019-01159-9

Retschnig, G., Williams, G.R., Mehmann, M.M., Yañez, O., de Miranda, J.R., Neumann, P., 2014. Sex-specific differences in pathogen susceptibility in Honey Bees (Apis mellifera). PLoS One 9, e85261. https://doi.org/10.1371/journal.pone.0085261

Ribbands, C.R., 1953. The behaviour and social life of honeybees., Bee Research Association Limited, London. ed.

Rollin, O., Bretagnolle, V., Decourtye, A., Michel, N., Vaissière, B.E., **Henry, M.**, 2013. Differences of floral
resource use between honey bees and wild bees in an intensive farming system. Agric. Ecosyst. Environ. 179, 78–86. https://doi.org/10.1016/j.agee.2013.07.007

Rollin, O., Bretagnolle, V., Fortel, L., Guilbaud, L., **Henry, M.**, 2015. Habitat, spatial and temporal drivers of diversity patterns in a wild bee assemblage. Biodivers. Conserv. 1–20. https://doi.org/10.1007/s10531-014-0852-x

Rollin, O., Pérez-Méndez, N., Bretagnolle, V., **Henry, M.**, 2019. Preserving habitat quality at local and landscape scales increases wild bee diversity in intensive farming systems. Agric. Ecosyst. Environ. 275, 73–80. https://doi.org/10.1016/j.agee.2019.01.012

Rondeau, G., Sánchez-Bayo, F., Tennekes, H.A., Decourtye, A., Ramírez-Romero, R., Desneux, N., 2014. Delayed and time-cumulative toxicity of imidacloprid in bees, ants and termites. Sci. Rep. 4. https://doi.org/10.1038/srep05566

Rumkee, J.C.O., Becher, M.A., Thorbek, P., Kennedy, P.J., Osborne, J.L., 2015. Predicting honeybee colony failure: using the BEEHAVE model to simulate colony responses to pesticides. Environ. Sci. Technol. 49, 12879–12887. https://doi.org/10.1021/acs.est.5b03593

Sánchez-Bayo, F., Wyckhuys, K.A.G., 2019. Worldwide decline of the entomofauna: A review of its drivers. Biol. Conserv. 232, 8–27. https://doi.org/10.1016/j.biocon.2019.01.020

Schatz, B., Drossart, M., **Henry, M.**, Geslin, B., Allier, F., Savajol, C., Gérard, M., Michez, D., 2021. Pollinator conservation in the context of global changes with a focus on France and Belgium. Acta Oecologica 112, 103765. https://doi.org/10.1016/j.actao.2021.103765

Schmickl, T., Crailsheim, K., 2007. HoPoMo: A model of honeybee intracolonial population dynamics and resource management. Ecol. Model. 204, 219–245. https://doi.org/10.1016/j.ecolmodel.2007.01.001

Schmickl, T., Crailsheim, K., 2002. How honeybees (Apis mellifera L.) change their broodcare behaviour in response to non-foraging conditions and poor pollen conditions. Behav. Ecol. Sociobiol. 51, 415–425. https://doi.org/10.1007/s00265-002-0457-3

Steffan-Dewenter, I., Potts, S.G., Packer, L., 2005. Pollinator diversity and crop pollination services are at risk. Trends Ecol. Evol. 20, 651–652.

Stillman, R.A., Railsback, S.F., Giske, J., Berger, U., Grimm, V., 2015. Making predictions in a changing world: the benefits of individual-based ecology. BioScience 65, 140–150. https://doi.org/10.1093/biosci/biu192

Tison, L., Hahn, M.-L., Holtz, S., Rößner, A., Greggers, U., Bischoff, G., Menzel, R., 2016. Honey Bees' behavior is impaired by chronic exposure to the neonicotinoid thiacloprid in the field. Environ. Sci. Technol. 50, 7218–7227. https://doi.org/10.1021/acs.est.6b02658

Torres, D.J., Ricoy, U.M., Roybal, S., 2015. Modeling honey bee populations. PLOS ONE 10, e0130966. https://doi.org/10.1371/journal.pone.0130966

van der Sluijs, J.P., Simon-Delso, N., Goulson, D., Maxim, L., Bonmatin, J.-M., Belzunces, L.P., 2013. Neonicotinoids, bee disorders and the sustainability of pollinator services. Curr. Opin. Environ. Sustain., Open issue 5, 293–305. https://doi.org/10.1016/j.cosust.2013.05.007

Vanbergen, A.J., the Insect Pollinators Initiative, 2013. Threats to an ecosystem service: pressures on pollinators. Front. Ecol. Environ. 11, 251–259. https://doi.org/10.1890/120126

Vandame, R., Meled, M., Colin, E., Belzunces, L.P., 1995. Alteration of the homing-flight in the honey bee *Apis mellifera* L. exposed to sublethal dose of deltamethrin. Environ. Toxicol. Chem. 14, 855–860.

Vereecken, N.J., Weekers, T., Leclercq, N., De Greef, S., Hainaut, H., Molenberg, J.-M., Martin, Y., Janssens, X., Noël, G., Pauly, A., Roberts, S.P.M., Marshall, L., 2021. Insect biomass is not a consistent proxy for biodiversity metrics in wild bees. Ecol. Indic. 121, 107132. https://doi.org/10.1016/j.ecolind.2020.107132

Vidau, C., Diogon, M., Aufauvre, J., Fontbonne, R., Viguès, B., Brunet, J.-L., Texier, C., Biron, D.G., Blot, N., El Alaoui, H., Belzunces, L.P., Delbac, F., 2011. Exposure to sublethal doses of fipronil and thiacloprid highly increases mortality of honeybees previously infected by *Nosema ceranae*. PLoS One 6, e21550. https://doi.org/10.1371/journal.pone.0021550

Wharton, K.E., Dyer, F.C., Huang, Z.Y., Getty, T., 2007. The honeybee queen influences the regulation of colony drone production. Behav. Ecol. 18, 1092–1099. https://doi.org/10.1093/beheco/arm086

Winfree, R., 2010. The conservation and restoration of wild bees. Ann. N. Y. Acad. Sci. 1195, 169–197. https://doi.org/10.1111/j.1749-6632.2010.05449.x

Winston, M.L., 1987. The biology of the honey bee. Harvard Univ Press, Cambridge, USA.

Winston, M.L., Dropkin, J.A., Taylor, O.R., 1981. Demography and life history characteristics of two honey bee races (*Apis mellifera*). Oecologia 48, 407–413.

Wintermantel, D., Odoux, J.-F., Decourtye, A., **Henry, M.**, Allier, F., Bretagnolle, V., 2020. Neonicotinoid-induced mortality risk for bees foraging on oilseed rape nectar persists despite EU moratorium. Sci. Total Environ.
704, 135400. https://doi.org/10.1016/j.scitotenv.2019.135400

Wojcik, V.A., Morandin, L.A., Davies Adams, L., Rourke, K.E., 2018. Floral resource competition between Honey Bees and wild bees: is there clear evidence and can we guide management and conservation? Environ. Entomol. 47. https://doi.org/10.1093/ee/nvy077

Zacepins, A., Brusbardis, V., Meitalovs, J., Stalidzans, E., 2015. Challenges in the development of Precision Beekeeping. Biosyst. Eng. 130, 60–71. https://doi.org/10.1016/j.biosystemseng.2014.12.001

Zattara, E.E., Aizen, M.A., 2021. Worldwide occurrence records suggest a global decline in bee species richness. One Earth 4, 114–123. https://doi.org/10.1016/j.oneear.2020.12.005
# 8. Annexes

## Liste des Annexes :

8.1. Article Henry et al. 2012, Science (Chapitre 4)

8.2. Article Henry et al. 2014, Nature Communications (Chapitre 4)

8.3. Article Henry et al. 2015, Proceedings of the Royal Society B (Chapitre 4)

8.4. Article Henry et al. 2017, Apidologie (Chapitre 5)

8.5. Article Henry et Rodet 2018, Scientific Reports (Chapitre 6)

8.6. Curriculum vitae détaillé
# A Common Pesticide Decreases Foraging Success and Survival in Honey Bees

Mickaël Henry,<fnref n="1" />,²* Maxime Béguin,²,³ Fabrice Requier,⁴,⁵ Orianne Rollin,²,⁶ Jean-François Odoux,⁵ Pierrick Aupinel,⁵ Jean Aptel,<fnref n="1" />,² Sylvie Tchamitchian,<fnref n="1" />,² Axel Decourtye²,⁶

Nonlethal exposure of honey bees to thiamethoxam (neonicotinoid systemic pesticide) causes high mortality due to homing failure at levels that could put a colony at risk of collapse. Simulated exposure events on free-ranging foragers labeled with a radio-frequency identification tag suggest that homing is impaired by thiamethoxam intoxication. These experiments offer new insights into the consequences of common neonicotinoid pesticides used worldwide.

Colony collapse disorder (CCD) is a recent, pervasive syndrome affecting honey bee (Apis mellifera) colonies in the Northern hemisphere, which is characterized by a sudden disappearance of honey bees from the hive (1). Multiple causes of CCD have been proposed, such as pesticides, pathogens, parasites, and natural habitat degradation (2–4). However, the relative contribution of those stressors in CCD events remains unknown. Some scientists and beekeepers suspect pesticides to hold a central place in colony-weakening processes (1) or at least in interaction with other stressors (5, 6). In modern cereal farming systems, honey bees are readily exposed to pesticides because they rely heavily on common blooming crops, such as oilseed rape (Brassica napus), maize (Zea mays), or sunflower (Helianthus annuus), that are now routinely treated against insect pests (3). Systemic pesticides in particular diffuse throughout all the tissues as plants grow up, eventually contaminating nectar and pollen (7). Foraging honey bees are therefore directly exposed, but so is the rest of the colony as returning foragers store or exchange contaminated material with hive conspecifics (7, 8). Those exposure pathways are of important concern, and pesticide manufacturers pay special attention to reduce nonintentional intoxications in field conditions. Pesticide authorization procedures now require running mortality surveys to ensure that doses encountered in the field remain below lethal levels for honey bees.

However, a growing body of evidence shows that sublethal doses—doses that do not entail direct mortality—still have the potential to induce a variety of behavioral difficulties in foraging honey bees, such as memory and learning

dysfunctions and alteration of navigational skills (9). Neonicotinoid pesticides used to protect crops against aphids and other sap-sucking insects are especially liable to provoke such behavioral troubles. They are highly potent and selective agonists of nicotinic acetylcholine receptors, which are important excitatory neurotransmitter receptors in insects (10, 11). Effects of sublethal neonicotinoid exposures in honey bees may include abnormal foraging activity (12–14), reduced olfactory memory and learning performance (15–17), and possibly impaired orientation skills (18). Yet, the consequences of such behavioral difficulties on the fate of free-ranging foragers and on colony dynamics are extremely difficult to assess in the field and remain poorly investigated.

In this study, we tested the hypothesis that a sublethal exposure to a neonicotinoid indirectly increases hive death rate through homing failure in foraging honey bees. We focused our attention on thiamethoxam, a recently marketed neonicotinoid substance (19) currently being authorized in an increasing number of countries worldwide for the protection of oilseed rape, maize, and other blooming crops foraged by honey bees. We proceeded in two steps. First, we assessed mortality induced by homing failure ($m_{hf}$) in exposed foragers. This was achieved by monitoring free-ranging honey bees with radio-frequency identification (RFID) tagging technology (14, 20). Second, we assessed the extent to which $m_{hf}$, in combination with natural

forager mortality, may upset colony dynamics. For that purpose, $m_{hf}$ was introduced into a model of honey bee population dynamics (21).

We used a custom-made RFID device (20) to monitor the fate of 653 individual free-ranging foragers in the course of four separate treatment-versus-control homing experiments (22). The study was conducted in an intensive cereal farming system of western France, as a part of the ECOBEE monitoring facility (Zone Atelier Plaine et Val de Sèvre, Centre d'Études Biologiques de Chizé) and in a suburban area in Avignon, southern France. To simulate daily intoxication events, foragers received a field-realistic, sublethal dose of thiamethoxam (a real dose of 1.34 ng in a 20-µl sucrose solution) and were released away from their colony with a microchip glued on their thorax (Fig. 1A). RFID readers placed at the hive entrance (Fig. 1B) were set to detect on a continual basis tagged honey bees going through the entrance. Mortality due to postexposure homing failure, $m_{hf}$, was then derived from the proportion of nonreturning foragers. To further discriminate $m_{hf}$ from other causes of homing failure in treated foragers—such as natural mortality, predation, or handling stress—we simultaneously released equal numbers of control foragers fed with an untreated sucrose solution. Hence, $m_{hf}$ was calculated as the proportion of nonretuning treated foragers relative to expectations given by the proportion of returning control foragers. Depending on the experiment, tagged honey bees were released up to 1 km away from their respective colony, a distance usually covered by foragers during normal foraging flights (23). Experiments were conducted on individuals from three different colonies (22).

Our strategy was not to get an estimate of $m_{hf}$ per se. Instead, we assessed its upper and lower bounds, depending on whether foragers were familiar or not with the foraging site in which they might get intoxicated. Indeed, one might expect that foragers familiar with the pathway back to the colony are less prone to homing failure than are unfamiliar foragers. Under field conditions, many foragers are probably familiar with the pathway back to the colony because they repeatedly forage on the same

<figure><img src="image_33.png" /><figcaption>**Fig. 1.** Honey bee RFID monitoring equipment. (A) A pollen-forager honey bee fitted with a 3-mg RFID tag. (B) A hive entrance equipped with RFID readers for detecting returning marked foragers.</figcaption></figure>

<footnote n="1">¹INRA (Institut National de la Recherche Agronomique), UR406 Abeilles et Environnement, F-84914 Avignon, France. ²UMT Protection des Abeilles dans l'Environnement, Site Agroparc, F-84914 Avignon, France. ³Association pour le Développement de l'Apiculture Provençale, F-13626 Aix-en-Provence, France. ⁴Centre d'Etudes Biologiques de Chizé, CNRS (USC-INRA 1339), UPR1934, F-79360 Beauvoir-sur-Niort, France. ⁵INRA, UE1255, UE Entomologie, F-17700 Surgères, France. ⁶Association de Co-ordination Technique Agricole, Site Agroparc, F-84914 Avignon, France.</footnote>

<footnote>*To whom correspondence should be addressed. E-mail: mickael.henry@avignon.inra.fr</footnote>
site (24). However, many others are unfamiliar, too. Those include young honey bees at the onset of foraging, scouting honey bees that look for new food sources, and foragers newly recruited by scouting bees on the basis of the dance information (25). Most importantly, systemic pesticides such as thiamethoxam are readily present in the nectar and pollen when flowering starts and receive the first visits of honey bees, hitherto unfamiliar with this newly available food source.

To account for individuals' past foraging experience, we conducted two distinct homing experiments. Experiment 1 simulated intoxication at a familiar foraging site, and experiment 2 at a random site regarding past foraging experience. These experiments were assumed to return the lower and upper bounds of $m_{hf}$, respectively.

In experiment 1, we referred to as “familiar” foragers those foragers for which we could make sure they covered at least once the pathway from the release site back to the colony. For that purpose, we selectively captured foragers returning to the colony with pollen loads from a known location and subsequently released them at that location. To ascertain pollen origin, we sowed beforehand a 1-ha field with scorpion weed *Phacelia tanacetifolia*, a highly attractive floral resource with bright blue pollen that is easily recognizable (26). Given that no other phacelia fields occurred in the area, we could ensure that phacelia-carrying honey bees came back from our experimental field. The colony was specifically placed 1 km away from the field for subsequent forager release (Fig. 2). In

experiment 2, we used the non-phacelia pollen foragers. They were released in equal groups at six sites equally spaced on a 1-km circle around the colony (Fig. 2). Following that design, release sites were considered as random locations regarding the past experience of foragers.

Both experiments 1 and 2 evidenced substantial mortality due to postexposure homing failure, $m_{hf}$, with the proportion of treated foragers returning to the colony being significantly lower than that of control foragers (exact binomial tests, $P = 0.033$ and $P < 0.001$, respectively) (Fig. 3 and table S1). Additionally, $m_{hf}$ was greater in treated foragers that tended to be unfamiliar with the foraging site, as indicated by their significantly lower homing proportions as compared with familiar foragers (exact binomial tests, $P < 0.001$). Experiments 1 and 2 returned $m_{hf}$ estimates of 0.102 and 0.316, respectively, potentially setting the lower and upper bounds for real $m_{hf}$ values. In other words, 10.2 to 31.6% of exposed honey bees would fail to return to their colony when foraging in treated crops on a daily basis. For the sake of comparison, foragers live ~6.5 days and therefore die at an average rate of $1/6.5 = 0.154$ individual day$^{-1}$ (27). Therefore, the probability that a forager would die because of homing failure during a day spent foraging on treated crops (up to 0.316) may attain twice the probability this same forager has to die naturally that day (~0.154).

Such an additional mortality might represent a heavy burden to bear for colonies exposed to treated crops in their environment. When implementing $m_{hf}$ into a honey bee population dynamics model (21), all the tested scenarios predicted a major deviation from the expected dynamic (Fig. 4). In our simulations, we considered the evolution of a typical colony during the first 3 months of a beekeeping season, encompassing the oilseed rape blooming period, which was April to May in our study area (22). At this time of the year, colonies emerge from the wintering period. Population size is rather low (<20,000 individuals) and gradually expands in order to rapidly increase food storage and ensure colony sustainability. The daily egg-laying rate of the queen is a critical parameter in this colony dynamic because it determines the daily egg-hatching rate and in turn the rate at which honey bees working in the hive will be replaced as they become themselves foragers. We simulated three scenarios with realistic levels of egg-laying rate (28), namely a rate allowing for a normal colony development (Fig. 4A), a rate ensuring equilibrium population (Fig. 4B), and a slightly deficient rate forcing the population to stabilize at a lower size (Fig. 4C). In each case, we also computed the expected trends if most foragers (90%) were exposed to nectar of treated oilseed rape each day and therefore had a natural mortality increased by a homing failure probability $m_{hf}$. Regardless of the queens' egg-laying rate, populations from colonies exposed to the treated nectar would follow a marked decline during the

<figure><img src="image_34.png" /><figcaption>**Fig. 2.** Study area and location of honey bee release sites relative to the colony hive in experiments 1 and 2.</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>**Fig. 3.** Cumulative homing probability of foragers released 1 km away from the hive. Temporal gaps denote the nighttime between the first and second days of release. (A) Homing experiment 1 was carried out with foragers familiar with the release site, and (B) experiment 2 with foragers released at random sites regarding their past experience. In both cases, treated honey bees that received a nonlethal dose of thiamethoxam returned to the hive in significantly lower proportions than did control honey bees (table S1).</figcaption></figure>
<figure><img src="image_36.png" /><figcaption>**Fig. 4.** Comparison of honey bee population dynamics between simulated colonies exposed to thiamethoxam (red lines) or not exposed (blue lines), following six demographic scenarios. *L* is the queen's daily laying rate (eggs per day). "Exp" is the proportion of foragers exposed to treated crops during the day. The nonexposed colony follows either (**A** and **D**) a normal development trajectory (at *L* = 2000), (**B** and **E**) an equilibrium dynamic (*L* = 1800), or (**C** and **F**) a slightly declining trajectory (*L* = 1600). Shaded areas delineate the exposure period (for example, oilseed rape). Pairs of trajectories in exposed colonies were obtained with the lower and upper bounds of homing failure mortality (0.102 and 0.316) in order to delineate the best and worst estimates for population dynamics, respectively. Dotted lines extend the declining trajectory expected for a sustained exposure. [Simulations derive from demographic models in (21)]</figcaption></figure>

blooming period and would hardly recover after-
ward (Fig. 4, A to C). When combined with
natural forager mortality, m<sub>hf</sub> raised total forager
death rate up to a point that could hardly be
compensated for by the rate at which new for-
agers are recruited. In the worse scenarios, pop-
ulations would fall down to 5000 individuals,
which is the lowest level one can usually ob-
serve in current beekeeping practices. With an
exposure rate reduced to 50% of foragers ex-
posed to treated nectar each day (Fig. 4, D to F),
the model still predicts a major deviation from
normal conditions.

In an attempt to verify the applicability of
these results to other contexts, we repeated two
additional experiments with two different colo-
nies (fig. S2 and table S1). In experiment 3,
we tested whether m<sub>hf</sub> was still significant when
exposure occurred in the least challenging sit-
uation: in the direct vicinity of the colony and
with honey bees familiar with the foraging site.
We repeated experiment 1 with phacelia foragers
captured from a beehive placed at the phacelia
field margin and released from inside the phacelia

field, only 70 m away. Homing failure (m<sub>hf</sub> = 0.061) (fig. S2A and table S1) was much reduced as compared with that of experiment 1 (m<sub>hf</sub> = 0.102) but was still significant (exact binomial test, *P* = 0.003). In experiment 4, we transposed experiment 2 into a different landscape. A beehive was placed in a suburban area in southern France, including a mosaic of mixed farming fields and orchards of moderate size. Foragers were released 1 km away at six equidistant sites. Homing failure (m<sub>hf</sub> = 0.098) (fig. S2B and table S1) was significant as well (exact binomial test, *P* = 0.029) but much smaller than in experiment 2 (m<sub>hf</sub> = 0.316).

Our study clearly demonstrates that expo-
sure of foragers to nonlethal but commonly
encountered doses of thiamethoxam can affect
forager survival, with potential contributions to
collapse risk. Furthermore, the extent to which
exposures affect forager survival appears de-
pendent on the landscape context and the prior
knowledge of foragers about this landscape.
Higher risks are observed when the homing task
is more challenging. As a consequence, impact

studies are likely to severely underestimate sub-
lethal pesticide effects when they are conducted
on honey bee colonies placed in the immediate
proximity of treated crops. This study raises im-
portant issues concerning exposed solitary bee
species, whose population dynamics are prob-
ably less resilient to forager disappearance than
are honey bee colonies.

**References and Notes**

1. B. P. Oldroyd, *PLoS Biol.* **5**, e168 (2007).

2. D. L. Cox-Foster et al., *Science* **318**, 283 (2007).

3. C. A. Mullin et al., *PLoS ONE* **5**, e9754 (2010).

4. D. Naug, *Biol. Conserv.* **142**, 2369 (2009).

5. C. Alaux et al., *Environ. Microbiol.* **12**, 774 (2010).

6. C. Vidau et al., *PLoS ONE* **6**, e21550 (2011).

7. A. Rortais, G. Arnold, M. P. Halm, F. Touffet-Briens, *Apidologie* **36**, 71 (2005).

8. C. H. Krupke, G. J. Hunt, B. D. Eitzer, G. Andino, K. Given, *PLoS ONE* **7**, e29268 (2012).

9. N. Desneux, A. Decourtye, J. M. Delpuech, *Annu. Rev. Entomol.* **52**, 81 (2007).

10. M. Tomizawa, J. E. Casida, *Annu. Rev. Entomol.* **48**, 339 (2003).

11. N. S. Millar, I. Denholm, *Invert. Neurosci.* **7**, 53 (2007).

12. E. C. Yang, Y. C. Chuang, Y. L. Chen, L. H. Chang, *J. Econ. Entomol.* **101**, 1743 (2008).

13. M. E. Colin et al., *Arch. Environ. Contam. Toxicol.* **47**, 387 (2004).

14. C. W. Schneider, J. Tautz, B. Grünewald, S. Fuchs, *PLoS ONE* **7**, e30023 (2012).

15. A. Decourtye et al., *Pestic. Biochem. Physiol.* **78**, 83 (2004).

16. A. Decourtye, E. Lacassie, M. H. Pham-Delègue, *Pest Manag. Sci.* **59**, 269 (2003).

17. A. Decourtye, J. Devillers, S. Cluzeau, M. Charreton, M. H. Pham-Delègue, *Ecotoxicol. Environ. Saf.* **57**, 410 (2004).

18. L. Bortolotti et al., *Bull. Insectology* **56**, 63 (2003).

19. P. Maienfisch et al., *Pest Manag. Sci.* **57**, 165 (2001).

20. A. Decourtye et al., *Ecotoxicology* **20**, 429 (2011).

21. D. S. Khoury, M. R. Myerscough, A. B. Barron, *PLoS ONE* **6**, e18491 (2011).

22. Materials and methods are available as supporting online material on *Science* Online.

23. I. Steffan-Dewenter, A. Kuhn, *Proc. R. Soc. B Biol. Sci.* **270**, 569 (2003).

24. R. Menzel, R. J. De Marco, U. Greggers, *J. Comp. Physiol. A Neuroethol. Sens. Neural Behav. Physiol.* **192**, 889 (2006).

25. J. R. Riley, U. Greggers, A. D. Smith, D. R. Reynolds, R. Menzel, *Nature* **435**, 205 (2005).

26. M. Henry et al., *Ecol. Modell.* **225**, 103 (2012).

27. R. Dukas, *Insectes Soc.* **55**, 252 (2008).

28. T. Schmickl, K. Crailsheim, *Ecol. Modell.* **204**, 219 (2007).

**Acknowledgments:** This study was funded by the European Community program (797/2004) for French beekeeping coordinated by the French Ministry of Agriculture (convention FranceAgriMer 11-45R). Special thanks go to L. Belzunces, J. L. Brunet, B. Vaissière, A. Maisonnasse, D. Fortini, Y. Le Conte, C. McDonnell, P. Jourdan, J. P. Vermandère, and V. Bretagnolle for valuable help and corrections, as well as three anonymous reviewers for useful comments on the manuscript. We are grateful to M. Grijolot for allowing experiments to take place in his fields. RFID devices were designed by Tag Tracing Solutions, Valence, France, and adapted to honey bees with the help of J. Fourrier, F. Brun, J. Devillers, M. Gauthier, and the ECOBEE monitoring facility. Data are in the supporting online material.

**Supporting Online Material**

www.sciencemag.org/cgi/content/full/science.1215039/DC1

Materials and Methods

Figs. S1 and S2

Table S1

Reference (29)

Database S1

10 October 2011; accepted 5 March 2012

Published online 29 March 2012;

10.1126/science.1215039
<figure><img src="image_37.png" /></figure>

ARTICLE

Received 25 Apr 2014 | Accepted 10 Jun 2014 | Published 10 Jul 2014

DOI: 10.1038/ncomms5359

# Pesticide risk assessment in free-ranging bees is weather and landscape dependent

Mickaël Henry<fnref n="1" />,², Colette Bertrand³,⁴, Violette Le Féon<fnref n="1" />,², Fabrice Requier⁵,⁶, Jean-François Odoux⁵, Pierrick Aupinel⁵, Vincent Bretagnolle⁶ & Axel Decourtye²,⁷,⁸

The risk assessment of plant protection products on pollinators is currently based on the evaluation of lethal doses through repeatable lethal toxicity laboratory trials. Recent advances in honeybee toxicology have, however, raised interest on assessing sublethal effects in free-ranging individuals. Here, we show that the sublethal effects of a neonicotinoid pesticide are modified in magnitude by environmental interactions specific to the landscape and time of exposure events. Field sublethal assessment is therefore context dependent and should be addressed in a temporally and spatially explicit way, especially regarding weather and landscape physiognomy. We further develop an analytical Effective Dose (ED) framework to help disentangle context-induced from treatment-induced effects and thus to alleviate uncertainty in field studies. Although the ED framework involves trials at concentrations above the expected field exposure levels, it allows to explicitly delineating the climatic and landscape contexts that should be targeted for in-depth higher tier risk assessment.

<footnote n="1">¹INRA, UR406 Abeilles et Environnement, F-84914 Avignon, France. ²UMT Protection des Abeilles dans l'Environnement, Site Agroparc, F-84914 Avignon, France. ³INRA, UR980 SAD Paysage, F-35042 Rennes, France. ⁴CNRS et Université de Rennes 1, UMR6553 Ecosystèmes Biodiversité Evolution, F-35042 Rennes, France. ⁵INRA, UE1255, UE Entomologie, F-17700 Surgères, France. ⁶Centre d'Etudes Biologiques de Chizé, CNRS (USC-INRA 1339), UPR1934, F-79360 Beauvoir-sur-Niort, France. ⁷Association de Coordination Technique Agricole, Site Agroparc, F-84914 Avignon, France. ⁸ITSAP-Institut de l'abeille, Site Agroparc, F-84914 Avignon, France. Correspondence and requests for materials should be addressed to M.H. (email: mickael.henry@avignon.inra.fr).</footnote>
The current assessment of new plant protection products on bees requires a tiered sequence, the first tier covering lethal toxicity trials under controlled laboratory conditions, whereas higher tiers imply surveys under (semi-) field conditions. The laboratory trials are aimed at determining toxicity end points based on median Lethal Dose (LD₅₀) estimation following acute oral or contact exposure. LD₅₀ is designed to be a repeatable mortality parameter, returned by easily reproducible testing schemes. Conversely, there is currently a lack of standardized procedure, both at first and higher tiers, to quantify sublethal effects of pesticides, that is, behavioural impairments that may originate from non-lethal exposures. This assessment shortcoming was further underlined by recent studies on honeybees (Apis mellifera) and bumblebees (Bombus terrestris) carried out in free-ranging conditions¹⁻⁴. Subsequently, the European governments have banned for 2 years, starting on December 2013, the use of three neonicotinoid pesticides as seed-coating treatment for some common flowering crops⁵, while calling for more studies in field-realistic conditions. A key issue currently uncovered by risk assessment schemes is whether the levels of pesticide residues found in nectar and pollen, and certified as sublethal in laboratory, can lead to adverse effects at the colony level on the field. To fill this gap, the European Food Safety Agency⁶ has identified a need to set standards for the measurement of sublethal effects in free-ranging bees, for example, using homing experiments at home range scale¹,⁷,⁸, although scaling up from individuals to colony remains a major conceptual and technical challenge⁹.

Homing flight performance has been considered as a candidate criterion to reveal pesticide sublethal effects in free-ranging honeybees because it may provide useful information on a range of sublethal parameters⁶. Successful homing flight is contingent to the proper integration of multiple physiological and cognitive functions, including among others, orientation, spatial memory, associative learning, energetic metabolism and muscular flight activity. Therefore, homing experiments cover a wide spectrum of biological mechanisms of pesticides are liable to interfere with (ref. 6). For instance, neonicotinoids at sublethal concentrations alter the neurophysiological properties of Kenyon cells, which are the major neuronal components of the honeybee mushroom bodies where multisensory information signals are integrated¹⁰. This in turn will have detrimental effects on cognitive tasks such as memory retrieval and orientation, as revealed by increased Homing Failure (HF)¹,⁷,⁸. Homing experiments are also relevant for field assessment because they are performed at the home range scale and because an increased probability of HF may be converted into a mortality parameter and implemented into exploratory models of colony population dynamics¹,¹¹,¹².

On the other hand, homing flights are conditioned by keystone environmental correlates of bee navigation, which are liable to interfere with the assessment of neonicotinoid sublethal effects in free-ranging conditions. Herein, Navigation Evironmental Correlate, thereafter abbreviated as NECs, refer to those extrinsic drivers of bee navigation and flight performances, comprising weather conditions (temperature, cloudiness or visibility of sun) and salient landmarks acting as visual cues for orientation¹³⁻²². Homing honeybees rely on a network of learned interconnected vectors for orientation¹⁸, conceptualized as memorized distances and directions among known locations. They further have the ability to derive novel shortcut vectors from the integration of two or more previously learned vector flights¹⁹. Honeybees retrieve the directional components of vector flights based on the sun as a celestial cue, but also on known landmarks whenever the sun compass or any spot of polarized light in the sky are masked by clouds¹⁶,¹⁹,²¹. Navigation in honeybees is further driven by temperature. Ambient temperature affects metabolism and flight

speed¹⁷, for instance with lower wingbeat frequency and higher energy expenditure per wingbeat as ambient temperature decreases from 37 °C to 19 °C (ref. 22). Overall, temperature jointly with the availability of the sun compass and of the visual landmarks as orientation cues are regarded as key NECs.

In this study, using HF as a paradigm for field assessment of pesticide risks to honeybees, we tested whether the neonicotinoid thiamethoxam interacts with NECs, potentially changing the latent sublethal effects among experimental contexts. First, we establish a generic dose-response function of HF, covering a broad range of sublethal exposure levels (non-exposed control bees and treatments from LD₅₀/12 to LD₅₀/2, with the LD₅₀ estimated at 5 ng per bee (ref. 23)). Second, we test whether the specific weather, as well as landscape context at the time and place of individuals' release, may modulate the expression of generic dose-response function through an interaction effect with thiamethoxam exposure.

As the dose-response function proved to be highly context dependant, we finally propose analytical solutions to explicitly handle the resulting uncertainty. Following the LD principle, we computed the ED required for reaching a given sublethal effect size in free-ranging bees across a range of contexts. The ED framework is designed (i) to assess how the risk to honeybees do vary with weather and landscape contexts and (ii) to explicitly delineate the areas where the effect is expected to be strongest and therefore should be targeted for in-depth higher tier risk assessment. This approach will help disentangle environment-induced effects from treatment-induced effects in field surveys and thus will contribute to alleviate uncertainty in the assessment of pesticide side effects to bees.

## Results

**HF increases with oral exposure doses.** We carried out 60 homing trials (see Methods) involving 978 honeybee foragers from nine different colonies, released at 22 distinct locations

<figure><img src="image_38.png" /></figure>

**Figure 1 | Generic dose-response function of the risk of homing failure in exposed honeybees.** (a) The predicted dose-response function (continuous line ± estimate s.e. as dashed lines) derived from a binomial logit LME model, specifying random intercepts for colony and release point effects, with *n* = 978 individuals split into 60 homing trials from 22 release points around nine colonies. Quartile boxes and dots overlaying predictions display observed homing failure proportions for trials with >15 individuals. (b) A posteriori statistical power analysis, appraising the ability of the LME to detect a significant homing failure dose-response with the current sample size. Statistical power was measured as the proportion of significant effects out of 1,000 simulated homing failure data sets randomly recomputed with a constrained effect size (percentage increase at 1 ng). The observed effect size (46% increase at 1 ng, filled dot) had a 99% chance of being detected given our sample size.
<table><caption>**Table 1 | Summary of the LME binomial models performed to assess the thiamethoxam dose-response function of honeybee homing failure and its interactions with NECs.**</caption><thead><tr><th>Model parameter</th><th>Complete model estimate ± s.e.</th><th>Z</th><th>P-value</th><th>Multimodel averaged estimate ± s.e.</th><th>Relative importance</th></tr></thead><tbody><tr><td>Intercept</td><td>-1.198 ± 0.317</td><td>-3.780</td><td><strong>&lt;0.001</strong></td><td>-1.151 ± 0.327</td><td></td></tr><tr><td>Dose</td><td>0.689 ± 0.194</td><td>3.546</td><td><strong>&lt;0.001</strong></td><td>0.665 ± 0.193</td><td>100%</td></tr><tr><td>Landmarks</td><td>-0.204 ± 0.300</td><td>-0.679</td><td>0.497</td><td>-0.162 ± 0.268</td><td>60%</td></tr><tr><td>Temperature</td><td>-1.486 ± 0.268</td><td>-5.548</td><td><strong>&lt;0.001</strong></td><td>-1.537 ± 0.261</td><td>100%</td></tr><tr><td>Dose × landmarks</td><td>1.024 ± 0.468</td><td>2.187</td><td><strong>0.029</strong></td><td>0.827 ± 0.472</td><td>42%</td></tr><tr><td>Dose × temperature</td><td>-1.197 ± 0.445</td><td>-2.692</td><td><strong>0.007</strong></td><td>-1.002 ± 0.441</td><td>100%</td></tr><tr><td>Landmarks × temperature</td><td>0.203 ± 0.568</td><td>0.357</td><td>0.721</td><td>0.203 ± 0.568</td><td>18%</td></tr><tr><td>Dose × landmarks × temperature</td><td>-1.953 ± 1.171</td><td>-1.669</td><td>0.095</td><td>-1.953 ± 0.171</td><td>18%</td></tr></tbody></table>

LME, Linear Mixed Effect model; NEC, navigation environmental correlate. Effect estimates were computed after explanatory variables were centred and standardized to a mean of 0 and a s.d. of 0.5 to allow for direct comparison among effects. Results are given both for the complete LME statistical model and for a multimodel inference procedure²⁴ to account for a greater range of alternative, more parsimonious combinations of explanatory variables with respect to their relative weight of statistical support. The relative importance measures each variable's occurrence frequency within the best candidate models, weighted by the models' respective statistical support. A relative importance of 100% indicates that the variable appears in each of the best models, and therefore receives maximal support as a potential correlate of homing failure. Bold entries indicate significant effects at a 0.05 level.

<figure><img src="image_39.png" /><figcaption>**Figure 2 | Observed cumulative homing probability as a function of time lapse since honeybee release at 1 km from the colony.** The interaction effects between treatment (control versus all doses confounded) and navigation environmental correlates (NECs) are depicted by successively splitting data based on (a) average temperature (28 °C) and (b and c) average landmark availability for orientation (4.2-km network of hedgerows and forest edges in the 1-km diameter area between the release point and the colony). Release time, that is, the origin of the temporal axis, was standardized among trials to 1630–1830 hours. Horizontal arrows indicate night time.</figcaption></figure>

1 km away and monitored using RFID (Radio-Frequency IDentificationentification) microchips. The median probability of HF ranged from 0.17 in control trials to 0.40 in higher exposure (> 0.5 ng per bee) trials. Controlling for the random effects of colonies and release points, the risk of HF, expressed as a 0-1 binary response variable, significantly increased with exposure dose in nanogram per bee (Linear Mixed Effect model (LME), n = 978, Z = 3.374, P < 0.001), following the inverse logit function in equation 1.

$$HF = 1/(1 + e^{-(0.48 \times \text{dose} - 1.53)}) \quad (1)$$

The dose-response function (Fig. 1a) returns a 46% increase of HF at a 1-ng dose (HF = 0.259) compared with control (0.177).

This effect size could be detected with a fairly high confidence level given our sample size (a *posteriori* statistical power = 99%, Fig. 1b). Yet, a substantial variability was observed among trials, with high HF values in some control treatments, as well as low HF values at higher doses. Such high level of variability might be partly accounted for by NECs, as detailed below.

**The dose-response function is context dependent.** As predicted, NECs interacted significantly with the thiamethoxam exposure (Table 1, Fig. 2), therefore altering the shape and steepness of the generic dose-response function of HF (Fig. 3). We focused our analyses on ambient temperature and availability of landmarks for orientation (total length of hedgerows and forest edges in the
1-km diameter area between the release point and the colony). Our measurements of solar radiations (W m⁻²), as an indication of sun compass availability, were strongly correlated with ambient temperature (r_s = 0.68, P < 0.001). When implemented into separate models, both temperature and solar radiations significantly influenced HF (LME; n = 978, Z = -5.55, P < 0.001

and n = 978, Z = -4.58, P < 0.001, respectively). But most importantly, they significantly interacted with thiamethoxam (n = 978, Z = -2.70, P = 0.007 and n = 978, Z = -2.55, P = 0.010, respectively). However, temperature enabled a better fit to HF data than did solar radiations (Akaike Information Criterion (AIC) of LMEs with dose interactions = 971.5 versus 983.6, respectively, with lower AIC indicating a better model fit), and was therefore considered as the most relevant weather-related NEC for investigating HF dose-response variations in our study.

The significant interactions (Table 1) are merely depicted by the dissociated curves of homing statistics as a function of time after release (Fig. 2). The *Control-Treatment* relative difference is negligible at higher temperatures and increases with decreasing temperature (Fig. 2a), as revealed by the negative *Dose × Temperature* interaction. At low temperatures (Fig. 2b,c), the homing gap between control and treated bees increases with landmarks availability. This is in accordance with the positive *Dose × Landmarks* interaction in model outputs (Table 1). Notwithstanding, higher landmarks availability (Fig. 2b) was associated with overall lower HF than in more open landscapes (Fig. 2c). All in all, temperature appeared as the proximal driver of HF dose-response variations because it had steeper interaction estimates and with a greater relative importance than landmarks availability (Table 1). However, landmarks availability was also an important environmental correlate of HF, judging from its pervasive occurrence in the set of best candidate models (Table 2), and therefore should be explicitly accounted for whenever possible as well. In this study, none of the best candidate models of HF reached a predominant weight of evidence (all wᵢ values ≪ 90%, Table 2), indicating that it is advisable to resort to multimodel averaging to derive more conservative model predictions²⁴.

**The ED addresses context dependency.** Multimodel conditional predictions of the dose-response function (Fig. 3) illustrate well the conclusion that NECs substantially modify the generic dose-response function of HF one may assess on the field. The dose-response curve is flattest at high temperatures in rather open landscapes (few landmarks available), and steepest at low

<figure><img src="image_40.png" /><figcaption>**Figure 3 | Conditional predictions of the thiamethoxam dose-response function of honeybee homing failure probability.** The predicted dose-response curves were computed for different combinations of temperature and landmarks availability, ranging from low landmark density (2 km long network of hedgerows and forest edges in the 1-km diameter area homing bees have to cover between the release point and the colony) to high landmark density (8 km). Intercepts show expected homing failure for non-treated bees in the context of the experiment. Intercepts therefore set the control value needed to calculate mortality rate due to homing failure, m<sub>HP</sub>, for honeybees exposed to a given dose.</figcaption></figure>

<table><caption>**Table 2 | Detailed set of all candidate explanatory LME models of homing failure in honeybees exposed to thiamethoxam, ranked by decreasing statistical support.**</caption><thead><tr><th>Model rank</th><th>Intercept</th><th>Dose</th><th>Landmarks</th><th>Temperature</th><th>Dose × landmarks</th><th>Dose × temperature</th><th>Landmarks × temperature</th><th>Dose × landmarks × temperature</th><th>k</th><th>logLik</th><th>AIC</th><th>ΔAIC</th><th>w<sub>i</sub></th></tr></thead><tbody><tr><td>1</td><td>-1.133</td><td>0.6580</td><td></td><td>-1.556</td><td></td><td>-0.9457</td><td></td><td></td><td>7</td><td>-477.952</td><td>969.9</td><td>0</td><td>0.278</td></tr><tr><td>2</td><td>-1.153</td><td>0.6639</td><td>-0.1413</td><td>-1.538</td><td>0.6763</td><td>-0.9876</td><td></td><td></td><td>9</td><td>-476.467</td><td>970.9</td><td>1.03</td><td>0.166</td></tr><tr><td>3</td><td>-1.198</td><td>0.6887</td><td>-0.2037</td><td>-1.486</td><td>1.0240</td><td>-1.1970</td><td>0.2027</td><td>-1.953</td><td>11</td><td>-474.737</td><td>971.5</td><td>1.57</td><td>0.127</td></tr><tr><td>4</td><td>-1.138</td><td>0.6558</td><td>-0.1460</td><td>-1.541</td><td></td><td>-0.9438</td><td></td><td></td><td>8</td><td>-477.786</td><td>971.6</td><td>1.67</td><td>0.121</td></tr><tr><td>5</td><td>-1.203</td><td>0.6673</td><td>-0.2818</td><td>-1.486</td><td>0.6821</td><td>-1.0260</td><td>0.5539</td><td></td><td>10</td><td>-476.106</td><td>972.2</td><td>2.31</td><td>0.088</td></tr><tr><td>6</td><td>-1.185</td><td>0.6682</td><td></td><td>-1.413</td><td></td><td></td><td></td><td></td><td>6</td><td>-480.375</td><td>972.7</td><td>2.85</td><td>0.067</td></tr><tr><td>7</td><td>-1.190</td><td>0.6598</td><td>-0.2868</td><td>-1.490</td><td></td><td>-0.9891</td><td>0.5792</td><td></td><td>9</td><td>-477.411</td><td>972.8</td><td>2.92</td><td>0.065</td></tr><tr><td>8</td><td>-1.206</td><td>0.6731</td><td>-0.1528</td><td>-1.391</td><td>0.6078</td><td></td><td></td><td></td><td>8</td><td>-479.108</td><td>974.2</td><td>4.31</td><td>0.032</td></tr><tr><td>9</td><td>-1.191</td><td>0.6665</td><td>-0.1442</td><td>-1.398</td><td></td><td></td><td></td><td></td><td>7</td><td>-480.203</td><td>974.4</td><td>4.50</td><td>0.029</td></tr><tr><td>10</td><td>-1.241</td><td>0.6743</td><td>-0.2481</td><td>-1.351</td><td>0.6160</td><td></td><td>0.3800</td><td></td><td>9</td><td>-478.918</td><td>975.8</td><td>5.93</td><td>0.014</td></tr><tr><td>11</td><td>-1.225</td><td>0.6680</td><td>-0.2340</td><td>-1.359</td><td></td><td></td><td>0.3772</td><td></td><td>8</td><td>-480.017</td><td>976.0</td><td>6.13</td><td>0.013</td></tr><tr><td>12</td><td>-1.159</td><td></td><td></td><td>-1.384</td><td></td><td></td><td></td><td></td><td>5</td><td>-486.442</td><td>982.9</td><td>12.98</td><td>&lt;0.001</td></tr><tr><td>13</td><td>-1.165</td><td></td><td>-0.1505</td><td>-1.369</td><td></td><td></td><td></td><td></td><td>6</td><td>-486.250</td><td>984.5</td><td>14.60</td><td>&lt;0.001</td></tr><tr><td>14</td><td>-1.197</td><td></td><td>-0.2380</td><td>-1.332</td><td></td><td></td><td>0.3586</td><td></td><td>7</td><td>-486.082</td><td>986.2</td><td>16.26</td><td>&lt;0.001</td></tr><tr><td>15</td><td>-1.232</td><td>0.6439</td><td></td><td></td><td></td><td></td><td></td><td></td><td>5</td><td>-497.628</td><td>1,005.3</td><td>35.35</td><td>&lt;0.001</td></tr><tr><td>16</td><td>-1.256</td><td>0.6498</td><td>-0.2266</td><td></td><td>0.6397</td><td></td><td></td><td></td><td>7</td><td>-495.868</td><td>1,005.7</td><td>35.83</td><td>&lt;0.001</td></tr><tr><td>17</td><td>-1.239</td><td>0.6421</td><td>-0.2281</td><td></td><td></td><td></td><td></td><td></td><td>6</td><td>-497.099</td><td>1,006.2</td><td>36.29</td><td>&lt;0.001</td></tr><tr><td>18</td><td>-1.209</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>4</td><td>-503.355</td><td>1,014.7</td><td>44.81</td><td>&lt;0.001</td></tr><tr><td>19</td><td>-1.216</td><td></td><td>-0.2285</td><td></td><td></td><td></td><td></td><td></td><td>5</td><td>-502.810</td><td>1,015.6</td><td>45.72</td><td>&lt;0.001</td></tr></tbody></table>

AIC, Akaike Information Criterion; LME, Linear Mixed Effect model.

AIC is an inverse indicator of model parsimony²⁴ considering fit (logLik = log-Likelihood) and complexity (k = number of parameters to be estimated in the candidate model). Estimates are given for each explanatory variable after those were centered to mean = 0 and standardized to s.d. = 0.5. Models ignoring the dose effect all received a virtually negligible statistical support (AIC weight of evidence w; <0.001, that is, the probability of model *i* as being the best model in the set is nearly 0). The top-model set with a ΔAIC <2 (AIC difference with the best candidate model) comprises four concurrent models (ranks 1-4, in bold) with a weight of evidence w, ranging from about 12 to 28%. All environmental correlates of honeybee navigation and their interactions with the dose effect are retained in at least one of the top models. The ΔAIC <2 cutoff rule was used to define the top-model set²⁴.
temperatures in rather complex landscapes (high landmarks density). Those conditional dose-response functions (Fig. 3) should be viewed as tentative predictions specific to our experimental context, and be further refined using additional data at the edges of temperature and landmark ranges. Homing experiments carried out with foragers familiar with the release point may also provide lower bound estimates of HF¹.

Because the spatial context, as revealed by landmarks density, is an important correlate of HF, one can derive spatially explicit maps of the expected HF for a given exposure level. However, given that consensus rarely emerges about field exposure levels²⁵, one may instead map the exposure level required to reach a given effect size. A critical pesticide effect is reached when HF adds up to background mortality up to a point that entails a 35% reduction of colony size during an exposure event⁶, for example, during a 30-day oilseed rape flowering period. Using previously published colony dynamics parameters¹¹, we found this critical slowdown is triggered at a mortality rate because of HF $m_{HF} > 0.15$, that is, if non-returning treated bees represent more than 15% of returning bees, one might expect from the control. For the sake of illustration, we therefore targeted a critical $m_{HF}$ value of 20% and mapped the corresponding required thiamethoxam

dose, called ED₂₀ (Fig. 4). All landscape units (pixels of the map) within the focus landscape (450-km², Fig. 4) were assigned an ED₂₀ > 0.320 ng per bee at 20 °C. This is greater than worst-case thiamethoxam residue intakes by foragers for 1h flying (estimated to be <0.276 ng per bee for 20% sugar content nectar²⁶), meaning that the critical HF is unlikely to be reached over the focus landscape.

## Discussion

Our experiment demonstrates that HF in treated bees is highly context dependent, with temperature as a major proximal correlate of homing performance, and landmarks availability as a secondary, although interacting environmental correlate. Previous literature review²⁷ have already shown that some pesticides (pyrethroids) have greater toxicity at low temperatures (for example, below 28–30 °C), which could be related with impairments of neural mechanisms of thermogenesis rather than with troubles of flight muscular coordination or metabolic rate. Neonicotinoids may, however, involve different toxicological pathways. In particular, we cannot exclude that neonicotinoids interfere with the honeybee orientation system.

<figure><img src="image_41.png" /><figcaption>**Figure 4** | Spatial assessment of the honeybee experimental homing failure, as delineated by the ED₂₀ (effective dose that leads to a 0.20 mortality rate due to homing failure). Mortality due to homing failure, $m_{HF}$, is given by the relative difference of control versus treated homing probabilities. Conditional predictions show the sharp transition from (a) low-effect conditions at 28 °C to (b) high-effect conditions at 20 °C. Thick lines show the 1-ng ED₂₀ limits, delineating the landscape contexts where honeybees exposed to a 1-ng dose will be subject to a mortality $m_{HF} > 0.20$. This area covers 34–75% of the whole landscape depending on the temperature (c), but may be much reduced if one focuses on lower exposure scenarios. Eventually, for doses <0.320 ng per bee at 20 °C, the landscape is expected to be exempt of critical homing failure risks (ED₂₀ landscape coverage = 0%). For the sake of comparison, horizontal lines indicate worst-case scenarios of thiamethoxam intake by honeybees for 1h flying, or to forage a day, based on 20–40% sugar nectars²⁶. In particular, the *per-hour* exposure scenarios (for example, 0.184–0.276 ng per bee for a 20% oilseed rape nectar) fall in the 0% range for ED₂₀ landscape coverage. (d) ED₂₀ conditional predictions use the local geographic information on the network of hedgerows and forest edges.</figcaption></figure>
<table><caption>**Table 3 | Sample size allocation among trials, colonies and experimental doses.**</caption><thead><tr><th></th><th>Control</th><th>0.42 ng</th><th>0.85 ng</th><th>1.19 ng</th><th>1.34 ng</th><th>2.39 ng</th></tr></thead><tbody><tr><td>Total number of honeybees</td><td>370</td><td>186</td><td>180</td><td>62</td><td>118</td><td>62</td></tr><tr><td>Total number of trials</td><td>22</td><td>12</td><td>12</td><td>4</td><td>6</td><td>4</td></tr><tr><td>Trials with >15 honeybees</td><td>18</td><td>5</td><td>7</td><td>1</td><td>6</td><td>3</td></tr><tr><td>Number of experimental colonies</td><td>9</td><td>6</td><td>6</td><td>2</td><td>1</td><td>2</td></tr></tbody></table>

A trial refers to a group of either control or treated bees, caught at one of the nine experimental colonies and released at one particular point 1km away. Releases were further split between two consecutive days to cover a wider range of weather conditions. Average conditions (± s.d.) for released bees were 28.6 ± 4.7 °C ambient temperature, 437.2 ± 132.3 W m⁻² solar radiations and 4.2 ± 2.6 km landmark availability. One or two treated trials were carried out at a time, and systematically coupled with a simultaneous control trial from the same colony.

<figure><img src="image_42.png" /><figcaption>**Figure 5 | Map of the study area and location of experimental colonies.** (a) Location of the study area in western France. (b) Location of the experimental colonies (red) and associated release points (blue) superimposed on the woody vegetation layer (Institut Géographique National) illustrating regional variations in landscape complexity. The solid black line delineates regional boundaries. (c) Details of one experimental site, emphasizing the contrasted landscape context between release points.</figcaption></figure>

Interestingly, the greatest treatment effect was measured during poor weather conditions in environments with a greater density of landmarks (Fig. 2b), suggesting that treated foragers had greater difficulties to recover and integrate the information on vector flights among learned landmarks and/or to derive novel shortcut vector flights¹⁸,¹⁹. This also matches the recent evidence from harmonic radar tracking surveys that honeybees experience difficulties in their homing flight pattern after non-lethal exposure to neonicotinoid pesticides²⁸. When forced to split their route flight back to the hive into a vector phase based on sun compass and a homing search phase from a known landmark, exposed honeybees were less efficient in completing the second phase, which involves retrieval of remote memory²⁸. This conforms well to our observation that higher landmark availability is associated with lower ED levels to reach the critical 20% HF (Fig. 4).

The spatially explicit predictions of ED₂₀ delineate consistent high-effect areas at medium temperatures (Fig. 4a). As temperature decreases, the high-effect area expands over most of the landscape, restricting low-risk areas to a few remote spots (Fig. 4b). However, further research is necessary to determine whether those low-effect refuges may be viewed as effect free for apiaries, or not. Indeed, those areas correspond to the most open landscapes (Fig. 4d) where both control and treated HF are high (Fig. 2c). They may thus rather be viewed as environments with greater challenge imposed upon honeybee navigation during poor weather conditions, where any treatment effect will be harder to detect. Notwithstanding, with the help of experts' assessment of the field exposure levels²⁹,³⁰ and how they compare with ED₂₀ values, this approach allows determining the proportion of the landscape at risk of detrimental HF in case of exposure.

We found the worst-case daily cumulative thiamethoxam residue intake (for example, 0.736–2.944 ng per bee for 20% sugar content nectar²⁶) entails a 70–80% landscape ED₂₀ coverage (Fig. 4c). However, the residue intake cumulated over an entire foraging day is not perceived as a valid, field realistic, exposure scenario for free-ranging bees, particularly in the case of transient behavioural observations such as homing tests⁶,²⁶. Alternatively, one should focus on an hourly basis, best suited for studying behavioural issues at the scale of foraging bouts. Here, the worst-case scenarios of residue intakes by foragers for 1h flying (estimated to be <0.276 ng per bee for 20% sugar content nectar²⁶) are unlikely to reach critical ED₂₀ levels at 20 °C in any landscape unit of the focus area (all ED₂₀ > 0.320 ng, Fig. 4c). However, given the rather narrow gap between those two end points, it might be critical to refine estimations of field dietary residue intakes and associated uncertainty levels, as well as the occurrence likelihood of worst-case scenarios in field conditions.

In conclusion, the assessment of sublethal effects in free-ranging honeybees may be highly context dependent. The strong context dependency helps understand the difficulty to establish consistent links between laboratory-based toxicological assays and field observations²⁹,³¹. This compromises the applicability of HF as a paradigmatic approach to the sublethal risk assessment in field conditions because neonicotinoids will not return repeatable, generic, dose-response functions. Instead, we argue that field assessment schemes should either (i) provide details about landscape and weather conditions, (ii) use thoroughly standardized landscape and weather conditions or (iii) cover an extensive range of contexts along with substantial sample sizes. In addition, large-scale epidemiological surveys aimed at exploring potential links between colony losses and phytosanitary practices
across regions or countries should whenever possible incorporate interacting climate or landscape covariates. In that respect, the ED framework, although involving trials at a range of concentrations above the expected field exposure levels, has the potential to become an operative tool to help toxicologists identify some specific contexts that should deserve particular attention in the course of higher tier risk assessment.

## Methods

**Homing experiments.** This homing study covers ongoing experiments from springs 2011 (Experiment 2 in Henry *et al.*¹) and 2012, using RFID individual monitoring³² coupled with landscape-scale characterization of salient landmarks as orientation cues³³. We analysed the homing statistics of 978 individual honeybees obtained from nine colonies in the course of 60 distinct homing trials. A homing trial was defined as a group of about 15–20 bees (mean and s.d. = 16.3 ± 2.9) released at one given site after receiving an acute oral thiamethoxam or control treatment (Table 3). Five thiamethoxam doses were chosen so as to cover a broad range of LD₅₀ fractions (from LD₅₀/20 to LD₅₀/2, with LD₅₀ = 5 ng per bee) and to increase the ability of the study to delineate a consistent dose-response function. Real experimental doses were subsequently titrated at 0.42–2.39 ng per bee (LD50/12 to LD50/2). This range lies above the worst-case thiamethoxam residue intakes by foragers expected for an hour flying (for example, 0.184–0.276 ng per bee with an average 20% sugar content nectar of treated oilseed rape²⁶). However, this did not compromise the study, which was first and foremost designed to establish a dose-response profile and its context dependency. The dose was administered to bees in 20 µl of a 50% sucrose solution. Bees were kept for a 40-min assimilation period before release. Their homing flight was monitored by RFID microships glued on their thorax and a series of ad-hoc readers placed at the colony entrance¹,³².

**Field experimental design.** *A. mellifera* ligustica strain colonies were headed by sister queens bred the previous summer during the same grafting session and wintered in nuclei, then set up on experimental location in April. Colonies were disposed >3 km apart from each other (nearest distance range = 3.7–10.4 km, Fig. 5), across the ECOBEE platform, western France³⁴. Homing experiments were performed after oilseed rape flowering period to avoid interferences with field treatments. Twenty-two release sites, located 1 km away from colonies (two to six per colony, Fig. 5) were chosen so as to maximize the overall variability of their landmark density values (total length of hedgerows and forest edges included in the 1 km diameter circular area homing bees need to cross over to get back to their colony). Treated and control groups were released simultaneously at each trial, and at a standardized time among trials (16:30–18:30 h).

**NECs.** Experiments were exclusively carried out under weather compatible with current foraging activity, that is, excluding rainy, cold (<15 °C) and windy (>15 km h⁻¹) conditions³⁵. Ambient temperature was measured at the point and time of release. Solar radiations (W m⁻², averaged hourly) were measured using a pyrometer placed in the centre of the study area, as an indication of cloudiness, and recorded on a continuous basis. Landmarks availability was measured as the total length (km) of hedgerows and forest edges (*sensu*³³) in the 1-km diameter area homing bees have to cover between the release point and the colony.

**Data analysis.** HF was defined as an absence of RFID record during the 7 days post-release, and was coded as a 0/1 binary response variable (return/no return, respectively). In a first step, the binary data were analysed in relation with exposure doses using a LME to establish a generic dose-response function. To conform the assumption of independent outcomes in the binary homing observations, the identity of experimental colonies and of release points were included as random grouping variables, along with the dose effect as fixed explanatory variable. A statistical power analysis was performed *a posteriori* to assess the ability of the LME to detect a significant HF dose-response with the current sample size. Statistical power was measured as the proportion of significant effects out of 1,000 simulated HF data sets recomputed from random binomial distributions with probability constrained on a range of potential effect sizes. In a second step, NECs were added to the LME model as additional correlates of HF, interacting with the dose effect. All the possible two- and three-way interactions among explanatory variables were considered within the frame of a multimodel inference procedure²⁴. The most plausible combinations of explanatory variables were determined using the AIC of model parsimony, considering model fit and complexity. We used a ΔAIC<2 cutoff rule to define the most parsimonious models and determine the relative importance and significance of candidate exploratory variables. Variables were centered to mean = 0 and standardized to s.d. = 0.5 before the analysis, so that their respective effects are given in s.d. units and therefore are readily interpretable in terms of effect size and comparable among each others.

## References

1. Henry, M. *et al.* A common pesticide decreases foraging success and survival in honey bees. *Science* **336**, 348–350 (2012).

2. Whitehorn, P. R., O'Connor, S., Wackers, F. L. & Goulson, D. Neonicotinoid pesticide reduces bumble bee colony growth and queen production. *Science* **336**, 351–352 (2012).

3. Schneider, C. W., Tautz, J., Grünewald, B. & Fuchs, S. RFID tracking of sublethal effects of two neonicotinoid insecticides on the foraging behavior of *Apis mellifera*. *PLoS ONE* **7**, e30023 (2012).

4. Gill, R. J., Ramos-Rodriguez, O. & Raine, N. E. Combined pesticide exposure severely affects individual—and colony—level traits in bees. *Nature* **491**, 105–108 (2012).

5. Cressey, D. Europe debates risk to bees. *Nature* **496**, 408–408 (2013).

6. EFSA. Guidance on the risk assessment of plant protection products on bees (*Apis mellifera*, *Bombus* spp. and solitary bees). *EFSA J.* **11**, 3295 (2013).

7. Matsumoto, T. Reduction in homing flights in the honey bee *Apis mellifera* after a sublethal dose of neonicotinoid insecticides. *Bull. Insectol.* **66**, 1–9 (2013).

8. Yang, E. C., Chuang, Y. C., Chen, Y. L. & Chang, L. H. Abnormal foraging behavior induced by sublethal dosage of imidacloprid in the honey bee (Hymenoptera: Apidae). *J. Econ. Entomol.* **101**, 1743–1748 (2008).

9. Cresswell, J. A meta-analysis of experiments testing the effects of a neonicotinoid insecticide (imidacloprid) on honey bees. *Ecotoxicology*, **20**, 149–157 (2011).

10. Palmer, M. J. *et al.* Cholinergic pesticides cause mushroom body neuronal inactivation in honeybees. *Nat. Commun.* **4**, 1634 (2013).

11. Henry *et al.* Response to comment on 'A common pesticide decreases foraging success and survival in honey bees'. *Science* **337**, 1453–1453 (2012b).

12. Cresswell, J. E. & Thompson, H. M. Comment on 'A Common Pesticide Decreases Foraging Success and Survival in Honey Bees'. *Science* **337**, 1453–1453 (2012).

13. Chittka, L., Geiger, K. & Kunze, J. The influences of landmarks on distance estimation of honey bees. *Anim. Behav.* **50**, 23–32 (1995).

14. Chittka, L., Kunze, J., Shipman, C. & Buchmann, S. L. The significance of landmarks for path integration in homing honeybee foragers. *Naturwissenschaften* **82**, 341–343 (1995).

15. Crailsheim, K., Riessberger, U., Blaschon, B., Nowogrodzki, R. & HRassnigg, N. Short-term effects of simulated bad weather conditions upon the behaviour of food-storer honey bees during day and night (*Apis mellifera carnaica* Pollmann). *Apidologie* **30**, 299–310 (1999).

16. Dyer, F. C. & Gould, J. L. Honey bee orientation: a backup system for cloudy days. *Science* **214**, 1041–1042 (1981).

17. HRassnigg, N. & Crailsheim, K. Metabolic rates and metabolic power of honey bees in tethered flight related to temperature and drag (Hymenoptera: Apidae). *Entomol. Gen.* **24**, 23–30 (1999).

18. Menzel, R. *et al.* A common frame of reference for learned and communicated vectors in honeybee navigation. *Curr. Biol.* **21**, 645–650 (2011).

19. Menzel, R. *et al.* Vector integration and novel shortcutting in honeybee navigation. *Apidologie* **43**, 229–243 (2012).

20. Riessberger, U. & Crailsheim, K. Short-term effect of different weather conditions upon the behaviour of forager and nurse honey bees (*Apis mellifera carnaica* Pollmann). *Apidologie* **28**, 411–426 (1997).

21. Towne, W. F. & Moscrip, H. The connection between landscapes and the solar ephemeris in honeybees. *J. Exp. Biol.* **211**, 3729–3736 (2008).

22. Woods, W. A., Heinrich, B. & Stevenson, R. D. Honeybee flight metabolic rate: does it depend upon air temperature? *J. Exp. Biol.* **208**, 1161–1173 (2005).

23. EFSA. Conclusion on the peer review of the pesticide risk assessment for bees for the active substance thiamethoxam. EFSA J. **11**, 3067 (2013).

24. Burnham, K. P. & Anderson, D. R. *Model Selection and Multimodel Inference: a Practical Information-Theoretic Approach* (Springer, 2002).

25. Osborne, J. L. Ecology: Bumblebees and pesticides. *Nature* **491**, 43–45 (2012).

26. EFSA. Statement on the findings in recent studies investigating sub-lethal effects in bees of some neonicotinoids in consideration of the uses currently authorised in Europe. EFSA J. **10**, 2752 (2012).

27. Belzunces, L. P., Tchamitchian, S. & Brunet, J.-L. Neural effects of insecticides in the honey bee. *Apidologie* **43**, 348–370 (2012).

28. Fischer, J. *et al.* Neonicotinoids interfere with specific components of navigation in Honeybees. *PLoS ONE* **9**, e91364 (2014).

29. Pilling, E., Campbell, P., Coulson, M., Ruddle, N. & Tornier, I. A four-year field program investigating long-term effects of repeated exposure of honey bee colonies to flowering crops treated with thiamethoxam. *PLoS ONE* **8**, e77193 (2013).

30. Krupke, C. H., Hunt, G. J., Eitzer, B. D., Andino, G. & Given, K. Multiple routes of pesticide exposure for honey bees living near agricultural fields. *PLoS ONE* **7**, e29268 (2012).

31. Cutler, G. C. & Scott-Dupree, C. D. Exposure to clothianidin seed-treated canola has no long-term impact on honey bees. *J. Econ. Entomol.* **100**, 765–772 (2007).
32. Decourtye, A. et al. Honeybee tracking with microchips: a new methodology to measure the effects of pesticides. *Ecotoxicology* **20**, 429–437 (2011).

33. Henry, M. et al. Spatial autocorrelation in honeybee foraging activity reveals optimal focus scale for predicting agro-environmental scheme efficiency. *Ecol. Model*, **225**, 103–114 (2012).

34. Odoux, J.-F. et al. ECOBEE: a tool for long-term bee colony monitoring at landscape scale in West European intensive agrosystems. *J. Apicult. Res.* **53**, 57–66 (2014).

35. Rollin, O. et al. Differences of floral resource use between honey bees and wild bees in an intensive farming system. *Agric. Ecosyst. Environ.* **179**, 78–86 (2013).

## Acknowledgements

This study was funded by the European Community program (797/2004) for French beekeeping coordinated by the French Ministry of Agriculture (TECHBEE and RISQAPI projects). Special thanks also go to J. Aptel, M. Béguin, D. Fortini, C. Le Mogne, L. Miot, O. Rollin, S. Tchamitchian, L. Tison, C. Toullet, B. Vaissière for valuable suggestions and help on the field, as well as O. Marloie for providing the pyrometer for solar radiation measurement. Experimental solutions were dosed by the GIRPA (Angers, France). RFID

devices were designed by Tag Tracing Solutions, Valence, France, and adapted to honey bees with the help of J. Fourrier, F. Brun, J. Devillers, M. Gauthier, the ECOBEE and LTER 'Zone Atelier Plaine & Val de Sèvre' monitoring facilities.

## Author contributions

M.H., C.B., V.L.F., V.B. and A.D. designed the experiments and data analyses; M.H., C.B., F.R., J.-F.O. and P.A. carried out the fieldwork; all authors contributed to the discussion and manuscript writing.

## Additional information

**Competing financial interests:** The authors declare no competing financial interests.

**Reprints and permission** information is available online at http://npg.nature.com/reprintsandpermissions/

**How to cite this article:** Henry, M. et al. Pesticide risk assessment in free-ranging bees is weather and landscape dependent. *Nat. Commun.* 5:4359 doi: 10.1038/ncomms5359 (2014).
PROCEEDINGS B

rspb.royalsocietypublishing.org

Research

<figure><img src="image_43.png" /></figure>

CrossMark

*click for updates*

**Cite this article:** Henry M, Cerrutti N, Aupinel P, Decourtye A, Gayrard M, Odoux J-F, Pissard A, Rüger C, Bretagnolle V. 2015 Reconciling laboratory and field assessments of neonicotinoid toxicity to honeybees.

*Proc. R. Soc. B* **282**: 20152110.

http://dx.doi.org/10.1098/rspb.2015.2110

Received: 2 September 2015

Accepted: 19 October 2015

**Subject Areas:**

environmental science, ecology

**Keywords:**

*Apis mellifera*, imidacloprid, oilseed rape, pesticides, thiamethoxam

**Author for correspondence:**

Mickaël Henry

e-mail: mickael.henry@paca.inra.fr

Electronic supplementary material is available at http://dx.doi.org/10.1098/rspb.2015.2110 or via http://rspb.royalsocietypublishing.org.

Reconciling laboratory and field assessments of neonicotinoid toxicity to honeybees

Mickaël Henry<sup>1,2,†</sup>, Nicolas Cerrutti<sup>2,3,†</sup>, Pierrick Aupinel<sup>4</sup>, Axel Decourtye<sup>2,5,6</sup>, Mélanie Gayrard<sup>3</sup>, Jean-François Odoux<sup>4</sup>, Aurélien Pissard<sup>5</sup>, Charlotte Rüger<sup>3</sup> and Vincent Bretagnolle<sup>7,8</sup>

European governments have banned the use of three common neonicotinoid pesticides due to insufficiently identified risks to bees. This policy decision is controversial given the absence of clear consistency between toxicity assessments of those substances in the laboratory and in the field. Although laboratory trials report deleterious effects in honeybees at trace levels, field surveys reveal no decrease in the performance of honeybee colonies in the vicinity of treated fields. Here we provide the missing link, showing that individual honeybees near thiamethoxam-treated fields do indeed disappear at a faster rate, but the impact of this is buffered by the colonies' demographic regulation response. Although we could ascertain the exposure pathway of thiamethoxam residues from treated flowers to honeybee dietary nectar, we uncovered an unexpected pervasive co-occurrence of similar concentrations of imidacloprid, another neonicotinoid normally restricted to non-entomophilous crops in the study country. Thus, its origin and transfer pathways through the succession of annual crops need be elucidated to conveniently appraise the risks of combined neonicotinoid exposures. This study reconciles the conflicting laboratory and field toxicity assessments of neonicotinoids on honeybees and further highlights the difficulty in actually detecting non-intentional effects on the field through conventional risk assessment methods.

**1. Introduction**

In the current context of global honeybee decline, much attention has been paid to evaluating the possible contribution of neonicotinoid insecticides to colony weakening and collapse [1–3]. These systemic insecticides, which now represent ca 30% of insecticide use worldwide [4], pose a particular risk for pollinators, because once the active substance has been taken up in the plant, its residues translocate to the pollen and nectar collected by foragers throughout flowering. However, after 15 years of active research on the side effects of neonicotinoids on bees, a gap has emerged between the results of toxicity assessments in the laboratory and in the field [1–3]. Artificial exposure in laboratory experiments typically consists of providing individuals with contaminated food and comparing relevant physiological and behavioural endpoints with healthy control groups [5]. This approach has led to the identification of a range of sublethal effects, i.e. adverse physiological or behavioural

<footnote>†These authors contributed equally to this study.</footnote>

<footnote n="1"><sup>1</sup>INRA, UR406 Abeilles et Environnement, 84914 Avignon, France</footnote>

<footnote n="2"><sup>2</sup>UMT Protection des Abeilles dans l'Environnement, Site Agroparc, 84914 Avignon, France</footnote>

<footnote n="3"><sup>3</sup>Terres Inovia, Centre de Grignon, Avenue Lucien Brétignières, 78850 Thiverval Grignon, France</footnote>

<footnote n="4"><sup>4</sup>INRA, UE1255, UE Entomologie, 17700 Surgères, France</footnote>

<footnote n="5"><sup>5</sup>Association de Coordination Technique Agricole, Site Agroparc, 84914 Avignon, France</footnote>

<footnote n="6"><sup>6</sup>ITSAP - Institut de l'Abeille, Site Agroparc, 84914 Avignon, France</footnote>

<footnote n="7"><sup>7</sup>Centre d'Etudes Biologiques de Chizé, UMR 7372, CNRS and Université de La Rochelle, 79360 Beauvoir-sur-Niort, France</footnote>

<footnote n="8"><sup>8</sup>LTER 'Zone Atelier Plaine and Val de Sèvre', Centre d'Etudes Biologiques de Chizé, CNRS, 79360 Villiers-en-Bois, France</footnote>
<figure><img src="image_44.png" /></figure>

**Figure 1.** Field exposure experimental design. The maps show the location of the experimental fields with oilseed rape grown from thiamethoxam-treated seeds (23 fields (total 153 ha) in 2013 and 18 fields (135 ha) in 2014). Dots show the positions of beehives fitted with RFID readers, spatially allocated in a way that covers a broad range of field exposure levels. Thiamethoxam field exposure was computed as the sum of all treated surfaces in the territory, with individual field surfaces moderated by an ordinary IDW interpolation. Surfaces of treated fields located farther away than the average honeybee foraging range (1 km) were down-weighted by a $1/d^2$ multiplicative coefficient, where *d* is the distance (km) of the field to the colony. The resulting field exposure values ranged from 1 to 63 (mean = 15.7 ± 16.8 (s.d.)). The most exposed colony (exposure value = 82) showed dramatically high mortality rates compared to the other colonies, but concomitantly developed foulbrood syndromes. It was therefore discarded from the analyses to avoid overstating the excess mortality due to field exposure.

changes measurable at levels well below the 50% lethal dose (LD₅₀) currently used as a legal reference threshold for plant protection product approval. Neonicotinoids have been found to have sublethal effects on mobility, orientation, foraging and learning performances [6–11]. However, it is not clear whether these endpoints are reflected in a corresponding effect at field level [12]. The various attempts made so far to appraise the possible consequences of honeybee exposure to treated fields under normal agricultural practices—hereafter called *field exposure*—have failed to detect any noteworthy change in colony performance [13–16]. Without formal identification of deleterious effects in real-field exposure conditions, policy decisions will remain controversial [1–3].

Various explanations have been put forward for the empirical mismatch between laboratory and field experiments. One may argue that laboratory experiments have been conducted at exposure levels much higher than would normally occur in the field [17], or that they have overstated the insecticide effects owing to acute, rather than chronic, exposure. Alternatively, effects in the field may be compensated for by the resilience of honeybee colonies owing to demographic regulation mechanisms and honey storage [1–3]. Exposure effects may also occur after a time lag [18], with colony weakening becoming apparent only later in the season or with lower colony survival the next winter. Space should be taken into account as well as time. For instance, homing failure after artificial exposure has been detected in free ranging bees at the foraging range scale (i.e. 1 km away from the colony [9,10]), whereas field exposure experiments have typically surveyed colonies located in the immediate vicinity of, or even inside, treated crops. Such uncertainty

requires a spatially explicit approach to reconcile laboratory and field exposure experiments.

Here, we show that a landscape-scale increase in exposure to treated oilseed rape fields does indeed entail higher individual mortality. This study follows the recommendations of ANSES, the French food safety agency, to reassess the possible side effects of thiamethoxam [9] under real agricultural usage conditions. We obtained authorization from the French ministry of agriculture to grow winter oilseed rape from seeds treated with thiamethoxam (Cruiser OSR® formulation), for experimental purposes. Thiamethoxam is currently prohibited at the national and European levels [19]. With the permission of the land-owning farmers, Cruiser rapeseed was sown in 2 consecutive years on a total of 288 ha (153 ha in 2013 and 135 ha in 2014, figure 1) within the LTER Zone Atelier Plaine & Val de Sèvre area, France. Using RFID technology (radio frequency identification [20,21]), we monitored the life histories of 6847 individual bees in relation to levels of thiamethoxam exposure from the oilseed rape fields, and questioned whether individual survival as well as colony dynamics would vary with field exposure.

## 2. Material and methods

### (a) Study design

The study was initially designed to produce a gradient of real-field exposure to oilseed rape grown from seeds treated with thiamethoxam (figure 1). However, an unexpected concomitant exposure to imidacloprid, another neonicotinoid insecticide, was detected at substantial levels both in the nectar of experimental oilseed rape treated with thiamethoxam, and in the dietary nectar
ingested by foragers (see Results). Therefore, the studied field exposure level referred to in this study actually represents a gradient of combined exposure to both neonicotinoid products.

The field exposure gradient was experimentally achieved by sowing oilseed rape seeds coated with the Cruiser® formulation (thiamethoxam content 280 g l⁻¹) in a total of 288 ha (153 ha in 2013 and 135 ha in 2014) in a portion of the study area and positioning colonies around at various distances and directions to cover a range of exposure levels to treated fields. The study involved local volunteer farmers and required a special derogation from the French ministry of agriculture due to the current prohibition of any neonicotinoid treatment on oilseed rape.

Eighteen standardized experimental colonies bred from sister queens were set up in 10-frame Dadant hives and monitored in the field for colony dynamics and state of health [22]. Hives were fitted with RFID readers ([21]; electronic supplementary material, figure S1) so as to monitor the life history of a total of 46 cohorts of 100–250 honeybees during the oilseed rape flowering period (one to two cohorts of just-emerged bees per colony and one cohort of adult foragers per colony). We noticed the most exposed colony (field exposure = 82 units, see below) became largely depopulated and showed dramatically high cohort mortality rates unequalled by the other colonies. However, it concomitantly developed foulbrood syndromes. Therefore, we discarded this colony from analyses to avoid any risk of overstating the excess mortality due to field exposure *per se*. Analyses comprise the lifelong monitoring of 6847 individual bees from 17 colonies (1638 bees from 17 cohorts tagged at the foraging stage and 5209 bees from 27 cohorts tagged just after emergence).

## (b) Computation and validation of the field exposure level

Colony field exposure was calculated as the sum of all treated field areas with an ordinary inverse distance weighted (IDW) interpolation (figure 1). Surfaces of treated fields located farther away than the average honeybee foraging range (1 km) were down-weighted by a $1/d^2$ multiplicative coefficient, where *d* is the distance (km) of the field to the focus colony. One field exposure unit is therefore virtually equivalent to 1 ha of treated oilseed rape within a 1-km distance from the colony, or e.g. 4 ha at a 2-km distance.

The resulting field exposure values ranged from 1 to 63 (mean = 15.7 ± 16.8 (s.d.)) for the monitored colonies. Special attention was paid to reducing variability due to landscape composition or configuration. The experiment was carried out in the western part of the study area (ca 150 km²) characterized by an open field landscape with few semi-natural elements and with homogeneously scattered oilseed rape fields typically accounting for 8–10% of total land cover. Landscape complexity, referring to the amount of non-cropped interstitial habitats around fields (hedgerows and forest edges total length [10]) averaged 6.6 ± 2.8 km km⁻² in the vicinity of the colonies. Total oilseed rape (including treated experimental fields) land cover within a 1-km radius averaged 29.2 ± 19.7 ha. We further ensured that field exposure variations were independent of both total oilseed rape land cover and landscape complexity (Pearson's product-moment correlation, *r* = 0.32, *p* = 0.19 and *r* = -0.04, *p* = 0.86, respectively).

The experimental field exposure gradient was validated by neonicotinoid multi-residual analysis of dietary nectar collected from foragers at each colony entrance (see Results). On three occasions with one-week intervals during oilseed rape blooming, 200 returning honeybees were collected when entering their hive and narcotized with ether in a cage. One by one, their abdomens were gently pressed until their crop nectar content was regurgitated. Nectar samples were then pooled in eppendorf tubes to

get an average value of neonicotinoid residual content per colony and sampling date. The correlative link between field exposure level and neonicotinoid dietary residues was statistically validated by zero-inflated generalized linear mixed models (ZI-GLMM) because an excess of zeroes was detected in the data distribution, thus simple GLMMs did not fulfil model residual normality and homoscedasticity requirements.

## (c) Colony dynamics and radio frequency identification individual monitoring

Colony monitoring comprised systematic inspections for diseases as well as measurements of adult population size, honey reserves and female (worker) and drone (male dispersers) brood production [22]. Measurements were taken at the onset of flowering (18 April 2013 and 25 March 2014), and during the first and the fourth week after the end of flowering. Initial colony state was independent of the field exposure level (Pearson's product-moment correlations, population size: *r* = -0.15, *p* = 0.53 honey reserves: *r* = -0.41, *p* = 0.090; total brood surface: *r* = 0.21, *p* = 0.40; worker brood surface: *r* = 0.22, *p* = 0.37; drone brood: *r* = -0.080, *p* = 0.742). Mean weekly changes of colony parameters were computed for the five-week period referred to as *during* flowering and the three-week period referred to as *after* flowering. Colony changes were analysed against field exposure level using linear models (LM).

The RFID individual lifelong monitoring includes 6847 individual bees (5209 bees tagged just after emergence and 1638 tagged at the foraging stage) assigned to a total of 44 cohorts belonging to 17 colonies. The RFID technology was used to assess three key individual life-history parameters in the context of field toxicology, namely mortality rate, frequency of flight activity and precocious behavioural maturation—with the precocious onset of foraging suspected as a possible mechanism for compensating forager excess mortality [23–25]. The fully detailed information about the RFID system provider, RFID tag characteristics and fixation on honeybees may be found in previously published material [9,10,21]. The cohorts of newly emerged worker bees with homogeneous age were obtained by caging a queen with preselected brood frames containing enough empty cells for eggs [26]. Therefore, each year, newly emerged cohorts were obtained from a single colony independent of the experimental ones, but still originating from sister queens within the same livestock, held at the UE Entomologie of INRA Magneraud. Newly emerged cohorts were tagged and released in colonies during the week preceding the onset of flowering, and with a 5-day interval whenever two distinct cohorts were introduced in the same colony. Returning adult foragers were captured on their way back at the entrance of colonies, RFID-tagged and reintroduced into colonies. The forager tagging sessions occurred ca 10–12 days after tagging of the just-emerged bees, i.e. the approximate average duration needed for newly emerged cohorts to perform first flights and become foragers themselves [27–29]. In doing so, newly emerged cohorts and forager cohorts could be simultaneously monitored during their foraging life stage and during the oilseed rape full blooming period. RFID readers record the tag signal of bees passing through the hive entrance and store detection events along with the date and time (±0.01 s). Detections were pooled to the nearest minute, assumed to be the minimal time required for a bee to perform a flight. Therefore, occasional multiple detections of the same individual in less than a minute were considered to be part of the same event.

The hypotheses of (i) increased mortality rate, (ii) precocious behavioural maturation, and (iii) decreased flight activity were investigated using, respectively, the date of the last RFID detection event, the date of the first event (newly emerged cohorts only) and the number of events per day. Tagged honeybees were recorded as alive until the day following last tag detection event. Cox proportional hazard survival models (Cox PH [30]) were then
applied on a daily basis to test whether field exposure was associ-
ated with an excess mortality compared to baseline mortality. Cox
PH models are semi-parametric analyses specifically designed to
test the effect of covariates on the time lapse before occurrence of
an event. We applied the Cox PH formula to compare time to last
event (increased mortality rate hypothesis) and time to first event
(precocious behavioural maturation hypothesis). The non-indepen-
dency of individuals from the same colony was accounted for by
specifying colonies' identity as a grouping cluster. Whenever ana-
lysed jointly in the same Cox PH model, cohort types (i.e. foragers
versus just-emerged bees) were specified as distinct strata due
to their different baseline mortality. Finally, variations of flight
activity (number of events recorded per day, normalized with a
log(value + 1) transformation) were investigated using GLMMs
with colony identity specified as a random grouping variable.
Random slopes as a function of date and maximal daily temperature
were also allowed to account for the temporal non-independency of
repeated measurements as well as the possibly contrasted weather
conditions throughout the study. Rainy days (more than 1 mm)
were discarded from the flight activity analyses. For all studied
life-history traits, inter-annual variations were also explicitly
tested, and whenever significant, possible latent interactions
with the field exposure level were also inspected. Analyses were
performed with the R software for statistical computing [31].

(d) Imidacloprid nectar contamination

To investigate the source of imidacloprid contamination in the
environment, we also performed neonicotinoid multi-residual
dosages on the treated oilseed rape nectar collected with
microcapillaries directly from flowers in a subset of nine thia-
methoxam experimental fields. Imidacloprid residuals were
detected at similar concentrations and frequency in flower
nectar samples (0.1-0.9 ppb in six out of nine samples) and
honeybee dietary nectar (0.1-1.0 ppb in 13 out of 17 colonies).

Given the substantial levels of imidacloprid residuals unex-
pectedly found in both oilseed rape flower nectar and in
honeybee dietary nectar, we extended the flower nectar survey
to as many non-experimental oilseed rape fields as possible
within the 450 km² study area. A total of 73 oilseed rape fields
were thus sampled for nectar between 15 and 24 April 2014.
Nectar was collected between 9.00 h and 19.00 h by gently insert-
ing a 5 µl glass microcapillary tube (Drummond Scientific,
Broomall, PA, USA) into randomly chosen open flowers [32]
until achieving a cumulative volume of 25 µl per field. Nectar
samples were collected beyond a 10 m buffer distance from the
field margin to avoid any edge effect.

Both flower nectar samples collected from oilseed rape fields
and dietary nectar samples collected from forager stomachs
were sent to the European Union reference laboratory for neo-
nicotinoid multi-residual analyses (ANSES, Sophia-Antipolis,
France). Residues were quantified (limit of detection = 0.1 ppb,
limit of quantification = 0.3 ppb) by liquid chromatography
with electrospray tandem mass spectrometry [33].

**3. Results and discussion**

In this study we found that individual honeybees near
thiamethoxam-treated fields do indeed disappear at a
faster rate, but the impact of this is buffered by the colonies'
demographic regulation response.

(a) Validation of the field exposure design

Thiamethoxam residues found in dietary nectar brought
back to colonies by foragers increased significantly with
experimental thiamethoxam field exposure, validating the

field exposure design (electronic supplementary material,
figure S2A). Thiamethoxam residues remained undetected
in dietary nectar in colonies with limited field exposure
(less than eight field exposure units) and ranged from 0.1
to 0.8 ppb in colonies with the highest field exposure
(8–63 units). Furthermore, when dietary nectar thia-
methoxam is regressed against field exposure, the intercept
is not significantly different from zero (electronic supplemen-
tary material, figure S2A), confirming that the least exposed
locations may be viewed as thiamethoxam-free environments
for a sound basis of comparison. However, we discovered an
unexpected and substantial concomitant exposure to imida-
cloprid (electronic supplementary material, figure S2) in the
dietary nectar samples (0.1–1.0 ppb in 13 out of the 17 sur-
veyed colonies). Those residual levels are high enough to
potentially exert side effects on bees [34,35]. In France, this
neonicotinoid insecticide is currently used as a seed dressing
treatment (Gaucho® formulation) for a range of non-entomo-
philous crops such as wheat, barley and sugar beet, but has
been prohibited for sunflower since 1999 and has never
been used on oilseed rape [36]. Imidacloprid residues cor-
related with our experimental field exposure level and
*a fortiori* with thiamethoxam residues (electronic supple-
mentary material, figure S2B and C). This unexpected
contamination forced us to reconsider the experimental thia-
methoxam field exposure as a combined (thiamethoxam and
imidacloprid) neonicotinoid field exposure.

(b) Honeybee survival and life histories in relation to
field exposure

Honeybees disappeared at a faster rate with increasing field
exposure and this excess mortality increased over time (elec-
tronic supplementary material, table S1). The baseline
mortality increased by 10.1% (CI95%= 3.0–17.7%) per
15 field exposure units (Cox PH survival analysis, n = 78 716
daily records for 6847 bees, z = 2.85, p = 0.004). However,
we detected a highly significant deviation from hazard pro-
portionality ($\chi^2$ = 258, p < 0.001), indicating that the excess
mortality was not stationary but actually increased with
time. We therefore reassessed the field exposure effect, expres-
sing it in interaction with time [30]; this indeed returned a
better model fit to the data (Cox PH, field exposure effect:
z = 3.52, p < 0.001). The level of excess mortality was now
raised by 5.6% per 15 field exposure units per additional
week, i.e. from an average 5.6% excess mortality at the onset
of flowering to 22.4% after three more weeks had passed.
This time-mediated survival pattern was consistent between
the two study years (Cox PH, inter-annual variations: z =
0.41, p = 0.68), and might be interpreted either as a conse-
quence of the accumulation of neonicotinoid residues in food
and hive materials over time, with a delayed effect [34], or as
the emergence of a chronic exposure effect [18]. The field
exposure effect appeared to be independent of the landscape
spatial context (electronic supplementary material, table S2),
contrary to what would have been expected from previous
studies [10,37].

To rule out possible confounding effects due to cohort
types, we recomputed separate non-stationary survival
models for foragers (n = 1638 bees from 20 cohorts) and
just-emerged bees (n = 5209 bees from 24 cohorts). We con-
firmed that field exposure accelerated the disappearance of
both forager and just-emerged bee cohorts (figure 2), but
<figure><img src="image_45.png" /></figure>

<figure><img src="image_46.png" /></figure>

**Figure 2.** Non-stationary excess mortality due to field exposure during oilseed rape flowering. The honeybee excess mortality level due to field exposure is not stationary, but steadily increases in time. (a) Non-stationary excess mortality due to field exposure in 1638 bees tagged at the foraging stage with RFID microchips. (b) Non-stationary excess mortality due to field exposure in 5209 bees tagged after emergence with RFID microchips. Shaded areas show the 95% confidence envelope of the Cox PH estimate of excess mortality, expressed in per cent of the baseline mortality, and indicating a significant field exposure effect whenever it is above zero. Excess mortality was sequentially reassessed along the temporal axis using a left-censoring procedure, i.e. discarding newly disappeared individuals at each time step. For the sake of comparison, dashed lines show the field exposure effect that would be found on the assumption of stationary excess mortality.

the excess mortality was 2.4 times higher for the former group (11.9% versus 4.9% excess mortality per 15 field exposure units and per week; electronic supplementary material, table S1). As an illustration of the non-stationary hazard during oilseed rape flowering, the excess mortality linked with field exposure increased steadily throughout the survey for both cohort types (figure 2). The rise in mortality is quite clear (figure 3a,b) when average survival below the field exposure threshold that entails no dietary residues of thiamethoxam (field exposure < 8 units, *n* = 24 cohorts from nine colonies) is compared with average survival above that threshold (more than 8 units, *n* = 20 cohorts from eight colonies). An *a posteriori* power analysis (figure 3c) established that the field exposure effect as revealed by Cox PH survival analyses was strong enough to be satisfactorily detected (statistical power > 80% [38]) with a subset of 12 out of our 17 surveyed colonies.

Orientation disorders reported by artificial exposure experiments [9,39,40] may contribute to explaining our field exposure results. Conversely, a recent treated versus control field exposure RFID survey reported no effect of foraging on thiamethoxam-treated oilseed rape [41]. However, this study was carried out at a very different spatial scale and

<figure><img src="image_47.png" /></figure>

<figure><img src="image_48.png" /></figure>

<figure><img src="image_49.png" /></figure>

**Figure 3.** Mean cumulative disappearance of monitored individuals under high versus low field exposure. Mean proportion of (a) undetected foragers and (b) undetected just-emerged bees over time with high (8–63 units, black dots) and low (less than eight units, white dots) field exposure values. The high non-detection levels during the first week correspond to the in-hive life stage before first exits in young bees. (c) A posteriori power analysis showing the ability (%) of the Cox PH survival analysis to detect the significant excess mortality due to field exposure out of random subsets of the 17 experimental colonies (*n* = 100 subsets for dots up to 14 colonies, and *n* = 50 and 16 subsets for 15 and 16 colonies, respectively). Dashed line shows the recommended 80% statistical power threshold.

with no molecular validation of dietary nectar contamination, therefore precluding any meaningful comparison. In particular, treated and control treatments from that study would both fall in the very initial part of the field exposure range we actually covered here, with fairly low scores (ca 1–2
field exposure units, with a single 2-ha treated field and 1-km isolation distances [41]). This study encompasses up to 63 field exposure units, and reports evidence of thiamethoxam residues only beyond eight field exposure units.

An increased forager mortality rate is expected to trigger the precocious recruitment of younger, less efficient, foragers as an adaptive compensatory mechanism [23–25]. We, however, found no evidence that the excess mortality due to field exposure was further echoed by side effects on age at first exit (electronic supplementary material, table S3) or daily rate of flight activity (electronic supplementary material, table S4). Although no significant modification of age at first exit was detected, we noticed a significant non-stationarity, i.e. a significant variation over time, in the rate at which just-emerged monitored bees performed their first flights in relation to field exposure (model 1 in electronic supplementary material, table S3). Up to the 10th day of monitoring, first flights occurred 19.6% earlier for an average field exposure of 15 units. But up to the 20th day of monitoring, first flights occurred on average 8.8% later in the honeybee life for an average field exposure of 15 units, therefore offsetting the apparent first flight precocity. This non-stationarity of age at first flight along the field exposure gradient suggests that a portion of the monitored young bees indeed became precocious foragers while the others were compelled to spend a longer lifetime as in-hive worker to cope with the increased nursing tasks. Though rather speculative, this tentative scenario is consistent with the observations made at colony level on brood production dynamics (electronic supplementary material, table S5 and figure S3), as detailed below.

### (c) Colony dynamics in relation to field exposure

In spite of the excess mortality measured at the individual scale, highly exposed colonies did not show altered performance *per se* in terms of population growth and honey and brood production. However, there was a change in the way reproductive effort was allocated between female (worker) brood and drone (male disperser) brood. The field exposure did not significantly affect week-by-week changes in the adult population during or after oilseed rape flowering (models 1 and 7 in the electronic supplementary material, table S5), or weekly honey storage or weekly brood production (models 2, 3, 8 and 9 in the electronic supplementary material, table S5), at least within the effect size resolution limits permitted by the study sample size (statistical power analysis, smallest detectable effects estimated to ±21–58% of the mean; electronic supplementary material, table S5). But the field exposure did trigger a significant change in the relative proportions of worker brood versus drone brood production (electronic supplementary material, figure S3). During flowering, the most exposed colonies tended to invest more in worker brood production at the expense of drone brood production (electronic supplementary material, figure S3A). Drone brood development was delayed in exposed colonies; after flowering, drone brood production followed the field exposure gradient, being significantly higher in the more exposed hives (electronic supplementary material, figure S3B).

Rather than a decline in colony performance strictly speaking, these patterns should be viewed as a by-product of colonies' demographic compensation and regulation of reproductive investment [42]. Drones are more costly to produce

and maintain than workers, among others because they do not participate to the foraging task force. Thus, colonies decrease drone production when foraging conditions are poor, either due to resource scarcity [42] or seemingly due to forager excess mortality. Drone production typically peaks in spring [22,42] when virgin queens are most abundant and then when drones are most likely to successfully mate and pass along genes to other colonies. Delayed drone production might somehow disrupt this biological synchrony, and should therefore be addressed in terms of mating success or fitness value of reared drones. This should be further explored along with colony performance metrics using large colony monitoring datasets for enhanced statistical power.

### (d) Imidacloprid nectar contamination

Likewise, more detailed studies on the environmental fate of neonicotinoid residues are urgently needed to properly control for potential confounding effects or synergistic effects between different active substances. Indeed, the concomitant occurrence of imidacloprid residues makes it difficult to assign the excess mortality to thiamethoxam alone or to a combined effect, as has been revealed in various combined exposure trials [43,44]. An increasing number of studies [45–47] report substantial contamination of soil and puddles by neonicotinoid residues that may be subsequently taken up by the next crop in the succession. To investigate the source of imidacloprid contamination of the dietary nectar, we collected and analysed oilseed rape floral nectar samples from 82 oilseed rape fields in the study area in 2014 (nine of the treated experimental fields and 73 additional fields). Imidacloprid was undetected in the nectar of 30 sampled fields (36.6%), and varied from 0.1 to 1.6 ppb (median = 0.4 ppb) in the remaining 52 fields (63.4%). Those results concur with the substantial re-uptake of neonicotinoid residues recently reported in pollen and nectar samples from wild flowers in field margins [48], sometimes at even higher concentrations than in the flowering crop nectar itself. Various hypotheses may be proposed for those observations, including the persistence and accumulation of neonicotinoid residues in the soil throughout one or more annual crop succession cycles, their lateral movement and leaching in adjacent slopes, or even possible contaminations via the seed-coating machinery [48]. The precise pathway by which imidacloprid used on wheat or barley can transfer to oilseed rape nectar later on therefore requires urgent clarification because it is liable to compromise any effort scientists and risk assessors make to reconcile the findings of laboratory and field exposure surveys.

## 4. Conclusion

Overall, our results lead to two main conclusions. First, we found that field exposure to thiamethoxam combined with concomitant imidacloprid contamination is associated with a significant excess mortality in free-ranging bees. This provides a strong and unprecedented link between predictions drawn from artificial exposure experiments [9–11,40] and evidence from real-field surveys. Second, colonies appeared to be able to compensate for the excess mortality so as to preserve unaltered performance in terms of population size and honey production. Instead, the most exposed colonies modified the timing of their reproductive investment, delaying
drone brood production in favour of increased worker brood production. We have now reconciled the conflicting laboratory and field assessments of neonicotinoid toxicity. It is thus urgent that risk assessors take into account the scientific evidence for behavioural disorders triggered by trace levels of neonicotinoids.

**Data accessibility.** Detailed data and statistics supporting the results are available in the electronic supplementary material.

**Authors' contribution.** M.H. and N.C. contributed equally to this work. All authors contributed to reviewing and editing the manuscript. N.C., V.B. and A.D. conceived and designed the study. N.C., J.F.O., P.A., C.R., M.G. and A.P. coordinated the fieldwork and collected the data. M.H., N.C., C.R., M.G. and V.B. analysed the data. M.H., N.C., A.D., J.F.O., P.A. and V.B. prepared the manuscript.

**Competing interests.** We declare we have no competing interests.

**Funding.** Funds for equipment and fieldwork were provided by Terres Inovia, the French Ministry of Agriculture and the European Community programme (797/2004) for French beekeeping coordinated by the French Ministry of Agriculture (RISQAPI project). ACTA and INRA SPE Division provided funds for field assistance.

**Acknowledgements.** Special thanks go to volunteer farmers involved in the study. T. Perrot, M. Liaigre, M. Gourrat, F. Allier and I. Badenhausser interviewed farmers, and R. Drieu helped with archiving farmer data. We thank A. Merrien, A. Pouzet and two anonymous reviewers for their useful comments on the manuscript, F. Requier for fruitful discussions on honeybee life-history traits and H. Coleman for language corrections. The Terres Inovia collaborators from INRA Magneraud as well as the INRA UE Entomologie technical staff provided extensive help in the field. RFID devices were developed and updated by E. Lewden from Tag Tracing Solutions (Valence, France). Neonicotinoid multi-residues were assayed by A.-C. Martel at the ANSES laboratory (Sophia-Antipolis, France).

## References

1. Godfray HJ, Blacquière T, Field LM, Hails RS, Petrokofsky G, Potts SG, Raine NE, Vanbergen AJ, McLean AR. 2014 A restatement of the natural science evidence base concerning neonicotinoid insecticides and insect pollinators. *Proc. R. Soc. B* **281**, 20140558. (doi:10.1098/rspb.2014.0558)

2. Sánchez-Bayo F. 2014 The trouble with neonicotinoids. *Science* **346**, 806–807. (doi:10.1126/science.1259159)

3. Goulson D, Nicholls E, Botías C, Rotheray EL. 2015 Bee declines driven by combined stress from parasites, pesticides, and lack of flowers. *Science* **347**, 1255957. (doi:10.1126/science.1255957)

4. Simon-Delso N et al. 2015 Systemic insecticides (neonicotinoids and fipronil): trends, uses, mode of action and metabolites. *Environ. Sci. Pollut. Res.* **22**, 5–34. (doi:10.1007/s11356-014-3470-y)

5. Desneux N, Decourtye A, Delpuech JM. 2007 The sublethal effects of pesticides on beneficial arthropods. *Annu. Rev. Entomol.* **52**, 81–106. (doi:10.1146/annurev.ento.52.110405.091440)

6. Decourtye A, Devillers J, Cluzeau S, Charreton M, Pham-Delègue MH. 2004 Effects of imidacloprid and deltamethrin on associative learning in honeybees under semi-field and laboratory conditions. *Ecotox. Environ. Safe* **57**, 410–419. (doi:10.1016/j.ecoenv. 2003.08.001)

7. El Hassani AK, Dacher M, Gary V, Lambin M, Gauthier M, Armengaud C. 2008 Effects of sublethal doses of acetamiprid and thiamethoxam on the behavior of the honeybee (*Apis mellifera*). *Arch. Environ. Contam. Toxicol.* **54**, 653–661. (doi:10.1007/s00244-007-9071-8)

8. Yang EC, Chuang YC, Chen YL, Chang LH. 2008 Abnormal foraging behavior induced by sublethal dosage of imidacloprid in the honey bee (Hymenoptera: Apidae). *J. Econ. Entomol.* **101**, 1743–1748. (doi:10.1603/0022-0493-101.6.1743)

9. Henry M, Béguin M, Requier F, Rollin O, Odoux J-F, Aupinel P, Aptel J, Tchamitchian S, Decourtye A. 2012 A common pesticide decreases foraging success and survival in honey bees. *Science* **336**, 348–350. (doi:10.1126/science.1215039)

10. Henry M, Bertrand C, Le Féon V, Requier F, Odoux J-F, Aupinel P, Bretagnolle V, Decourtye A. 2014 Pesticide risk assessment in free-ranging bees is weather and landscape dependent. *Nat. Commun.* **5**, 4359. (doi:10.1038/ncomms5359)

11. Pisa LW et al. 2015 Effects of neonicotinoids and fipronil on non-target invertebrates. *Environ. Sci. Pollut. Res.* **22**, 68–102. (doi:10.1007/s11356-014-3471-x)

12. Thompson HM, Maus C. 2007 The relevance of sublethal effects in honey bee testing for pesticide risk assessment. *Pest. Manag. Sci.* **63**, 1058–1061. (doi:10.1002/ps.1458)

13. Cutler GC, Scott-Dupree CD. 2007 Exposure to clothianidin seed-treated canola has no long-term impact on honey bees. *J. Econ. Entomol.* **100**, 765–772. (doi:10.1603/0022-0493(2007)100:765:ETCSCH)2.0.C0;2)

14. Pilling E, Campbell P, Coulson M, Ruddle N, Tornier I. 2013 A four-year field program investigating long-term effects of repeated exposure of honey bee colonies to flowering crops treated with thiamethoxam. *PLoS ONE* **8**, e77193. (doi:10.1371/journal.pone.0077193)

15. Cutler GC, Scott-Dupree CD, Sultan M, McFarlane AD, Brewer L. 2014 A large-scale field study examining effects of exposure to clothianidin seed- treated canola on honey bee colony health, development, and overwintering success. *PeerJ* **2**, e652. (doi:10.7717/peerj.652)

16. Rundlöf M et al. 2015 Seed coating with a neonicotinoid insecticide negatively affects wild bees. *Nature* **521**, 77–80. (doi:10.1038/ nature14420)

17. Carreck NL, Ratnieks FLW. 2014 The dose makes the poison: have 'field realistic' rates of exposure of bees to neonicotinoid insecticides been overestimated in laboratory studies? *J. Api. Res.* **53**, 607–614. (doi:10.3896/IBRA.1.53.5.08)

18. Dively GP, Embrey MS, Kamel A, Hawthorne DJ, Pettis JS. 2015 Assessment of chronic sublethal effects of imidacloprid on honey bee colony health. *PLoS ONE* **10**, e0118748. (doi:10.1371/journal.pone.0118748)

19. Erickson B. 2013 Regulation: Europe bans three neonicotinoids linked to honeybee population declines. *Chem. Eng. News Archive* **91**, 11. (doi:10.1021/cen-09118-notw9)

20. Streit S, Bock F, Pirk CWW, Tautz J. 2003 Automatic life-long monitoring of individual insect behaviour now possible. *Zoology* **106**, 169–171. (doi:10.1078/0944-2006-00113)

21. Decourtye A, Devillers J, Aupinel P, Brun F, Bagnis C, Fourrier J, Gauthier M. 2011 Honeybee tracking with microchips: a new methodology to measure the effects of pesticides. *Ecotoxicology* **20**, 429–437. (doi:10.1007/s10646-011-0594-4)

22. Odoux J-F, Aupinel P, Gateff S, Requier F, Henry M, Bretagnolle V. 2014 ECOBEE: a tool for long-term bee colony monitoring at landscape scale in West European intensive agrosystems. *J. Api. Res.* **53**, 57–66. (doi:10.3896/IBRA.1.53.1.05)

23. Beshers SN, Fewell JH. 2001 Models of division of labor in social insects. *Ann. Rev. Entomol.* **46**, 413–440. (doi:10.1146/annurev.ento.46.1.413)

24. Leoncini I, Causer D, Robinson GE, Conte YL. 2004 Worker-worker inhibition of honey bee behavioural development independent of queen and brood. *Insect. Soc.* **51**, 392–394. (doi:10.1007/s00040-004-0757-x)

25. Perry CJ, Savik E, Myerscough MR, Barron AB. 2015 Rapid behavioral maturation accelerates failure of stressed honey bee colonies. *Proc. Natl Acad. Sci. USA* **112**, 3427–3432. (doi:10.1073/pnas.1422089112)

26. Williams GR et al. 2013 Standard methods for maintaining adult *Apis mellifera* in cages under *in vitro* laboratory conditions. *J. Api. Res.* **52**, 1–36. (doi:10.3896/IBRA.1.52.1.04)

27. Seeley TD. 1982 Adaptive significance of the age polyethism schedule in honeybee colonies. *Behav. Ecol. Sociobiol.* **11**, 287–293. (doi:10.1007/ BF00299306)

28. Winston ML, Punnett EN. 1982 Factors determining temporal division of labor in honeybees. *Can. J. Zool.* **60**, 2947–2952. (doi:10.1139/z82-372)

29. Capaldi EA et al. 2000 Ontogeny of orientation flight in the honeybee revealed by harmonic radar. *Nature* **403**, 537–540. (doi:10.1038/3500056)
30. Fox J. 2002 Cox proportional-hazards regression for survival data. In *An R and S-PLUS companion to applied regression* (ed. J Fox), pp. 1-20. Thousand Oaks, CA: Sage.

31. R Development Core Team. 2008 *R: a language and environment for statistical computing*. Vienna, Austria: R Foundation for Statistical Computing.

32. Corbet SA, Cuthill I, Fallows M, Harrison T, Hartley G. 1981 Why do nectar-foraging bees and wasps work upwards on inflorescences? *Oecologia* **51**, 79-83. (doi:10.1007/BF00344656)

33. Martel A-C, Mangoni P, Gastaldi-Thiery C. 2013 Determination of neonicotinoid residues in nectar by liquid chromatography coupled with tandem mass spectrometry (LC-MS/MS). *EuroReference J.* **11**, 18-21.

34. Rondeau G, Sánchez-Bayo F, Tennekes HA, Decourtye A, Ramírez-Romero R, Desneux N. 2014 Delayed and time-cumulative toxicity of imidacloprid in bees, ants and termites. *Sci. Rep.* **4**. (doi:10.1038/srep05566)

35. Prisco GD, Cavaliere V, Annoscia D, Varricchio P, Caprio E, Nazzi F, Gargiulo G, Pennacchio F. 2013 Neonicotinoid clothianidin adversely affects insect immunity and promotes replication of a viral pathogen in honey bees. *Proc. Natl Acad. Sci. USA* **110**, 18 466-18 471. (doi:10.1073/pnas.1314923110)

36. Maxim L, van der Sluijs JP. 2007 Uncertainty: cause or effect of stakeholders' debates?: Analysis of a case study: the risk for honeybees of the insecticide

Gaucho®. *Sci. Total Environ.* **376**, 1-17. (doi:10.1016/j.scitotenv.2006.12.052)

37. Park MG, Blitzer EJ, Gibbs J, Losey JE, Danforth BN. 2015 Negative effects of pesticides on wild bee communities can be buffered by landscape context. *Proc. R. Soc. B* **282**, 20150299. (doi:10.1098/rspb.2015.0299)

38. Sokal RR, Rohlf FJ. 1994 *Biometry: the principles and practice of statistics in biological research*, 3rd revised edn. New York, NY: W. H. Freeman & Co Ltd.

39. Fischer J, Müller T, Spatz A-K, Greggers U, Grünewald B, Menzel R. 2014 Neonicotinoids interfere with specific components of navigation in honeybees. *PLoS ONE* **9**, e91364. (doi:10.1371/journal.pone.0091364)

40. Karahan A, Çakmak I, Hranitz JM, Karaca I, Wells H. In press. Sublethal imidacloprid effects on honey bee flower choices when foraging. *Ecotoxicology*. (doi:10.1007/s10646-015-1537-2)

41. Thompson H, Coulson M, Ruddle N, Wilkins S, Harkin S. In press. Thiamethoxam: assessing flight activity of honeybees foraging on treated oilseed rape using RFID technology. *Environ. Toxicol. Chem.* (doi:10.1002/etc.3183)

42. Boes K. 2010 Honeybee colony drone production and maintenance in accordance with environmental factors: an interplay of queen and worker decisions. *Insect. Soc.* **57**, 1-9. (doi:10.1007/s00040-009-0046-9)

43. Gill RJ, Ramos-Rodriguez O, Raine NE. 2012 Combined pesticide exposure severely affects individual- and colony-level traits in bees. *Nature* **491**, 105-108. (doi:10.1038/nature11585)

44. Johnson RM, Dahlgren L, Siegfried BD, Ellis MD. 2013 Acaricide, fungicide and drug interactions in honey bees (*Apis mellifera*). *PLoS ONE* **8**, e54092. (doi:10.1371/journal.pone.0054092)

45. Samson-Robert O, Labrie G, Chagnon M, Fournier V. 2014 Neonicotinoid-contaminated puddles of water represent a risk of intoxication for honey bees. *PLoS ONE* **9**, e108443. (doi:10.1371/journal.pone.0108443)

46. Douglas MR, Rohr JR, Tooker JF. 2015 Neonicotinoid insecticide travels through a soil food chain, disrupting biological control of non-target pests and decreasing soya bean yield. *J. Appl. Ecol.* **52**, 250-260. (doi:10.1111/1365-2664.12372)

47. Schaafsma A, Limay-Rios V, Baute T, Smith J, Xue Y. 2015 Neonicotinoid insecticide residues in surface water and soil associated with commercial maize (corn) fields in southwestern Ontario. *PLoS ONE* **10**, e0118139. (doi:10.1371/journal.pone.0118139)

48. Botías C, David A, Horwood J, Abdul-Sada A, Nicholls E, Hill E, Goulson D. 2015 Neonicotinoid residues in wildflowers, a potential route of chronic exposure for bees. *Environ. Sci. Technol.* **49**, 12 731-12 740. (doi:10.1021/acs.est.5b03459)
# Predictive systems models can help elucidate bee declines driven by multiple combined stressors

Mickaël HENRY<fnref n="1" />,<fnref n="2" />, Matthias A. BECHER<fnref n="3" />, Juliet L. OSBORNE<fnref n="3" />, Peter J. KENNEDY<fnref n="3" />, Pierrick AUPINEL<fnref n="4" />, Vincent BRETAGNOLLE<fnref n="5" />,<fnref n="6" />, François BRUN<fnref n="7" />, Volker GRIMM<fnref n="8" />, Juliane HORN<fnref n="8" />, Fabrice REQUIER<fnref n="1" />,<fnref n="2" />

Received 30 March 2016 – Revised 4 August 2016 – Accepted 22 September 2016

**Abstract** – Bee declines are driven by multiple combined stresses, making it exceedingly difficult to identify experimentally the most critical threats to bees and their pollination services. We highlight here the too often ignored potential of mechanistic models in identifying critical stress combinations. Advanced bee models are now available as open access tools and offer an unprecedented opportunity for bee biologists to explore bee resilience tipping points in a variety of environmental contexts. We provide general guidelines on how to run bee models to help detect a priori critical stress combinations to be targeted in the field. This so-called funnel analysis should be performed in tight conjunction with the recent development of large-scale field monitoring programs for bee health surveillance.

*Apis mellifera* / field monitoring program / honeybees / mechanistic modeling / agent-based models

## 1. INTRODUCTION

The early reviews on the causes of pollinator declines (Potts et al. 2010; Vanbergen and the Insect Pollinators Initiative 2013) highlighted the potential role of multiple stressor interactions and synergies. Recently, Goulson et al. (2015) updated the state of knowledge on combined stressors on bees, reviewing evidence from in vitro and in vivo experiments. They concluded that disentangling the effects of multiple combined stressors is

exceedingly difficult, implying that evidence will not appear rapidly enough given the urgency of the situation. Here, we would like to draw bee biologists' attention to an under-used shortcut that may help bypass some technical barriers, namely the in silico pathway (Evans et al. 2013; Stillman et al. 2015). Mechanistic ecological modeling has matured over the last decade and leads to more predictive and flexible models, which are increasingly used for solving a wide range of applied problems (Stillman et al. 2015). We provide some guidelines on how to run bee mechanistic models to help detect critical stress combinations from large-scale monitoring program datasets. We concur with Goulson et al. (2015) that it is necessary to develop pollinator monitoring programs, but

<footnote n="1">¹INRA, UR406 Abeilles et Environnement, F-84914, Avignon, France</footnote>

<footnote n="2">²UMT Protection des Abeilles dans l'Environnement, Site Agroparc, F-84914, Avignon, France</footnote>

<footnote n="3">³Environment and Sustainability Institute, University of Exeter, Penryn Campus, Penryn, Cornwall TR10 9FE, UK</footnote>

<footnote n="4">⁴INRA, UE1255, UE Entomologie, F-17700, Surgères, France</footnote>

<footnote n="5">⁵Centre d'Etudes Biologiques de Chizé, UMR 7372, CNRS & Université de La Rochelle, F-79360, Beauvoir-sur-Niort, France</footnote>

<footnote n="6">⁶LTER Zone Atelier Plaine & Val de Sèvre, Centre d'Etudes Biologiques de Chizé, CNRS, F-79360, Villiers-en-Bois, France</footnote>

<footnote n="7">⁷ACTA INRA, UMR 1248 AGIR, F-31326, Castanet Tolosan cedex, France</footnote>

<footnote n="8">⁸UFZ, Helmholtz Centre for Environmental Research – UFZ, Permoserstr. 15, DE-04318, Leipzig, Germany</footnote>

<footnote>Corresponding author: M. Henry, mickaël.henry@paca.inra.fr</footnote>

<footnote>Manuscript editor: David Tarpy</footnote>
the field design of those programs should, whenever possible, target the capture of information on critical combinations of stressors predicted in silico.

Ecological data collected at landscape scale (i.e., hundreds of km²) under question-oriented monitoring surveys has propelled a reappraisal of environmental stressors in the observed decline of wild bees and honeybees. For instance, in the ECOBEE program (Odoux et al. 2014), a strong focus has been developed toward a better understanding of the relationships between honeybee colony traits and the environment. This has demonstrated that food availability appears as one of the main drivers of honeybee colony dynamics (Requier et al. 2014), and temporal limitations in food availability between mass flowering bloom events (Bretagnolle and Gaba 2015) generate strong cascading effects on colony health (Requier et al. 2016; see also Horn et al. 2016). Still, other stressors can act synergistically so that their respective effect on bee survival may be hard to discriminate in the field (Henry et al. 2014, 2015), and there is accumulating evidence that combined stresses actually exacerbate mortality in honeybees and wild bees. Excess mortality has been empirically evidenced mostly with two-way stressor interactions, including parasites × pathogens (Nazzi et al. 2012; Retschnig et al. 2014), poor nutrition × pathogens (Di Pasquale et al. 2013), pesticides × pathogens (Alaux et al. 2010; Aufauvre et al. 2012; Doublet et al. 2015; Fauser-Misslin et al. 2014), pesticides × pesticides (Vidau et al. 2011; Gill et al. 2012; Johnson et al. 2013; Henry et al. 2015), and pesticides × weather or landscape context combinations (Henry et al. 2014; Park et al. 2015). See also Collison et al. (2015) for a thorough quantitative review of pesticides × pathogens case studies.

However, the number of possible combinations increases exponentially as more candidate stressors are considered, eventually making it prohibitively expensive to cover the whole range of possibilities through traditional in vitro and in vivo approaches (Goulson et al. 2015). Alternatively, the in silico pathway provides a useful framework for combining and assessing complex stressor interactions in bees (Barron 2015). The idea of running bee models to simulate in silico

any multiple-way combination of stressors has recently emerged (Becher et al. 2013), stemming from the current flourishing of literature on bee mechanistic models addressing specific issues (e.g., honeybee colonies [Schmickl and Crailsheim 2007; Becher et al. 2010, 2014; Khoury et al. 2011, 2013; Torres et al. 2015], bumblebee colonies [Bryden et al. 2013], solitary bees [Everaars and Dormann 2014]). The in silico pathway may be used to investigate tipping points under combined stressors and has provided significant improvement on the understanding of bee mortality (Henry et al. 2012; Perry et al. 2015; Horn et al. 2016; Rumkee et al. 2015). Tipping points (Figure 1) refer here to stress levels beyond which a small variation translates into a drastic constraint on the sustainability of the whole system, i.e., a steep increase in the probability to reach a given endpoint threshold (typically the bee population or colony collapse).

Some of these open source bee models are expected to be used by an increasing number of bee biologists for the study of combined stresses (Barron 2015) and recalibrated for a variety of specific bioclimatic contexts and research issues. In that respect, we believe that a coherent research agenda is urgently needed, along with a unified terminology to facilitate communication among

<figure><img src="image_50.png" /></figure>

**Figure 1.** Predictive systems models may be used to investigate tipping points in bee populations under environmental stresses. Tipping points denote the critical stress levels beyond which a slight increase translates into a steep change in the state of the whole system, e.g., a drastic rise of the probability to reach a given monitored endpoint threshold (*double arrows*). In honeybees, typical endpoints of interest are colony collapse or low population size (e.g., <4000 individuals at the end of the season; Becher et al. 2014; Horn et al. 2016; Rumkee et al. 2015).
users. We do not plead however for an in silico approach as a substitute for empirical studies. Instead, it should be viewed as a tool to guide empirical studies throughout the immense scope of possibilities down to a reduced subset of plausible and testable predictions, that is to say, *funneling* candidate scenarios into the conception of both targeted small-scale experiments and of large field monitoring programs (e.g., the ECOBEE monitoring survey in France [Odoux et al. 2014], the German bee monitoring project [Genersch et al. 2010], the US national surveys [van Engelsdorp et al. 2008], the EPILOBEE European survey [EPILOBEE Consortium et al. 2016]).

We called this approach *funnel analysis* and provide here general guidelines for using it, with special emphasis on the well-documented honeybee, though the same procedure applies to models targeting other bee species. We first describe the procedure in the *Material and methods* section, along with an appropriate terminology. Second, by using BEEHAVE as a case study honeybee system model (Becher et al. 2014), we tested whether demographic tipping points may be identified when simulating individual stresses on important colony traits taken separately (e.g., forager mortality, larvae mortality, *Varroa mite* infestation, foraging activity) and whether combining those stresses would hasten the tipping point, as is expected from current hypotheses on the multifactorial nature of bee declines. We finally discuss the forthcoming opportunities for bee colony model users to coordinate their efforts for a better understanding of the environmental drivers of bee declines.

## 2. MATERIAL AND METHODS

### 2.1. Prerequisite for the multistress funnel analysis

Models, even if they are complex, are simplified representations. Therefore, to base inference about real bees on model output, all model assumptions and their justifications need to be documented so that they can be scrutinized. For example, the *Overview, Design concepts, Details* (ODD; Grimm et al. 2006, 2010) standard documentation of the BEEHAVE honeybee colony

model (Becher et al. 2014) allowed a European expert group to evaluate all model assumptions (European Food Safety Authority 2015). Furthermore, models should be rich enough in structure and processes so that they can be empirically validated at different hierarchical levels (DeAngelis and Wolf 2003; Grimm et al. 2005). For example, Horn et al. (2016) used the BEEHAVE model to identify time periods within a foraging season where honeybee colonies are most sensitive to gaps in forage supply. They then trace these findings back to specific changes in colony structure and dynamics, which can be tested in subsequent targeted small-scale experiments. Finally, to enable testing of a model under different environmental settings, effects of the most important environmental drivers should emerge from first principles instead of imposing them via empirical rules. For honeybee models, this means making colony dynamics dependent on weather and use the availability, distribution, and foraging of nectar and pollen so that the same model can be used in different landscapes (e.g., Becher et al. 2014; Horn et al. 2016; Rumkee et al. 2015).

If a model fulfills these requirements, it can be used for what we here refer to as funnel analysis (Figure 2). For such analyses, one needs to (i) implement appropriate climate and/or landscape *input* model parameters for the bioclimatic context of interest and (ii) choose a relevant endpoint threshold among the bee model *outputs*. Simulating combined stresses on bees requires that users carefully decide what endpoint critical values are to be considered for drawing conclusions from simulations. In the case of honeybees, the typical endpoint of interest for bee biologists and beekeepers is colony survival over a season. Other endpoints may include colony size at the end of the year, peak colony size, economic value of honey stores or pollination services, or production of queens and offspring in wild bees. In general, endpoints should be used that can demonstrate thresholds and tipping points (Figure 1) as this translates the output of even complex models into a format that can directly inform management and policy development (Stillman et al. 2016).

### 2.2. Essentials of the multistress funnel analysis

The multistress funnel analysis we propose here (Figure 2) is based on a mixture of *implicit* and
<figure><img src="image_51.png" /><figcaption>**Figure 2.** The multistressor funnel analysis applied to honeybee colony dynamic models. In silico virtual colony (Schmickl and Crailsheim 2007; Becher et al. 2010, 2014; Khoury et al. 2011, 2013; Torres et al. 2015) simulations should sequentially drop combined stressor scenarios until delivering a limited subset of plausible scenarios testable as a part of field monitoring programs (van Engelsdorp et al. 2008; Genersch et al. 2010; Odoux et al. 2014). See *Material and methods* section for step-by-step details.</figcaption></figure>

explicit representations of stressor effects on core honeybee colony traits. Few stressors are explicitly parameterized in current models, i.e., explicitly taken into account via specific parameters and equations. This is the case with *Varroa mites* and the deformed wing virus (DWV) whose infestation levels are tunable as a part of the input parameters in BEEHAVE virtual colonies (Becher et al. 2014). Distance and quantity of nectar and pollen available are also explicitly defined in BEEHAVE as well as the effects of lack of pollen and nectar on the brood (Becher et al. 2014; Horn et al. 2016). Still, many candidate honeybee stressors have to be implicitly implemented by simulating their effect on core colony life history traits, i.e., those colony traits that most influence the fate of the colony (step A1 in Figure 2). For instance, the excess mortality due to homing failure in intoxicated forager bees was implicitly implemented by inflating the foraging mortality parameter (Henry et al. 2012; Becher et al. 2014). Core colony traits may be identified by means of model sensitivity

analyses (Saltelli et al. 2008; Thiele et al. 2014). Such analyses can comprise factorial designs of stressors which are then condensed into variance partitioning methods (Horn et al. 2016) or via screening for the most sensitive parameters and then exploring, via global analysis, the relative importance of these processes and their interactions. Typical core colony traits concern the baseline demographic functioning of the colony (e.g., brood mortality, in-hive workers, and foragers mortality; Schmickl and Crailsheim 2007; Becher et al. 2014; Rumkee et al. 2015).

Therefore, the funnel analysis approach we recommend here actually consists of simulating variations of one or multiple core colony traits and identifying which scenarios may lead to the endpoint thresholds of interest (steps A2 and A3 in Figure 2) and then which ones are realistic enough, in terms of stress realism and parsimony, to be tested in the field (steps A4 and B in Figure 2). Those steps are further detailed below and illustrated in Figure 2:
(A) From the full scope of possibilities down to the most plausible subset of critical stress combinations.

(A1) *Stressor effect assignment*: identifying from literature and expert knowledge a qualitative network of effects linking the candidate stressors (e.g., *nosema* spores, neonicotinoid pesticides, floral resource scarcity) with the affected core colony traits (e.g., foraging intensity, brood mortality, in-hive workers, and foragers mortality).

(A2) *Heuristic simulation of effects on core colony traits*: simulating all possible combinations of core colony trait variations. Trait variations may be either persistent throughout the season or restricted to a particular period. Likewise, the magnitude of trait variations may be allowed to vary within a realistic range. One may compute broad variation scenarios (e.g., low, medium, or large effect size) or series of scenarios with finer resolution depending on computational performance.

(A3) *Screening candidate scenarios for specific goals and endpoint criteria*: The candidate scenarios that result in a high probability of colonies reaching the predetermined endpoint threshold of interest are singled out. The threshold probability for retaining and rejecting scenarios is decided before the analyses, based on goals for protection of colonies or risk deemed acceptable. One example of specific goals being applied is in the pesticide regulatory arena. In line with the European Food Safety Authority guidance for acceptable honeybee protection goals, Rumkee et al. (2015) proposed the lethal imposed stress ($LIS_x$) framework for colonies, analogous to the standard lethal dose ($LD_{50}$) currently used for pesticide risk assessment on an individual basis. The $LIS_{50}$ (or $LIS_{10}$) indicates the level of imposed stress resulting in 50% (or 10%) probability of colonies reaching a chosen endpoint threshold. Using such a framework would screen out the scenarios that predict less than a 50% (or 10%) chance of reaching the chosen threshold. Alternatively, a cutoff rule may be based on the statistical power of the envisioned field monitoring, i.e., the finest detectable

change of endpoint threshold occurrence frequency given actual sample size.

(A4) *Screening candidate scenarios for realism and parsimony criteria*: Within the subset of candidate scenarios meeting the endpoint criteria, some are likely to be more realistic than others considering stressor temporal or spatial occurrence. Expert knowledge is required at this step to decide which scenario should receive priority given field realism and parsimony of stressor combinations. Experts may set aside candidate scenarios with unrealistically high stress levels or with stressors occurring at inappropriate periods of the season. For instance, some experts in a given eco-region may be primarily concerned by the effects of increased forager mortality due to neonicotinoid residuals during oilseed rape blooms and of inhibited pollen foraging activity due to the presence of attacking yellow-legged hornets at hive entrance. Those stresses occur, respectively, in spring (Henry et al. 2015) and in early fall (Rome et al. 2015). Experts may therefore recommend focusing on the corresponding simulations with acute forager mortality in spring and reduced foraging activity in fall in order to further assess potential associations with other stressors throughout the season, such as pathogens or climate hazards. Conversely, experts may judge it unlikely that all life stages (eggs, larvae, pupae, in-hive workers, established foragers) be simultaneously subject to acute mortality. For the sake of parsimony, they may recommend selecting simulations with only one or two life stages being simultaneously affected by acute mortality levels.

(B) Priority field experiments.

The most plausible and parsimonious scenarios may finally be addressed in the field, either by suggesting new hypotheses or novel dataset cross-comparisons to be investigated in already existing large-scale monitoring programs. Field (in vivo) monitoring surveys should be aimed at confirming the links between endpoint threshold occurrence and core colony trait alterations. Laboratory (in vitro) experiments may further help untangle the explicit links with specific stressors. Regardless
of the chosen monitoring design, time series ana-
lytical approaches would be requested, whereby
correlative links would be tested between the pres-
ence or magnitude of candidate stressors, the core
colony traits at a given time of the season, and the
subsequent state of colonies later on (endpoint
threshold reached or not).

## 2.3. Assessing combined effects of multiple stressors

To reveal critical stressor combinations, predictive
systems models should be able to reveal combined
effects, i.e., they should predict that multiple stressors
produce a greater effect than individual stressors taken
separately. To test this assertion, we compared in silico
the fate of honeybee colonies with a series of stresses set
at the tipping point and simulated either separately or in
combinations. Although purely illustrative, this simple
example shows how independent stressors can combine
their respective effect to promote colony failure.

*Initial settings* We used the BEEHAVE model (Becher et al. 2014) with Horn's et al. (2016) initial colony settings, and we implemented a typical weather and landscape dataset from the ECOBEE colony monitoring area (Odoux et al. 2014). We allowed a baseline infestation with 100 virus-infected *Varroa* mites on the first day of simulation as default value. Current beekeeping practices were also simulated with the ad hoc options (default settings for *Varroa* treatment and honey harvests, see Becher et al. 2014).

*Endpoint* Colony sizes smaller than 4000 adult bees at the end of the year (Becher et al. 2014; Horn et al. 2016; Rumkee et al. 2015) were deemed insufficient for sustainable honey production and survival throughout the beekeeping season without human intervention. We therefore used this demographic threshold as an endpoint value to illustrate the combined stress analysis. Colony simulations were run for 3 years or until the endpoint threshold was reached.

*Core colony traits and tipping points* We focused on four core colony traits that are thought to have profound impact on demographic trajectories: forager mortality, larvae mortality, virus-infected *Varroa* mites, and reduction of the foraging activity. We performed a sensitivity analysis to identify critical tipping points

associated with those core colony traits taken individu-
ally. Default model parameters were gradually increased
(larvae and forager mortality rate, *Varroa* infestation) or
decreased (maximal foraging distance allowed) to de-
lineate the endpoint threshold occurrence probability
with increasing stress values. First, exploratory simula-
tions were run by series of 10 repetitions, until detecting
the range of stress values leading to low (n = 1) and
high (n = 10) endpoint threshold occurrence. The range
was then systematically covered by 20 series of 10
simulation repetitions for a better resolution. Simulation
outputs were finally fitted with a standard logistic func-
tion by using a generalized linear model (GLM) for
binomial data (*glm* function in R software for
statistical computing, R Core Team 2014). The tipping
point level was taken to be the highest stress value the
simulated colonies could afford without reaching the
endpoint threshold.

*Stress implementation* We tested the hypothesis that combining two or more stresses at their tipping point (the others being set to the default value) would precipitate colonies to endpoint, which would normally not occur with a single stressor. All multiple-way tipping point combinations (n = 16) were simulated 20 times each, totaling 320 simulated colonies. The probability to reach the endpoint threshold out of the 20 simulations was computed and compared by using simple a posteriori pairwise exact binomial tests among stress combinations (*binom.test* function, R Core Team 2014).

3. RESULTS

*Core colony traits and tipping points* A tipping point pattern was successfully established when tuning each of the four studied core colony traits (Figure 3) and satisfactorily modeled by using logistic functions (Table I). Estimated tipping points were 1.42 for forager mortality rate per 10⁵s foraging (default = 1), 0.06 for daily larval mortality rate (default = 0.01), 1142 for initial number of infected mites (default = 100), and 3 km for maximal distance flown by foragers on a single day (default limited by duration of clem-ent weather for foraging).

*Stress implementation* Of the 320 simulated colo-
nies covering the 16 possible tipping point
<figure><img src="image_52.png" /><figcaption>**Figure 3.** Assessment of tipping points associated with four core colony traits. The probability to reach the endpoint threshold (<4000 adult honeybees) dramatically increases beyond a critical stress value referred to as tipping points (vertical line). Dots and bars stand for the mean ± SD endpoint threshold occurrence per series of 10 replicated colony simulations. Trends are depicted by a logistic fit to simulation data (thick lines and 95 % confidence limit as shaded areas, see Table I for details).</figcaption></figure>

combinations, 14 collapsed or fell below the end-point threshold in the course of the first year, 120 in the second year, and 15 in the third year. Simulation outputs merely illustrate the combined stressor hypothesis developed in current scientific

literature (Figure 4). Taken alone, stressors at tipping point levels did not precipitate simulated colonies to reach the endpoint threshold. For instance, colonies with a 42 % increase in forager mortality were already on the brink of reaching

<table><caption>**Table 1.** Logistic model statistics for predicting of the probability to reach the endpoint threshold as a function of stress level imposed on four core colony traits.</caption><thead><tr><th>Core colony trait value x</th><th>α</th><th>β</th><th>Chi-squared statistics for deviance</th></tr></thead><tbody><tr><td>Forager mortality rate (per 10<sup>5</sup> s foraging)</td><td>9.77 ± 1.73</td><td>-17.49 ± 3.12</td><td>κ = 64.0, p < 0.001</td></tr><tr><td>Larval mortality rate (per day)</td><td>132.4 ± 26.6</td><td>-12.31 ± 2.50</td><td>κ = 49.7, p < 0.001</td></tr><tr><td>Varroa infestation (nb. of infected mites)</td><td>0.0011 ± 0.00015</td><td>-4.12 ± 0.63</td><td>κ = 115.8, p < 0.001</td></tr><tr><td>Foraging activity (maximal daily foraging distance, km)</td><td>-4.67 ± 1.67</td><td>8.43 ± 3.19</td><td>κ = 26.4, p < 0.001</td></tr></tbody></table>

The predicted endpoint curves (Figure 3) are of the form $p_{Endpoint} = 1 / (1 + e^{-(αx+β)})$, with x the focus colony trait value and α and β the logistic model estimates returned by GLMs
<figure><img src="image_53.png" /><figcaption>**Figure 4.** Comparison of combined stress scenarios for endpoint criteria. For each of the 16 simulated combined stress scenarios (*left panel*), the risk to reach the critical endpoint threshold (*right panel*) is given as a probability value (±SD). Probability values that are identified by distinct *letters* in the right panel are significantly different from each other (a posteriori exact binomial tests among scenarios of combined stresses, see text).</figcaption></figure>

the endpoint threshold, but regulation mechanisms in the colony compensated these forager losses. However, the two-way combinations were associated with a range of significantly higher endpoint threshold occurrences (10 to 90 %, Figure 4), confirming the model ability to identify sharp combined effects contrasting with separate effects of stresses taken alone. With three or more combined stresses, the endpoint is almost inevitably reached (>95 %, Figure 4).

*Toward funnel analysis* This demonstrates that tuning only a couple of core honeybee colony parameters can result in a broad range of collapse scenarios. In practice, the funnel analysis is indeed intended to single out which of those scenarios should receive priority in question-oriented field monitoring programs. From this simplistic hypothetical example (Figure 4), the funnel process

would first screen out the scenarios that had no chance to lead to the endpoint of interest (n = 6 scenarios). Then, experts may decide to screen out scenarios with a 10 % endpoint probability cutoff rule (n = 1 scenario), given that this approximates their maximal field monitoring statistical resolution (see also the LIS₅₀ framework in Rumkee et al. 2015). They may further decide to screen out scenarios with altered foraging activity (n = 5 scenario) because they may consider that this level of stress—maximal daily foraging distance limited down to 3 km per forager only—is unrealistically restrictive. Having reduced candidate scenarios from 16 down to 4, experts may finally adjust their field study so as to target the possible environmental drivers of larval mortality, which delivers the highest endpoint threshold probabilities when combined with *Varroa* mites.
# 4. DISCUSSION

The idea underlying the suggested multistress funnel analysis is to funnel candidate scenarios into explicit field monitoring programs while screening for endpoint criteria and for realism and parsimony of candidate stressor combinations. Having confirmed that bee systems models can reveal critical multistress combinations, we recommend that future research investigates the ability of funnel analyses to effectively end up in a small subset of candidate scenarios. There is no guarantee that it will invariably succeed, but it will at least help identify the most influential factors and discover unexpected ones in the studied systems.

Sensitivity to certain stressors is likely to depend on eco-region, weather, landscape structure, and management, so that focusing monitoring and experiments on the wrong factors in a certain eco-region might lead to enormous waste of resources and time. The ultimate aim of model-based funnel analysis is to screen the thousands of conceivable scenarios of combined stressors down to a short list of plausible scenarios humanly manageable in small-scale experiments and field monitoring programs. In that respect, we fully support Goulson et al.'s (2015) call for developing pollinator monitoring programs, but we further encourage bee biologists to conceive or adjust those programs considering mechanistic model predictions of what might be the most critical stress combinations for the considered endpoint. With this in mind, we expect higher colony resilience in social bees owing to buffering mechanisms as a consequence of colony living and division of labor (Henry et al. 2015; Rundlöf et al. 2015). This makes social bees also a testbed for testing theories and approaches regarding resilience and tipping points. With managed honeybee colonies, for example, we can, in contrast to full ecosystems, manipulate buffer mechanisms and study recovery rates and the ability of colonies to absorb disturbances and maintain

functioning. This is also in line with the plea for a more mechanistic and predictive “systems ecology” (Evans et al. 2013). Solitary wild bees are likely to be more susceptible and hence might provide even more sensitive model organisms to test combined stresses (European Academies Science Advisory Council 2015; Rundlöf et al. 2015), but with different resilience mechanisms.

Bee biologists are now at a crossroads of action and reflection. Although efforts should be spent in the development of operative solutions for sustainable beekeeping and maintaining viable wild bee populations (Goulson et al. 2015), disentangling the respective effects of the many candidate stressors at the roots of bee decline is still required. Simulating virtual bee systems is one promising tool in that respect. Mechanistic models of complex biological systems cannot be viewed as perfectly reflecting reality. However, even though a given model may be lacking empirical support for a specific context of interest, we would not discourage bee biologists from investigating the model simulation outputs as a part of their research, providing they keep in mind that the outputs are conditional predictions, i.e., predictions that are dependent on the assumptions the model is based on. Model assumptions, formalized by a series of parameters and equations, will gain accuracy as developers and users gain empirical knowledge. At the same time, the next steps in model development will be (i) to develop specific modules for interaction effects among stressors, as is currently evidenced in laboratory but still poorly considered in predictive systems models and (ii) to move beyond the dynamics of the colony to the dynamics of the actual populations, so that the densities of managed colonies or native species can be predicted over time and space. Using funnel analysis to match model simulations with empirical work will lead to deeper insight into the mechanisms behind broader scale patterns of colony and even species loss.
## ACKNOWLEDGMENTS

This work was supported by the European Community program (797/2004) for French beekeeping coordinated by the French Ministry of Agriculture (RISQAPI project). JLO, MAB, and PK were funded by grants from the Biotechnology and Biological Sciences Research Council, UK [BB/J014915/1; BB/K014463/1]. We thank Axel Decourtye and Jean-François Odoux for useful discussions on honeybee colony modeling and field monitoring design, as well as two anonymous reviewers for constructive comments on the manuscript.

**Author contributions** Mickaël Henry, Matthias A. Becher, Juliet L. Osborne, Peter J. Kennedy, Pierrick Aupinel, Vincent Bretagnolle, François Brun, Volker Grimm, Juliane Horn, and Fabrice Requier contributed to the design and writing of the study. FR and JH performed the simulations.

### Les modèles prédictifs de systèmes peuvent contribuer à élucider les phénomènes de déclin des abeilles soumises à des stress multiples

*Apis mellifera* / programme de surveillance / abeilles mellifères / modélisation mécaniste / modèles agent-basés

### Vorhersagemodelle können helfen, den Rückgang der Bienen aufgrund multipler Stressoren aufzuklären

*apis mellifera* / feldmonitoring / honigbienen / mechanistische modellierung / agentenbasierte modelle

## REFERENCES

Alaux, C., Brunet, J.-L., Dussaubat, C., Mondet, F., Tchamitchan, S., Cousin, M., Brillard, J., Baldy, A., Belzunces, L.P., Conte, Y.L. (2010) Interactions between *Nosema* microspores and a neonicotinoid weaken honeybees (*Apis mellifera*). *Environ. Microbiol.* **12**, 774–782

Aufauvre, J., Biron, D.G., Vidau, C., Fontbonne, R., Roudel, M., Diogon, M., Viguès, B., Belzunces, L.P., Delbac, F., Blot, N. (2012) Parasite-insecticide interactions: a case study of *Nosema ceranae* and fipronil synergy on honeybee. *Sci. Rep.* **2**, 326

Barron, A.B. (2015) Death of the bee hive: understanding the failure of an insect society. *Curr. Opin. Insect Sc.* **10**, 45–50

Becher, M.A., Hildenbrandt, H., Hemelrijk, C.K., Moritz, R.F.A. (2010) Brood temperature, task division and

colony survival in honeybees: a model. *Ecol. Model.* **221**, 769–776

Becher, M.A., Osborne, J.L., Thorbek, P., Kennedy, P.J., Grimm, V. (2013) Towards a systems approach for understanding honeybee decline: a stocktaking and synthesis of existing models. *J. Appl. Ecol.* **50**, 868–880

Becher, M.A., Grimm, V., Thorbek, P., Horn, J., Kennedy, P.J., Osborne, J.L. (2014) BEEHAVE: a systems model of honeybee colony dynamics and foraging to explore multifactorial causes of colony failure. *J. Appl. Ecol.* **51**, 470–482

Bretagnolle, V., Gaba, S. (2015) Weeds for bees? A review. *Agron. Sustain. Dev.* **35**, 891–909

Bryden, J., Gill, R.J., Mitton, R.A.A., Raine, N.E., Jansen, V.A.A. (2013) Chronic sublethal stress causes bee colony failure. *Ecol. Lett.* **16**, 1463–1469

Collison, E., Hird, H., Cresswell, J., Tyler, C. (2015) Interactive effects of pesticide exposure and pathogen infection on bee health—a critical analysis. *Biol. Rev. in press*, doi: 10.1111/brv.12206

R Core Team (2014) R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna

DeAngelis, D.L., Wolf, M. (2003) In praise of mechanistically rich models, in: Canham, C.D., Cole, J.J., Lauenroth, W.K. (Eds.), Models in ecosystem science, Princeton University Press, Princeton, pp. 63–82

Di Pasquale, G., Salignon, M., Le Conte, Y., Belzunces, L.P., Decourtye, A., Kretzschmar, A., Suchail, S., Brunet, J.-L., Alaux, C. (2013) Influence of pollen nutrition on honey bee health: do pollen quality and diversity matter? *PLoS One* **8**, e72016

Doublet, V., Labarussias, M., de Miranda, J.R., Moritz, R.F.A., Paxton, R.J. (2015) Bees under stress: sublethal doses of a neonicotinoid pesticide and pathogens interact to elevate honey bee mortality across the life cycle. *Environ. Microbiol.* **17**, 969–983

EPILOBEE Consortium, Chauzat, M.-P., Jacques, A., Laurent, M., Bougeard, S., Hendrikx, P., Ribière-Chabert, M. (2016) Risk indicators affecting honeybee colony survival in Europe: one year of surveillance. *Apidologie*, **47**, 348–378.

European Academies Science Advisory Council. (2015) Ecosystem Services, Agriculture and Neonicotinoids, EASAC policy report 26, Halle, Germany

European Food Safety Authority. (2015) Statement on the suitability of the BEEHAVE model for its potential use in a regulatory context and for the risk assessment of multiple stressors in honeybees at the landscape level. *EFSA Journal* **13**, 4125:4216

Evans, M.R., Bithell, M., Cornell, S.J., Dall, S.R.X., Diaz, S., et al. (2013) Predictive systems ecology. *Proc. R. Soc. B* **280**, 20131452

Everaars, J., Dormann, C.F. (2014) Simulation of oliary (non-*Apis*) bees competing for pollen, in Devillers, J. (Ed.), *Silico Bees*, CRC Press, Taylor, Francis Group, Boca Raton, pp. 209–268
Fauser-Misslin, A., Sadd, B.M., Neumann, P., Sandrock, C. (2014) Influence of combined pesticide and parasite exposure on bumblebee colony traits in the laboratory. *J. Appl. Ecol.* **51**, 450–459

Genersch, E., von der Ohe, W., Kaatz, H., Schroeder, A., Otten, C., et al. (2010) The German bee monitoring project: a long term study to understand periodically high winter losses of honey bee colonies. *Apidologie* **41**, 332–352

Gill, R.J., Ramos-Rodriguez, O., Raine, N.E. (2012) Combined pesticide exposure severely affects individual- and colony-level traits in bees. *Nature* **491**, 105–108

Goulson, D., Nicholls, E., Botías, C., Rotheray, E.L. (2015) Bee declines driven by combined stress from parasites, pesticides, and lack of flowers. *Science* **347**, 1255957

Grimm, V., Revilla, E., Berger, U., Jeltsch, F., Mooij, W.M., Railsback, S.F., Thulke, H.-H., Weiner, J., Wiegand, T., DeAngelis, D.L. (2005) Pattern-oriented modeling of agent-based complex systems: lessons from ecology. *Science* **310**, 987–991

Grimm, V., Berger, U., Bastiansen, F., Eliassen, S., Ginot, V., et al. (2006) A standard protocol for describing individual-based and agent-based models. *Ecol. Model.* **198**, 115–126

Grimm, V., Berger, U., DeAngelis, D.L., Polhill, J.G., Giske, J., Railsback, S.F. (2010) The ODD protocol: a review and first update. *Ecol. Model.* **221**, 2760–2768

Henry, M., Béguin, M., Requier, F., Rollin, O., Odoux, J.-F., Aupinel, P., Aptel, J., Tchamitchian, S., Decourtye, A. (2012) A common pesticide decreases foraging success and survival in honey bees. *Science* **336**, 348–350

Henry, M., Bertrand, C., Le Féon, V., Requier, F., Odoux, J.-F., Aupinel, P., Bretagnolle, V., Decourtye, A. (2014) Pesticide risk assessment in free-ranging bees is weather and landscape dependent. *Nat. Commun.* **5**, 4359

Henry, M., Cerrutti, N., Aupinel, P., Decourtye, A., Gayrard, M., Odoux, J.-F., Pissard, A., Rüger, C., Bretagnolle, V. (2015) Reconciling laboratory and field assessments of neonicotinoid toxicity to honeybees. *Proc. R. Soc. B* **282**, 20152110

Horn, J., Becher, M.A., Kennedy, P.J., Osborne, J.L., Grimm, V. (2016) Multiple stressors: using the honeybee model BEEHAVE to explore how spatial and temporal forage stress affects colony resilience. *Oikos* **125**, 1001–1016

Johnson, R.M., Dahlgren, L., Siegfried, B.D., Ellis, M.D. (2013) Acaricide, fungicide and drug interactions in honey bees (*Apis mellifera*). *PLoS One* **8**, e54092

Khoury, D.S., Myerscough, M.R., Barron, A.B. (2011) A quantitative model of honey bee colony population dynamics. *PLoS One* **6** e18491

Khoury, D.S., Barron, A.B., Myerscough, M.R. (2013) Modelling food and population dynamics in honey bee colonies. *PLoS One* **8**, e59084

Nazzi, F., Brown, S.P., Annoscia, D., Del Piccolo, F., Di Prisco, G., Varricchio, P., Della Vedova, G., Cattonaro, F., Caprio, E., Pennacchio, F. (2012) Synergistic

parasite-pathogen interactions mediated by host immunity can drive the collapse of honeybee colonies. *PLoS Pathog.* **8**, e1002735

Odoux, J.-F., Aupinel, P., Gateff, S., Requier, F., Henry, M., Bretagnolle, V. (2014) ECOBEE: a tool for long-term bee colony monitoring at landscape scale in west European intensive agrosystems. *J. Apic. Res.* **53**, 57–66

Park, M.G., Blitzer, E.J., Gibbs, J., Losey, J.E., Danforth, B.N. (2015) Negative effects of pesticides on wild bee communities can be buffered by landscape context. *Proc. R. Soc. B* **282**, 20150299

Perry, C.J., Søvik, E., Myerscough, M.R., Barron, A.B. (2015) Rapid behavioral maturation accelerates failure of stressed honey bee colonies. *Proc. Natl. Acad. Sc. U.S.A.* doi: 10.1073/pnas.1422089112

Potts, S.G., Biesmeijer, J.C., Kremen, C., Neumann, P., Schweiger, O., Kunin, W.E. (2010) Global pollinator declines: trends, impacts and drivers. *Trends Ecol. Evol.* **25**, 345–353

Requier, F., Odoux, J.-F., Tamic, T., Moreau, N., Henry, M., Decourtye, A., Bretagnolle, V. (2014) Honey bee diet in intensive farmland habitats reveals an unexpectedly high flower richness and a major role of weeds. *Ecol. Appl.* **25**, 881–890

Requier, F., Odoux, J.F., Henry, M., Bretagnolle, V. (2016) The carry-over effects of spring pollen shortage negatively impact the colony dynamics and survival of managed honeybees. *J. Appl. Ecol, in press*

Retschnig, G., Williams, G.R., Mehmann, M.M., Yañez, O., de Miranda, J.R., Neumann, P. (2014) Sex-specific differences in pathogen susceptibility in honey bees (*Apis mellifera*). *PLoS One* **9**, e85261

Rome, Q., Muller, F.J., Touret-Alby, A., Darrouzet, E., Perrard, A., Villemant, C. (2015) Caste differentiation and seasonal changes in *Vespa velutina* (Hym.: Vespidae) colonies in its introduced range. *J. Appl. Entomol.* **139**, 771–782

Rumkee, J.C.O., Becher, M.A., Thorbek, P., Kennedy, P.J., Osborne, J.L. (2015) Predicting honeybee colony failure: using the BEEHAVE model to simulate colony responses to pesticides. *Environ Sci. Technol.* **49**, 12879–12887

Rundlöf, M., Andersson, G.K.S., Bommarco, R., Fries, I., Hederström, V., et al. (2015) Seed coating with a neonicotinoid insecticide negatively affects wild bees. *Nature* **521**, 77–80

Saltelli, A., Ratto, M., Andres, T., Campolongo, F., Cariboni, J., Gatelli, D., Saisana, M., Tarantola, S. (2008) Global sensitivity analysis: the primer, John Wiley Sons, Chichester.

Schmickl, T., Crailsheim, K. (2007) HoPoMo: a model of honeybee intracolonial population dynamics and resource management. *Ecol. Model.* **204**, 219–245

Stillman, R.A., Railsback, S.F., Giske, J., Berger, U., Grimm, V. (2015) Making predictions in a changing world: the benefits of individual-based ecology. *Bio-Science*, **65**, 140–150

Stillman, R.A., Wood, K.A., Goss-Custard, J.D. (2016) Deriving simple predictions from complex models to
support environmental decision-making. Ecol. Model.,
10.1016/j.ecolmodel.2015.04.014. in press

Thiele, J.C., Kurth, W., Grimm, V. (2014) Facilitating pa-
rameter estimation and sensitivity analysis of agent-
based models: a cookbook using NetLogo and “R”.
JASSS **17**, 11. doi:10.18564/jasss.2503

Torres, D.J., Ricoy, U.M., Roybal, S. (2015) Modeling honey bee populations. PLoS One **10**, e0130966

van Engelsdorp, D., Hayes, J., Jr., Underwood, R.M., Pettis, J. (2008) A survey of honey bee colony losses

in the U.S., fall 2007 to spring 2008. PLoS One 3,
e4071

Vanbergen, A.J., the Insect Pollinators Initiative. (2013)
Threats to an ecosystem service: pressures on pollina-
tors. Front. Ecol. Environ. **11**, 251–259

Vidau, C., Diogon, M., Aufauvre, J., Fontbonne, R., Viguès, B., et al. (2011) Exposure to sublethal doses of fipronil and thiacloprid highly increases mortality of honeybees previously infected by *Nosema ceranae*. PLoS One **6**, e21550
SCIENTIFIC REPORTS

<figure><img src="image_54.png" /></figure>

OPEN

# Controlling the impact of the managed honeybee on wild bees in protected areas

Mickaël Henry & Guy Rodet

Received: 15 December 2017

Accepted: 1 June 2018

Published online: 18 June 2018

In recent years, conservation biologists have raised awareness about the risk of ecological interference between massively introduced managed honeybees and the native wild bee fauna in protected natural areas. In this study, we surveyed wild bees and quantified their nectar and pollen foraging success in a rosemary Mediterranean scrubland in southern France, under different conditions of apiary size and proximity. We found that high-density beekeeping triggers foraging competition which depresses not only the occurrence (–55%) and nectar foraging success (–50%) of local wild bees but also nectar (–44%) and pollen (–36%) harvesting by the honeybees themselves. Overall, those competition effects spanned distances of 600–1.100 m around apiaries, i.e. covering 1.1–3.8km² areas. Regardless the considered competition criterion, setting distance thresholds among apiaries appeared more tractable than setting colony density thresholds for beekeeping regulation. Moreover, the intraspecific competition among the honeybees has practical implications for beekeepers. It shows that the local carrying capacity has been exceeded and raises concerns for honey yields and colony sustainability. It also offers an effective ecological criterion for pragmatic decision-making whenever conservation practitioners envision progressively reducing beekeeping in protected areas. Although specific to the studied area, the recommendations provided here may help raise consciousness about the threat high-density beekeeping may pose to local nature conservation initiatives, especially in areas with sensitive or endangered plant or bee species such as small oceanic islands with high levels of endemism.

In spite of the well-established environmental threats and economic difficulties affecting beekeeping sustainability, the amount of commercial honeybee (Apis mellifera L.) colonies has globally increased by 45% over the last half century, which denotes the rise of beekeeping industry in the era of economic globalization and growing worldwide human population¹. In the meantime, seasonal migration of beehives (also called transhumance) has become a common practice in professional beekeeping where intensive agricultural systems have reduced floral resources and honey yields²,³. Often, these landscapes can no longer sustain apiaries all year round. Instead, beekeepers periodically move their hives into protected natural areas, with potential consequences for the integrity of the native pollinators' interaction networks⁴⁻⁶.

A growing body of literature reports examples of massively introduced honeybees having detrimental effects on the abundance and diversity of local wild bees⁴,⁷,⁸. However, the findings are often inconsistent between different studies because the underlying ecological processes are not well understood and obviously depend on the local situation, including region⁹, habitat type and season⁵,¹⁰⁻¹², or the degree of dietary specialisation in the wild bee population¹³. In the absence of clear evidence of the underlying ecological processes, and without specific guideline or legislation on this issue, protected land managers remain unaware of the potential threat high-density beekeeping poses to their nature conservation efforts. Uncertainty also fuels the debate among bee biologists¹⁴⁻¹⁶.

There are two main ecological processes by which massively introduced managed honeybees may compete with local wild bee populations and affect their spatial patterns of occurrence⁴. In interference competition¹⁷, the superior competitor directly deters others by physical interference. Interference competition, also termed interspecific displacement, is seldom reported to be the main driver of species occurrence in foraging bees¹⁸. Aggressive interaction may also be more often triggered by wild bees than honeybees¹⁹. In exploitative competition, the superior competitor indirectly alters other competitors' fitness or abundance by monopolizing and depleting available resources. Although exploitative competition is expected to be the most common driver of

INRA, UR406 Abeilles et Environnement, Avignon, F-84000, France. Correspondence and requests for materials should be addressed to M.H. (email: mickael.henry@inra.fr)
species occurrence through space and time in the context of competition among foraging bees²⁰, it has hardly ever been evidenced as such because it requires fine assessments of nectar and pollen resource availability. Instead, most studies have investigated side effects of honeybee-induced competition through the study of wild bee flower visitation rates²¹, body size²² or reproductive success¹⁹,²³. Some studies have also reported long term wild bee declines and pollination network deficiencies that are possibly concomitant with the introduction of beekeeping in protected natural areas or other sensitive environments such as small oceanic islands with high levels of endemism²⁴⁻²⁶. In such large-scale correlative studies, it is however difficult to firmly establish the possible contribution of honeybee-induced competition.

In this study, conducted in a protected area exploited by managed honeybees, we first report evidence of wild bee occurrence patterns that are consistent with competition with the introduced honeybees. By “bee occurrence” we mean foraging intensity expressed as flower visitation rate. Second, we show that this pattern is mediated by exploitative competition for floral resources, based on measurements of harvested material in honeybee and wild bee nectar crops and pollen sacs. These estimates of foraging success reflect both interspecies competition between managed honeybees and wild bees inhabiting the protected area, and intraspecific competition among the honeybees themselves. Finally, we provide a rationale for guiding wild pollinator conservation policy in protected areas; we identify operational metrics and thresholds to guide decision-making by conservationists concerned by the impact of beekeeping on wild bee populations.

Sampling was carried out in a 5,700 ha area of Mediterranean scrubland, the Côte Bleue coastal area, southern France, of which 3,400 ha have protection status under Europe's Natura-2000 programme. During the rosemary (*Rosemarinus officinalis* L., Lamiaceae) blooming period (c. 1 month, in March-April), up to 830 honeybee colonies are deposited in 28 apiaries – a density of >14 colonies.km⁻². For comparison, the average national density is 2.5 colonies.km⁻² in France²⁷. During the 2015 and 2016 rosemary honey flows, a total of 180 honeybee and wild bee samples were collected at 60 sampling sites chosen in such a way as to cover the broadest possible range of distances from apiaries (850 ± 830 m (s.d.), range = [10 m, 4000 m]) and of colony density scores (128 ± 78 colonies, range = [12, 287]; see Methods for the density computation). Distance and density scores were analysed separately in relation to wild bee occurrence and foraging success, because land managers may use either colony distance or density metrics to regulate beekeeping, e.g. by setting distance or density thresholds. Wild bee occurrence was assessed as the number of foraging individuals counted on standardized amounts of flowering rosemary shrubs. The nectar and pollen foraging success was assessed in both wild bees and the honeybee. A nectar foraging success index was derived from measurements of harvested nectar in bees' crop, collected with capillary glass tubes after applying a gentle pressure on their abdomen. The pollen foraging success index was obtained from measurements of pollen loads in their pollen collection apparatus.

## Results and Discussion

Overall, our results indicate that the wild bee occurrence pattern in the study area is compatible with competitive exclusion by managed honeybees. As a first striking feature, the study revealed that median honeybee occurrence was 15.3 and 12.9 times greater than that of wild bees in 2015 and 2016, respectively (58.5 vs. 3.8 and 78.8 vs. 6.1 visits per 100 units of rosemary flowering volume, respectively, see Methods). Those ratio cover more than an order of magnitude, mirroring the findings of other studies in mass-flowering areas where there is migratory beekeeping⁵,²⁸,²⁹. Most importantly, the wild bee occurrence rate was affected by the presence of managed bees. It decreased significantly with closeness to an apiary (Table 1, Supplementary Fig. S1) and with increasing honeybee colony density (Supplementary Table S1). In both cases, however, the pattern only emerged after a one-year time lag. For instance, median wild bee occurrence values dropped from 10.4 individuals at 1 km or more from the nearest apiary down to 4.4 individuals at shorter distances. This suggests that shorter distances to, or higher densities of, honeybee colonies cause wild bees to disperse further afield and/or depress their fitness²³, with consequences for the next generation's local occurrence rate. Indeed, as most wild bee species are solitary and univoltine, i.e. with a single generation per year, their current occurrence pattern typically depends on the previous year's nesting conditions, particularly the floral resources available for provisioning nests³⁰,³¹.

Exploitative competition, rather than interference competition, is the most obvious underlying process for this effect on wild bee occurrence. First, despite many thousand observations of foraging bees in the field, aggressive interactions between foraging honeybees and wild bees have seldom been witnessed (M.H. and G.R., personal observations), just as in other studies on the subject¹⁸,¹⁹. Secondly, the large-bodied bees (e.g. *Anthophora, Bombus* or *Xylocopa* species, Supplementary Table S2), which are more likely to physically outcompete honeybees, were in fact more prone to competitive exclusion, as they were found farther away from the apiaries (Table 1, Supplementary Fig. S2). Although rather weak (12% body size decrease around apiaries), this finding is consistent with the exploitative competition hypothesis. Large-bodied bees need more pollen and nectar, and will therefore suffer the effects of exploitative competition before small-bodied ones⁶. Larger bees are also more mobile and can easily disperse away from apiaries to forage and nest in low-competition areas³²,³³. Finally, exploitative competition was directly evidenced by the decrease in nectar foraging success with decreasing distance from apiaries, in both honeybees and wild bees (Table 1, Supplementary Fig. S1). This was revealed by nectar crop measurements carefully controlled for body size variations (see Methods). Exploitative competition for pollen, however, was only evidenced at the intraspecific level. Honeybee pollen foraging success was significantly lower where apiaries were closer (Table 1, Supplementary Fig. S1), while no significant pattern emerged in wild bees. Although speculative at this point, it is possible that the spatial rearrangement of foraging wild bees along the apiary distance gradient has successfully achieved an Ideal Free Distribution (*sensu*³⁴) for pollen harvesting, whereby the lower competitors adjust their distribution so as to balance costs and benefits of foraging and to equalize foraging success³⁵.

Still, further data on pollen and nectar availability were needed to disentangle the effect of intraspecific competition from a possible behavioural trade-off, in honeybees, between distance and harvest. We therefore collected additional field data (see Methods) on nectar and pollen availability in rosemary flowers and performed
<table><caption>**Table 1.** Effect of increasing distance to the nearest apiary on bee occurrence and foraging success. Wild bee occurrence in foraging surveys is better explained by the previous year's apiary distances (inter-annual scale) than by current year distances (annual scale). Analogous statistics for colony density effects are shown in Supplementary Table S1. *All models are LMMs, except wild bee foraging occurrence: Zero-Inflated GLMM (negative-binomial family distribution and log-link function); †Estimates stand for changes per apiary distance unit (km, with log-correction in wild bees); ‡AIC weight of evidence in favour of the apiary distance effect being a better predictor than the colony density effect. The AIC weight ω is shown only when at least one of the two candidate predictors has a significant effect (see Supplementary Table S1 for the colony density effect).</caption><thead><tr><th>Bee occurrence and foraging response variables*</th><th>Sample size (Nb of sites)</th><th>Intercept</th><th>Estimates†</th><th>Statistics</th><th>P-value (effect sign)</th><th>AIC weight (ω)<sup>t</sup></th></tr></thead><tbody><tr><td colspan="7"><strong>Wild bees</strong></td></tr><tr><td>Wild bee foraging occurrence, inter-annual scale (Foraging intensity for 100 flowering volume units)</td><td>180 (60)</td><td>2.63 ± 0.22</td><td>0.46 ± 0.171</td><td>z = 2.71</td><td><strong>0.006 (+)</strong></td><td>78.4%</td></tr><tr><td>Wild bee foraging occurrence, annual scale (Foraging intensity for 100 flowering volume units)</td><td>180 (60)</td><td>2.68 ± 0.20</td><td>0.26 ± 0.13</td><td>z = 1.92</td><td>0.055</td><td></td></tr><tr><td>Mean nectar foraging success (Standardized nectar crop content)</td><td>82 (35)</td><td>23.67 ± 2.61</td><td>8.43 ± 3.85</td><td>t = 2.19</td><td><strong>0.033 (+)</strong></td><td>>99%</td></tr><tr><td>Mean pollen foraging success (Pollen load score)</td><td>78 (39)</td><td>29.16 ± 4.50</td><td>3.10 ± 3.80</td><td>t = 0.82</td><td>0.42</td><td></td></tr><tr><td>Body size, inter-annual scale (Body length, mm)</td><td>220 (44)</td><td>12.84 ± 0.36</td><td>1.02 ± 0.37</td><td>t = 2.75</td><td><strong>0.006 (+)</strong></td><td>>99%</td></tr><tr><td>Body size, annual scale (Body length, mm)</td><td>220 (44)</td><td>12.83 ± 0.26</td><td>1.18 ± 0.37</td><td>t = 3.22</td><td><strong>0.001 (+)</strong></td><td>>99%</td></tr><tr><td colspan="7"><strong>Honeybees</strong></td></tr><tr><td>Mean nectar foraging success (Nectar crop content, μl)</td><td>144 (49)</td><td>4.37 ± 0.54</td><td>1.64 ± 0.32</td><td>t = 5.17</td><td><strong>&lt;0.001 (+)</strong></td><td>>99%</td></tr><tr><td>Mean pollen foraging success (Pollen load score)</td><td>106 (44)</td><td>1.58 ± 0.29</td><td>0.66 ± 0.22</td><td>t = 3.02</td><td><strong>0.004 (+)</strong></td><td>>99%</td></tr></tbody></table>

a confirmatory path analysis³⁶,³⁷ to reconstruct the most plausible chain of causal links with honeybee metrics (Fig. 1). The best causal scenario behind the honeybee intraspecific competition involves apiary distance and colony density as joint drivers of honeybee foraging intensity. Greater honeybee foraging intensity led in turn to lower pollen and nectar availability in rosemary flowers. This result validates the hypothesis of intraspecific exploitative competition among honeybees for pollen and nectar where colony density is high. Distance to nearest apiary was still the key driver of competition in the studied scrubland (see AIC weights in Table 1) and should therefore take precedence over density score as a criterion for regulating beekeeping.

Practical thresholds emerged for apiary distance vs. competition effect size, paving the way for pragmatic decision-making by managers of protected land who are concerned with wild bee conservation (Fig. 2). We have expressed competition effect size (%) as the percentage changes in bee foraging success and occurrence with distance from nearest apiary. Effect sizes were recursively recomputed by adjusting the distance limit between closer and more distant sites, from the first (150 m) to the third (1,200 m) quartiles of distances covered by the study (Fig. 2). We detected marked effect size peaks which showed that competition was highest within a certain distance from an apiary and relaxed beyond that distance. The competition effect was operative (i) 600 m away from the apiary with a 50% decrease in wild bee nectar foraging success, (ii) 900 m away with a 55% decrease in wild bee occurrence, (iii) 1,100 m away with a 44% decrease in honeybee nectar foraging success and (iv) 1,200 m away or more with a 36% decrease in honeybee pollen foraging success (Fig. 2). The honeybee distance thresholds fell within the median foraging ranges of 1–2 km usually reported in the literature³⁸⁻⁴⁰. It also fits distances at which detrimental effects of apiaries have been reported on native bumblebee foraging behaviour (until 1,200 m in⁴¹).

We argue that protected land managers could use honeybee intraspecific competition as a criterion to guide their regulation of mass-flowering resource exploitation by beekeepers. Competition effect size measured in wild bees peaked at shorter distances from apiaries (600–900 m, Fig. 2) than honeybee intraspecific competition (Nectar: 1.1 km, pollen: >1.2 km, Fig. 2). Therefore, any buffer distance rule derived from the latter metrics will be more conservative regarding wild bee protection. It will also inform managers about the carrying capacity of the area and help optimize beekeeping honey flows. For instance, in the Mediterranean scrubland area we studied, with an average apiary size of 30.9 ± 21.8 (s.d.) colonies, intraspecific competition for nectar is theoretically relaxed beyond the effect size peak distance of 1.1 km. Thus the area theoretically reaches saturation with a distance of 1.1 × 2 = 2.2 km between neighbouring apiaries, i.e. ~1 apiary per 3.8 km² (that is 0.26 apiary/km²). Managers wishing to allocate half an area to wild bee conservation under relaxed competition will then need an apiary load twice as low, i.e. ~1 apiary per 7.6km² (that is 0.13 apiary/km² or 3.1 km spacing between apiaries).

Overall, those distance-based recommendations do not return honeybee colony densities that are fundamentally different from others found in literature. Considering our average apiary size of about 30 colonies, the “saturation” (0.26 apiary/km²) and “half-saturation” (0.13 apiary/km²) apiary loads are equivalent to 7.8 and 3.9 colonies/km², respectively. Interestingly, the latter value, which stands for a 50% land-sharing between honey production and wild bee conservation, is only slightly more permissive than previous recommendations (3.1 colonies/km² in⁴² or 3.5 colonies/km² in⁶). But most importantly, the distance-based thresholds we suggest herein might help inform the recent debate on the effectiveness of beekeeping regulation in protected areas with respect to honeybee foraging range¹⁴,¹⁵. Given the high mobility of foraging honeybees, protected areas may be theoretically exploited from off-site apiaries located up to 10 km away, potentially making local regulation rules inefficient¹⁴. On the other hand, such extreme foraging ranges are thought to be uncommon, and competition effects is rather expected to range around the mean honeybee foraging distances, i.e. about 1 km away from apiaries¹⁵. Here, we provide empirical support to the view that local beekeeping regulation will indeed benefit wild bee
<figure><img src="image_55.png" /><figcaption>**Figure 1.** Reconstruction of the most plausible causal chain scenario behind honeybee intraspecific competition. Arrows indicate causal links among variables. P-values and signs stand for link significance and effect direction. In a first tentative path model scenario, apiary distance was positioned as the only proximal driver of honeybee occurrence (foraging intensity) in sampling sites, affecting in turn nectar and pollen availability in flowers. However, a significant effect of colony density scores, conditional on distance, was detected in the causal chain (dashed arrow, see Supplementary Table S5 for detailed statistics), making the conditional independence requirements close to being rejected (d-separation test of deviation from conditional independences, $C = 19.22$, $df = 12$, $P = 0.083$). A joint distance-density effect path model was also computed, whereby apiary distance and colony density jointly influenced honeybee foraging intensity (see Methods for the computation of detrended colony density). The joint distance-density effect scenario satisfactorily met the conditional independence requirements ($C = 8.75$, $df = 10$, $P = 0.56$) and was far better supported than the distance scenario, given the much greater AICc weight of evidence ($\omega > 99\%$ in favour of the joint effect scenario, Supplementary Table S5).</figcaption></figure>

conservation in protected areas. We based our reasoning on functional aspects of foraging ecology, rather than on any rough guess of the actual honeybee foraging range.

Managers of protected land should apply the distance-based thresholds only where natural mass-flowering resources are over-exploited, with a view to gradually reducing existing beekeeping pressure. Managed honeybees should not be introduced into pristine areas or areas with sensitive or endangered plant or bee species⁴, such as in small oceanic islands with high levels of endemism²⁴⁻²⁶,⁴³. In continental environments, managers of protected areas may also constrain the location of apiaries by applying the distance-based regulation threshold around habitats or microsites of special conservation interest that are identified to host threatened or emblematic plant or pollinator species. Thresholds should be carefully re-evaluated for each situation concerned. In addition to threshold-based regulation, land managers could envision periodic break years to temporarily halt competition disturbance regime and boost resilience in wild bee populations.

Regulating colony density will also benefit honeybees themselves. Beekeepers may not perceive a substantial honey yield decrease under high-density management, owing to the colonies' internal regulatory processes. But exploitative competition may constrain honeybee foragers' lifespan and trigger a cascade of problems that eventually leads to colony weakening or collapse later in the season, long after the migration period has ended³.

The honeybee tends to take precedence over wild bees as a target species in conservation programs because for policy makers and land managers the honeybee is an emblematic pollinator species⁷. Some European countries, including France, have agro-environmental schemes that envision subsidising beekeepers to set up apiaries in natural areas. We believe these agro-environmental schemes should, on the contrary, reward beekeepers who make joint efforts with farmers to maintain their apiaries in agro-ecosystems all year round. Protected land managers and beekeepers should realise that mass-flowering resources in natural areas are shared resource systems. If beekeepers exploit them independently, according to their own self-interest, they have a high chance of acting contrary to the common good of all users, others beekeepers and wild bees, by depleting or spoiling that resource through their cumulative actions. This is called the tragedy of the commons⁴⁴.
<figure><img src="image_56.png" /><figcaption>**Figure 2.** Beekeeping-induced competition as a function of distance from nearest apiary in (a) wild bees and (b) honeybees. The competition effect size (percentage decrease in foraging intensity and foraging success close to apiaries) was recomputed by varying the distance between closer and farther sites, from the first (150 m) to the third (1,200 m) quartiles of distances covered in the study. Trends are depicted by LOESS local regression fits and 95% confidence envelopes (solid lines and coloured areas respectively). For each competition metric, the panels on the right show the distribution of values (quartile boxes) for sites located closer to vs. farther away from the nearest apiary. The AIC weight $\omega$ gives the probability that competition is better accounted for by a two-step threshold effect model (closer-$v$s. - farther binary distance variable) rather than a progressive effect model (continuous distance variable). Thresholds emerged ($\omega >> 50\%$) at 600 m for wild bee nectar foraging success, at 900 m for wild bee occurrence and at 1100 m for honeybee nectar foraging success. The progressive effect model was better supported for honeybee pollen foraging success ($\omega << 50\%$) as effect size did not peak but steadily increased until 1200 m. Effect sizes for wild bee occurrence were based on the inter-annual scale (Table 1). See Supplementary Table S6 for detailed data and sample sizes at each distance class and AIC model selection statistics.</figcaption></figure>

## Methods

The hypotheses of altered wild bee occurrence and depressed foraging success around apiaries were investigated by means of honeybee and wild bee occurrence surveys and foraging success assessment in a range of sampling sites located at different distances from apiaries in a protected natural area.

**Study area and sampling design.** The Côte Bleue study area is a 5.700-ha protected Mediterranean scrub-land, dominated by Kermes oak *Quercus coccifera* and Rosemary *Rosmarinus officinalis*. The area is managed by the French coastal protection agency (Conservatoire National du Littoral) in partnership with the French National Forest Office (Office national des forêts, ONF). Each year, ONF allocates 28 sites in the Côte Bleue to accredited beekeepers so that the location and size of each apiary is thoroughly registered. Those apiary
registration data were systematically ground-truthed at each study year. Beekeepers typically set up apiaries for rosemary honey production during ca. 1 month early in the season, in between Mars and April.

During the 2015 and 2016 rosemary blooms, a total of 180 wild bee field samples were carried out in 60 structurally similar scrubland sampling sites chosen so as to cover broad gradients of beekeeping intensity in terms of nearest apiary distance and size. The sampling design was however constrained by inter-annual fluctuations in the presence of apiaries. We varied accordingly sampling site locations between years (Supplementary Fig. S3) in order to maintain a balanced site allocation across beekeeping gradients. Data analysis was therefore based on a generalised mixed model framework to account for the resulting spatial and temporal nested design (see *Data Analysis*).

Sampling site selection was first driven by apiary location, and then by local rosemary flowering cover. The area was subdivided into four main contiguous sectors with different access trails. In each sector, a first series of sites were chosen for their proximity to the main apiaries, i.e. a few tens of meters away, corresponding to the smallest foraging range usually reported in literature for wild bees⁴⁵. Then, a second series of sites were located as far as possible from apiaries (1 to 4 km). Finally, additional sites were chosen at intermediary positions on the way, so that the overall design covered a wide range of distances (850 m ± 830 m (s.d.)) with balanced site numbers among sectors. Sites were defined as an area of 50 m in diameter with a widespread, uniform, rosemary distribution – totalling on average 10% to 15% of soil cover. Rosemary is by far the dominant flowering resource throughout the scrubland, and particularly during the beekeeping migration period when secondary resources (mainly *Thymus vulgaris*, *Cistus albidus*, *Cistus salviifolius*, *Helianthemum marifolium*, *Reseda phyteuma*), have not bloomed yet. We still controlled for food resource conditions by avoiding sites with conspicuous or unusually dense patches of secondary resources. Accordingly, only 3.5% of sampled bees were captured on the secondary floral resources, which we considered to be too restrictive to deserve specific analyses.

All sampling sites within a given sector were typically processed during the same half-day. No sample was performed under rainy, windy or cold (<12°C) weather. All sectors were therefore covered in 2-days sessions. The rosemary blooming period was covered with a total of six sessions each year. Depending on local phenological fluctuations, individual sites were visited on average 3.0 ± 1.6 (s.d.) times. During a given session, we randomized the sector visit order so that sample locations were independent from time of the day (varying from 9:00 to 18:00, solar time). Finally, we assumed that our sampling design was poorly affected by possibly undetected large apiaries located outside of the boundaries of the protected area. The area is bordering the Mediterranean sea on the south part, and is surrounded by a large (155 km²) water body, the Pond of Berre, on its West and Northern parts. The immediately adjacent lands on the North and East parts of the area are densely urbanised and are not suitable for professional apiaries.

**Wild bee occurrence.** Field samples consisted in the joint assessment of the local rosemary floral resource availability and of the wild bee occurrence on those flowers (foraging intensity expressed as flower visitation rate). Fifteen flowering rosemary shrubs were carefully inspected, starting from the largest one in the site, and then moving step by step to the nearest neighbours. Contiguous shrubs with undistinguishable, coalescing crowns were treated as a single individual. Shrub flowering volume estimates were derived from the three shrub crown dimensions (length, width and height) rounded to the nearest 50 cm, which we defined as the smallest tractable volume dimension unit. If necessary, shrubs with irregular shapes were sized in two or more steps. Shrub flowering volumes were further weighed by a floribondity multiplicative coefficient, reflecting the percentage (±10%) of open flowers relatively to the expected maximal number of flowers branches may actually bear (up to 80 flowers per 20 cm). For instance, if a shrub is visually estimated to be at 50% of its flowering potential, its crown volume is then corrected to half its actual size. Finally, we tallied foraging honeybees and wild bees while inspecting shrubs, and computed the corresponding numbers of bees per unit of flowering rosemary volume to serve as a foraging intensity measurement. The same two observers performed all those field estimates. Preliminary blind comparisons led to consistent and highly correlated estimates between observers. Nevertheless, an observer was kept unchanged throughout each sampling session to avoid biases. The bee survey routine included net captures for assessing wild bee and honeybee individual foraging success.

**Individual Foraging success.** Foraging success was assessed in female wild bees and honeybees using nectar crop content and pollen load measurements. A variable transect walk (*sensu*¹⁶) was performed to collect all bees with nets. observers walked at moderate speed among flowering rosemary shrubs and collected bees during a minimum of 20 minutes.person, which was usually sufficient to capture at least 10 foraging honeybees. We did not constrain samples with a minimum number of wild bee captures because those were much less abundant. Captured individuals were maintained in a cooler at about 4°C to slow down their metabolism before being processed. Throughout the study, we favoured non-invasive methods. Once processed and identified to family or genus, wild bees were released right on the capture site. Only few specimens per morphotype were collected for subsequent identification to species. The resulting preliminary checklist for the study area is provided in Supplementary Table S2.

Nectar foraging success was assessed by measuring the nectar volume stored in their crop stomach at the time of capture, also termed field nectar load⁴⁷,⁴⁸. Bees were first narcotized for few seconds with CO₂ (ProFlora U500 Cylinder, JBL GmbH & co, Neuhofen, Germany). A gentle dorso-ventral pressure was then applied on their abdomen until their crop nectar content was regurgitated. The extracted nectar was capillary-collected (10 µl micro-pipettes Ringcaps, Hirshmann Laborgeräte GmbH & co, Eberstadt, Germany) for volume measurement with an estimated ± 0.05 µl resolution. We then used a refractometer (REF108, Index Instruments Ltd., Cambridgeshire, England) to ascertain the presence of high concentrations of sugar. On very rare occasions, crop content was identified to be mostly water, with null or low (<5%) sugar concentrations. Bee water foragers were discarded from the nectar foraging databases. Finally, we measured (to the nearest mm) wild bee body length, from head to