Impact of pesticide use on wild bee distributions across the United States

Abstract

The decline of many wild bee species has major consequences for pollination in natural and agro-ecosystems. One hypothesized cause of the declines is pesticide use; neonicotinoids and pyrethroids in particular have been shown to have pernicious effects in laboratory and field experiments, and have been linked to population declines in a few focal species. We used aggregated museum records, ecological surveys and community science data from across the contiguous United States, including 178,589 unique observations from 1,081 bee species (33% of species with records in the United States) across six families, to model species occupancy from 1995 to 2015 with linked land use data. While there are numerous causes of bee declines, we discovered that the negative effects of pesticides are widespread; the increase in neonicotinoid and pyrethroid use is a major driver of changes in occupancy across hundreds of wild bee species. In some groups, high pesticide use contributes to a 43.3% decrease in the probability that a species occurs at a site. These results suggest that mechanisms that reduce pesticide use (such as integrative pest management) can potentially facilitate pollination conservation.

Main

There are widespread reports of bee declines in Europe and North America, but the status of most species is poorly known^1,2,3,4,5,6. Insect pollination, largely from wild and managed bees, benefits ~75% of crop species worldwide and 88% of flowering plant species¹. Further, the majority of crop pollination is provided by wild pollinators worldwide, and wild pollinators can enhance yields regardless of managed bee abundances^7,8. Overall, this suggests that the decline of wild pollinators will have strong detrimental effects on pollination services, with both economic and ecological consequences.

The major drivers of wild bee declines include climate change⁵, land use change and habitat loss⁹, disease and pathogens¹⁰, dietary stress, and pesticide use^9,11. Many of these factors are primarily associated with agricultural intensification (Extended Data Fig. 1)¹². First, agricultural intensification generally reduces the diversity of floral and nesting resources available¹³. This is particularly the case in crop monocultures that do not provide resources to pollinators, such as grain monocultures¹⁴. Second, agricultural intensification can increase the exposure of wild bees to combinations of pesticides¹⁵. Among these pesticides, two notable classes are neonicotinoids and pyrethroids. The usage of these compounds is widespread across the United States, with neonicotinoid usage having rapidly increased since their introduction in the mid-1990s (Supplementary Figs. 1–3), while the usage of other classes of insecticides declined over the same time period (Extended Data Fig. 2)¹⁶. Evidence from both laboratory and field experiments has demonstrated that both types of pesticide are harmful to individual bees; thus, exposure to the high quantities typically used in agriculture could potentially cause population declines^17,18,19. Neonicotinoids are a class of neuro-active insecticides that target the central nervous system²⁰. They are either sprayed or applied as soil drenches or seed treatments. When used as a seed treatment, they are systemic, expressed throughout plant tissues including pollen and nectar. The effects of neonicotinoids are typically sublethal and chronic, and therefore difficult to detect under typical regulatory studies²¹. Because of their chronic and sublethal effects, the European Union banned neonicotinoids in 2018 (with a moratorium since 2013)²². Pyrethroids are synthetic, modified versions of pyrethrins and target the closure of voltage-gated sodium channels in axonal membranes of insects²³. Lastly, many farmers, especially those with large acreages of crops that rely or benefit from animal pollination, use managed honeybees (Apis mellifera) for crop pollination. There is evidence that managed honeybees can negatively impact wild bees through competition for floral resources and disease transmission²⁴.

While there are multiple avenues by which agricultural intensification might harm wild bee populations, identifying mechanistic pathways has been difficult. This is partly owing to large data gaps^25,26, even in relatively well-sampled areas such as the United States. Similarly, lack of available data on historical application of different pesticides (both spatial and temporal) has hindered our ability to evaluate their full effects on communities. Indeed, most of the evidence that links pesticide use to bee declines comes from experimental or small-scale observational studies which may not reflect large-scale (for example, continent-wide) patterns. The majority of these studies were performed on a handful of species²⁷ (most notably the western honeybee Apis mellifera) and therefore may not be representative of the many other wild bee species. However, recent studies from the United Kingdom have demonstrated that population-level extinction rates at a country level are associated with the use of neonicotinoid seed treatment¹¹. Developing effective policies to protect wild bees requires understanding the causes of decline, and many of the purported causal factors are interlinked. Understanding the impacts of pesticides in relation to other potential drivers of decline is critical to sustainable management of ecosystems and food systems, and has important ecological and economic implications worldwide.

Here we address these challenges using one of the largest databases of bee records for the contiguous United States ever compiled; these were aggregated from museum specimens, surveys and community science observations²⁵. The United States has the largest number of described bee species on the planet with 3,594 species (including Hawaii and Alaska)²⁶, 17% of known species, and a large proportion of agriculture is intensive agriculture with recent land use change, making it an ideal place to understand the impacts of intensive agriculture on bee diversity. Our analysis included occurrence records for 1,081 bee species across the following families: 220 species from Andrenidae, 284 from Apidae (excluding Apis mellifera), 69 from Colletidae, 221 from Halictidae, 278 from Megachilidae and 9 from Melittidae (Fig. 1d). We combined multispecies occupancy models with tools of causal inference to estimate the effect of (1) pesticide use, (2) animal-pollinated agriculture (that is, agriculture that benefits from animal pollination) and (3) honeybee colonies on wild bee distributions across the contiguous United States. We did this via three types of multispecies occupancy model, each run independently on each bee family. We focused our analysis on crops that benefit from animal pollination because farms growing these may import honeybee colonies which may change their impacts on wild bees (Supplementary Fig. 4a vs b). Further, we expect that crops that benefit from animal pollination can provide additional resources known to attract wild bees. Finally, our analysis focuses specifically on the effects of neonicotinoids and pyrethroids because these classes of insecticides are thought to be linked to declines of wild pollinators¹, and increasing use of these compounds has increased risk to bees¹⁶.

