Gabriela Flores-Vargas
Darren R. Korber1
Jordyn Bergsveinson- 1Food and Bioproduct Sciences, University of Saskatchewan, Saskatoon, SK, Canada
- 2Watershed Hydrology and Ecology Research Division, Environment and Climate Change Canada, Saskatoon, SK, Canada
The effects of sub-minimum inhibitory concentrations (sub-MICs) of antibiotics on aquatic environments is not yet fully understood. Here, we explore these effects by employing a replicated microcosm system fed with river water where biofilm communities were continuously exposed over an eight-week period to sub-MIC exposure (1/10, 1/50, and 1/100 MIC) to a mix of common antibiotics (ciprofloxacin, streptomycin, and oxytetracycline). Biofilms were examined using a structure–function approach entailing microscopy and metagenomic techniques, revealing details on the microbiome, resistome, virulome, and functional prediction. A comparison of three commonly used microbiome and resistome databases was also performed. Differences in biofilm architecture were observed between sub-MIC antibiotic treatments, with an overall reduction of extracellular polymeric substances and autotroph (algal and cyanobacteria) and protozoan biomass, particularly at the 1/10 sub-MIC condition. While metagenomic analyses demonstrated that microbial diversity was lowest at the sub-MIC 1/10 antibiotic treatment, resistome diversity was highest at sub-MIC 1/50. This study also notes the importance of benchmarking analysis tools and careful selection of reference databases, given the disparity in detected antimicrobial resistance genes (ARGs) identity and abundance across methods. Ultimately, the most detected ARGs in sub-MICs exposed biofilms were those that conferred resistance to aminoglycosides, tetracyclines, β-lactams, sulfonamides, and trimethoprim. Co-occurrence of microbiome and resistome features consistently showed a relationship between Proteobacteria genera and aminoglycoside ARGs. Our results support the hypothesis that constant exposure to sub-MICs antibiotics facilitate the transmission and promote prevalence of antibiotic resistance in riverine biofilms communities, and additionally shift overall microbial community metabolic function.
1. Introduction
Over the last decade, many studies have examined anthropogenic systems with high rates of antibiotic resistance events, colloquially-known as “Antimicrobial resistance (AMR) hotspots” (Karkman et al., 2018; Kraemer et al., 2019; Kunhikannan et al., 2021; Yadav and Kapley, 2021). Monitoring these hotspots (i.e., clinical settings, wastewater treatment plants, and pharmaceutical manufacturing plants) is critical to understand mechanisms of acquired resistance and the prevalence, maintenance, and transmission risk of pathogenic bacteria. The linkage between AMR hotspots and increased dissemination of antimicrobial resistance genes (ARGs) and bacteria (ARB) is clear (Munk et al., 2022); however, relatively little is known of AMR abundance and dynamics downstream of hotspots where sub-minimum inhibitory concentration (sub-MIC) of antibiotics may still exert selective pressure on environmental microbiomes (Ebmeyer et al., 2021).
Exposure to sub-MICs of antibiotics have been shown to promote the upregulation of mutagenesis and DNA transfer events which can lead organisms within a mixed community to acquire ARGs (Chow et al., 2021). There is evidence that sub-MICs of antibiotics induce the expression of genes coding for virulence factors (Laureti et al., 2013) which can proliferate via horizontal gene transfer (HGT) mechanisms. HGT events themselves, along with associated integron and transposon recombination, have been shown to be induced or mediated by the bacterial SOS response (Baharoglu et al., 2010). Indeed, sub-MICs levels of aminoglycosides, fluoroquinolones and β-lactams have reported to directly increase mutation and conjugation rates via the SOS response in Escherichia coli, Staphylococcus aureus, and Vibrio cholerae (Baharoglu et al., 2013; Andersson and Hughes, 2014).
To date, most studies examining sub-MICs effects on microorganism behavior have focused on single-species of pathogenic bacteria, or at best, the same family of bacteria (Davies et al., 2006; Laureti et al., 2013; Abe et al., 2020). Multispecies experimental data is scarce, but recently, studies have surveyed antibiotic concentrations in environmental microbiomes and their relation to AMR (Danner et al., 2019; Chow et al., 2021; Sanchez-Cid et al., 2022). In freshwater aquatic systems, which are the primary receiving environments of wastewater treatment plants (WWTP) discharge and agricultural run-off, biofilm communities are abundant. Biofilms are known to support the maintenance and dissemination of antibiotic resistance mechanisms, genes and a diversity of organisms (Abe et al., 2020). The physical nature of biofilm architecture and lifecycle confer considerable protection from surrounding chemical stressors because of the production of extracellular polymeric substances (EPS). This EPS matrix, and the close proximity of constituent microorganisms, allows for adaptive tolerance mechanisms to be supported as well as the inter-organism exchange of genetic material by various means (Flemming et al., 2022). Thus, naturally occurring aquatic biofilms could act as potential reservoirs and maintenance environments of AMR in a myriad of conditions (Balcázar et al., 2015; Guo et al., 2018; Flores-Vargas et al., 2021; Matviichuk et al., 2022). Studies examining biofilms response to high concentrations of antibiotics have often focused on microbial communities in close proximity to WWTP effluents (Aubertheau et al., 2017; Lépesová et al., 2018; Chonova et al., 2019; Maestre-Carballa et al., 2019; Petrovich et al., 2019; Matviichuk et al., 2022). However, a clear understanding of biofilm community response to constant, sub-MICs of antibiotics and their role in promoting or maintaining AMR in supposed “non-impacted” environments is still lacking (Cairns et al., 2018). Whole community metagenomic analysis provides for necessary data to understand the transfer and prevalence of AMR in the environment, which are fundamental dynamics aimed to perform future environmental risk assessments. To determine the potential risk that low antibiotic concentrations present in aquatic biofilm systems, the range of antibiotic concentrations that have an effect in shaping naturally occurring microbial communities must be defined, as must the response of the communities’ tendency to the disseminate ARGs and ARB. Herein we characterize the response of natural riverine biofilm communities under selective presence of various sub-MIC antibiotics using metagenomic sequencing and biofilm architecture analysis, to determine whether exposure influences both community structure and function, and the abundance of ARGs and AMR-related functions.
2. Materials and methods
2.1. Experimental design
Microcosm experiments for biofilm development were performed in rotating annular bioreactor (RAB) systems which have been previously described in detail (Lawrence et al., 2000, 2004). Natural river water from the South Saskatchewan River (Saskatoon, Canada) was collected weekly and used as inoculum to establish the microbial biofilms, and as a source of carbon and nutrients for biofilm growth and development. Physicochemical parameters of the South Saskatchewan River were previously described (Lawrence et al., 2007a); where water was collected within city limits, yet upstream of both the municipal WWTP facility and sources of agricultural discharges. River biofilms were provided a one-week establishment phase without antibiotics to initiate growth on the surface of the 12 removable polycarbonate strips (1 × 11 cm) within each RAB. The 12 RABs were run in parallel, and each treatment was replicated (n = 3). RAB dose treatments were provided by the direct and constant addition of an antibiotic cocktail of three commonly used antibiotics from different drug classes at sub-MIC dosages of 1/10, 1/50, and 1/100 their respective MIC level, resulting in μg/L levels within the RAB system. Antibiotics concentration and degradation within the RAB system were not assessed after their addition.
MIC values for the biofilm communities were determined based on breakpoint reports of EUCAST (European Committee on Antimicrobial Testing EUCAST) and CLSI (Clinical Laboratory Standards Institute; CLSI, 2018; EUCAST, 2019): ciprofloxacin (MIC 0.5 mg/L; 1/10: 50 μg/L, 1/50: 10 μg/L, 1/100: 5 μg/L), streptomycin (MIC 512 mg/L 1/10: 51,200 μg/L, 1/50: 10,240 μg/L, 1/100: 5120 μg/L), and oxytetracycline (MIC 125 mg/L 1/10: 12,500 μg/L, 1/50: 2,500 μg/L, 1/100: 1,250 μg/L) were used in this study (Supplementary Table S1). The selection of these three antibiotic drug classes was based on the list of antimicrobial susceptibility surveillance criteria by the WHO, where streptomycin and ciprofloxacin are listed as critically important and oxytetracycline is considered as highly important (WHO, 2017). Additionally, ciprofloxacin is an antibiotic routinely found in WWTP discharge; and oxytetracycline and streptomycin are broad-spectrum antibiotics widely used in Canadian and global livestock operations (Kraemer et al., 2019).
Control RABs were operated with river water alone. Biofilms were grown under treatment and control conditions for 8 weeks following a one-week establishment period, after which strips were removed from each RAB, immediately frozen and stored at −80°C for subsequent molecular analyses. All analyses (see below) were conducted on subsamples from randomly selected biofilm strips from each RAB replicate.
2.2. Microscopic analysis
2.2.1. Confocal laser scanning microscopy
Coupon pieces of 1 cm2 were excised from a randomly selected strip of each RAB and stained (below) prior to observation using a Confocal Laser Scanning Microscope Nikon Eclipse LV 110 DU and C2 camera with water-immersible lenses (10×, 40×, 60×; Nikon, Chiyoda, Tokyo, Japan). Biofilm architecture was quantified using a three-channel procedure (Neu et al., 2001; Dynes et al., 2006), where signals at green (excitation 488 nm, emission 522/32 nm), red (excitation 568 nm, emission 605/32 nm) and far-red channels (excitation 647 nm, emission 680/32 nm) were obtained. Coupons of 1 cm2 colonized by biofilm communities were directly stained with SYTO 9 (Molecular Probes, Eugene, OR, United States; excitation wavelength 488 nm, emission wavelength 522–532 nm) to detect nucleic acids of bacteria. EPS components of the biofilm matrix were visualized using three fluorescent fluor-conjugated lectin-binding dyes at 1 mg/mL: Triticum vulgaris-TRITC (TRITC: tetramethyl rhodamine isothiocyanate; excitation 568 nm, emission 605/32 nm; Sigma Chemicals, St. Louis, MI, United States) with polymer binding specificity for N-acetylglucosamine residues and oligomers; Arachis hypogaea-FITC (FITC: fluorescein isothiocyanate; excitation 485/495 nm, emission 510/600 nm; Sigma Chemicals, St. Louis, MI, United States) with polymer binding specificity for galactose and N-acetylglucosamine; and Canavalia ensiformis-FITC, also known as “Concanavalin A” (excitation 495/500 nm; emission 495/519 nm; Sigma Chemicals, St. Louis, MI, United States) with polymer binding specificity for mannose and glucose residues. Additionally, chlorophyl autofluorescence from algal and cyanobacteria cells (excitation 647 nm, emission 680/32 nm) was detected in the far-red channel (Lawrence et al., 2004; Dynes et al., 2006). Five areas per coupon were randomly selected for Z-stack image scanning with a slice interval of 5 μm (40×) and 1 μm (60×). CLSM image sequences were collected as previously described by Lawrence et al. (2005) and used for image analysis via ImageJ (Schneider et al., 2012) to define biofilm depth (or thickness) and architecture, and biomass of bacteria, EPS and autotrophs.
2.2.2. Protozoa enumeration
Enumeration of protozoan organisms were also considered for overall biofilm community structure. Weekly counts were performed in water-immersed 1 cm2 coupons from each RAB using a light microscope and protozoa or micrometazoa were identified based on morphology and manually counted (Packroff et al., 2002).
2.3. Metagenomic analysis
2.3.1. DNA extraction and whole genome sequencing
For each RAB, biofilm material from five randomly selected strips were recovered. Cellular biomass was aseptically obtained using a cell scraper (08-100-241; Fisher Scientific, Pittsburgh, PA) and the collected material centrifuged for 5 min at 9,000 × g to separate the water phase (Lawrence et al., 2009; Bergsveinson et al., 2020) and concentrate biofilm material to 300 mg wet weight. Total DNA was extracted using the Mag-Bind® Universal Pathogen Kit (Omega Bioservices, Norcross, GA, United States) and concentration and yield of DNA (5–10 ng/μL) was measured with the QuantiFluor dsDNA System on a Quantus Fluorometer (Promega, Madison, WI, United States). Whole genome libraries targeted for prokaryotic organisms were constructed using the KAPA Biosystems HyperPlus Kit following manufacturer’s instructions (Kapa Biosystems, Wilmington, MA, United States). Briefly, DNA was fragmented, ends were repaired, 3′ adenylated, and ligated to adapters. The resulting adapter-ligated libraries were PCR-amplified, Illumina indexes added, and pooled for multiplexed sequencing on an Illumina HiSeq4000/X10 platform (Illumina, San Diego, CA, United States) using the paired-end 150 bp run format, yielding 29 to 49 million reads per sample (average of 36 million reads; Supplementary Table S2).
2.3.2. Data normalization
Paired-end raw reads were quality filtered with Trimmomatic v0.36 (Bolger et al., 2014) using the following parameters: HEADCROP:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15. Adapters were removed using the TruSeq3 adapter sequence file as reference. All samples were subsampled to 25 million reads (Mreads) using the seqkt tool package.1 Briefly, FASTQ reads depth from each sample was rarefied to the lowest number of Mreads/Sample, remaining reads were randomly discarded. This method ensured a uniform population diversity to perform comparative analyses (Connelly et al., 2017; Chekabab et al., 2020). Trimmed and subsampled FASTQ reads were then analyzed by alignment to different reference databases, as described below.
2.3.3. Taxonomic profiling
Identification of bacteria and calculation of relative abundance at the genera level was assessed using three metagenome taxonomic profiling pipelines: CosmosID v2.0 (CosmosID Inc., Rockville, MD), Kraken v2 (Wood et al., 2019) and MetaPhlAn v3.0 (Segata et al., 2012). CosmosID is a web-based platform that utilizes data mining k-mer-based algorithms (Genius software) and Genbook, a high-performance curated comparator database, and includes references from NCBI-RefSeq, CARD, ARDB, VFDB, IMG and DDBJ and contains over 15,000 bacterial, 5,000 viral, 250 protists and 1,500 fungal species (Connelly et al., 2017; Junqueira et al., 2017; Yan et al., 2019; Chekabab et al., 2020). Taxonomic classification with the Kraken pipeline was performed with subsampled, trimmed FASTQ files using default parameters with the “—paired” option to indicate paired read files. Reads were mapped against the standard Kraken2 database comprised of complete genomes based on the NCBI RefSeq database (downloaded February 2022; O’Leary et al., 2016) for bacteria, viruses, fungi, protozoa and archaea. MetaPhlAn analysis was performed using the default parameters in addition to the “—ignore_eukaryotes” and “-t rel_ab_w_read_stats” options, the latter employed to estimate the number of reads per clade and obtain absolute abundance outputs for subsequent diversity indexes analysis. Microbial composition was compared against the CHOCOPhlAn_201901 database. Taxonomic abundance matrices from CosmosID, Kraken and MetaPhlAn analysis of each sample were prepared for analysis and visualization with all microbial taxa above 0.01% relative abundance included.
2.3.4. Resistome analysis
The resistome profiles of ARGs from all samples were identified and compared through three different databases: CosmosID, the Comprehensive Antibiotic Resistance Database (CARD, v.3.0.0; Jia et al., 2016), and ResFinder v. 4.1.11 (Bortolaia et al., 2020). Metagenomes were first assembled using MEGAHIT v. 1.2.9 (Li et al., 2015) and contigs annotated with prokka v.1.14.5 (Seemann, 2014; Supplementary Table S2). Quality-assembled contigs were mapped against CARD and ResFinder databases. The RGI (Resistance Gene Identifier) application v 6.0.2 in bwt mode for metagenomic reads was used for CARD analysis. Analysis using CosmosID was performed by mapping the trimmed and concatenated paired reads (.fastq) against the platform’s database using default settings and using the filtering threshold to ensure high confidence. Output of CosmosID included abundance counts table estimates of ARGs downloaded for further data analysis. The virulome of each sample was additionally profiled by comparing contig sequences against the Virulence Factor Database (VFDB; Liu et al., 2019) using default parameters and cut-off threshold of 90% identity and coverage. The ABRicate tool v.1.0.0 (Seemann, 2019) was used against the built-in ResFinder, CARD, and VFDB databases to corroborate the identity of all ARG/virulence factors identified.
2.3.5. Functional prediction
For functional prediction of the metabolic pathways influenced by sub-MIC antibiotics, the HUMAnN v.3.1.1 pipeline (Beghini et al., 2021)2 was used. Subsampled, trimmed paired reads were concatenated into a single file per sample and mapped against the ChocoPhlAn (microbial database) and the full UniRef (version: uniref90_201901b) reference databases for gene family identification. Default parameters were used with the “—taxonomic-profile” option, since the previously generated taxonomic profile output from MetaPhlAn was used as a custom taxonomic profile for the analysis. HUMAnN’s default curated metabolic pathway database, MetaCyc (Caspi et al., 2018) was used for pathway annotation. The abundance of each pathway (path_abundance.tsv file) of all samples was merged into a single table using the “humann_join_tables” option. The list of each identified pathway is represented in reads per kilobase (RPK) units with the abundance of each pathway categorized into bacterial organisms (in HUMAnN referred to as “stratifications”). After filtering out unidentified functional pathways, a pathway abundance table was used for subsequent multivariate and statistical analyses. To generate coverage information of calculated functional genes and ARG sequences, trimmed paired reads (.fastq) were mapped against the contigs of each metagenome assembly using Bowtie2 v. 2.4.4 (Langmead and Salzberg, 2012) with -bowtie2-build and default settings. Then, Samtools v.1.15.1 (Li et al., 2009) was used to visualize a file report of sequence coverage across samples.
2.4. Statistical analysis
All statistical analyses were performed in R version 4.2.0 (R Core Team, 2020). One-way ANOVA was used when normality assumptions were met (tested using Tukey HSD, for instance, for microscopic analyses). For non-parametric data, Kruskal-Wallis was applied to untransformed data to test for significant differences (p < 0.05) in biofilm composition, biomass abundance percentage, and diversity indices. p-values were adjusted using the Benjamini-Hochberg method to reduce false positive results (padj < 0.05; Benjamini and Hochberg, 1995). Diversity indices were calculated based on the untransformed abundance data of bacterial biofilm communities at genus and species levels. Species richness (Chao1), Shannon and Simpson indices were estimated within samples to describe α-diversity (Supplemental Figures S3, S4).
Comparisons between control and sub-MIC treatment-exposed communities were assessed through permutational analyses of variance based on Bray–Curtis dissimilarity indices (PERMANOVA, 999 permutations) calculated from Hellinger-transformed abundance count tables of the microbiome or resistome using vegdist function with vegan package v. 2.6.2 (Oksanen et al., 2020). Pairwise comparisons between control and treatments were performed to determine significance (p < 0.05). Bray-Curtis distance matrices were used to calculate β-diversity between samples. β-diversity was visualized with non-metric multidimensional scaling plots (nMDS) to ordinate microbiome and resistome data. Significant dissimilarity of ordination between groups was assessed using the Analysis of Similarities statistic (ANOSIM). To screen for significant differences (log2-fold changes) in gene and functional pathway abundance, the DESeq2 v.1.35 (Love et al., 2014) package was used, where the contrast function was employed to extract values per sub-MIC antibiotic condition.
The Bray-Curtis distance matrices resulting from the nMDS ordination were visualized for correlation patterns between the microbiome and resistome through the procrustes function in the vegan package. The Spearman’s rank correlation via Mantel test was used to compare the similarity of the microbiome and resistome databases, and to assess the relationship between the microbiome and the resistome in the presence of sub-MICs antibiotics, where correlations considered to be significant (p < 0.05) and with a dissimilarity coefficient > 0.75 were visualized in network analysis using Hmisc v.4.7-2 and igraph package v.1.3.5.
3. Results
3.1. Structural biofilm composition
Riverine biofilm structure was altered following sub-MICs antibiotic exposure. CLSM stacked images from microscopic analysis showed differences in biofilm architecture across treatments (Figure 1). Overall, sub-MICs antibiotic exposure decreased the complexity of biofilm communities, with the sub-MIC 1/10 treatment having the most pronounced effects on overall biofilm architecture, including biofilm thickness, and biomass percentage of bacteria, autotroph, and EPS composition. While pennate diatoms were observed in all samples, cyanobacteria (magenta fluorescence) were only recorded in control biofilms (Figure 1).
Figure 1. CLSM stacked image results of biofilms grown without or in the presence of sub-MICs antibiotic exposure after 8 weeks of development. Samples were stained by (A) Triticum vulgaris—TRITC lectin and SYTO9 showing EPS (red), bacteria (yellow), algae (blue) and cyanobacteria (magenta), (B) Arachis hypogaea—FITC and (C) Canavalia ensiformis—FITC lectins showing EPS (green), algae (blue), and cyanobacteria (magenta). Scale bar indicates 50 μm.
Biofilm thickness was calculated based on the fluorescence emitted by the lectin probes: C. ensiformis-FITC, A. hypogaea-FITC and T. vulgaris-TRITC (Figure 2) which target the EPS matrix associated with diverse microbial community members and are commonly used as indicators of changes within biofilm communities (Lawrence et al., 2007b, 2009). Decreased biofilm thickness was consistently observed with increasing strength of antibiotic treatment. The sub-MIC 1/10 treatment resulted in significantly decreased biofilm thickness relative to controls and the lowest (sub-MIC 1/100) treatment (F = 12.3, p < 0.05; Figure 2). sub-MIC 1/50 and sub-MIC 1/100 treatments resulted in thinner biofilms compared to control samples, though not significantly so.
Figure 2. Biofilm thickness of microbial communities under sub-MICs antibiotic exposure was based on fluorescence of lectin-binding specificity to EPS glycoconjugate residues. Mean values are displayed with biological (n = 3) and technical (n = 5) replicates (**p < 0.05). Edges represent quartile values.
The EPS characterized by the glycoconjugate distribution showed residues of mannose, glucose and N-acetylglucosamine through all biofilm samples, which are carbohydrate residues typically found in freshwater biofilms (Dynes et al., 2009; Flemming et al., 2022). The EPS biomass percentage, as measured by the fluorescent signal from the three lectin probes, significantly shifted (F = 8.8, p < 0.05) between samples (Figure 3B). Canavalia ensiformis-FITC and A. hypogaea-FITC signal percentage increased in the presence of sub-MICs antibiotics, whereas T. vulgaris-TRITC lectin binding to N-acetylglucosamine residues and oligomers decreased. Exopolymer composition throughout all samples was dominated by C. ensiformis-FITC binding to mannose and glucose residues, with an increase in fluorescence (binding specificity) after sub-MICs antibiotic exposure, and most notably observed in the sub-MIC 1/10 treatment.
Figure 3. Proportional abundance of biofilm architecture elements exposed to sub-MICs antibiotics. (A) Relative biomass abundance of bacteria, EPS and autotroph composition, and (B) Relative lectin-binding specificity of three probes used to characterize the glycoconjugate composition of the EPS matrix. Proportional abundance is the coverage percentage recorded of the fluorescent signal of each biofilm structure element from Z-stack images. Mean values are displayed with biological (n = 3) and technical (n = 5) replicates.
Proportional biomass abundance changes in the bacterial, autotroph and glycoconjugate composition (EPS matrix) were calculated based on analysis of the fluorescence of the three-channel images (Figure 3A). The autotroph proportional biomass amount (from signals of both prokaryotic and eukaryotic autotroph species) was significantly different between treatments (F = 4.9, p < 0.05) and was highest in the sub-MIC 1/100 treatment (36%), followed by the sub-MIC 1/50 treatment (33%), control samples (27%) and the sub-MIC 1/10 treatment (17%). Overall, EPS biomass significantly decreased (F = 12.2, p < 0.05) after sub-MICs antibiotic exposure relative to control values (63%), with the lowest biomass EPS composition detected in the sub-MIC 1/50 treatment (Figure 3A). Interestingly, CLSM image analysis did not show significant differences in the amount of bacterial biomass, although highest proportion was observed in the sub-MIC 1/10 treatment (25%), followed by the sub-MIC 1/50 treatment (17%), the sub-MIC 1/100 treatment (13%), and lastly the control samples (10%).
The cumulative number of protozoa was also tracked (Supplementary Figure S1), revealing that ciliates dominated control biofilms and overall protozoa counts significantly decreased (x2 = 7.8, p < 0.05) in the presence of sub-MICs antibiotics. However, no correlations, or more specific protozoan predation or grazing pattern interactions between biofilms developed under antibiotic concentrations (sub-MIC 1/10, 1/50 and 1/100), were detected. Despite the heterogeneous biofilm communities and thickness range, overall visualization of the EPS matrix confirmed that control biofilms contained thicker EPS than sub-MICs treated communities, yet conversely, bacterial biomass increased with sub-MICs antibiotic treatments. These results, in combination with decreased protozoan and cyanobacteria, indicate a selective effect of sub-MICs antibiotics on biofilm communities.
3.2. Microbiome of riverine biofilm communities
The bacterial community composition was analyzed by using three commonly utilized bioinformatic tools for metagenome taxonomic profiling: CosmosID, Kraken and MetaPhlAn. Across the three approaches, the taxonomic profile of the biofilm shifted in response to sub-MICs antibiotic treatment, with Actinobacteria and Proteobacteria being the predominant phyla groups, constituting more than 70% of overall community across all samples (Figure 4). Interestingly, the proportion of Proteobacteria was lowest in the control samples and increased along with sub-MIC antibiotics concentration. In contrast, Actinobacteria species were consistently reduced in the presence of sub-MIC antibiotics, as were abundances of Bacteroidetes, Cyanobacteria and Planctomycetes. The phylum Firmicutes was only detected via Kraken, with an average relative abundance of 30% and little variation across samples.
Figure 4. Comparison of taxonomic classification tools. Relative bacterial sequence abundances for (A) CosmosID, (B) Kraken, and (C) MetaPhlAn tools at the phylum level (n = 3).
While there were notable differences in the taxonomic profiles produced by the three different identification databases, most identified taxa were common across tools. A total of 11 different classes (Supplementary Figure S2) and 20 order groups (Figure 5) were identified with >1% relative abundance across samples. At the class level, bacterial abundance changes were mostly associated with Actinobacteria and Alphaproteobacteria, the former decreasing at the highest treatment concentration (sub-MIC 1/10). All tools showed similar abundance trends at the order level for Burkholderiales, Corynebacteriales, Sphingomondales, and Rhizobiales. The Bacilli class and Bacillales order were only observed in Kraken, constituting 25%–30% of relative abundance.
Figure 5. Comparison of taxonomic classification tools. Relative bacterial sequence abundances for (A) CosmosID, (B) Kraken, and (C) MetaPhlAn tools at the order level (>1%; n = 3).
At the level of genus and species, profile differences were observed across the identification tools. CosmosID analysis detected 84 different genera, and 100 bacterial species; MetaPhlAn identified 31 genera and 39 species; and Kraken detected 251 genera and 2,572 species (Figure 6B). Comparison of these tools showed that MetaPhlAn was less sensitive for identifying bacteria at species level (most groups were categorized as unclassified), whereas Kraken showed high sensitivity to species strains. α-diversity showed markedly different trends according to the taxonomic pipeline used. Genus-level richness was determined to be higher in control biofilms relative to sub-MIC samples, except for MetaPhlAn results; however, significant differences were only observed between the 1/100 antibiotic treatment and control samples (Kruskal-Wallis x2 = 8.4, p < 0.05, Kraken, Supplementary Figure S4). Species level-richness was consistently lowest in the sub-MIC 1/100 treatment, although the analysis did not demonstrate significant differences across samples (Supplementary Figure S5).
Figure 6. Bacterial composition at the genus level according to three taxonomic databases of riverine biofilm communities cultivated under sub-MICs antibiotic exposure. (A) nMDS two-dimensional plot of Bray-Curtis similarities showing β-diversity identified in the three taxonomic databases: CosmosID (84 = genera), Kraken (251 = genera) and MetaPhlAn (31 = genera). (B) Venn diagram displaying identified total genera groups across taxonomic databases. Overlap indicates the number of shared taxa among databases. (C) nMDS plot of Bray-Curtis similarities of significant genera (p < 0.001) detected at each taxonomic database. Ellipses shape is defined by covariance of each group, and ellipses centroid represents the group mean.
Despite the observed differences, Shannon and Simpson indices calculated from all databases demonstrated higher genera diversity of biofilms at the sub-MIC 1/50 treatment compared to the sub-MIC 1/10 condition. Species-level diversity yielded consistently similar diversity distributions, but just for the Shannon index. Overall, calculated genus and species-level α-diversities from CosmosID and MetaPhlAn were more similar to each other than to the Kraken results. For example, diversity of control biofilms relative to sub-MIC treatments were less diverse in MetaPhlAn and CosmosID analyses, yet control samples in the Kraken analysis, the most species-sensitive package, were the most taxonomically diverse, with significant differences when compared to sub-MIC 1/10 treated communities (x2 = 7.8, p < 0.05; Supplementary Figure S4).
The Bray-Curtis distance matrices of identified bacteria at genus level were used to generate nMDS plots to assess the dissimilarities between biofilm microbial communities developed under different sub-MIC antibiotic exposures (Figure 6A). The three taxonomic tools displayed significant differences across sub-MIC antibiotic treated samples (CosmosID: R2 = 0.577, p = 0.017; Kraken: R2 = 0.609, p = 0.040; MetaPhlAn: R2 = 0.537, p = 0.018). Of the total 296 unique genus groups identified across taxonomic analyses, only 20 were shared across the three databases (Figure 6B). Interestingly, seven genera were shown to be significant across all three taxonomic tools (PERMANOVA, F = 3.6, p < 0.001, Bray-Curtis distance), suggesting their abundances drove the ordination pattern observed among biofilm communities (Figure 6C) Mycobacterium and Nitrospira genera were associated with control samples, whereas Flavobacterium, Methylorubrum and Aminobacter genera were associated with biofilm communities cultivated under sub-MIC 1/10 and sub-MIC 1/50 conditions.
Correlation analysis between the taxonomic databases indicates that identified genera from the CosmosID and MetaPhlAn databases shared the most similarities in taxonomic composition when examining genus-level abundance (Spearman’s rho = 0.895; Supplementary Figure S7A). Comparison of the genera abundance from the three databases was significant (p < 0.05), but Mantel results from CosmosID vs. Kraken (Spearman’s rho = 0.711) and MetaPhlAn vs. Kraken (Spearman’s rho = 0.644) were only moderately correlated.
3.3. Biofilm community resistome and virulome trends
Three different antibiotic resistance databases were compared to determine the resistome of riverine biofilm communities under sub-MIC antibiotic exposure. Across all biofilm community samples, a total of 92 ARGs were detected via CosmosID, 93 ARGs in CARD, and 67 ARGs in the ResFinder database (Figure 7B). Normalized abundance of ARGs determined from the three database analyses demonstrated a marked increase in the copy numbers of ARG in samples under sub-MIC exposure relative to control samples (Supplementary Figure S3).
Figure 7. ARGs composition according to the resistome database in riverine biofilm communities under sub-MICs antibiotic exposure. (A) Relative abundance classified by drug class resistance identified in each database. (B) Venn diagram displaying total identified number of ARGs. Overlap indicates the number of shared ARGs across databases. (C) ARGs with significant fold-changes (p < 0.05) between treatments. The normalized abundance heatmap is scaled to each gene (Z-score) after rlog-transforming counts.
Resistome composition from the CosmosID pipeline consisted of 33 different types of ARGs corresponding to seven drug resistance classes: aminoglycosides, β-lactams, macrolides, phenicols, sulfonamides, tetracyclines and trimethoprim. ARGs from the aminoglycoside drug class represented the majority of the resistome, followed by tetracycline and β-lactam genes (Figure 7C). Results from the CARD database detected 24 different types of ARGs which were grouped into nine drug resistance classes: aminoglycosides, β-lactams, cyclic peptides, macrolides, penams, phenicols, sulfonamides, tetracyclines and trimethoprim. CosmosID, relative abundances showed ARGs from aminoglycoside drug class to be the predominant group across control and sub-MICs exposed samples; however, macrolide and sulphonamide ARGs were the next most-abundant genes (Figure 7A). The ResFinder database identified 25 different types of ARGs across samples (Figure 7A) within eight different drug classes: aminoglycosides, β-lactams, macrolides, phenicols, rifamycins, sulfonamides, tetracyclines, and trimethoprim. However, no ARGs in the ResFinder database with an identity above the 90% identity cut-off threshold were detected in any of the control biological replicates.
Several ARGs were detected in at least two of the resistome databases (Figure 7B), with all three pipelines predicting the same seven ARGs under sub-MIC conditions (Supplementary Table S3). Aminoglycoside resistance was present at greater than 70% relative abundance through all samples according to CosmosID and CARD databases and were the predominant ARGs group in the ResFinder analysis. The aph6-Id aminoglycoside gene was prevalent throughout all samples and across all three pipelines. When examining ARGs that were more abundant in sub-MIC exposures relative to control samples (log2FoldChange, padj < 0.05), all three gene resistance databases predicted 18 common ARGs consisting of the classes aminoglycosides, β-lactams, macrolides, sulphonamides, tetracyclines and trimethoprim (Figure 7C). Most notably, gene copies of aadA5, drfA1, ereA, and ermF were highest at the sub-MIC 1/50 treatment, whereas sul2 was highest at the sub-MIC 1/100 treatment.
We further measured the richness of the resistome in all samples classified by gene, drug class, and mechanism of resistance (Supplementary Figure S6). Across the three profiled pipelines, microbial communities developed under sub-MICs antibiotic exposure displayed higher absolute numbers of ARGs with gene richness significantly greater than control samples (H=, p < 0.05). Annotations performed by the CosmosID pipeline showed that the seven drug classes and four mechanisms of antibiotic resistance were significantly increased in sub-MIC treatments relative to controls. The CARD pipeline revealed significant differences between sub-MIC 1/10 and 1/100 treatments relative to control samples only at the mechanism level, but not at the gene or drug class level. Results from the ResFinder pipeline showed significant differences between resistome richness of the sub-MIC 1/50 treatment relative to control samples at the level of individual ARG, drug class, and mechanisms of resistance (Supplementary Figure S6). ARGs from the CARD database displayed a significant correlation with the ResFinder database (Spearman’s rho = 0.556, p = 0.004). However, the relationships between the other databases showed non-significant and weak associations with identified ARGs (Supplementary Figure S7B).
The virulome of riverine biofilm communities was also analyzed, since it is known that pathogenic or non-pathogenic bacteria can acquire “accidental virulence” as result of the advantageous conditions within biofilm niches (Andersson et al., 2018). Virulence genes associated with conjugative transfer, efflux pumps, DNA-transfer primases, integrons, mating-pair-formation, resolvases, transposons, and multidrug factors, were detected across the three sub-MIC antibiotic conditions (Figure 8; Supplementary Figure S8). Virulome composition was shown to be significantly different across sub-MIC conditions (x2 = 8.9, p < 0.05). A total of 89 genes associated with virulence were identified across all samples, among which transposon genes (i.e., TpnA) and multidrug resistance genes (i.e., sul1 and qacEdelta1) were not only observed in control samples but were also the most abundant type of virulence factors according to the VFDB database (Figure 8). Samples from sub-MIC 1/10 and 1/50 treatments showed very similar virulome compositions. However, extended-spectrum β-lactamase (ESBL) genes, such as blaOXA-1 and blaOXA-2, were associated with sub-MIC 1/50 and 1/100 antibiotic samples. Moreover, ssb, a gene coding for single-stranded DNA-binding protein, was only detected in one control replicate.
Figure 8. Differential abundance of normalized virulence factor comparisons with significant differences between sub-MICs antibiotic treatments (p < 0.05) using Benjamini-Hochberg correction for multiple comparisons.
3.4. Biofilm community functional capacity
The prediction of metabolic pathways in riverine biofilms revealed a total of 106 functional pathways present throughout samples. “Unidentified” pathways and those that confer functionality to non-bacterial organisms (such as plants or mammals) were manually removed prior to statistical analysis and visualization. Differential gene abundance changes of 99 pathways were analyzed via DESeq2, yielding 28 significant functional pathways related to sub-MIC antibiotic treatments relative to control samples (Figure 9). In general, detected bacterial metabolic pathways were related to biosynthesis and nucleotide salvage. Following the nomenclature by the MetaCyc database, these biosynthesis pathways included the super-classes for nucleoside and nucleotide biosynthesis (PWY-5686, PWY-6126, PWY-6125, PWY-6277, PWY-6609), phosphate-related pathways (NONOXIPENT-PWY, PENTOSE-P-PWY), cell wall biosynthesis (PWY-6385, PEPTIDOGLYCANSYN-PWY), fatty acid and lipid biosynthesis (PHOSLIPSYN-PWY, FASYN-ELONG-PWY, PWY-6282), secondary metabolite biosynthesis (NONMEVIPP-PWY, PWY-6703), Cofactor, carrier and vitamin biosynthesis (PANTO-PWY, COBALSYN-PWY, 1CMET2-PWY, COA-PWY, HEME-II), and aromatic compound biosynthesis (ARO-PWY).
Figure 9. Functional profiling of riverine biofilm communities under sub-MICs antibiotic exposure. Y-axis represents metabolic pathways annotated by MetCyc from HUMAnN analysis; and X-axis represents significant differential abundance sub-MICs antibiotic condition relative to control samples using DESeq2 (p < 0.05).
Among genes associated to superclasses for amino-acid biosynthesis, L-isoleucine (PWY-3001) and L-arginine (PWY-5154) were more abundant under sub-MIC antibiotic exposure compared to control samples; yet L-histidine (HISTSYN-PWY), L-arginine (ARGSYN-PWY, ARGSYNBSUB-PWY) and L-lysine (DAPLYSINESYN-PWY) displayed modest decreases in differential gene abundance (log2 Fold-change < −1., p < 0.05) in sub-MIC antibiotic conditions, except for sub-MIC 1/100 (Figure 9). Interestingly, genes related to arginine biosynthesis via acetyl-L-citrulline showed two log2 Fold-change higher abundance in the three sub-MIC antibiotic treatments; however, genes related to arginine biosynthesis via the acetyl cycle and the L-ornithine pathway were less abundant (−0.5 log2 Fold-change) in the sub-MIC 1/10 and sub-MIC 1/50 treatments, respectively. Overall, genes of functional pathways involved in biosynthesis showed negative differential abundances under exposure to sub-MIC treatments, with abundance fold-change varying across treatment. Several genes involved in functional pathways showed significant relative abundance differences (p < 0.05) between sub-MIC 1/10 (23 pathways) and sub-MIC 1/50 (22 pathways) treatments relative to controls, although the sub-MIC 1/100 treated samples showed fewer significantly differentiated metabolic pathways (13 pathways).
Nucleoside and nucleotide biosynthesis pathways demonstrated lower relative gene abundances (log2 Fold-change −2.40 to −0.23., p < 0.05) at all sub-MIC treatments in comparison to control samples. In contrast, two cell wall biosynthesis pathways showed higher abundances across the three sub-MIC antibiotic treatments (log2 Fold change1.9–2.2, p < 0.05) relative to controls. Notable degradation pathways were assigned to the general fatty-acid pathway (FAO-PWY), which was present in sub-MIC 1/10 and 1/50 treatments, and the urea cycle (PWY-4984). The urea cycle and fatty acid (PWY-6282, FASYN-ELONG-PWY) pathways had positive abundance shifts for sub-MIC 1/10 and sub-MIC 1/50 treatments but negative shifts for the sub-MIC 1/100 treatment relative to control biofilms. Notably, urea degradation is commonly associated to Nitrospira spp.; however, this genus was most strongly associated with the control biofilm condition, suggesting nitrogen and urea cycling are not clearly impacted or linked to sub-MIC treatments. Functional pathways with the most pronounced abundance shift between treatments were illustrated according to the contributing bacterial genera (Supplementary Figure S9). The bacterial species associated with detected functional pathways included the genera Afipia, Aminobacter, Blastomonas, Bosea, Bradyrhizobium, Caulobacter, Flavobacterium, Sphingomonas, Methylibium, Methylorubrum, Mycobacterium, Mycobacteroides, and Variovorax.
3.5. Microbiome and resistome relationship
The co-occurrence relationship between the microbiome and the resistome composition was assayed through the Procrustes analysis by using the Mantel test and Procrustes test for statistical validation (Supplementary Figure S10). Both tests showed significant correlation between the microbiome and resistome dissimilarity matrices (Mantel r = 0.652, p = 0.001; Procrustes M2 = 0.257, r = 0.861, p = 0.001), thus demonstrating that the bacterial and resistome composition were closely associated. Since the Procrustes analysis only illustrates overall correlations, a network analysis was performed to further examine individual co-occurrence patterns between the microbiome and resistome in biofilm communities after sub-MIC antibiotic exposure (Figure 10). The network analysis was based on strong and significant correlations (Spearman’s rho > 0.75, p < 0.05) between a correlation matrix constructed with 23 ARGs and 20 bacterial taxa (genus level).
Figure 10. Network analysis showing correlation patterns between the microbiome and resistome after sub-MICs antibiotics exposure. Edges (lines) connecting nodes mean a strong correlation (Spearman’s rho > 0.75, p < 0.05), negative indicated in red and positive in gray. The size of each node is proportional to the number of connections (degree; n = 12).
Additionally, co-occurrence analysis indicating resistome and microbiome patterns for each of the sub-MIC treatments (1/10, 1/50, and 1/100) was performed (Supplementary Figure S11). For details about the networks’ topological properties, such as the number of edges between nodes of each network, see Supplementary Tables S5, S6. Based on the connections between nodes, the most densely connected nodes were regarded as indicators of co-occurrence of ARGs. Microbiome composition results showed that Aminobacter, Bosea, Flavobacterium, Sphingomonas and Variovorax displayed the highest co-occurrence (4 nodes) to microbial and resistome nodes. Aminoglycoside ARGs (aph3”-Ib, aadA2, rrsB) demonstrated co-occurrence to Proteobacteria genera, but also to Actinobacteria (aadA, rrsA, rrsB), Bacteroidetes (aph3”-Ib) and Nitrospirota (rrsB; Figure 10). Indeed, Proteobacteria genera, such as Bradyrhizobium, Methylorubrum, Novosphingobium, Sphingobium, and Variovorax, correlated with aph3”-Ib, aph6-Id, aadA-, which were prevalent genes. The resistome composition displayed the highest co-occurrence to microbiome nodes for the sub-MIC 1/50 treatment with 23 ARGs, followed by the sub-MIC 1/100 treatment with 11 ARGs, and the sub-MIC 1/10 treatment with six ARGs. Aminoglycoside ARGs (aph6-Id, rrsA, rrsB) were present in all sub-MIC-exposed samples and correlated with Proteobacteria, as well as ARGs from other drug classes, except in control biofilms where these ARGs correlated only within the resistome composition (Supplementary Figure S11). The ARGs from the resistome nodes with the highest correlation to microbial nodes were aph3”-Ia and blaOXA-2. In the sub-MIC 1/10 treatment condition, Aminobacter, Bosea, Caulobacter and Sphingorhabdus genera showed a positive correlation to sul1.
This pattern indicates that Proteobacteria such as Caulobacter, Bosea, Novosphingobium, Hydrogenophaga, Methylibium, and more specifically Variovorax may act as potential reservoirs for aminoglycoside, β-lactam and sulfonamide resistance factors.
4. Discussion
4.1. Sub-MICs antibiotic affects the structure of riverine biofilm communities
Biofilm communities favor the maintenance and transmission of antibiotic resistance due to their diverse composition, close cell proximity, structural stability, and physical protection against chemical and physical stressors (Balcázar et al., 2015; Abe et al., 2020; Flores-Vargas et al., 2021; Matviichuk et al., 2022). In this study, heterogeneity between and within communities was observed, with up to a 54-μm range in biofilm thickness observed in samples from the same treatment (Figure 2). However, the effects of sub-MIC antibiotic treatments were most notable on biofilm communities developed under the sub-MIC 1/10 condition, where significantly thinner biofilms compared to other treatments (p < 0.05) were observed. Thinner biofilms have previously been noted to form in the presence of chemical stressors (Romero et al., 2018); therefore, thicker biofilms would seemingly proliferate under normal (non-stressed) growing-conditions. Thicker biofilms imply more diffusional distance within the matrix biomass for chemicals (i.e., antibiotics) to react with, thus initially providing more time (and reduced effective compound concentration) for microbial cells within the biofilm to adapt to stressor conditions. Furthermore, thicker biofilms imply greater protection provided by the volume of the EPS matrix, which is associated to increased cell density and closer proximity of cells (physical juxtapositioning), and consequently higher likelihood of cell–cell communication (quorum sensing) and genetic exchange which would facilitate HGT events (Dynes et al., 2009; Abe et al., 2020).
Carbohydrates are the structural skeleton of EPS matrix and can range from oligosaccharides through to polymers of up to 100 kD in size (Hanlon et al., 2006). Although polysaccharide composition in EPS of biofilms is highly heterogenous, most consist of a N-acetylglucosamine core (Flemming et al., 2022). The role of N-acetylglucosamine in EPS has mostly been studied in clinical monoculture biofilms such as Pseudomonas aeruginosa (Reichhardt et al., 2020), and is also known for being one of the major components of cell wall in bacteria (Gonzalez-Martinez et al., 2014). Result from our CLSM approach targeting glycoconjugate residues show a clear reduction of EPS polysaccharides, which could translate into lower signaling and hindered genetic exchange (Zhang et al., 2015). Since among biofilm elements, cyanobacteria and green algae exude large quantities of polysaccharides and amino-acids among other organic compounds (Battin et al., 2016; Delattre et al., 2016), it is reasonable to assume that the decrease of autotrophic organisms in biofilms treated with sub-MIC antibiotics, particularly at sub-MIC 1/10, had repercussions involving EPS accumulation (Figures 1, 3).
Riverine biofilms play a key role in river biogeochemical organic and inorganic cycles (Romaní et al., 2014) and EPS can be used as carbon source and energy reserves for bacterial heterotroph community and macrofauna (Pierre et al., 2012; Delattre et al., 2016), thus EPS reduction could hinder the further development of aquatic microorganisms. Additionally, EPS accumulation is crucial for adhesion to substrates (Fulaz et al., 2019), meaning that a reduction of EPS could contribute to earlier and/or higher frequency of detachment of biofilms, and further promote the spread of ARB and elements containing ARGs to and throughout the surrounding environment.
Biofilm formation in response to sub-MIC antibiotics remains contradictory even in single, clinically relevant bacterial isolates. For example, a commonly cited study demonstrates sub-MICs of aminoglycosides increased biofilm formation in P. aeruginosa and E. coli (Hoffman et al., 2005). However, a reduction in biofilm formation was reported for sub-MIC of β-lactams in Staphylococcus sp. (Cerca et al., 2005), macrolides in Mycobacterium avium (Carter et al., 2004), and mupirocin in P. aeruginosa (Horii et al., 2003). A more recent study documented how the presence of several antibiotics at sub-MIC levels affected S. aureus biofilm formation (Majidpour et al., 2017), where azithromycin and vancomycin increased biofilm formation, but linezolid, cefazolin, and clarithromycin inhibited it. A similar response was observed in Staphylococcus epidermidis biofilms where 0.5 MIC ciprofloxacin reduced biofilm thickness, but 0.5 and 0.25 MIC tigecycline increased biofilm growth, as well as the expression of icaA, altE, and sigB genes involved in polysaccharide intercellular adhesion (Szczuka et al., 2017). These studies are consistent with our microbiome-level results, given that sub-MIC antibiotic cocktail mix treatments containing ciprofloxacin in range between 5 to 50 μg/L decreased biofilm structural complexity. These findings demonstrate that effect on biofilm formation and its constituent members is likely antibiotic-dependent; thus, antibiotics affect diverse bacteria within biofilms differently. Indeed, the bioavailability and chemical state of antibiotics would also be presumed to influence the severity of response and/or impact on the biofilm, though these parameters are not within the scope of the current study.
In addition to chemical and physical stressors, biofilms are also exposed to biological stressors that include nutrient availability, protozoan predation and viral lysis (Andersson et al., 2018). Different protozoa have the potential to influence the proportion of bacterial species via selective or non-selective grazing, and have been shown to promote the occurrence of HGT events during short time frames (Sun et al., 2018). Protozoa were first detected in the fourth week of treatment exposure, with counts increasing slightly over the following weeks of the experiment, in alignment with the expected biofilm’s structural maturation throughout the eight-week period. As the biofilm matures, increased thickness and biomass would be expected to result in increased grazing pressure by protozoan organisms (Lawrence et al., 2009). Protozoa counts decreased in the presence of sub-MIC antibiotics, though we could not correlate the sub-MIC treatment to the protozoa composition or predation-interactions. It is possible that the observed reduction of protozoa organisms in biofilms under sub-MIC exposure was due to the limited availability of autotrophic bacteria recorded under sub-MIC antibiotic conditions, indicating an indirect effect through trophic interactions (Figure 1; Supplementary Figure S1). It has been previously reported that high protozoan predation on bacteria was associated with persistence of Mycobacterium and Rickettsia spp. (Andersson et al., 2018). Mycobacterium spp. were prevalent taxa throughout the biofilm communities of this study; thus, it is possible that predatory protozoa in these samples would be associated with, and contribute to the maintenance of, Mycobacterium spp.
Ultimately, reduced EPS results in an increase in the proportion of unprotected bacterial biomass but a reduction in protozoan and cyanobacteria (autotrophic bacteria). Loss of diversity and decreased protection from the EPS matrix signify a broader ecological impact of sub-MIC antibiotic exposure that could further accelerate the acquisition or expression of resistance functions of bacteria species within biofilm communities.
4.2. Microbiome composition shifts due to sub-MIC antibiotic exposure
Microbial composition analysis revealed that Proteobacteria species were favored by sub-MICs antibiotic presence, and overall, a considerable shift in bacterial community composition was observed at sub-MIC 1/10. Relative abundance profiles in this study (Figure 4) were in accordance with previous taxonomic results noted for South Saskatchewan River biofilms by Lawrence et al. (2020), where the dominant phyla (Proteobacteria, Actinobacteria and Cyanobacteria) comprised 63% of total bacterial community membership. Proteobacteria, particularly Betaproteobacteria, are implicated in more readily attaching to surfaces during initial biofilm development (Luo et al., 2020). Proteobacteria predominance throughout treated biofilm communities also signals their adaptability to low or sub-MICs antibiotics (Manaia, 2017) supporting the idea that species within this phylum possess antibiotic resistance mechanisms or potential for acquired resistance.
Additionally, river biofilms influenced by nearby run-off from WWTPs also showed Proteobacteria to be the dominant phyla (Guo et al., 2020; Matviichuk et al., 2022). Our results are in agreement with Zou et al. (2018) where oxytetracycline residues selected for Proteobacteria abundance. In contrast, Subirats et al. (2018) did not detect significant variations in bacterial abundance of biofilms treated with ciprofloxacin, erythromycin, sulfamethoxazole, diclofenac, and methylparaben, even though these antibiotics are known to inhibit bacterial growth at sub-MIC levels. Notably, that experimental design comprised 28 days, instead of the 56-day period used for this study. Thus, it is possible in Subirats et al. that the EPS matrix provided protection during early sub-MICs antibiotic exposure. However, the prolonged period of antibiotic exposure in our study might indicate the scale at which ecological impacts may be observed due to resultant shifts in microbial community composition.
nMDS analysis demonstrated that sub-MIC 1/10 and sub-MIC 1/50 treatments most profoundly shaped microbiome composition at the phyla and genus level. For example, Brevundimonas and Methylorubrum (Proteobacteria) as well as Flavobacterium (Bacteoridetes) relative abundance displayed significant correlation to biofilms grown in sub-MIC 1/50 antibiotic treatment; interestingly the two latter genera also demonstrated co-occurrence to Variovorax in the network analysis (Figure 10). Nitrospira (Nitrospirota phylum) revealed to be strongly associated with control biofilms (Figure 6) and showed co-occurred with several ARGs and genera in sub-MIC 1/50 treated biofilms. The prevalence of Mycobacterium spp. despite reduced abundance (up to 50%) in the presence of sub-MIC antibiotics (relative to control biofilms) presumably indicates that Actinobacteria were sensitive to the effects of antibiotics at these concentrations, or that antibiotics hindered the development of these genera.
Calculation and interpretation of α-diversity and species richness proved to be very dependent on diversity indices, taxonomic level, and databases (Supplementary Figures S4, S5). Bacterial α-diversity data of biofilms in the South Saskatchewan River is limited, however results from CosmosID and Kraken databases demonstrated that biofilms in this study were similar in diversity and community distribution to previous microcosm characterizations (Romaní et al., 2014). Furthermore, diversity distribution was similar between the sub-MIC 1/10 and sub-MIC 1/50 samples, which were generally less diverse than controls. This is not unexpected as it has been previously recorded that bacterial richness, diversity and composition are disturbed by the presence of antibiotics such as ciprofloxacin (Gonzalez-Martinez et al., 2014) and oxytetracycline (Zhang et al., 2013).
In line with this observation, it is frequently concluded that a decrease in diversity is an indicator of altered microbial taxonomy and by extension, functional environmental health functions and services. However, this data cautions against using α- or β-diversity as a sole metric or predictor of ecological impact or health in face of environmental stresses, given the dependence on database health, data normalizations and granularity of taxa level (Shade, 2017).
4.3. Comparisons of microbial reference databases
To our knowledge, this is the first study to evaluate different commonly used resistome and taxonomic databases for environmental aquatic microbiomes; and thus can serve to provide guidance for future study considerations and best practices. CosmosID, Kraken and MetaPhlAn are commonly-used analysis pipelines that were all separately applied in this study to profile taxa of biofilm communities (Couto et al., 2018). These analysis tools and pipelines notably influenced observations, though consensus data of bacterial abundance was achieved. These differences in data outputs are a result of the underlying databases used for alignment. Specifically, Kraken uses genomes based on the NCBI RefSeq database consisting of over 128,299 organisms (O’Leary et al., 2016); MetaPhlAn compares sequences against the CHOCOPhlAn_201901 database v.4.0 (Segata et al., 2012), containing over 771,500 metagenomic assembled genomes; and CosmosID aligns against the GenBook database (150,000 microbial genomes and 5,000 viral species) using k-mer-based alignment (Yan et al., 2019; Zaouri et al., 2020). In the present study, MetaPhlAn was less able to identify sequences to the species level (20 of 59 species were categorized as unclassified), whereas Kraken assigned the most sequences to species level.
Similar to taxonomic analysis, three different commonly used resistome databases were evaluated for sensitivity and accuracy of ARG annotation of metagenomes: CosmosID (5,500 ARGs), CARD (1,600 manually-curated ARG sequences) and ResFinder (3,000 ARG). Again, the selection of database directly affected final observations; however, all three analysis streams displayed the trend of sub-MIC antibiotics selecting for an increase in ARG abundance within biofilms. Correlation analysis revealed that CARD and ResFinder databases were most similar in their output and performance (Supplementary Figure S7B).
Despite the known challenges of comparing metagenomic data between studies (Couto et al., 2018), the dependence of data interpretation and conclusions on database quality, environmental metagenomics studies frequently present data produced from a favored tool or database. By performing analysis using one method and failing to compare the output against other available tools, particularly in the field of environmental AMR, researchers’ risk either under- or over-interpreting their results. Here we conclude that the CARD database and associated tools are currently the most reliable resistome database due to its high level of curation and maintenance, and that the CosmosID database is a reliable source for taxonomic profiling of aquatic environmental microbiomes, especially if the study seeks to identify Gram-negative or pathogenic bacteria within samples (Couto et al., 2018).
4.4. Resistome composition enhanced after sub-MIC antibiotic exposure
Resistome profiling of riverine biofilms developed under a sub-MICs antibiotic cocktail mix of ciprofloxacin, oxytetracycline and streptomycin revealed an abundance prevalence and dominance for aminoglycoside ARGs (Supplementary Figure S3). Resistome distribution classified by drug resistance and mechanism was consistently similar across the three resistome pipelines, with sub-MICs antibiotic treated biofilms yielding greater number of resistome elements in comparison to control biofilms (Supplementary Figure S6). These findings are consistent with resistome abundance observed in biofilms influenced by hospital wastewater (Petrovich et al., 2020) and WWTP discharge (Matviichuk et al., 2022). Microbiome composition aligned with resistome results in that the highest richness of taxa/genes were registered for the sub-MIC 1/50 treatment in comparison to control biofilms. Overall, aph6-Id, aph3”-Ib, aadA2 tetC, tetG, sul1, and sul2 were highly abundant and/or prevalent across all biofilms under sub-MICs exposure (Supplementary Table S3).
It has been previously reported that aminoglycoside resistance genes co-select with other ARGs such as those encoding β-lactamases and are often located within different mobile genetic elements (MGE) such as integrons and plasmids (Beceiro et al., 2013). Since ARGs are often co-located with pathogen virulence factors in biofilm organisms, we also screened for the occurrence of virulence genes. Virulence refers not only the ability of bacteria to cause disease in the host (i.e., degree of pathogenicity) but also the ability of an organism to infiltrate and colonize a host (Schroeder et al., 2017). TpnA transposase and sul1 sulfonamide genes were notably prevalent across all treatments including control biofilms according to the VFDB database (Supplementary Figure S8). TpnA has been reported to commonly occur in drinking water systems (Brumfield et al., 2020), whereas sul1 and sul2 have been recorded to be prevalent ARGs in influent and effluent riverine biofilms (Auguet et al., 2017).
Interestingly, relative abundance of intI1, a class I integron gene usually located in MGE such as transposons and plasmids (Gillings et al., 2015) and linked to several resistance genes including sul1 (Subirats et al., 2018; Cheng et al., 2020), was enhanced in all sub-MICs antibiotic treated samples. The sul1 ARG is of significant clinical relevance (Petrovich et al., 2020) and its abundance has been found to be greater in extracellular DNA than intracellular DNA in estuarine biofilms (Guo et al., 2018), supporting the idea that mobilization of this gene plays an important role in resistance transmission. The prevalence of intI1, sul1 and tpnA genes across biofilms under sub-MICs antibiotic exposure might be result of their mobilization across bacteria, particularly Proteobacteria species, as HGT dissemination events are more common between phylogenetically related bacteria (Subirats et al., 2018). Given this increase in suspected HGT events and the fact that constant sub-MICs exposure of antibiotics can select for tolerance mutations (Schuster et al., 2022) it is plausible to assume that presence of both virulence and ARGs would result in multiple drug resistance within constituent microorganism (Schroeder et al., 2017). However, there is still not a clear understanding of the co-selection relationship(s) between virulence and resistance factors, nor virulence and antibiotic concentrations (Beceiro et al., 2013; Matviichuk et al., 2022). Therefore, tracking AMR in environmental microbiomes is an essential step required to anticipate the spread and risk of ARGs/ARB transmission.
Environmental concentrations of ciprofloxacin in European WWTP effluent released into surface freshwater were higher than estimated and found to promote resistance in 5% of bacterial species in contact with the antibiotic (Emara et al., 2023). The same study reported a similar pattern for tetracycline, trimethoprim, ofloxacin, and norfloxacin. Recently, ciprofloxacin (among other antimicrobials) has been measured in the range of 542–32,800 ng/L in effluents close to two Canadian pharmaceutical manufacturing facilities (Kleywegt et al., 2019). It has also been determined in France that river biofilms exposed to WWTP discharge with oxytetracycline concentrations ≥ to 108 μg/L represented an environmental risk for aquatic bacteria (Matviichuk et al., 2022). These in-situ concentrations are within the range used in our present study (Supplementary Table S1), with values covering 1/10, 1/50, and 1/100 of MIC and ranging from 5,000 to 50,000 ng/L (below MICs of 500,000 ng/L), supporting the notion that global freshwater microbial communities are currently exposed to sub-MICs of antibiotics, with potential functional and structural impacts such as those observed in this study.
4.5. Broader microbial functional impacts
Significant changes in metabolic pathways found within the whole biofilm community under sub-MICs antibiotic exposure were also observed. In general, genes belonging to metabolic pathways related to biosynthesis processes were differentially abundant in treatment biofilms. For example, two peptidoglycan pathways were similarly abundant in each of the three sub-MIC antibiotic treatments; suggesting that even at the lowest concentration (5–5,120 μg/L) antibiotic presence can increase peptidoglycan biosynthesis. This upregulation of cell wall biosynthesis is likely in response to sub-MICs antibiotics exerting selective pressure on microbial communities, and in particular, the abundant group of Proteobacteria (Gram-negative bacteria) present. β-lactam antibiotics can interfere with the bacterial cell wall biosynthesis by inactivation of penicillin-binding proteins (PBPs), enzymes requisite for synthesis of the peptidoglycan layer (Pazda et al., 2019). β-lactam resistance mechanisms against Gram-negative bacteria mostly consist in β-lactamases, followed by permeability alterations and extrusion by efflux pumps (Beceiro et al., 2013); blaOXA-10 and blaOXA-24 have further been implicated in virulence properties by changing peptidoglycan composition (Fernández et al., 2012). Results from our functional prediction demonstrate that sub-MICs antibiotics enhanced the abundance of peptidoglycan-related metabolism, with a similar pattern observed for the abundance of blaOXA-1 and blaOXA-2 genes.
Functional metagenomic analysis also revealed shifts in gene abundance related to amino acid biosynthesis pathways at sub-MIC 1/10 and 1/50 conditions, including those responsible for synthesis of arginine, lysine, histidine, and isoleucine (Figure 9). Genes associated with lysine biosynthesis have been shown to be required in bacteria for protein synthesis, cell wall biosynthesis and generation of diaminopimelate (DAP), the latter being a key cell wall component in many Gram-negative bacteria (Grundy et al., 2003). Nucleotide and nucleoside biosynthesis genes consistently showed reduced (negative fold-change) gene abundance in treated samples relative to control samples. These results indicate potential disruptions in bacterial growth, since genes associated with these pathways are key in synthesizing precursor metabolites for central anabolic processes (Caspi et al., 2018). Considering that antibiotic resistance often causes fitness costs (Andersson and Hughes, 2012), shifts in genes associated with nucleoside and nucleotide pathways and overall growth rates could be related to the increased replication demands from more abundant MGE such as plasmids or integrons.
Interestingly, no other metabolic pathway genes associated with stress response or resistance displayed significant differential abundance in response to sub-MICs antibiotic exposure. In addition, no significant variations in carbon-related gene abundance were observed. This was not surprising given that carbon biosynthesis would largely result from autotrophic and/or microalgae organisms (non-prokaryotic) processes which were excluded from bacterial gene abundance screen regarding metabolic pathway analyses. Overall, sub-MIC antibiotic exposure seemingly influences cell wall biosynthesis and processes related to bacterial growth regulation. As biofilms are key participants in biogeochemical organic and inorganic cycles, including primary productivity, alterations to biofilm growth dynamics can potentially lead to downstream impacts on these environmental functions.
4.6. Microbiome and resistome linkages post-antibiotic exposure
Network analysis patterns from biofilm communities at different sub-MIC antibiotic conditions demonstrated that the overall resistome composition correlated to Proteobacteria genera.
Proteobacteria genera such as Bradyrhizobium, Methylibium, Methylorubrum, Novosphingobium, Sphingobium and Variovorax showed the highest correlations (number of edges) to the resistome (Figure 10; Supplementary Figure S11) suggesting that the presence of these genera strongly contribute to shaping the microbial composition and that they may play important ecological roles in riverine biofilms.
With respect to ARGs composition, aminoglycosides (aadA3, aph3” Ib, aph6-Id, rrsA), β-lactam (blaOXA-2, blaOXA-21), and sulfonamide (sul1) were the most prevalent across biofilms under sub-MICs antibiotic exposure (Figure 10). The aph6-Id gene displayed co-occurrence to aph3”-Ib in all sub-MICs treatments, and correlation with other ARGs increased with the concentration of sub-MICs antibiotics. Indeed, the co-occurrence of Proteobacteria, particularly Variovorax, suggests that this bacterial genera could act as potential host of aminoglycosides, β-lactams and sulfonamides and trimethoprim ARGs.
Moreover, Proteobacteria taxa, Mycobacterium and Mycobacteroides (Actinobacteria) showed correlation to blaOXA-2, sul1, and tetG genes. Occurrence of the sul1 gene was negatively correlated to Nitrospira, showing different observations from that of Matviichuk et al. (2022) who reported that sul1 was correlated with Bacteroidetes and Nitrospirae phyla from samples upstream of a WWTP, but with no correlation in downstream biofilms. Though co-occurrence patterns between resistome and microbiome across the different sub-MICs antibiotic conditions were observed, this approach and dataset is limited in its ability to ascribe with certainty the abundance of resistance genes to specific bacterial taxa (Supplementary Figure S11). However, the totality of our results suggests that predominance of Proteobacteria could serve as an indicator of AMR risk and aid in defining strategies for aquatic environmental surveillance efforts, given that the abundance of this taxa group is strongly correlated with ARG abundance.
5. Conclusion
This study demonstrated that the combined presence of sub-MICs oxytetracycline, ciprofloxacin and streptomycin in riverine biofilm communities decreased microbial diversity and selected for Proteobacteria genera including Methylorubrum and Variovorax, particularly at the sub-MIC 1/10 antibiotic level. In contrast, the relative abundance and richness of resistome composition increased at the sub-MIC 1/50 antibiotic level, suggesting greater opportunities for HGT events and transfer of MGE occurred under this condition. At sub-MIC 1/50 antibiotic conditions, exposure did not affect EPS matrix thickness as significantly as observed for the 1/10 treatment level, thereby potentially enabling enhanced cell–cell communication (quorum sensing) between biofilm microorganisms. Accordingly, there are a range of biofilm impacts and lifestyle trade-offs, with respect to structural stability and genetic diversity across the range of sub-MICs concentrations examined here. Our results have important implications for the overall functional capacity, ecosystem services, and broader ecological roles that sub-MICs antibiotic impacted biofilms play in aquatic ecosystems. Altogether, these observations confirm that the presence of environmentally relevant concentrations of antibiotics promote and accelerate the acquisition of resistance in aquatic microbial communities.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://www.ncbi.nlm.nih.gov/, BioProject ID: PRJNA944808.
Author contributions
GF-V wrote the manuscript and performed experiments and data analysis. JB wrote and reviewed manuscript and designed and supervised research. DK wrote and reviewed manuscript and designed and supervised research. All authors contributed to the article and approved the submitted version.
Funding
Financial support for this research was provided by Environment and Climate Change Canada and the University of Saskatchewan. Also, support of the Natural Sciences and Engineering Research Council of Canada (NSERC; grant to DK) is acknowledged. This work was also supported by Nous remercions le Conseil de Recherches en Sciences Naturelles et en Génie du Canada (CRSNG) de son soutien.
Acknowledgments
The authors thank Alixandra Schebel for assisting on sample collection and set-up of the microcosm experiment. The authors also thank George D. W. Swerhone and Yuriy Kryachko for helping with the microscopic analysis.
Conflict of interest
The authors declare that the research was conducted in the absence of any financial relationship that could be construed as potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2023.1194952/full#supplementary-material
Footnotes
References
Abe, K., Nomura, N., and Suzuki, S. (2020). Biofilms: hot spots of horizontal gene transfer (HGT) in aquatic environments, with a focus on a new HGT mechanism. FEMS Microbiol. Ecol. 96:31. doi: 10.1093/femsec/fiaa031
Andersson, A., Ahlinder, J., Mathisen, P., Hägglund, M., Bäckman, S., Nilsson, E., et al. (2018). Predators and nutrient availability favor protozoa-resisting bacteria in aquatic systems. Sci. Rep. 8:8415. doi: 10.1038/s41598-018-26422-4
Andersson, D. I., and Hughes, D. (2012). Evolution of antibiotic resistance at non-lethal drug concentrations. Drug Resist. Updat. 15, 162–172. doi: 10.1016/j.drup.2012.03.005
Andersson, D. I., and Hughes, D. (2014). Microbiological effects of sublethal levels of antibiotics. Nat. Rev. Microbiol. 12, 465–478. doi: 10.1038/nrmicro3270
Aubertheau, E., Stalder, T., Mondamert, L., Ploy, M. C., Dagot, C., and Labanowski, J. (2017). Impact of wastewater treatment plant discharge on the contamination of river biofilms by pharmaceuticals and antibiotic resistance. Sci. Total Environ. 579, 1387–1398. doi: 10.1016/j.scitotenv.2016.11.136
Auguet, O., Pijuan, M., Borrego, C. M., Rodriguez-Mozaz, S., Triadó-Margarit, X., Giustina, S. V. D., et al. (2017). Sewers as potential reservoirs of antibiotic resistance. Sci. Total Environ. 605-606, 1047–1054. doi: 10.1016/j.scitotenv.2017.06.153
Baharoglu, Z., Bikard, D., and Mazel, D. (2010). Conjugative DNA transfer induces the bacterial SOS response and promotes antibiotic resistance development through integron activation. PLoS Genet. 6, e1001165–e1001110. doi: 10.1371/journal.pgen.1001165
Baharoglu, Z., Krin, E., and Mazel, D. (2013). RpoS plays a central role in the SOS induction by sub-lethal aminoglycoside concentrations in Vibrio cholerae. PLoS Genet. 9, e1003421–e1003411. doi: 10.1371/journal.pgen.1003421
Balcázar, J. L., Subirats, J., and Borrego, C. M. (2015). The role of biofilms as environmental reservoirs of antibiotic resistance. Front. Microbiol. 6, 1–9. doi: 10.3389/fmicb.2015.01216
Battin, T. J., Besemer, K., Bengtsson, M. M., Romani, A. M., and Packmann, A. I. (2016). The ecology and biogeochemistry of stream biofilms. Nat. Rev. Microbiol. 14, 251–263. doi: 10.1038/nrmicro.2016.15
Beceiro, A., Tomás, M., and Bou, G. (2013). Antimicrobial resistance and virulence: a successful or deleterious association in the bacterial world? Clin. Microbiol. Rev. 26, 185–230. doi: 10.1128/CMR.00059-12
Beghini, F., McIver, L. J., Blanco-Míguez, A., Dubois, L., Asnicar, F., Maharjan, S., et al. (2021). Integrating taxonomic, functional, and strain-level profiling of diverse microbial communities with bioBakery 3. Elife 10:e65088. doi: 10.7554/eLife.65088
Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. Ser. B 57, 289–300. doi: 10.1111/j.2517-6161.1995.tb02031.x
Bergsveinson, J., Roy, J., Maynard, C., Sanschagrin, S., Freeman, C. N., Swerhone, G. D. W., et al. (2020). Metatranscriptomic insights into the response of river biofilm communities to ionic and nano-zinc oxide exposures. Front. Microbiol. 11, 1–15. doi: 10.3389/fmicb.2020.00267
Bolger, A. M., Lohse, M., and Usadel, B. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120. doi: 10.1093/bioinformatics/btu170
Bortolaia, V., Kaas, R. S., Ruppe, E., Roberts, M. C., Schwarz, S., Cattoir, V., et al. (2020). ResFinder 4.0 for predictions of phenotypes from genotypes. J. Antimicrob. Chemother. 75, 3491–3500. doi: 10.1093/jac/dkaa345
Brumfield, K. D., Hasan, N. A., Leddy, M. B., Cotruvo, J. A., Rashed, S. M., Colwell, R. R., et al. (2020). A comparative analysis of drinking water employing metagenomics. PLoS One 15, 1–27. doi: 10.1371/journal.pone.0231210
Cairns, J., Ruokolainen, L., Hultman, J., Tamminen, M., Virta, M., and Hiltunen, T. (2018). Ecology determines how low antibiotic concentration impacts community composition and horizontal transfer of resistance genes. Commun. Biol. 1:35. doi: 10.1038/s42003-018-0041-7
Carter, G., Young, L. S., and Bermudez, L. E. (2004). A subinhibitory concentration of clarithromycin inhibits Mycobacterium avium biofilm formation. Antimicrob. Agents Chemother. 48, 4907–4910. doi: 10.1128/AAC.48.12.4907-4910.2004
Caspi, R., Billington, R., Fulcher, C. A., Keseler, I. M., Kothari, A., Krummenacker, M., et al. (2018). The MetaCyc database of metabolic pathways and enzymes. Nucleic Acids Res. 46, D633–D639. doi: 10.1093/nar/gkx935
Cerca, N., Martins, S., Sillankorva, S., Jefferson, K. K., Pier, G. B., Oliveira, R., et al. (2005). Effects of growth in the presence of subinhibitory concentrations of dicloxacillin on Staphylococcus epidermidis and Staphylococcus haemolyticus biofilms. Appl. Environ. Microbiol. 71, 8677–8682. doi: 10.1128/AEM.71.12.8677-8682.2005
Chekabab, S. M., Lawrence, J. R., Alvarado, A., Predicala, B., and Korber, D. R. (2020). A health metadata-based management approach for comparative analysis of high-throughput genetic sequences for quantifying antimicrobial resistance reduction in Canadian hog barns. Comput. Struct. Biotechnol. J. 18, 2629–2638. doi: 10.1016/j.csbj.2020.09.012
Cheng, D., Hao Ngo, H., Guo, W., Wang Chang, S., Duc Nguyen, D., Liu, Y., et al. (2020). Contribution of antibiotics to the fate of antibiotic resistance genes in anaerobic treatment processes of swine wastewater: a review. Bioresour. Technol. 299:122654. doi: 10.1016/j.biortech.2019.122654
Chonova, T., Kurmayer, R., Rimet, F., Labanowski, J., Vasselon, V., Keck, F., et al. (2019). Benthic diatom communities in an alpine river impacted by waste water treatment effluents as revealed using DNA metabarcoding. Front. Microbiol. 10, 1–17. doi: 10.3389/fmicb.2019.00653
Chow, L. K. M., Ghaly, T. M., and Gillings, M. R. (2021). A survey of sub-inhibitory concentrations of antibiotics in the environment. J. Environ. Sci. 99, 21–27. doi: 10.1016/j.jes.2020.05.030
CLSI (2018). Performance standards for antimicrobial susceptibility testing. 28th ed. CLSI supplement M100. Wayne, PA:
Connelly, S., Bristol, J. A., Hubert, S., Subramanian, P., Hasan, N. A., Colwell, R. R., et al. (2017). SYN-004 (ribaxamase), an oral beta-lactamase, mitigates antibiotic-mediated dysbiosis in a porcine gut microbiome model. J. Appl. Microbiol. 123, 66–79. doi: 10.1111/jam.13432
Couto, N., Schuele, L., Raangs, E. C., Machado, M. P., Mendes, C. I., Jesus, T. F., et al. (2018). Critical steps in clinical shotgun metagenomics for the concomitant detection and typing of microbial pathogens. Sci. Rep. 8:13767. doi: 10.1038/s41598-018-31873-w
Danner, M. C., Robertson, A., Behrends, V., and Reiss, J. (2019). Antibiotic pollution in surface fresh waters: occurrence and effects. Sci. Total Environ. 664, 793–804. doi: 10.1016/j.scitotenv.2019.01.406
Davies, J., Spiegelman, G. B., and Yim, G. (2006). The world of subinhibitory antibiotic concentrations. Curr. Opin. Microbiol. 9, 445–453. doi: 10.1016/j.mib.2006.08.006
Delattre, C., Pierre, G., Laroche, C., and Michaud, P. (2016). Production, extraction and characterization of microalgal and cyanobacterial exopolysaccharides. Biotechnol. Adv. 34, 1159–1179. doi: 10.1016/j.biotechadv.2016.08.001
Dynes, J. J., Lawrence, J. R., Korber, D. R., Swerhone, G. D. W., Leppard, G. G., and Hitchcock, A. P. (2006). Quantitative mapping of chlorhexidine in natural river biofilms. Sci. Total Environ. 369, 369–383. doi: 10.1016/j.scitotenv.2006.04.033
Dynes, J. J., Lawrence, J. R., Korber, D. R., Swerhone, G. D. W., Leppard, G. G., and Hitchcock, A. P. (2009). Morphological and biochemical changes in Pseudomonas fluorescens biofilms induced by sub-inhibitory exposure to antimicrobial agents. Can. J. Microbiol. 55, 163–178. doi: 10.1139/W08-109
Ebmeyer, S., Kristiansson, E., and Larsson, D. G. J. (2021). A framework for identifying the recent origins of mobile antibiotic resistance genes. Commun. Biol. 4, 8–10. doi: 10.1038/s42003-020-01545-5
Emara, Y., Jolliet, O., Finkbeiner, M., Heß, S., Kosnik, M., Siegert, M. W., et al. (2023). Comparative selective pressure potential of antibiotics in the environment. Environ. Pollut. 318:120873. doi: 10.1016/j.envpol.2022.120873
EUCAST (2019). The European committee on antimicrobial susceptibility testing. Breakpoint tables for interpretation of MICs and zone diameters. Version 9.0. Available at: http://www.eucast.org.
Fernández, A., Pérez, A., Ayala, J. A., Mallo, S., Rumbo-Feal, S., Tomás, M., et al. (2012). Expression of OXA-type and SFO-1 β-lactamases induces changes in peptidoglycan composition and affects bacterial fitness. Antimicrob. Agents Chemother. 56, 1877–1884. doi: 10.1128/AAC.05402-11
Flemming, H.-C., van Hullebusch, E. D., Neu, T. R., Nielsen, P. H., Seviour, T., Stoodley, P., et al. (2022). The biofilm matrix: multitasking in a shared space. Nat. Rev. Microbiol. 21, 70–86. doi: 10.1038/s41579-022-00791-0
Flores-Vargas, G., Bergsveinson, J., Lawrence, J. R., and Korber, D. R. (2021). Environmental biofilms as reservoirs for antimicrobial resistance. Front. Microbiol. 12:766242. doi: 10.3389/fmicb.2021.766242
Fulaz, S., Vitale, S., Quinn, L., and Casey, E. (2019). Nanoparticle–biofilm interactions: the role of the EPS matrix. Trends Microbiol. 27, 915–926. doi: 10.1016/j.tim.2019.07.004
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
Gonzalez-Martinez, A., Rodriguez-Sanchez, A., Martinez-Toledo, M. V., Garcia-Ruiz, M.-J., Hontoria, E., Osorio-Robles, F., et al. (2014). Effect of ciprofloxacin antibiotic on the partial-nitritation process and bacterial community structure of a submerged biofilter. Sci. Total Environ. 476-477, 276–287. doi: 10.1016/j.scitotenv.2014.01.012
Grundy, F. J., Lehman, S. C., and Henkin, T. M. (2003). The L box regulon: lysine sensing by leader RNAs of bacterial lysine biosynthesis genes. Proc. Natl. Acad. Sci. 100, 12057–12062. doi: 10.1073/pnas.2133705100
Guo, X. P., Sun, X. L., Chen, Y. R., Hou, L., Liu, M., and Yang, Y. (2020). Antibiotic resistance genes in biofilms on plastic wastes in an estuarine environment. Sci. Total Environ. 745:140916. doi: 10.1016/j.scitotenv.2020.140916
Guo, X., Yang, Y., Lu, D., Niu, Z., Feng, J., Chen, Y., et al. (2018). Biofilms as a sink for antibiotic resistance genes (ARGs) in the Yangtze estuary. Water Res. 129, 277–286. doi: 10.1016/j.watres.2017.11.029
Hanlon, A. R. M., Bellinger, B., Haynes, K., Xiao, G., Hofmann, T. A., Gretz, M. R., et al. (2006). Dynamics of extracellular polymeric substance (EPS) production and loss in an estuarine, diatom-dominated, microalgal biofilm over a tidal emersio—immersion period. Limnol. Oceanogr. 51, 79–93. doi: 10.4319/lo.2006.51.1.0079
Hoffman, L. R., D’Argenio, D. A., MacCoss, M. J., Zhang, Z., Jones, R. A., and Miller, S. I. (2005). Aminoglycoside antibiotics induce bacterial biofilm formation. Nature 436, 1171–1175. doi: 10.1038/nature03912
Horii, T., Morita, M., Muramatsu, H., Muranaka, Y., Kanno, T., and Maekawa, M. (2003). Effects of mupirocin at subinhibitory concentrations on flagella formation in Pseudomonas aeruginosa and Proteus mirabilis. J. Antimicrob. Chemother. 51, 1175–1179. doi: 10.1093/jac/dkg226
Jia, B., Raphenya, A. R., Alcock, B., Waglechner, N., Guo, P., Tsang, K. K., et al. (2016). CARD 2017: expansion and model-centric curation of the comprehensive antibiotic resistance database. Nucleic Acids Res. 45, D566–D573. doi: 10.1093/nar/gkw1004
Junqueira, A. C., Ratan, A., Acerbi, E., Drautz-Moses, D. I., Premkrishnan, B. N. V., Costea, P. I., et al. (2017). The microbiomes of blowflies and houseflies as bacterial transmission reservoirs. Sci. Rep. 7:16324. doi: 10.1038/s41598-017-16353-x
Karkman, A., Do, T. T., Walsh, F., and Virta, M. P. J. (2018). Antibiotic-resistance genes in waste water. Trends Microbiol. 26, 220–228. doi: 10.1016/j.tim.2017.09.005
Kleywegt, S., Payne, M., Ng, F., and Fletcher, T. (2019). Environmental loadings of active pharmaceutical ingredients from manufacturing facilities in Canada. Sci. Total Environ. 646, 257–264. doi: 10.1016/j.scitotenv.2018.07.240
Kraemer, S. A., Ramachandran, A., and Perron, G. G. (2019). Antibiotic pollution in the environment: from microbial ecology to public policy. Microorganisms 7:180. doi: 10.3390/microorganisms7060180
Kunhikannan, S., Thomas, C. J., Franks, A. E., Mahadevaiah, S., Kumar, S., and Petrovski, S. (2021). Environmental hotspots for antibiotic resistance genes. Microbiology 10:e1197. doi: 10.1002/mbo3.1197
Langmead, B., and Salzberg, S. L. (2012). Fast gapped-read alignment with bowtie 2. Nat. Methods 9, 357–359. doi: 10.1038/nmeth.1923
Laureti, L., Matic, I., and Gutierrez, A. (2013). Bacterial responses and genome instability induced by subinhibitory concentrations of antibiotics. Antibiotics 2, 100–114. doi: 10.3390/antibiotics2010100
Lawrence, J. R., Chenier, M. R., Roy, R., Beaumier, D., Fortin, N., Swerhone, G. D. W., et al. (2004). Microscale and molecular assessment of impacts of nickel, nutrients, and oxygen level on structure and function of river biofilm communities. Appl. Environ. Microbiol. 70, 4326–4339. doi: 10.1128/AEM.70.7.4326-4339.2004
Lawrence, J. R., Paule, A., Swerhone, G. D. W., Roy, J., Grigoryan, A. A., Dynes, J. J., et al. (2020). Microscale and molecular analyses of river biofilm communities treated with microgram levels of cerium oxide nanoparticles indicate limited but significant effects. Environ. Pollut. 256:113515. doi: 10.1016/j.envpol.2019.113515
Lawrence, J. R., Swerhone, G. D. W., Kuhlicke, U., and Neu, T. R. (2007a). In situ evidence for microdomains in the polymer matrix of bacterial microcolonies. Can. J. Microbiol. 53, 450–458. doi: 10.1139/W06-146
Lawrence, J. R., Swerhone, G. D. W., and Neu, T. R. (2000). A simple rotating annular reactor for replicated biofilm studies. J. Microbiol. Methods 42, 215–224. doi: 10.1016/S0167-7012(00)00195-0
Lawrence, J. R., Swerhone, G. D. W., Topp, E., Korber, D. R., Neu, T. R., and Wassenaar, L. I. (2007b). Structural and functional responses of river biofilm communities to the nonsteroidal anti-inflammatory diclofenac. Environ. Toxicol. Chem. 26, 573–582. doi: 10.1897/06-340r.1
Lawrence, J. R., Swerhone, G. D. W., Wassenaar, L. I., and Neu, T. R. (2005). Effects of selected pharmaceuticals on riverine biofilm communities. Can. J. Microbiol. 51, 655–669. doi: 10.1139/w05-047
Lawrence, J. R., Zhu, B., Swerhone, G. D. W., Roy, J., Wassenaar, L. I., Topp, E., et al. (2009). Comparative microscale analysis of the effects of triclosan and triclocarban on the structure and function of river biofilm communities. Sci. Total Environ. 407, 3307–3316. doi: 10.1016/j.scitotenv.2009.01.060
Lépesová, K., Kraková, L., Pangallo, D., Medveďová, A., Olejníková, P., Mackuľak, T., et al. (2018). Prevalence of antibiotic-resistant coliform bacteria, enterococcus spp. and staphylococcus spp. in wastewater sewerage biofilm. J. Glob. Antimicrob. Resist. 14, 145–151. doi: 10.1016/j.jgar.2018.03.008
Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., et al. (2009). The sequence alignment/map format and SAMtools. Bioinformatics 25, 2078–2079. doi: 10.1093/bioinformatics/btp352
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
Liu, B., Zheng, D., Jin, Q., Chen, L., and Yang, J. (2019). VFDB 2019: a comparative pathogenomic platform with an interactive web interface. Nucleic Acids Res. 47, D687–D692. doi: 10.1093/nar/gky1080
Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Gen. Biol. 15, 550. doi: 10.1186/s13059-014-0550-8
Luo, X., Xiang, X., Yang, Y., Huang, G., Fu, K., Che, R., et al. (2020). Seasonal effects of river flow on microbial community coalescence and diversity in a riverine network. FEMS Microbiol. Ecol. 96:fiaa132. doi: 10.1093/femsec/fiaa132
Maestre-Carballa, L., Lluesma Gomez, M., Angla Navarro, A., Garcia-Heredia, I., Martinez-Hernandez, F., and Martinez-Garcia, M. (2019). Insights into the antibiotic resistance dissemination in a wastewater effluent microbiome: bacteria, viruses and vesicles matter. Environ. Microbiol. 21, 4582–4596. doi: 10.1111/1462-2920.14758
Majidpour, A., Fathizadeh, S., Afshar, M., Rahbar, M., Boustanshenas, M., Heidarzadeh, M., et al. (2017). Dose-dependent effects of common antibiotics used to treat Staphylococcus aureus on biofilm formation. Iran. J. Pathol. 12, 362–370. doi: 10.30699/ijp.2017.27993
Manaia, C. M. (2017). Assessing the Rrsk of antibiotic resistance transmission from the environment to humans: non-direct proportionality between abundance and risk. Trends Microbiol. 25, 173–181. doi: 10.1016/j.tim.2016.11.014
Matviichuk, O., Mondamert, L., Geffroy, C., Gaschet, M., Dagot, C., and Labanowski, J. (2022). River biofilms microbiome and resistome responses to wastewater treatment plant effluents containing antibiotics. Front. Microbiol. 13:5206. doi: 10.3389/fmicb.2022.795206
Munk, P., Brinch, C., Møller, F. D., Petersen, T. N., Hendriksen, R. S., Seyfarth, A. M., et al. (2022). Genomic analysis of sewage from 101 countries reveals global landscape of antimicrobial resistance. Nat. Commun. 13, 7251–7216. doi: 10.1038/s41467-022-34312-7
Neu, T. R., Swerhone, G. D. W., and Lawrence, J. R. (2001). Assessment of lectin-binding analysis for in situ detection of glycoconjugates in biofilm systems. Microbiology 147, 299–313. doi: 10.1099/00221287-147-2-299
O’Leary, N. A., Wright, M. W., Brister, J. R., Ciufo, S., Haddad, D., McVeigh, R., et al. (2016). Reference sequence (RefSeq) database at NCBI: current status, taxonomic expansion, and functional annotation. Nucleic Acids Res. 44, D733–D745. doi: 10.1093/nar/gkv1189
Oksanen, J., Blanchet, F. G., Friendly, M., Kindt, R., Legendre, P., McGlinn, D., et al. (2020). Vegan: Community ecology package. R package. Available at: https://cran.r-project.org/package=vegan
Packroff, G., Lawrence, J. R., and Neu, T. R. (2002). In situ confocal laser scanning microscopy of protozoans in cultures and complex biofilm communities. Acta Protozool. 41, 245–253.
Pazda, M., Kumirska, J., Stepnowski, P., and Mulkiewicz, E. (2019). Antibiotic resistance genes identified in wastewater treatment plant systems – a review. Sci. Total Environ. 697:134023. doi: 10.1016/j.scitotenv.2019.134023
Petrovich, M. L., Ben Maamar, S., Hartmann, E. M., Murphy, B. T., Poretsky, R. S., and Wells, G. F. (2019). Viral composition and context in metagenomes from biofilm and suspended growth municipal wastewater treatment plants. J. Microbial. Biotechnol. 12, 1324–1336. doi: 10.1111/1751-7915.13464
Petrovich, M. L., Zilberman, A., Kaplan, A., Eliraz, G. R., Wang, Y., Langenfeld, K., et al. (2020). Microbial and viral communities and their antibiotic resistance genes throughout a hospital wastewater treatment system. Front. Microbiol. 11, 1–13. doi: 10.3389/fmicb.2020.00153
Pierre, G., Graber, M., Rafiliposon, B. A., Dupuy, C., Orvain, F., De Crignis, M., et al. (2012). Biochemical composition and changes of extracellular polysaccharides (ECPS) produced during microphytobenthic biofilm development (Marennes-Oléron, France). Microb. Ecol. 63, 157–169. doi: 10.1007/s00248-011-9959-8
R Core Team (2020). R: A language and environment for statistical computing Vienna: R Foundation for Statistical Computing. Available at: https:// www.R-project.org/
Reichhardt, C., Jacobs, H. M., Matwichuk, M., Wong, C., Wozniak, D. J., and Parsek, M. R. (2020). The versatile Pseudomonas aeruginosa biofilm matrix protein CdrA promotes aggregation through different extracellular exopolysaccharide interactions. J. Bacteriol. 202, e00216–e00220. doi: 10.1128/JB.00216-20
Romaní, A. M., Borrego, C. M., Díaz-Villanueva, V., Freixa, A., Gich, F., and Ylla, I. (2014). Shifts in microbial community structure and function in light- and dark-grown biofilms driven by warming. Environ. Microbiol. 16, 2550–2567. doi: 10.1111/1462-2920.12428
Romero, F., Sabater, S., Timoner, X., and Acuña, V. (2018). Multistressor effects on river biofilms under global change conditions. Sci. Total Environ. 627, 1–10. doi: 10.1016/j.scitotenv.2018.01.161
Sanchez-Cid, C., Guironnet, A., Keuschnig, C., Wiest, L., Vulliet, E., and Vogel, T. M. (2022). Gentamicin at sub-inhibitory concentrations selects for antibiotic resistance in the environment. ISME Commun. 2, 1–9. doi: 10.1038/s43705-022-00101-y
Schneider, C. A., Rasband, W. S., and Eliceiri, K. W. (2012). NIH image to ImageJ: 25 years of image analysis. Nat. Methods 9, 671–675. doi: 10.1038/nmeth.2089
Schroeder, M., Brooks, B. D., and Brooks, A. E. (2017). The complex relationship between virulence and antibiotic resistance. Genes 8:39. doi: 10.3390/genes8010039
Schuster, D., Axtmann, K., Holstein, N., Felder, C., Voigt, A., Färber, H., et al. (2022). Antibiotic concentrations in raw hospital wastewater surpass minimal selective and minimum inhibitory concentrations of resistant Acinetobacter baylyi strains. Environ. Microbiol. 24, 5721–5733. doi: 10.1111/1462-2920.16206
Seemann, T. (2014). Prokka: rapid prokaryotic genome annotation. Bioinformatics 30, 2068–2069. doi: 10.1093/bioinformatics/btu153
Seemann, T. (2019). Abricate. Github. Available at: https://github.com/tseemann/abricate
Segata, N., Waldron, L., Ballarini, A., Narasimhan, V., Jousson, O., Huttenhower, C., et al. (2012). Metagenomic microbial community profiling using unique clade-specific marker genes. Nat. Methods 9, 811–814. doi: 10.1038/nmeth.2066
Shade, A. (2017). Diversity is the question, not the answer. ISME J. 11, 1–6. doi: 10.1038/ismej.2016.118
Subirats, J., Timoner, X., Sànchez-Melsió, A., Balcázar, J. L., Acuña, V., Sabater, S., et al. (2018). Emerging contaminants and nutrients synergistically affect the spread of class 1 integron-integrase (intI1) and sul1 genes within stable streambed bacterial communities. Water Res. 138, 77–85. doi: 10.1016/j.watres.2018.03.025
Sun, S., Noorian, P., and McDougald, D. (2018). Dual role of mechanisms involved in resistance to predation by protozoa and virulence to humans. Front. Microbiol. 9:1017. doi: 10.3389/fmicb.2018.01017
Szczuka, E., Jabłońska, L., and Kaznowski, A. (2017). Effect of subinhibitory concentrations of tigecycline and ciprofloxacin on the expression of biofilm-associated genes and biofilm structure of Staphylococcus epidermidis. Microbiology 163, 712–718. doi: 10.1099/mic.0.000453
WHO (2017). Integrated surveillance of antimicrobial resistance in foodborne bacteria: Application of a one health approach. Available at: https://apps.who.int/iris/bitstream/handle/10665/255747/9789241512411-eng.pdf;jsessionid=1710DEAA0E355CBCB99450559C0DD0C8?sequence=1%0Ahttp://apps.who.int/iris/bitstream/10665/255747/1/9789241512411-eng.pdf?ua=1%0Ahttp://apps.who.int/iris/bitstream/10665/.
Wood, D. E., Lu, J., and Langmead, B. (2019). Improved metagenomic analysis with kraken 2. Genome Biol. 20:257. doi: 10.1186/s13059-019-1891-0
Yadav, S., and Kapley, A. (2021). Antibiotic resistance: global health crisis and metagenomics. Biotechnol. Rep 29:e00604. doi: 10.1016/j.btre.2021.e00604
Yan, Q., Mi Wi, Y., Thoendel, M. J., Raval, Y. S., Greenwood-Quaintance, K. E., Abdel, M. P., et al. (2019). Evaluation of the CosmosID bioinformatics platform for prosthetic joint-associated sonicate fluid shotgun metagenomic data analysis. J. Clin. Microbiol. 57, 1–13. doi: 10.1128/JCM.01182-18
Zaouri, N., Jumat, M. R., Cheema, T., and Hong, P.-Y. (2020). Metagenomics-based evaluation of groundwater microbial profiles in response to treated wastewater discharge. Environ. Res. 180:108835. doi: 10.1016/j.envres.2019.108835
Zhang, R. Y., Neu, T. R., Bellenberg, S., Kuhlicke, U., Sand, W., and Vera, M. (2015). Use of lectins to in situ visualize glycoconjugates of extracellular polymeric substances in acidophilic archaeal biofilms. J. Microbial. Biotechnol. 8, 448–461. doi: 10.1111/1751-7915.12188
Zhang, Y., Xie, J., Liu, M., Tian, Z., He, Z., van Nostrand, J. D., et al. (2013). Microbial community functional structure in response to antibiotics in pharmaceutical wastewater treatment systems. Water Res. 47, 6298–6308. doi: 10.1016/j.watres.2013.08.003
Keywords: antibiotics, antibiotic resistance, biofilms, environmental resistome, microbiome, sub-MIC
Citation: Flores-Vargas G, Korber DR and Bergsveinson J (2023) Sub-MIC antibiotics influence the microbiome, resistome and structure of riverine biofilm communities. Front. Microbiol. 14:1194952. doi: 10.3389/fmicb.2023.1194952
Edited by:
Nona Sheila Romualdo Agawin, University of the Balearic Islands, SpainReviewed by:
Maialen Barret, Université de Toulouse, FranceJingfeng Wang, Tianjin Institute of Environmental and Operational Medicine, China
Copyright © 2023 Flores-Vargas, Korber and Bergsveinson. 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(s) 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.
*Correspondence: Jordyn Bergsveinson, am9yZHluLmJyb2FkYmVudEBlYy5nYy5jYQ==