Managing the evolution of herbicide resistance.

BACKGROUND: Understanding and managing the evolutionary responses of pests and pathogens to control efforts is essential to human health and survival. Herbicide-resistant (HR) weeds undermine agricultural sustainability, productivity and profitability, yet the epidemiology of resistance evolution - particularly at landscape scales - is poorly understood. We studied glyphosate resistance in a major agricultural weed, Amaranthus tuberculatus (common waterhemp), using landscape, weed and management data from 105 central Illinois grain farms, including over 500 site-years of herbicide application records.
RESULTS: Glyphosate-resistant (GR) A. tuberculatus occurrence was greatest in fields with frequent glyphosate applications, high annual rates of herbicide mechanism of action (MOA) turnover and few MOAs field(-1) year(-1) . Combining herbicide MOAs at the time of application by herbicide mixing reduced the likelihood of GR A. tuberculatus.
CONCLUSIONS: These findings illustrate the importance of examining large-scale evolutionary processes at relevant spatial scales. Although measures such as herbicide mixing may delay GR or other HR weed traits, they are unlikely to prevent them. Long-term weed management will require truly diversified management practices that minimize selection for herbicide resistance traits.

INTRODUCTION

Widespread, rapid evolution of herbicide‐resistant (HR) weeds is destabilizing weed management in commercial agriculture. This problem has reached epidemic levels, driving up crop production costs, decreasing farm profitability,1, 2 and forcing farmers out of business in some cases (Norsworthy JK, private communication, 2014). Current efforts to impede the evolution and spread of HR weeds focus on diversifying the herbicide mechanisms of action (MOAs) (i.e. herbicidal active chemistries that affect and disrupt specific target sites or biochemical processes3) to which weed populations are exposed.4, 5 Herbicide rotations and sequences (also called cycling) reduce selective pressure on individual target sites by employing herbicides with different MOAs in successive growing seasons or within the same season, respectively, whereas herbicide mixing exposes weeds to multiple MOAs simultaneously.4 Drug mixture and rotation strategies used to combat drug‐resistant infectious agents are based on the same principles; the evolution of drug and herbicide resistance are parallel evolutionary processes.6

Although herbicide rotation and mixing are commonly recommended for preventing HR weeds, almost no quantitative data exist on the effectiveness of either strategy under conditions relevant to production farming7 (see the supporting information). Current understanding is derived from simulation models, laboratory studies and small‐scale‐plot studies focused on proximate instead of ultimate causes of resistance. There is a large gap between the urgency and scale of the problem and the scope of relevant knowledge. Given the epidemic proportions of the problem, we believe it is time to take an epidemiological approach to understanding it.

Our goal was to identify risk factors associated with resistance to the broad‐spectrum herbicide glyphosate in Amaranthus tuberculatus (Moq.) Sauer var. rudis (Sauer) Costea & Tardif (common waterhemp) in field crops (see the supporting information). In recent decades, A. tuberculatus has expanded from its native range in Midwest US riparian areas to become a highly competitive weed in arable systems of this region and beyond. It also has become emblematic of HR weeds in general;8 numerous populations of A. tuberculatus possess multiple herbicide resistance mechanisms that, in the most serious cases, severely limit the number of remaining effective herbicide options.9

Herbicide resistance evolves within a complex of environmental and management factors acting within agroecosystems. To detect risk factors for herbicide resistance evolution that apply broadly to real‐world production systems, as opposed to narrowly defined experimental systems, we took an epidemiological approach. We compiled A. tuberculatus glyphosate resistance frequency, soil, landscape, farm management and weed data, and almost 500 site‐years of field‐level herbicide application records from 105 central Illinois grain farms managed by a single custom retail applicator from 2004 to 2010 (Fig. 1). In addition to deciphering environmental and management signals linked to herbicide resistance evolution, our analysis examined putative causal relationships between past herbicide use and contemporary glyphosate resistance frequencies. Because producers have questioned whether they will be able to manage resistance effectively if their neighbors do not, we addressed spatial contagion in a separate analysis.

Figure 1

Study field relative locations within Fayette and Effingham counties, Illinois. A more detailed map is available in supporting information Fig. S1.

MATERIALS AND METHODS

Data collection and processing

Fields selected for inclusion in the study were identified by a commercial custom applicator that provides services for several hundred fields in central Illinois. The applicator released field‐specific management histories for 2004–2011, and provided access for our team to collect seeds, soil samples and other site data. The fields in the study were distributed over an area of approximately 300 miles2 (777 km2). While small on a continental scale, this is an area of high strategic importance and provided the opportunity to study glyphosate resistance evolution in waterhemp as it was unfolding.

