High-throughput functional screening reveals low frequency of antibiotic resistance genes in DNA recovered from the Upper Mississippi River

Antibiotic-resistant bacteria were first reported in the 1940s and have become an increasing public health concern due to the diversity of resistant species, as well as an increasing number of antibiotics against which resistance has been observed (Levy & Marshall 2004; Capita & Alonso-Calleja 2013). Acquisition of antibiotic resistance is primarily attributed to the misuse and overuse of antibiotics in medicine and animal husbandry (Finch 2004; Bywater 2005). However, it is also well known that environmental bacteria produce antibiotics, and naturally possess intrinsic resistance mechanisms (Martínez 2008). While the role of antibiotics in inhibiting growth of microbial competitors is well understood, several hypotheses suggest that low concentrations of antibiotics may additionally serve as signaling molecules (Fajardo & Martínez 2008) or, in some cases, as nutrient sources (Dantas et al. 2008). Furthermore, genes that confer resistance to elevated concentrations of antibiotics may also have drastically different metabolic or ecological functions in the environment (Martinez et al. 2009).

Some, if not all, antibiotic resistance genes (ARGs) have evolved chromosomally to either confer resistance to an antimicrobial compound or to perform a separate function, which may confer antibiotic resistance in another organism (D'Costa et al. 2006). However, due to anthropogenic impacts, including the reintroduction of antibiotics into the environment, many of these genes have also been incorporated into mobile genetic platforms (e.g. plasmids and transposons) making their spread to pathogens more efficient via horizontal gene transfer (HGT) (Alonso et al. 2001; LaPara et al. 2011). The spread of ARGs among pathogenic bacteria as well as in wastewater communities has received considerable attention, but only in the last decade or so have studies focused on environmental reservoirs, especially surface waters, as important reservoirs for ARGs (Baquero et al. 2008; LaPara et al. 2011). Recently, bacterial communities in water environments including rivers, streams, and lakes have been implicated as important vehicles for the retention and transfer of ARGs to human pathogens (Baquero et al. 2008; Lupo et al. 2012).

Although river water samples have been shown to harbor lower concentrations of antibiotics than associated sediment samples (Kim & Carlson 2007), sub-inhibitory concentrations of antibiotics have been shown to support the selection of resistant phenotypes (Gullberg et al. 2011). Maintenance of antibiotic resistance generally confers a fitness cost to the organism. This suggests that a lack of selective pressure would select for reversal of resistance (i.e. sensitivity). However, compensatory evolution and/or mitigation of fitness cost through community-level interactions has been shown to drastically slow reversal among bacterial communities (Andersson & Hughes 2010). Owing to their persistence and potential for HGT between pathogens and non-pathogens, ARGs themselves, separate from antibiotics, have been suggested to be emerging contaminants in surface waters (Pruden et al. 2006; Martinez 2009).

Aquatic systems, in particular rivers, have recently been suggested as drivers for the spread of ARGs due to constant mixing of the bacterial community as well as a multitude of anthropogenic impacts (Baquero et al. 2008; Taylor et al. 2011). The spread of ARGs as a result of this mixing may result in the transfer of novel resistance genes from indigenous bacteria to potential human pathogens that, in turn, transfer the newly acquired resistance to the clinical setting (Poirel et al. 2002, 2005; Wright 2010). Furthermore, selective pressures from the contribution of antibiotics from agricultural, industrial, and municipal runoff, as well as wastewater effluent, may promote the emergence of new resistant phenotypes, or ARGs may be indirectly selected for in response to adaptations to other nutrient or chemical contamination from anthropogenic sources (Martinez 2009).

One of the major difficulties in evaluating the resistome of environmental communities – the suite of ARGs present in the community (D'Costa et al. 2006) – is the lack of culturability of >99% of environmental species (Amann et al. 1995). To circumvent this obstacle, recent studies assessing the prevalence and distribution of ARGs in environmental samples have relied on culture-independent quantitative polymerase chain reaction (qPCR) assays (Huerta et al. 2013; Marti et al. 2013). However, functional metagenomic screening of large clone libraries has been recently suggested to characterize antimicrobial resistance frequencies of environmental samples (Martínez & Osburne 2013). This approach has been taken using bacterial artificial chromosome libraries to assess antibiotic resistance in soil communities (Riesenfeld et al. 2004), and more recently, fosmid libraries were constructed and screened to assess levels of antibiotic resistance in river sediments (Amos et al. 2014). The latter system is very amenable to manipulation of genetic material as there is less of a requirement for isolation of very large DNA fragments from environmental samples.

In this study, we used a high-throughput fosmid library screening method to assess the frequencies of antibiotic resistance to a β-lactam (ampicillin), a cephalosporin (cephalothin), an aminoglycoside (kanamycin), and a tetracycline (tetracycline) in bacterial communities throughout the Mississippi River in Minnesota. Fosmid libraries consisting of approximately 10,000 clones were constructed from each of 11 sites sampled in the summers of 2011 and 2012. Physicochemical, nutrient, land, and bacterial community taxonomic data were also collected to determine how these factors were related to antibiotic resistance frequencies. We initially hypothesized that chemical parameters would have a more significant relationship with antibiotic resistance than community composition, in large part due to the persistence of a taxonomically and functionally conserved, core microbial community throughout the study area (Staley et al. 2013, 2014). Results of this study reveal how anthropogenic chemical inputs and community structure presumably influence the distribution of ARGs in a major river ecosystem.

Eight sampling sites were selected along the main branch of the Mississippi River in Minnesota from the headwaters at Lake Itasca to near the southern border at La Crescent (Figure 1). In addition, the Minnesota, St Croix, and Zumbro Rivers were sampled. The 11 sites were selected to cover the length of the Mississippi River throughout the state and characterize communities in pristine (forested), agricultural, and urban-developed areas as well as in the major confluent rivers. Each site was sampled once between May and July in both 2011 and 2012. At each site, 40 L of water was collected from the surface, approximately 1.8 m from the shoreline in two 20 L carboys and transported back to the laboratory. Water temperature and pH were also recorded at the time of sampling, and rainfall up to 3 days prior to sampling was obtained from the website http://www.wunderground.com.

Additional 1 L samples were also collected for nutrient analysis in sterile amber bottles. Determination of the concentrations of ammonium, colorimetric nitrite/nitrate (NO₂/NO₃), orthophosphate, total phosphorus, total dissolved solids (TDS), and total organic carbon was performed at the Research Analytical Laboratory at the University of Minnesota (Saint Paul) via standard methods (http://ral.cfans.umn.edu/types-of-analysis-offered/water/). Land cover data were extrapolated from the 2006 National Land Cover Database (Fry et al. 2011) by overlaying a map of hydrologic unit code (HUC) boundaries at a scale of 1:250,000 using ArcGIS (Esri, Redlands, CA, USA). Maps were obtained from the US Geological Survey (http://water.usgs.gov/maps.html). Major land cover types (forested, developed, or agricultural) were assigned based on percentage area within the HUC boundary (Table 1).

Samples were either processed immediately or stored at 15 °C for <24 h before filtration. Water was filtered as previously described (Staley et al. 2013). Briefly, water was pre-filtered through 90 mm diameter P5 filters (Whatman Inc., Piscataway, NJ, USA), and microorganisms were concentrated on 142 mm diameter, 0.45 μm polyethane-sulfonate filters (Pall Co., Port Washington, NY, USA) followed by elutriation in pyrophosphate buffer (0.1% sodium pyrophosphate buffer, pH 7.0, 0.2% Tween 20). Cell pellets (six per sample, each representing 6–7 L of water) from cell suspensions were stored at −80 °C.

One cell pellet per sample was shipped on dry ice to the Clemson University Genomics Institute (CUGI) (http://www.genome.clemson.edu/) for fosmid construction. DNA from each of the samples was extracted using the Metagenomic DNA Isolation Kit for Water (Epicentre Biotechnologies, Madison, WI, USA) followed by end-repair/phosphorylation. DNA fragments between 35–50 kb were size-selected by pulsed-field gel electrophoresis and ligated into pCC2FOS (Epicentre Biotechnologies, Madison, WI, USA). Ligated fosmids were transduced into Escherichia coli DH10B by λ phage at CUGI. Fosmid libraries for each site contained a minimum of 50,000 clones and were shipped back to the laboratory on dry ice as glycerol stocks.

Fosmid libraries were diluted to 2.5 colony forming units (CFU) μL⁻¹ and 1 mL aliquots were plated on 20 × 20 cm Luria Bertani (LB) agar plates containing 12.5 μg mL⁻¹ chloramphenicol (CAM). Colonies (approximately 10,000 per library per site per year, Table 1) were transferred to 384-well plates containing Hogness modified freezing media (Yan et al. 2007) with 12.5 μg mL⁻¹ CAM using the QBot colony picking robot (Genetix, Sunnyvale, CA, USA). Fosmid libraries were stored at −80 °C.

Fosmid libraries were thawed at room temperature immediately prior to functional screening. Functional screening was performed by plating libraries on 20 × 20 cm Müller-Hinton plates (Himedia, Mumbai, India) amended with 7 μg mL⁻¹ CAM and the experimentally determined minimal inhibitory concentrations (MICs) of antibiotic – 20 μg mL⁻¹ ampicillin (AMP), 35 μg mL⁻¹ cephalothin (CET), 15 μg mL⁻¹ kanamycin (KAN), or 10 μg mL⁻¹ tetracycline (TET; see below). Plating was performed using a flame-sterilized 384-well replicator (Boekel Scientific, Feasterville, PA, USA) and up to six 384-well plates were stamped on each 20 × 20 cm plate. For each set of plates, a negative control (E. coli DH10B containing a fosmid without insert) was also streaked. Plates were incubated overnight (16–18 h) at 37 °C. Resistant isolates were determined as those that formed opaque colonies at least 1 mm in diameter and reported as a percentage of the total clone library for each sample.

MICs of clones were determined as the lowest concentrations of antibiotic that reproducibly (triplicate cultures) prevented growth of a control strain (E. coli DH10B containing pCC2FOS without insert) grown overnight at 37 °C in 5 mL LB broth with agitation at 250 rpm. Concentrations of antibiotic were adjusted in 5 μg mL⁻¹ increments until MICs were established and all plates were amended with 7 μg mL⁻¹ CAM. All antibiotics were obtained from Sigma–Aldrich (St Louis, MO, USA).

DNA was extracted from two separate cell pellets using the Metagenomic DNA Isolation Kit for Water (Epicentre Biotechnologies). The V6 hypervariable region of the 16S rRNA gene was amplified using barcoded 967F/1046R primers as described previously (Sogin et al. 2006; Staley et al. 2013) and purified using the QiaQuick^® Gel Extraction Kit (Qiagen, Valencia, CA, USA) according to the manufacturer's instructions. Purified amplicons were pooled in equal amounts for sequencing. Amplicons originating from one cell pellet were sequenced on the Illumina MiSeq platform (2 × 150 read length) at the University of Minnesota Genomics Center (UMGC, Minneapolis). For replication and to enable greater sample multiplexing, DNA from the second cell pellet was sequenced in duplicate by UMGC on the HiSeq2000 (2 × 100 read length), as cross-platform variation in data has been previously reported (Caporaso et al. 2012). All sequences were deposited in the National Center for Biotechnology Information Sequence Read Archive under accession number SRP018728.

All sequence processing was performed using Mothur ver. 1.29.2 (Schloss et al. 2009). Sequences were trimmed to 100 bp and paired-end aligned using fastq-join (Aronesty 2013). Sequences were quality trimmed using a window of 50 bp and an average quality score of 35. In addition, singleton sequences, those containing an ambiguous base, homopolymers >8 bp, and sequences that did not have 100% identity to primer and barcode sequences were removed. Samples were aligned against the SILVA reference database ver. 102 (Pruesse et al. 2007) and subjected to a 2% precluster (Huse et al. 2010; Kunin et al. 2010). Chimeric sequences were removed using UCHIME (Edgar et al. 2011), samples were normalized to 25,717 sequence reads per sample by random subsampling, and operational taxonomic units were assigned at 97% similarity using the furthest-neighbor algorithm. Taxonomic classification was also performed against the Ribosomal Database Project database ver. 9 (Cole et al. 2009).

To compare taxonomic data with antibiotic resistance frequencies, the abundances of phyla were averaged among triplicates. Analysis of variance (ANOVA) and Spearman rank correlations relating physicochemical data, land cover, taxonomic data, and antibiotic resistance frequencies were performed using SPSS software ver. 19 (IBM, Armonk, NY, USA). Redundancy analysis was performed using XLSTAT (Addinsoft, Belmont, MA, USA). All statistics were evaluated at α = 0.05.

The frequency of antibiotic resistance to AMP, CET, and KAN was observed throughout the study area during both years. However, frequencies were extremely low, never with frequencies ≥1.0% observed among all fosmid libraries (Table 1). Resistance to TET was not observed among any of the fosmid libraries screened. Annual differences were observed in the frequency of resistance to both CET and KAN, but not AMP. In 2011, the frequency of CET resistance was significantly higher compared to 2012 (P = 0.019), but resistance to KAN was lower (P = 0.008). However, throughout the study, the frequencies of CET and KAN resistance were positively correlated (r = 0.426, P = 0.048).

Relationships between resistance frequencies and climatic and physicochemical parameters were investigated to determine if these parameters were associated with increased resistance frequencies. Few relationships were significant, although the ammonium concentration was significantly correlated with KAN resistance frequency (r = 0.617, P = 0.002), and the concentration of TDS was negatively correlated with CET resistance frequency (r = −0.449, P = 0.036). These were among the only parameters, in addition to temperature and total phosphorus concentration, which differed significantly (P < 0.05) between years (Table 2).

The surrounding land cover types of the basins under study were further interrogated to determine if specific land cover types (developed, agricultural, or forested cover) were related to increases in resistance to specific antibiotics. The frequency of resistance to CET was significantly positively correlated with the percentage of surrounding forested cover (r = 0.444, P = 0.039) and inversely correlated with agriculturally associated cover (r = −0.435, P = 0.043). Forested area also co-varied significantly with TDS concentration (r = −0.518, P = 0.014). No co-variation was observed between other parameters that were significantly associated with antibiotic resistance frequency and surrounding land cover.

Since runoff from a variety of non-point sources might be associated with the increase in the frequency of ARGs, due potentially to introduction of non-indigenous, resistant taxa, the bacterial community was characterized via 16S rRNA sequencing of the V6 hypervariable region. The community was found to be comprised primarily of α-, β-, and γ-proteobacteria, Actinobacteria, and Bacteroidetes, among 35 phyla identified (Table 2). No differences in the relative abundance of these groups were observed between years of study at α = 0.05. Among all the phyla identified, only the relative abundances of Planctomycetes (r = −0.524, P = 0.012) and candidate division OD1 (r = 0.456, P = 0.033) were significantly correlated with the frequency of AMP resistance. Resistance to CET or KAN was not associated with any of the phyla identified.

Lastly, redundancy analysis was performed incorporating all of the parameters investigated by traditional correlation analysis except rainfall, to simplify the model (Figure 2). The results of the redundancy analysis generally corroborated the results of the traditional analyses, showing weak to moderate, positive associations between ammonium and forested land cover with resistances of KAN and CET, respectively. Similarly, TDS concentration was shown to be negatively associated with the frequency of CET resistance. Interestingly, however, there was an apparent weak, but positive, association between the relative abundance of γ-proteobacteria and AMP resistance frequency, which was not captured by traditional analysis.

ARGs are now recognized as emerging contaminants of surface waters and pose a risk of conferring antibiotic resistance to human pathogens (Martinez 2009), and recent research has suggested that resistant bacteria are able to enter drinking water supplies indicating a further risk of these bacteria spreading to the human food chain (Walsh et al. 2011). The Mississippi River is used as a source of drinking water for more than 50 cities, affecting millions of people, so it is critical to protect this ecosystem from the introduction of ARGs, among other pollutants.

In this study, functional metagenomic screening of fosmid libraries revealed extremely low frequencies of ARGs in the Upper Mississippi River. The frequencies of antibiotic resistance reported here are similar to those found in a recent study of river sediments using functional metagenomic screening (Amos et al. 2014). However, our results are in contrast to a previous study that employed direct plating of water collected from the Mississippi River near Minneapolis, MN, USA (Ash et al. 2002). In the prior study, a mean concentration of 3.09 log₁₀ CFU mL⁻¹ was observed for two water samples analyzed, and 19.7 and 23.7% of isolates were resistant to AMP. Furthermore, among AMP-resistant isolates (n = 72), 38% were resistant to CET and 28% were resistant to amoxicillin. The discrepancy in the frequency of antibiotic resistance is likely due to differences in methodology, where direct plating may have selected for rare, antibiotic-resistant species. Conversely, since the community in the Mississippi River is comprised primarily of a small number of highly abundant species (Staley et al. 2013), it is likely that these species were over-represented in the fosmid libraries and did not possess ARGs. It is also possible, based on the nature of fosmid library construction, that ARGs, which are less abundant versus genes encoding metabolic enzymes, were diluted out of the fosmid libraries in favor of more abundant housekeeping genes. In addition, some genes incorporated into the fosmid libraries may simply have not been expressed by the E. coli host strain used due potentially to orientation in the fosmid vector.

To our knowledge, this study is among the first to examine the relationships between nutrient concentrations, which are potentially related to anthropogenic impacts, and antibiotic resistance frequency. While most parameters measured were poorly related to resistance frequency, ammonium and TDS concentrations were found to be significantly associated with resistance to CET and KAN. While CET resistance was observed to be higher in 2011, so were TDS concentrations, which were negatively correlated with this resistance frequency. Similarly, in 2012, KAN resistance was higher, but so were ammonium concentrations, and the two were negatively correlated. However, concentrations of CET and KAN resistance were also positively correlated, suggesting that if nutrient concentrations play a role in selecting for ARGs, the dynamics are likely complex and require further study. In addition to responses directly related to nutrient concentrations, it is also possible that biofilms potentially present on TDS particles contribute to increased frequencies of antibiotic resistance as these communities are thought to harbor high densities of antibiotic-resistant phenotypes (Marti et al. 2013).

Interestingly, both traditional correlation and redundancy analyses revealed that forested land cover was significantly associated with antibiotic-resistant phenotypes among fosmid libraries, while agricultural cover was negatively correlated with CET resistance, and no significant relationships were observed for developed land cover. Only two sites (Itasca and La Crescent) were classified as forested areas, and the Itasca site is relatively shallow and potentially more highly influenced by soil communities with intrinsic resistance mechanisms (D'Costa et al. 2006). Factors such as canopy cover that would reduce exposure to UV light may also serve to protect resistant species, but conclusive determination of factors affecting this finding remains to be studied. This result is in contrast to a prior study conducted on the South Platte River Basin, where a positive correlation between the capacities of upstream wastewater treatment plants and animal feeding operations to the frequency of sulfonamide resistance was observed using qPCR targeting sul1 (Pruden et al. 2012). Furthermore, the prior study found that distribution of TET resistance, targeting tet(W), was independent of land use. Another study investigating tertiary-treated wastewater and surface waters in the Duluth-Superior Harbor identified approximately 20-fold higher concentrations of genes encoding TET resistance (tet(A), tet(W), and tet(X)) and a gene encoding the integrase of class 1 integrons (intI1) by qPCR in wastewater compared to surface waters (LaPara et al. 2011). The lack of detection of TET in this study is surprising and may potentially be a result of the method used. However, it has been shown (De Francesco et al. 2010) that qPCR may overestimate the frequency of antibiotic resistance due, at least in part, to detection of heteroresistant organisms.

A prior study implicated the phyla Actinobacteria and Firmicutes as responsible for the transport and dissemination of ARGs in Mediterranean water reservoirs (Huerta et al. 2013). Similarly, in a study of wastewater, the phyla Bacteroidetes and Firmicutes, as well as the β-proteobacteria, were positively correlated with antibiotic-resistant populations (Novo et al. 2013). In the present study, the abundance of only one candidate division, OD1, was positively correlated with resistance to AMP, and this division has been reported to be present at relatively low abundance in the Mississippi River (Staley et al. 2013). Given the relatively low frequency of resistance detection and the potential restriction of the fosmid libraries to more abundant species, it is possible to suggest that minority members of the community in the water column are responsible for ARG dissemination, and their detection may have been limited in this study. Such minority members of the community may be harbored in greater densities in sediments and biofilms where ARGs may be more abundant (Marti et al. 2013), and antibiotic-resistant communities in these reservoirs have been suggested to protect planktonic bacteria from perturbation due to antibiotic release (Baquero et al. 2008). The interchange of ARGs and resistant phenotypes between the water column and sediment or biofilm communities will require further study.

While significant associations were observed between nutrient concentrations, land, and community structure on the distribution of ARGs in this ecosystem, elucidation of the mechanisms behind these associations could not be determined in the current study. Similarly, the extent to which ARGs identified in the Mississippi River may pose a risk to human health, either by entering drinking water or via HGT to pathogens, could not be determined based on functional metagenomic screening. However, this study provides important insights into which factors may play a role in HGT, as well as ARG selection and dissemination. To inform us further on this issue, more focused genetic and genomic studies will be necessary to determine such mechanisms. Importantly, this study provides insight into community-level distributions of ARGs and identifies factors that may be mechanistically important to these distributions.

In this study, a functional metagenomic screening strategy revealed that the frequency of resistance to AMP, CET, KAN, and TET were low-to-non-detectable. Evaluation of nutrient concentrations, land cover, and taxonomic composition of the river community suggest that ammonium and TDS may be involved in a complex selection process for ARGs and that much of the resistance observed may be of natural origin. Furthermore, we suggest that minority species, potentially contributed in low numbers from sediment and biofilm reservoirs, may be the primary carriers of ARGs in this riverine system.

Funding for this project was provided by the American Recovery and Reinvestment Act of 2009 and the Minnesota Environment and Natural Resources Trust Fund as recommended by the Legislative-Citizen Commission on Minnesota Resources. This work was carried out in part using computing resources at the University of Minnesota Supercomputing Institute.