A Glyphosate-Based Herbicide Cross-Selects for Antibiotic Resistance Genes in Bacterioplankton Communities.

Agrochemicals often contaminate freshwater bodies, affecting microbial communities that underlie aquatic food webs. For example, the herbicide glyphosate has the potential to indirectly select for antibiotic-resistant bacteria. Such cross-selection could occur if the same genes (encoding efflux pumps, for example) confer resistance to both glyphosate and antibiotics. To test for cross-resistance in natural aquatic bacterial communities, we added a glyphosate-based herbicide (GBH) to 1,000-liter mesocosms filled with water from a pristine lake. Over 57 days, we tracked changes in bacterial communities with shotgun metagenomic sequencing and annotated metagenome-assembled genomes (MAGs) for the presence of known antibiotic resistance genes (ARGs), plasmids, and resistance mutations in the enzyme targeted by glyphosate (enolpyruvyl-shikimate-3-phosphate synthase; EPSPS). We found that high doses of GBH significantly increased ARG frequency and selected for multidrug efflux pumps in particular. The relative abundance of MAGs after a high dose of GBH was predictable based on the number of ARGs in their genomes (17% of variation explained) and, to a lesser extent, by resistance mutations in EPSPS. Together, these results indicate that GBHs can cross-select for antibiotic resistance in natural freshwater bacteria. IMPORTANCE Glyphosate-based herbicides (GBHs) such as Roundup formulations may have the unintended consequence of selecting for antibiotic resistance genes (ARGs), as demonstrated in previous experiments. However, the effects of GBHs on ARGs remain unknown in natural aquatic communities, which are often contaminated with pesticides from agricultural runoff. Moreover, the resistance provided by ARGs compared to canonical mutations in the glyphosate target enzyme, EPSPS, remains unclear. Here, we performed a freshwater mesocosm experiment showing that a GBH strongly selects for ARGs, particularly multidrug efflux pumps. These selective effects were evident after just a few days, and the ability of bacteria to survive and thrive after GBH stress was predictable by the number of ARGs in their genomes and, to a lesser extent, by mutations in EPSPS. Intensive GBH application may therefore have the unintended consequence of selecting for ARGs in natural freshwater communities.

INTRODUCTION

Glyphosate-based herbicides (GBHs) are by far the most extensively used weed-killers worldwide, especially since the introduction of transgenic glyphosate-resistant crops in the 1990s (1, 2). Glyphosate residues can spread widely and accumulate in soil, water, and plant products, raising concerns over human and environmental health (3). A recent systematic review and risk analysis concluded that glyphosate poses a moderate to high risk to freshwater biodiversity in 20 of the countries investigated (4). Some of the highest aquatic concentrations of glyphosate were found in countries with the largest production of genetically engineered glyphosate-tolerant crops globally, including the United States, Brazil, and Argentina (2, 4).

Although designed to control weed growth, glyphosate may also affect microorganisms that use the herbicide’s molecular target, the enzyme enolpyruvyl-shikimate-3-phosphate synthase (EPSPS), to synthesize aromatic amino acids (5). The EPSPS is classified into four classes according to mutations in the enzyme active site that confer differential sensitivities to glyphosate (6). In bacteria, EPSPS classes I and II, which are respectively sensitive and tolerant to glyphosate are the most frequently found, while classes III and IV are rarer and both confer glyphosate resistance (6). The EPSPS class II sequence isolated from a strain of Agrobacterium tumefaciens is used as the transgene conferring tolerance in most commercially available glyphosate-resistant crops (7, 8).

Experiments conducted in diverse environments, such as soil and freshwater (9–11) and the bee gut microbiome (12), have shown that bacterial taxa from natural ecosystems vary in their sensitivity to glyphosate. Some of this variation is explained by the distribution of different EPSPS classes. However, while strains with EPSPS class I are known to be sensitive, they have also been observed to tolerate glyphosate through unknown mechanisms (12), indicating that additional EPSPS-independent glyphosate resistance mechanisms exist in nature.

Studies with bacterial cultures have shown increased resistance to antibiotics after exposure to high concentrations of glyphosate and other herbicides (13–17). In the presence of glyphosate, the expression of membrane transporters may confer resistance to glyphosate and antibiotics simultaneously (18). Specifically, multidrug efflux pumps have been experimentally shown to confer resistance to both glyphosate and antibiotics, presumably by exporting a variety of small molecules (13, 14). This is an example of cross-resistance, a mechanism of indirect selection through which one resistance gene or biochemical system confers resistance to other antimicrobial agents (19, 20).

Direct selection of antibiotic resistance occurs when bacteria are exposed to an antibiotic agent and mutations conferring resistance to this agent are selected (21). In contrast, indirect selection for antibiotic resistance occurs in the absence of the antibiotic, either via cross- or co-resistance (19, 20). Cross-resistance occurs when the same gene confers resistance to multiple antibiotic agents, while co-resistance occurs when a resistance gene is genetically linked to another gene that is not necessarily an antibiotic resistance gene (ARG) but that is under positive selection.

Most studies of cross-resistance induced by herbicides have focused on bacterial isolates in laboratory experiments (13–16, 22). A recent study showed that herbicide selection increases the prevalence of ARGs in soil bacterial communities, using observational and experimental field data (23). However, we still lack evidence on aquatic communities, which are of particular interest because herbicides often reach water bodies through leaching, runoff, and spray drift from agricultural fields (4, 24). Moreover, the extent of direct selection on EPSPS mutations compared to indirect selection on ARGs is unclear. In a previous study, we used 16S ribosomal gene amplicon sequencing to assess how the composition of freshwater bacterioplankton communities responds to a GBH applied alone or in combination with a widely used neonicotinoid insecticide (11). As part of the same experiment, we also showed that phytoplankton and zooplankton communities responded strongly to high doses of GBH (25, 26). Here, we expand on our previous work and investigate the effects of the GBH on ARG frequencies in aquatic bacterial communities, using the same outdoor array of experimental ponds (Fig. 1A).

FIG 1