Fig. 1: Counties across the contiguous United States differ in their pesticide application, percentage of agriculture, use of honeybee colonies and expected wild bee richness.

a, The distributions of pesticide use (2013–2015), modified from ref. ⁶⁵. Pesticide use was calculated by dividing the kg of compound used by the LD₅₀ (for honeybees), averaged over a 3-year period, summed across compounds and standardized to county area. b, Animal-pollinated agriculture (2008), modified from ref. ³³. The percent of the county that is animal pollinated removes the following crops: corn, wheat, rice, soybean, sorghum, barley and oat. Here we include all crops that would present a high exposure to pesticides for wild bees. c, Honeybee colonies (2012), modified from ref. ⁵¹. The number of honeybee colonies comes from the Census of Agriculture of farms that sell more than US$1,000. Values for a–c are log transformed and scaled. d, Expected wild bee species richness. Wild bee expected richness was compiled from ref. ²⁵ and divided by area of the county and log transformed (raw expected richness is presented in Supplementary Fig. 8).

We analyse individual species’ occupancy trends through space and time using a multispecies framework. In these models, individual species’ effects are drawn from distributions whose parameters are informed by data from all species. We analyse by family to reduce the computational burden. In addition, different families analysed separately provide semi-independent validation and, furthermore, bee species’ responses may be taxonomically grouped. We combined Melittidae and Colletidae (the two smallest families) because Melittidae had too few species to model on its own. In total, we ran 15 multispecies occupancy models (5 family groups × 3 effects). We modelled occupancy across all of the counties that fall within a species geographic range for 7 time periods, each lasting 3 years, from 1995 to 2015. We terminated in 2015 because the USGS Pesticide National Synthesis Project does not provide data on seed-coated neonicotinoids beyond 2015²⁸. We treated each of the 3 years within each time period as an opportunity for a potential visit by someone collecting bees. For each of these potential visits, we inferred species absences for a particular species at a site if another species within the same genus was observed at that site in that year^4,29,30. We obtained data on pesticide use application for every county and every year from the USGS Pesticide National Synthesis Project³¹ (Fig. 1a). Because the application of neonicotinoids and pyrethroids are correlated across space and time, we aggregated all active compounds (for both neonicotinoids and pyrethroids), controlling for their toxicity by weighting each by its median lethal dose (LD₅₀) as measured on honey bees (from the US Environmental Protection Agency (EPA) ECOTOX Database³²). In doing this, we are implicitly assuming that the relative LD₅₀ of different pesticides for honey bees is representative of the relative LD₅₀ of those same pesticides for native bees, but we recognize that LD₅₀s probably vary among wild bee species. We obtained honeybee colony data (the number of colonies per county across multiple years) from the National Agricultural Statistics Service (NASS) Census of Agriculture (Fig. 1c). While this census ignores small farm operations by only tracking farms that sell more than US$1,000 per year, it includes most large honeybee operations. Finally, we obtained agricultural distribution data from the Crop Data Layer³³ (Fig. 1b) and climatological data from CHELSA³⁴. Because we modelled occupancy across space and time, estimated effects can be driven by changes in predictors across both space and time.

To select predictors for the occupancy models, we relied on structural causal model analysis, which uses directed acyclic graphs (DAGs). DAGs are a visual representation of the presumed relationships between predictor variables and can help identify potential controlling variables in the context of multiple regression³⁵. In a DAG analysis, we first construct a diagram that includes all potential relationships between all potential predictor variables. This is done using previous knowledge about the system. For example, the amount of pesticide applied in a county probably depends on the amount of agriculture in that county. Next, we compared candidate DAGs using methods that leverage the statistical correlations among our predictor variables to arrive at a best DAG (Extended Data Fig. 1). Once we settle on a DAG that both represents potential relationships between variables and is consistent with the data, we identify the minimal set of variables needed to estimate the effect of our predictor(s) of interest. This analysis allows us to block confounding paths using the ‘backdoor criterion’. Namely, we identify which variables can confound the effect of the predictor(s) of interest. The main result is to avoid overcontrolling by including mediator or collider variables. While including many variables in an analysis can increase the predictive ability of a model, it may not allow for inference on effect sizes of the predictor of interest. Instead, by carefully considering potential confounders, colliders and mediators, we determine which minimal set of variables is needed to correctly estimate effect size(s) of the predictor(s) of interest. These steps are conducted before a regression analysis. Once a predictor set is chosen, we run the multispecies occupancy models. While a full overview of DAG analysis is beyond the scope of this paper, these steps have been summarized for ecological approaches in ref. ³⁶, have been identified as crucial for robust inferences in occupancy models in ref. ³⁷ and are explained in ref. ³⁵.

Results

