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

Wild and managed bees provide essential pollination services for ~75% of crop species and 88% of flowering plants, yet many wild bee species are declining across Europe and North America. Multiple factors have been implicated, including climate change, land-use change and habitat loss, pathogens and disease, nutritional stress, and pesticide exposure—pressures often associated with agricultural intensification. Neonicotinoids and pyrethroids have expanded in use in the United States since the mid-1990s, with evidence of sublethal and chronic effects on bees that may scale to population impacts. However, large-scale, multispecies assessments across broad spatial and temporal extents remain limited, and most prior work focuses on a few model taxa (e.g., Apis mellifera). This study asks how pesticide use, the extent of animal-pollinated agriculture, and honey bee colony numbers relate to wild bee occupancy across the contiguous United States, aiming to disentangle these effects using causal inference and multispecies occupancy modeling applied to one of the largest compilations of U.S. bee records.

Prior research documents widespread pollinator declines and highlights the importance of wild pollinators for crop yield independent of managed honey bee abundance. Laboratory and field studies show neonicotinoids and pyrethroids impair navigation, foraging, longevity, and reproduction in bees, with sublethal, chronic effects that are difficult to detect in regulatory tests; the EU instituted restrictions and a 2018 ban on key neonicotinoids. Observational studies at national scales (e.g., England) have linked neonicotinoid use to increased extinction rates in wild bees, and landscape-level pesticide risk negatively affects bumble bees. Additional literature indicates managed honey bees can compete with wild bees and transmit diseases. Despite this, most large-scale inferences have been constrained by data gaps in bee occurrences and historical, spatially explicit pesticide applications, and by reliance on a few focal species. The present work builds on this foundation by integrating extensive occurrence records with county-level pesticide use, agricultural, climatic, and honey bee data, guided by structural causal models.

Study scope and data: The analysis focused on the contiguous United States from 1995 to 2015, modeling occupancy for 1,081 wild bee species across six families (Andrenidae 220, Apidae 284 excluding Apis mellifera, Colletidae 69, Halictidae 221, Megachilidae 278, Melittidae 9). Occurrence data (museum specimens, surveys, community science) were compiled from prior work (GBIF and SCAN), cleaned for taxonomy and errors, restricted to the contiguous U.S., and reduced to unique species × county × year presences. Records before 1994 or after 2016 were removed, species with <10 unique observations or present in <3 years or <3 counties were excluded, and genera required ≥3 species. This yielded 634,597 unique presence combinations and 178,589 unique observations within the 1995–2015 analysis frame. Species ranges were approximated as the set of counties intersecting the convex hull around occurrence points to define plausible sites. Predictors: Pesticide use (USGS Pesticide National Synthesis Project) was compiled annually at county level for neonicotinoids (acetamiprid, clothianidin, dinotefuran, imidacloprid, thiamethoxam, thiacloprid) and pyrethroids (cyfluthrin, cypermethrin, permethrin, tefluthrin, tralomethrin, fenvalerate, deltamethrin, cyhalothrin-gamma, resmethrin, fluvalinate-tau). Active ingredients were weighted by honey bee contact LD50 (EPA ECOTOX), standardized to ng/bee, summed across compounds, divided by county area, averaged over 3-year intervals, log-transformed, and scaled. Because neonicotinoid and pyrethroid use were correlated, models considered both combined and, in sensitivity analyses, separated. Honey bee colonies: County-level counts came from USDA NASS Census of Agriculture (years 2002, 2007, 2012, 2017), divided by county area, log-transformed and scaled; values were carried to adjacent non-census years to match 3-year occupancy intervals. Animal-pollinated agriculture (APA): From USDA Cropland Data Layer (CDL; nationwide coverage in 2008), APA was computed as the proportion of county area in crops requiring/benefiting from animal pollination. Three definitions were tested: crops requiring pollination; crops using managed pollinators; and crops attractive to pollinators. Proportions were log-transformed and scaled. Because CDL was only fully available in 2008, APA was treated as time-invariant in the main models. Land cover: NLCD data (2001–2016) were used to quantify agricultural cover and assess correlations between natural/semi-natural habitat loss and pesticide use/APA, which were weak (|r| ~ 0.16–0.21); land conversion was not included as a mediator in the primary causal sets. Climate: CHELSA monthly precipitation and maximum temperature at ~1 km resolution were aggregated to county-level yearly summaries (max annual temperature; total precipitation), averaged over each 3-year interval, and scaled. Causal analysis and model design: Directed acyclic graphs (DAGs) were constructed to represent hypothesized relationships among variables (climate, topography/soil, APA, honey bees, pesticide use, floral/nest resources, bee occupancy). Using dagitty/ggdag, minimal adjustment sets were determined for estimating each exposure's causal effect while avoiding overcontrol of mediators/colliders: Model 1 (pesticide exposure) adjusted for APA and honey bees; Model 2 (APA exposure) adjusted for climate (temperature—linear and quadratic—and precipitation); Model 3 (honey bee exposure) adjusted for APA. Occupancy modeling: Multispecies occupancy-detection models were fit separately for each family (Colletidae and Melittidae combined due to small sample size), across seven 3-year eras: 1995–1997, 1998–2000, 2001–2003, 2004–2006, 2007–2009, 2010–2012, 2013–2015. Detection sub-intervals were 1 year. Presence-only data were leveraged by inferring non-detections: within a species’ range, a year–county was treated as a potential visit if another species in the same genus was observed that year in that county; non-detections were inferred for congeners not recorded. Three models differed by occupancy covariates: Model 1 included pesticide, APA, honey bees, county area, and species random intercepts; Model 2 included APA, temperature (linear and quadratic), precipitation, area, and species effects; Model 3 included honey bees, APA, area, and species effects. Species-specific slopes for covariates were drawn from normal hyperdistributions to estimate family-level means and among-species variation. Detection probability included era trend, species random intercept, and site×era random effect. Inference and computation: Models were implemented in JAGS via runjags with uninformative priors, 3 chains, 100,000 iterations, burn-in 1,000, thinning 100, yielding 300 posterior samples; convergence assessed by chain inspection and Gelman–Rubin R < 1.1. Sensitivity analyses varied the closure interval (2–5 years), APA definitions, and separated pesticide classes; results were qualitatively robust.

