Occurrence of an herbicide-resistant plant trait in agricultural field margins.

Agricultural environments allow study of evolutionary change in plants. An example of evolution within agroecological systems is the selection for resistance to the herbicide glyphosate within the weed, Conyza canadensis. Changes in survivorship and reproduction associated with the development of glyphosate resistance (GR) may impact fitness and influence the frequency of occurrence of the GR trait. We hypothesized that site characteristics and history would affect the occurrence of GR C. canadensis in field margins. We surveyed GR occurrence in field margins and asked whether there were correlations between GR occurrence and location, crop rotation, GR crop trait rotation, crop type, use of tillage, and the diversity of herbicides used. In a field experiment, we hypothesized that there would be no difference in fitness between GR and glyphosate-susceptible (GS) plants. We asked whether there were differences in survivorship, phenology, reproduction, and herbivory between 2 GR and 2 GS populations of C. canadensis in agrestal and ruderal habitats. We found that geographic location was an important factor in the occurrence of GR C. canadensis in field margins. Although not consistently associated with either glyphosate resistance or glyphosate susceptibility, there were differences in phenology, survivorship, and herbivory among biotypes of C. canadensis. We found equal or greater fitness in GR biotypes, compared to GS biotypes, and GR plants were present in field margins. Field margins or ruderal habitats may provide refugia for GR C. canadensis, allowing reproduction and further selection to occur as seeds recolonize the agrestal habitat. Agricultural practices may select for ecological changes that feed back into the evolution of plants in ruderal habitats.

Introduction

