Risk of potential pesticide use to honeybee and bumblebee survival and distribution: A country‐wide analysis for The Netherlands

Bees play an important role in natural ecosystems and the world's food supply. In the past decades, bee abundance and diversity have declined globally. This has resulted in decreased pollination services for natural ecosystems and the agricultural sector at the field scale. One of the causes of the decline in bee abundance and diversity is the use of pesticides. Linking pesticide use, land use and bee presence could provide crucial insights into areas, and pesticides that pose a significant threat to the abundance and diversity of bees. Obtaining actual figures of farmer pesticide use is rarely possible. Therefore, we designed a method to study the effects of potential pesticide use on the survival and distribution of honeybees and bumblebees.


| INTRODUC TI ON
Bees are important organisms for ecosystems, food production and the economy. They serve as pollinators, playing a role in the abundance, distribution and diversity of flowering plants. Some of these plants are cultivated by humans. Worldwide, 35% of all cultivated crops depend on pollinators (Klein et al., 2007). With a contribution of 153 billion euros, pollinators represent a share of 9.5% of the total value of the global food production (Gallai, Salles, Settele, & Vaissière, 2009). Both honeybees and wild bees play an important role in plant pollination with wild pollinator communities often outperforming managed honeybees in pollination success (Garibaldi et al., 2013). While beekeepers can compensate actively for yearly honeybee losses, this type of population recovery is not possible for wild bees. Flowering plants that depend on pollination by bees could become the subject of local or even global extinction as the abundance and diversity of bees decline. This process has been observed in Britain and the Netherlands (Biesmeijer et al., 2006). The declines of pollinator abundance and diversity have been reported to occur around the globe (Potts et al., 2010(Potts et al., , 2016Vanbergen,2013) having pesticide use, habitat loss, parasite load and changes in climatic conditions as some of the main contributing factors for the declines.
The risk of pesticides to bee health has been studied extensively for individual active ingredients and various active ingredients have been found to be hazardous for bees. For instance, imidacloprid (neonicotinoid) has been shown to affect both honeybees and wild bees (Goulson, 2013;Rundlöf et al., 2015;Woodcock et al., 2016).
Although the effects of field-level active ingredient concentrations are often sublethal, chronic exposure has been shown to negatively affect foraging behaviour in bumblebees (Gill, Ramos-Rodriguez, & Raine, 2012) and honeybees (Schneider, Tautz, Grünewald, & Fuchs, 2012). More recently, field-realistic exposure to neonicotinoids was shown to reduce honeybee health in corn-growing regions (Tsvetkov et al., 2017). A parallel study by Woodcock et al. (2017) found negative effects of neonicotinoid seed coatings on honeybee health in winter-sown oilseed rape, but also found the reproduction of wild bees to be negatively correlated with neonicotinoid residues in wild bee nests. Importantly, the findings of the latter two studies were based on field-based experiments.
Honeybees are considered highly sensitive to pesticide pressure (Porrini et al., 2003). The number of pesticide products in beeswax has been found to have a negative effect on honeybee colony survival (Traynor et al., 2016). The median lethal dose (LD 50 ) values for this bee species are known for a broad range of chemicals applied in crops. A recent review found high variability in LD 50 values among other bee species for multiple active ingredients (Arena & Sgolastra, 2014), thus assessing the risk of pesticide use for each bee species separately seems critical to investigate the species-specific impact. A recent study showed that pesticide risk is buffered by surrounding natural habitat (Park, Blitzer, Gibbs, Losey, & Danforth, 2015), indicating that a habitat specific approach is needed. Pesticide risk assessments for bees can be done using a hazard quotient (HQ , EPPO, 2010). It has been shown that by using the treated area combined with application rates of active ingredients and the LD 50 values, the risk of bee mortality can be predicted (Mineau et al., 2008).
The loss of natural habitat has been identified as one of the causes for the decline in bee abundance and diversity (Winfree, Aguilar, Vázquez, LeBuhn, & Aizen, 2009). However, the magnitude of the effect was not found to be large, and only significant in areas with very little natural habitat remaining. A review by Potts et al. (2010) confirmed this finding, but also noted that some studies found the opposite effect for some bee species. For example, honeybee (Apis meliffera) abundance was higher close to urban habitats on a European scale (Carré et al., 2009). This abundance is logically caused by activities of beekeepers, making a comparison between the abundance of honeybees to wild bee species difficult. However, the study by Carré et al. also found positive effects of urban habitats on wild bee abundance, noting that this could be the result of suitable nesting habitats in bare grounds for certain wild bee species.
A more recent study underlined the importance of the historical structure of habitats, showing that bumblebees were more susceptible to habitat fragmentation in areas with high historical species spillover potential (i.e., the potential to disperse from one land class to another) and additionally, the edge density between managed and natural areas was found to have a positive effect on bumblebee species richness (Aguirre-Gutiérrez et al., 2015). Given the relatively fragmented state of (semi-)natural areas in urbanized regions, it is pivotal to identify vegetation structures that positively affect bee abundance and diversity to allow for more effective conservation measures. This is of importance as different vegetation with varied structures may render a diverse set of feeding and nesting resources for pollinators (Proctor, Nol, Burke, & Crins, 2012). One of the technologies that is increasingly used to study and map vegetation structure is light detection and ranging (LiDAR). This method allows for three-dimensional mapping of vegetation structure in high resolution (Davies & Asner, 2014;Jetz et al., 2016;Simonson, Allen, & Coomes, 2014).
Changes in climatic conditions may have consequences for the phenology of plants and pollinators, which could result in the desynchronization of the two (Goulson, Nicholls, Botias, & Rotheray, 2015;Pyke, Thomson, Inouye, & Miller, 2016). The changing climate could also lead to range shifts of plants and pollinators, possibly leading to spatial mismatches (Pyke et al., 2016;Schweiger et al., 2010). Predicting future pollinator distributions is crucial to assess the consequences of changes in climatic conditions. In the past century, the average temperature increased by 1.7°C in the Netherlands (Ligtvoet, vanMinnen, & Franken, 2013). A recent study by Aguirre-Gutiérrez, Kissling, et al. (2017) found an increasing variable importance of temperature in predicting bee distributions in this country over the studied 1951-2014 period. This underlines the possibility of severe consequences for bee abundance and diversity as climatic conditions change.
In this study, we developed potential pesticide risk maps based on the number of allowed active ingredients on all agricultural parcels in the Netherlands. The maps were analysed to test the effects of potential pesticide use on honeybee colony survival and bumblebee presence. Given the harmful effects of various pesticides on bees reported in recent studies, a negative effect of potential pesticide risk on honeybee colony survival and bumblebee presence was hypothesized. Consequently, honeybee colony losses were expected to be higher in regions that pose a relatively higher potential pesticide risk compared with regions that pose a relatively lower risk.
For bumblebees, the presence in these high-risk areas was expected to be lower compared with low-risk areas. We further analysed the effects that land use, vegetation structure and climate may have on honeybee colony survival and bumblebee distributions. We hypothesized that these factors play an important role in shaping survival and distribution of honey and bumblebees as they may facilitate or hamper access to feeding and nesting resources.

| ME THODS
We modelled the potential pesticide risk for honeybees and bumble-

| Potential pesticide use data
The potential pesticide use was based on the allowed pesticide products for each crop. A list of all allowed products for 2017 was obtained from the College voor de Toelating van Gewasbeschermingsmiddelen (CTGB; www.ctgb.nl), the pesticide regulatory organ in the Netherlands. The list consisted of 2,556 files with information on the crops and pesticide application rates associated with it. If a maximum application amount per year or per cultivation cycle was available, it was used as the maximum allowed application quantity. If this quantity was not available, the regular advised amount was used as the maximum allowed application quantity. The average application amount of the crops F I G U R E 1 Schematic illustration of the workflow. Data concerning allowed pesticide use associated with each crop were linked to parcel data including crop specifications (top). On every agricultural parcel, the number of allowed risky active ingredients was counted, resulting in the risk rasters (bottom). The risk value associated with a honeybee colony or bumblebee presence point was calculated by taking the average risk value of the raster cells within a 3-km buffer identified per pesticide product was used to calculate the hazard quotient. The final dataset contained the maximum allowed amount of pesticide product use per hectare for every crop and the concentration of the product.

| Active ingredient sensitivity and crops attracting bees
The median lethal dose values for honeybees and bumblebees were obtained from Sanchez-Bayo and Goka (2014; Appendix S7). A list of crops attracting honeybees was obtained from the CTGB (Appendix S6). Such a list was not available for bumblebees, and we created it using several sources from Koppert Biological Systems (www. koppe rt.com), Klein et al. (2007)

| Pesticide risk map model
The risk assessment for active ingredients was done using a hazard quotient (HQ). The method was based on an insecticide loading assessment method presented in Gillespie et al. (2017). The hazard quotient was calculated as follows: For every crop that attracts the specific bee group, the allowed pesticide use was analysed. The chemicals posing a possible risk to bees (HQ > 50; EPPO, 2010) were counted on every agricultural parcel.
If no median lethal dose value was available for a certain chemical and bee group, this chemical was not included in the assessment. If a parcel was specified as "organic", the risk value of the parcel would be zero.
However, if on that organic parcel a crop was cultivated in which the application of the active ingredient spinosad was allowed, a risk value of one was assigned to the parcel, spinosad is the only risky chemical that is allowed in organic farming. This workflow described the number of risky chemicals that may be used per agricultural parcel, and based on it, we created a spatially explicit raster map layer, with each raster cell holding the potential risk value of the corresponding cultivated crop (Figure 1) TA B L E 1 Description of environmental variables used in the initial GLM(M) models risky chemicals that could be applied in the raster cell. We created the "risk" raster layer at a fine spatial resolution of 100 × 100 m. LiDAR-derived vegetation structure information was obtained from Aguirre-Gutiérrez, WallisDeVries, et al. (2017)

| Bioclimatic variables
We obtained precipitation and temperature data for the period 2000 to 2015 from the Dutch meteorological institute, KNMI (https ://data.knmi.nl/datasets). These data were used to generate monthly bioclimatic variables (Fick & Hijmans, 2017) representing ecologically meaningful descriptions of precipitation and temperature for the Netherlands with a spatial resolution of 100 × 100 m (Table 1).

| Statistical analysis
To investigate the contribution of the potential pesticide risk to honeybee winter colony losses (survival/death), we applied generalized linear mixed models (GLMM with a binomial error structure). In the GLMM's, we included the beekeeper as a random factor to account for the repeated measurements with beekeepers (most beekeepers had multiple colonies). As explanatory variables, we included the potential pesticide risk, land use, vegetation structure and the bioclimatic variables. In total, we selected 21 explanatory variables ( We first split the 21 explanatory variables into two classes: land use and other variables. We then combined the most important variables that resulted from the modelling selection protocol for each of the two classes. These most important variables were then used for another modelling selection run. Finally, the most parsimonious models were selected to calculate the average final model. We modelled the bumblebee distributions using generalized linear models (GLM) with a data train/test ratio of 0.75/0.25. As true absences were not available, we generated 1,000 pseudo-absences and created 100 models. The best models were selected based on the highest AUC (area under the curve; Hanley & McNeil, 1982). The used 24 predictor variables (Table 1) had a VIF < 5 (Appendix S4). All analyses were implemented in R version 3.4.0 using the "lme4" and "MuMIn" packages.

| Pesticide risk maps
The analysis of the potential pesticide risk model using the allowed pesticide data resulted in the identification of 46 crops with potential risk for honeybees and 23 crops with potential risk for bumblebees (Table 2; Appendices S2 and S3). This translated into a larger area with potential risk for honeybees compared with bumblebees.
The active ingredients responsible for most hectares with potential risk were lambda-cyhalothrin for honeybees and abamectin for bumblebees ( Figure 2). Most parcels showed a change in risk value from 2015 to 2016, both for honeybees and bumblebees (Figure 3; Appendix S1). Median potential risk value 1 1

TA B L E 2 Summary of Appendices S2 and S3; Crops and potential risk
Note: Potential risk was assigned if the specific crop attracted honeybees or bumblebees and the hazard quotient associated with the chemical and the allowed application rate on the crop exceeded 50.

| Honeybee colony survival
The average model for 2015 showed negative effects of canopy density between 0.5 and 2 m, greenhouses, annual mean temperature and mean temperature of the wettest quarter on honeybee colony survival (Table 3)

| Bumblebee distribution
The GLM models showed negative effects of the potential pesticide risk factors for 2015 and 2016, but non-significant (p = .07 and p = .12, respectively, Table 4

| D ISCUSS I ON
Recent studies investigating the impact of agricultural pesticide use have shown a negative effect on honeybee health in field-realistic settings (Tsvetkov et al., 2017;Woodcock et al., 2017). Since spatial distributions of bees vary per species, it is of critical importance to study pesticide risk to bees in a spatial and species-specific manner by incorporating knowledge concerning local agricultural circumstances.
In this study, we estimated and mapped the potential pesticide risk to honeybees and bumblebees. This was done by combining data concerning the allowed pesticide use per crop and pesticide product, the Dutch parcel registry specifying the location of all agricultural parcels and the associated cultivated crops, and the honeybee or bumblebee specific sensitivities to pesticides. We analysed if and to what extent the risk of this potential pesticide use, land use, vegetation structure and climate may drive survival of honeybee colonies and the distribution of bumblebees in the Netherlands. We find that while the variable importance of the modelled potential pesticide risk seems to be low, the identified non-significant negative effects should be studied further using knowledge concerning actual pesticide use data.
Additionally, including more potential risk maps of other years than 2015 and 2016 could be vital in the understanding of the development of pesticide pressure. Finally, we find a strong indication that pesticide risk from greenhouses should be included in such future studies.

| Potential pesticide risk for honeybees and bumblebees
For honeybees, more active pesticide ingredients were identified as a potential risk compared with bumblebees. This could have been caused by several factors. Firstly, bumblebees were less sensitive to certain active ingredients than honeybees. For example, the median F I G U R E 2 Area with potential risk per active ingredient for honeybees and bumblebees in 2015 and 2016. The numbers above the bars represent the number of crops in which the specific risky active ingredients can be used, separated for honeybees (pink) and bumblebees (red) lethal dose by contact for deltamethrin was 0.024 µg for honeybees and 0.28 µg for bumblebees. This resulted in deltamethrin not having an HQ value higher than 50 on any crop for bumblebees. The second factor was the absence of median lethal dose data for bumblebees. For example, the value for esfenvalerate was not available for bumblebees, while this chemical was hazardous to honeybees (0.026 µg). Thirdly, some crops were identified as honeybee-attracting crops, but not as bumblebee-attracting crops (Appendix S6). The crop groups with most potential risk hectares were fruits (for honeybees and bumblebees) and potatoes (for honeybees). The active ingredient group representing most risk was pyrethroids for honeybees (e.g., deltamethrin and lambda-cyhalothrin). For bumblebees, the potential risk was represented by only four chemicals, of which abamectin caused the largest risk area. Imidacloprid (neonicotinoid) can be used in the cultivation of both apples and pears (spray application) and was identified as a potential risk for both honeybees and bumblebees.
Unfortunately, the automated identification of the allowed application rates from the product manuals for the active ingredient thiamethoxam failed. Thiamethoxam is allowed in potatoes, and its HQ value would have exceeded 50 for honeybees. This means that the risk value assigned to potatoes was underestimated with a potential risk value of one.
The implemented full GLMM models using the honeybee colony survival data indicated non-significant negative effects of potential pesticide risk. Showing an overall low variable importance in the full GLMM models, the potential pesticide risk factors were not included in the average models. Therefore, this study cannot confirm the findings of Traynor et al. (2016), where a negative effect of the number of pesticide products in beeswax on honeybee colony survival was found. For the bumblebee presence data, the implemented GLM models showed non-significant negative effects of potential pesticide risk.
Due to the assumptions made in this study, the results should be interpreted with caution. Several important factors have been simplified during the potential pesticide risk mapping process. These factors will be described below.
Firstly, the mode of application and information concerning whether a pesticide is systemic or not was not included in the potential pesticide risk model. Secondly, the risk of exposure to pesticides will change over time. The number of allowed pesticide products will change every year, for the CTGB may alter the product regulations. Additionally, chemicals will degrade and leach over time, every chemical and environment will have specific rates at which these processes occur. Moreover, we ignored regulations concerning pesticide application timings.
The actual pesticide use will depend on several conditions, such as the crop type that is cultivated or the presence or absence of pests that need to be controlled. Therefore, the usage is subjected to local changes. However, data concerning local actual pesticide use are not (publicly) available. Moreover, the data that are publicly available are very limited, making the linkage of pesticide use of individual active ingredients to crop types virtually impossible.
The foraging of bees on flowering weeds in crops that do not attract bees was not included in this study. We assumed that bees only forage on parcels containing bee attracting crops (Appendix S6). However, a recent study found a risk of exposure to pesticides through weeds present in crops unattractive to bees (Simon-Delso, Furthermore, all flowering crops were assumed equal (e.g., same size, nutritional value), crop attractiveness was assumed to be binary (either attractive or non-attractive for bees), the varying routes of exposure ignored (the median lethal dosages by contact were used to calculate the hazard quotients), the toxicity of active ingredient combinations was assumed to be additive and the food intake rates of honeybees and bumblebees were assumed to be equal.
The provided risk maps can be the first step towards a risk assessment based on local pesticide pressure. Increasing the accuracy of the approach used in this study could be accomplished in various ways. Firstly, the inclusion of a risk assessment incorporating seed treatments, for this application method has been found to be harmful in several cases, for example with imidacloprid (Girolami et al., 2009). This could be achieved by adding a toxicity exposure ratio to the mapping approach and using it with soil and seed treatments while retaining the hazard quotient method for spray applications (EPPO, 2010). Adding this method would also introduce a speciesspecific estimation of the amount of food ingested by bees, which could have an impact on the risk assessment. Secondly, the potential pesticide use data should be specified for every year separately, following the changes in the CTGB's allowed pesticides database.
Thirdly, as noted previously, the inclusion of pesticide risk posed by greenhouses, for the runoff water from this source has been found to contain pesticides (Haarstad, Bavor, & Roseth, 2012;Tamis, van't Zelfde, & Vijver, 2016). Fourthly, incorporating knowledge concerning the actual pesticide use in the Netherlands is crucial to grasp the local variations in the risk posed by pesticides. Data concerning average actual pesticide use by farmers per province or municipality could give a rough estimation of these local variations. Ultimately, a public system registering the actual used amounts of pesticides per agricultural parcel would be a significant improvement. Such a public system should be justified by the benefits for human and animal health that could result from improved risk assessments and mappings. All the previously mentioned assessment improvements focus on one single niche of pesticide use; the agricultural sector. Including an estimation of the pesticide use of households and the government should improve the risk assessment, especially in urban areas.

| Effects of land use and vegetation structure on honeybees and bumblebees
One

| Effects of climate on the survival and distribution of honeybees and bumblebees
In both the honeybee average models and the bumblebee models, climatic factors were of relatively high importance. For honeybees, the annual mean temperature showed negative effects, which can be interpreted to be in line with a recent study of Switanek, Crailsheim, Truhetz, and Brodschneider (2017). This study found an increase in honeybee winter colony mortality when weather conditions in the preceding year were warmer and drier. For the honeybee data of 2015, the most important factor was the mean temperature of the wettest quarter, showing a negative effect on honeybee colony survival. This again seems to be in line with the latter study. In the bumblebee models, the temperature annual range factor had a significant negative effect on bumblebee presence, which could be an indication for a sensitivity to temperature extremes.

| Concluding remarks
In this study, the potential risk of pesticide use to honeybees and bumblebees was mapped and analysed. The presented potential pesticide risk maps could aid in the conservation of wild bee species and the prevention of honeybee colony losses. The maps could help identify regions with relatively high pesticide pressure in a speciesspecific manner, enabling conservation actions on a local scale. This could result in a lower local pesticide pressure for the bee species in question, while minimizing economic damages, since enforced pesticide regulations could be tailored to local high-risk areas. Such actions could help restore threatened bee species and lower honeybee colony losses, which would benefit pollinator dependent farmers and plant species by the increase of pollination services.

ACK N OWLED G EM ENTS
We thank the pollinator group of Naturalis Biodiversity Center for their support and comments during the study. Additionally, we thank Maarten van 't Zelfde for his help and support and Menno Reemer for the collaboration with EIS-Nederland.

DATA AVA I L A B I L I T Y S TAT E M E N T
The vegetation structure data are available at the Dryad Digital Repository: https ://doi.org/10.5061/dryad.fg1r11d.
F I G U R E 5 Variable importance in the GLM models for the 2015 and 2016 bumblebee data