Rising insecticide potency outweighs falling application rate to make US farmland increasingly hazardous to insects

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 insect decline. Nevertheless, answers to basic questions about the spatiotemporal patterns of insecticide use remain elusive, due in part to the inherent complexity of insecticide use, and exacerbated by the dispersed nature of the relevant data, divided between several government repositories. Here, we integrate these public datasets to generate county-level annual estimates of total ‘insect 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 drivers of the observed changes, we divide insect toxic load into the components of extent (area treated) and intensity (application rate x potency). We show that while contact-based insect toxic load remained relatively steady over the period of our analysis, oral-based insect toxic load increased roughly 9-fold, with reductions in application rate outweighed by disproportionate increases in potency (toxicity/kg) and increases in extent. This pattern varied markedly by region, with the greatest increases seen in Heartland and Northern Great Plains regions, likely driven by use of neonicotinoid seed treatments in corn and soybean. In this “potency paradox,” US farmland has become more hazardous to insects despite lower volumes of insecticides applied, raising serious concerns about insect conservation and highlighting the importance of integrative approaches to pesticide use monitoring. Significance statement Previous analyses disagree about whether US insecticide use is increasing or decreasing, a question of significant importance given the putative role of insecticides in recent insect declines. We integrated information from multiple national databases to estimate ‘insect toxic load’ (represented as honey bee lethal doses) of the agricultural insecticides applied in each US county from 1997 to 2012, and factors responsible for its change. Across the US, insect toxic load – calculated on the basis of oral toxicity – increased 9-fold. This increase was due to increases in the potency (toxicity/kg) of insecticides applied and in the area treated; the volume of insecticides applied declined. Toxic load increased most dramatically in regions where neonicotinoid seed treatments for field crops are commonly used.