We found that the mean effect of pesticide on species occupancy across all families was negative (Fig. 2a). This effect was strongest (95% Bayesian credible intervals do not overlap zero) for Andrenidae and Apidae, and Colletidae/Melittidae, and negative (but 95% Bayesian credible intervals overlap zero) for Halictidae and Megachilidae. For the different families, this translates into the following declines in mean occupancy probability (for a corresponding increase in pesticide use from zero to the maximum observed value in our dataset): 43.3% decline in Apidae, 28.9% in Andrenidae, 23% in Colletidae and Melittidae, 19% in Halictidae and 0.4% in Megachilidae (Fig. 3). These findings were qualitatively unchanged under multiple possible classifications of animal-pollinated agriculture. Namely, it did not matter whether we considered animal-pollinated agriculture those crops that (1) need pollination, (2) attract pollinators or (3) use managed pollinators (Extended Data Fig. 3a,d). Similarly, our findings were robust to our chosen duration of the occupancy interval (Supplementary Fig. 5) and whether we modelled the neonicotinoids or pyrethroids together or separately (Extended Data Fig. 4). When we aggregated these species-level effects within each genus, we found that effects varied by genus (Fig. 4a, and Supplementary Figs. 6a and 7a), with estimated effects of pesticide ranging from a 54% decline to a 62% increase (for a corresponding increase in pesticide use from zero to the maximum observed value in our dataset; Extended Data Fig. 5). We found that the mean effect of animal-pollinated agriculture was positive after controlling for climate (via temperature and precipitation) for Andrenidae, Colletidae and Melittidae, and Megachilidae (Fig. 2b), and had no effect for other groups. This effect was only consistent for Andrenidae; when we used a broader definition of animal-pollinated agriculture, credible intervals overlapped with zero for all other families (Extended Data Fig. 3b,e). This effect also varied by genus (Fig. 4b, and Supplementary Figs. 6b and 7b), with the estimated effect of animal-pollinated agriculture ranging from a 23.5% decline to a 401% increase for a corresponding increase in the percent animal-pollinated agriculture from zero to the maximum observed value in our dataset (Extended Data Fig. 6). Finally, to quantify potential effects of honey bees (Fig. 2c, Model 3), we controlled for the distribution of animal-pollinated agriculture. The mean effect of honeybee colonies on wild bee occupancy was not distinguishable from zero for all groups (Fig. 2c).

Fig. 2: Family-level mean effect sizes on occupancy.

a, Model 1. The mean effect size for pesticide use is strongly negative across five families. b, Model 2. The mean effect size of animal-pollinated agriculture is largely positive. c, Model 3. The mean effect size of honey bees is mixed across families. Asterisks mean that the 95% credible interval does not overlap with zero and error bars represent 95% Bayesian credible intervals. Sample sizes vary for each family but are the same for each model (Extended Data Table 1). Photos obtained from iNaturalist taken by Sam Droege under a CC-0 license. In order: Megachile fortis, Agapostemon angelicus, Colletes willistoni, Bombus griseocollis, Andrena polemonii. The family Apidae excludes Apis mellifera.

Fig. 3: Mean expected occupancy decreases as neonicotinoid and pyrethroid use increases across bee families.

The black lines represent the mean occupancy, the shaded grey lines are the 95% credible intervals. Occupancy probability was estimated using the posterior distribution of the mean intercept and mean slope effect estimated for pesticide use while keeping all other predictors at the mean value. Pesticide use was varied from the minimum to the maximum ever observed. Pesticide use is the combined sum of kg weighted by LD₅₀ and area, log transformed and scaled. Sample sizes vary for each family (Extended Data Table 1).

Fig. 4: Genus-level effect sizes for genera that contain at least ten species in our dataset.

a, Model 1. The mean effect of pesticide use is strongly negative across most genera. b, Model 2. The mean effect of animal-pollinated agriculture is largely positive. c, Model 3. The mean effect of honey bees is mixed across genera. Points denote mean effects and bars 95% Bayesian credible intervals. Points and bars are coloured grey if the credible intervals overlap zero, blue if they do not overlap zero and the mean effect is positive, and orange if they do not overlap zero and the mean effect is negative. We present coefficients for all genera of bees in Supplementary Fig. 6. Sample sizes vary for each family but are the same for each model (Extended Data Table 1).

Discussion

Across the contiguous United States, we found that higher pesticide use resulted in lower occupancy of wild bees. Here, pesticides were aggregated across neonicotinoid and pyrethroid compounds, hence we interpret the pesticide use effect as a measure of pesticide intensity by county. The negative effect of pesticide use was consistent across all five families of bees and for three families, the credible intervals did not contain zero. We also found negative effects of pesticides for important crop pollinator groups. For example, we found that bumble bees (genus Bombus) are negatively affected by pesticide use. Bumble bees are important pollinators of tomatoes, eggplants, peppers and melons, among other crops³⁸. Similarly, the genus Andrena was negatively impacted by pesticide use, and bees in this group are important pollinators of apples and a wide variety of native plants.

Laboratory and field studies have demonstrated that realistic levels of exposure to some types of pesticide have negative effects on bees. These negative effects include impairment of navigation, reduced foraging success, reduced longevity and reductions in reproductive health, among others¹. Some of these studies have been conducted across multiple countries^18,39 and provide a foundation for larger-scale future experiments. At country-wide scales, observational studies have demonstrated that these negative effects of pesticides do scale up to population-level trends¹¹. Our study adds additional such evidence by documenting decreasing species occupancy across hundreds of species, a finding that is consistent with experimental studies that have evaluated landscape-level pesticide risk for bumble bees³⁹. Further, we identify that certain groups, such as bumble bees, may be sensitive to pesticides, which could be driving stronger declines in those groups. This finding for bumble bees, in particular, is consistent with laboratory studies that have shown that they are particularly sensitive to some neonicotinoids such as acetamiprid and imidacloprid, and combinations of compounds such as clothianidin and propiconazole^27,40,41.

