Reduction of Crop Diversity Does Not Drive Insecticide Use

Reduction of crop diversity, on farm and landscape levels, has been suggested as a factor that leads to lower biodiversity in agricultural systems, and thus makes them more prone to pest damage. To determine whether the relationship between increase of corn production and reduction of crop diversity and the proportion of cropland treated with insecticide found previously in the Midwestern states is universally applicable, we applied spatial panel model analysis to USDA Agricultural Census data for 1997, 2002, 2007, and 2012 using county-level data for the entire continental USA. The 7 Midwestern states and the remaining 41 states were analyzed separately. Simpson’s diversity index was used as the metric for crop diversity. We also examined the effect of temperature, the proportion of all cropland that is corn and average crop market value as additional predictor variables. The results show that expansion of corn production, together with the market value of crop products, and accumulated amount of heat could be key factors driving the increase of acres treated with insecticide. This phenomenon is observed in the entire continental USA, indicating that it was not a reduction of crop diversity in general, but specifically the predominance of corn that is likely driving increased insecticide use in the Midwestern states in recent years.


Introduction
Large-scale application of insecticides can detrimentally impact human health and ecosystem functions (Chagnon et al., 2015;Whiting, et al., 2014).Despite pesticide regulations and improvements in integrated pest management (IPM) (EPA, 2014;Goulson, 2013), the proportion of county cropland area treated with insecticides in the continental US has increased over the past four USDA Agricultural Census periods (2002( -2012, Figure , Figure 1).Insecticide use intensity can be affected by a variety of biological and economic factors: climate, soil, population dynamics of pest and their natural enemies, crop yield, market value, and the farmer's perception of pest risk (Waterfield & Zilberman, 2012).
The expansion of monoculture, a key feature of modern intensified agriculture, for agricultural production is recognized as an essential factor contributing to increased insect pest pressure due to lower populations of natural enemies of pests and higher crop patch connectivity (Chaplin-Kramer et al., 2011;Letourneau et al., 2010;Lin et al., 2011;Matson et al., 1997;Tscharntke et al., 2007).Two recent studies that applied different statistical methods on the same variable sets have drawn controversial conclusions regarding the effect of landscape simplification on insecticide use intensity, using the proportion of cropland treated with insecticides as a proxy for "use intensity" (Larsen, 2013;Meehan et al., 2011).While Larsen (2013) used ordinary least squares and Meehan et al. (2011) used spatial error models, they both found that in 7 Midwest states, i.e.Illinois, Indiana, Iowa, Michigan, Minnesota, Ohio, and Wisconsin, 1) the proportion of corn acreage and 2) the proportion of fruits and vegetables in total harvested cropland of a county both significantly correlated to increased insecticide use in recent years.To understand the relationship between landscape simplification in cropping systems and percentage of cropland treated with insecticides, a relevant metric is needed for quantifying the state of simplification.Species diversity indices are commonly used to examine variance of species abundance and distribution, and can represent the diversity of crop landscapes in both homogeneous and heterogeneous areas (Gardiner et al., 2009).
Current methods for investigating insecticide usage on crops typically use either pounds of insecticide applied per acre of cropland (Vialou et al., 2008), or percentage of cropland treated with insecticide (Fausti et al., 2012;Larsen, 2013;Meehan et al., 2011;Tscharntke et al., 2007).National Agricultural Statistics Service (NASS) conducts an annual agricultural survey each June based on random sampling of growers and a Census of Agriculture (COA) every five years based on a comprehensive survey of all farm and ranch operators.The annual survey data provides records of amount of insecticides applied for each crop at the state level, and the COA data contains records of total crop acres treated with insecticides at the county level.The former method that only collects chemical use data ignores important properties such as active ingredient content and application rate, making it problematic to assess changes in insecticide use over time.Since reliable data for estimating chemical characteristics is not available at the county-level for the entire continental US, we developed a quantitative model using proportion of county cropland treated with insecticides as the metric of insecticide use.
Panel data refers to a dataset where a group of individual entities are surveyed repeatedly over several periods of time (Frees, 2004).Depending on the nature of the study and data availability, the entity of interest can be persons, zip codes, cities, regions, counties, states, etc.With the increasing availability of census data, panel data analysis is being used by agricultural economists and other social scientists with increasing frequency (Baylis et al., 2011).The classical panel data model is able to control for the heterogeneity across the entity of interest and/or time period (Baltagi, 2011).When the panel data of interest are geographical units, the existence of spatial autocorrelation can induce bias.In this study, we applied a spatial panel model that controls for both heterogeneity and spatial autocorrelation (Anselin, 1998(Anselin, , 2003;;Baltagi, 2008;Baylis et al., 2011).The objective of the study is to investigate whether crop diversity leads to the increased insecticide use throughout the USA by using the panel data of insecticide use, crop market value, crop composition, total acreages of cropland, and temperature.

Materials and Methods
Panel data for model analysis were mainly obtained from COA data, which is conducted every 5 years, and provides detailed farm and ranch data at the county level.Data used in this study includes 1997, 2002, 2007, and 2012 census years, and includes the following county-level variables: land area treated with insecticides, market value of all crops, total land area, total agricultural land area, total harvested area, and harvested acres of specific crops.All are in acres, as given the the COA, except for crop value which were converted to US dollars/hectare to be consistent with previous studies (Larsen, 2013;Meehan et al., 2011).
Crops that were produced for different purposes, e.g.processing tomato and fresh market tomato, were aggregated to single crop type.Crop landscape diversity is based on the crop specific acres harvested data.Numerous diversity indices have been created to measure the richness and abundance of species in an ecosystem.Shannon's diversity index and Simpson's diversity index, which differ in their theoretical foundations, are the most common (Morris et al., 2014;Wagner, et al., 2000).Both indices increase as the diversity increases, and Simpson's diversity Index was chosen in this study because it is more sensitive to abundant species (Begon, et al., 2005): (1) Where, D is the diversity index with no units, N is the total number of cropland acres, and n is the acres of each crop type.The data were available only for major crops, and we assumed that rare crops will not substantially affect the total acreages treated with insecticide.
The market value variable for the spatial panel data model was calculated using the data extracted from the COA database.The data cleaning procedures applied included removing records with extremely high values, and counties with variable values which fluctuated more than 10-fold during the census period.
In terms of crop diversity measures, NASS Cropland Data Layers (CDL) maps derived from satellite images and ground truths contain very detailed spatial information of crop distribution patterns (Boryan et al., 2011; USDA/NASS, 2015b).While the basic spatial unit of COA data is the county-level, the CDL maps represent crop distribution at a scale down to either 56-or 30-meter resolution.CDL maps have equivalent temporal resolution to COA data, but the available years are more limited.However, since CDL maps for most areas are only available for 2007 and 2012, we decided to use CDL as a reference to verify the diversity measure derived from COA data.We downloaded CDL data for the 7 Midwestern states to verify our evaluation of crop diveristy based on COA data.Landscape metrics, which is a set of spatial pattern indices designed to measure the spatial configuration of ecological habitats and other categorical maps, have been used to measure the status of landscape diversity of natural or cultivated landscapes (Bharath et al., 2012;Li & Wu, 2004).CDL maps that contain information of crop land distribution were used to calculate the landscape metrics.Percentage of Like Adjacencies (PLADJ) and Core Area Percent of Landscape (CPLAND) of the most widely cultivated crop in each county were selected as indicators of landscape simplification and were calculated using the R package SDMTools (McGarigal et al., 2002;Van Der Wal et al., 2014).
The number of heat units accumulated during a specific period is an essential factor pest population dynamics (Latham & Mills, 2011).We used accumulated growing degree-days (GDD) above 9.475 °C, the average base developmental threshold temperature among 287 pest arthropods listed in the Insect Development Database (IDD) (NAPPFAST, 2015;Nietschke et al., 2007).The result is as a measure of the heat available to drive the development and reproduction of arthropod pests.The annual accumulation is the sum of GDD for each day in a calendar year.The basic equation of GDD for each calendar day is: (2) Where, T MAX is the daily maximum air temperature, T MIN the daily minimum air temperature, and T BASE is the temperature below which the biological processes do not progress (McMaster & Wilhelm, 1997).While the IDD includes data for 500 arthropod species, those that are generally considered beneficial, such as parasitoids, were excluded from the calculation.The average base temperature we obtained is slightly lower than the 9.6 °C estimation of Nietschke et al. (2007).The difference might be generated from a difference in calculation methods, or our exclusion of non-pest arthropods.
County-level daily maximum and minimum air temperature data were obtained from the North America Land Data Assimilation System (NLDAS) Daily Air Temperatures and Heat Index website (CDC, 2012(CDC, ) for 1997(CDC, , 2002(CDC, , 2007(CDC, , 2010(CDC, , and 2011.Since the 2012 data is not available yet, we used the average of the same date in 2010 and 2011 as the proxy of 2012 data. The advantages of using the panel data model rather than a single year model are increasing the size of the dataset, which is more informative, and an ability to control of individual heterogeneity.One common issue in agricultural research is that it is impossible to collect data for all potential explanatory variables, e.g.farmer behavior, policy effects, etc.The effect of the unobserved variables can be included in a linear model as the intercept of the regression line.There are two standard approaches for estimating the intercept: random effects ( 1) and fixed effects (Baltagi, 2011).In this analysis, fixed effects were chosen over random effects based on the test suggested by Baltagi et al. (2003).
As shown in Figure 2, counties with a high proportion of insecticide treated acres tend to be clustered spatially, which indicates there might be spatial dependence.Spatial dependence amongst the data violates the independent error term assumption of ordinary least square regression.Two common approaches to address the spatial dependency issue are to include a spatial lag parameter in the dependent variable and spatial autoregressive disturbance (Elhorst, 2014;LeSage & Pace, 2009).A static (not accounting for temporal autocorrelation) spatial panel model takes into account the spatial lag of the dependent variable, the spatial autocorrelation of the error term, and the specific effects (whether fixed or random) (Baltagi, 2008;Kapoor et al., 2007;LeSage & Pace, 2009;Millo & Piras, 2012).The general form of the model based on Baltagi et al. (2003) and Kapoor et al. (2007) includes a spatial lag of the dependent variable and spatial autoregressive disturbances as follows (Millo & Piras, 2012): Where, y is a vector of the dependent variables of NT × 1 dimension, X is a matrix of explanatory variables of NT × k dimension, I T is an identity matrix of T dimension, W N is the N × N spatial weights matrix, and λ is the corresponding spatial parameter.The disturbance vector is: Where ι T is a T × 1 vector of ones, I N is an N × N identity matrix, μ is a vector of spatially independently time-invariant individual specific effects, and ε is a vector of spatially autoregressive disturbance with the following component: Where, ρ is the spatial autoregressive parameter, W N is the spatial weighted matrix, and Several spatial weight matrices can be used in the spatial panel model.The spatial weights matrix based on inverse distance was applied in this study because it accounts for both global (the entire modeling area) and local (neighboring counties) effects (Kostov, 2012).The spatial panel model was estimated using splm package in R.

Spatial Panel Model
Variables of interest included 1) proportion of county cropland area treated with insecticides, 2) SDI or proportion of corn in county cropland, 3) crop market value per hectare, and 4) accumulated growing degree-days (GDD), using 9.475 ºC as the average lower developmental threshold for 287 pest arthropods listed in the IDD.These variables were tested using the classical panel model without spatial accommodation.The residuals were then tested for spatial autocorrelation with Moran's I test and the result indicated significant spatial dependence that needed to be remedied (Blanchet et al., 2008).Counties where total cropland was less than 1% of the county area, and counties that did not have data for all variables in all years were excluded.County-level regions analyzed numbered 567 and 1,396 counties in the seven Midwestern states and the remaining 41 continental US states, respectively.
We applied the fixed effect spatial panel model that accounts for both spatial lag and spatial error effects using the splm package in R (Millo & Piras, 2012).The estimated parameters of explanatory variables and impact/effect of each variable of the county-level spatial panel model for the 7 Midwestern states and the 41 remaining states are summarized in Tables 1 and 2, respectively.
While discussed in greater detail in the Methods section, it is important to consider these results in light of the assumptions of the spatially and time-specific fixed effects models.The spatially fixed effect model (i.e., fixed effect on spatial entities, which are counties in this study) assumed there are some attributes of the counties that don't change over time but arbitrarily correlate with other independent variables-whether observable or unobservable, such as natural, cultural, and economical factors of a given county.
In this study, the county-specific fixed effects can include any unobserved or unmeasurable factors.The time-specific fixed model (i.e., fixed effect on time periods, which are years in this study) produces time-specific fixed effects that do not vary among the counties.Year-specific fixed effects can include climate, policy, or market events that affect cropping and thus insecticide use.
Results of the county fixed effect model show that SDI is significantly negatively related to increasing insecticide use in the seven Midwestern states (Table 1A).The λ coefficient at the lagged spatial dependent variable, when positive, indicates spatial dependences among the counties.The ρ coefficient of the spatial autocorrelation of the error term can be seen as a short-term unobservable disturbance.County-or year-specific fixed effects are similar to a dummy variable that absorbs the factors not included among the model variables.
The direct impact is used to test whether an explanatory variable has a significant effect on the dependent variable within a particular county, while the indirect impact is used to test whether the explanatory variable of neighboring counties has a significant effect on the dependent variable in the county (Elhorst, 2014;Millo & Piras, 2012).In the 41 states model, all λ and ρ coefficients are small and positive, which indicate spatial dependency does exist (Table 2).The coefficient of SDI in both the county-fixed-effect and the year-fixed-effect models are positive, indicating that county-specific, time-invariant, and unobserved variables are likely to play more important roles in insecticide use than cropland diversity.
The proportion of corn acres, on the other hand, is associated with higher proportion of insecticide treated acres.
The seven Midwestern states analyzed in this model are all of the major corn producing states in the USA (Figure 4), and characterized by low crop diversity.This result represents patterns similar to those found in Meehan et al. (2011).Crop diversity which is dominated by corn is likely to induce more insecticide use in Midwestern cropland.When we changed the model input to the data of the 41 other continental USA states, however, the coefficients of both SDI and corn proportion are positive.This indicates that both variables are positively related to an increase in insecticide use in these states.Not surprisingly, in counties where corn was not a dominant crop, corn proportion was not the major contributor of neither crop diversity nor simplification (Figure 4).
Based on data from insecticide phenology studies (e.g., McMaster & Wilhelm, 1997;Latham & Mils, 2011), temperature is expected to be positively related to pest pressure, and therefore insecticide use.However, when the year-specific fixed effect is controlled for, the growing degree-day variable is either not significant or has a negative coefficient in the Midwestern states model data (Table 1).The negative coefficient of the model variable in this case may indicate that spatial variation of growing degree-days in these Midwestern counties might not be great enough to lead to statistically significant differences in pest population development.The climatic homogeneity of this 7-state region, compared with the entire USA for example, is further driven by the fact that most of the cropland, including corn fields with high insecticide treated rate, are found in the southern and middle portion of this area (Figure 2).As discussed below, it is also likely that other year-specific factor(s), such as Bt corn resistance development, may have overshadowed the role of temperature in driving the changes of insecticide use intensity in this region.Unlike in the Midwest states model, the coefficients of accumulated growing degree-days variables are positive in the 41 states model, and the effect is higher in the year-fixed effect model.
The market value variables are significant in all models.In the county-specific fixed model, we assumed that each county has some non-random characteristics that do not change over time.The positive coefficient of the market value variable indicates that in a year when market value of crop is higher, the proportion of cropland treated with insecticide is also higher.The coefficients are larger in year-specific fixed effect models than in the county-specific fixed model.This indicates that when year-specific factor(s) are excluded, the variation of prices among different counties have a stronger effect on insecticide use.For the 41 states model, the effect of value variables in the county-fixed-effect model is higher than the coefficient in year-fixed-model, which emphasizes both that high value crops tend to receive more insecticide treatment and that each county tends to have its own insecticide use characteristics.
In the 41 states model with corn proportion variables, all estimated coefficients indicate a positive relationship between corn proportion and insecticide use (Table 2).The coefficients of corn proportion variables in these models are much smaller than in the Midwestern states models, but the impact is still relatively high when compared to other variables.Figure 4 shows that most corn production counties in the US has more than 40% of cropland dominated by corn fields, which indicates that corn production counties tend to have more proportion of cropland treated with insecticide.

Discussion
The purpose of this study is to understand the relationship between landscape simplification and insecticide application rates.Similar to Meehan et al. (2011) and ( 2013), we also considered that in the Midwestern states, the proportion of corn acreage in cropland could be a relevant metric of cropland simplification.Additionally, we adopted SDI as a metric that negatively correlated crop diversity with the corn proportion in both the Midwestern model and the 41 states model.The initial model variable statistics reveal a high correlation between corn proportion and cropland simplification in Midwestern states, and low correlation in the other 41 states.Figure4 shows that when the proportion of corn acres in a county was small, there was no apparent relationship between corn proportion and diversity; however, when the proportion of corn acres exceeded a level of about 40% crop diversity was negatively related to corn proportion.
The panel data models used in this study show that reduced crop diversity is not necessarily related to higher proportion of acres treated with insecticide among cropping systems with a wide range of geographic and climatic conditions.Insecticide use is dependent on the conditions of other potential drivers, however the increases in corn planted area over the past decade (Landis, Gardiner, van der Werf, & Swinton, 2008) might have led to increasing insecticide treated acres.
There are two likely reasons why corn proportion is positively related to increasing insecticide treatment area: Bacillus thuringiensis (Bt) resistance and biofuel-driven increases in corn acreage and value.Bt corn is genetically modified to produce a protein that is toxic to various insects (Benbrook, 2012).Bt crops have been planted in the US since 1996, and represented about 65% of corn and 75% of cotton planted in 2011 (Benbrook, 2012).The evolution of resistance in the corn rootworm to the widely planted Bt corn since the mid-2000s is likely to necessitate increased chemical insecticide use in affected areas (Gassmann et al., 2011).In addition, there is some evidence that Bt corn fields may support reduced populations of non-target insects, including natural enemies of pests (Marvier et al., 2007).During the four NASS census periods included in this study, the proportion of corn acreage in the seven Midwestern states increased from 38% to 43%, and the proportion of corn acres in the remaining 41 states fluctuated between 11% and 13%.This means that most corn acres in the US were cultivated in the Midwest, a region with relatively homogeneous physical environmental and cropping patterns.The increase in corn cultivation between 2002 and 2007 has been widely attributed to biofuel policy (Landis et al., 2008).The model coefficients obtained show that corn proportion is positively related to increased insecticide use in the seven Midwestern states, which supports the analysis of Meehan et al. (2011) using 1997 CDL data.In addition, the model of the remaining 41 states also shows a positive, though less strong, relationship between the proportion of corn acres and insecticide use.Our data support the conclusion that one impact of this increase in corn acres has been an increase in insecticide-treated acres.This also raises concerns about agricultural landscape simplification (monoculture), pest management, and a range of other ecological issues that may be indirectly impacted (Landis et al., 2008;Meehan et al., 2011;Wright & Wimberly, 2013).
Given its toxicity to major corn pests, a general consensus in the literature is that planting Bt corn would reduce pounds of insecticide used per acre (Gassmann et al., 2011).Benbrook (2012) also argues that only a small proportion of corn acres were treated with insecticide prior to Bt corn adoption, which means that the proportion of corn acres treated with insecticide should have been consistently low (.COA data shows the increases of counties with a high proportion of cropland treated with insecticide in the 7 Midwestern states (Figure 4).This supports Benbrook's statement that even in the early stage (1997) of Bt adoption, corn production was not characterized by a high proportion of insecticide treatment.Bt resistance at the field-level has also been documented (Landis et al., 2008b).Fausti et al. (2012) analyzed insecticide use data in corn-producing states with high Bt corn adoption rates, and suggests that farmers might have increased insecticide use on corn fields where corn rootworm began evolving resistance to Bt corn in the mid-2000s.The timing of the increase in corn-derived ethanol demand for biofuel production, has resulted in Bt corn accounting for most of the new planted acreage.With the evolution of Bt resistance, it is widely expected that corn production will continue to require increasing insecticide use as corn takes over other cropland and natural habitats (Wright & Wimberly, 2013).The corn acreage analyzed in this study includes sweet corn, which consists of primarily non-Bt varieties.However, since it is mainly produced in warmer states including Florida and California, and has high value and cosmetic requirement for fresh market, sweet corn also is also expected to be characterized by relatively high insecticide applications.
Unlike the correlation with corn cropland proportion, the relationship of SDI and insecticide use differ between the 7-state and 41-state models: while county-level crop diversity was negatively related to an increase in insecticide use in the Midwestern states, crop diversity in the other states was positively related to higher insecticide use.While this is contrary to the general consensus that landscape simplification should lead to increasing insecticide use, increasing crop diversity is also expected to affect a variety of factors which may tend to either increase or decrease pest pressure.For example, while increasing complexity decreases patch size it also increases areas of refugia where pests (particularly those with a large number of crop and non-crop host plants) may survive insecticide treatments and subsequently re-colonize fields.

Conclusions
The interactions among organisms and their environments is complex and viable IPM strategies for commercial-scale agricultural production in the US have not yet developed (Schellhorn et al., 2015;Tscharntke et al., 2007).The conversion of natural habitats to agricultural use and the intensification of crop systems have reduced biodiversity in both native plants and crops.The pattern of plant functional traits, however, can vary dramatically depending on the crop species composition (Lin et al., 2011).Although increasing crop diversity at the landscape level could ameliorate losses of natural plant functional traits, this system still needs further investigation.The level of insecticide use in each county and year is likely driven by (1) crop value, (2) pest pressure (and the various factors that affect it such as length of growing season, resistance, natural enemies, and others), and (3) other social/economic factors that affect growers' pest management decisions.California, Washington, and Florida, for example, are areas with high crop diversity which include higher proportions of fresh-market vegetables and fruits that have high cosmetic integrity requirements and market values.However, one or more factors may be the dominant drivers locally where their influence exceeds that of the other factors.
In the continental US, regions with a high proportion of cropland treated with insecticides at the county-level were distributed mainly in the coastal areas, southern states and Midwestern states.Insecticide use in the southeast was generally higher than in the northeast, probably due to higher temperatures; and higher in the West Coast relative to similar latitudes on the East Coast, probably due to higher crop values.The coastal and southern areas are also characterized by high Simpson's diversity index, and the Midwest, on the contrary, is characterized by large-scale monoculture with a high proportion of corn.Our spatial panel model demonstrates these patterns quantitatively and reveals the relationship between proportion of corn acres and higher insecticide use.Understanding the ways in which crop and overall landscape diversity and factors other than crop value and pest pressure affect insecticide use are key factors in developing effective pest management strategies.

Figure 1 .
Figure 1.Proportion of county-level cropland treated with insecticide.(A) The average in the 7 Midwestern states has increased during each of the four census periods.(B) The average of all county-level regions in contiguous USA is lower than the Midwest average, and has also increased

Figure 2a .
Figure 2a.Distribution of panel model variables (1997).All values are based on NASS census data.(A) Proportion of county cropland that was treated with insecticides; (B) County level Simpson's diversity index; (C) Annual growing degree days (GDD); (D) Crop market value (total value divided by harvested acres acres as given in Census, and converted to hectares) in US dollars/hectare at county level; (E) Proportion of corn acreages of county cropland

Figure 3 .
Figure 3. Scatter plot of Simpson's diversity Index (SDI) based on Census of Agriculture (COA) data compared with landscape simplification measures based on CDL maps.(A) The relationship between the county level SDI and percentage of like adjacencies (PLADJ) index.(B) The relationship between the county level SDI and core area percent of landscape (CPLAND)

Figure 4 .
Figure 4. Scatter plot of Simpson's diversity Index (SDI) based on Census of Agriculture (COA) data and the proportion of corn acreage at the county level for the continental states

Table 1A
Note.Prop.Corn = proportion of corn in the county cropland.

Table 2A .
Coefficients and estimated impacts from spatial panel model (41 States: individual effects) using either SDI (A) and proportion of corn acreage (B) Note.See Table1footnote for abbreviations.Number of counties in the 41 States: 1,396.