Experimental area and design. (A) Aerial photograph of the Large Experimental Array of Ponds (LEAP) at Gault Nature Reserve in Mont Saint-Hilaire (Canada). The laboratory facility and inflow reservoir, where water from our source lake was redirected to before filling the mesocosms, can be seen at the top of the photograph. Our source lake, Lake Hertel, is located upstream (not shown in the photograph). (B) Schematic representation of the subset of mesocosms selected for metagenomic sequencing in this study. A total of eight ponds were sampled 11 times over the course of the 8-week experiment, which was divided in two phases: phase I (6 weeks) and phase II (2 weeks). Phase I included two pulse applications (doses) of GBH, with three target glyphosate concentrations (0, 0.5, and 15 mg/liter). In phase II, all ponds except for two controls, shown in gray, received a higher dose of glyphosate (40 mg/liter). Phase I included four control ponds (gray and yellow), while phase II only included two controls (gray). Note that yellow ponds only received GBH in phase II. Nutrients were also added to ponds to reproduce mesotrophic or eutrophic conditions, represented by circles and squares (target phosphorus concentrations are indicated). TP, total phosphorus.

To test the extent to which GBH (in the form of Roundup) cross-selects for ARGs in complex aquatic communities over time, we exposed freshwater mesocosms to two glyphosate concentrations for 6 weeks (0.3 and 15 mg/liter; phase I) and to a higher dose for the next 2 weeks (40 mg/liter; phase II) (Fig. 1B). We sequenced metagenomes from each mesocosm and reconstructed metagenome-assembled genomes (MAGs) of bacteria, which were annotated according to their taxonomy, presence of ARGs, plasmids, and resistance mutations in the EPSPS enzyme (see Materials and Methods). We hypothesize that the frequency of ARGs would increase after exposure to a high concentration of glyphosate and that efflux pumps are among the main resistance mechanisms promoted by GBH. We also expect that MAGs encoding many ARGs and/or the resistant classes of the EPSPS gene will be the most likely to survive and proliferate after GBH exposure. Consistent with these expectations, we find that high doses of GBH (15 and 40 mg/liter glyphosate) cross-select for ARGs, particularly multidrug efflux pumps. These results demonstrate that severe contamination of aquatic systems with GBH could indirectly select for antibiotic resistance.

RESULTS

GBH treatment changes community composition and increases ARG frequency.

To assess how GBH treatments affected bacterioplankton community composition over time, we built principal response curves (PRCs). The response variables used for PRCs were either the estimated relative abundance of MAGs or the summed abundance of MAGs grouped at more inclusive taxonomic levels (phylum or class). In phase I of the experiment, two pulses of a GBH were applied to reach concentrations of 0.3 mg/liter or 15 mg/liter glyphosate. The first GBH pulse at the highest concentration changed the MAG composition irreversibly for the duration of the experiment (see Fig. S1A in the supplemental material). However, when these MAGs were grouped at more inclusive taxonomic levels, we observed a recovery of community composition 20 to 30 days after the first pulse (Fig. S1B and C). The second pulse in phase I had a weaker effect that depended on the taxonomic resolution, with the most pronounced changes visible at the class level. In phase II, when a single dose of 40 mg/liter glyphosate was applied to all mesocosms except for the phase II controls, a strong effect was observed in all ponds regardless of taxonomic resolution. These results are broadly consistent with our previous 16S rRNA gene amplicon sequencing from the same experiment, which showed community resilience to GBH pulses at higher taxonomic levels only (11).

Principal response curves (PRCs) of the experimental treatment effect on the composition of MAGs (A) and MAGs grouped at the phylum (B) and class (C) level. Treatment effect is shown on the left y axis, while taxon scores (proportional to their contribution to the treatment effect) are shown on the right y axis. Dashed vertical lines indicate the application of phase I glyphosate pulses and the solid vertical line the phase II glyphosate pulse. Glyphosate concentration of pulses applied in phase I (dose 1 and dose 2) are indicated by the legend, while in phase II all treatments received 40 mg/liter of glyphosate, except the phase II controls. Treatment effect zero is equivalent to the low-nutrient control phase II pond. (A) Only taxon scores higher than 0.1 are shown, the axis explains 16.1% of total variance, F = 14.6 and P = 0.011 by PERMUTEST; (B) 51.5% of total variance explained, F = 44.4 and P = 0.024 by PERMUTEST; (C) 31.8% of total variance explained, F = 34.3 and P = 0.008 by PERMUTEST. Download FIG S1, PDF file, 0.07 MB.

To test the effects of GBH on ARG frequency over the experiment, we tracked variation in the number of metagenomic reads mapped to the Comprehensive Antibiotic Resistance Database (CARD), here referred to as ARG reads, and in the counts of unique ARGs over time, both normalized by the total number of reads in each sample (Fig. 2). Among the GBH treatments applied in phase I, only the highest concentration increased ARG frequencies over time, either when measured as the number of unique ARGs (GAM F = 15.65, P < 0.001; Table 1 and Fig. S2) or as the number of ARG reads (GAM F = 15.78, P < 0.001; Table 1 and Fig. S2). The concordance of these two metrics suggests that the effect of GBH on ARGs was not due to a few highly responsive resistance genes but to multiple unique genes. The first GBH pulse produced a more prominent effect than the second pulse in phase I. After each pulse, the ARG frequencies returned close to their baseline, tracking with the community recovery observed at phylum and class taxonomic levels but not at the level of MAGs (Fig. S1). In phase II, the single dose of 40 mg/liter glyphosate triggered an increase in ARG frequencies across all treated ponds (Fig. 2). ARG frequencies increased over time due mainly to the phase II GBH pulse (Table 1 and Fig. S2). Nutrient enrichment produced a weak but significant effect only when considered alone, not in interaction with time (Table 1). Overall, these results support the hypothesis that the GBH treatment has the most dominant and strongest positive effect on ARG frequencies over time. The observation that ARG frequencies recover to baseline concordantly with higher taxonomic units but not the finest units (MAGs) suggests that ARGs are carried by a wide range of distantly related bacteria rather than specific species or strains.

FIG 2

ARG frequencies increase in GBH treatments over time. (A and B) Number of unique ARGs per million metagenomic reads (A) and number of metagenomic reads mapped to ARGs per million metagenomic reads (B) vary according to treatment and time. Dashed vertical lines indicate the application of phase I GBH pulses and solid vertical line the phase II pulse. The color code refers to the target glyphosate concentrations in phase I (pulse 1 and pulse 2), while in phase II all treated ponds received a target of 40 mg/liter glyphosate.