Pesticide use had a greater impact on native bees than the presence of honeybee colonies or the type of agriculture (animal pollinated), where the effect of animal-pollinated agriculture was largely positive and the effect of honeybee colonies was indistinguishable from zero. However, we note that the honeybee data that we used in this study (the Census of Agriculture for honey-producing colonies) were coarse both spatially and temporally, and do not account for the transient nature of honeybee hives that results from their transportation around the continent during the growing season. Further, our quantification of honeybee colonies at a county level may be too coarse a spatial scale to enable detection of their full effects. The positive (but not significant) association between honey bees and some wild bees may be a consequence of beekeepers targeting areas that are also good for wild bees such as those with ample floral resources. We also note that we did not explicitly quantify the effects of habitat conversion. Land cover and pesticide intensity in our model vary across space and through time, and we do not differentiate between those types of variation, hence our estimate should not be interpreted as necessarily indicative of temporal trends only. While conversion of natural and semi-natural areas into urban or industrial areas also affects wild bees⁴², we found that such habitat conversion was only weakly correlated with pesticide use and the fraction of animal-pollinated agriculture. Despite the lack of data, we want to highlight that this study was only possible due to the collection, curation and release of pesticide use data by the USGS, as well as the large-scale mobilization of bee occurrence data. Extension of this approach to other parts of the world requires both spatially explicit pesticide application records, honeybee colony location data, crop usage data and new bee inventories, particularly in undersurveyed areas.

While the use of occupancy models on presence-only data is becoming more commonplace^{2,3,5,6,29,43,44,45}, doing so requires making assumptions about how the data have been collected. In our case, we assumed that detection of some species provides information about the detection probability, and thus potential presence of other species (namely, those in the same genus). While this assumption may not always be valid (particularly when observers are only collecting a single species), recent studies have found that occupancy models are fairly robust to the inclusion of incomplete checklists⁴⁶ and yield unbiased results as long as a minimum of ~50% of the collection events target multiple species³⁰. Model caveats aside, there are no obvious mechanisms that would cause our use of presence-only data to generate statistical correlations between bee occupancy and pesticide application. Throughout the text, we have used the word ‘driver’ to link pesticides to changes in bee occupancy. We have used the strongest methods available to justify the use of such a term. However, the results of our models depend on the causal assumptions we made using the DAGs and, while our data were consistent with our DAGs, other hypotheses that include additional variables and/or links between variables might generate different outcomes.

Declining wild bee population health is troubling, particularly because wild bees are important pollinators of many crops and are essential to maintaining stability of our ecosystems in the United States. Our findings suggest that pesticide use may lower wild bee occupancy. These trends do not necessarily reflect population-level trends, particularly because changes in occupancy do not necessarily capture changes in abundance. Nevertheless, changes in occupancy only occur when a particular species either colonizes a previously vacant site or is extirpated at a site. Thus, declines in occupancy should cause concern, as the underlying patterns in abundance required to generate changes in occupancy are probably even more concerning. Analysing abundance patterns (not possible with available data) may reveal more troubling or wider-spread declines. Due to their high importance in agricultural settings and in pollinating native plants, identifying strategies that protect wild bees is necessary to ensure that these populations in the United States survive into the future. Techniques that enhance integrated pest and pollinator management (IPPM)⁴⁷, which aims to co-manage both pests and pollinators, can help minimize trade-offs of pest control and pollination services. For example, increasing habitat via hedgerows has been shown to be a cost-effective method to increase pest control and pollination services⁴⁸.

Methods

Bee records

We acquired North American occurrence records comprising six bee families (Andrenidae, Apidae, Colletidae, Halictidae, Megachilidae and Melittidae) from ref. ²⁵. In ref. ²⁵, records were downloaded from the Global Biodiversity Information Facility and Symbiota Collections of Arthropods Network, and were cleaned to ensure that (1) taxonomic names are correct and (2) erroneous records are removed. Further, ref. ²⁵ removed records for the western honeybee (Apis mellifera) and restricted records to only those within the contiguous United States. Ref. ²⁵ presented a total number of 1,923,814 occurrence records for 3,158 bee species from 1700 to 2021 (Supplementary Fig. 8). Of these, many are multiple observations of the same species on the same date at the same geographic location. Because the raw specimen data we used are presence-only records and because there is no measure of effort, the number of observations at a site is not a direct measure of abundance, hence treating it as such may not be appropriate. Further, because occupancy analyses use presence/absence data, rather than abundance, we reduced abundance measures to binary 0/1 (absence/presence). This resulted in a total of 634,597 unique species × date × location presence combinations (times and locations are explained below). Because our emphasis is on post-1990s pesticide application, we removed any observation from before 1994 or after 2016. Further, we removed any species that had less than 10 unique observations (unique dates and locations). We also removed species that were present in less than 3 years or 3 counties. We also removed all A. angelicus records and kept only A. texanus records because females cannot be morphologically differentiated between the two. Finally, our models only included genera for which we had data for at least three species, which resulted in a total of 1,081 bee species across the following families: Andrenidae 220 species, Apidae 284 species, Colletidae 69 species, Halictidae 221 species, Megachilidae 278 species and Melittidae 9 species.

Species ranges

For each species, we constructed plausible species ranges by identifying all counties that intersected a convex hull drawn around all observations of that species. This resulted in the plausible set of sites where that species may be found. While this may include counties where a species would probably not be found, subsetting to an approximate species range reduces the bias that might occur if one instead modelled every species at every site^4,30.