We quantified the glyphosate resistance frequencies of A. tuberculatus at each field in the greenhouse as the proportion of plants grown from field‐collected seeds that survived glyphosate exposure. Two people each sampled seeds from ten plants per field (N = 20 plants field−1) in September 2010. The total seed produced by each sample plant was added to one bulk seed sample per field. Thus, the sample was representative of the frequency of glyphosate resistance traits contributed to the soil seedbank in 2010 rather than the frequency of those traits in the sampled parent plants. We sampled plants along transects according to the following level of priority: infield (highest), head‐row (intermediate), field border (lowest). The majority of plants were sampled from within the fields. All plants sampled, including the minority that came from field borders, were exposed to any pollen‐bearing plants within the field.10 Because we sampled seeds at the end of the growing season, the sample weed populations were potentially enriched for resistance traits against the MOAs that had already been used that year (i.e. plants resistant to the herbicides applied were most likely to survive and produce seed for collection). This is a potential source of bias, but this method allowed us to capture a representative sample of the seeds contributing to the next year's weed population. Seeds were threshed, pooled from within each field and stratified at 4 °C for at least 6 weeks to break dormancy. Seeds were germinated on moistened filter paper in petri dishes incubated in a growth chamber set to 12 h days at 35 °C/15 °C day/night temperatures. After 48 h, germinated seedlings were transplanted to a growth medium of 3:1:1:1 commercial potting mix:soil:peat:sand in 3.8 cm × 21 cm cone‐tainers. Seedlings were grown in a greenhouse with 16 h days at 28–30 °C/24–26 °C day/night temperatures, and with supplemental lighting to maintain a minimum of 800 µmol m−2 s−1. Twenty‐one uniform plants 4–5 cm tall were sprayed with a commercial formulation of glyphosate at 1260 g ae ha−1 using a moving‐nozzle spray chamber.9 This glyphosate rate was chosen on the basis of preliminary experiments to distinguish effectively known glyphosate‐resistant (GR) and sensitive populations in a high‐throughput greenhouse screen. Plants were visually evaluated for glyphosate resistance 14 days after treatment and rated as sensitive (dead) or resistant (alive). A. tuberculatus populations previously determined to be sensitive or resistant to glyphosate were included in all herbicide application runs as controls, and most field populations were evaluated in at least two separate runs (yielding a sample of 42 plants field−1 for most fields). Glyphosate resistance frequencies in each field were estimated using a logistic regression with site as a fixed effect and the number of resistant versus sensitive plants per field as the binomial response.

Data used as covariates in the analyses were compiled from several sources. Farm management records were obtained from the custom retail applicator. Landscape attributes, including field area, field perimeter length, distance to nearest road, field margin composition, presence of and proximity to waterways, etc., were obtained from analyses of aerial imagery in Google Earth. Weed community composition was assessed and soil samples collected in April 2010. Soil samples were analyzed for chemical and physical properties by Great Lakes Laboratories.

To quantify within‐year herbicide diversity, we calculated the total unique herbicide MOAs year−1, as well as the mean and maximum MOAs application−1 year−1, each averaged across years. The distinctions among these are important. The mean MOAs year−1 measures average herbicide diversity per growing season. In contrast, the per‐application metrics quantify the diversity of herbicides applied simultaneously through herbicide mixing. We calculated the proportion of consecutive 2 year sequences in which glyphosate was applied in both years (i.e. a ‘run’ of glyphosate years), the overall proportion of years that glyphosate was applied and the mean number of glyphosate applications per year. Herbicide MOA temporal diversity was quantified with Whittaker's β W,11 Harrison's β H1,12 and our own index of herbicide turnover across growing seasons (‘herbicide turnover index’, based on Diamond and May's13 and Wilson and Shmida's14 indices of species turnover and β diversity). Turnover is the time‐weighted mean proportion of herbicide MOAs that change through gains or losses between data‐years for each field. It is calculated as

where n is the number of data‐years for a field, n − 1 is therefore the number of 2 year sequences of data for the field, t is the index of the 2 year data sequence (e.g. 2004 and 2005, or 2004 and 2006 if data from 2005 are absent), g and l are the number of herbicide MOAs gained and lost during sequence t respectively, S is the total number of MOA between both years, and Δt is the number of years between year 1 and year 2 in data sequence t. All metrics of herbicide MOA diversity were calculated both for the years 2004–2006 within the 55 fields with complete data and in the 76 fields with at least 3 years of data from 2004 to 2010. Our metrics of MOA mixtures, β diversity and turnover weighted all MOAs equally and were not specific to inclusion of glyphosate. Thus, the concepts of herbicide rotation and turnover that we employed did not consider whether glyphosate or any other particular MOA was rotated from year to year. They were based solely on the number or proportion of MOAs that changed annually.