TABLE 1

Summary of GAMs showing the effect of GBH on ARG frequencies in phase I only and in both phases

Response variable/adjusted R2 and predictor	Factora	Estimate (SE) or EDF		t value or F		P valueb
		Phase I	Both phases	Phase I	Both phases	Phase I	Both phases
Unique ARG counts per million metagenomic reads [log10(x + 1)], adjusted R2 = 65.1% (phase I, n = 64)/74.6% (both phases, n = 88)
Parametric terms
Treatment	Control phase I	−0.001 (±0.006)	0.018 (±0.006)	−0.1	2.9	0.890	0.005*
	Glyphosate 0.3 mg/liter	0.003 (±0.006)	0.023 (±0.006)	0.5	3.8	0.633	<0.001*
	Glyphosate 15 mg/liter	0.035 (±0.006)	0.040 (±0.006)	6.2	6.7	<0.001*	<0.001*
Nutrient	High nutrient	-0.012 (±0.004)	-0.011 (±0.004)	−3.0	-2.7	0.005*	0.009*
Smooth terms
ti(day)		1.0	6.8	0.02	2.71	0.903	0.010*
ti(day, by=treatment)	Control phase I	1.0	3.2	0.03	10.04	0.861	<0.001*
	Glyphosate, 0.3 mg/liter	1.0	4.7	0.57	7.47	0.453	<0.001*
	Glyphosate, 15 mg/liter	3.9	4.3	15.65	5.01	<0.001*	0.002*
ti(day, by=nutrient)	High nutrient	1.0	1.0	0.60	0.25	0.444	0.620
ARG reads per million metagenomic reads [log10(x + 1)], adjusted R2 = 66.6% (phase I, n = 64)/77.3% (both phases, n = 88)
Parametric terms
Treatment	Control phase I	−0.028 (±0.105)	0.322 (±0.103)	−0.3	3.1	0.790	0.003*
	Glyphosate, 0.3 mg/liter	0.013 (±0.105)	0.320 (±0.103)	0.1	3.1	0.899	0.003*
	Glyphosate, 15 mg/liter	0.678 (±0.105)	0.804 (±0.103)	6.5	7.8	<0.001*	<0.001*
Nutrient	High nutrient	−0.197 (±0.074)	−0.185 (±0.073)	−2.7	−2.6	0.010*	0.013
Smooth terms
ti(day)		1.0	6.7	0.21	2.77	0.648	0.009*
ti(day, by=treatment)	Control phase I	1.0	3.5	0.11	12.27	0.737	<0.001*
	Glyphosate, 0.3 mg/liter	1.0	2.6	0.47	8.93	0.497	<0.001*
	Glyphosate, 15 mg/liter	3.9	4.6	15.78	6.90	<0.001	<0.001*
ti(day, by=nutrient)	High nutrient	1.0	1.0	3.92	1.52	0.053	0.222

When factor is absent, the respective predictor variable is continuous (“day”).

Asterisks indicate significant P values after Bonferroni correction (<0.0125).

For each predictor of the model, when it is a parametric term we report the respective parameter estimate with standard errors (SE) and t value; when it is a smooth term, we report the effective degrees of freedom (EDF) and F statistic. Smooths terms are described as the mgcv syntax, and “ti()” phrases are tensor product interactions. P values are reported for each predictor, and reports of significant factors after Bonferroni correction (P < 0.0125) are highlighted in boldface. A Gaussian residual distribution was used.

Glyphosate increases ARG frequencies in experimental ponds. GAMs illustrating the time-dependent effect of GBH and nutrient treatments on unique ARGs in phase I (A), in both phase I and II (B), on ARG reads in phase I (C), and in both phase I and phase II (D). Dashed vertical lines indicate the application of phase I GBH pulses and solid vertical line the phase II pulse. Glyphosate acid concentration of pulses applied in phase I (dose 1 and dose 2) are indicated in the legend, while in phase II all treatments received 40 mg/liter except the control phase II. Shades indicate a confidence interval of 95%. Download FIG S2, PDF file, 0.04 MB.

GBH selects for specific gene functions, including antibiotic efflux.

To assess how GBH affected known gene functions beyond ARGs in the bacterial communities, we built PRCs based on SEED annotations of genes in the metagenomes. The PRCs revealed a clear effect of GBH on the composition of gene functions (Fig. S3), similar to the treatment effect detected in the PRCs for community composition at the phylum and class levels (Fig. S1B and C). In phase I, the first pulse of 15 mg/liter glyphosate induced greater deviations from controls than the second pulse. In phase II, all ponds receiving 40 mg/liter glyphosate deviated from the controls. Resistance to antibiotics is among the functions positively affected by GBH treatment, as indicated by the positive scores of the SEED subsystems “virulence, disease, and defense,” at level 1 (Fig. S3A) and “resistance to antibiotics and toxic compounds” at level 2 (Fig. S3B). Table S1 shows the complete list of PRC scores for all SEED subsystems at levels 1 and 2. Membrane transporters (level 1, Fig. S3A), such as the ATP-binding cassette (ABC) transporters (level 2, Fig. S3B), are among the positively selected functions. These genes could plausibly change cell permeability to various molecules, including glyphosate.

PRC scores from functional annotations shown in Fig. S3. Download Table S1, XLSX file, 0.01 MB.

PRCs of the experimental treatment effect on the composition of gene functional profiles predicted from metagenomic reads grouped according to SEED subsystem level 1 (A) and level 2 (B). Treatment effect is shown on the left y axis, while scores of genes (proportional to their contribution to the treatment effect) are shown on the right y axis. Dashed vertical lines indicate the application of phase I glyphosate pulses and solid vertical line the phase II glyphosate pulse. Glyphosate concentration of pulses applied in phase I (dose 1 and dose 2) is indicated by the legend, while in phase II all treatments received 40 mg/liter glyphosate except the phase II controls. Treatment effect zero is equivalent to the low-nutrient control phase II pond. Function of resistance to antibiotics is highlighted in red according to how it is named in SEED subsystem level 1 (50.9% of total variance explained, F = 43.1 and P = 0.023 by PERMUTEST) (A) and SEED subsystem level 2 (33.1% of total variance explained, F = 25.8 and P = 0.027 by PERMUTEST) (B), where only scores with absolute values larger than 0.05 are reported (all scores are shown in Table S1). Download FIG S3, PDF file, 0.06 MB.