Agriculture and land use change

We obtained land cover data from the National Land Cover Database (NLCD), which provides nationwide data on land cover at a 30 m² resolution for every 2–3 years from 2001 to 2016. The NLDC provides data for 16 land use categories. Of relevance to us is the category of ‘Cultivated crops’ (hereafter agriculture), which contains areas used for the production of annual crops such as corn, soybeans and vegetables, perennial woody crops (such as orchards and vineyards), and land actively tilled⁴⁹. For each county, we calculated the proportion of the county covered in agriculture. While this database does not cover every year in our study, we interpolate the data within years. This is appropriate given that the amount of agricultural cover generally did not change appreciably within counties between 1995 and 2015.

With the NLCD, we also evaluated whether the loss of natural habitat was correlated with pesticide use or animal-pollinated agriculture. We defined as ‘Natural or semi-natural area’ the following land cover categories: Deciduous forest, Evergreen forest, Mixed forest, Woody wetlands, Emergent herbaceous wetlands, Shrub/Scrub, Grassland/Herbaceous, Sedge/Herbaceous and Pasture/Hay. We tested whether changes in natural or semi-natural areas were associated with the amount of pesticide use or agricultural cover. We found low correlations between the loss of natural areas and the amount of pesticide applied (between −0.157 and −0.206). Therefore, we do not consider further the effect of land use change as it relates to the effect of pesticide use or animal-pollinated agriculture.

Type of agriculture

We obtained crop cover data from the Crop Data Layer (CDL), which provides nationwide agricultural land cover data at a 30 m² resolution. Because the amount of agricultural cover did not change appreciably through time, we only used agricultural cover data from 2008 (the first year to cover all states) to quantify crop types for all years of the analyses. The CDL provides information on 137 types of cover³³. For each crop category, we also divided crops into ‘Non-animal-pollinated agriculture (NAPA)’ and ‘Animal-pollinated agriculture (APA)’. We considered three ways to categorize animal and non-animal-pollinated agriculture. First, we divided crops depending on whether they require pollination (presented in the main text); second, we divided crops on the basis of whether they use managed pollinators; and third, we divided crops depending on whether they attract pollinators (even if they do not require pollinators). For crops that do not require pollination, we included the following crops: corn, wheat, rice, soybean, sorghum, barley and oat. All of the other types of crop were considered in the APA category. Second, we used the USDA ‘Attractiveness of Agricultural Crops to Pollinating Bees for the Collection of Nectar and/or Pollen, 2017’. Here we separated crops on the basis of whether they used managed pollinators. This is because even though a crop may require animal pollination, it may not be supplemented by managed pollinators. The crops included in the category of not using managed pollinators were: barley, dry beans, chickpeas, corn, cotton, eggplants, garlic, grapes, hops, lentils, lettuce, oats, olives, oranges, peas, pistachios, potatoes, rice, rye, sorghum, soybeans, sugarbeets, sugarcane, sweet potatoes, tobacco, triticale, vetch, walnuts, winter wheat, durum wheat, spring wheat, turnips, celery, mustard, citrus, pecans, millet and flaxseed. All other crops were considered to supplement pollination via managed pollinators (Extended Data Fig. 7c). We also selected crops that do not necessarily require pollination but may attract pollinators and, therefore, expose wild bees to pesticides. In the ‘Non-attractive to pollinators agriculture’ category we included: barley, grapes, hops, lentils, oats, onions, pistachios, rice, rye, sugarcane, triticale, walnuts, winter wheat, durum wheat, spring wheat, pecans and millet. All other crops were considered to attract pollinators (Extended Data Fig. 7d). These additional categorizations of ‘animal-pollinated agriculture’ did not result in substantially different model estimates (Extended Data Fig. 3). For each county, we calculated the proportion of the county that is APA. Proportions were scaled and log transformed to ensure normality.

Pesticide use

We obtained pesticide use data from the USGS Pesticide National Synthesis Project, which provides national data on pesticide use for each county for every year from 1992 to 2021. The pesticide use data are provided as kg of active ingredient used in each county for 448 active ingredient types³¹. From these active ingredients, we considered neonicotinoids and pyrethroids because (1) they are known to be highly toxic to honey bees (via LD₅₀ measurements), (2) the usage of neonicotinoids increased during the study period, and the usage of pyrethroids remained constant, while the amount of kg used of other classes of insecticides decreased (Extended Data Fig. 2) and (3) because neonicotinoids and pyrethroids have been previously identified as problematic due to their increased toxicity when combined with fungicides¹. Specifically, we included in our analyses the following neonicotinoids: acetamiprid, clothianidin, dinotefuran, imidacloprid, thiamethoxam and thiacloprid; and the following pyrethroids: cyfluthrin, cypermethrin, permethrin, tefluthrin, tralomethrin, fenvalerate, deltamethrin, cyhalothrin-gamma, resmethrin and fluvalinate-tau. Because pesticides vary by orders of magnitude in toxicity, we weight each chemical by its contact LD₅₀, as measured on honey bees (following the procedure in ref. ⁵⁰), before we sum chemicals together. While LD₅₀ measurements collected on honey bees may differ from LD₅₀s for wild bees, such measurements for native bees generally do not exist. Further, we expect that LD₅₀ measured on honey bees should broadly reflect toxicity of different active ingredients for native bees. We acquired LD₅₀ estimates from the EPA ECOTOX Database³². For all results, we standardized the mean response to be in units of ng per bee, and we removed any LD₅₀ estimate that could not be standardized in this way. We only used LD₅₀ from studies that estimated acute toxicity (studies whose duration was less than 4 days) because they characterized most studies. We removed any study that did not report the duration. For each compound and using only contact exposure, we calculated the mean LD₅₀ across studies for each compound (Supplementary Table 1). We only used contact LD₅₀ because this was the most commonly tested across compounds. We set to zero each county × year × compound combination that did not have any pesticide application data. To calculate the total pesticide use, we summed the kg of pesticides used per year per county across LD₅₀ weighted compounds. We then divided the total amount of pesticide use by the area of the county to normalize for county size. Because we estimated occupancy at 3-year intervals (see below), we averaged 3 years of pesticide use for every county. To ensure the data were normally distributed, we log transformed the weighted sums and converted values into standard normal deviance to facilitate comparisons between different predictors. Because the application amounts of pyrethroids and neonicotinoids are highly correlated (correlation per year varies from between 0.30 when neonicotinoids were rarely applied in 1996 to 0.94 in 2008), we ran models where pesticides of both classes were combined or separated. Animal-pollinated agriculture was not highly correlated with pesticide use across all of the counties modelled (Supplementary Fig. 9)