Pesticides: Across all families, the mean effect of pesticide use (combined neonicotinoids and pyrethroids, LD50-weighted) on occupancy was negative. Family-level declines in mean occupancy probability, corresponding to an increase in pesticide use from the dataset minimum to maximum, were: Apidae −43.3%, Andrenidae −28.9%, Colletidae/Melittidae −23.0%, Halictidae −19.0%, Megachilidae −0.4%. For three families (Andrenidae, Apidae, Colletidae/Melittidae) 95% credible intervals for the mean effect did not overlap zero; for Halictidae and Megachilidae they included zero. Genus-level effects varied widely, from approximately −54% to +62% changes in occupancy over the observed pesticide gradient, indicating heterogeneous sensitivities. When separating classes, neonicotinoids showed strongly negative mean effects for Andrenidae and Apidae, while pyrethroids showed strongly negative mean effects for Apidae, Colletidae/Melittidae, and Halictidae. Animal-pollinated agriculture (APA): After controlling for climate, mean APA effects were generally positive for Andrenidae, Colletidae/Melittidae, and Megachilidae, and not different from zero for other groups; consistency across alternative APA definitions was limited, with only Andrenidae remaining consistently positive. Genus-level APA effects ranged from −23.5% to +401% across the observed APA gradient, suggesting that some genera benefit substantially from APA-associated resources, while others do not. Honey bees: After controlling for APA, mean effects of honey bee colonies on wild bee occupancy were not distinguishable from zero at the family level. Robustness: The negative pesticide effect held across alternative APA categorizations, occupancy interval lengths (2–5 years), and when modeling pesticide classes together or separately. Notable sensitive groups included Bombus and Andrena, both showing negative responses to pesticide intensity.

The analysis demonstrates that county-level pesticide intensity—captured by LD50-weighted neonicotinoid and pyrethroid use—is a widespread negative correlate and inferred driver of reduced wild bee occupancy across the contiguous United States. This signal is consistent across five major bee families and aligns with experimental literature documenting sublethal and chronic impacts on individuals that scale to population and community patterns. Important crop-pollinating genera such as Bombus (pollinators of tomatoes, peppers, melons) and Andrena (pollinators of apples and native flora) were negatively affected, indicating potential risks to both agricultural production and native plant reproduction. In contrast, the presence of honey bee colonies showed no detectable mean effect at county scale after adjustment, potentially due to data coarseness and the transient nature of hive movements; positive associations in some cases may reflect shared selection of resource-rich areas by beekeepers and wild bees. APA tended to have positive or neutral effects, consistent with added floral resources in animal-pollinated cropping systems, though responses varied by genus and definition of APA. The findings suggest that pesticide management may yield larger benefits for wild bee persistence than manipulating honey bee stocking densities, and that integrated pest and pollinator management can reduce trade-offs between pest control and pollination services. While occupancy integrates both spatial and temporal variation and cannot isolate purely temporal trends, the observed declines in occupancy imply underlying abundance reductions sufficient to drive local extirpations or reduced colonization, warranting conservation concern.

Using one of the most extensive compilations of U.S. bee records linked to spatially explicit pesticide and agricultural data, this study provides continent-scale evidence that increased use of neonicotinoids and pyrethroids is a major driver of reduced occupancy across hundreds of wild bee species. Effects were strongest for Apidae and Andrenidae, with heterogeneous responses among genera. Animal-pollinated agriculture often had neutral to positive associations with occupancy, while honey bee colonies showed no consistent effect at county scale. These results reinforce experimental findings and highlight pesticide reduction as a key lever for sustaining wild pollinators. Future work should (1) obtain finer-resolution, temporally explicit honey bee colony data and seasonal dynamics; (2) improve pesticide application records, including seed treatments post-2015 and co-application with fungicides; (3) develop independent, spatially resolved floral and nesting resource datasets; (4) incorporate abundance data where feasible to complement occupancy; and (5) extend this framework to other regions with improved inventories to inform global policy and IPPM strategies such as habitat enhancement (e.g., hedgerows) and reduced pesticide reliance.

Occupancy models infer presence while accounting for imperfect detection but do not measure abundance; declines in occupancy may underestimate or mischaracterize abundance changes. Presence-only historical data required assumptions about inferring non-detections from congeneric detections; while prior work suggests robustness when many events include multiple species, bias may persist when sampling targets single species or identification varies among genera. Honey bee colony data were coarse (county-level, census years only, lacking seasonal/transhumance dynamics), potentially obscuring localized competitive or disease effects. Pesticide toxicity weighting relied on honey bee contact LD50 values as proxies for native bees, omitting oral/chronic endpoints and species-specific sensitivities; LD50 data for wild bees are scarce. The pesticide metric aggregates across compounds and space-time, conflating spatial and temporal variation; seed-coated neonicotinoid data were unavailable beyond 2015. APA was treated as time-invariant due to CDL limitations, and habitat conversion effects were not modeled explicitly; although correlations with pesticide use and APA were weak, unmeasured habitat mediators could remain. County-level aggregation may mask fine-scale exposure and habitat heterogeneity. Causal inference depends on DAG assumptions; omitted variables or alternative structures could alter estimated effects.