To assess the effects of GBH on ARGs at a higher level of resolution, we built another set of PRCs based on ARG profiles predicted from reads mapping to CARD. The resulting PRC plot showed a prominent effect of the first and second pulses of 15 mg/liter glyphosate in phase I (Fig. 3). In phase II, the GBH had an effect in all treatments that received a last pulse (40 mg/liter glyphosate). This result is consistent with the greater effect of the large phase II pulse compared to smaller phase I pulses on total ARG frequencies (Fig. 2 and Fig. S2). The two principal resistance mechanisms of the ARGs annotated by CARD are antibiotic efflux and antibiotic inactivation (shown in blue and red text in Fig. 3). Genes encoding antibiotic efflux functions were more often found with positive PRC scores (Fisher’s exact test, P = 0.013), suggesting that they tend to be selected more often than other ARGs in the presence of GBH. This result supports the hypothesis that membrane transporters used for antibiotic efflux also play a role in exporting glyphosate from bacterial cells.

FIG 3

GBH skews composition of ARGs in favor of antibiotic efflux pumps. Principal response curves (PRCs) illustrating divergence (relative to controls) in the composition of ARGs in response to GBH exposure. The left y axis represents the magnitude or ARG compositional response, while the right y axis represents individual gene scores (i.e., relative contribution to overall compositional changes). Gene names (ARO) are color-coded based on their mechanism of resistance. Dashed vertical lines indicate the timing of GBH pulses in phase I, and the solid vertical line represents the pulse in phase II. The zero line (y = 0) represents the low-nutrient control pond from both phases I and II. The PRC explains 30% of the total variance (F = 22.8, P = 0.024 by PERMUTEST permutation test for redundancy analysis). Treatments and time interactively explain 74.8% of the variance, while 25% is explained by time alone.

Connecting resistance genes to genomes and plasmids.

Thus far, our results have only considered ARGs outside the context of the bacterial genomes or plasmids in which they occur. On average, 71% (±3%; range, 45 to 94%; Table S2) of ARG reads (those mapping to CARD) across all samples also mapped to MAGs, implying that MAGs captured a large fraction of ARG reads in the metagenomes. We identified putative plasmids in 390 MAGs, with an average of 43 plasmid contigs per MAG (minimum, 1; maximum, 520; standard errors [SE],  3.5; Table S3). However, only 27 plasmid contigs in total were annotated with ARGs. Out of a total of 188 MAGs with ARGs, only 24 (13%) of them had at least one ARG identified in a potential plasmid. Although some ARGs are clearly found on plasmids, ARGs are more associated with genomes than with MAG plasmids in our study.

Metagenomic sample information, summary of RGI output for hits above mapping threshold (MAPQ of >10 and minimum of 50% gene coverage) and proportion of sample reads mapped to CARD (ARG reads) that mapped back to MAGs. Download Table S2, XLSX file, 0.01 MB.

MAG information, predicted EPSPS amino acid sequence, summary of ARGs, and plasmids. For each predicted EPSPS sequence, the putative classification regarding glyphosate resistance is shown. The number of potential plasmid contigs and how many of these had ARGs annotated is also shown. The number of ARGs annotated to MAG contigs (total RGI strict hits) is provided in the last column. Download Table S3, XLSX file, 0.1 MB.

Of the 426 total MAGs, only 20 recruited 100 or more ARG reads, and the classification of their EPSPS genes varied (Fig. S4 and S5). To visualize which ARGs were more abundant in GBH treatments and in which MAGs they were found, we examined the frequency of metagenomic reads mapped to ARGs according to their antibiotic resistance ontology (ARO) classification (top graphs in Fig. S4 and S5) as well as the proportion of these reads that were mapped to MAGs (bottom graphs in Fig. S3 and S5). These visualizations confirmed that efflux pumps (e.g., mex genes) increased in frequency in response to GBH. The relative abundance of mex genes is strongly associated with a Pseudomonas putida MAG (Fig. S4, bottom right) but is sometimes also associated with other MAGs, such as Aeromonas veronii (Fig. S4), Oxalobacteraceae, and Azospirillum (Fig. S5). Thus, it is likely that GBH selects for efflux pump genes in multiple different genomic backgrounds.

Metagenomic reads mapped to ARGs classified according to their ARO (top graph) and ARG reads mapped to MAGs (bottom graph) in low-nutrient ponds. MAG identities are followed by their finest taxonomic assignment (o, order; f, family; g, genus; s, species). Only alignments with MAPQ of >10 were tallied. Dashed vertical lines represent phase I GBH, and solid vertical lines are phase II pulses (all at 40 mg/liter glyphosate). Download FIG S4, PDF file, 0.05 MB.

Metagenomic reads mapped to ARGs classified according to their ARO (top graph) and ARG reads mapped to MAGs (bottom graph) in high-nutrient ponds. MAG identities are followed by their finest taxonomic assignment. Only alignments with MAPQ of >10 were tallied. Dashed vertical lines represent phase I GBH pulses, and solid vertical lines are phase II pulses (all at 40 mg/liter glyphosate). Download FIG S5, PDF file, 0.05 MB.

The number of ARGs in a MAG predicts its frequency after severe GBH exposure.

Collectively, our results suggest an important role for ARGs, and efflux pumps in particular, in allowing bacterioplankton to survive and grow in the presence of a GBH. We next asked about the importance of ARGs relative to genetic variation in the glyphosate target enzyme EPSPS. Based on known sequence variation in the EPSPS gene, we were able to classify MAGs as putatively glyphosate resistant, sensitive, or unclassified. We also defined a MAG's antibiotic resistance potential as the number of ARGs they contain (i.e., the number of resistance gene identifier [RGI] strict hits). These definitions are genetic predictions based on the presence or absence of resistance-associated genes in MAGs, but they do not represent confirmed phenotypes. We nonetheless tested the extent to which these genomic features were predictive of a MAG's average relative abundance across ponds at the end of the experiment, after receiving 40 mg/liter glyphosate in phase II. We found that MAGs encoding more unique ARGs tended to have higher relative abundance after receiving the GBH pulse in phase II (Fig. 4A, Table 2). This effect of antibiotic resistance potential was highly significant (multiple linear regression model, t = 9.53, P < 0.001; Table 2) and was not observed in control ponds that did not receive the phase II pulse (Table S1 and Fig. S6A; t = 2.26, P = 0.025; not significant after Bonferroni correction). The relative abundance of MAGs at the end of the experiment in these control ponds was predicted by their relative abundance in phase I (40% of variance explained; Table S4), consistent with temporal autocorrelation (e.g., due to fluctuations in species abundances unrelated to experimental treatments). In contrast to the strong effect of ARGs on predicting MAG relative abundance after GBH stress (17% of variance explained; Table 2), EPSPS classification explained only 2% of the variation in both phase II treatment and control ponds.