Statistical analyses

Classification and regression tree analyses

We used classification and regression trees (CART)15, 16 to identify potentially important relationships between GR A. tuberculatus presence or proportional frequency and 66 environmental, soil, landscape, weed community and management variables (supporting information Table S1). CART is a non‐parametric machine learning technique that recursively partitions a response variable into increasingly homogeneous subgroups using both continuous and categorical independent variables as splitting criteria. Variables that most effectively minimize the heterogeneity of the response at each node are retained in the final tree, which is then pruned using the one‐SE rule15 to maximize parsimony. The tenfold cross‐validated CART fitting procedure was bootstrapped using 30 repeats. This protected against selecting local model solutions based on initial conditions in the recursive partitioning algorithm. See the documentation for the CARET package in R for more details.17 Our use of CART helped to focus our attention on herbicide diversification and rotation practices in subsequent analyses. The CART analyses of management variables used 76 fields with at least 3 years of management records. All other CART analyses used 105 fields. Comparisons among models with different sets of independent variables were based on the percentage of variance explained by each model relative to the null, single‐node model with zero splits. The management‐based models for both proportional and binary resistance response variables had the highest R 2 values and thus were pursued in the subsequent logistic and GLMM model ranking analyses.

Model building and ranking

Because the CART models pointed to herbicide diversification practices as important, we hypothesized that management prior to the first detection of GR A. tuberculatus would be predictive of the presence/absence or prevalence of resistance years later. We fitted and ranked generalized linear mixed models (GLMMs) of resistance probability (glyphosate resistance presence/absence field−1) and the proportion of GR plants site−1 separately for eight measures of herbicide diversity. The GLMMs had binomial errors and random slopes and intercepts within spatial clusters of sites to control for spatial correlations. We then ranked the models using Akaike weights (w) and AICc18 to identify the predictors and random effects (or lack of random effects) best supported by the data for each response variable (proportion of resistant plants or probability of resistance 4–6 years later). Random slopes and intercepts were retained in the proportional resistance models, but no random effects were supported in the binary (probability of resistance) models, which were thus reduced to generalized linear models (logistic regressions). In these models we used data from 40 sites with complete management records for the period 2004–2006, within 5 km of the nearest neighboring site, and identified using hierarchical clustering as belonging to a spatial cluster of three or more sites.

To verify that these results were independent of proximity effects, we calculated pairwise differences in glyphosate resistance (ΔGR), maximum MOAs application−1 (2004–2006 data) (ΔMOA) and intersite distance for all site pairs. We then used semi‐partial correlations19 to quantify the variance in ΔGR uniquely explained by (a) ΔMOA after controlling for distance effects and (b) distance after controlling for ΔMOA. See the supporting information for further details on model fitting and semi‐partial correlations. The data reported in this paper and the analytical code necessary to reproduce the analyses are archived in the supporting information.

RESULTS AND DISCUSSION

Management practices were most predictive of glyphosate resistance among 66 environmental, soil, landscape, weed and management variables analyzed (supporting information Table S2). Frequency of glyphosate resistance among seeds produced within a given field was highest (supporting information Table S2; Fig. 2) in fields with at least 57% turnover in herbicide MOAs year−1 (i.e. rotation) and at least 0.93 glyphosate applications year−1 (mean GR = 18% of plants field−1). Likewise, the presence of GR A. tuberculatus in fields was associated with high annual turnover in herbicide MOAs, using fewer herbicide MOAs year−1 and using glyphosate in over 75% of years (supporting information Table S2; Fig. 3). Contrary to prevailing wisdom, both binary and proportional models indicated that herbicide turnover increased the frequency of resistance. Logistic regression models supported this result, with positive or non‐significant effects of several metrics of herbicide MOA turnover and temporal β diversity on proportion of GR plants (Fig. 4).

Figure 2

Final regression tree of proportional herbicide resistance. R 2 = 0.27. HR is the mean proportion of herbicide‐resistant plants per field at each node. The percentage of fields at each node is shown below the node. The bold text indicates the splitting criteria at each node. Fields where the criterion is true are moved down the tree to the right to the next node or terminal leaf. For example, the first split indicates that if the herbicide turnover index is at least 0.57, move to the right. If it is less than 0.57, move to the left. Observations were weighted by the number of years of management data available. N = 76 fields.

