County-level analysis reveals a rapidly shifting landscape of insecticide hazard to honey bees (Apis mellifera) on US farmland

Each year, millions of kilograms of insecticides are applied to crops in the US. While insecticide use supports food, fuel, and fiber production, it can also threaten non-target organisms, a concern underscored by mounting evidence of widespread decline of pollinator populations. Here, we integrate several public datasets to generate county-level annual estimates of total ‘bee toxic load’ (honey bee lethal doses) for insecticides applied in the US between 1997–2012, calculated separately for oral and contact toxicity. To explore the underlying components of the observed changes, we divide bee toxic load into extent (area treated) and intensity (application rate x potency). We show that while contact-based bee toxic load remained relatively steady, oral-based bee toxic load increased roughly 9-fold, with reductions in application rate outweighed by disproportionate increases in potency (toxicity/kg) and extent. This pattern varied markedly by region, with the greatest increase seen in Heartland (121-fold increase), likely driven by use of neonicotinoid seed treatments in corn and soybean. In this “potency paradox”, farmland in the central US has become more hazardous to bees despite lower volumes of insecticides applied, raising concerns about insect conservation and highlighting the importance of integrative approaches to pesticide use monitoring.

have not been problematic enough to justify treatment in low-value, large-acreage commodity crops such as wheat, soybean, and corn. For those fields that are treated, the intensity of insecticide use varies as a function of the application rate (kg active ingredient [AI]/ha) of products applied and the potency (toxicity/kg AI) of those products toward insects, typically measured via the LD 50 (the dose required to kill 50% of a population). These factors jointly determine insect toxic load, the total number of insect lethal doses applied in a given area 18,19 . Because the honey bee (Apis mellifera) is the standard terrestrial insect in regulatory pesticide tests and so has the most comprehensive toxicity data available, thus far, calculations of this or similar metrics have depended on honey bee data 18,19 . We have followed this convention and will henceforth refer to this value as "bee toxic load. " Past analyses of insecticide use on US farmland have arrived at divergent narratives of overall trends in insecticide use. For example, analyses based on farmer survey data from the US Department of Agriculture (USDA) show long-term national declines in the application rate of insecticides in corn and cotton crops, and corresponding reductions in the overall weight of insecticides applied 17,20 . A common interpretation of these trends is that introduction of transgenic insect-resistant (Bt) crops starting in the mid-1990s led to a reduction in insecticide application 17,20 . Conversely, data from the US Agricultural Census indicate that the extent of US insecticide use has increased since 1997, precisely in those regions dominated by Bt corn production [21][22][23] . The hypothesized cause of this increase in extent is the widespread adoption of neonicotinoid seed treatments 22,24 . The first analysis did not consider insecticide potency, while the second did not account for intensity. As these examples illustrate, the use of different metrics and data sources leads researchers to conflicting interpretations of overall trends in insecticide use and their associated causes.
Importantly, describing aggregate trends in pesticide use is not very meaningful without a consideration of the potency of diverse active ingredients 18,25 . In the US alone, there are more than 1000 pesticide active ingredients registered for use 26 , which vary by orders of magnitude in potency to target and non-target organisms 27 . Aggregate trends in application rate or total weight applied can mask important shifts in the use of different pesticide products over time. For insecticides, this is particularly problematic given evidence of a long-term trend toward development of active ingredients with increasing potency toward insects 28 .
Two recent analyses examined trends in bee toxic load at the national scale in Great Britian 18 and the US 19 . Goulson et al. 18 found that in Great Britain, from 1990 to 2015, the quantity (kg) of pesticide applied decreased, but the amount of treated agricultural land (ha) increased and the toxic load to bees increased 6-fold. DiBartolomeis et al. 19 demonstrated that in the US, from 1992 to 2014, there was an overall increase in bee toxic load by 4-fold (calculated using honey bee contact LD 50 ) or 48-fold (calculated using honey bee oral LD 50 ). DiBartolomeis et al. further found that neonicotinoids in corn and soy crops were the primary drivers of this increase.
Building on these national analyses, the research presented here characterizes for the first time the spatial dimension of changes in bee toxic load, by estimating this value and its contributing elements of insecticide extent and intensity at the county scale. In the US, there is considerable regional variation in insecticide use patterns due to the large size of the country, the concentration of different crop production systems in particular areas, the divergent insecticide regimes associated with different crops, and regional variation in pest pressure 23,29,30 . Analyses at regional to local scales are vital to target insecticide mitigation efforts and to facilitate research relating spatial patterns in insecticide use to spatial patterns in important outcomes such as crop production, insecticide resistance, and pollinator decline. Differentiating components of insecticide extent and intensity further elucidates syndromes of insecticide use in agricultural landscapes, with potential implications for both pest management and conservation of non-target species. Here, we integrate data from several national databases on the components of insecticide extent and intensity ( Fig. 1) to characterize spatiotemporal patterns in bee toxic load and its components from 1997-2012 at the county scale. This time period was chosen in part for pragmatic reasons relating to data availability, but perhaps more importantly it corresponds with the window in which several significant technological and policy changes occurred in US agricultural pest management: the passage of the Food Quality Protection Act in 1996, the introduction of transgenic Bt crops starting in 1996, and the widespread use of neonicotinoid seed treatments starting in the mid-2000s. Our goals were to (i) assess trends in bee toxic load on a contact and oral basis from 1997 to 2012, nationally and for counties in the contiguous US (ii) identify the relative contribution of extent and intensity responsible for these changes over the same time period, and (iii) describe regional variation in these trends among agricultural production regions.