FIG 4

Antibiotic resistance potential predicts MAG relative abundance after severe GBH stress. (A) Boxplots show a positive correlation between MAG abundance in phase II and their potential for antibiotic resistance. Each dot represents a MAG that is color-coded based on the predicted resistance of their EPSPS. A slight offset on the x axis (jitter) was introduced to facilitate data visualization. See Table 2 for regression coefficients. (B) Regression tree confirms the significance of the correlation seen in panel A, particularly for antibiotic efflux genes. Two other factors were also included and have small effects on MAG relative abundance in phase II: the EPSPS classification and the average abundance of MAGs in phase I.

TABLE 2

Multiple linear regression model and variance partitioning of MAGs abundance in phase II in treatment mesocosms

Predictor	Estimate (SE)	t value	P value	Explained variance (%)
EPSPS classification				2
Sensitive	0.002 (±0.127)	0.02	0.987
Resistant	0.413 (±0.133)	3.11	0.002*
MAG antibiotic resistance potential	0.496 (±0.052)	9.53	<0.001*	17
MAG mean abundance in phase I treatment mesocosms (log10)	0.178 (±0.066)	2.69	0.007	1; residuals, 79

P values are reported for each predictor, asterisks indicate significant P values after Bonferroni correction (P < 0.0125), and reports of significant factors are highlighted in boldface. Adjusted R-squared equals 21.1% for MAG persistence in treatments (n = 426; F = 29.5). The response variable was MAG average abundance (log10) in phase II treatment mesocosms.

Multiple linear regression model and variance partitioning of MAGs abundance in phase II in control mesocosms. P values are reported for each predictor. Asterisks indicate significant P values after Bonferroni correction (P < 0.0125), and reports of significant factors are highlighted in boldface. Adjusted R2 equals 43.2% for MAG abundance in controls as response variable (n = 425, F = 78.7). Download Table S4, XLSX file, 0.01 MB.

MAG mean relative abundance in controls of phase II as a function of antibiotic resistance potential (or the amount of ARGs annotated to their genomes) and the classification of EPSPS enzyme (resistant, sensitive, or unclassified). (A) Series of boxplots show the absence of correlation between MAGs abundance in phase II and their potential for antibiotic resistance. Each dot represents a MAG that is color-coded according to the potential resistance of their EPSPS. To facilitate visualization, a small amount of random variation (jitter) was added so dots would not overlap. Table 2 reports statistics of a linear model that tested how MAG abundance in phase II controls could be explained by EPSPS classification, antibiotic resistance potential, and MAG abundance in phase I. (B) Regression tree with MAG abundance in controls of phase II as the response variable and the following predictors: the EPSPS enzyme classification, the number of ARGs classified as antibiotic efflux, antibiotic inactivation or target alteration, and the MAG relative abundance in phase I. Download FIG S6, PDF file, 0.07 MB.

To further explore these results, we used a regression tree analysis to identify the drivers of MAG abundance at the end of phase II. Instead of combining the three major classes of ARGs (antibiotic target alteration, antibiotic inactivation, and antibiotic efflux), we used each as a separate predictor in the regression tree. The first division splits MAGs with at least one antibiotic efflux gene (Fig. 4B, node 7), which were, on average, more abundant post-GBH pulse than those without efflux genes (Fig. 4B, node 2). Among MAGs with efflux genes, the more genes they had, the higher their abundance. Among MAGs without antibiotic efflux genes, the EPSPS classification was an important driver of their abundance, followed by the MAG's average abundance in phase I. In the absence of a GBH pulse in phase II, the primary driver of MAG abundance in phase II controls was their mean relative abundance in phase I (Fig. S6). Control pond regression trees also included a split between resistant/sensitive and unclassified EPSPS, which is difficult to interpret biologically and likely attributable to noise. This could also explain why 2% of the variation in MAG relative abundance in control ponds was explained by EPSPS class. Together, these results indicate that a bacterial genome's ARG coding potential is predictive of its ability to persist in the face of GBH stress, more so than the class of EPSPS enzyme it encodes.

DISCUSSION

Our mesocosm experiment used deep metagenomic sequencing to assess the effect of a GBH on microbial genes and genomes in seminatural freshwater bacterial communities. We show that exposure to Roundup at high concentrations (15 mg/liter and 40 mg/liter glyphosate) changes community composition and increases the frequency of ARGs in freshwater bacterioplankton. Consistent with our previous 16S rRNA gene amplicon sequencing from the same experiment (11), we found that more inclusive taxonomic groupings were resilient to GBH pulses, whereas the precise MAG composition never recovered. Moreover, we show that the abundance of MAGs after severe contamination (40 mg/liter glyphosate) was predictable based on the number of ARGs encoded, and such successful MAGs tended to have at least one antibiotic efflux gene annotated in their genomes. The effect of GBH on ARGs is likely due to cross-resistance, since the multidrug efflux pumps that rise in frequency in response to GBH could transport glyphosate in addition to antibiotics (18). Alternatively, co-resistance could play a role if GBH selects for glyphosate-resistant bacterial genomes or genetic elements (rather than specific genes) that happen also to encode ARGs. While we cannot exclude a role for co-resistance entirely, the cross-resistance model is more plausible since efflux genes are strongly affected, likely in multiple independent genomic backgrounds. As discussed in detail below, direct selection on the EPSPS locus appears to be weak, implying that ARGs are unlikely to achieve high frequency due to genetic linkage with resistant EPSPS alleles.