Honeybee hives

The NASS produces surveys and reports in their Bee and Honey Program⁵¹. Since 2002, NASS has conducted the Census of Agriculture, which surveys any farm that reports revenue in excess of US$1,000. In this census, NASS collects information on the number of colonies, the number of farms, the amount of honey produced and the total sale of honey products. Thes data are available for the years 2002, 2007, 2012 and 2017. We divided the total number of colonies used by the area of the county to normalize for county size and then log transformed and scaled the data. Because these data lack seasonal information, they do not allow us to incorporate information about when the hives were actually present in farms.

Temperature and precipitation

We obtained climatic variables (both temperature and precipitation) from the CHELSA high-resolution climate data³⁴. Through the CHELSA time series, we were able to obtain monthly precipitation and maximum temperature at 30 arc-seconds resolution (~1 km). To calculate the yearly maximum temperature for every county, we extracted the maximum temperatures for the entire year and averaged the values within each county. Of these, we selected the maximum value per year. To calculate total precipitation, we extracted the precipitation for every month of the year, and averaged it within each county and across every year. We averaged both maximum temperature and mean precipitation across the number of time visits in a time interval to ensure consistency in the sensitivity analysis. We then scaled temperature and precipitation across all counties and years.

Data processing

Because pesticide use data are reported at county level, we calculated every other environmental predictor at county level and report each county as a ‘site’. We applied occupancy-detection models to historical data using the approach of ref. ³⁰. This approach requires specifying time intervals over which to estimate occupancy, and sub-intervals over which to estimate detection. We used 3 years as the time interval to estimate occupancy, and 1 year as the time interval to estimate detection. As a result, we did not model seasonal and within-year variation in occupancy. Recent studies have shown that trends inferred by occupancy models are robust to the choice of closure period (that is, the time period at which occupancy is estimated)⁴⁵ and spatial grain^5,45. We conducted a sensitivity analysis where we considered occupancy intervals from 2 to 5 years; results were qualitatively unchanged from those in the main text that use a 3-year interval (Supplementary Figs. 5, 10 and 11). We inferred that a visit to a site that might have detected a given species (for example, a potential of 0 or 1) occurred only when at least one other species in the same genus as that species was observed in the same year and county, and we only did this for sites that fall within that species’ range. Because of this, we excluded any species where there are fewer than three species in a genus. The application of occupancy models to presence-only data is a growing field^{2,3,5,6,29,43,44,45}. These studies have successfully leveraged the assumption that the presence of another species can inform non-detections for related species. Occupancy models can be robust to the breaking of this assumption, depending on the underlying data. For example, ref. ⁴⁶ showed that inferences from occupancy models are robust to the combination of complete and incomplete checklist data. Further, ref. ³⁰ found that using occupancy models that inferred non-detections on the basis of other species’ presences performed relatively well when applied to datasets comprising as much as 50% of sampling events where a collector was only seeking a single species.

However, it remains to be shown how sensitive these methods are to the taxonomic scope over which inferences about non-detections are extrapolated. We inferred non-detections across species at the genus, subfamily or family level. When inferring non-detections at the family level, we found that genera that are harder to identify tended to show declining trends (results not shown); however, when we inferred non-detections at the genus level, trends were consistent with those reported by experts. Our approach implicitly assumes that search effort, collection, identification and digitization of a species are comparable for different species within the same genus. While this may not be the case individually across all species, our analysis does allow for some genera to be under-identified or underdigitized compared with other genera via the non-detection imputation. We modelled each bee family separately, except for Colletidae and Melittidae, which we modelled together due to their containing a lower number of species.

Because not all families were observed in every county, we only modelled county × time interval combinations that had the presence of a species from a given family. Consequently, we modelled Andrenidae in 626 counties, Apidae in 1,763 counties, Halictidae in 1,046 counties, Megachilidae in 970, and Colletidae and Melittidae in 500 counties (Extended Data Fig. 8).