Results
trends in bee toxic load. Between 1997 and 2012 bee toxic load across the contiguous US increased more than 9-fold on an oral LD 50 basis and was roughly constant on a contact LD 50 basis, despite an overall decline in the weight of insecticides applied ( Fig. 2; note that the weight of active ingredients, not formulations, are used in our analyses). Toxic load was spatially heterogenous (Fig. 3); in 2012, the top 10% of counties accounted for 55% and 48% of contact and oral toxic load, respectively. From 1997 to 2012, oral toxic load increased in 87% of counties and decreased in 13% of counties (median change: +15-fold, IQR: +3.2-fold to +103-fold) while contact toxic load increased in 51% of counties and decreased in 49% of counties (median change: +5%, IQR: −56% to +142%).
In 2012, we estimate that agricultural insecticides were applied to ~5% of the land area of the contiguous US (26% of cropland), at an average intensity of 5 billion honey bee contact LD 50 s per treated ha and 16 billion honey bee oral LD 50 s per treated ha ( Table 1). components of bee toxic load. The main contributor to the increase in oral toxic load was insecticide intensity: a 16-fold increase in oral potency far surpassed a 64% decline in application rate (Fig. 4). The increase in insecticide potency was associated with the shifting composition of insecticide classes, away from organophosphates and toward pyrethroids and neonicotinoids (Figs. 2, S1); by 2012 neonicotinoids alone contributed 98% of oral toxic load. While contact potency also increased from 1997 to 2012 (+76%), this change was roughly counterbalanced by a corresponding decline in application rate.
Changes in extent also occurred over the study period and influenced overall toxic load; while total cropland declined slightly (−12%), the proportion of cropland treated with insecticides increased by 78%, from 15% of cropland treated in 1997 to 26% in 2012.
Regional variation in bee toxic load and its components. The temporal dynamics of bee toxic load varied across the 9 major agricultural production regions (Figs. [3][4][5], which vary in crop composition (Fig. S2). Oral toxic load increased significantly in all regions except the Basin and Range (Table S1), with reductions in application rate overwhelmed by increases in potency and cropland treated (Fig. S3). The most dramatic increases occurred in the Heartland (121-fold increase) and the Northern Great Plains (53-fold increase) (Figs. [4][5]. This pattern was underscored by time series cluster analysis, where the majority of total variation was captured in the node dividing the Heartland + Northern Great Plains cluster from the other seven regions (Fig. 5b). The cause of this pattern, evident when oral toxic load is parsed by chemical class, was the substantial application of neonicotinoid insecticides in the Heartland and Northern Great Plains, beginning in 2006 and increasing exponentially through 2012 (Fig. 5c). This likely reflects the dominance of corn and soybean in the Heartland and the shift to greater area under these crops in the Northern Great Plains (Fig. S2). www.nature.com/scientificreports www.nature.com/scientificreports/ Conversely, changes in contact toxic load were more modest and variable, with most regions experiencing no significant trends (Table S1). The exceptions included the Eastern Uplands and Northern Crescent, where significant declines were observed (−38% and −10%, respectively), driven by application rate, and to a lesser extent, declines in total cropland (Table S1, Fig. S3). However, contact toxic load increased more than 3-fold in the Northern Great Plains (Table S1, Fig. S3), a pattern associated with a 2.7-fold increase in cropland treated and a 3.8-fold increase in contact potency (Fig. 4). In the cluster analysis, the majority of variation was captured in the node dividing the Northern Great Plains region from all other regions (Fig. S4). Notably, organophosphate toxic load increased in the Northern Great Plains from 2006 to 2012, in contrast to all other regions where the contribution of this chemical class was stable or declining (Fig. S4). . Toxic load is normalized by the total area (ha) of each county, to account for variable county size. Note that California data excludes seed-applied insecticides. Color scales are based on Jenks natural breaks, which maximize similarity within classes and minimize similarity between classes, adjusted for (b,d) to center on zero. Black outlines show USDA Farm Resource Regions. Maps were created using the R statistical language version 3.6.1 (https://www.R-project.org) using the data generated in this study.

Contact intensity (bil bee LD 50 s/treated-ha)
Oral intensity (bil bee LD 50  www.nature.com/scientificreports www.nature.com/scientificreports/ Sensitivity analysis. We performed several analyses to characterize the robustness of our results to uncertainty in the underlying data sources and interpolation of missing values (Supplemental Material). We found that interpolation of insecticide use, crop area, and treated area in counties with missing values in our source data had only a minor effect on the overall results ( Figures S5-S7). Furthermore, while the LD 50 values for some of the insecticides used in our analysis were derived from sources with low and medium certainty, the vast majority of bee toxic load was derived from insecticides with LD 50 values reported in high quality sources (Fig. S8). Finally, while we used the more conservative USGS 'low' pesticide estimate in our analyses, repeating the analyses with the 'high' estimate did not change the direction or significance of national or regional trends. That said, numerical estimates differed (Fig. S9, Table S2). For example, at the national scale in 2012 the high estimate was 39% and 7% greater than the low estimate for contact and oral toxic load, respectively.

Discussion
By creating a novel dataset synthesizing land use, insecticide use, and honey bee toxicity, we characterized spatiotemporal patterns in the total inputs of insecticides to cropland in the contiguous US at the county scale. We found that bee toxic load from 1997 to 2012 was relatively constant when calculated using contact LD 50 values, but increased nine-fold when calculated using oral LD 50 values, with considerable variation among regions. Bee toxic load increased despite a decrease in average application rate, which emphasizes the importance of accounting for potency when assessing trends in insecticide use, as recently argued for herbicide use 25,31 . The pattern of divergence between decreasing application rate and increasing extent and potency, combined with significant variation across agricultural regions, helps to explain and reconcile the conflicting narratives of insecticide dependency that characterize recent research on US insecticide use trends 17,[20][21][22][23]31 . While this analysis does not allow us to directly evaluate exposure of pollinators to insecticides, the county-level dataset provided here can be used to test relationships between insecticide use and pollinator populations and to prioritize areas for in-depth risk assessments. Furthermore, these results can inform the collection of data on insecticide use and recommendations to reduce exposure of pollinators to insecticides. www.nature.com/scientificreports www.nature.com/scientificreports/ DiBartolomeis et al. 19 found a greater increase in overall oral toxic load, of 48-fold compared to our calculation of 9-fold. This difference is likely due to the use of a slightly different time period, and the inclusion of an environmental half-life in the toxic load calculations, which increases the relative contribution of long-lasting insecticides. While we appreciate the value of incorporating degradation, we chose not to do so because pesticide half-lives depend on the matrix considered (e.g. soil vs. plant surface vs. water), and the pesticide source data do not include information on application method that could inform the selection of a value for half-life. For example, the Pesticide Properties Database reports a soil half-life under field conditions of 174 days for the neonicotinoid imidacloprid, while the residual half-life of the same active ingredient in the plant matrix is 4.9 days 32 . Nonetheless, it is interesting to note that the inclusion of degradation by DiBartolomeis et al. only magnified the trend toward increasing toxic load, suggesting that our analysis may be conservative.
We found that the increase in oral toxic load was particularly acute in the Heartland and Northern Great Plains regions, which showed a 121-fold and 53-fold increase, respectively. We attribute this pattern to the increasing use of neonicotinoid seed treatments in corn and soybean (see Figs. 2c, 5c, S9). In the Heartland, >90% of cropland (excluding pasture) in 1997 and 2012 was devoted to these two crops (Fig. S2), while the Northern Great Plains experienced significant land use change over the study period in which corn and soybean displaced other crops and conservation land 33 . Neonicotinoids accounted for the overwhelming majority of oral toxic load by 2012, and previous research 24 showed that virtually all neonicotinoid use in corn and soybean is via seed application. These results are consistent with DiBartolomeis et al., which also found that neonicotinoid use in corn and soy was the primary driver of the observed increase in oral toxic load at the national scale. The map was created using the R statistical language version 3.6.1 (https://www.R-project.org). Hierarchical clustering (b) grouped regions with similar patterns of toxic load using a Euclidean distance matrix and Ward's linkage method. The y-axis distance between two nodes and their nearest common node is inversely proportional to similarity. Hence, the majority of variation among regions is captured in the first split that separates the Heartland and Northern Great Plains from the other seven regions. The contribution of different chemical classes to the overall toxic load pattern in each region is depicted with stacked bar plots (c).

Scientific RepoRtS |
(2020) 10:797 | https://doi.org/10.1038/s41598-019-57225-w www.nature.com/scientificreports www.nature.com/scientificreports/ The important role of neonicotinoid seed treatments in this analysis underscores another source of conflicting interpretations of US insecticide use trends -inconsistent or incomplete data collection. USDA started collecting data on seed treatments only in 2015 34 , and the data are not yet publicly available. The USGS Pesticide National Synthesis Project included seed-applied pesticides through 2014 (except in California), then stopped when the data provider excluded them 35 . Analyses derived from the Census of Agriculture suggest that seed treatments are at least partly captured by that dataset 22,23 , although the survey language is ambiguous 36 . Consequently, our estimate of insecticide extent is likely an underestimate. Our results show that analyses of insecticide use that do not account for seed treatments are likely missing a major component of bee toxic load, highlighting a need for more consistent and comprehensive data collection. For similar reasons, our findings on insecticide use trends for counties in California should be interpreted cautiously given the lack of seed treatment data in that state. Similarly, our results point to a need for more detailed investigation of the combined agronomic and socioeconomic drivers of neonicotinoid seed treatments, which do not seem to be cleanly related to pest pressure or field-level economics 24,30,37,38 .
Recommendations to farmers for reducing exposure of pollinators and natural enemies to insecticides have historically emphasized avoiding contact exposure, for example by spraying at night or not spraying during bloom 6 . Our finding that oral toxic load -driven by the use of systemic neonicotinoids -has surpassed contact toxic load suggests that this guidance may need to be amended. Systemic insecticides are taken up by plant tissues where they can remain active, sometimes for long periods. Indeed, a major route of exposure for foraging bees to neonicotinoids and other pesticides is through nectar and pollen obtained from flowering weeds in agricultural areas, where the weeds have taken up residues from the soil [39][40][41] . Seed treatments likely exacerbate this route of exposure, since neonicotinoids can remain present in the soil for years 42 ; indeed, recent studies identified imidacloprid residues in plants and bee colonies, though imidacloprid use had been discontinued in these study regions in previous years 39,43 . It is sometimes possible to adjust management practices to avoid non-target exposure to systemic pesticides, but this requires longer-term planning than mitigation at the time of spraying. In apples, for example, spraying neonicotinoids 5-10 days prior to bloom ensures levels are sufficiently low in nectar and pollen to protect pollinators and natural enemies 44 . Furthermore, in a recent study of authorized uses of neonicotinoids in agricultural systems in France, in 96% of cases, an effective alternative was available, and in 78% of cases, an effective non-chemical alternative was available (such as physical or biocontrol methods), suggesting it is possible to alter practices to continue managing pests but reduce off-target effects 45 .
Importantly, our approach is not a formal risk assessment, nor does it fit neatly into existing risk assessment frameworks like that of the US-EPA 46 . The formal concept of risk-the likelihood of adverse outcome-is defined as the product of hazard and exposure. Bee toxic load corresponds to the hazard component of this definition, but we do not attempt to estimate or model pesticide fate in the environment or exposure. Formal risk assessments are also, by necessity, narrow in scope, typically evaluating single compounds with respect to specific receptor organisms over a narrow range of use scenarios. Our analysis incorporated >100 active ingredients applied in varied ways to dozens of crops, and the source data does not include information on application timing and method that would be necessary to derive a mechanistic understanding of fate and exposure. Rather, our approach integrates across compounds and use scenarios to estimate cumulative, landscape-scale insecticide hazard to bees. In this sense, we present an bee "hazard cup", analogous to the "risk cup" paradigm used in human risk assessment 47 and echoing Berenbaum's 48 extension of this paradigm to honey bee risk assessment. While there is clearly variation across insect species in toxicity, honey bees do not appear to be consistently more or less sensitive to insecticides than other insect species 49 , and so "bee toxic load" may serve as a starting point for a general assessment of the toxicity of insecticides to insect populations in a given landscape. Future extensions of this framework could seek to integrate additional layers of complexity, for example synergy between jointly applied insecticide active ingredients or between insecticides and fungicides.
The development of a county-level measure of bee toxic load, aggregated across insecticides and national in scope, advances our ability to evaluate the impacts of insecticide use on target and non-target insect populations. Previous landscape-scale studies seeking to relate insecticide use patterns to invertebrate outcomes have relied on the percentage of land in agriculture 50 , the percentage of land used to produce particular crops 51 , the percentage of land treated 22,52 , the weight of insecticides applied 53 , or weight of specific compounds 14,54 . By integrating information from national databases to summarize insecticide intensity, extent, and bee hazard (toxic load), our study provides a novel dataset that can be used in the US to investigate the potential impacts of insecticide use on diverse insect populations and dependent wildlife. Indeed, a recent study in the Netherlands created a "potential pesticide risk map" based on pesticides allowed to be applied on crops, and found a negative association between increasing pesticide toxicity in the landscape and both honey bee colony survival and wild bumble bee distributions 55 ; the approach we demonstrate allows for a more empirical method for generating such maps, accounting for actual insecticide use practices, which vary regionally. While our estimates use data from 1997-2012, it is expected that the overall use patterns have remained consistent in recent years (if anything, use of neonicotinoid seed treatments in field crops has increased 56 ,) and many datasets documenting insect declines encompass the time period in our study 12,13,15 . Furthermore, these national maps of bee toxic load, in combination with maps examining insect abundance and ecosystem services 57 , can be used to identify regions that should be prioritized for conservation (areas where insect abundance, diversity, or ecosystem services are high and insecticide hazard is low) or mitigation (areas where abundance of imperiled species, such as monarch butterflies, is high 58 and insecticide hazard is also high). While there are clearly limitations to this analysis, it complements more traditional 'bottom-up' approaches to ecotoxicology by enabling 'top-down' research at the landscape scales at which insect populations and ecosystem services are structured 59 . (2020) 10:797 | https://doi.org/10.1038/s41598-019-57225-w www.nature.com/scientificreports www.nature.com/scientificreports/

Materials and Methods
This study required synthesizing data on insecticide use, land use, and toxicity from multiple national databases to generate a novel dataset documenting county-level spatiotemporal patterns in bee toxic load and its components. All data processing and analysis were performed in the R statistical language, version 3.6.1 60 . The code used to generate and analyze the dataset is documented on the following website: https://land-4-bees.github.io/ bee_tox_county. insecticide use data. Data on the weight of insecticides applied by county and year were downloaded from the US Geological Survey (USGS) Pesticide National Synthesis Project 61 . This dataset reports pesticide use (kg active ingredient applied in each county) for >400 of the most commonly applied herbicides, insecticides, and fungicides on agricultural crops. The estimates are derived from proprietary farmer surveys, except for California, where estimates are drawn from the state's Department of Pesticide Regulation 61 . The dataset includes seed-applied pesticides through 2014, except in California, where they have never been included. USGS provides two pesticide estimates; we used the more conservative 'low' estimate (data presented in the text and figures) but repeated our analyses with the 'high' estimate (data presented in the Supplementary Materials) to ensure that changes over time and comparisons among regions were robust to differences in pesticide estimation. We focused on the 138 insecticides in the dataset, classified as such based on the Pesticide Properties Database (PPDB) 32 and the Insecticide Resistance Action Committee (IRAC) 62 . Following the custom of the US Environmental Protection Agency (US-EPA) we excluded petroleum oils and distillates 4 .
Land use data. The US Census of Agriculture is conducted every five years. We used this dataset as a source of county-level data for the years 1997, 2002, 2007, and 2012 on (i) cropland area treated with at least one insecticide, (ii) total cropland area, and (iii) total land area in each county. Data were downloaded from the developer page of the Quick Stats database maintained by the National Agricultural Statistics Service 63 . We focused on area in cropland rather than pastureland because USGS data indicates that only ~1% of insecticides are applied to pastureland. Missing values occurred in the dataset when data were withheld to avoid disclosing information about individual farm operations; this situation occurs mostly in cases where agricultural acreage in a county is very low. In those cases, we imputed missing values with the average of values in other years for that county if possible, or with zeros if values for the county were missing in all years 23 . Finally, counties were grouped into regions using a system developed by US Department of Agriculture ('Farm Resource Regions') to capture geographic variation in farming systems related to cropping patterns, similarities in soil and climate, and economic characteristics 29 . toxicity data. We considered toxicity on the basis of contact and oral exposure. Contact exposure occurs when an insect encounters the insecticide on the outside of its body, for example if it is directly sprayed. Oral exposure occurs when an insect ingests the insecticide, for example when feeding on a recently treated plant. We conducted our analyses separately with each toxicity measure because of the difficulty of predicting the dominant mode of exposure from insecticide use data alone (particularly given that the data do not include method of application).
Honey bee acute toxicity values were compiled mainly from US-EPA's ECOTOX database 64 and the PPDB. The ECOTOX database was queried in July 2017 and values were recorded from the PPDB in June 2018. Both sources include toxicity data generated during pesticide regulatory procedures as well as from other sources (e.g. scientific journal articles). We gave preference to values generated in the regulatory processes of the US-EPA and the European Food Safety Authority (EFSA) because they tend to be generated using standardized testing protocols.
To compile honey bee LD 50 values for the 138 insecticides present in the USGS dataset, we first searched ECOTOX for all lab-generated LD 50 estimates for Apis mellifera. We then processed the data by recoding exposure types into 'contact' , 'oral' , or 'other' and standardizing units where possible into µg/bee, consulting original sources as needed to verify extreme values and to fill in missing information. The following criteria were used to select records: (i) exposure time 4 days or less 65 , (ii) contact or oral exposure, (iii) tests on adult bees, and (iv) non-zero LD 50 estimate in units of μg/bee. We related this cleaned dataset to USGS pesticide names on the basis of CAS numbers. Similarly, we generated a dataset comprising honey bee LD 50 values (acute contact and oral) from the PPDB, along with the source code identifying high-quality values coming from the EU regulatory process (code ' A5').
To generate a consensus list of contact and oral LD 50 values for all insecticides reported in the USGS dataset (Supplemental Data 1) we used the following procedure, adapted from previous research 66 : (1) If point estimate(s) were available from US or EU regulatory bodies, we used those values, taking a geometric mean when estimates were available from both sources; otherwise, (2) If point estimate(s) were available from other sources in ECOTOX or PPDB, we used those values, taking a geometric mean when estimates were available from multiple sources; otherwise, (3) If unbounded estimate(s) were available from US or EU regulatory bodies (i.e. "greater than" or "less than" some value), we used the minimum (for "less than") or the maximum (for "greater than"); otherwise, (4) If unbounded estimate(s) were available from other sources in ECOTOX or PPDB, we used the minimum (for "less than") or the maximum (for "greater than"); otherwise, (5) We used the median toxicity value for the insecticide mode-of-action group, with group defined by the Insecticide Resistance Action Committee. (6) In rare cases (n = 1/138 compounds for contact toxicity and 8/138 compounds for oral toxicity) after this procedure we were still left without a toxicity estimate for a particular compound. In those cases, we used the median for all insecticides.
Scientific RepoRtS | (2020) 10:797 | https://doi.org/10.1038/s41598-019-57225-w www.nature.com/scientificreports www.nature.com/scientificreports/ Data synthesis and analysis. We synthesized the data sources described above to generate a novel dataset describing bee toxic load and its contributors (Supplemental Data 2). First, we calculated the total contact and oral bee toxic load for each county-year combination from 1997 to 2012, using the following equation adapted from previous research 25  where weight is the total weight (kg) of insecticide active ingredient (ai) applied in a county in a particular year, and toxicity is the contact or oral acute toxicity to honey bees (LD 50 in μg/bee) for each ai. There were a minority of counties for which insecticide use data were missing in particular years. We used linear interpolation to fill in missing values for kg applied, contact toxic load, and oral toxic load. If agricultural insecticide use was never reported in a particular county (this occurred in a few, highly urban counties like New York City), we assumed it was zero. Bee toxic load in a county is a function of the area of land treated with insecticides (extent) and the aggregate toxicity of the insecticides applied per unit area (intensity) (Fig. 1). To better elucidate bee toxic load over the study period, we calculated elements of insecticide extent and intensity for those years in which land use data were available from the agricultural census (1997,2002,2007,2012). Data on county-level land use were joined with data on insecticide use and toxic load on the basis of county FIPS codes (Supplemental Data 3). Once joined, we identified and calculated the contributors to bee toxic load as defined in Fig. 1.
To test for monotonic trends in bee toxic load at the national and regional scales we used the non-parametric Mann-Kendall test. To further analyze change over time, and to compare the relative magnitude of change across indicators with disparate units, for each of our insecticide indicators we calculated fold-change from 1997 to 2012 as the response ratio: Value 2012 /Value 1997 , so that a value of one represents no change, two represents a doubling, and one half represents a decline of 50%.
To analyze regional differences in the temporal dynamics of toxic load, we performed hierarchical time series clustering using the R package 'dtwclust' 67 . Contact and oral toxic load were calculated for years 1997-2012 each USDA Farm Resource Region. Because regions differed by as much as two orders of magnitude in absolute toxic load, we converted the absolute toxic load within each region to fold-change relative to 1997, as described above, enabling more informative clustering driven by patterns of relative toxic load. We then performed hierarchical time series clustering using a Euclidean distance matrix and Ward's linkage method, with oral and contact toxic load analyzed separately.

Data availability
The R code and workflow associated with this paper is documented on the following website: https://land-4bees.github.io/bee_tox_county. The datasets generated through our analyses and their associated metadata are included in Supplemental Information.