There are several explanations for the observation that ARG frequencies and higher taxonomic units, phyla and classes but not MAGs, concurrently increase and then recover to baseline following GBH pulses. One technical explanation is that ARG families are defined at levels of genetic similarity more in line with phyla or classes than with finer taxonomic levels. It is also possible that the recovery of ARG frequencies to baseline after phase I pulses could be due to gene loss and horizontal transfer events (as discussed further below). However, even if such events occur, they are not sufficient to obscure the predictability of successful MAGs after a severe GBH pulse, which requires relatively tight linkage of ARGs with genomes.

An association between glyphosate, ARGs, and mobile genetic elements was previously found in soil microbiomes, as demonstrated in a recent study combining experimental microcosms and environmental data from agricultural field sites in China (23). Through laboratory assays of three bacterial strains, the authors quantified the conjugation frequency of a multidrug resistance plasmid induced by glyphosate and further investigated changes in cell membrane permeability. They detected a significant increase in conjugation frequency and augmented cell membrane permeability in the presence of glyphosate, suggesting that glyphosate stress increases membrane permeability, thereby promoting plasmid movement. Here, we provide additional support for the hypothesis that cell membrane permeability is altered in the presence of a GBH, as demonstrated by the selection of membrane transport mechanisms, such as ABC transporters (27), among the annotated gene functions most responsive to our GBH treatments. In contrast, although we did not quantify the frequency of conjugation in our experiment, we did identify some ARGs located on putative plasmids that were present in a small fraction of the bacterial community. Of the MAGs encoding ARGs, only 13% contained a predicted plasmid-encoded ARG. It is possible that unassembled plasmids or plasmids not associated with MAGs could harbor ARGs. Including such plasmids would not be expected to change our main conclusion that ARGs are more predictive of MAG frequency post-GBH exposure than EPSPS. In addition to plasmids, other mechanisms also contribute to horizontal gene transfer between bacteria, such as phage-mediated transduction and transformation (28), and future studies could test how these processes can be affected by GBH stress.

Strikingly, antibiotic resistance potential, particularly the presence of antibiotic efflux genes, was more important than the EPSPS classification in explaining variation in MAG abundance in phase II after a high-GBH pulse. This evidence of cross-resistance in seminatural communities may help explain why, in previous experiments also performed with complex communities, bacterial strains with the sensitive EPSPS-encoding gene were resistant to glyphosate, as was the case in two strains of Snodgrassella alvi in the bee gut microbiome (12). Although EPSPS alleles were weakly predictive of MAG relative abundance after the phase II GBH pulse, their effects were clearly secondary to the strong effects of ARGs. Computational gene annotations of both ARGs and resistant or sensitive EPSPS have limitations because they are based on sequence similarity, not on phenotypic measurements. Therefore, we cannot exclude a role for EPSPS alleles in conferring GBH resistance in nature, but their effects were small in our experiment. Together, our results strongly suggest that ARGs (and efflux pumps in particular) are more relevant to glyphosate resistance in nature than mutations in the glyphosate target enzyme.

Our study aligns with previous single-strain laboratory evidence that antibiotic resistance may enhance bacterial survival in the presence of pesticides. Laboratory assays of bacterial isolates showed that exposure to agrochemicals accelerated the rate of antibiotic resistance evolution (15, 16). In other studies, depending on the combination of herbicide and antibiotics tested, herbicides increased or decreased antibiotic susceptibility of bacterial strains (13, 14). In our experiment, GBH pulses caused a general increase in ARG frequency in the bacterioplankton community, suggesting an overall positive effect of GBH on resistance to antibiotics. This does not exclude the possibility that GBH could also select against specific ARGs or mutations, making bacteria less resistant to certain antibiotics. As a metagenomic study, we only measured genetic correlates of resistance, and future experiments will be needed to examine resistance phenotypes.

Our study shows that antibiotic efflux is a major mechanism of antibiotic resistance that is cross-selected by GBH stress, which also corroborates previous laboratory assays. It has been shown that the targeted deletion of efflux pump genes can neutralize the increased tolerance to kanamycin and ciprofloxacin in Escherichia coli and Salmonella enterica serovar Typhimurium in the presence of GBH (13, 14). Here, we provide evidence that efflux pumps also provide resistance to both glyphosate and certain antibiotics in a more natural and complex system. Whether all efflux pumps are equally capable of transporting various molecules out of the cell remains to be seen, and other resistance mechanisms could also play a role. Resistance to both antibiotics and GBH could also be modulated by changing the expression of efflux pumps. While we used a metagenomic approach to show how GBH affects ARG frequencies in a community, this does not exclude the possibility of changes in gene expression, which could be tracked using metatranscriptomics in future experiments.

It should be noted that we used a commercial Roundup formulation of the herbicide glyphosate, which includes other constituents that may also influence microbial communities and cellular physiology. For example, the surfactant polyethoxylenamine (POEA) has produced negative effects on Vibrio fischeri at lower concentrations than pure glyphosate acid (29). However, given that our results are in general agreement with previous soil experiments using pure glyphosate (23), we believe that our findings are at least in part attributable to an effect of glyphosate itself. Furthermore, regardless of whether it is glyphosate or other constituents of GBH that drive cross-selection of ARGs, assessing the risks associated with commercial formulations is ecologically relevant, as these formulations are used in agriculture fields and lawns (30).

On an applied level, the safety assessment process for pesticides such as glyphosate, currently based on toxicity to model eukaryotic organisms (31, 32), could also consider the potential effects on bacterioplankton and ARGs. Our results highlight the role of GBH contamination as an indirect selective pressure favoring ARGs in natural communities. Although glyphosate concentrations as high as the ones inducing this effect (i.e., 15 mg/liter and 40 mg/liter) are rarely found in nature, there are reports of even higher glyphosate levels during the rainy season close to agricultural fields, as observed in Argentina (105 mg/liter), for example (4). Additionally, currently regulated acceptable concentrations of glyphosate in freshwaters in the United States and Canada for short-term exposure (1 to 4 days) are close to the concentrations used in our experiment (49.9 mg/liter [32] and 27 mg/liter [31]). Here, we have shown that ARG frequencies can rise dramatically just a few days after GBH treatment at such doses. The extent to which these ARGs, and the bacteria that carry them, can be mobilized across aquatic ecosystems, and from these ecosystems into animals and humans, remains to be seen.

MATERIALS AND METHODS

Experimental design.