Introduction 60
Insects are the most diverse and abundant class of animals on earth, with an estimated 5.5 million 61 species that dominate animal biomass in many ecosystems (1,2). Given their ubiquity, it is not 62 surprising that insect populations serve in key roles as both friend and foe to human societies. 63 This is particularly true in agriculture, where farmers seek to manage populations of insect pests 64 to produce essential food, fuel and fiber, a task in which insecticides play an important role. In 65 the United States, agriculture accounts for ~57% of insecticide weight applied and ~85% of area 66 treated, and so constitutes the single largest contributor to insecticide use (3,4). However, since 67 at least the 1960's it has been widely recognized that insecticide application can also negatively 68 affect non-target species, including populations of insect pollinators and natural enemies that 69 serve to support crop production (5,6). There has been a concomitant effort to reduce reliance on 70 insecticides and/or minimize their non-target effects, including state and federal regulation, 71 Integrated Pest Management (IPM) programs, and the development of alternative pest-72 management technologies and production/marketing systems such as organic farming (7). 73 Nonetheless, studies suggest recent and widespread declines in insect abundance, 74 diversity, and range (8)(9)(10), and insecticide use has been identified as a likely contributor along 75 with habitat loss, species introductions, and climate change (11). In the US, declines have been 76 documented in populations of several wild bee species, butterflies in Ohio and lowland 77 California, and the migratory monarch butterfly (12)(13)(14)(15), while beekeepers sustain losses of > 78 40% of their managed honey bee colonies annually (16). There is also evidence that populations 79 of some insect pests are declining (17,18). 80 To understand trends in pest management and potential impacts on target and non-target 81 insect species, we argue it is important to disentangle four distinct dimensions of agricultural 82 insecticide use, organized into two broad categories (Figure 1), extent and intensity. At the 83 landscape scale, the extent of insecticide use is a function of the area devoted to cropland as well 84 as the proportion of that cropland treated with insecticides. Historically in the US, a majority of 85 cropland has not been treated with insecticides (19), likely because insect pest populations have 86 not been problematic enough to justify treatment in low-value, large-acreage commodity crops 87 such as wheat, soybean, and corn. For those fields that are treated, the intensity of insecticide use 88 varies as a function of the application rate (kg/ha) of products applied and the potency 89 (toxicity/kg) of those products toward insects, typically measured via the LD50 (the dose required 90 to kill 50% of a population). These factors jointly determine insect toxic load, the total number of 91 insect lethal doses applied in a given area (as described in 20). 92 Past analyses of insecticide use on US farmland, focusing on only a limited subset of the 93 contributors to insect toxic load, have arrived at divergent narratives of overall trends in 94 insecticide use. For example, analyses based on farmer survey data from the US Department of 95 Agriculture (USDA) show long-term national declines in the application rate of insecticides in 96 corn and cotton crops, and corresponding reductions in the overall weight of insecticides applied 97 (19,21). A common interpretation of these trends is that introduction of transgenic insect-98 resistant (Bt) crops starting in the mid-1990s led to a reduction in insecticide application (19, 21). 99 Conversely, data from the US Agricultural Census indicate that the extent of US insecticide use 100 has increased since 1997, precisely in those regions dominated by Bt corn production (22-24). 101 The hypothesized cause of this increase in extent is the widespread adoption of neonicotinoid 102 seed treatments (23, 25). The first analysis did not consider insecticide potency, while the second 103 did not account for intensity. As these examples illustrate, the use of different metrics and data 104 sources leads researchers to conflicting interpretations of overall trends in insecticide use and 105 their associated drivers. 106 Importantly, describing aggregate trends in pesticide use is not very meaningful without a 107 consideration of the potency of diverse active ingredients (20, 26). In the US alone, there are 108 more than 1000 pesticide active ingredients registered for use (27), which vary by orders of 109 magnitude in potency to target and non-target organisms (28). Aggregate trends in application 110 rate or total weight applied can mask important shifts in the use of different insecticide products 111 over time. For insecticides, this is particularly problematic given evidence of a long-term trend 112 toward development of active ingredients with increasing potency toward insects (29). 113 Finally, insecticide use trends can vary substantially depending on the spatial scale 114 considered. Analyses at the national scale may obscure important local variation, particularly in 115 the US given the large size of the country, the concentration of different crop production systems 116 in particular areas, the divergent insecticide regimes associated with different crops, and regional 117 variation in pest pressure (24,30,31). Analyses at regional to local scales can help to target 118 insecticide mitigation efforts and facilitate research relating spatial patterns in insecticide use to 119 spatial patterns in important outcomes such as crop production, insecticide resistance, and 120 pollinator decline. Given limitations of existing data sources, most previous studies at the county 121 scale have focused solely on the extent of insecticide use, missing potentially important changes 122 in intensity. 123 Here, we integrate data from several national databases on the components of insecticide 124 extent and intensity ( Figure 1) to characterize spatiotemporal patterns in insect toxic load and its 125 drivers from 1997-2012 at the county scale. We use the honey bee (Apis mellifera) as a 126 representative insect species to calculate insect toxic load because it is the standard terrestrial 127 insect in regulatory pesticide tests, and so has the most comprehensive toxicity data available. 128 Analyses based on both contact toxicity (exposure to the outside of the body) and oral toxicity 129 (exposure via ingestion) captures uncertainty and variability in exposure pathways. Our goals 130 were to i) assess trends in insect toxic load on a contact and oral basis from 1997 to 2012, 131 nationally and for counties in the contiguous US, ii) identify the relative contribution of factors 132 (extent and intensity) responsible for these changes over the same time period, and iii) describe 133 regional variation in these trends among agricultural production regions. 134 135

Trends in insect toxic load 137
Between 1997 and 2012 insect toxic load across the contiguous US increased more than 9-fold 138 on an oral LD50 basis and was roughly constant on a contact LD50 basis, despite an overall 139 decline in the weight of insecticides applied ( Figure 2). Toxic load was spatially heterogenous 140 ( Figure 3); in 2012, the top 10% of counties accounted for 55% and 48% of contact and oral 141 toxic load, respectively. From 1997 to 2012, oral toxic load increased in 87% of counties and 142 decreased in 13% of counties (median change: +15-fold, IQR: +3.2-fold to +103-fold) while 143 contact toxic load increased in 51% of counties and decreased in 49% of counties (median 144 change: +5%, IQR: -56% to +142%). 145 In 2012, we estimate that agricultural insecticides were applied to ~5% of the land area of 146 the contiguous US (26% of cropland), at an average intensity of 5 billion honey bee contact 147 LD50s per treated ha and 16 billion honey bee oral LD50s per treated ha (Table 1). 148 149