Figure 3

Final classification tree of the presence/absence of the resistance trait within the field. R 2 = 0.53, accuracy = 0.60. At each node, HR is the mean probability of the glyphosate resistance trait occurring within a field, followed by the percentage of fields represented at the node or leaf. For example, the first node includes 100% of the farms in the dataset, of which 49% have glyphosate resistance. The first split is based on manure use. The label indicates that if manure is applied, go to the left. The leaf to the left includes 18% of the fields in the dataset that together have a mean resistance probability of 0.15. Note: for all other splits, if the splitting criterion is true, observations are included in the branch or leaf to the right of the split. Observations were weighted by the number of years of management data available. N = 105 fields.

Figure 4

Logistic regressions of the per capita probability of glyphosate resistance (proportion of plants resistant) versus three indices of herbicide MOA β diversity. Each metric was calculated for the 3 year period 2004–2006 (n = 55 sites; left column) and for all available years for sites with at least 3 years of data (n = 76 sites; right column). The three metrics are Harrison's β 1 (β H1), Whittaker's β (β W) and our own herbicide turnover index (see methods).

Herbicide mixing was strongly linked with reduced selection for glyphosate resistance. A field with a mean herbicide complexity of 2.5 MOAs application−1 in 2004–2006 was 83 times less likely to produce GR A. tuberculatus seeds 4–6 years later than one with a mean complexity of 1.5 MOAs application−1 (Fig. 5A) (odds ratio = 0.012; ΔAICc = 4.2; w = 0.83). Similarly, GR seed relative abundance decreased as the maximum MOAs application−1 year−1 increased after controlling for random spatial effects (Fig. 5B) (odds ratio = 0.0195; ΔAICc = 16.0; w > 0.99). A. tuberculatus seeds in fields with maximum annual application mixtures of three MOAs were 51 times less likely to be glyphosate resistant than those from fields with two MOAs per application. Model rankings are shown in supporting information Tables S3 and S4.

Figure 5

GR seed production of A. tuberculatus in 2010 versus herbicide mixture complexity prior to detection of resistance in Illinois. (A) Probability that resistance > 0 versus the mean MOAs application−1 year−1. The solid black curve shows the predicted probabilities from a logistic regression with resistance status as the binary response variable. The dashed line shows the predictions excluding sites with MOAs > 2.5 from the analysis. (B) Proportion of GR seed produced per field versus the maximum MOAs application−1 year−1. Colors correspond to cluster identities in supporting information Fig. S1. N = 40 in both panels.

After partialling out distance effects, ΔMOA explained 19.4% of the variance in ΔGR (r = −0.440, P < 0.0001, N = 780), while distance explained less than 1% of the variance after controlling for ΔMOA (r = −0.035, P = 0.3235, N = 780). It is most probable that management differences rather than spatial contagion were the primary drivers of differences in glyphosate resistance frequencies between fields. This suggests that individuals may be able to mitigate GR A. tuberculatus through aggressive management even if nearby fields become highly resistant.

Many variables that were not retained in our best‐supported models quantified the environmental context in which glyphosate resistance evolved (supporting information Table S2). These included the proximity of the nearest field ‘infected’ by glyphosate resistance, the presence of adjacent watercourses (a common route for A. tuberculatus seed dispersal) and waterhemp population density. The strong management signal that we detected in our analysis was context‐invariant within the study region, the epicenter of A. tuberculatus glyphosate resistance evolution in the northern Corn Belt. Use of herbicide mixtures was effective in limiting glyphosate resistance regardless of spatial contagion or weed community, soil or landscape characteristics. We could not have evaluated the relative strengths of these signals without using a data‐mining approach in an epidemiological context.

The central tenets of pesticide (and antibiotic) resistance management are that (a) increased exposure to a selective agent elevates the probability that resistance to that agent will evolve and become enriched,4 and that (b) resistance alleles carry fitness costs that lead to their depletion in the absence of the selective agent.20 Rotating through mechanisms of action reduces exposure to each agent individually, but depletion will only occur if the fitness costs of resistance are high relative to the rotation interval. Although some resistance traits do have high fitness costs (e.g. triazine resistance in A. hybridus and Chenopodium album 21), glyphosate resistance in the related weed A. palmeri comes with little or no fitness cost.22, 23 Our data suggest that glyphosate resistance costs are likely too low for MOA rotation to reduce GR in A. tuberculatus. Because A. tuberculatus forms large, dormant seedbanks, standing allele frequencies for GR will be maintained over time and enriched each time glyphosate is applied.