An 8-week mesocosm experiment was conducted at the Large Experimental Array of Ponds (LEAP) facility (Fig. 1A), located at McGill University’s Gault Nature Reserve (QC, Canada), from 17 August (day 1) to 12 October (day 57) 2016, as previously described (11, 25, 26). Pond mesocosms were filled with 1,000 liters of water and planktonic communities from Lake Hertel (45°32′ N, 73°09′ W). Lake water was passed through a coarse sieve to prevent fish introduction while retaining lake bacterioplankton, zooplankton, and phytoplankton, whose responses to experimental treatments have been described in previous studies (11, 25, 26).

Figure 1B illustrates the experimental design of a subset of eight treatments selected for the metagenomic sequencing analyses reported here (see reference 25 for a full description of all treatments at the LEAP facility in 2016). The eight ponds were sampled at 11 time points throughout phases I and II of the experiment. In phase I (days 1 to 44), all ponds received nutrient inputs biweekly, simulating mesotrophic or eutrophic lake conditions with additions of a concentrated nutrient solution. Four ponds were treated with a GBH to reach target concentrations of 0.3 or 15 mg/liter of the active ingredient (glyphosate; acid equivalent), while the other four were kept as control ponds. The GBH was applied in two pulses in phase I, at days 6 and 33. In phase II (days 45 to 57), two control ponds (here referred to as control phase I) and the four treatment ponds received one pulse of the GBH at a higher dose (40 mg/liter glyphosate) on day 44, while the other two control ponds (here referred to as control phase II) received no pulse.

Target doses of the active ingredient were calculated based on the glyphosate acid content in Roundup super concentrate grass and weed control (reg. no. 22759; Bayer), the formulation used for the experiment. We used a commercial formulation to mimic environmental contamination and because the costs of using pure glyphosate salt would be prohibitive in a large-scale field experiment. Treatments are referred to by their glyphosate acid concentration to allow comparison with other formulations. Nutrients were added in the form of nitrate (KNO3) and phosphate (KH2PO4 and K2PO4), with target concentrations of 15 μg P/liter and 231 μg N/liter in the low-nutrient (mesotrophic) treatment ponds and 60 μg P/liter and 924 μg N/liter in the high-nutrient (eutrophic) treatment ponds. The concentrated nutrient solution had an N:P molar ratio of 33, comparable to our source lake. As reported in previous studies (11, 25), target doses of glyphosate acid and nutrients were achieved reasonably well, although glyphosate accumulated over phase I, reaching higher concentrations than intended after the second dose.

DNA extraction and metagenomic sequencing.

The eight experimental ponds were sampled for bacterioplankton DNA at 8 time points during phase I (days 1, 7, 15, 30, 35, 38, 41, and 43) and 3 time points during phase II (days 45, 49, and 57). Water samples were collected with 35-cm-long integrated samplers (2.5-cm-diameter PVC tubing) at multiple locations in the same pond and stored in 1-liter dark Nalgene bottles at 4°C until being filtered within 4 h. We filtered 250 mL of each sample onsite through 0.22-μm-pore-size Millipore hydrophilic polyethersulfone membranes 47 mm in diameter (Sigma-Aldrich, St. Louis, USA). Filters were stored at –80°C until DNA extraction.

We extracted DNA from a total of 88 filter samples using the PowerWater DNA isolation kit (MoBio Technologies Inc.) by following the manufacturer’s guidelines. Shotgun metagenomic sequencing was performed using the Illumina HiSeq 4000 technology with 100-bp paired-end reads. Libraries were prepared with 50 ng of DNA using the NEBNext Ultra II DNA library prep kit for Illumina (New England Biolabs) per the manufacturer’s recommendations and had an average fragment size of 390 bp.

Metagenomic read trimming, functional annotation, and ARG inference from metagenomic reads.

We removed Illumina adapters and quality filtered metagenomic reads using Trimmomatic (33) in the paired-end mode. We used FragGeneScan (34) for gene prediction from trimmed metagenomic reads and annotated predicted genes with SEED subsystems (35). To identify known ARGs in the metagenomic reads, we used the resistance gene identifier (RGI) “bwt” function that maps FASTQ files of reads passing quality control to CARD (36) using Bowtie2 (version 2.4) as an aligner (37). Only alignments with mapping quality (MAPQ) higher than 10 and gene coverage of 50% were retained. To calculate the proportion of metagenomic reads mapped to CARD that have been assembled and binned to genomes, we extracted reads that aligned to CARD using SAMtools (38) and mapped them to MAGs using Bowtie2 (37). Table S2 in the supplemental material shows the total number of reads by sample after trimming and a summary of the RGI output by sample for hits with minimum gene coverage of 50% and average MAPQ of >10.

Metagenomic de novo coassembly, binning, dereplication, and curation of MAGs.

We organized the data set into eight sets of metagenomes, each of them containing samples of the same mesocosm pond (Fig. 1B) from multiple time points. We coassembled reads from each of the 8 time series using MEGAHIT v1.1.1 (39), with a minimum contig length of 1 kbp. We used anvi’o v5.1 (40) to profile contigs, to identify genes using Prodigal v2.6.3 (41) and HMMER v3.2.1 (42), to infer the taxonomy of genes with Centrifuge v1.0.4 (43), to map metagenomic reads to contigs using Bowtie2 v2.4.2 (37), and then to estimate depth of read coverage across contigs. Finally, we used anvi’o to cluster contigs according to their sequence composition and coverage across samples with the automatic binning algorithm CONCOCT (44), and we manually refined the bins (n = 830) using the anvi’o interactive interface, as suggested by the developers (40), by removing splits that diverged in the differential coverage and/or tetranucleotide frequency of most splits in the same bin.

We dereplicated bins as described previously (45). In summary, we calculated the Pearson correlation coefficient between the relative abundance (i.e., the mean coverage calculated by the function “anvi-summarize” within anvi’o) for each pair of bins in the metagenomic samples, using the “cor” function in R (46) and the average nucleotide identity (ANI) of bins affiliated with the same phylum using NUCmer (47). Taxonomy assignment of redundant bins was done using CheckM (48). Bins with a Pearson correlation coefficient above 0.9 and ANI of 98% or more were considered redundant. In a total of 830 bins obtained before performing the dereplication, we found 607 nonredundant bins, of which 426 were classified as MAGs, as they had at least 70% completeness and no more than 10% redundancy (Table S2). We then created a nonredundant genomic database of these 426 MAGs to which we mapped metagenomic reads to calculate the relative abundance of each MAG across the different samples. Here, we define a MAG's relative abundance as the number of metagenomic reads recruited to a MAG divided by the total number of metagenomic reads in a given sample.

