Original Research ARTICLE
Human Activity Determines the Presence of Integron-Associated and Antibiotic Resistance Genes in Southwestern British Columbia
- 1Department of Pathology & Laboratory Medicine, The University of British Columbia, Vancouver, BC, Canada
- 2BC Centre for Disease Control Public Health Laboratory, Vancouver, BC, Canada
- 3Provincial Laboratory for Public Health, Edmonton, AB, Canada
- 4Department of Laboratory Medicine and Pathology, Faculty of Medicine & Dentistry, University of Alberta, Edmonton, AB, Canada
- 5Laboratory Services, Public Health Ontario, Toronto, ON, Canada
- 6National Microbiology Laboratory, Public Health Agency of Canada, Winnipeg, MB, Canada
- 7Coastal Genomics, Inc., Burnaby, BC, Canada
- 8Pacific Biological Station, Nanaimo, BC, Canada
- 9Department of Pathology, Sidra Medical and Research Center, Doha, Qatar
The dissemination of antibiotic resistant bacteria from anthropogenic sources into the environment poses an emerging public health threat. Antibiotic resistance genes (ARGs) and gene-capturing systems such as integron-associated integrase genes (intI) play a key role in alterations of microbial communities and the spread of antibiotic resistant bacteria into the environment. In order to assess the effect of anthropogenic activities on watersheds in southwestern British Columbia, the presence of putative antibiotic resistance and integrase genes was analyzed in the microbiome of agricultural, urban influenced, and protected watersheds. A metagenomics approach and high-throughput quantitative PCR (HT qPCR) were used to screen for elements of resistance including ARGs and intI. Metagenomic sequencing of bacterial genomic DNA was used to characterize the resistome of microbial communities present in watersheds over a 1-year period. There was a low prevalence of ARGs relative to the microbial population (<1%). Analysis of the metagenomic sequences detected a total of 60 elements of resistance including 46 ARGs, intI1, and groEL/intI1 genes and 12 quaternary ammonium compounds (qac) resistance genes across all watershed locations. The relative abundance and richness of ARGs was found to be highest in agriculture impacted watersheds compared to urban and protected watersheds. A downstream transport pattern was observed in the impacted watersheds (urban and agricultural) during dry months. Similar to other reports, this study found a strong association between intI1 and ARGs (e.g., sul1), an association which may be used as a proxy for anthropogenic activities. Chemical analysis of water samples for three major groups of antibiotics was below the detection limit. However, the high richness and gene copy numbers (GCNs) of ARGs in impacted sites suggest that the effects of effluents on microbial communities are occurring even at low concentrations of antimicrobials in the water column. Antibiotic resistance and integrase genes in a year-long metagenomic study showed that ARGs were driven mainly by environmental factors from anthropogenized sites in agriculture and urban watersheds. Environmental factors such as land-use and water quality parameters accounted for 45% of the variability observed in watershed locations.
Antibiotic resistance is recognized as a major emerging health threat worldwide. While not a new phenomenon, factors such as the global population growth, overuse, and limited development of novel antibiotics have increased overall the morbidity and mortality as well as the costs of treating of bacterial diseases (Spellberg et al., 2008; Colomer-Lluch et al., 2011).
Globally, the annual consumption of antibiotics is estimated to be 70 billion standard units/year for human use (Van Boeckel et al., 2014) and 63,151 ± 1,560 tons/year for livestock (Van Boeckel et al., 2015). Within the next 15 years, usage is predicted to increase by 30% and 67% for human and veterinary purposes, respectively (Gelbrand et al., 2015). Currently, mortality rates attributable to antimicrobial resistance are estimated to be 700,000 deaths annually (O’Neill, 2014). The same report suggests that this rate will increase to 10 million deaths/year by 2050; antimicrobial resistance will become a leading cause of death world-wide (O’Neill, 2014). While this number is being debated (de Kraker et al., 2016), an upward trend will undoubtedly continue as a consequence of the use and misuse of antibiotics.
It is estimated that between 75% and 90% of antibiotics are poorly absorbed by either humans or animal hosts and are excreted, unaltered, in feces or urine (Karthikeyan and Meyer, 2006; Sarmah et al., 2006; Bockelmann et al., 2009). Thus, spillage of antibiotics and their metabolites into the environment is widespread; contamination hotspots include hospital sewage discharges, health care facilities, and community wastewater treatment plants (WWTPs), pharmaceutical industry facilities, and confined animal feeding operations (Pruden et al., 2006; Berendonk et al., 2015). Terrestrial and aquatic ecosystems such as soil, rivers, streams, watersheds, groundwater, and sediments are the recipients of both antibiotic residues and antibiotic-resistant bacteria (Pruden et al., 2006; Baquero et al., 2008). While some antibiotics may degrade quickly, others accumulate in the soil or sediments and persist in the environment for a longer period of time (Shi et al., 2014; Xiong et al., 2015a). Surface water remains the main vehicle of dissemination of antibiotic resistant bacteria, antibiotic residues, and antibiotic resistance gene (ARG) elements into the environment (Heberer, 2002; Marti et al., 2014).
Microbial acquisition of ARGs occurs through a variety of gene transfer systems using genetic elements such as conjugative plasmids, transposons, and integrons (Bennett, 2008; Garriss et al., 2009). Moreover, horizontal gene transfer enables genes to move from one bacterial cell to another and from one ecosystem to another (Bennett, 2008). The role of elements such as phages and integrons in the spread of resistance appears to be significant in freshwater ecosystems (Lupo et al., 2012) and has been well documented using culture-independent approaches (Marti et al., 2013; Stalder et al., 2014; Xiong et al., 2014, 2015b). In particular, integrons are considered the main agents of bacterial evolution due to their role in the dissemination of ARGs, development of multiple drug resistance, and their ability to add gene structures into bacterial genomes (Mazel, 2006; Joss et al., 2009; Cambray et al., 2010; Gillings, 2014). These assembly platforms containing gene cassettes with a variety of functions, mainly antimicrobial resistance, are composed of a gene encoding integrase (intI), a recombination site (attI), and a promoter (PC) located upstream of the gene cassette (Fluit and Schmitz, 2004; Gillings, 2014). Among these elements, intI genes are considered the signatures to differentiate integron classes (Sá et al., 2010). In this context, there are at least five intI classes (Ravi et al., 2014); the first three classes of intI (intI1, intI2, and intI3) are mainly linked with horizontal gene transfer among multiple bacterial species (Shibata et al., 2003; Fluit and Schmitz, 2004; Ramirez et al., 2005; Xu et al., 2007; Gillings et al., 2008, 2015; Wright et al., 2008; Partridge et al., 2009; Rodriguez-Minguela et al., 2009; Barkovskii et al., 2010; Ravi et al., 2014), while the other two classes, intI4 and intI5, are chromosomal integrases harbored by Vibrio isolates conferring resistance to trimethoprim (Diaz-Mejia et al., 2008; Ravi et al., 2014). It is hypothesized that ancestral integrons may have not harbored multiple ARGs compared to modern integrons, which encode resistance to quaternary ammonium compounds/disinfectants, sulfonamides, and aminoglycosides and exemplify the increase in antibiotic usage (Gillings et al., 2009; Ravi et al., 2014; Adesoji et al., 2017) in a wide variety of environments.
Several metagenomic studies have been conducted in aquatic ecosystems analyzing the resistome of microbial communities from sediments, WWTPs, and effluents (Uyaguari et al., 2011; Gomez-Alvarez et al., 2012; Zhu et al., 2013; Baumlisberger et al., 2015). Only a few of these studies, however, have focused on ARGs or elements of resistance in surface water or watersheds (Li B. et al., 2015; Staley et al., 2015). Since watersheds are a key source of drinking water, they may represent a pivotal role in the spread of antibiotic resistance. As watersheds are impacted by land use and anthropogenic activities, the water they provide flows into other larger aquatic ecosystems such as rivers, lakes, basins, estuaries, and salt water bodies. In the present study, we collected raw surface water samples from three watersheds in southwestern British Columbia, Canada, over a 1-year period. Each watershed was characterized by different land use (agricultural, urban, and protected). Metagenomic sequencing of the bacterial fraction and high-throughput quantitative PCR were used to characterize and quantify elements of resistance in study water samples. Understanding the interplay between land use, watersheds, and the dissemination of ARGs or genetic elements associated to resistance, is important to understanding and planning to prevent future impacts to both the public health and the environment.
Materials and Methods
In the present study, we combined shotgun metagenomics and high-throughput quantitative PCR (HT qPCR) to detect, characterize, and quantify ARGs and elements of resistance in watersheds with different land use. Figure 1 summarizes the experimental design used in the study.
FIGURE 1. Experimental design used in the characterization and quantitation of antibiotic resistance genes from agriculture, urban, and protected watersheds in southwestern British Columbia.
A case-control design was used to characterize spatial and temporal distribution of representative antibiotic resistance and integrase genes in watersheds. Each sampling study site represented different land-use: agriculture (Agricultural Upstream or AUP site located in a forested area with minimal housing, Agricultural Polluted or APL site receiving effluents from multiple farms, and Agricultural Downstream or ADS site, 2.5 km from APL and fed by waters from APL and another tributary with both agricultural and some urban influence); urban (Urban Polluted or UPL site characterized by residential development in a mountainous forest area, and urban downstream or UDS site, 1 km from UPL and running through residential, forested, and park areas); and non-impacted (Protected Upstream or PUP site, in a protected watershed that provides drinking water for a large community, and Protected Downstream or PDS site, 16 km from PUP). This last site was located after the drinking water supply (raw reservoir water) passed through a 9-km pipe (2 m in diameter; Supplementary Figure S1). Agricultural, urban and protected watersheds are not connected; they were separated by a distance of at least 63 km. A total of 89 samples were collected during a 13-month period (April 2012 to April 2013). Samples related to urban watersheds were obtained in a 12-month period. All 40-l samples were collected in sterile plastic carboys. Samples were pre-filtered in situ using a 105-μm spectra/mesh polypropylene filter (SpectrumLabs, Rancho Dominguez, CA, United States) to remove larger debris. The filtrate was kept on ice, and transported to the British Columbia Centre for Disease Control Public Health Laboratory (BCCDC PHL) for processing. All samples were kept at 4°C until analysis. An additional 250-ml subsample for each of the 89 samples was collected in amber high-density polyethylene sterile bottles without prefiltration, transported to the laboratory and stored at -80°C for further analysis for antibiotics.
Water quality parameters were measured in situ using a YSI Professional Plus handheld multiparameter instrument (YSI Inc., Yellow Springs, OH, United States). Physico-chemical parameters of watersheds included: temperature (°C), dissolved oxygen (mg/l), specific conductivity (μS/Cm), total dissolved solids (mg/l), salinity (PSU), pressure (mmHg), and pH. Turbidity (NTU) was measured using a VWR turbidity meter model No. 66120-200 (VWR, Radnor, PA, United States). Water flow data (m3/s) was determined in situ using a Swoffer 3000 current meter (Swoffer Instruments, Seattle, WA, United States). Total coliform and Escherichia coli counts were determined using the Colilert-24 testing procedure (IDEXX Laboratories, Westbrook, ME, United States). Laboratory analysis included dissolved chloride (mg/l) and Ammonia (mg/l) using automated colorimetry (SM-4500-Cl-G) and spectrophotometry (SM-4500NH3G). Other nutrients such as chlorophyll a (Welschmeyer, 1994), orthophosphate (Murphy and Riley, 1962), nitrite, and nitrate (Wood et al., 1967) were also analyzed.
Filtration and DNA Extraction
The bacterial fraction of each water sample was captured by passing water first through a 1 μm-Envirochek HV sampling capsules (Pall Corporation, Ann Harbor, MI, United States), followed by further filtration using a 0.22-μm 142 mm Supor-200 membrane disc filters (Pall Corporation, Ann Harbor, MI, United States) that retained bacterial cells as previously described (Uyaguari-Diaz et al., 2015, 2016). Filters were cut into small strips (1 cm × 1 cm) using sterile scissors and placed into 50 ml sterile centrifuge tubes (VWR, Radnor, PA, United States). Bacterial-sized cells captured on the Supor-200 membrane disc filters were then washed with 15 ml of 1× phosphate buffered saline (PBS) and 0.01% Tween (pH 7.4). Mechanical procedures involving vigorous vortexing for 20 min in 50 ml tubes vortex adaptors and centrifugation (3,300 × g, 15 min at 4°C) were used to remove and further concentrate cells. The supernatant was removed up to the 5 ml mark, and then bacterial cells were resuspended and aliquots of 1 ml distributed into sterile 1.7 ml microcentrifuge tubes. Bacterial cells were further pelleted down at 10,000 × g, 10 min at 4°C, supernatant was removed up to ∼200 μl mark. A total of five cell aliquots per sample were then stored at -80°C for DNA extraction. DNA was extracted from multiple concentrated cell aliquots using the PowerLyzer PowerSoil DNA kit (MoBio, Carlsbad, CA, United States) by following the manufacturer’s instructions. DNA was precipitated using 10% 3M sodium acetate and 2× 100% ethanol, and 5 μl of 5 μg/μl linear acrylamide, washed with 1 ml of 70% ice-cold ethanol, and eluted in 34 μl of 10 mM Tris solution. DNA concentration and purity was assessed with Qubit fluorometer (Life Technologies, Carlsbad, CA, United States) and NanoDrop spectrophotometer (NanoDrop technologies, Inc., Wilmington, DE, United States), respectively.
To characterize microbial communities and subsequently identify antibiotic resistance elements in samples, libraries were prepared using Nextera XT DNA Sample Preparation kit (Illumina, Inc., San Diego, CA, United States). Sequencing libraries were generated using 1 ng of DNA, according to the manufacturer’s instructions, with a gel size selection modification using Coastal Genomics’ Ranger Technology (Uyaguari-Diaz et al., 2015). Metagenomic sequencing was performed on a MiSeq platform (Illumina, Inc., San Diego, CA, United States) using an Illumina MiSeq V2 2 × 250 bp paired-end reagent kit. Seventeen water samples were multiplexed on each MiSeq cartridge. To balance the diversity and account for batch effect in each sequence run, we included five samples from agricultural, five samples from urban, five samples from protected watersheds, one negative or background control, and one mock community (Peabody et al., 2015). Genomic DNA from known bacterial species were pooled in equal molar amounts and processed like indicated in the Nextera XT DNA sample preparation kit (Illumina, Inc., San Diego, CA, United States). Raw bacterial reads were created and are available as part of a large-scale Watershed Metagenomics project, BioProject ID: 2878401.
Trimmomatic version 0.32 (Bolger et al., 2014) was used to remove adapters using the sequences packaged with the A5-Miseq assembly pipeline (Coil et al., 2015). Low quality and short sequences (<75 nt) were discarded. Sequences were then assembled using PandaSeq (Masella et al., 2012); unassembled pairs were retained. De novo assembly was conducted on the assembled and unassembled pairs using MEGAHIT (Li D. et al., 2015) and contigs shorter than 200 nt were discarded. Nucleotide sequences were aligned against the comprehensive antibiotic resistance database (CARD; McArthur et al., 2013) and Integrall (Moura et al., 2009) using BLAST (Altschul et al., 1990).
Probe and Primer Design
Fully annotated genes from CARD encoding resistance to known antibiotics were used to design primers and probes for qPCR using Primer3Plus software (Untergasser et al., 2007). Specific primers for intI1-3 genes were selected from the literature (Barraud et al., 2010). The primers were first checked for sensitivity and specificity by SYBR green-based PCR (data not shown). TaqMan probes (Life Technologies, Carlsbad, CA, United States) were then designed and validated for the primers to provide suitable sensitivity and specificity (data not shown). Primers with low sensitivity and specificity were excluded from the HT qPCR. Assays with cycle threshold (Ct) values greater than 35 were not included in the HT qPCR run using the BioMark system (Fluidigm Corporation, South San Francisco, CA, United States). Table 1 summarizes primers and probes used in this study. All probes used a 5′ 6-FAM dye with an internal ZEN quencher and 3′ Iowa Black fluorescent quencher (Life Technologies, Carlsbad, CA, United States).
Genomic DNA from either agricultural or urban impacted watersheds was used as template to generate amplicons for each ARG (Table 1). Non-template controls used nuclease free-water (Promega Corporation, Fitchburg, WI, United States). DNA from purified strains of E. coli ATCC 25922, and E. coli strains JVC1076, JVC1170, and JM109 (generous gift from Drs. Davies, Miao, and Villanueva, The University of British Columbia) were used as amplification controls for 16S rRNA gene, intI1, intI2, and intI3 genes, respectively. PCR conditions were defined as follows: 94°C for 5 min, followed by 35 cycles of 94°C for 30 s, 55°C (50°C were used for intI genes, Table 1) for 45 s, 72°C for 1 min, and a final extension at 72°C for 10 min. Amplicons were visualized on a 1.5% agarose gel. PCR amplicons were purified with a QIAQuick PCR Purification Kit (Qiagen Sciences, Inc., Germantown, MD, United States) according to the manufacturer’s instructions. The purified amplicons were ligated into pCR2.1-TOPO cloning vectors (Invitrogen, Carlsbad, CA, United States) and transformed into One Shot E. coli DH5 α-T1R competent cells following the manufacturer’s protocol. Transformants were grown overnight at 37°C in lysogeny broth with 50 μg/ml of kanamycin. Plasmids were extracted and purified using PureLink Quick Plasmid Miniprep Kit (Life Technologies, Carlsbad, CA, United States), and quantified using a Qubit dsDNA high sensitivity kit on a Qubit 3.0 fluorometer (Life Technologies, Carlsbad, CA, United States). Plasmids for each ARG and integrase gene class were end-sequenced using an ABI 3130xl Genetic Analyzer (Life Technologies, Carlsbad, CA, United States) with M13 forward primer (-20) (5′-GTAAAACGACGGCCAG-3′) and M13 reverse primer (5′- CAGGAAACAGCTATGACC-3′) using BigDye Terminator version 3.1 cycle sequencing kit (Applied Biosystems, Warrington, United Kingdom). The resultant set of DNA sequences for each ARG and integrase classes were searched against the GenBank database using blastx with default settings.
Plasmid DNA harboring amplicon gene standards were linearized by digestion with the BamHI endonuclease (Life Technologies, Carlsbad, CA, United States). Serial dilutions (n = 6) of the linearized plasmids were multiplexed and used as templates to generate standard curves. Additional bioinformatic analysis screened for potential primer collision cases against the BLAST nt database using up to 1000 nt as maximum collision distance (Supplementary Spreadsheet File S1). Estimates of bacteria, antibiotic resistance, and integrase gene were determined using 16S rRNA, ARGs and intI gene fragments, respectively (Table 1). GCNs per ml of sample were calculated as previously described (Ritalahti et al., 2006).
High-Throughput Multiplex Quantitative Polymerase Chain Reaction
DNA extracts from watershed samples were diluted 10-fold and 1.25 μl of DNA from each sample was pre-amplified with low-concentration primer pairs (0.2 μM) corresponding to all assays (Table 1) in a 5 μl reaction volume using TaqMan Preamp Master Mix (Life Technologies, Carlsbad, CA, United States) according to the BioMark protocol (Fluidigm Corporation, South San Francisco, CA, United States). Unincorporated primers were removed using ExoSAP-IT High-Throughput PCR Product Clean Up (MJS BioLynx Inc., Brockville, ON, Canada) and samples were diluted 1:5 in DNA Suspension Buffer (TEKnova, Hollister, CA, United States).
The pre-amplified products were run on the BioMark system (Fluidigm Corporation, South San Francisco, CA, United States) using 96.96 dynamic arrays. Five μl of 10× assay mix (9 μM primers and 2 μM probes) were loaded to the assay inlet, while 5 μl of sample mix (2× TaqMan Mastermix; Life Technologies, Carlsbad, CA, United States), 20× GE Sample Loading Reagent, nuclease-free water, and 2.25 μl of pre-amplified DNA were loaded to each sample inlet of the array following the manufacturer’s recommendations. After mixing the assays and samples into the chip with an IFC controller HX (Fluidigm Corporation, South San Francisco, CA, United States), qPCR was performed with the following conditions: 50°C for 2 min, 95°C for 10 min, followed by 40 cycles of 95°C for 15 s and 60°C for 1 min. Samples were run in quadruplicates for all six standards, for the 89 environmental samples and for the one no-template control.
Analysis of Antibiotic Residues
A subset of water samples (n = 19) stored at -80°C was sent to a commercial laboratory to screen for antibiotic residues. This included 13 water samples from APL and one sample from the other sites (AUP, ADS, UPL, UDS, PUP, and PDS). Selection of antibiotic residues for testing was based on the major groups of ARGs found through the CARD database and included: ampicillin (β-lactam group); sulfamethoxazole (sulfonamide group); and chlortetracycline, doxycycline, oxytetracycline, and tetracycline (tetracycline group). Residues of ampicillin and sulfamethoxazole were analyzed by LC/MS following the methods described in EPA 549.2 (EPA, 1997), while the tetracycline group was tested by LC/MS–MS following a modified EPA 549.2 method (EPA, 1997).
To estimate the percentages of ARGs and intI, we used the number of contigs matching CARD and Integrall divided by the total number of contigs per site (Table 2 and Supplementary Spreadsheet File S2). In order to estimate the “true” resistance potential in the microbial community, a manual curation was conducted on these matching hits. Mechanisms such as mutation based resistance genes, target proteins for resistance genes, and regulators for resistance genes were not included in downstream analysis. GCNs and metadata were transformed using log10 function for analysis. Longitudinal analysis (PROC mixed with repeated measures) was conducted on the qPCR data for ARGs and intI to detect differences between watershed sites over time. Replicate values of all antibiotic resistance and integrase genes for each site were averaged and these values were introduced into the model. Tukey’s test was used to determine time effect among the different sites. Spearman’s rank correlation analysis was also conducted among ARGs, intI and water quality parameters. Factor analysis employed a parsimax oblique rotation that included metadata and the ratio of each ARG normalized by the 16S rRNA gene (as estimated by HT qPCR platform) measured in watershed locations. Statistical analyses were performed using Statistical Analysis System (SAS, version 9.4 for Windows). A p-value of 0.05 was assumed for all tests as a minimum level of significance. Additional taxonomic characterization was conducted with Taxonomer (Flygare et al., 2016) and microbial community diversity and richness indices were calculated using PAST (version 3.18; Hammer et al., 2001).
TABLE 2. Year-long summary statistics of bacterial metagenomes in watershed locations and relative abundance (expressed as %) of antibiotic resistance genes using 2 databases: CARD and Integrall.
Results and Discussion
Eighty-nine surface water samples among seven sampling sites located in three watersheds of differing land-use in southwestern British Columbia were analyzed. Samples were collected at regular intervals over the course of one year. Bacterial metagenomic sequencing generated 36 Gb for this dataset. CARD and Integrall were selected over other databases due to their characteristics of available ARG spectrum data and resources including functionality and features (Xavier et al., 2016). A total of 3.95 million contigs (≥200 bp) out of 150 million reads were assembled and analyzed using CARD and Integrall.
From the total sequencing reads identified by the databases (CARD and Integrall) as elements of resistance (Table 2), we observed a large proportion (∼96%) of them as housekeeping function-associated genes such as elongation factors (55.4%), DNA-directed RNA polymerase (23.8%), DNA topoisomerases/gyrases (10.6%), resistance nodulation cell division (3.2%), and transport systems (2.0%) (Supplementary Figure S2B). Although these functions can be targeted by antibiotics such as elfamycin, daptomycin/rifampin, and fluoroquinolones, ethambutol, or other drugs (Supplementary Figure S2A), identification of these genes is not in itself an indication of resistance; the universality of these structural genes as cellular protection mechanisms or regulators has been widely documented in a variety of environments (Yu and Zhang, 2012; Ma et al., 2014; Baumlisberger et al., 2015; Chen et al., 2015). Thus, we further manually curated both databases and found that the remaining percentage (∼4%) of these genes encoded resistance to at least one antibiotic. We also included tetracycline resistance genes as this antibiotic group acquires new genes through mobile plasmids or transposons (Agerso and Sandvang, 2005; Roberts and Schwarz, 2009). After the curation process, using CARD, we identified a total of 12,044 hits from 353 contigs as elements of resistance, with grouping into 46 ARG types (conferring resistance to at least one antibiotic). Using Integrall, we observed a total of 6,688 hits with 1200 contigs mostly associated to intI1 and 12 qac genes. These contigs per watershed location, time point, and database (CARD and Integrall) are provided in Supplementary Spreadsheet File S2. The results from this study describe only genes with resistance potential within microbial communities in watersheds.
Main Classes of Antibiotics and Mechanisms of Antibiotic Resistance
A more diverse group of elements of resistance conferring resistant to aminoglycosides, β-lactams, sulfonamides, tetracyclines, and multiple antibiotics was observed in agriculture impacted sites (APL and ADS) compared to AUP, urban, and protected watersheds (Figure 2). Aminoglycoside resistance genes were only detected in agricultural influenced watersheds and ranged from 7.9% (APL) up to 37.4% (UDS) (Figure 2A). Aminoglycoside acetyltransferases and dihydropteroate synthase were the most common mechanisms conferring resistance to aminoglycosides. They were found in agricultural impacted sites with average values of 14.9% and 13.3%, respectively (Figure 2B). Other aminoglycoside resistance mechanism present in agricultural impacted watersheds included streptomycin phosphotransferase (APL: 1.4% and ADS: 11.6%; Figure 2B). Genes involved in resistance to β-lactam antibiotics made up at least 2.3% of ARGs in the APL and were detected in all sites (Figure 2A). Although β-lactamases were present in all watersheds, cephalosporinases were only observed in APL and ADS with values of 1.4% and 6.1%, respectively. In the same context, genes conferring resistance to sulfonamides were only observed in APL and ADS with values of 4.5% and 31.3%, respectively. Moreover, relative abundance of tetracycline resistance genes had values of 0.4% in APL and 2.0% in ADS. Interestingly, multiple resistance as indicated by the integron-integrase genes class 1 were present in agricultural and urban influenced watersheds with average values of 42.5% and 23.4%, respectively. Furthermore, qac genes were observed in higher relative abundance in impacted watersheds (APL: 48.1% and UPL: 71.6%).
FIGURE 2. Composition of antibiotic resistance gene categories assigned by CARD in watershed locations: (A) Group of antibiotic, and (B) Mechanism of action. AUP, agricultural upstream site; APL, agricultural polluted; ADS, agricultural downstream; UPL, urban polluted; UDS, urban downstream; PUP, protected upstream; PDS, protected downstream.
In urban impacted watersheds a higher relative abundance of tetracycline resistance genes was observed compared to other watersheds. In this context, UDS had a higher relative abundance (53.8%) of tetracycline resistance genes compared to other watershed locations. Furthermore, a resistance mechanism associated to chloramphenicol acetyltransferase (1.1%; Figure 2B) could be only detected in UPL. Finally, using metagenomics we could only identify ARGs to β-lactams in protected watershed sites (Figure 2). The ubiquity of β-lactams in protected watersheds may reflect microbial communities unique to this soil/sediment as previously demonstrated in low impact or undisturbed environments (Allen et al., 2009; Agga et al., 2015) or may reflect biofilm formation (Kostakioti et al., 2013; Balcazar et al., 2015) within the pipe carrying water from a drinking water catchment to the community (samples were collected that flowed out of this pipe). Overall, our results confirmed the prevalence of antibiotic resistance determinants to these major classes of antibiotics (β-lactams, aminoglycosides, sulfonamides, and tetracyclines) and elements of resistance such as intI1 and qac genes in anthropogenic-influenced environments (Szczepanowski et al., 2009; Zhu et al., 2013; Wichmann et al., 2014; Baumlisberger et al., 2015). Gene richness and phylogenetic information associated with these mobile genetic elements are described below.
Diversity and Distribution of Antibiotic Resistance Genes
From the bacterial DNA metagenomes using CARD and Integrall, we identified a total of 46 ARGs, intI1, and groEL/intI1 genes and 12 qac genes associated with the different watersheds distributed in 28 bacterial families. Figure 3 shows the distribution of ARGs and elements of resistance in the watershed locations with the phylogenetic affiliation at family level based on original input from the databases.
FIGURE 3. Number of records of integrase and antibiotic resistance genes in (A) agricultural, (B) urban DNA, and (C) protected watershed locations. AUP, agricultural upstream site; APL, agricultural polluted; ADS, agricultural downstream; UPL, urban polluted; UDS, urban downstream; PUP, protected upstream; PDS, protected downstream. ∗Single records have been filtered in agricultural watershed location for visualization purposes. For further details refer to Supplementary Supplementary Spreadsheet File S2.
Forty-seven ARGs including intI1 and groEL/intI1 genes and five qac genes were observed in agriculture-influenced watersheds and distributed among 24 bacterial families. A large number of these ARGs and elements of resistance were observed in APL (n = 41) and ADS (n = 28), while less richness of ARGs genes was found in AUP (n = 7; Figure 3A). In AUP, intI1 were observed in 71% of the ARGs, with 53% associated to members of the Pseudomonadaceae family. A small number of genes conferring resistance to β-lactams (n = 3) and aminoglycosides (2) were observed distributed among Pseudomonadaceae, Enterobacteriaceae, Caulobacteraceae, and Ralstoniaceae families. Moreover, qac genes in AUP were observed only in members of the family Pseudomonadaceae (8.3%).
In APL, intI1 genes represented a large proportion of the elements of resistance (48%) and were distributed mainly between Enterobacteriaceae (42.1%), Pseudomonadaceae (27.5%), uncultured bacteria (10%), Aeromonadaceae (4.1%), Moraxacellaceae (3.3%), and plasmids (3.1%), (Figure 3A). Other representative ARGs in APL included aadA genes (15.3%), mainly aadA1 with 8.1%, sul1 (10.0%), ampC (4.4%), aac (2.2%), and strA and sul2 with 2.0% each. In ADS, sul1 (30.6%) distributed mostly between Enterobacteriaceae (55.5%) and Pseudomonadaceae (22.2%), was the ARG with the highest relative abundance followed by intI1 (19.7%). Other abundant genes in ADS included strB (9.5%) and aadA1 (8.8%) with the family Enterobacteriaceae representing 44.4% of these mobile genetic elements. Although detected in low relative ratios, tet(W) was found in agriculture impacted watersheds (APL and ADS) with average values of 1.3% and affiliated to members of the family Bifidobacteriaceae. Overall, we observed a wide range of ARGs in both agriculture-impacted sites, APL (n = 41) and ADS (n = 28) (Figure 3A). The relative higher abundance of genes such as intI1 and presence of sul1 in APL and ADS is consistent with previous studies reporting the prevalence of these genes in agriculture impacted environments (Byrne-Bailey et al., 2009; Gaze et al., 2011; Stalder et al., 2012; Wang N. et al., 2014; Gillings et al., 2015).
Urban sites were less diverse in ARGs compared to APL and ADS. In UPL and UDS, a total of 12 and 9 ARGs were observed, respectively. From the sequencing data, we observed that urban sites shared at least 2 ARG types between them (intI1 and blaIMP genes). Similar to the agriculture impacted sites, UPL and UDS had high relative abundances of intI1 genes (∼23.4%), which may be an indication of the prevalence of these genes in anthropogenic impacted environments (Stokes and Gillings, 2011; Stalder et al., 2012, 2014). In this context, integrase gene class 1 in UPL and UDS was distributed in similar percentages among members of Enterobacteriaceae (50%), followed by Pseudomonadaceae, Moraxellaceae, Alcaligenaceae, and Corynebacteriaceae with values of 10% each. Moreover, quaternary ammonium compound resistance genes were observed in lower proportion of 2.6% in UDS (affiliated only with Pseudomonadaceae) compared to 71.6% in UPL. Within this latter watershed, Enterobacteriaceae represented 36.8%, followed by Pseudomonadaceae and uncultured bacteria with 25% and 23.5%, respectively. We also noted the same pattern of low qac in downstream sites compared to sites directly impacted by anthropogenic activities (urban and agricultural). In this context, quaternary ammonium compound resistance genes have been reported from municipal wastewater, hospital, and livestock environments in members of Pseudomonadaceae (Wang et al., 2007; Romao et al., 2011; Zou et al., 2014; Wan and Chou, 2015), Enterobacteriaceae (Azadpour et al., 2015), Vibrionaceae (Kazama et al., 1999) as well as Enterococcaceae and Staphylococcaceae, indicating this element is present in Gram-negative and Gram-positive bacteria (Kazama et al., 1998).
Members of Pseudomonadaceae family were linked to blaIMP gene in UPL and UDS (Figure 3B). Moreover, tet genes were found in UPL to be affiliated with the family Streptomycetaceae (1.1%). Within the same watershed we also found a chloramphenicol acetyltransferase gene (catB10) associated with Pseudomonadaceae (1.1%). Members of this family and specifically the genus Pseudomonas carrying this gene have been reported in human-impacted environments (Roberts and Schwarz, 2016). Most of ARGs in UDS consisted of tet genes (60%) and were distributed in four types of ARGs: tet(O) (Campylobacteriaceae, 33.3%), tet(32), and tet(Q) with 23.8% each, both present in Clostridiaceae and Bacteroidaceae, respectively, and finally tet(W) (Bifidobacteriaceae, 19.1%) (Figure 3B). blaIMP gene was also observed in ADS and associated to Pseudomonadaceae (17.1%). More generally, in urban sites we observed determinants of antibiotic resistance such as intI, tet, and blaIMP genes, reported to be associated with wastewater and surface water receiving such effluents (Pellegrini et al., 2009; Zhang et al., 2009; Makowska et al., 2016).
When compared to agricultural and urban locations, a low richness of ARGs (three types) was observed in the protected watershed. In PUP, the blaIMP gene represented 60% of ARGs from Pseudomonadaceae. While the blaIMP gene was also observed in urban sites and agricultural watersheds, significantly different IMP variants (Docquier et al., 2003; Pellegrini et al., 2009) may have been present in each site. We could not type the exact IMP variant based on the recovered sequences. The remaining 40% of ARGs detected in PUP were blaOXA genes and distributed evenly between Pseudomonadaceae and Ralstoniaceae (Figure 3C). A class B β-lactamase LRA from uncultured bacteria was the only ARG detected in PDS using metagenomics. In this context, blaLRA gene encoding resistance to ampicillin and some members of the cephalosporin structural class such as cefamandole, ceftazidime, and cefoxitin, has been reported in remote Alaskan soil not associated with human activity (Allen et al., 2009). Although richness of ARGs in protected watersheds was low, we hypothesize that the ubiquity of β-lactamases may be related to ancestral/natural β-lactamases present in soil/sediments (Allen et al., 2009; Agga et al., 2015), wildlife (Guenther et al., 2011; Stedt et al., 2015), or biofilms (Schwartz et al., 2003; Balcazar et al., 2015). It is important to note that integrase genes were not detected in any of these protected watershed locations using a metagenomics approach.
Across watersheds, predominant groups of bacteria included Actinobacteria, Firmicutes, Bacteroidetes, and Proteobacteria (Supplementary Spreadsheet File S3) as previously reported (Peabody et al., 2017). It is important to note that within this latter phyla, Enterobacteriacea and Pseudomonadaceae represented on average only 1.1% and 2.2% of the total microbial community, respectively. Nevertheless, when combined 39 ARGs and elements of resistance were mostly distributed between Enterobacteriaceae (24 gene types) and Pseudomonadaceae (26 gene types) families. Moreover, high peaks of relative abundance were observed after T8 in ADS and APL for both families (Supplementary Figure S3). APL had higher relative abundance of Enterobacteriaceae than AUP (p = 0.0424), PUP (p = 0.0178), and PDS (p < 0.0001), while that for Pseudomonadaceae this difference was only detected between APL and PDS (p-value = 0.0021). No other significant differences were detected. These fluctuations over time are also reflected in APL and ADS as members of Pseudomonadaceae were among the top 50 OTUs observed in these locations (Supplementary Spreadsheet File S3). Overall, agriculture-impacted sites were found to have a greater richness of ARGs compared to urban and protected sites. Our metagenomics approach detected intl1, a potential indicator of gene cassettes in only the impacted sites, while neither intI2 nor intI3 were detected in any of the study sites. Although most of the bacterial community was captured across watersheds as estimated by rarefaction analysis, diversity, and richness indices (Supplementary Spreadsheet File S3), detection of functional genes including ARGs by metagenomics is impacted by the relative abundance of other microbes and sequencing depth; whereas the use of quantitative PCR is indifferent to the community structure. Thus, to complement our findings with metagenomics in watersheds further, we used high-throughput qPCR (HT qPCR) analysis.
Detection and Quantitation of Antibiotic Resistance Genes
To confirm and quantify ARG prevalence in the study watersheds, primers and probes were designed for 60 elements of resistance, based on results from this metagenomics study. Although only class 1 of intI genes was detected by metagenomics, we incorporated two other main classes of integrase genes (intI2 and intI3) that have been reported in the literature (Barraud et al., 2010). Because class 1 integron carry qac and multiple ARGs (Gaze et al., 2011), intI1 was also used as indicator of qac gene activity.
Samples collected from the same sampling locations, but not part of this study, were used to validate qPCR assays for the ARGs. The final HT qPCR panel included a total of 10 ARG, three classes of intI, and 16S rRNA gene primers and probes (Table 1) which were run on DNA templates from all study samples. The 16S rRNA gene was also included to estimate bacterial counts. A heat map of ARGs normalized by bacterial counts (as estimated by the 16S rRNA gene) shows the relative abundance and distribution of each gene in watershed locations over time (Figure 4). Most of the ARGs were estimated to be present in an order of magnitude of 1 × 10-2 (Figure 4). These findings are in agreement with other studies conducted in agriculture- and urban-influenced environments (Marti et al., 2013; Li B. et al., 2015; Sui et al., 2016). Genes such as tet(32), tet(Q), and tet(W) were found in agriculture-impacted sites and UDS and had ratios of up to 3.2 × 10-1 (relative to the 16S rRNA gene). We hypothesize that these values are related to agricultural discharges and in the UDS site, runoff from storm water could explain the high ratio of tetracycline resistance genes. These gene ratios were also high in agriculture-impacted sites. When ratios were compared across study sites, genes with more distinctive patterns were observed. For instance, genes such as aacA1, aadA1, strA, strB, sul1, and sul2 showed high ratios in the agriculture-impacted sites, while relative concentrations as high as 3.2 × 10-1 and 1.5 × 10-1 were observed for tet(32) and tet(Q), respectively in the urban-impacted watersheds. Gene tet(W) was prevalent in both urban- and agriculture-impacted sites compared to protected watersheds. We also observed that the intI1 gene was present over the whole year and in all sites (April 2012 to April 2013). Class 1 integrase genes had ratios ranging from 3.4 to 3.7 × 10-2 in ADS and APL, respectively, while those for AUP had a ∼5.8 × 10-3 relative concentration. Ratios detected in agriculture-impacted sites were relatively high compared to urban-impacted (1.1–1.9 × 10-2) and protected watersheds (1.6–7.7 × 10-3). Class 2 integrase genes were detected in agricultural sites (1.6–1.9 × 10-3) and to a lesser degree in urban sites ranging from 7.0 × 10-5 (UPL) to 1.4 × 10-4 (UDS), while in protected sites, intI2 was only detected once in PUP, at a low relative concentration (8.0 × 10-5). Although intermittently observed in the protected watershed locations, integrase class 3 was detected in all watershed locations. On average concentrations of intI3 were one order of magnitude higher in downstream sites of impacted watersheds (ADS and UDS) compared to the other sites. The prevalence of certain elements of resistance over time may serve as proxies of anthropogenic impacts, as proposed by other studies (Chen et al., 2013; He et al., 2014; Gillings et al., 2015). To elucidate differences over time of ARGs and intI in watershed sites, we conducted a longitudinal analysis in terms of sample volume and total extracted bacterial biomass.
FIGURE 4. Heat map depicting the ratio of antibiotic resistance genes and 16S rRNA gene over time (as estimated by high-throughput quantitative PCR platform) in watershed sites. AUP, agricultural upstream site; APL, agricultural polluted; ADS, agricultural downstream; UPL, urban polluted; UDS, urban downstream; PUP, protected upstream; PDS, protected downstream. Months are represented by T1 through T13 corresponding to the long-year study (April 2012 to April 2013).
Absolute numbers generated from standard curves in the HT qPCR platform were normalized per ml of sample (volume) and ng of water DNA (biomass) as described by Ritalahti et al. (2006). Due to the multicopy nature of the 16S rRNA gene, a factor of 4.3 was used to normalize bacterial counts (Lee et al., 2009). Figure 5 and Supplementary Figure S4 depict GCNs per ml of water sample and ng DNA, respectively, for all watershed study sites. Longitudinal analysis revealed striking differences between watershed sites in terms of 16S rRNA gene, ARGs, and intI genes. Copy numbers of the 16S rRNA gene in orders of magnitude of 105 (volume) were detected and found to change over time within the watersheds. When comparing sites, we observed average quantities of the 16S rRNA gene with 7.34 × 105 GCN/ml sample in PDS, followed by ADS and APL locations with 5.67 × 105 and 4.17 × 105 GCN/ml of sample, respectively. It is probable that values observed in agriculture impacted watersheds are associated with farm discharges. We propose that the high GCN in PDS may be due to biofilm in the pipe where study samples were collected (Uyaguari-Diaz et al., 2016) as noted above.
FIGURE 5. Gene copy numbers of antibiotic resistance genes, integrase gene classes 1, 2, and 3, and 16S rRNA gene per ml of environmental water sample over time in the different locations. AUP, agricultural upstream site; APL, agricultural polluted; ADS, agricultural downstream; UPL, urban polluted; UDS, urban downstream; PUP, protected upstream; PDS, protected downstream. Red dotted lines represent means for a specific time point. Number on the lower right represents p-value from the PROC mixed with repeated measures. Statistical significance was set at the 0.05 level.
In the agriculture-impacted sites, significantly higher GCNs of aacA1, aadA1, strA, strB, sul1, sul2, intI classes 1–3 were detected over time compared to the protected watersheds (Figure 5). Genes such as tet(32), tet(Q), and tet(W) were more abundant in both the agriculture- and urban-impacted sites compared to the protected and upstream sites. During this yearlong study, integrase classes 1, 2, and 3 appeared to be predominant in APL and ADS. Moreover, the means of intI genes differed (p < 0.0001) over time per watershed location. In agricultural sites, intI class 1, 2, and 3 genes were detected in at least one order of magnitude higher than in urban or protected watersheds, with values ranging from 104 to 105, 102 to 104, and 102 to 103 GCNs/ml of sample for intI1, 2, and 3, respectively. Integrase class 1 genes were detected in all sites; GCNs were lower in the protected watersheds compared to the agriculture- and urban-impacted watersheds. Integrase class 2 was mostly observed in APL and ADS over time, while intI3 was detected in higher GCNs in agricultural and urban sites. It was observed that in some sampling events intI2 (T1) or intI3 (T10–T12) were detected in PDS, but in significantly lower GCNs compared to the impacted sites. We also observed that most of these significant changes in absolute GCNs per ml of water occurred around October (T7), which coincides with the beginning of the rainy season and subsequently, greater microbial transport. For instance, higher GCNs per ml of water (after T7) were detected for genes such as aadA1, strA, strB, tet(32), tet(Q), tet(W), intI2, and intI3 (Figure 5). Although patterns and predictions of ARGs may differ among aquatic systems and gene types, significant increases in absolute GCNs of ARGs have been reported associated to rainfall and seasonality (Di Cesare et al., 2015, 2017). On the contrary and within this context, GCNs of blaIMP, were higher in urban-impacted watersheds from T1 to T7 (Figure 5); followed by the onset of the rainy season, a slight decline in GCNs of blaIMP was detected and persisted over time. Although AUP is not directly impacted by agriculture activities, the presence of some residences near the sampling site may have some influence in the blaIMP gene patterns observed.
An analysis of GCNs expressed per ng DNA (Supplementary Figure S4) revealed seasonal downstream transport of ARGs in the sampling sites. We observed that these intermittent changes occurred at different times of the year. From T1 through T8 (April to November 2012), high GCNs of aacA1, sul1, sul2, intI1, and intI2 were detected in ADS compared to APL. During the same time period (T1–T8), high GCNs per ng DNA of blaIMP, tet(32), and tet(Q) were observed in UDS compared to UPL (Supplementary Figure S4). Moreover, a downstream transport of genes such as aadA1, strA, strB, and tet(W) was observed in both agricultural and urban downstream locations. Low GCNs of ARGs were detected during T4 and T5 in all watershed locations (Figure 5 and Supplementary Figure S4), the driest months (July and August) in southwestern British Columbia. This pattern of seasonal variation was similar for most ARGs described (Supplementary Figure S4). After T8 and with the onset of the rainy season, higher GCNs were observed. No transport patterns, however, were observed downstream from the impacted sites. Seasonal changes have been documented in other freshwater ecosystems studies (Knapp et al., 2012). Overall, some gene patterns in terms of seasonality were observed between the dry (time points 1–7, corresponding to spring and summer) and rainy season (time points 8–13, corresponding to fall and winter; Figure 5); more time series studies across seasons are needed to confirm this observation.
Exploratory Factor Analysis of Watershed Microbiomes
An exploratory factor analysis was conducted using the ratios of ARGs to the 16S rRNA gene with physicochemical and biological parameters from each study watershed. Both orthogonal and oblique rotations were conducted. Figure 6 depicts an oblique parsimax rotation that best fits all variables assessed in this study. Factor 1 or “anthropogenic stressors” accounted for 45.2% of the observed variability, while that factor 2 or “aerobic conditions” accounted for 36.0% of the observed variability. Four clusters can visually be identified (Figure 6): APL and ADS; AUP and UPL; UDS; PUP; and PDS. Integrase genes class 1, sul1, sul2, strA, strB, aacA1, aadA1, and tet(W) appear to be driven by anthropogenic stressors in agriculture impacted sites (APL and ADS), UDS, and to a lesser extent UPL. Moreover, nutrients associated with land-use such as orthophosphate, nitrate, nitrite, and ammonia (Fraterrigo and Downing, 2008; Barnes and Raymond, 2010; Soranno et al., 2015) were observed within the same quadrant as these ARGs. In additional models, where metadata was excluded from the analysis or where only nutrients were incorporated into the model, intI1 and sul1 clustered within the same quadrant as anthropogenic perturbations (data not shown). Furthermore, a significant positive correlation (r = 0.8078, p < 0.0001) was detected between these two genes (Supplementary Figure S5). This finding is in agreement with previous studies linking intI1 and sul1 to anthropogenic activities in freshwater ecosystems (Pruden et al., 2012; Wang F.H. et al., 2014; Gillings et al., 2015). Note that loading of traditional markers of water quality (total coliform and E. coli) also fell within the same quadrant as agricultural impacted sites. While the presence of these indicator organisms does not necessarily indicate the presence of harmful bacteria in water, their counts in APL and ADS may be heavily associated with farm discharges, effluents from human septic systems, as well as from non-point source fecal contamination such as wildlife or human recreational activities (Harmel et al., 2010; Walters et al., 2011).
FIGURE 6. Factor analysis of antibiotic resistance genes/16S rRNA gene and environmental variables observed over time in watershed locations. AUP, agricultural upstream site; APL, agricultural polluted; ADS, agricultural downstream; UPL, urban polluted; UDS, urban downstream; PUP, protected upstream; PDS, protected downstream. Factor 1 and Factor 2 represent environmental stressors and aerobic conditions, respectively. AMR, percentage of antibiotic resistance genes found in metagenomic sequences (based on CARD and Integrall); Chl a, chlorophyll a; Cl-, dissolved chloride; DO, dissolved oxygen; Eco, E. coli counts; GC, percentage of guanine-cytosine content; NH3-N, ammonia; NO2, nitrite; NO3, nitrate; pH, potential of hydrogen; PO4, orthophosphate; Sal, salinity; SPC, specific conductivity; TC, total coliform counts; TDS, total dissolved solids; Temp, temperature; Turb, turbidity; WFR, water flow rate. Blue dashed lines represent factor loading values for antibiotic resistance genes and integron-integrase genes.
On the other hand, ratios of tet(Q) appeared to be related to aerobic conditions as observed in AUP, UPL, and PUP (Figure 6). A significant positive correlation (p < 0.0001) were detected between this gene and dissolved oxygen (DO). tet(Q) has primarily been associated with anaerobic bacteria, and in a lesser extend to aerobic and facultative anaerobic bacteria (Chung et al., 1999; Roberts and Schwarz, 2009, 2016). It is possible that the latter groups of bacteria carrying tet(Q) may have driven the relative abundance of this gene in these watersheds. Whereas intI2 was positively correlated with DO (r = 0.4214, p < 0.0001) but rarely detected in non-impacted sites (Figure 5), we hypothesize that their relative abundance in agricultural watersheds aligns with the presence of aerobic bacterial hosts of integron genes rather than anthropogenic stressors or other environmental factors. Furthermore, integrase class 2 and tet(Q) were positively correlated among them (r = 0.5450, p < 0.0001). It is possible that aerobic conditions favored the ratio of these genes compared to the overall genes in a microbial community.
When comparing ratios of tet(32) and intI3, they were not affected by environmental stressors or aerobic conditions (Supplementary Figures S6, S7). Instead, natural factors may influence their occurrence in the studied watersheds. For instance, tet(32), previously documented in only anaerobic bacteria from clinical samples (Melville et al., 2001; Warburton et al., 2009) has recently been reported to be widespread in human, animal, and environmental resistomes (Pal et al., 2016). Another example is the wide distribution of the lesser known integrase gene, intI3, in natural environments (Barkovskii et al., 2010; Uyaguari et al., 2013). It is important to mention that biofilms in PDS (as noted above) may have also had an effect on the microbial community structure and thus of the naturally occurring ARGs. The role of biofilms as reservoirs of ARGs has been well documented as has their influence on increasing resistance of bacteria to antibiotics (10–1000 times compared to free-living bacteria; Gilbert et al., 2002; Sabater et al., 2007; Balcazar et al., 2015). Finally, the ratio of blaIMP seems as well to be aligned with the natural occurrence of bacteria in watersheds. This metallo-β-lactamase gene conferring resistance to carbapenems was positively correlated (p ≤ 0.0233) with parameters such chlorophyll a and water flow rate (Figure 6 and Supplementary Figures S5–S7). Although an inverse relationship between chlorophyll a and water flow rate has been widely reported in freshwater ecosystems (Reynolds, 2000; Salmaso and Braioni, 2008; Lucas et al., 2009; Li et al., 2010; Desortova and Puncochar, 2011), this condition only applied in the more impacted watersheds (APL, ADS, UPL, and UDS). In this study, we observed a positive correlation between both parameters in lesser-impacted environments such as AUP or non-impacted watersheds such as PUP. No association could be determined in PDS, perhaps due to the influence of pipe-associated biofilms on this site. Another correlation (p = 0.0001) was observed between temperature and chlorophyll a, but this factor did not seem to have a major effect (p = 0.2733) on the population density harboring the blaIMP gene. It should be noted that absolute GCNs of blaIMP in urban, AUP, and protected watersheds were detected in at least one order of magnitude higher than APL or ADS (Figure 5 and Supplementary Figure S4). Additional analysis of phylogenetic information associated to blaIMP genes suggested that species of Pseudomonas (Figure 3) are the most probable bacterial host for this gene. Members of this genus are important phytopathogens and are also opportunistic agents of human infections (Tripathy et al., 2007). They have been reported in low nutrient or oligotrophic environments, urban environments, colonizing biofilms, and plumbing structures (Mena and Gerba, 2009). These observations may explain the relative abundance of blaIMP in non-agricultural impacted sites.
The model described in this study included only three of the major factors that may explain the overall variability of ARGs in watershed locations. Besides land-use practices, additional factors may have influenced the pattern of bacterial communities and ARGs in watersheds, such as seasonal conditions, water flow (i.e., flow in ADS is regulated by gated dams located 8.7 km further downstream), and indirect human interventions (i.e., water collected in PDS passes through a pipe) (Knapp et al., 2012; Zhang et al., 2015).
Analysis of Antibiotic Residues in Freshwater Samples
To further understand the impact of antibiotics on the aquatic environment, we screened for antibiotic residues using a subset of water samples. The selection of antibiotic metabolites screened was derived from the metagenomics analysis supported with information on the most commonly employed antibiotics for agricultural and human purposes in Canada (Government of Canada, 2015, 2016; Agunos et al., 2016). Information from antibiotics used on farms was not available for this study area (Henrich et al., 2015) compared to other studies (Li et al., 2013; Zhu et al., 2013; Zhao et al., 2016). Detection limits for the analytical methods used were relatively low for each antibiotic as follows: ampicillin (0.020 μg/l), sulfamethoxazole (0.0050 μg/l), chlortetracycline (0.025 μg/l), doxycycline (0.050 μg/l), oxytetracycline (0.010 μg/l), and tetracycline (0.025 μg/l). On analysis, none of the three groups of antibiotics were detectable in the subset of water samples and as such, further testing on the entire dataset (n = 89) was discontinued. For various reasons, it is probable that the amounts of antimicrobials used in this part of Canada are a smaller fraction of those used in other countries (i.e. China; Yin et al., 2013; Mao et al., 2015). Other factors such as a rapid degradation (including environmentally relevant conditions; Soeborg et al., 2004; Sanderson et al., 2005; Zhang et al., 2013; Li et al., 2016), formation of metabolites (even at a higher concentration than the parent molecule; Halling-Sorensen et al., 2002; Kolpin et al., 2002; Bonvin et al., 2013; Jiang et al., 2013), and distribution (into sediments rather than the water column; Li et al., 2012; Liang et al., 2013; Chen and Zhou, 2014; Xu et al., 2014) may have resulted in these antibiotics being below the detection limit. On the other hand, it is known that low concentrations of antibiotics or their metabolites have been associated with selectivity for antibiotic resistant bacteria (Gullberg et al., 2011; Sandegren, 2014). The observation of high richness of ARGs in the agriculture-impacted sites compared to the urban (∼1:7 ratio) and protected (∼1:11 ratio) sites suggests that the ARGs were more likely derived from bacteria from the effluents rather than the de novo acquisition of resistance in the naturally occurring bacteria in the water. Finding greater GCNs of ARGs quantified in agriculture and urban impacted sites compared to non-impacted watershed locations also supports this.
The total number of sequences associated with ARGs from bacterial reads were low in aquatic environments (<1%). The metagenomics approach identified a total of 46 different ARGs, one integrase class type 1, a groEL/intI1 gene, and 12 qac genes across all sampling sites. Agriculture-impacted sites contained a higher richness of ARGs compared to the urban or protected environments. Thirteen genes using HT qPCR were further screened to quantify GCNs of ARG in the study sites. We included integrase gene classes 1, 2, and 3 due to their relevance to mobile genetic elements and ARGs. We detected higher GCNs of ARGs in agriculture and urban impacted sites than in protected sites. A downstream transport pattern was identified for most of the ARGs during the dry season, while these differences became undetectable with the onset of higher precipitation in the study area. Genes such as aacA1, aadA1, strA, strB, sul1, sul2, and intI2 were more prevalent in agricultural sites, while tet(32) and tet(Q) had a higher prevalence in the urban sites. Genes tet(W) and intI1 were prevalent in both urban and agricultural settings. Moreover, an exploratory factor analysis found that there were three major contributors/drivers of ARGs in the watershed study sites: anthropogenic stressors (45.2%), aerobic conditions (36.0%), as well as natural occurrence (18.8%). The inability to detect antibiotics in the water suggests that the ARGs may have come from organisms in effluents from impacted sites. This is consistent with the resilience/stability of antibiotic resistant organisms even after they enter the environment. Although the occurrence of ARGs in these sites was low compared to the bacterial population, high richness and GCNs in agricultural sites and to a lesser extent in urban sites demonstrates the influence of anthropogenic activities on the aquatic environment.
MU-D conceived, designed, and performed the experiments and wrote the manuscript. MaC led the bioinformatic analyses and designed the experiments. ZL designed and performed the experiments. KC, MiC, SL, and WB performed the experiments. MN performed the size selection of sequencing libraries. DD and WH provided additional bioinformatic analyses. MU-D, MaC, KM, JI-R, PT, and NP designed the experiments, contributed analysis tools, guided the analyses, and aided in interpretations. All authors contributed to final revisions of the manuscript.
This work was funded by a Genome British Columbia and Genome Canada grant (LSARP-165WAT). MU-D was supported by a Mitacs Accelerate Fellowship (IT04560).
Conflict of Interest Statement
MN holds shares of Coastal Genomics, a privately owned British Columbia company offering the Ranger Technology used in this study.
The other authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
We would like to thank Belinda Wong, Frankie Tsang, and the staff members of the Environmental Microbiology Laboratory at the BCCDC Public Health Laboratory for their support. We also thank Brian Auk, Mark McCabe, and staff in the Molecular Microbiology and Genomics Program at the BCCDC Public Health Laboratory for their sequencing expertise. We would like to thank Jared R. Slobodan (Coastal Genomics) for his assistance with gel size selection. Thanks to Dr. Robert Balshaw (BCCDC/UBC), Min-Kuang Lee (BCCDC), David Andrade Laborde (Ryder System, Inc.), and Guillermo Baños Cruz (University of Guayaquil) for providing statistical advice. We would also like to thank Drs. Julian Davies, Vivian Miao, and Ivan Villanueva at The University of British Columbia for providing integrase control strains.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2018.00852/full#supplementary-material
FIGURE S1 | Sampling locations with land cover: (A) Agricultural watershed, (B) Urban watershed, and (C) Protected watershed.
FIGURE S2 | Pie chart depicting the relative abundance of antibiotic resistance gene categories identified by CARD in watershed locations: (A) Group of antibiotic, and (B) Mechanism of action. This figure includes all contigs identified by CARD that were assigned as housekeeping function-associated genes (i.e., elongation factors, DNA-directed RNA polymerase and DNA topoisomerases/gyrases, resistance nodulation cell division, and transport systems).
FIGURE S3 | Line plots of Enterobacteriaceae and Pseudomonadaceae observed in watershed locations. AUP: agricultural upstream site; APL: agricultural polluted; ADS: agricultural downstream; UPL: urban polluted; UDS: urban downstream; PUP: protected upstream; PDS: protected downstream. The x-axis represents sample collection time (T1–T13) and the y-axis shows percentages of relative abundance compared to the microbial community using shotgun metagenomics.
FIGURE S4 | Gene copy numbers of antibiotic resistance genes, integrase gene classes 1, 2, and 3, and 16S rRNA gene per ng of DNA over time in watershed locations. AUP, agricultural upstream site; APL, agricultural polluted; ADS, agricultural downstream; UPL, urban polluted; UDS, urban downstream; PUP, protected upstream; PDS, protected downstream. Red dotted lines represent mean values for a specific time point. Number on the lower right represents p-value from the PROC mixed with repeated measures. Statistical significance was set at the 0.05 level.
FIGURE S5 | Heat map showing the Spearman’s rank correlation analysis between integrase and antibiotic resistance genes and water quality parameters.
FIGURE S6 | Factor analysis of antibiotic resistance genes/16S rRNA gene and environmental variables observed over time in watershed locations. AUP, agricultural upstream site; APL, agricultural polluted; ADS, agricultural downstream; UPL, urban polluted; UDS, urban downstream; PUP, protected upstream; PDS, protected downstream. Factor 1 and Factor 3 represent environmental stressors and natural occurrence, respectively. AMR, percentage of antibiotic resistance genes found in metagenomic sequences (based on CARD); Chl a, chlorophyll a; Cl-, dissolved chloride; DO, dissolved oxygen; Eco, E. coli counts; GC, percentage of guanine-cytosine content; NH3-N, ammonia; NO2, nitrite; NO3, nitrate; pH, potential of hydrogen; PO4, orthophosphate; Sal, salinity; SPC, specific conductivity; TC, total coliform counts; TDS, total dissolved solids; Temp, temperature; Turb, turbidity; WFR, water flow rate. Blue dashed lines represent factor loading values for antibiotic resistance genes and integron-integrase genes.
FIGURE S7 | Factor analysis of antibiotic resistance genes/16S rRNA gene and environmental variables observed over time in watershed locations. AUP, agricultural upstream site; APL, agricultural polluted; ADS, agricultural downstream; UPL, urban polluted; UDS, urban downstream; PUP, protected upstream; PDS, protected downstream. Factor 2 and Factor 3 represent aerobic conditions and natural occurrence, respectively. AMR: percentage of antibiotic resistance genes found in metagenomic sequences (based on CARD); Chl a: chlorophyll a; Cl-: dissolved chloride; DO: dissolved oxygen; Eco: E. coli counts; GC: percentage of guanine-cytosine content; NH3-N: ammonia; NO2: nitrite; NO3: nitrate; pH: potential of hydrogen; PO4: orthophosphate; Sal: salinity; SPC: specific conductivity; TC: total coliform counts; TDS: total dissolved solids; Temp: temperature; Turb: turbidity; WFR: water flow rate. Blue dashed lines represent factor loading values for antibiotic resistance genes and integron-integrase genes.
SPREADSHEET FILE S1 | Summary of Primer pair collision detection.
SPREADSHEET FILE S2 | Number of contigs matching CARD and Integrall.
SPREADSHEET FILE S3 | Microbial community analysis.
Adesoji, A. T., Ogunjobi, A. A., and Olatoye, I. O. (2017). Characterization of integrons and sulfonamide resistance genes among bacteria from drinking water distribution systems in Southwestern Nigeria. Chemotherapy 62, 34–42. doi: 10.1159/000446150
Agerso, Y., and Sandvang, D. (2005). Class 1 integrons and tetracycline resistance genes in Alcaligenes, Arthrobacter, and Pseudomonas spp. isolated from pigsties and manured soil. Appl. Environ. Microbiol. 71, 7941–7947. doi: 10.1128/AEM.71.12.7941-7947.2005
Agga, G. E., Arthur, T. M., Durso, L. M., Harhay, D. M., and Schmidt, J. W. (2015). Antimicrobial-Resistant bacterial populations and antimicrobial resistance genes obtained from environments impacted by livestock and municipal waste. PLoS One 10:e0132586. doi: 10.1371/journal.pone.0132586
Agunos, A., Leìger, D., Gow, S., Carson, C., Deckert, A., Irwin, R., et al. (2016). “Antimicrobial use monitoring in Canadian broiler flocks – Results from the CIPARS farm surveillance program,” in Proceedings of the Sixty-Fifth Western Poultry Disease Conference, ed. D. Frame (Vancouver, BC: American College of Poultry Veterinarians), 9–11. Available at: https://aaap.memberclicks.net/assets/WPDC/wpdc_2016_final.pdf
Allen, H. K., Moe, L. A., Rodbumrer, J., Gaarder, A., and Handelsman, J. (2009). Functional metagenomics reveals diverse beta-lactamases in a remote Alaskan soil. ISME J. 3, 243–251. doi: 10.1038/ismej.2008.86
Azadpour, M., Nowroozi, J., Goudarzi, G. R., and Mahmoudvand, H. (2015). Presence of qacEDelta1 and cepA genes and susceptibility to a hospital biocide in clinical isolates of Klebsiella pneumoniae in Iran. Trop. Biomed. 32, 109–115.
Barkovskii, A. L., Green, C., and Hurley, D. (2010). The occurrence, spatial and temporal distribution, and environmental routes of tetracycline resistance and integrase genes in Crassostrea virginica beds. Mar. Pollut. Bull. 60, 2215–2224. doi: 10.1016/j.marpolbul.2010.08.016
Barnes, R. T., and Raymond, P. A. (2010). Land-use controls on sources and processing of nitrate in small watersheds: insights from dual isotopic analysis. Ecol. Appl. 20, 1961–1978. doi: 10.1890/08-1328.1
Barraud, O., Baclet, M. C., Denis, F., and Ploy, M. C. (2010). Quantitative multiplex real-time PCR for detecting class 1, 2 and 3 integrons. J. Antimicrob. Chemother. 65, 1642–1645. doi: 10.1093/jac/dkq167
Baumlisberger, M., Youssar, L., Schilhabel, M. B., and Jonas, D. (2015). Influence of a non-hospital medical care facility on antimicrobial resistance in wastewater. PLoS One 10:e0122635. doi: 10.1371/journal.pone.0122635
Bennett, P. M. (2008). Plasmid encoded antibiotic resistance: acquisition and transfer of antibiotic resistance genes in bacteria. Br. J. Pharmacol. 153(Suppl. 1), S347–S357. doi: 10.1038/sj.bjp.0707607
Berendonk, T. U., Manaia, C. M., Merlin, C., Fatta-Kassinos, D., Cytryn, E., Walsh, F., et al. (2015). Tackling antibiotic resistance: the environmental framework. Nat. Rev. Microbiol. 13, 310–317. doi: 10.1038/nrmicro3439
Bockelmann, U., Dorries, H. H., Ayuso-Gabella, M. N., Salgot de Marcay, M., Tandoi, V., Levantesi, C., et al. (2009). Quantitative PCR monitoring of antibiotic resistance genes and bacterial pathogens in three European artificial groundwater recharge systems. Appl. Environ. Microbiol. 75, 154–163. doi: 10.1128/AEM.01649-08
Bonvin, F., Omlin, J., Rutler, R., Schweizer, W. B., Alaimo, P. J., Strathmann, T. J., et al. (2013). Direct photolysis of human metabolites of the antibiotic sulfamethoxazole: evidence for abiotic back-transformation. Environ. Sci. Technol. 47, 6746–6755. doi: 10.1021/es303777k
Byrne-Bailey, K. G., Gaze, W. H., Kay, P., Boxall, A. B., Hawkey, P. M., and Wellington, E. M. (2009). Prevalence of sulfonamide resistance genes in bacterial isolates from manured agricultural soils and pig slurry in the United Kingdom. Antimicrob. Agents Chemother. 53, 696–702. doi: 10.1128/AAC.00652-07
Chen, B., Liang, X., Huang, X., Zhang, T., and Li, X. (2013). Differentiating anthropogenic impacts on ARGs in the Pearl River Estuary by using suitable gene indicators. Water Res. 47, 2811–2820. doi: 10.1016/j.watres.2013.02.042
Chen, K., and Zhou, J. L. (2014). Occurrence and behavior of antibiotics in water and sediments from the Huangpu River. Shanghai, China. Chemosphere 95, 604–612. doi: 10.1016/j.chemosphere.2013.09.119
Chen, L. X., Hu, M., Huang, L. N., Hua, Z. S., Kuang, J. L., Li, S. J., et al. (2015). Comparative metagenomic and metatranscriptomic analyses of microbial communities in acid mine drainage. ISME J. 9, 1579–1592. doi: 10.1038/ismej.2014.245
Chung, W. O., Young, K., Leng, Z., and Roberts, M. C. (1999). Mobile elements carrying ermF and tetQ genes in gram-positive and gram-negative bacteria. J. Antimicrob. Chemother. 44, 329–335. doi: 10.1093/jac/44.3.329
Di Cesare, A., Eckert, E. M., Rogora, M., and Corno, G. (2017). Rainfall increases the abundance of antibiotic resistance genes within a riverine microbial community. Environ. Pollut. 226, 473–478. doi: 10.1016/j.envpol.2017.04.036
Di Cesare, A., Eckert, E. M., Teruggi, A., Fontaneto, D., Bertoni, R., Callieri, C., et al. (2015). Constitutive presence of antibiotic resistance genes within the bacterial community of a large subalpine lake. Mol. Ecol. 24, 3888–3900. doi: 10.1111/mec.13293
Diaz-Mejia, J. J., Amabile-Cuevas, C. F., Rosas, I., and Souza, V. (2008). An analysis of the evolutionary relationships of integron integrases, with emphasis on the prevalence of class 1 integrons in Escherichia coli isolates from clinical and environmental origins. Microbiology 154(Pt 1), 94–102. doi: 10.1099/mic.0.2007/008649-0
Docquier, J. D., Riccio, M. L., Mugnaioli, C., Luzzaro, F., Endimiani, A., Toniolo, A., et al. (2003). IMP-12, a new plasmid-encoded metallo-beta-lactamase from a Pseudomonas putida clinical isolate. Antimicrob. Agents Chemother. 47, 1522–1528. doi: 10.1128/AAC.47.5.1522-1528.2003
EPA (1997). Method 549.2. Determination of Diquat and Paraquat in Drinking Water by Liquid-Solid Extraction and High Performance Liquid Chromatography with Ultraviolet Detection, ed. J. W. M. A. W. J. Bashe. Cincinnati, OH: U.S. Environmental Protection Agency.
Flygare, S., Simmon, K., Miller, C., Qiao, Y., Kennedy, B., Di Sera, T., et al. (2016). Taxonomer: an interactive metagenomics analysis portal for universal pathogen detection and host mRNA expression profiling. Genome Biol. 17:111. doi: 10.1186/s13059-016-0969-1
Gaze, W. H., Zhang, L., Abdouslam, N. A., Hawkey, P. M., Calvo-Bado, L., Royle, J., et al. (2011). Impacts of anthropogenic activity on the ecology of class 1 integrons and integron-associated genes in the environment. ISME J. 5, 1253–1261. doi: 10.1038/ismej.2011.15
Gilbert, P., Maira-Litran, T., McBain, A. J., Rickard, A. H., and Whyte, F. W. (2002). The physiology and collective recalcitrance of microbial biofilm communities. Adv. Microb. Physiol. 46, 202–256. doi: 10.1016/S0065-2911(02)46005-5
Gillings, M. R., Gaze, W. H., Pruden, A., Smalla, K., Tiedje, J. M., and Zhu, Y. G. (2015). Using the class 1 integron-integrase gene as a proxy for anthropogenic pollution. ISME J. 9, 1269–1279. doi: 10.1038/ismej.2014.226
Gillings, M. R., Krishnan, S., Worden, P. J., and Hardwick, S. A. (2008). Recovery of diverse genes for class 1 integron-integrases from environmental DNA samples. FEMS Microbiol. Lett. 287, 56–62. doi: 10.1111/j.1574-6968.2008.01291.x
Gillings, M. R., Xuejun, D., Hardwick, S. A., Holley, M. P., and Stokes, H. W. (2009). Gene cassettes encoding resistance to quaternary ammonium compounds: a role in the origin of clinical class 1 integrons? ISME J. 3, 209–215. doi: 10.1038/ismej.2008.98
Gomez-Alvarez, V., Revetta, R. P., and Santo Domingo, J. W. (2012). Metagenomic analyses of drinking water receiving different disinfection treatments. Appl. Environ. Microbiol. 78, 6095–6102. doi: 10.1128/AEM.01018-12
Government of Canada (2015). Canadian Integrated Program for Antimicrobial Resistance Surveillance (CIPARS) Annual Report 2012. Antimicrobial use in Animals. Guelph, ON: Public Health Agency of Canada, 49.
Guenther, S., Ewers, C., and Wieler, L. H. (2011). Extended-spectrum beta-lactamases producing E. coli in wildlife, yet another form of environmental pollution? Front. Microbiol. 2:246. doi: 10.3389/fmicb.2011.00246
Gullberg, E., Cao, S., Berg, O. G., Ilback, C., Sandegren, L., Hughes, D., et al. (2011). Selection of resistant bacteria at very low antibiotic concentrations. PLoS Pathog. 7:e1002158. doi: 10.1371/journal.ppat.1002158
Halling-Sorensen, B., Sengelov, G., and Tjornelund, J. (2002). Toxicity of tetracyclines and tetracycline degradation products to environmentally relevant bacteria, including selected tetracycline-resistant bacteria. Arch. Environ. Contam. Toxicol. 42, 263–271. doi: 10.1007/s00244-001-0017-2
Harmel, R. D., Karthikeyan, R., Gentry, T., and Srinivasan, R. (2010). Effects of agricultural management, land use, and watershed scale on E. coli concentrations in runoff and streamflow. Trans. ASABE 53, 1833–1841. doi: 10.13031/2013.35809
He, L. Y., Liu, Y. S., Su, H. C., Zhao, J. L., Liu, S. S., Chen, J., et al. (2014). Dissemination of antibiotic resistance genes in representative broiler feedlots environments: identification of indicator ARGs and correlations with environmental variables. Environ. Sci. Technol. 48, 13120–13129. doi: 10.1021/es5041267
Henrich, N., Holmes, B., and Prystajecky, N. (2015). Looking upstream: findings from focus groups on public perceptions of source water quality in British Columbia. Canada. PLoS One 10:e0141533. doi: 10.1371/journal.pone.0141533
Jiang, H., Zhang, D., Xiao, S., Geng, C., and Zhang, X. (2013). Occurrence and sources of antibiotics and their metabolites in river water, WWTPs, and swine wastewater in Jiulongjiang River basin, south China. Environ. Sci. Pollut. Res. Int. 20, 9075–9083. doi: 10.1007/s11356-013-1924-2
Joss, M. J., Koenig, J. E., Labbate, M., Polz, M. F., Gillings, M. R., Stokes, H. W., et al. (2009). ACID: annotation of cassette and integron data. BMC Bioinformatics 10:118. doi: 10.1186/1471-2105-10-118
Kazama, H., Hamashima, H., Sasatsu, M., and Arai, T. (1999). Characterization of the antiseptic-resistance gene qacE delta 1 isolated from clinical and environmental isolates of Vibrio parahaemolyticus and Vibrio cholerae non-O1. FEMS Microbiol. Lett. 174, 379–384.
Knapp, C. W., Lima, L., Olivares-Rieumont, S., Bowen, E., Werner, D., and Graham, D. W. (2012). Seasonal variations in antibiotic resistance gene transport in the Almendares river, havana, cuba. Front. Microbiol. 3:396. doi: 10.3389/fmicb.2012.00396
Kolpin, D. W., Furlong, E. T., Meyer, M. T., Thurman, E. M., Zaugg, S. D., Barber, L. B., et al. (2002). Pharmaceuticals, hormones, and other organic wastewater contaminants in U.S. streams, 1999–2000: a national reconnaissance. Environ. Sci. Technol. 36, 1202–1211. doi: 10.1021/es011055j
Kostakioti, M., Hadjifrangiskou, M., and Hultgren, S. J. (2013). Bacterial biofilms: development, dispersal, and therapeutic strategies in the dawn of the postantibiotic era. Cold Spring Harb. Perspect. Med. 3:a010306. doi: 10.1101/cshperspect.a010306
Li, B., Yang, Y., Ma, L., Ju, F., Guo, F., Tiedje, J. M., et al. (2015). Metagenomic and network analysis reveal wide distribution and co-occurrence of environmental antibiotic resistance genes. ISME J. 9, 2490–2502. doi: 10.1038/ismej.2015.59
Li, D., Liu, C. M., Luo, R., Sadakane, K., and Lam, T. W. (2015). MEGAHIT: an ultra-fast single-node solution for large and complex metagenomics assembly via succinct de Bruijn graph. Bioinformatics 31, 1674–1676. doi: 10.1093/bioinformatics/btv033
Li, F. P., Zhang, H. P., Zhu, Y. P., Chen, L., and Zhao, J. F. (2010). “Spatial and temporal dynamics in the relationship of phytoplankton biomass and limnological variables in a small artificial lake,” in Proceedings of the 2nd International Symposium on Aqua Science, Water Resource and Low Carbon Energy, Sanya, 29–32. doi: 10.1063/1.3456271
Li, L., Sun, J., Liu, B., Zhao, D., Ma, J., Deng, H., et al. (2013). Quantification of lincomycin resistance genes associated with lincomycin residues in waters and soils adjacent to representative swine farms in China. Front. Microbiol. 4:364. doi: 10.3389/fmicb.2013.00364
Li, W., Shi, Y., Gao, L., Liu, J., and Cai, Y. (2012). Occurrence of antibiotics in water, sediments, aquatic plants, and animals from Baiyangdian Lake in North China. Chemosphere 89, 1307–1315. doi: 10.1016/j.chemosphere.2012.05.079
Li, Y., Wang, H., Liu, X., Zhao, G., and Sun, Y. (2016). Dissipation kinetics of oxytetracycline, tetracycline, and chlortetracycline residues in soil. Environ. Sci. Pollut. Res. Int. 23, 13822–13831. doi: 10.1007/s11356-016-6513-8
Liang, X., Chen, B., Nie, X., Shi, Z., Huang, X., and Li, X. (2013). The distribution and partitioning of common antibiotics in water and sediment of the Pearl River Estuary. South China. Chemosphere 92, 1410–1416. doi: 10.1016/j.chemosphere.2013.03.044
Lucas, L. V., Thompson, J. K., and Brown, L. R. (2009). Why are diverse relationships observed between phytoplankton biomass and transport time? Limnol. Oceanogr. 54, 381–390. doi: 10.4319/lo.2009.54.1.0381
Lupo, A., Coyne, S., and Berendonk, T. U. (2012). Origin and evolution of antibiotic resistance: the common mechanisms of emergence and spread in water bodies. Front. Microbiol. 3:18. doi: 10.3389/fmicb.2012.00018
Ma, L., Li, B., and Zhang, T. (2014). Abundant rifampin resistance genes and significant correlations of antibiotic resistance genes and plasmids in various environments revealed by metagenomic analysis. Appl. Microbiol. Biotechnol. 98, 5195–5204. doi: 10.1007/s00253-014-5511-3
Makowska, N., Koczura, R., and Mokracka, J. (2016). Class 1 integrase, sulfonamide and tetracycline resistance genes in wastewater treatment plant and surface water. Chemosphere 144, 1665–1673. doi: 10.1016/j.chemosphere.2015.10.044
Marti, E., Jofre, J., and Balcazar, J. L. (2013). Prevalence of antibiotic resistance genes and bacterial community composition in a river influenced by a wastewater treatment plant. PLoS One 8:e78906. doi: 10.1371/journal.pone.0078906
Masella, A. P., Bartram, A. K., Truszkowski, J. M., Brown, D. G., and Neufeld, J. D. (2012). PANDAseq: paired-end assembler for illumina sequences. BMC Bioinformatics 13:31. doi: 10.1186/1471-2105-13-31
McArthur, A. G., Waglechner, N., Nizam, F., Yan, A., Azad, M. A., Baylay, A. J., et al. (2013). The comprehensive antibiotic resistance database. Antimicrob. Agents Chemother. 57, 3348–3357. doi: 10.1128/AAC.00419-13
Melville, C. M., Scott, K. P., Mercer, D. K., and Flint, H. J. (2001). Novel tetracycline resistance gene, tet(32), in the Clostridium-related human colonic anaerobe K10 and its transmission in vitro to the rumen anaerobe Butyrivibrio fibrisolvens. Antimicrob. Agents Chemother. 45, 3246–3249. doi: 10.1128/AAC.45.11.3246-3249.2001
Moura, A., Soares, M., Pereira, C., Leitao, N., Henriques, I., and Correia, A. (2009). INTEGRALL: a database and search engine for integrons, integrases and gene cassettes. Bioinformatics 25, 1096–1098. doi: 10.1093/bioinformatics/btp105
O’Neill, J. (2014). Review on Antimicrobial Resistance. Antimicrobial Resistance: Tackling a Crisis for the Health and Wealth of Nations. London, 5–6. Available at: http://www.jpiamr.eu/wp-content/uploads/2014/12/AMR-Review-Paper-Tackling-a-crisis-for-the-health-and-wealth-of-nations_1-2.pdf
Partridge, S. R., Tsafnat, G., Coiera, E., and Iredell, J. R. (2009). Gene cassettes and cassette arrays in mobile resistance integrons. FEMS Microbiol. Rev. 33, 757–784. doi: 10.1111/j.1574-6976.2009.00175.x
Peabody, M. A., Caravas, J. A., Morrison, S. S., Mercante, J. W., Prystajecky, N. A., Raphael, B. H., et al. (2017). Characterization of Legionella species from watersheds in British Columbia, Canada. mSphere 2:e00246-17. doi: 10.1128/mSphere.00246-17
Peabody, M. A., Van Rossum, T., Lo, R., and Brinkman, F. S. (2015). Evaluation of shotgun metagenomics sequence classification methods using in silico and in vitro simulated communities. BMC Bioinformatics 16:363. doi: 10.1186/s12859-015-0788-5
Pellegrini, C., Mercuri, P. S., Celenza, G., Galleni, M., Segatore, B., Sacchetti, E., et al. (2009). Identification of bla(IMP-22) in Pseudomonas spp. in urban wastewater and nosocomial environments: biochemical characterization of a new IMP metallo-enzyme variant and its genetic location. J. Antimicrob. Chemother. 63, 901–908. doi: 10.1093/jac/dkp061
Pruden, A., Arabi, M., and Storteboom, H. N. (2012). Correlation between upstream human activities and riverine antibiotic resistance genes. Environ. Sci. Technol. 46, 11541–11549. doi: 10.1021/es302657r
Pruden, A., Pei, R., Storteboom, H., and Carlson, K. H. (2006). Antibiotic resistance genes as emerging contaminants: studies in northern Colorado. Environ. Sci. Technol. 40, 7445–7450. doi: 10.1021/es060413l
Ramirez, M. S., Vargas, L. J., Cagnoni, V., Tokumoto, M., and Centron, D. (2005). Class 2 integron with a novel cassette array in a Burkholderia cenocepacia isolate. Antimicrob. Agents Chemother. 49, 4418–4420. doi: 10.1128/AAC.49.10.4418-4420.2005
Ravi, A., Avershina, E., Ludvigsen, J., L’Abee-Lund, T. M., and Rudi, K. (2014). Integrons in the intestinal microbiota as reservoirs for transmission of antibiotic resistance genes. Pathogens 3, 238–248. doi: 10.3390/pathogens3020238
Ritalahti, K. M., Amos, B. K., Sung, Y., Wu, Q., Koenigsberg, S. S., and Loffler, F. E. (2006). Quantitative PCR targeting 16S rRNA and reductive dehalogenase genes simultaneously monitors multiple Dehalococcoides strains. Appl. Environ. Microbiol. 72, 2765–2774. doi: 10.1128/AEM.72.4.2765-2774.2006
Roberts, M. C., and Schwarz, S. (2009). “Tetracycline and chloramphenicol resistance mechanisms,” in Antimicrobial Drug Resistance: Mechanisms of Drug Resistance, ed. D. L. Mayers (Totowa, NJ: Humana Press), 183–193. doi: 10.1007/978-1-59745-180-2_15
Roberts, M. C., and Schwarz, S. (2016). Tetracycline and phenicol resistance genes and mechanisms: importance for agriculture, the environment, and humans. J. Environ. Qual. 45, 576–592. doi: 10.2134/jeq2015.04.0207
Rodriguez-Minguela, C. M., Apajalahti, J. H., Chai, B., Cole, J. R., and Tiedje, J. M. (2009). Worldwide prevalence of class 2 integrases outside the clinical setting is associated with human impact. Appl. Environ. Microbiol. 75, 5100–5110. doi: 10.1128/AEM.00133-09
Romao, C., Miranda, C. A., Silva, J., Mandetta Clementino, M., de Filippis, I., and Asensi, M. (2011). Presence of qacEDelta1 gene and susceptibility to a hospital biocide in clinical isolates of Pseudomonas aeruginosa resistant to antibiotics. Curr. Microbiol. 63, 16–21. doi: 10.1007/s00284-011-9934-0
Sá, L. L., Fonseca, E. L., Pellegrini, M., Freitas, F., Loureiro, E. C., and Vicente, A. C. (2010). Occurrence and composition of class 1 and class 2 integrons in clinical and environmental O1 and non-O1/non-O139 Vibrio cholerae strains from the Brazilian Amazon. Mem. Inst. Oswaldo Cruz 105, 229–232. doi: 10.1590/S0074-02762010000200021
Sabater, S., Guasch, H., Ricart, M., Romani, A., Vidal, G., Klunder, C., et al. (2007). Monitoring the effect of chemicals on biological communities. The biofilm as an interface. Anal. Bioanal. Chem. 387, 1425–1434. doi: 10.1007/s00216-006-1051-8
Salmaso, N., and Braioni, M. G. (2008). Factors controlling the seasonal development and distribution of the phytoplankton community in the lowland course of a large river in Northern Italy (River Adige). Aquat. Ecol. 42, 533–545. doi: 10.1007/s10452-007-9135-x
Sanderson, H., Ingerslev, F., Brain, R. A., Halling-Sorensen, B., Bestari, J. K., Wilson, C. J., et al. (2005). Dissipation of oxytetracycline, chlortetracycline, tetracycline and doxycycline using HPLC-UV and LC/MS/MS under aquatic semi-field microcosm conditions. Chemosphere 60, 619–629. doi: 10.1016/j.chemosphere.2005.01.035
Sarmah, A. K., Meyer, M. T., and Boxall, A. B. (2006). A global perspective on the use, sales, exposure pathways, occurrence, fate and effects of veterinary antibiotics (VAs) in the environment. Chemosphere 65, 725–759. doi: 10.1016/j.chemosphere.2006.03.026
Schwartz, T., Kohnen, W., Jansen, B., and Obst, U. (2003). Detection of antibiotic-resistant bacteria and their resistance genes in wastewater, surface water, and drinking water biofilms. FEMS Microbiol. Ecol. 43, 325–335. doi: 10.1111/j.1574-6941.2003.tb01073.x
Shi, H., Yang, Y., Liu, M., Yan, C., Yue, H., and Zhou, J. (2014). Occurrence and distribution of antibiotics in the surface sediments of the Yangtze Estuary and nearby coastal areas. Mar. Pollut. Bull. 83, 317–323. doi: 10.1016/j.marpolbul.2014.04.034
Shibata, N., Doi, Y., Yamane, K., Yagi, T., Kurokawa, H., Shibayama, K., et al. (2003). PCR typing of genetic determinants for metallo-beta-lactamases and integrases carried by gram-negative bacteria isolated in Japan, with focus on the class 3 integron. J. Clin. Microbiol. 41, 5407–5413. doi: 10.1128/JCM.41.12.5407-5413.2003
Soeborg, T., Ingerslev, F., and Halling-Sorensen, B. (2004). Chemical stability of chlortetracycline and chlortetracycline degradation products and epimers in soil interstitial water. Chemosphere 57, 1515–1524. doi: 10.1016/j.chemosphere.2004.09.020
Soranno, P. A., Cheruvelil, K. S., Wagner, T., Webster, K. E., and Bremigan, M. T. (2015). Effects of land use on lake nutrients: the importance of scale, hydrologic connectivity, and region. PLoS One 10:e0135454. doi: 10.1371/journal.pone.0135454
Spellberg, B., Guidos, R., Gilbert, D., Bradley, J., Boucher, H. W., Scheld, W. M., et al. (2008). The epidemic of antibiotic-resistant infections: a call to action for the medical community from the Infectious Diseases Society of America. Clin. Infect. Dis. 46, 155–164. doi: 10.1086/524891
Stalder, T., Barraud, O., Casellas, M., Dagot, C., and Ploy, M. C. (2012). Integron involvement in environmental spread of antibiotic resistance. Front. Microbiol. 3:119. doi: 10.3389/fmicb.2012.00119
Stalder, T., Barraud, O., Jove, T., Casellas, M., Gaschet, M., Dagot, C., et al. (2014). Quantitative and qualitative impact of hospital effluent on dissemination of the integron pool. ISME J. 8, 768–777. doi: 10.1038/ismej.2013.189
Staley, C., Gould, T. J., Wang, P., Phillips, J., Cotner, J. B., and Sadowsky, M. J. (2015). High-throughput functional screening reveals low frequency of antibiotic resistance genes in DNA recovered from the Upper Mississippi River. J. Water Health 13, 693–703. doi: 10.2166/wh.2014.215
Stedt, J., Bonnedahl, J., Hernandez, J., Waldenstrom, J., McMahon, B. J., Tolf, C., et al. (2015). Carriage of CTX-M type extended spectrum beta-lactamases (ESBLs) in gulls across Europe. Acta Vet. Scand. 57:74. doi: 10.1186/s13028-015-0166-3
Stokes, H. W., and Gillings, M. R. (2011). Gene flow, mobile genetic elements and the recruitment of antibiotic resistance genes into Gram-negative pathogens. FEMS Microbiol. Rev. 35, 790–819. doi: 10.1111/j.1574-6976.2011.00273.x
Sui, Q., Zhang, J., Chen, M., Tong, J., Wang, R., and Wei, Y. (2016). Distribution of antibiotic resistance genes (ARGs) in anaerobic digestion and land application of swine wastewater. Environ. Pollut. 213, 751–759. doi: 10.1016/j.envpol.2016.03.038
Szczepanowski, R., Linke, B., Krahn, I., Gartemann, K. H., Gutzkow, T., Eichler, W., et al. (2009). Detection of 140 clinically relevant antibiotic-resistance genes in the plasmid metagenome of wastewater treatment plant bacteria showing reduced susceptibility to selected antibiotics. Microbiology 155(Pt 7), 2306–2319. doi: 10.1099/mic.0.028233-0
Tripathy, S., Kumar, N., Mohanty, S., Samanta, M., Mandal, R. N., and Maiti, N. K. (2007). Characterisation of Pseudomonas aeruginosa isolated from freshwater culture systems. Microbiol. Res. 162, 391–396. doi: 10.1016/j.micres.2006.08.005
Untergasser, A., Nijveen, H., Rao, X., Bisseling, T., Geurts, R., and Leunissen, J. A. (2007). Primer3Plus, an enhanced web interface to Primer3. Nucleic Acids Res. 35, W71–W74. doi: 10.1093/nar/gkm306
Uyaguari, M. I., Fichot, E. B., Scott, G. I., and Norman, R. S. (2011). Characterization and quantitation of a novel beta-lactamase gene found in a wastewater treatment facility and the surrounding coastal ecosystem. Appl. Environ. Microbiol. 77, 8226–8233. doi: 10.1128/AEM.02732-10
Uyaguari, M. I., Scott, G. I., and Norman, R. S. (2013). Abundance of class 1-3 integrons in South Carolina estuarine ecosystems under high and low levels of anthropogenic influence. Mar. Pollut. Bull. 76, 77–84. doi: 10.1016/j.marpolbul.2013.09.027
Uyaguari-Diaz, M. I., Chan, M., Chaban, B. L., Croxen, M. A., Finke, J. F., Hill, J. E., et al. (2016). A comprehensive method for amplicon-based and metagenomic characterization of viruses, bacteria, and eukaryotes in freshwater samples. Microbiome 4:20. doi: 10.1186/s40168-016-0166-1
Uyaguari-Diaz, M. I., Slobodan, J. R., Nesbitt, M. J., Croxen, M. A., Isaac-Renton, J., Prystajecky, N. A., et al. (2015). Automated gel size selection to improve the quality of next-generation sequencing libraries prepared from environmental water samples. J. Vis. Exp. 2015:e52685. doi: 10.3791/52685
Van Boeckel, T. P., Brower, C., Gilbert, M., Grenfell, B. T., Levin, S. A., Robinson, T. P., et al. (2015). Global trends in antimicrobial use in food animals. Proc. Natl. Acad. Sci. U.S.A. 112, 5649–5654. doi: 10.1073/pnas.1503141112
Van Boeckel, T. P., Gandra, S., Ashok, A., Caudron, Q., Grenfell, B. T., Levin, S. A., et al. (2014). Global antibiotic consumption 2000 to 2010: an analysis of national pharmaceutical sales data. Lancet Infect. Dis. 14, 742–750. doi: 10.1016/S1473-3099(14)70780-7
Walters, S. P., Thebo, A. L., and Boehm, A. B. (2011). Impact of urbanization and agriculture on the occurrence of bacterial pathogens and stx genes in coastal waterbodies of central California. Water Res. 45, 1752–1762. doi: 10.1016/j.watres.2010.11.032
Wan, M. T., and Chou, C. C. (2015). Class 1 integrons and the antiseptic resistance gene (qacEDelta1) in municipal and swine slaughterhouse wastewater treatment plants and wastewater-associated methicillin-resistant Staphylococcus aureus. Int. J. Environ. Res. Public Health 12, 6249–6260. doi: 10.3390/ijerph120606249
Wang, C., Cai, P., Guo, Y., and Mi, Z. (2007). Distribution of the antiseptic-resistance genes qacEDelta1 in 331 clinical isolates of Pseudomonas aeruginosa in China. J. Hosp. Infect. 66, 93–95. doi: 10.1016/j.jhin.2007.01.012
Wang, F. H., Qiao, M., Lv, Z. E., Guo, G. X., Jia, Y., Su, Y. H., et al. (2014). Impact of reclaimed water irrigation on antibiotic resistance in public parks. Beijing, China. Environ. Pollut. 184, 247–253. doi: 10.1016/j.envpol.2013.08.038
Wang, N., Yang, X., Jiao, S., Zhang, J., Ye, B., and Gao, S. (2014). Sulfonamide-resistant bacteria and their resistance genes in soils fertilized with manures from Jiangsu Province. Southeastern China. PLoS One 9:e112626. doi: 10.1371/journal.pone.0112626
Warburton, P., Roberts, A. P., Allan, E., Seville, L., Lancaster, H., and Mullany, P. (2009). Characterization of tet(32) genes from the oral metagenome. Antimicrob. Agents Chemother. 53, 273–276. doi: 10.1128/AAC.00788-08
Wood, E. D., Armstron, F. A. J., and Richards, F. A. (1967). Determination of nitrate in sea water by cadmium-copper reduction to nitrite. J. Mar. Biol. Assoc. U.K. 47, 23–31. doi: 10.1017/S002531540003352X
Wright, M. S., Baker-Austin, C., Lindell, A. H., Stepanauskas, R., Stokes, H. W., and McArthur, J. V. (2008). Influence of industrial contamination on mobile genetic elements: class 1 integron abundance and gene cassette structure in aquatic bacterial communities. ISME J. 2, 417–428. doi: 10.1038/ismej.2008.8
Xavier, B. B., Das, A. J., Cochrane, G., De Ganck, S., Kumar-Singh, S., Aarestrup, F. M., et al. (2016). Consolidating and exploring antibiotic resistance gene data resources. J. Clin. Microbiol. 54, 851–859. doi: 10.1128/JCM.02717-15
Xiong, W., Sun, Y., Ding, X., Zhang, Y., and Zeng, Z. (2014). Antibiotic resistance genes occurrence and bacterial community composition in the Liuxi River. Front. Environ. Sci. 2:61. doi: 10.3389/fenvs.2014.00061
Xiong, W., Sun, Y., Ding, X., Wang, M., and Zeng, Z. (2015a). Selective pressure of antibiotics on ARGs and bacterial communities in manure-polluted freshwater-sediment microcosms. Front. Microbiol. 6:194. doi: 10.3389/fmicb.2015.00194
Xiong, W., Sun, Y., Zhang, T., Ding, X., Li, Y., Wang, M., et al. (2015b). Antibiotics, antibiotic resistance genes, and bacterial community composition in fresh water aquaculture environment in China. Microb. Ecol. 70, 425–432. doi: 10.1007/s00248-015-0583-x
Xu, J., Zhang, Y., Zhou, C., Guo, C., Wang, D., Du, P., et al. (2014). Distribution, sources and composition of antibiotics in sediment, overlying water and pore water from Taihu Lake. China. Sci. Total Environ. 49, 267–273. doi: 10.1016/j.scitotenv.2014.07.114
Yu, K., and Zhang, T. (2012). Metagenomic and metatranscriptomic analysis of microbial community structure and gene expression of activated sludge. PLoS One 7:e38183. doi: 10.1371/journal.pone.0038183
Zhang, H. P., Chen, R. H., Li, F. P., and Chen, L. (2015). Effect of flow rate on environmental variables and phytoplankton dynamics: results from field enclosures. Chin. J. Oceanol. Limnol. 33, 430–438. doi: 10.1007/s00343-015-4063-4
Zhang, X. X., Zhang, T., Zhang, M., Fang, H. H., and Cheng, S. P. (2009). Characterization and quantification of class 1 integrons and associated gene cassettes in sewage treatment plants. Appl. Microbiol. Biotechnol. 82, 1169–1177. doi: 10.1007/s00253-009-1886-y
Zhang, Y., Xu, J., Zhong, Z., Guo, C., Li, L., He, Y., et al. (2013). Degradation of sulfonamides antibiotics in lake water and sediment. Environ. Sci. Pollut. Res. Int. 20, 2372–2380. doi: 10.1007/s11356-012-1121-8
Zhao, Q., Wang, Y., Wang, S., Wang, Z., Du, X. D., Jiang, H., et al. (2016). Prevalence and abundance of florfenicol and linezolid resistance genes in soils adjacent to swine feedlots. Sci. Rep. 6:32192. doi: 10.1038/srep32192
Zhu, Y. G., Johnson, T. A., Su, J. Q., Qiao, M., Guo, G. X., Stedtfeld, R. D., et al. (2013). Diverse and abundant antibiotic resistance genes in Chinese swine farms. Proc. Natl. Acad. Sci. U.S.A. 110, 3435–3440. doi: 10.1073/pnas.1222743110
Zou, L., Meng, J., McDermott, P. F., Wang, F., Yang, Q., Cao, G., et al. (2014). Presence of disinfectant resistance genes in Escherichia coli isolated from retail meats in the USA. J. Antimicrob. Chemother. 69, 2644–2649. doi: 10.1093/jac/dku197
Keywords: antibiotic resistance genes, watersheds, metagenomics, high throughput screening, quantitative PCR, land-use, antibiotics
Citation: Uyaguari-Díaz MI, Croxen MA, Luo Z, Cronin KI, Chan M, Baticados WN, Nesbitt MJ, Li S, Miller KM, Dooley D, Hsiao W, Isaac-Renton JL, Tang P and Prystajecky N (2018) Human Activity Determines the Presence of Integron-Associated and Antibiotic Resistance Genes in Southwestern British Columbia. Front. Microbiol. 9:852. doi: 10.3389/fmicb.2018.00852
Received: 09 December 2017; Accepted: 13 April 2018;
Published: 01 May 2018.
Edited by:Karla B. Heidelberg, University of Southern California, United States
Reviewed by:Erica Marie Hartmann, Northwestern University, United States
Amelie Garenaux, University of California, Irvine, United States
Copyright © 2018 Uyaguari-Díaz, Croxen, Luo, Cronin, Chan, Baticados, Nesbitt, Li, Miller, Dooley, Hsiao, Isaac-Renton, Tang and Prystajecky. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.