Herbicide mixes will delay resistance evolution, temporarily

Mixing herbicide MOAs was effective against GR A. tuberculatus because mixing strategies are not reliant on fitness costs to drive depletion. Rather, mixing depletes any resistance alleles by decreasing survival probabilities of all individuals. This strategy will work only if multiple mix components individually suppress a given weed species. One challenge moving forward is finding effective mixes against weeds that have already evolved resistance to many of the previously effective herbicide options. This will be true even after new crops are brought to market that contain multiple herbicide tolerance traits. Another pitfall of mixing is that it may increase the potential for selecting resistance mechanisms that confer broad cross‐resistance, particularly when herbicides are applied at reduced or marginal rates.24 Most GR A. tuberculatus (including samples in the study area that were evaluated) contain EPSPS gene amplification as the resistance mechanism.8 Because amplification of the EPSPS gene specifically affects glyphosate, this type of resistance is probably less likely to arise from selection by herbicides with different mechanisms of action than, for example, metabolic‐based resistance.7 These caveats aside, our results demonstrate that mixing was an effective resistance mitigation strategy for GR A. tuberculatus in a large area located within the heart of US corn/soybean production.

One obvious concern about using herbicide mixes is the increased quantity of herbicidal compounds required. The approach could increase chemical costs and exposure for growers and applicators and elevate risks of environmental damage. Some have advocated that mixed herbicides be applied at lower rates than when applied individually,4 but low‐dose mixtures can potentially increase the risk of non‐target‐site resistance and cross‐resistance evolution.24, 25 Although our analysis does not specifically address application rates, we suggest caution when considering application rate reductions for herbicide mixtures.

Herbicide mixes are not a permanent solution to the problem of target‐site resistance; herbicidal mixtures may delay evolution of resistance, but they do not prevent it.4, 6, 24 HR weed evolution is inevitable, enabled by the current, nearly exclusive reliance on herbicides and the greatly reduced discovery and commercialization of new herbicide chemistries over the last two decades. Long‐term, cost‐effective, environmentally sound weed management will require truly diversified management practices that minimize selection for herbicide resistance traits.5 Combining chemical, cultural, physical and biological tactics can provide cost‐effective weed management while reducing reliance on herbicides.26

The striking parallels between pesticide‐resistant weeds and insects and drug‐resistant pathogens stem from their identical underlying evolutionary causes, which explains the similarities in their responses to management.6 This is why, for example, antibiotic cycling in hospitals is unlikely to suppress evolution and spread of antibiotic‐resistant microbes,27 while mixtures or ‘cocktails’ of multiple drugs administered simultaneously can suppress HIV and drug‐resistant bacterial infections.6, 28 These measures work because the probability of target‐site resistance to multiple drug MOAs is the product of the individual resistance probabilities. By killing all bacterial or viral particles in the host, the evolution of resistance traits is halted.6 Likewise, MOA rotation strategies can fail by allowing resistance to evolve if the cost of resistance is low relative to the rotation interval.27 If the fitness cost of a resistance trait is high, depletion will occur rapidly in the absence of the selective agent. If it is small, depletion will take a long time. In a wind‐pollinated weed such as A. tuberculatus, low‐cost resistance traits can persist in individuals growing outside crop fields (e.g. in field margins) with minimal depletion and thus can survive a multiyear herbicide MOA rotation schedule. Even if the other agents in the rotation are 100% effective, the resistance trait can reinoculate the field.

We have a high level of confidence in our conclusions, given (a) the broad inference domain of our study, (b) the consistent management signals we detected through noisy, real‐world data and (c) that in an analogous problem from the field of health care the same strategy of administering multiple MOA drug mixtures has repeatedly proven to be effective at managing drug‐resistant infections.6 We will encounter resistance evolution repeatedly in natural systems managed for human benefit. Sustainable stewardship of these systems will depend on recognizing that we are always applying selective pressures, and that management responses need to grow from our understanding of applied evolution.

Supplementary Materials

Click here for additional data file.

Table S1. Variables used in CART analyses

Click here for additional data file.

Table S2. CART Model Selection

Click here for additional data file.

Table S3. Generalized linear model rankings for proportional resistance

Click here for additional data file.

Table S4. Generalized linear model rankings for binary resistance

Click here for additional data file.

Figure S1. Study field locations within Fayette and Effingham Counties, Illinois, USA.

Click here for additional data file.

Code Supplement S1. Complete data and analytical code for the analyses reported here are available for download. To access them, download the compressed zip file linked to the article online, unzip the folder, and follow the directions in the README.txt file.

Click here for additional data file.