We matched the temporal scale of the occurrence data and all of the predictors. For pesticide use, temperature and precipitation, this process was fairly simple because all of these datasets were reported yearly. Therefore, we averaged them over the period over which we estimated occupancy (3 years in the main results and 2–5 years in the sensitivity analysis). For the number of honeybee colonies, we had to rely on the NASS census of agriculture data which are only available for the census years 2002, 2007, 2012 and 2017. Therefore, for the other periods, we interpolated values from the census year to any previous year for a given time period. For example, using 3-year time intervals, we had seven time periods (1995–1997, 1998–2000, 2001–2003, 2004–2006, 2007–2009, 2010–2012 and 2013–2015); the first 4 time intervals would use honeybee data from 2002, while interval 5 would use those from 2007, 6 would use those from 2012, and 7 would use those from 2017. We followed a similar process for the agricultural data from NLDC which is available from 2001 to 2016 at varying time intervals. The only predictor that we did not vary over time was animal-pollinated agriculture because this relied on the Crop Data Layer that was only available from 2008 onwards. Therefore, we treated this as a time-independent variable.

Causal analysis

To infer the causal effect of pesticide use, animal-pollinated agriculture and honeybee usage, we used DAGs to find the minimal set of adjustment variables, given each exposure set. First, we generated the DAG (Extended Data Fig. 1a). We used the following logic when building the DAG: we assumed that (1) conversion to agriculture has decreased the amount of floral and nesting resources available to wild bees⁵², (2) the expansion of agricultural landscapes has led to an increase in insecticide use⁵³, (3) animal-pollinated crops also require the supplementation of honey bees for pollination⁵⁴, (4) the presence of honey bees can affect the amount of pesticide used because one of the mitigation strategies proposed by the EPA to reduce the risk of pesticide to honey bees is to encourage the reduction of pesticide use via Managed Pollinator Protection Plans⁵⁵. Furthermore, the location of agriculture in the contiguous United States is largely informed by both climate and topographical and soil characteristics. Climate (both temperature and precipitation) also affects the availability of floral resources for bees. These relationships represent the initial DAGs we created (Supplementary Fig. 4a,b).

Methods that help build a DAG use the statistical correlations among predictors to support inclusion or exclusion of potential links between variables. For example, in our data, we observed that climate variables are correlated with honeybee colony abundance, suggesting that a link should connect these two predictors. This correlation might arise because climate variables directly impact honeybee colonies or, alternatively, it could be a consequence of a third variable that mediates this effect (for example, climate might impact the types of crop that are grown, which could subsequently impact whether or not honey bees are imported). Indeed, we found that including animal-pollinated agriculture as a mediator between climate and honeybee colony abundance is supported. Construction of our DAG proceeded by using this type of logic to determine which connections should be drawn between all of our variables (Supplementary Fig. 4). We then used our final DAG to identify the combinations of predictors needed to test our various hypotheses in our occupancy models. These models differed depending on which ‘exposure’ variable we were interested in (for example, pesticide vs animal-pollinated agriculture vs honeybee colonies). While we would ideally use DAG E in Supplementary Fig. 4, we were not able to obtain floral resource data that are independent from agriculture (for example, ref. ⁵ estimated floral resources on the basis of the Crop Data Layer), hence we used DAG D. Fortunately, for the exposure variables we were interested in here, the final model for each remained the same whether we used DAG C, D or E.

We used the R packages dagitty⁵⁶ and ggdag⁵⁷ to evaluate the minimum adjustment variables to estimate the causal effect of (1) pesticide use, (2) animal-pollinated agriculture and (3) honey bees, on wild bee occupancy (Extended Data Fig. 1). To estimate the effect of pesticide use, we included both animal-pollinated agriculture and honeybee surveys (Model 1: Extended Data Fig. 1b); for animal-pollinated agriculture, we only included climate in the model (Model 2: Extended Data Fig. 1c); and for honeybee usage, we only included animal-pollinated agriculture (Model 3: Extended Data Fig. 1d).

Occupancy models

Occupancy models aim to estimate the occupancy, denoted ψ, of a given species i on a site j in an era k by accounting for imperfect detection. They do this by incorporating multiple visits to the same sites. We developed three multispecies occupancy models for bee occurrence records in the contiguous United States that estimate the effect of pesticides, animal pollinated agriculture or honey bees on species occupancy.

For all occupancy models developed here, the probability that a species i is detected at site j in era k is drawn from a Bernoulli distribution with probability p_ijkz_ijk,

$${y}_{ijk} \sim {\rm{Bernoulli}}(\,{p}_{ijk}{z}_{ijk})$$

(1)

where p_ijk is the detection probability and z_ijk is the true (but unknown) occupancy state. If a site is occupied, then z_ijk = 1, and z_ijk = 0 otherwise. The true occupancy state is also drawn from a Bernoulli distribution with success probability ψ_ijk.

$${z}_{ijk} \sim {\rm{Bernoulli}}({\psi }_{ijk})$$

(2)

Both occupancy probability, ψ, and detection probability, p, can be formulated as functions of covariates, and we did this in three different ways for ψ.

Specifically, we used the following models:

Model 1

$$\begin{array}{ll}\psi [i,j,k]=&{\rm{expit}}\left({\beta }_{0}+\right.\\ &{\beta }_{{\rm{species}}}[i]+\\ &{\beta }_{{\rm{area}}}\times {\rm{area}}[\;j]+\\ &{\beta }_{{\rm{pesticide}}}[i]\times {\rm{pesticide}}[\;j,k]+\\ &{\beta }_{{\rm{honeybee}}}[i]\times {\rm{honeybee}}[\;j,k]+\\ &\left.{\beta }_{{\rm{apa}}}[i]\times {\rm{apa}}[\;j]\right)\end{array}$$

(3)

Model 2