Identification of ARGs, EPSPS, and plasmids in MAGs.

We annotated ARGs within MAG contigs with the RGI “main” function, which compares predicted protein sequences from contigs to the CARD protein reference sequence data. Within RGI, we used the BLAST (49) alignment option and the strict algorithm (excluding nudge of loose hits to strict hits) for low-quality contigs (<20,000 bp). The RGI low sequence quality option uses Prodigal anonymous mode (41) for the prediction of open reading frames, supporting calls of partial ARGs from short or low-quality contigs.

To identify EPSPS sequences from MAG contigs, we first used anvi’o to predict amino acid sequences of the nonredundant MAGs with the flag “report-aa-seqs-for-gene-calls” of the function “anvi-summarize.” Gene calls of all the MAGs were concatenated, conserving the original split names, and transformed into a fasta file. We then searched the predicted amino acid sequences against a custom database with sequences of the EPSPS enzyme, using BLASTp (49) and a minimum E value of 1e−5. After selecting the gene call with the best match (i.e., lowest E value) to an EPSPS sequence in each of the 426 MAGs, we used the EPSPSClass web server (6) to classify the retrieved sequences according to resistance to glyphosate. Sequences were classified as EPSPS class I, class II, or class IV if they contained all the amino acid markers from the respective reference, i.e., if the percent identity was equal to 1, and classified as class III when they contained at least one complete motif out of 18 of the resistance-associated sequences, as explained in reference 6. MAGs whose EPSPS sequences did not match these criteria of having at least one motif of class III or 100% identity to class I, II, or IV, or those in which no predicted amino acid sequence matched a known EPSPS sequence, were set as unclassified (roughly 27% of MAGs). EPSPS sequences matching class I were considered putative sensitive and those with at least one motif of class III or matching class II as putative resistant. No sequences were found that matched class IV.

To identify potential plasmid contigs assembled to MAGs, we used the plasmid classifier PlasClass (50). We counted all contigs classified as plasmids with a minimum of 70% probability as well as how many of these potential plasmid contigs were annotated with ARGs through RGI. Table S2 summarizes MAG information, including the predicted EPSPS sequence found in the genome, the EPSPS classification, the number of estimated plasmid contigs, and how many of them contained ARG sequences.

Statistical analyses.

All statistical analyses were conducted in R v.4.0.2 (46). Time series of (log-transformed) ARG counts and ARG reads per million metagenomic reads were modeled using additive models (GAM) using the “mgcv” R package (51). We used GAMs to account for nonlinear relationships among the response variable and the predictors. Some predictors (nutrient and herbicide treatment levels) were coded as ordered factors; Table 1 lists all factors and predictors of the model. We built the models using the “gam” function and assessed significance of effects with the “summary.gam” function. We validated the models with the “gam.check” function, inspecting the distribution of model residuals, comparing fitted and observed values, and checking if the basis dimension (k) of smooth terms was large enough.

We used PRCs to test for the effect of treatments on the composition of MAGs, ARGs, and gene functional profiles over time. PRCs are a special case of partial redundancy analysis (pRDA) used in temporal experimental studies where treatments and the interaction between treatment and time are used as explanatory variables (52). Time is the covariable (or conditioning variable) whose effect is portioned out, and the response variable is the matrix containing compositional data (taxon or gene family relative abundances). We built PRCs using relative abundances of predicted genes grouped according to the SEED subsystem levels 1 and 2. In a more focused analysis, we built a PRC for the matrix of ARGs found in each sample, i.e., metagenomic reads mapped to each ARG from the CARD reference classified according to their antibiotic resistance ontology (ARO). The matrices were transformed using the Hellinger transformation (53). The PRC diagram displays the treatment effect on the y axis, expressed as deviations from the experimental controls at each time point. It also shows species scores on the right y axis, which here can be interpreted as the contribution of each function or gene to the treatment response curves. We assessed the significance of the first PRC axis by permuting the treatment label of ponds while keeping the temporal order, using the “permute” package (54), followed by a permutation test (999 permutations) using the “vegan” package (55). For the PRC based on ARG composition, we tested if the distribution of PRC positive and negative scores was different among the resistance mechanisms of the identified ARGs using the “fisher.test” function in the “stats” package in R (46).

To test if MAG abundance in phase II glyphosate treatments was correlated with their antibiotic resistance potential, we built a multiple linear regression with the “lm” function of the R package “stats” (46). The response variable was the average relative abundance of a MAG in glyphosate-treated ponds in phase II, averaging abundance across all ponds in which a MAG was found in phase II. The three predictors were the MAG's antibiotic resistance potential (defined as the number of RGI strict hits found in the MAG), the average MAG relative abundance in the same ponds of phase I, and their EPSPS sequence classification (resistant, sensitive, or unclassified). To assess the relative contribution of the different predictors to MAG survival in phase II, we performed a variance partitioning analysis with the “varpart” function of the R package “vegan” (55). Finally, to visualize the hierarchy among predictors, we constructed a conditional inference regression tree. Response variable and predictors were the same as those described above, except that instead of grouping all ARG hits, we transformed them into three variables according to their function: antibiotic target alteration, antibiotic inactivation, or antibiotic efflux. The regression tree was fitted with the “ctree” function in the R package “party” (46). As a negative control, we repeated the same analyses for MAGs found in control ponds of phase II.

As multiple predictors were tested, we used a conservative Bonferroni correction for the additive and linear models, whereby the P value significance threshold of 0.05 was divided by the number of statistical tests.

Graphs and heatmaps for time series data visualization were built using the functions “geom_point” and “geom_tile,” respectively, in the R package “ggplot2” (56).

Data accessibility.

Sequence data of the 88 metagenomic samples were submitted to NCBI SRA (BioProject no. PRJNA767443, accession numbers SRR16126824–SRR16126911), and the genomes of 426 predicted MAGs have been deposited and associated with the same BioProject (BioSample accession numbers are in Table S3).