Herbicides impose selection pressure, influencing the evolution of plant species in managed landscapes. Anthropogenically mediated selection pressures are associated with examples of contemporary evolution and allow investigation of ecological feedbacks (Palumbi 2001; Ashley et al. 2003; Rice and Emery 2003; Carroll et al. 2005, 2007; Hendry et al. 2006; Strauss et al. 2006). Such contemporary evolution can occur rapidly and feedback to influence ecological processes (Shefferson and Salguero‐Gómez (2015). A high‐profile example of contemporary evolution has been the increase in herbicide resistance in response to the selection pressures of agricultural management systems (Neve et al. 2014). Insight into the evolutionary dynamics of plants may come from understanding the adaptive evolution of herbicide‐resistant weeds (Murray et al. 2012). This study explores the ecology of herbicide resistance in the absence of selection pressure. We ask what happens when glyphosate‐resistant (GR) weeds move out of agricultural fields into surrounding nonagricultural settings. Modern agriculture is dependent upon glyphosate (Young 2006), and there are 32 species of GR weeds worldwide (Heap 2015). The ecology and evolution of many of these GR biotypes in the absence of herbicide selection pressure are unknown.

The ecological implications and dynamics of glyphosate‐resistant populations depend upon the fitness of evolved GR biotypes. If there is a fitness cost associated with resistance, removal of the selection pressure of herbicide may allow the regression of resistance (Prather et al. 2000). However, if a plant only expresses resistance in the presence of the herbicide, such as the upregulation of genes responsible for resistance, then no change in fitness may occur in the absence of the herbicide or selection pressure. The costs of resistance on fitness may be difficult to detect (Bergelson and Purrington 1996; Hails 2000). Mixed results from hundreds of studies attempting to measure the cost of resistance in the absence of selection pressure may be partially attributable to genetic linkages between genes responsible for compensatory fitness effects and genes responsible for resistance (Bergelson and Purrington 1996). Mixed results may also be a result of limited experimental design and small sample sizes and may be avoided by following robust design recommendations (Vila‐Aiub et al. 2009). Also, the trade‐offs involved in fitness costs can have complex ecological interactions (Strauss et al. 2002).

In order to examine the potential of GR biotypes to persist and feedback into further evolutionary processes, a large‐scale survey of glyphosate resistance was undertaken in the field margins of soybean fields with the GR crop trait, relating the GR trait in C. canadensis to field environmental factors. This study tested the hypothesis (Ho1) that the occurrence of GR C. canadensis in field margins was correlated with environmental factors associated with previous management in the field interiors. The long distance dispersal ability of C. canadensis may decrease any correlation between field history of herbicide use and the ability to predict occurrence and dynamics of GR populations within the field or margins (Davis et al. 2009b). Typically, however, 99% of seed fall 100 m or less from the mother plant (Dauer et al. 2007). Other research has looked at local land use in relation to the occurrence of resistant C. canadensis in perennial crop and noncrop areas and found little or no relationship (Hanson et al. 2009). However, other work found that management factors, such as crop rotation, tillage, soil properties, crop residue cover, and geography, were associated with the occurrence of GR C. canadensis within crop fields (Davis et al. 2009b). To further investigate the potential persistence and fitness of GR biotypes, we conducted an experimental field study, which focused on survival and reproduction, testing the hypothesis (Ho2) that there would be no difference in fitness between GR and glyphosate‐susceptible (GS) C. canadensis. We asked whether there were differences in survivorship, phenology, reproduction, and herbivory of two populations of GR and two populations of GS C. canadensis when each was planted in a ruderal, old‐field habitat and an agrestal, soybean habitat.

Materials and Methods

Study species

Conyza canadensis L. Cronq. (formerly Erigeron canadensis L., common names horseweed, marestail, and Canada fleabane) was the first eudicot species to display glyphosate resistance (GR) and has evolved resistance on multiple, independent occasions (Yuan et al. 2010; Okada et al. 2013). To date, 37 GR biotypes have been identified, more than any other GR weed (Heap 2015). Conyza canadensis may reduce crop yields by 90% in some cases (Bruce and Kells 1990), has a wide geographic distribution from latitudes N 55° to S 45°, and tolerates a wide range of conditions (Weaver 2001). The plant is native to North America, although it is now globally dispersed. Within the last 350 years since its introduction to Europe, the plant has become naturalized and abundant (Thébaud and Abbott 1995).

Previous studies on the fitness of GR C. canadensis have been conflicting and inconclusive (Zelaya et al. 2004; Shrestha et al. 2007; Davis et al. 2009a; Shrestha et al. 2010; Grantz et al. 2008). It is unclear whether or not any differences in fitness may be related to the GR trait. The mechanism(s) of resistance is(are) not well understood (Shaner 2009), although evidence implicates a temperature‐dependent process involving a tonoplast‐membrane pump, which sequesters glyphosate into the central plant vacuole (Ge et al. 2010, 2011). There may also be other changes in the way that xenobiotics or plant secondary compounds are transported or stored in GR C. canadensis as compared to GS plants.

Regional survey (Tests of Ho1)

Field sites and survey methods

This study uses field sites and data collected as part of the Benchmark Study, a large‐scale field study of the management practices and impacts of cropping systems using the GR crop trait (Givens et al. 2009, 2011; Owen et al. 2011; Shaw et al. 2011; Weirich et al. 2011a,b; Wilson et al. 2011; Gibson et al. 2013). The Benchmark Study is a long‐term multistate project that seeks to collect data on weed dynamics from 156 sites from six US states (Illinois, Indiana, Iowa, Mississippi, Nebraska, and North Carolina) in GR crop systems to advance stewardship of the GR crop trait technology. Within the Benchmark Study, environmental factors, relevant for the success of managing crop systems with the GR crop trait (crop type – corn, cotton, and soy; crop rotation; rotation of the GR crop trait; herbicide diversity; and tillage), were identified and studied beginning in 2006. Fields were divided in half and managed according to the grower's (AG1) or the researcher's (AG2) recommendations to provide a comparison of methods in order to determine the sustainable practices for minimizing weed population shifts and the development of herbicide resistance (Wilson et al. 2011; Gibson et al. 2015). Detailed methods can be found in Shaw et al. (2011).

Seed collection

Seeds were collected from Benchmark Study field margins at 17 sites in three US states in order to examine the occurrence of the GR trait and the relationship to the selection pressures in field interiors (Fig. 1). These sites were chosen based upon 2006 field survey data that showed the presence of C. canadensis in field interiors. Seeds were collected from up to 11 plants in field margins which had not been exposed to herbicide spray or drift (Table S1). Seeds from each plant were kept separate so that the inference of GR status could be made for individuals through tests of their progeny. In order to assign a biologically meaningful variable to geographic location, field sites were identified by USDA plant hardiness zone (Cathey 1990). USDA plant hardiness zones reflect the average minimum frost‐free days for a given location and are therefore a direct measure of temperature tolerances for plant distributions that may affect the occurrence of resistance.

Figure 1

Map of field margins indicating sites of seed collection, labeled by percentage of GR C. canadensis plants collected.

Spray trials for resistance screening

Seeds from individual “mother” plants were germinated and grown in the greenhouse for GR screening in January 2009. At the one‐leaf stage, seedlings were transplanted into 10 cm × 10 cm pots, two plants per pot, for a total of eight plants per herbicide treatment. Additional seeds were germinated, and plants were transplanted for a second run. After transplant, plants were fertilized every three days and were grown at a temperature range of 18 to 32°C, with supplemental light providing 120 to 140 μmol m−2 s−1 photosynthetic photon flux for a 16‐h photoperiod.

Resistance screening of greenhouse plants was conducted in a spray chamber. Plants were sprayed with glyphosate (Roundup WeatherMax®, Monsanto Company, 800 North Lindbergh Boulevard, St. Louis, MO) once they reached the 5 cm diameter rosette stage. Herbicide was applied uniformly with a moving nozzle sprayer, delivering 187 L ha−1 at 276 kPa (R and D Sprayers, Opelousas, LA). The herbicide rates used were as follows: 1X the labeled rate, 0.84 kg ae ha−1, and 3X the labeled rate, 2.52 kg ae ha−1. At 7, 14, and 21 days after treatment (DAT), plants were visually assessed for signs of herbicide injury, similar to the methods of Davis et al. (2009a). After 21 DAT visual ratings, plant biomass was harvested at the soil surface, oven dried at 55°C, and weighed. Biomass reductions for the treatments, each averaged by mother plant, were calculated as a percentage of untreated controls. Previously screened populations of known GR and GS plants were also treated and visually rated for reference in the statistical tests.

Quantification of glyphosate resistance

The whole‐plant glyphosate screen data were used to determine the occurrence of GR C. canadensis in 2008 field margins by assigning a resistance score using a principle component analysis (PCA). Each mother plant was also classified as GR or GS based on progeny response to the spray test using a discriminant analysis. All data analysis was performed using SAS® software (Version 9.2, 2002–2008, SAS Institute Inc., Cary, NC). Variables were visual ratings at each DAT and herbicide treatment rate (7 DAT 1X, 7 DAT 3X, 14 DAT 1X, 14 DAT 3X, 21 DAT 1X, and 21 DAT 3X) and percent reduction of dry weight for 1X and 3X treatments, reported for each mother plant per run, as an average of four replicates per run. PCA was computed using a variance–covariance cross‐products matrix, considering that all variables were percentages, of the same scale and magnitude (Naik and Khattree 1996). The number of significant principle components to retain for interpretation and further analysis was determined using parallel analysis (Franklin et al. 1995). Each mother plant was assigned a score for axis loadings within the PCA. PCA scores were averaged by field, and these scores were then used for further analyses as a determination of the level of GR.

The visual rating at 14 DAT for the 1X rate was selected as the most informative classification variable and used in a training set by which to run a discriminant analysis (DA). DA is a type of canonical variates analysis, where a set of linear functions of the predictor variables is constructed and used to classify unknown samples (Gittins 1985). The control plants, two GR and two GS C. canadensis populations, were used as the subjects of the training set. The 14 DAT 1X variable was selected because it had a high component loading for PCA axis 1 and was indicated as a significant predictor of glyphosate resistance for runs 1 and 2 in the DA. Although classification based upon all six of the visual ratings yielded approximately the same results, selection of only one variable avoided the problems of multicollinearity arising from using all the variables in one analysis. Using progeny response to the spray tests, mother plants were categorized as GR or GS. The percentage of GR mother plants per field was quantified for runs 1 and 2. In order to obtain a single quantification of percentage of resistant plants per field for subsequent analysis, the percentage of resistant plants for runs 1 and 2 were averaged and coded as DA percent resistance.

The occurrence of resistance in field margins was tested for correlations with environmental variables (test of Ho1). The two dependent variables, mother plant PCA axis 1 score and DA estimate of percent resistance for plants in 2008 field margins, are referred to collectively in the discussion as measures of occurrence of GR C. canadensis. A mixed model analysis (Littell et al. 2006) was used to test the two dependent variables by association with environmental variables in years 2006 to 2008. The Kenward–Rogers approximation procedure (Littell et al. 2006) was used for determining the correct degrees of freedom for the estimates. Several one‐way mixed model analyses were conducted to overcome the limitations of low sample size (N = 17 fields) on modeling measures of resistance with independent variables: USDA plant hardiness zone (zones 4, 5, and 6), crop rotations (0 – no, 1 – yes), rotation of the GR crop trait (0 – no, 1 – yes), crop type (corn or soybean in 2006, 2007, and 2008), tillage (0 – no‐till, 1 – minimum or conventional), herbicide selection pressure (total number of herbicide applications (low: 1–2 times, high: 3–5 times), number of times glyphosate was applied (low: 0–1, high: 2–3), and number of times herbicides with other modes of action were applied (low: 0–1, high: 2–5) in years 2006, 2007, and 2008. The covariance structure used in the analysis was unstructured or variance components, based on the structure type returning the lowest AIC fit statistic. The relationship between measures of GR C. canadensis resistance in 2008 field margins and continuous variables, for example latitude and longitude, was examined using Pearson's correlation analysis providing a test of our hypothesized relationship (Ho1) between GR status and environmental factors.

Potential fitness costs of GR and GS biotypes (tests of Ho2)

Study site and experimental design

The study examined differences in performance between four populations with known resistance status in agrestal and ruderal habitats providing a test of our hypothesis (Ho2) that there would be no difference in fitness between GR and glyphosate‐susceptible (GS) C. canadensis biotypes. In the summer of 2008, 100 plants from germinated seed of the 2 GR and 2 GS control populations used in the previously described spray test were transplanted into 5‐cm‐diameter tubes at the one‐leaf stage and allowed to establish for two weeks. Individuals were then transplanted one meter apart in a random placement in two habitats, a soybean field (agrestal habitat, 5 rows of 40 plants) and a second‐year successional old field (ruderal habitat, 4 rows of 50 plants), similar to a field margin in plant community and competition, (37.698295°N, −89.242460°W) in June. A week prior to field transplant, soybeans with the GR crop trait were planted in 76‐cm rows using a no‐till planter in a two‐year no‐till soybean field. Fifty C. canadensis seedlings from each population were planted midway between soybean rows. Competing vegetation was pulled from the site within 8 cm of the transplant to facilitate establishment and supplemental water was provided for two weeks. Once seedlings reached a height of approximately 10 cm, the soybean field was sprayed with a backpack sprayer and handheld boom using a 1X labeled rate of 870 g ae ha−1 (22 fl oz/acre) of glyphosate (Roundup WeatherMax®, Monsanto Company) to remove competition and simulate normal management. The seedlings were covered with 470‐ml plastic cups, and herbicide was allowed to dry on surrounding plant leaves before cups were removed. Climate data were monitored using an on‐site weather station (data not presented; WARM 2012).

Following field establishment, surrogate measures of fitness were quantified for individual plants. Plants were surveyed on days 29, 62, 101, 133, 168, 198, 398, and 470 after planting. Data were collected on survivorship, height, rosette diameter, number of leaves, life stage, and herbivory at each survey date. Percentage of leaves with signs of herbivory was calculated, after methods of Girdler and Radtke (2006). Digital images were taken of individual herbivores within 5 cm of leaf damage for later identification. Flowering stalks were covered with polyester mesh bags to prevent seed dispersal. After 80% of capitulae on each individual had matured, or first frost, flowering stalks were collected for further measures. Three capitulae from each flowering individual were randomly selected, and mature seeds were counted. When available, three samples of 1–1.5 mg of seeds from each flowering individual were weighed on a scale accurate to the nearest μg, and a mass per seed value was obtained. Six plants in the ruderal environment were still living at the end of the first growing season, and data collection continued the following season.

Statistical analyses

Statistical analyses were conducted to determine whether there were differences in surrogate measures of fitness between populations and habitats and to test Ho2. All statistical analyses were completed with SAS 9.2 software (SAS Institute Inc., Cary, NC). To determine the differences between habitat and population and habitat by population interactions on survivorship, time to bolting, and time to flowering, Wilcoxon's tests in survivorship (failure time) analysis were used to test for homogeneity between groups. Significance of multiple comparisons was obtained using Z‐statistics (Fox 1993). Phenological stage at the time of death was analyzed with a mixed model analysis to test for differences in population within each habitat. Maximum percentage of herbivory for each individual per individual's lifetime was analyzed using a nonparametric test (Kruskal Wallis H test). Comparisons between populations and habitats were made by calculating Z critical values, taking into account family‐wise error (post hoc test, Siegel and Castellan 1988). The variable leaf number was cube‐root‐transformed to meet statistical assumptions. Plant height, leaf number, and plant diameter were tested for differences between population, day, and population by day interaction, up to day 101. Following day 101, there were too few cases for repeated measures analysis. Least squares means comparisons were made for population by day interactions in the agrestal habitat at days 62 and 101 for plant height, leaf number, and plant diameter. The reproductive measures, capitulae count per plant, average seed count per capitulum were tested, and average weight per seed was tested for differences by population in each habitat. Degrees of freedom were calculated using Satterthwaite or Kenward–Rogers methods, as appropriate (Littell et al. 2006).

Results

Regional survey (Tests of Ho1)

Occurrence of GR in field margins varied across a wide geographic gradient. Field sites ranged from approximately 37.29 to 42.82°N and −101.62 to −85.15°W (Fig. 1; Table S2). Management practices within fields included rotation of crops (65%), rotation of the GR crop trait (35%), and tillage, which varied by year (2006 – 18%, 2007 – 29%, and 2008 – 12%).

Principal component analysis and discriminate analysis

Plants were classified as GR or GS using PCA and DA. PCA axis 1 explained 81% of the variability in the data, while PCA axis 2 explained an additional 10% (data not shown). Parallel analysis (Franklin et al. 1995) determined that only PCA axis 1 scores should be retained for interpretation and further analysis, although axis 2 scores were used to produce scatter plots for visual examination (Fig. 2A). Component loadings for axis 1 indicated that visual ratings at 21 DAT and 14 DAT for the 1X rate treatment had the most influence on percentage of explained variability (Eigenvectors 0.52 and 0.50, respectively). In order of decreasing importance were 14 DAT 3X rate, 7 DAT 1X, 7 DAT 3X, and 21 DAT 3X (Eigenvectors 0.37, 0.35, 0.35, and 0.30). The least explanatory variables were the percentage of dry weight reductions for both 1X and 3X rates (Eigenvectors 0.004 and 0.002). By examining the placement of the known resistant and susceptible reference plants within the scatter plots, it was determined that a more negative score indicated a higher degree of herbicide resistance for mother plants as detected by the progeny spray tests (Fig. 2A).

Figure 2

Occurrence of GR C. canadensis in 2008 field margins by USDA plant hardiness zone and management factors. (A) Principle component analysis (PCA) of data from spray test by field margin. USDA plant hardiness zones 4 and 5 were not different from each other (P = 0.38), but were different from 6 (P = 0.002, P = 0.004, respectively). (B) Mixed model analysis of discriminant analysis (DA) percent resistance for GR C. canadensis in 2008 field margins by USDA plant hardiness zone. (C) The relationship of occurrence of GR C. canadensis in 2008 field margins to tillage (no‐tillage or minimum/conventional tillage) in 2007. (D) The number of herbicide applications with mode of action (MOA) other than glyphosate for the field half managed according to grower practice (AG1) in 2007 by percentage of GR C. canadensis in 2008 field margins. Low herbicide diversity is 0 to 1 applications; high, 2 to 5 applications.

Once the DA was used to classify the progeny as glyphosate resistant or glyphosate susceptible, the percentage of resistant mother plants by field margin was determined to range from 10 to 100% (mean value of 46%) (Table S2). Wilks' lambda indicated a significant difference between the GR and GS classified plants, based upon visual ratings for 14 DAT at the 1X rate (Run 1: Λ = 0.03, F 1,2 = 64.67, P = 0.02; Run 2: Λ = 0.0004, F 1,2 = 4455.68, P = 0.0002).

Mixed model and correlation analyses of dependent variables with GR status (Test of Ho1)

The occurrence of resistance for both dependent variables (DA percentage of GR in field margins and PCA axis 1 scores) was more strongly associated with geographic location than any of the other variables analyzed (Table S3). GR C. canadensis was more common and more variable in field margins at the lowest latitudes below 40° (Fig. S1), as reflected by relationships with PCA axis 1 scores. The latitudinal patterns were reflected in the analysis of USDA plant hardiness zones, with the most southern zone, zone 6, having the greatest occurrence of GR C. canadensis in 2008 field margins, while zones 4 and 5 were not different (Fig. 2B). There was no effect of percent clay, sand, or silt in the soil (data not presented) or the management factors of crop rotation, rotation of GR crop trait, or crop type within the field interiors on occurrence of GR C. canadensis in 2008 field margins, although tillage in 2007 was used more often in field interiors with a greater occurrence of resistance in field margins in 2008 (Fig. 2C). No relationship between in‐field herbicide usage for years 2006 to 2008 and occurrence of GR C. canadensis in 2008 field margins was detected until 2007, and this was only for herbicide diversity (the use of herbicides with alternative modes of action to glyphosate). Response of both dependent variables indicated that occurrence of resistance was greatest in 2008 field margins where more applications of herbicides with alternative modes of action had occurred in the field interior during the previous growing season in 2007 (Fig. 2D).

Potential fitness costs of GR and GS biotypes (Tests of Ho2)

Survivorship, bolting time, and flowering time

There was a significant interaction between habitat type and population for survivorship, bolting, and flowering (Table 1). The log‐rank test showed similar significance (not shown). The resulting survivorship curves were similar to a Deevey (1947) type I curve, where mortality rates increase with increasing age of individuals (Fig. 3A). The GR1 population in the ruderal habitat had the greatest survivorship of all eight population‐by‐habitat combinations (all comparisons: Z > 1.65, P < 0.05, not shown individually). Ruderal populations of GR2, GS1, and GS2 were not different from one another (GR2 vs. GS1: Z = 0.48, P = 0.31; GR2 vs. GS2: Z = 0.002, P = 0.50; GS1 vs. GS2: Z = 0.49, P = 0.31). Survivorship of all populations in the agrestal habitat was similar (all comparisons: Z < 0.79, P > 0.21, not shown individually). Population GR1 in the agrestal habitat had similar survivorship to populations GS2 and GR2 in the ruderal habitat (GR1 vs. GS2: Z = 1.34, P = 0.09; GR1 vs. GR2: Z = 1.32, P = 0.09), and survivorship of the GR2 population did not differ by habitat (GR2 vs. GR2: Z = 1.29, P = 0.10).

Table 1

Failure‐time analysis of survivorship and advancement of phenological stages. Wilcoxon's test statistic (χ ) is reported

	df	χ 2	P
(a) Time to death, N = 400
Habitat	1	19.91	<0.0001
Population	3	4.06	0.25
Habitat × population	7	27.03	0.0003
(b) Time to bolt, N = 400
Habitat	1	69.83	<0.0001
Population	3	0.36	0.95
Habitat × population	7	73.80	<0.0001
(c) Time to flower, N = 400
Habitat	1	9.00	0.003
Population	3	6.41	0.09
Habitat × population	7	19.01	0.01

df, degrees of freedom.

Figure 3

Failure‐time analysis of survivorship, time to bolt, and time to flower of GR and GS populations in agrestal and field margin (ruderal) habitats: (A) Survivorship curve of individuals by population and habitat. (B) Survivorship curve of individuals that did not reach bolting stage by population and habitat. (C) Survivorship curve of individuals that did not reach flowering stage by population and habitat. Abbreviations in the legend indicate plants from the four populations: GR1 – GR population 1, GR2 – GR population 2, GS1 – GS population 1, GS2 – GS population 2. Sampling times: day 1, June 28, 2008; day 29, July 24, 2008; day 62, August 26, 2008; day 101, October 4, 2008; day 133, November 5, 2008; day 168, December 10, 2008; day 198, January 9, 2009; day 398, July 28, 2009; and day 470, October 8, 2009.

Although analysis of time to bolt indicated a significant interaction with stratification by habitat and population, there were no differences between populations in each habitat (within agrestal habitat, all comparisons: Z < 0.77, P > 0.22; within ruderal habitat, all comparisons: Z < 1.23, P > 0.11; between habitats, all comparisons: Z > 3.04, P < 0.001) (Fig. 3B). Bolting occurred later in the ruderal habitat in all populations, and fewer individuals reached the bolting stage. While there were no differences between populations within habitat types for bolting time, there were significant differences in flowering time (Fig. 3C).

Both GR populations in both habitats were observed to flower earlier than the GS populations; however, this observed trend was not significant in the survivorship analysis. Flowering of the two GR populations in the agrestal habitat were different (GR1 vs. GR2: Z = 2.05, P = 0.02). While the occurrence of flowering began at the same time, GR2 had 12% more flowering individuals than GR1 by the end of the study (Fig. GS2). Agrestal GR1 began flowering earlier but had fewer flowering individuals than agrestal GS1 (Z = 2.01, P = 0.02). Agrestal GR1 had a similar number of flowering individuals as agrestal GS2 (Z = 0.40, P = 0.34). Agrestal GR2 flowered earlier and had more individuals reach the flowering stage than agrestal GS2 (Z = 1.65, P = 0.05). Although agrestal GR2 began flowering earlier than ruderal GR2, the final number of flowering individuals was similar, and the resulting survivorship curves were not different (Z = 1.41, P = 0.08). Agrestal GR2 had many more flowering individuals than ruderal GR1 (Z = 2.39, P = 0.008), GS1 (Z = 2.93, P = 0.002), and GS2 (Z = 2.91, P = 0.002). Agrestal GS1 had the most flowering individuals of all other populations in the two habitats, with 22% of individuals reaching the reproductive stage. Agrestal GS1 plants flowered later with many more flowering individuals than agrestal GR1 plants and all ruderal populations except ruderal GR2. Although ruderal GR2 began flowering earlier than agrestal GS1 individuals, the resulting survivorship curves were not different (Z = 1.36, P = 0.09). Despite differences in number of flowering individuals by population, there were no differences between populations for life stage of individuals at the time of death in each habitat (Agrestal: F 3,173 = 0.36, P = 0.78; ruderal: F 3,124 = 0.28, P = 0.84).

Maximum percentage of herbivory

Herbivory damage was quantified to test for differences between populations and habitats (Appendix 1). Data were not normally distributed, and Kruskal–Wallis tests of maximum percentage of herbivory recorded over the entire life span of each individual was different by population in the ruderal habitat (χ 2 = 8.63, df = 3, P = 0.03) but not in the agrestal habitat (χ 2 = 1.83, df = 3, P = 0.61). Post hoc tests indicated that the only difference in percentage of herbivory was between GR1 and GR2, with GR2 experiencing more herbivory than GR1 in the ruderal habitat (Fig. S3).

Reproductive allocation

None of the reproductive measures were significantly different between populations in the ruderal habitat, possibly because of the low number of plants reaching reproductive maturity (number of individuals: GR1 = 3, GR2 = 5, GS1 = 2, GS2 = 2) (Table S4). Of the plants that reached the flowering stage, the GR biotype, GR1, had the highest capitulae count in the agrestal habitat, although not significantly different from GR2 (Fig. 4A). However, once capitulae production was normalized by the total number of plants in each population, including those that did not reproduce, GR2 had almost 1000 more capitulae per plant than GR1 (Fig. 4B). Population GS1 produced among the lowest capitulae counts and the lowest number of seeds per capitulum (Fig. 4C) of all the populations in the agrestal habitat. Although two GS1 plants reached the reproductive stage in the ruderal habitat, no mature capitulae were produced. There was no difference in weight per seed by population in either habitat (agrestal: N = 26, mean = 0.053 mg ± SE 0.002; ruderal: N = 10, mean = 0.021 mg ± 0.008). The result with the greatest implications was produced by calculating the total number of seeds produced by population. The GR2 population produced over 300,000 seeds more than any of the other populations in the agrestal habitat (Fig. 4D).

Figure 4

Comparison of reproductive measures of glyphostate‐resistant (GR) and glyphosate‐susceptible (GS) populations in the agrestal habitat: (A) Capitulae count by number of flowering individuals; number of individuals are GR1 = 3, GR2 = 10, GS1 = 11, GS2 = 6, (B) Capitulae count normalized by total number of individuals, N = 50, (C) Average number (N = 3) of seeds per capitulum per individual; number of individuals are GR1 = 3, GR2 = 9, GS1 = 10, GS2 = 5. Mean values with the same letter are not significantly different (α = 0.05). (D) Total number of seeds produced by population as a derived measure from total number of flowering plants × average number of capitulae per plant × average number of seeds per capitulum. Abbreviations for the four populations are: GR1 – GR population 1, GR2 – GR population 2, GS1 – GS population 1, GS2 – GS population 2.

Perennation

Conyza canadensis may exhibit a biennial form (USDA 2013). However, the overwintering plants within the region of the study sites are thought to germinate in the fall. The plants in this study were germinated from seed and transplanted into the field early in the summer, and therefore, overwintering was unexpected. Of 200 plants in the ruderal habitat, six survived into the winter as rosettes or bolted individuals (GR1: N = 4; GR2: N = 0; GS1: N = 1; GS2: N = 1). The following spring, four individuals remained (GR1: N = 3; GR2: N = 0; GS1: N = 1; GS2: N = 0), and, of those, three reached flowering stage in the fall (all GR1 population, data not presented). Following damage to the apical meristem at first frost, 14 bolting‐stage and one flowering‐stage individuals in the agrestal habitat produced adventitious rosettes at the base of the plant (GR1: N = 7; GR2: N = 4; GS1: N = 2; GS2: N = 2), although the agrestal plants did not survive the winter. To our knowledge, overwintering of plants that have survived a summer season and the occurrence of the production of adventitious rosettes has not been previously documented in the literature.

Discussion

This study suggests that biotypes of GR C. canadensis may persist in the ruderal habitats surrounding agricultural fields. These biotypes have the potential to reproduce and disperse into the agricultural field and feed back into the further evolution of the species. The patchiness and gradients of the agrestal–ruderal habitat, coupled with the refugia of the field margin, may interact with other biotic factors, such as climate and management intensity, and impact the evolution of GR biotypes of C. canadensis.

Our study supports our hypothesis (Ho1) that there is a strong relationship between environmental factors, especially geographic location, and selection pressure on the evolution of GR biotypes. Occurrence of resistance in field margins was related to geographic location and cropping systems utilizing the GR crop trait. The evolution of GR biotypes depended upon environmental factors related to location across a wide geographic area. Other studies have also noted trends in the geographic distribution of in‐field herbicide‐resistant C. canadensis (Barnes et al. 2004; Davis et al. 2008, 2009b; Hanson et al. 2009; Ge et al. 2011; Okada et al. 2013), although there have been no previous analyses of GR C. canadensis in field margins with the GR crop trait. In a survey of 400 growers concerning the prevalence of GR weeds in the northern “corn belt” or the southern “cotton belt” of the USA in the spring of 2006, 24% of northern and 39% of southern growers suspected problems with GR weeds; although 67% of southern growers had grown continuous GR crops for 3 to 5 years, while only 27% of northern growers had used no rotation of the GR trait (Foresman and Glasgow 2008). Foresman and Glasgow (2008) pointed out that it remains to be seen whether or not these differences in distribution are a product of land use and herbicide selection pressure or larger regional effects associated with geography, or both. However, the present study suggests an overarching effect of geography and that changes in management may have been implemented to reduce levels of resistance – an effect rather than the cause of resistance.

Linked to geography, temperature has been identified as a proximal factor in the susceptibility of C. canadensis to glyphosate. One mechanism of resistance is vacuolar sequestration; although GR and GS C. canadensis plants may uptake the same amount of herbicide into cells, twenty‐four hours following application, GR plants had sequestered 85% of glyphosate in vacuoles compared to only 15% in GS plants (Ge et al. 2010, 2011). Sequestration of glyphosate protects the chloroplasts from the herbicide mode of action. However, a brief period of cold acclimation, 12/8°C day and night temperatures for a minimum of 7 days, prior to application may increase the susceptibility of plants that exhibit a high level of resistance in warmer temperatures. The likely transport mechanism is a tonoplast‐membrane pump. The transport of glyphosate across the tonoplast into the vacuole is inhibited in cold conditions, leading to plant death as glyphosate enters the plant cell (Ge et al. 2011). Therefore, patterns of resistance associated with wide geographic range could be a combination of management practices and plant physiological tolerances.

While we did not see a fitness cost associated with GR, supporting our hypothesis (Ho2) that GR does not impose a fitness penalty, a small sample size makes it difficult to draw generalizations. However, previous studies comparing fitness among GS and GR biotypes have also been inconclusive. Populations of GR C. canadensis in California had earlier phenology (Shrestha et al. 2010) and more rapid growth in early phenological stages than GS plants (Shrestha et al. 2007, 2010; Grantz et al. 2008). One of our GR populations showed early phenology; however, other measures of growth and reproduction for this GR population were not different from either of the susceptible populations. Other studies have found few or no differences in growth and reproduction between GR and GS C. canadensis (Zelaya et al. 2004; Davis et al. 2009a).

This study suggests that habitat type (ruderal vs. agrestal) may influence the timing of phenological stages. Although it is advantageous for C. canadensis rosettes to bolt quickly to emerge from the canopy of competing vegetation (Regehr and Bazzaz 1979), competition has previously been linked to delayed onset of bolting (Thébaud et al. 1996). Several studies have documented earlier bolting in GR compared with GS C. canadensis biotypes, and it is unclear whether this change in phenology will persist under high levels of interspecific competition in a ruderal habitat dominated by graminoid species. If rapid bolting leads to early reproduction, reproductive success may be enhanced through avoidance of unfavorable conditions following the peak growing season, that is frost kill (Thébaud et al. 1996). In addition to the effect of altered competitive interactions as characteristic of different habitats, differences in rate of phenological development may influence plant–herbivore interactions (Jarzomski et al. 2000). In Europe, where C. canadensis is a relatively new invader, insect and mollusk herbivory was greatest in habitats with the highest densities of graminoid species (Prieur‐Richard et al. 2000; Prieur‐Richard et al. 2002). Although habitat type may affect the level of herbivory and thus survivorship, differences between biotypes are unlikely, as C. canadensis appears to be relatively herbivore tolerant. Other studies have shown that herbivory does not affect bolting time, flowering time, or reproductive output of C. canadensis (Thébaud et al. 1996; Abhilasha and Joshi 2009). Our study suggests that of the two GR biotypes studied, the GR biotype that experiences greater herbivory also has higher reproductive output.

Adaptive change in terms of phenological change has implications for the spread of GR individuals. It is unclear if early flowering might be a survivorship advantage, but similar shifts in phenology have been seen in other studies of C. canadensis (Alcorta et al. 2011). A dramatic change in phenology may lead to reproductive isolation between biotypes and further evolution. Even without disjunct phenology, C. canadensis is primarily autogamous, and the estimated outcrossing potential is only 4% (Regehr and Bazzaz 1979). However, in a plant that produces up to 200,000 seeds, 4% could significantly affect directional selection within a population under extreme pressure. The GR trait has been shown to spread by pollen (Zelaya et al. 2004), and actual outcrossing rate may be higher than 4% with increased stand density or large populations of insect floral visitors. The significantly higher numbers of seed produced by one of the GR populations in this study and the characteristics of the GR trait in C. canadensis suggest a rapid increase in the occurrence of GR individuals can occur. Additional increases in the frequency of GR individuals may occur with climate change, in the case of warmer spring temperatures at the time of herbicide application or milder winters which could select for perennation. The potential interaction between herbicide resistance‐driven evolution and climate change‐driven evolution (Shaw and Etterson 2012) may magnify phenological shifts.

Conclusions

We conclude that even in the absence of glyphosate, there was no consistent evidence that the GR trait will decline in frequency within populations of C. canadensis. Continued use of glyphosate without additional diversity in land management may select for GR C. canadensis with traits conferring increased fitness over GS plants. In this way, it may be possible that glyphosate resistance could drive other advantageous traits to fixation.

Data accessibility

The permanent archive of the data used in this study is available from the Dryad Digital Repository [to be completed upon publication].

Conflict of Interest

None declared.

Table S1. Additional information about Benchmark Study locations and samples taken from field margins to test for relationships between the occurrence of GR C. canadensis, field history, and location.

Table S2. Descriptive statistics for independent variables, Principal Component Analysis (PCA) axis 1 scores and Discriminant Analysis (DA) percentage of glyphosate‐resistant (GR) C. canadensis plants in field margins, and environmental variables associated with field sites (n = 17).

Table S3. Summary of correlation and mixed model analyses of independent variables as related to Principal Component Analysis (PCA) axis 1 scores and Discriminant Analysis (DA) percentage of glyphosate‐resistant (GR) C. canadensis plants in field margins. Significance is indicated by *P < 0.05, **P < 0.01 and ***P < 0.001.

Table S4. General linear mixed model analyses, testing measures of reproductive output by population in the agrestal and ruderal habitats. Num, Den df, numerator degrees of freedom and denominator degrees of freedom, respectively.

Figure S1. Correlation analysis of principle component analysis (PCA) axis 1 scores of GR C. canadensis in 2008 field margins by field with latitude, where the most negative scores are correlated with lower latitudes, indicating greater and more variable occurrence of resistance at lower latitudes.

Figure S2. Percentage of individuals reaching the flowering stage by population and habitat, based on survivorship curve data for flowering stage.

Figure S3. Maximum estimate of herbivory (percentage of leaves damaged/individual) for individuals in the ruderal habitat. Mean values with the same letter are not significantly different (α = 0.05).

Click here for additional data file.