$$\begin{array}{ll}\psi [i,j,k]=&{\rm{expit}}\left({\beta }_{0}+\right.\\ &{\beta }_{{\rm{species}}}[i]+\\ &{\beta }_{{\rm{area}}}\times {\rm{area}}[\;j]+\\ &{\beta }_{{\rm{apa}}}[i]\times {\rm{apa}}[\;j]+\\ &{\beta }_{{\rm{temp}}}[i]\times {\rm{temp}}[\;j,k]+\\ &{\beta }_{{\rm{temp2}}}[i]\times {\rm{temp}}{[\;j,k]}^{2}+\\ &\left.{\beta }_{{\rm{prec}}}[i]\times {\rm{prec}}[\;j,k]\right)\end{array}$$

(4)

Model 3

$$\begin{array}{ll}\psi [i,j,k]=&{\rm{expit}}\left({\beta }_{0}+\right.\\ &{\beta }_{{\rm{species}}}[i]+\\ &{\beta }_{{\rm{area}}}\times {\rm{area}}[\;j]+\\ &{\beta }_{{\rm{honeybee}}}[i]\times {\rm{honeybee}}[\;j,k]+\\ &\left.{\beta }_{{\rm{apa}}}[i]\times {\rm{apa}}[\;j]\right)\end{array}$$

(5)

where the expit function is defined as:

$${\rm{expit}}(x)=\frac{{e}^{\,x}}{1+{e}^{\,x}}$$

(6)

Here, β₀ denotes mean occupancy, β_species[i] denotes a species-specific random intercept, β_area denotes a fixed effect of area on occupancy, to account for the fact that some counties are larger than others. β_apa[i], β_temp[i], β_prec[i], β_honeybee[i] and β_pesticide[i] denote species-specific linear effects of animal-pollinated agricultural cover (apa), temperature, precipitation, honeybee colonies, and pesticide use. β_temp2[i] denotes a species-specific quadratic effect of temperature.

We assumed that species-specific slopes and intercepts are normally distributed about some mean. Specifically,

$$\begin{array}{ll}{\beta }_{{\rm{species}}}[i]& \sim {\mathcal{N}}(0,{\sigma }_{\beta {\rm{species}}})\\ {\beta }_{{\rm{apa}}}[i]& \sim {\mathcal{N}}({\mu }_{\beta {\rm{apa}}},{\sigma }_{\beta {\rm{apa}}}),\\ {\beta }_{{\rm{temp}}}[i]& \sim {\mathcal{N}}({\mu }_{\beta {\rm{temp}}},{\sigma }_{\beta {\rm{temp}}}),\\ {\beta }_{{\rm{temp2}}}[i]& \sim {\mathcal{N}}({\mu }_{\beta {\rm{temp2}}},{\sigma }_{\beta {\rm{temp2}}}),\\ {\beta }_{{\rm{prec}}}[i]& \sim {\mathcal{N}}({\mu }_{\beta {\rm{prec}}},{\sigma }_{\beta {\rm{prec}}}),\\ {\beta }_{{\rm{honeybee}}}[i]& \sim {\mathcal{N}}({\mu }_{\beta {\rm{honeybee}}},{\sigma }_{\beta {\rm{honeybee}}})\\ {\beta }_{{\rm{pesticide}}}[i]& \sim {\mathcal{N}}({\mu }_{\beta {\rm{pesticide}}},{\sigma }_{\beta {\rm{pesticide}}}),\end{array}$$

(7)

where μ_βapa, μ_βtemp, μ_βtemp2, μ_βprec, μ_βhoneybee and μ_βpesticide denote the mean effect of each corresponding predictor across species, and σ terms denote the variances about these means.

In all of the above models, we modelled detection probability as

$$\begin{array}{ll}{\rm{logit}}(p[i,j,k])=&{\alpha }_{0}+\\ &{\alpha }_{{\rm{era}}}\times k+\\ &{\alpha }_{{\rm{species}}}[i]+\\ &{\alpha }_{{\rm{site.era}}}[\;j,k]\end{array}$$

(8)

where α₀ denotes the mean detection probability, α_era denotes a slope of detection through time, α_species denotes a species-specific random intercept of detection, and α_site,era[j, k] denotes a site-specific random effect that is era specific. This latter term allows detection to vary relatively independently across sites and between eras. Specifically, we assumed

$$\begin{array}{ll}{\alpha }_{{\rm{species}}}[i]& \sim {\mathcal{N}}(0,{\sigma }_{\alpha {\rm{species}}})\\ {\alpha }_{{\rm{site.era}}}[\;j,k]& \sim {\mathcal{N}}({\mu }_{\alpha {\rm{site.era}}},{\sigma }_{\alpha {\rm{site.era}}}).\end{array}$$

(9)

We fitted 15 occupancy models, each of the types of model presented above (1, 2 and 3) for all five families of bees.

We ran models in JAGS⁵⁸ and assessed model convergence both by visually inspecting chains and checking the Gelman–Rubin statistic (we ensured that $\hat{R}$ was <1.1 for all parameters). We used flat, uninformative priors for all parameters and ran models for 100,000 iterations, with a burn-in of 1,000 iterations and thinning every 100 iterations across 3 chains, resulting in a total of 300 posterior samples. For all analyses, we used R (v.4.3.2)⁵⁹. For spatial manipulations, we used the package sf (v.1.0-15)⁶⁰; for data manipulation and visualization, we used tidyverse (v.2.0.0)⁶¹ and data.table (v.1.15.0)⁶²; and for running models, we used runjags (v.2.2.2-4)⁶³.

Reporting summary

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