Drivers of insect toxic load 150
The main driver of the increase in oral toxic load was insecticide intensity: a 16-fold increase in 151 oral potency far surpassed a 64% decline in application rate ( Changes in extent also occurred over the study period and influenced overall toxic load; 158 while total cropland declined slightly (-12%), the proportion of cropland treated with insecticides 159 increased by 78%, from 15% of cropland treated in 1997 to 26% in 2012. 160 161 162 163

Regional variation in insect toxic load and its drivers 164
The temporal dynamics of insect toxic load varied across the 9 major agricultural production 165 regions (Figures 3-5). Oral toxic load increased significantly in all regions except the Basin and 166 Range (Table S1) Conversely, changes in contact toxic load were more modest and variable, with most 175 regions experiencing no significant trends (Table S1). The exceptions included the Eastern 176 Uplands and Northern Crescent, where significant declines were observed (-38% and -10%, 177 respectively), driven by application rate, and to a lesser extent, declines in total cropland (Table  178 S1, Figure S3). Conversely, contact toxic load increased more than 3-fold in the Northern Great 179 Plains (Table S1, Figure S3), a pattern driven by a 2.7-fold increase in cropland treated and a  insecticides with LD50 values reported in high quality sources ( Figure S8). Finally, while we 194 used the more conservative USGS 'low' pesticide estimate in our analyses, repeating the 195 analyses with the 'high' estimate did not influence our conclusions regarding national and 196 regional trends, although numerical estimates for contact and oral toxic load differed ( Figure S9, 197 By creating a novel dataset synthesizing land use, insecticide use, and honey bee toxicity, we 201 characterized spatiotemporal patterns in the total inputs of insecticides to cropland in the 202 contiguous US. We found that insect toxic load from 1997 to 2012 was relatively constant when 203 calculated using contact LD50 values, but increased nine-fold when calculated using oral LD50 204 values, consistent with a recent analysis from the United Kingdom (20). Insect toxic load 205 increased despite a decrease in average application rate, which emphasizes the importance of 206 accounting for potency when assessing trends in insecticide use, as recently argued for herbicide 207 use (26, 32). The pattern of divergence between decreasing application rate and increasing extent 208 and potency, combined with significant variation across agricultural regions, helps to explain and 209 reconcile the conflicting narratives of insecticide dependency that characterize recent research on 210 US insecticide use trends (19, 21-24, 32). 211 Insect toxic load -calculated on an oral LD50 basis -increased in all regions, but the 212 increase was particularly acute in the Heartland and Northern Great Plains regions. We attribute 213 this pattern to the increasing use of neonicotinoid seed treatments in corn and soybean. In the 214 Heartland, > 90% of cropland (excluding pasture) in 1997 and 2012 was devoted to these two 215 crops ( Figure S2 Analyses derived from the Census of Agriculture suggest that seed treatments are at least partly 226 captured by that dataset (23, 24), although the survey language is ambiguous (36). Consequently, 227 our estimate of insecticide extent is likely an underestimate. Our results show that analyses of 228 insecticide use that do not account for seed treatments are likely missing a major component of 229 insect toxic load, highlighting a need for more consistent and comprehensive data collection. 230 Similarly, our results point to a need for more detailed investigation of the combined agronomic 231 and socioeconomic drivers of neonicotinoid seed treatments, which do not seem to be cleanly 232 related to pest pressure or field-level economics (25, 31). 233 Recommendations to farmers for reducing exposure of pollinators and natural enemies to 234 insecticides have historically emphasized avoiding contact exposure, for example by spraying at 235 night or not spraying during bloom (e.g. 6). Our finding that oral toxic load -driven by the use of 236 systemic neonicotinoids -has surpassed contact toxic load suggests that this guidance may need 237 to be amended. Systemic insecticides are taken up by plant tissues and can remain active for long 238 periods. Indeed, a major route of exposure for foraging bees to neonicotinoids and other 239 pesticides through nectar and pollen obtained from flowering weeds in agricultural areas, where 240 the weeds have taken up residues from the soil (37)(38)(39). Seed treatments likely exacerbate this 241 route of exposure, since neonicotinoids can remain present in the soil for years; indeed, recent 242 studies identified imidacloprid residues in plants and bee colonies, though imidacloprid use had 243 been discontinued in these study regions in previous years (37,40). It is sometimes possible to 244 adjust management practices to avoid non-target exposure to systemic pesticides, but this 245 requires longer-term planning than mitigation at the time of spraying. In apples, for example, 246 spraying neonicotinoids 5-10 days prior to bloom ensures levels are sufficiently low to protect 247 pollinators and natural enemies (41). 248 Importantly, our approach is not a formal risk assessment, nor does it fit neatly into 249 existing risk assessment frameworks (e.g. 42). The formal concept of risk--the likelihood of 250 adverse outcome--is defined as the product of hazard and exposure. Insect toxic load corresponds 251 to the hazard component of this definition, but we do not attempt to estimate or model pesticide 252 fate in the environment or exposure. Formal risk assessments are also, by necessity, narrow in 253 scope, typically evaluating single compounds with respect to specific receptor organisms over a 254 narrow range of use scenarios. In contrast, our approach integrates across compounds and use 255 scenarios to estimate cumulative, landscape-scale insecticide hazard to insects in general, 256 represented by a honey bee proxy. In this sense, we present an insect "hazard cup", analogous to 257 the "risk cup" paradigm used in human risk assessment (43) and echoing Berenbaum's (44)  258 extension of this paradigm to honey bee risk assessment. 259 The development of a county-level measure of insect toxic load, aggregated across 260 insecticides and national in scope, advances our ability to evaluate the impacts of insecticide use 261 on target and non-target insect populations. Previous large-scale studies have relied on the 262 percentage of land in agriculture (45), the percentage of land used to produce particular crops 263 (46), the percentage of land treated (23, 47), the weight of insecticides applied (48), or weight of 264 specific compounds (14,49). By developing a pipeline that integrates information from national 265 databases to summarize insecticide intensity, extent, and insect hazard (toxic load), our study 266 provides a common and repeatable framework that can be used across the US to investigate the 267 potential impacts of insecticide use on diverse insect populations and dependent wildlife. This study required synthesizing data on insecticide use, land use, and toxicity from multiple 279 national databases to generate a novel dataset documenting county-level spatiotemporal patterns 280 in insect toxic load and its drivers. All data processing and analysis was performed in the R 281 statistical language, version 3.6.1 (53). The code used to generate and analyze the dataset is 282 stored on GitHub and the workflow will be made publicly available upon publication.

Toxicity data 315
We considered toxicity on the basis of contact and oral exposure. Contact exposure occurs when 316 an insect encounters the insecticide on the outside of its body, for example if it is directly 317 sprayed. Oral exposure occurs when an insect ingests the insecticide, for example when feeding 318 on a recently treated plant. We conducted our analyses separately with each toxicity measure 319 because of the difficulty of predicting the dominant mode of exposure from insecticide use data 320 alone (particularly given that the data do not include method of application). 321 Honey bee acute toxicity values were compiled mainly from US-EPA's ECOTOX 322 database (58) and the PPDB. The ECOTOX database was queried in July 2017 and values were 323 recorded from the PPDB in June 2018. Both sources include toxicity data generated during 324 pesticide regulatory procedures as well as from other sources (e.g. scientific journal articles). We 325 gave preference to values generated in the regulatory processes of the US-EPA and the European 326 Food Safety Authority (EFSA) because they tend to be generated using standardized testing 327 protocols. 328 To compile honey bee LD50 values for the 138 insecticides present in the USGS dataset, 329 we first searched ECOTOX for all lab-generated LD50 estimates for Apis mellifera. We then 330 processed the data by recoding exposure types into 'contact', 'oral', or 'other' and standardizing 331 units where possible into ug/bee, consulting original sources as needed to verify extreme values 332 and to fill in missing information. The following criteria were used to select records: i) exposure 333 time 4 days or less (e.g. acute, after 59), ii) contact or oral exposure, iii) tests on adult bees, and 334 iv) non-zero LD50 estimate in units of µg/bee. We related this cleaned dataset to USGS pesticide 335 names on the basis of CAS numbers. Similarly, we generated a dataset comprising honey bee 336 LD50 values (acute contact and oral) from the PPDB, along with the source code identifying 337 high-quality values coming from the EU regulatory process (code 'A5'). 338 To generate a consensus list of contact and oral LD50 values for all insecticides reported 339 in the USGS dataset (Supplemental Data 1) we used the following procedure (adapted from 60): 340 (1) If point estimate(s) were available from US or EU regulatory bodies, we used those 341 values, taking a geometric mean when estimates were available from both sources; 342 otherwise, 343 (2) If point estimate(s) were available from other sources in ECOTOX or PPDB, we used 344 those values, taking a geometric mean when estimates were available from multiple 345 sources; otherwise, 346 (3) If unbounded estimate(s) were available from US or EU regulatory bodies (i.e. "greater 347 than" or "less than" some value), we used the minimum (for "less than") or the maximum 348 (for "greater than"); otherwise, 349 (4) If unbounded estimate(s) were available from other sources in ECOTOX or PPDB, we 350 used the minimum (for "less than") or the maximum (for "greater than"); otherwise, 351 (5) We used the median toxicity value for the insecticide mode-of-action group, with group 352 defined by the Insecticide Resistance Action Committee. 353 (6) In rare cases (n = 1/138 compounds for contact toxicity and 8/138 compounds for oral 354 toxicity) after this procedure we were still left without a toxicity estimate for a particular 355 compound. In those cases, we used the median for all insecticides. 356 357

Data synthesis and analysis 358
We synthesized the data sources described above to generate a novel dataset describing insect 359 toxic load and its contributors (Supplemental Data 2). First, we calculated the total contact and where weight is the total weight (kg) of insecticide active ingredient (ai) applied in a county in a 364 particular year, and toxicity is the contact or oral acute toxicity to honey bees (LD50 in µg/bee) 365 for each ai. There were a minority of counties for which insecticide use data were missing in 366 particular years. We used linear interpolation to fill in missing values for kg applied, contact 367 toxic load, and oral toxic load. If agricultural insecticide use was never reported in a particular 368 county (this occurred in a few, highly urban counties like New York City), we assumed it was 369 zero. 370 Insect toxic load in a county is a function of the area of land treated with insecticides 371 (extent) and the aggregate toxicity of the insecticides applied per unit area (intensity) (Figure 1). 372 To better elucidate drivers of insect toxic load over the study period, we calculated elements of 373 insecticide extent and intensity for those years in which land use data were available from the 374 agricultural census (1997,2002,2007,2012). Data on county-level land use were joined with 375 data on insecticide use and toxic load on the basis of county FIPS codes (Supplemental Data 3). 376 Once joined, we identified and calculated the contributors to insect toxic load as defined in 377 Figure 1. 378 To test for monotonic trends in insect toxic load at the national and regional scales we 379 used the non-parametric Mann-Kendall test (after 26). To further analyze change over time, and 380 to compare the relative magnitude of change across indicators with disparate units, for each of 381 our insecticide indicators we calculated fold-change from 1997 to 2012 as the response ratio: 382 Value2012/Value1997, so that a value of one represents no change, two represents a doubling, and 383 one half represents a decline of 50%. 384 To analyze regional differences in the temporal dynamics of toxic load, we performed 385 hierarchical time series clustering using the R package 'dtwclust' (61). Contact and oral toxic 386 load were calculated for years 1997-2012 for each USDA Farm Resource Region. Because 387 regions differed by as much as two orders of magnitude in absolute toxic load, we converted the 388 absolute toxic load within each region to fold-change relative to 1997, as described above, 389 enabling more informative clustering driven by patterns of relative toxic load. We then 390  Table 1. Estimates of 2012 insect toxic load and its contributors associated with insecticide use in the nine USDA Farm Resource 424 Regions of the contiguous US. 425

Sensitivity analysis
We performed several analyses to characterize the robustness of our results to uncertainty in the underlying data sources. One source of uncertainty in our estimates stemmed from interpolating missing values for insecticide use, crop area, and treated area in counties with missing values in our source data. However, the sensitivity analysis revealed that interpolation had a minor effect on estimates of insect toxic load ( Figures S5-S7). At the national scale, pesticide use was interpolated for fewer than 2% of counties in each year, and these counties contributed less than 0.05% to weight applied, contact toxic load, and oral toxic load. Interpolation was more frequent for treated cropland; nonetheless cropland treated was interpolated for < 6% of county-year combinations and these interpolated values contributed < 1% to total cropland treated. Across data items, interpolation tended to be rare in the regions driving our results (Heartland, Northern Great Plains), and more common in regions with generally little cropland (e.g. Basin and Range; Figures S5-S7).
Another source of uncertainty in our analyses concerned the honey bee LD50 values underlying the translation of insecticide use into insect toxic load. For each US county, we calculated the percentage of insect toxic load contributed by compounds with 'high', 'medium', and 'low' uncertainty in their LD50 values. Uncertainty was considered low if LD50 values were derived from US or EU regulatory procedures, medium if LD50 values were compound-specific but derived from other scientific sources, and high if LD50 values were estimated by a median value for the class or insecticides as a whole. In each year of the dataset active ingredients with highquality values (those generated as part of the US or EU regulatory processes) accounted for the vast majority of insecticide weight (79-86%), contact toxic load (91-98%), and oral toxic load (98-100%; Figure S7). Low-quality values (those missing and so estimated from class medians) contributed modestly to insecticide weight (1.5-15%) but very little to contact toxic load (< 1%) and oral toxic load (< 3%), likely reflecting that active ingredients without honey bee toxicity data are more likely to be found in classes with generally low toxicity.
To test the sensitivity of our findings to different methods of insecticide estimation, we repeated our main analyses using the less-conservative USGS 'high' estimates of insecticide use ( Figures  S3, S8, S9; Table S2). Qualitative patterns in insect toxic load and its drivers were extremely similar for the two estimates, supporting our conclusions regarding overall trends and regional variability in insecticide indicators. However, there were quantitative differences. For example, the 2012 estimate of insect toxic load at the national scale based on the 'high' USGS estimate was 39% higher on a contact-toxicity basis and 7% higher on an oral-toxicity basis than the equivalent value using the 'low' estimate (Table S2).   Figure S3. Change in the drivers of insect toxic load from 1997 to 2012 for all agricultural production regions, based on the USGS low (left) and high (right) pesticide use estimates. Foldchange is calculated as a response ratio: Value2012/Value1997, so that a value of one represents no change, two represents a doubling, one half represents a decline of 50% (values are presented on a log scale).

Figure S4. Contact toxic load by region and chemical class.
Toxic load time series were constructed for each of the nine USDA Farm Resource Regions. Hierarchical clustering (left) 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 Northern Great Plains from the other eight regions. The contribution of different chemical classes to the overall toxic load pattern in each region is depicted with stacked bar plots (right). Figure S5. Contribution of counties with interpolated insecticide values to total counties (n_interp), mass applied (kg_low), contact toxic load (ct_tox), and oral toxic load (or_tox) for each of nine agricultural regions and all regions combined. Figure S6. Contribution of counties with interpolated cropland area values to total counties (n_interp) and cropland (crop_ac) for each of nine agricultural regions and all regions combined. Figure S7. Contribution of counties with interpolated treated area values to total counties (n_interp) and treated cropland (trt_ins_ac) for each of nine agricultural regions and all regions combined. Figure S8. Contribution of LD50 values of low, medium, and high uncertainty to national estimates of insect toxic load. LD50 values are considered low-uncertainty if they are derived from US or EU regulatory procedures, medium-uncertainty if they are compound-specific but derived from the general scientific literature, and high-uncertainty if they are based on class or insecticide median values. Figure S9. National analysis repeated with the USGS 'high' insecticide use estimate. * California pesticide use data excludes products applied as seed treatments so insecticide toxic 6 load may be underestimated in these regions, which contain California counties. 7 8 9 Table S2. Estimates of 2012 insect toxic load and its contributors for agricultural production regions of the U.S. associated with 10 insecticide use, based on the less-conservative 'high' estimate from the USGS National Pesticide Synthesis Project. 11  * California pesticide use data excludes products applied as seed treatments so insecticide intensity and toxic load may be 14 underestimated in these regions, which contain California counties 15 16