River Biofilms Microbiome and Resistome Responses to Wastewater Treatment Plant Effluents Containing Antibiotics.

Continuous exposure to low concentrations of antibiotics (sub-minimal inhibitory concentration: sub-MIC) is thought to lead to the development of antimicrobial resistance (AMR) in the environmental microbiota. However, the relationship between antibiotic exposure and resistance selection in environmental bacterial communities is still poorly understood and unproven. Therefore, we measured the concentration of twenty antibiotics, resistome quality, and analyzed the taxonomic composition of microorganisms in river biofilms collected upstream (UPS) and downstream (DWS) (at the point of discharge) from the wastewater treatment plant (WWTP) of Poitiers (France). The results of statistical analysis showed that the antibiotic content, resistome, and microbiome composition in biofilms collected UPS were statistically different from that collected DWS. According to Procrustes analysis, microbial community composition and antibiotics content may be determinants of antibiotic resistance genes (ARGs) composition in samples collected DWS. However, network analysis showed that the occurrence and concentration of antibiotics measured in biofilms did not correlate with the occurrence and abundance of antibiotic resistance genes and mobile genetic elements. In addition, network analysis suggested patterns of co-occurrence between several ARGs and three classes of bacteria/algae: Bacteroidetes incertae sedis, Cyanobacteria/Chloroplast, and Nitrospira, in biofilm collected UPS. The absence of a direct effect of antibiotics on the selection of resistance genes in the collected samples suggests that the emergence of antibiotic resistance is probably not only due to the presence of antibiotics but is a more complex process involving the cumulative effect of the interaction between the bacterial communities (biotic) and the abiotic matrix. Nevertheless, this study confirms that WWTP is an important reservoir of various ARGs, and additional efforts and legislation with clearly defined concentration limits for antibiotics and resistance determinants in WWTP effluents are needed to prevent their spread and persistence in the environment.

Introduction

Antibiotics (ATBs) are now considered to be ubiquitous contaminants in all kinds of environmental matrices (sediments, soils, water, plants and animals) and present worldwide (Li et al., 2012; Gothwal and Shashidhar, 2015; Aus Der Beek et al., 2016; Patel et al., 2019). Their main origin is related to their excretion and persistence in human or animal feces after administration. In France, there is a whole system of wastewater collection and treatment, either in wastewater treatment plants (WWTP) or, in the case of animal waste (such as manure and slurry) in pits, digesters or methanation units. But all of these treatment methods used are not adapted/optimized to remove antibiotics and other pharmaceuticals (Luo et al., 2014; Kim and Zoh, 2016). Treated effluents, therefore, contain antibiotic residues (and/or their transformation products) which are released into the environment when the effluent is discharged or reused (agricultural spreading). In rivers, these discharges lead to a constant presence of residual concentrations of ATBs and, consequently, exposure of naturally occurring aquatic bacteria that live in structured communities called biofilms (Besemer, 2016).

It has already been found that continuous exposure of natural biofilms to low concentrations of ATBs can lead to various consequences, such as changes in taxonomic composition (Cairns et al., 2018; Danner et al., 2019; Tian et al., 2020), effects on metabolic activity (Crabbé et al., 2019), or even cell death (Friman et al., 2015; Gebreyohannes et al., 2019). However, the most important concern relates to the emergence of antimicrobial resistance (AMR) (Gullberg, 2014; Flemming et al., 2016; Khan et al., 2017). Continuous exposure to low concentrations (sub-minimal inhibitory concentration: sub-MIC) can’t only increase the relative number of antibiotic-resistant bacteria, but also select them for resistance, promote adaptive evolution, and even have unexpected consequences for trophic levels and natural microbial ecosystem functioning (Friman et al., 2015; Hathroubi et al., 2017; Danner et al., 2019).

Although knowledge about the occurrence, transmission, evolution, acquisition, selection, preservation, etc., of AMR and all related theoretical mechanisms have advanced in the last decade, especially through the “One Health” approach, they are still incomplete (Singer et al., 2016; Bengtsson-Palme et al., 2018). The main research has been conducted in the medical/clinical field and has linked the spread of AMR to antimicrobial misuse and overuse as a major factor in the maintenance of AMR. In contrast, the natural variability of resistance genes (RGs), mobile genetic elements (MGEs) and the frequency of their transmission in different types of natural environments are not sufficiently studied (Zhuang et al., 2021). A recent survey of AMR-related research funding showed that only 3% of projects involved research at the environmental scale (Singer et al., 2016). Under natural conditions, the relationship between antibiotic exposure and resistance selection in bacterial communities is less understood. However, awareness of the residual presence of ATBs raises many questions about the role of the environment in the emergence, dissemination and/or maintenance of AMR. Furthermore, it has been shown that resistant bacteria in natural ecosystems have the potential to transfer resistance to pathogenic species through horizontal gene transfer (Stanton et al., 2020). Moreover, there is no simple solution/technology that could be effective in an emergency to clean up all natural environments and stop an explosive spread of AMR (Larsson et al., 2018; Stanton et al., 2020). Thus, repeated/continuous release of antibiotics and resistant bacteria could ultimately pose a significant risk to public health if it makes the environment an uncontrollable reservoir of AMR. The objective of this work was to determine whether current antibiotic contamination of the environment contributes to the spread of AMR in rivers. Therefore, this study examined the relationship between antibiotic concentrations, taxonomic composition and the abundance of antibiotic resistance markers (genes and MGEs) in microbial communities during a 1-year field study.

Materials and Methods

Sampling Area and Sample Collection

Natural epilithic biofilm was collected from the Clain river located in the Midwest of France (Supplementary Figure 1). Sampling was conducted within the city of Poitiers at two sites: 1st—upstream (UPS) of the Poitiers WWTP (approximately 2 km from the WWTP discharge area) and 2nd—downstream (DWS) of the WWTP (approximately 5 m from the treated wastewater discharge point). The Poitiers’ WWTP has a nominal capacity of 152,500 inhabitant equivalents (per capita) and treats an average of 20,797 m3 of urban wastewater per day (in 2019) (more information provided in the Supplementary Information). Briefly, sterile rocks were placed in the river 5 months before the experiment to allow the formation of natural epilithic biofilms from the river’s microbial communities. At the end of this step, allowing the biofilm to establish, develop, and be exposed, the stones were collected to obtain biofilm samples. A 1-year sampling campaign was conducted from January to December 2018. Beginning from January, several rocks were gradually collected each month, which corresponded to a single bulk sample. Biofilm was carefully scraped off the rocks using a sterile toothbrush that was constantly rinsed with MilliQ grade water. The obtained solution was then carefully centrifuged to remove excess water. The resulting “fresh” biofilms were stored in bottles at −80°C until analyzed.

Sample’s Preparation and Quantification of Antibiotics

Twenty investigated antibiotics—enrofloxacin (ENR), erythromycin (ERY), levofloxacin (LVX), norfloxacin (NFX), oxytetracycline (OTC), sulfamethazine (SMZ), sulfamethoxazole (SMX), sulfaquinoxaline (SQX), trimethoprim (TMP), ciprofloxacin (CPR), metronidazole (MTZ), roxithromycin (RXM), clarithromycin (CLR), flumequine (FMQ), enoxacin (ENX), tylosin tartrate (T-T), midecamycin (MED), spiramycin (SPR), josamycin (JOS), azithromycin (AZM) (Supplementary Table 1) were purchased from Sigma Aldrich (Buchs, Switzerland). LC/MS-grade methanol, acetonitrile and formic acid were purchased from Carlo Erba Reagents (Val-de-Reuil, France), ultrapure water (>18.2 MΩ cm–1) was prepared with the Milli-Q IQ 7000 system (Millipore SAS, Molsheim, France). Individual raw solutions (200 mg/L) were prepared for each antibiotic by dissolving the standard powder in methanol and stored in the dark at −20°C. The standard solution mixtures were prepared at three concentrations (10, 5, 1 mg/L) by diluting the raw solutions with water and then stored in the dark at 4°C. Antibiotics were measured in all the collected biofilms by the procedure described in Aubertheau et al. (2017).

DNA Extraction and Resistome Analysis

Total bacterial DNA was extracted from freeze-dried biofilms using the Fast DNA® SPIN Kit for Feces (MP Biomedicals, Illkirch, France) according to the manufacturer’s protocol. DNA concentrations were measured using a NanoDrop™ OneC Microvolume UV-Vis Spectrophotometer (Thermo Fisher Scientific, United States). Class 1, 2 and 3 integrons were measured by qPCR analysis using the MX3005P real-time detection system (Stratagene®) according to the Taqman method described by Stalder et al. (2014). The rest of RGs and MGEs (91 in total) were analyzed using nanoliter-scale quantitative PCR. Between them 66 ARGs, 5 multidrug efflux pumps, 6 heavy metal resistance genes (MRGs), 3 disinfectant RGs and 11 MGEs. The list of RGs and their corresponding primer sequences are given in Supplementary Table 5. The analysis methodology was performed using 96.96 BioMark® Dynamic Array for Real-Time PCR (Fluidigm Corporation, San Francisco, CA, United States). A threshold cycle (Ct) values were extracted using BioMark Real-Time PCR analysis software. The detection limit was set at 20, and the normalized gene abundance was calculated using the 16S rRNA gene abundance according to the following formula: 2{−[. The analysis procedure is described in more detail by Buelow et al. (2018, 2020).

Microbial Community Analysis

Part of each sample was sent to a European sequencing platform (MACROGEN, EUROPE B.V., Amsterdam, Netherlands[1]) to perform microbial community analysis. There, total DNA was extracted using the PowerSoil® DNA Isolation Kit (QBioGene, Solon, OH, United States) following the manufacturer’s instructions. Sequencing samples were prepared according to Illumina 16S Metagenomic sequencing library protocols. The hypervariable V3-V4 16S region of the bacterial gene was targeted by the universal bacterial primers PCR 341F and 805R (Herlemann et al., 2011). Gene amplification and sequencing [using the Illumina platform (Illumina Inc., San Diego, CA, United States)] were performed by Macrogen Ltd. (Seoul, Korea), with the Herculase II Fusion DNA Polymerase Nextera XT Index kit V2 (library kit) and following the 16S metagenomics sequencing library preparation protocol[2].

Assembly was performed using the FLASH (Fast Length Adjustment of SHort Reads) software tool to merge paired-end reads. The raw data were processed using CD-HIT-OTU. Short reads were filtered out and ultra-long tails were trimmed. The filtered reads were clustered at 100% identity (using CD-HIT-DUP) and chimeric reads were identified. Secondary clusters were then joined to primary clusters and noise sequences from these clusters were removed. Finally, the remaining representative reads from non-chimeric clusters were clustered using a greedy algorithm into an operational taxonomic unit (OTU) (at the species level). OTU clustering was used in the de novo OTU selection method (CD-HIT). Taxonomic distribution and diversity statistics were performed using Qiime 2 version 2019.10 Taxonomic software using the NCBI database (BLAST, NCBI_16S_20210518). Classified OTUs were determined at ≥ 98% sequence homology and converted to percentages (relative abundance) to determine the representation of each microbe among treatments. OTUs with a relative abundance of less than 0.001% were excluded. Diversity alpha (Shannon, Simpson, and Chao1), sparsity analysis, and beta diversity index (Bray–Curtis similarity) were evaluated using Qiime2 software (version 2019.10).

Statistical Analysis

Discriminant function analysis (DFA) was used to examine whether antibiotic concentrations and the abundance of resistance genes differ significantly between the two sites studied (UPS and DWS). This analysis was performed using XLSTAT. Non-metric multidimensional scaling (NMDS), Goodness-of-fit test, Mantel test, Adonis test and Procrustes analysis were conducted based on Bray–Curtis dissimilarity matrices on RStudio (version 4.0.5) using package “Vegan” for ecological statistical analysis. Based on the NMDS results, Procrustes analysis was performed using the protest function with 999 permutations to verify the significance. To assess the effect of different ATB concentrations on RGs abundance, the quantitative data were converted into binary data. To do this, we compared the obtained concentrations with the minimal selective concentration (MSC) and the minimal inhibitory concentration (MIC), and according to whether this concentration corresponded to one or another of the mutant selection window (MSW) hypothesis ranges, this value was converted to 1, or 0. Thus, three data sets were prepared that separated samples with concentrations less than MSC, greater than MIC and those that fell within MSW. These data were then used to construct a network analysis. The network analysis between ATBs and RGs was based on the dissimilarity coefficient calculated from the binary data using XLSAT. Similarly, Pearson’s dissimilarity coefficient from the quantitative data was calculated between the most abundant classes of bacteria and RGs. The network visualization, as well as the calculation of the degree of nodes and edges betweenness, were performed in Cytoscape 3.8.2 software, where only the dissimilarity coefficients > 0.7 were displayed.

Results and Discussion

Temporal Dynamics of Antibiotics, Resistome and Microbiome in the Clain River Biofilms

The results of the 1-year survey showed the presence of ATBs from different families in the biofilms. Between 19 and 20 antibiotics (out of 20 analyzed) were detected in Clain river biofilms collected upstream (UPS) and downstream (DWS) of the WWTP, respectively (Figure 1A and Supplementary Figure 2). Only CPR was not detected at the UPS samples. The values obtained show that the total median ATB concentrations DWS exceed the corresponding concentrations measured UPS (220 vs. 851 ng/g of biofilm). Statistical comparison of the values by DFA showed that this difference is significant (Wilks’ lambda = 0.000002770 and p-value = 0.005—Table 1) between the biofilms collected from both sites. The difference is mainly due to twelve of the 20 ATBs studied, whose concentrations were found to be significantly different (Supplementary Table 2). The results also revealed differences in contamination depending on the sampling periods at the same site. In DWS biofilms, CLR and SPR concentrations increased in February and March, while AZM was not detected in the same period, although it was the main ATB detected in other sampling months. At UPS, the distribution of ATBs is much more variable each month, with no trend or pattern. In November the concentration of almost all ATBs was significantly higher than during other months. Some previous works had already pointed out the presence of similar concentrations of pharmaceuticals (including AZM, CPR, CLR, LVX, NFX, SMX and TMP) in river biofilms (Aubertheau et al., 2017; Auguet et al., 2017; Carles et al., 2021; Valdés et al., 2021) and also observed variations in the distribution of these compounds over time (Valdés et al., 2021).

FIGURE 1

Dynamic of panel (A) ATBs concentration in ng/g of biofilm, of panel (B) cumulative abundance of the resistome, grouped into gene classes and of panel (C) most abundant bacterial phyla in biofilms collected UPS and DWS of the WWTP. Due to analytical problems results for October (resistome composition) and for September (taxonomic composition) are not included.

TABLE 1

Results of the statistical comparison of the concentration of ATBs and normalized abundance of RGs measured in biofilms collected UPS and DWS from the WWPT.

	Comparison of ATBs concentration	Comparison of normalized abundance of resistance genes and MGEs
Wilks’ lambda	0.000002770	0.00000648
F (Observed value)	22,566.915	10,283.206
F (Critical value)	246.464	245.950
p-value	0.005	0.008
Alpha	0.05	0.05

It is well known that ecosystems affected by human activities are generally the most polluted by ARGs (Chu et al., 2018; Zhao et al., 2018; González-Plaza et al., 2019; Buelow et al., 2020; Kunhikannan et al., 2021). Zhuang et al. (2021) conducted a large analysis of scientific publications on the distribution and diversity of ARGs published over the past 30 years (from 1990 to 2020). The results of the nearly 10,000 publications they analyzed showed that the most contaminated habitats for ARGs (after hospitals and farms) were sewage treatment plants and surface water, with little difference between them. The resistome study of UPS samples showed that even without direct exposure to WWTP discharge, river biofilms contain significant amounts of RGs. Anthropogenic expansion means that even areas upstream of cities are likely to be impacted (to a lesser extent) by different discharges (e.g., agricultural, on-site sanitation) that maintain/contribute to the presence of RGs in river biofilms. If the ARGs genetic diversity is close between UPS and DWS, the location directly downstream of discharge clearly maintains a greater potential of selection pressure. Thus, river biofilms are reservoirs with the potential to acquire and maintain resistance to antibiotics, which, however, is characterized by their exposure degree.

The study of the resistome-mobilome showed that the differences in the abundance of resistance genes and MGEs measured in biofilms collected from both sites are statistically significant (Wilks’ lambda = 0.00000648, p-value = 0.008) (Table 1). The unidimensional test of equality showed (Supplementary Table 3) that the normalized abundance of 31 of the 57 RGs and MGEs, which were quantified in more than three samples at both sites, differed significantly between the biofilms collected UPS and DWS. The normalized abundances, shown in Figure 1B and Supplementary Table 7, demonstrate the difference in gene abundance distribution (in %) between the two sites. More tetracycline, beta-lactam and efflux RGs were detected in UPS but fewer MGEs, streptogramin, sulfonamides and polymyxin compared to DWS. The results also show that the gene distributions evolve slightly by periods of a few months over the 1-year monitoring period. Thus, at the UPS site, a significant proportion of streptogramin and trimethoprim RGs is observed from January to March, then their abundance decreases strongly in favor of tetracycline RGs from July to September, but from October to December their proportion increases again. At the DWS site, resistome analysis shows a significant proportion of trimethoprim and tetracycline RGs from January to March, then, from July to September, their abundance decreases in favor of methicillin and polymyxin RGs, after which methicillin abundance decreases, and beta-lactams increase in November and December.

To complete the monitoring of biofilms, the taxonomic composition of biofilms was analyzed in parallel with the above-mentioned analyses. The results show very little variation in taxonomic composition (by Phylum level) between samples collected at the same site for 1 year. At the UPS site, the changes were associated mainly with the number of Cyanobacteria/Chloroplast and Nitrospirae (the number of both decreased in December) and Verrucomicrobia (increased in December). At the DWS site in July and August there was a decrease in the number of Nitrospirae, but an increase in Proteobacteria. The taxonomy diversity also reveals varieties of different species and much larger numbers of Cyanobacteria/Chloroplast and Verrucomicrobia in the UPS-collected samples (Figure 1C). Biofilms collected DWS contain visibly more Nitrospirae (Nitrospira) and Proteobacteria (Alpha, Beta-, Delta- and Gamma- proteobacteria), as well as Bacteroidetes (Bacteroidetes_incertae_sedis, Flavobacteria and Sphingobacteriia) (Supplementary Figure 3). The phylum Proteobacteria consists of Gram-negative bacteria that include a wide variety of pathogens. The abundance of this phylum in WWTPs discharges, in sediments and activated sludge is common (Liang et al., 2009; Atashgahi et al., 2015; Zhang et al., 2019), since Proteobacteria adapted to low-nutrient conditions and are very common in the mammalian gastrointestinal tract (Moon et al., 2018). The presence of Nitrospirae can be explained by the fact that effluents of WWTPs are rich in nitrate, as nitrification is a common process in the biological removal of nitrogen in WWTPs (McCarty, 2018). These results correlate with previous studies that analyzed microbial communities of river bacterioplankton in Australia, plankton with sediments in an urbanized river of Brazil and biofilm collected in drinking water distribution system in China where Cyano-, Proteo -, and Actinobacteria were dominant phyla (Carney et al., 2015; Köchling et al., 2017; Chen et al., 2020). A metagenomic study of artificially grown biofilms in a Chinese urban river showed the same results (Cai et al., 2018).

The taxonomic analysis also pointed out 17 different genera of pathogenic or opportunistic pathogens at both sites, as shown in Supplementary Figure 4. Aeromonas, Clostridium and Pseudomonas spp. were found to be the most abundant potentially pathogenic genera in studied samples. The number of pathogenic bacteria can be particularly significant DWS from WWTPs, which are important sources of pathogen releases since traditional treatment usually has low efficiency for pathogen removal (Aghalari et al., 2020; Mao et al., 2021). In addition, biofilm formation is characteristic of most bacterial pathogens (for about 80%) (Balcázar et al., 2015; Khatoon et al., 2018; Chen et al., 2020).

Incidence of Season and Exposure

A NMDS based on Bray–Curtis dissimilarity matrices was carried out to examine the dissimilarities in ATBs concentration, resistome abundance and microbial community composition among the biofilms collected over 1 year (samples were grouped into three seasons). The resulting ordination plots demonstrate that the differences in microbial taxonomic composition (Figure 2A) and ATBs content (Figure 2B) between the biofilms for both studied sites are statistically significant (confirmed by Goodness-of-fit and Adonis tests). This finding reveals that exposure to WWTP effluent alters taxonomic diversity and increases antibiotic contamination of the biofilm. However, likely, some of the species identified are not living bacteria, but rather remnants of dead bacteria dumped by the WWTP.

FIGURE 2

Bray–Curtis based NMDS plots showing difference between square root transformed data of panel (A) taxonomic composition (Goodness-of-fit: R2 = 0.61, p-value < 0.001; Adonis test: R2 = 0.53, p-value = 0.002; Stress = 0.087), panel (B) ATBs composition (Goodness-of-fit: R2 = 0.83, p-value < 0.001; Adonis test: R2 = 0.68, p-value < 0.001; Stress = 0.045) and panel (C) resistome composition (Goodness-of-fit: R2 = 0.68, p-value < 0.001; Adonis test: R2 = 0.51, p-value < 0.001; Stress = 0.013) between studied biofilms by sampling site and season [samples were divided into three corresponding seasons: summer (June—August), autumn (September—November), and winter (December—February)]. Due to analytical problems results for September (taxonomic composition) and October (resistome composition) are not included. The samples were named by the month of sampling, but they actually reflect the 30 days prior to the sampling date.

Resistome composition also significantly differed by location (Figure 2C). In addition, Goodness-of-fit and Adonis tests revealed statistically significant changes in resistome and taxonomic composition in the biofilms of both sites (UPS and DWS) as a function of the season (p-value = 0.025/0.005 for resistome and p-value = 0.014/0.019 for microbiome collected UPS and DWS, respectively). At the same time, seasonal changes in ATB concentration were observed only in samples collected DWS (p-value = 0.003). Note that the NMDS dissimilarity plots shown in Figure 2 do not give a good indication of the effect of seasons on the distribution of samples, since these are two-dimensional plots.

This finding is consistent with other studies showing changes in river bacterial communities with changes in water temperature (Kaevska et al., 2016; Wang et al., 2021). Several studies have also reported seasonal variations in resistome composition in untreated hospital wastewater in France and Turkey (Buelow et al., 2020; Gönder et al., 2021), in sewage wastewater in Germany (Caucci et al., 2016), in WWTP effluents and river water in Poland (Harnisz et al., 2020). This effect may be related to the concentration of ATBs in the wastewater, which was also reported to be seasonally dependent (Zhang et al., 2015; Gönder et al., 2021). For example, the seasonality of ARGs in wastewater observed by Caucci et al. (2016) was related to the level of antibiotic prescribing in the sampling region. Seasonality in ATB use has already been reported in Europe and the United States, where the highest ATB consumption occurred during the cold period (winter and spring) (Elseviers et al., 2007; Suda et al., 2014).

What Potential Risk Does River Biofilm Contamination Pose?

The permanent coexistence of antibiotics (and probably other co-selectors), microorganisms (including pathogens) and antibiotic resistance determinants (ARGs, MGEs, etc.) makes biofilms potentially favorable micro-environments for the selection, co-selection, transmission and spread of RGs or resistant bacteria. To assess the extent of ATB concentrations in biofilms, the concentrations found were compared with the corresponding predicted no-effect concentrations for the resistance selection (PNEC-MIC) (Figure 3). The comparison showed that the concentrations measured in biofilms could theoretically cause the maintenance and enrichment of existing resistant strains and even lead to the selection of new resistant mutants.

FIGURE 3

A spider chart showing the range (from min to max) of antibiotics concentrations measured in biofilms collected upstream (UPS) and downstream (DWS) from the WWTP compared to PNEC-MIC.

Many studies confirm that anthropogenic pressure contributes to the increase of antibiotic resistance markers (genes and bacteria) in the environment. However, selection and maintenance are not necessarily due to the effect of ATB presence alone. The enrichment of ARGs can be due to cross-resistance effects or to the presence of co-selectors, which are widely present in anthropized ecosystems. Indeed, several studies have proven that MRGs and MGEs can contribute to the maintenance and proliferation of ARGs (Wales and Davies, 2015; Pal et al., 2017; Partridge et al., 2018; Yang et al., 2018; Glibota et al., 2020). Based on this, several options for the enrichment of ARBs in human-exposed environments are possible: (1) selection under high antibiotic pressure or (2) enrichment by MGEs carrying RGs to various antimicrobial agents (including antibiotics) and selective pressure from a cocktail of contaminants (i.e., co-selection).

Some previous works had already pointed out the presence of pharmaceuticals (including AZM, CPR, CLR, LVX, NFX, SMX and TMP) (Aubertheau et al., 2017; Auguet et al., 2017; Carles et al., 2021; Valdés et al., 2021) and MGEs (integrons) in similar samples (Aubertheau et al., 2017; Auguet et al., 2017; Chen et al., 2020; Li and Zhang, 2020). However, it was not clear how the presence of integrons can be linked to the presence of ATBs, as integrons also contain genes other than ARGs or encoding proteins of unknown function.

Relationship Between Resistome and Antibiotics

Procrustes analysis was performed to find out whether the occurrence of RGs is related to the presence of ATBs. This assay was used to see if the total number of RGs was related to the amount of antibiotics measured in biofilms. For this purpose, the results of NMDS ordination based on the dissimilarity matrix were used. The Mantel test was also performed to confirm the results of the previous one. As shown in Figure 4A, the results of the Mantel and Procrustes analyses showed no correlation between the matrices. In contrast, for samples collected DWS (Figure 4B), the resistome dissimilarity matrix has a significant positive relationship with ATBs (M2 = 0.45, r = 0.74, p-value = 0.016 for Procrustes analysis; r = 0.44, p-value = 0.034 for Mantel analysis).

FIGURE 4

NMDS-based procrustean analysis showing correlation between resistome and ATB content (A,B) and resistome and microbial community (C,D) for both studied sites. M2, sum of squares, r, correlation in a symmetric Procrustes rotation, P, p-value (significance).

However, while Procrustes analysis can provide a global view of the overall correlations between ATBs and resistome, it cannot clarify the relationships between individual variables. To examine the relationships between each ATB and its corresponding RG, we assessed the correlation using a network analysis, which can be an effective complement to general correlation information obtained by Procrustes analysis. The patterns of co-occurrence were investigated using network analysis based on strong (r > 0.7) correlations. To better study the effect of different ATBs concentration on AMR development, the idea of the MSW hypothesis was applied (Drlica and Zhao, 2007; Gullberg et al., 2011). For this purpose, ATB concentrations were sorted into three groups: those less than the MSC, those between the MSC and the MIC (which means that they fall within the MSW range), and those greater than the MIC. Accordingly, a biofilm network analysis was performed for each of the three cases (Supplementary Figures 5–7).

To simplify the reading of complex network graphs with multiple connections that overlap each other, the results, we are most interested in, have been summarized in Table 2. A synthesis of the results presented in this table shows that the correlation between ATBs and RGs can be divided into several groups. For the UPS-collected biofilms four groups can be distinguished. The first group (LVX, AZM and ENR) corresponds to the ATBs that have been found to correlate with at least one of the associated RG in each measured concentration range. The second group (NFX, FMQ and CLR) corresponds to the ATBs that are present in two concentration ranges (at < MSC and in MSW) and only when present in MSW they were found to be associated with the relevant RGs. The third group (T-T, RXM, ERY and SMX) corresponds to the ATBs that are quantified in all samples only at concentrations < MSC and were associated with several RGs. The fourth group consists of ATBs such as TMP, OTC and MTZ that have not been correlated with any of the associated RGs, although they have been detected in concentrations consistent with those that may be responsible for resistance selection. Finally, several of the listed ATBs (except SPR, NFX, TMP, OTC and MTZ) were found correlated with a group of the same RGs, without matching logic.

TABLE 2

Summary table of correlation between resistome and ATBs in concentrations corresponding to the three ranges of the MSW hypothesis.

Antibiotic concentration/Related RGs	UPS			DWS
	< MSC	MSC-MIC	> MIC	< MSC	MSC-MIC	> MIC
CPR	–	–	–	–	+	+
	–	–	–	–	qnrA	qnrA
ENR	+	–	+	+	–	+
	qnrA, qnrS	–	qnrB, qnrC, mfsA	qnrA, qnrS, mfsA	–	qnrA, qnrS, mfsA
FMQ	+	+	–	–	+	–
	–	qnrB, qnrC, mfsA	–	–	qnrS	–
LVX	–	+	+	–	+	+
	–	qnrA, qnrS	qnrA, qnrS	–	qnrA	qnrA
NFX	+	+	+	–	+	+
	–	mfsA		–	qnrA	qnrA, mfsA
AZM	–	+	+	–	+	+
	–	ermC, ermX, macB, mefA_10	ermC, mefA_10	–	ermC	ermC
CLR	+	+	–	–	+	+
	–	macB, mefA_10, ermC	–	–	–	–
ERY	+	–	–	–	+	+
	ermX, macB	–	–	–	–	–
RXM	+	–	–	+	–	+
	ermC, ermX, macB, mefA_10	–	–	mefA_10	–	–
SPR	+	+	–	–	+	+
	–	–	–	–	ermC, ermX, ermY, macB	–
T-T	+	–	–	+	–	–
	ermC, ermX, macB, mefA_10	–	–	–	–	–
SMX	+	–	–	+	–	–
	sul1, sulA	–	–	sul1, sulA	–	–
MTZ	+	+	–	–	+	–
	–	–	–	–	–	–
OTC	+	+	–	+	–	–
	–	–	–	tetB, tetQ	–	–
TMP	+	+	+	+	+	+
	–	–	–	–	–	–

The first column lists all antibiotics with already defined concentration ranges for the mutant selection window (MSC and MIC). There are two rows opposite each antibiotic. The top row shows the concentration range in which the antibiotic was measured in the biofilm (highlighted in green), and the bottom row shows whether the network analysis showed a link between that antibiotic and the corresponding RGs, and if so, with which gene.

For the biofilms collected DWS, the distribution of ATBs by groups is slightly different. CPR, NFX, LVX and AZM were found in concentrations falling in two out of three ranges (MSW and > MIC) and in both cases, they were associated with at least one corresponding RG. RXM, SMX, and OTC were found to be associated with the corresponding RGs at concentrations < MSC, while FMQ and SPR at concentrations < MIC. ENR was detected in concentrations either below or above the MSW, but in all cases was associated with several fluoroquinolones’ RGs. The remaining ATBs (CLR, ERY, TMP, MTZ and T-T) were not correlated with any of the corresponding RGs in any case.

The results of the network analysis (DWS) show correlations, but these are more related to co-occurrence than to pressure/selection mechanisms (based on the MSC and MIC). Thus, RXM and SPR correlate with corresponding RGs at concentrations < MSC or in the MSC-MIC range (in the case of SPR), but no longer correlate at concentrations > MIC, where resistant strains generally dominate. Most likely, the explanation can be found in the hypothesis that the genes measured in the DWS-collected biofilms originate from genetic residues captured/trapped by the biofilm matrix. Indeed, DWS biofilms contain MGEs and cell-free DNA that play an important role in horizontal gene transfer (HGT) and the spread of ARGs in bacterial communities (Woegerbauer et al., 2020). In addition, because the MSC and MIC values were derived from assays performed on pure cultures, these values were compared with the predicted no-effect concentrations for resistance selection (PNECR) published by Murray et al. (2020). In this paper, the authors presented data on experimentally derived selective effect concentrations of several antibiotics (five of the 20 analyzed in our study) from studies of reduced growth of a complex bacterial community in wastewater and changes in the abundance of class 1 integrons. For two antibiotics (TMP and CPR), the MSC and MIC values we used were lower than the PNECR values published by Murray et al. (2020) (Supplementary Table 4). For the other antibiotics (ERY, CLR, and AZM), the PNECR values are much higher than the MIC. However, when looking at Table 2, for two antibiotics (ERY, CLR) no correlations were found with the corresponding genes at concentrations greater than MIC, while for AZM correlations were observed at concentrations > MSC and > MIC. Thus, based on the literature data (obtained from both pure cultures and multispecies communities), it can be concluded that there is no correlation between the concentrations of various antibiotics theoretically favorable for the selection of ARGs/MGEs in river biofilms, suggesting that the occurrence of ARGs and MGEs is not related to the antibiotic concentrations measured in river biofilms (no matter exposed or not to the WWTP discharges).

The emergence of resistance genes may be related to their persistence, since resistance acquired by co-selection can persist even at low concentrations of ATB (which has already been observed in various freshwater ecosystems). There is even considerable evidence in the literature for the persistence of ARGs in environments without ATB or located far from sources of antibiotic release and therefore containing very low concentrations (Andersson and Hughes, 2011; Proia et al., 2016; Lopatkin et al., 2017).

Relationship Between Resistome and Microbiome

The relationship between resistome and microbiome was also investigated. The results of the Procrustean analysis revealed a significant correlation pattern of resistome with the microbial communities of the biofilms collected DWS of the WWTP (M2 = 0.24, r = 0.87, p-value = 0.003 for Procrustes analysis; r = 0.44, p-value = 0.032 for Mantel analysis), whereas no correlation was observed in samples collected UPS (Figures 4C,D).

The relationship analysis between resistome and microbiome was supplemented by network analysis, as previously done for resistome and ATB concentrations. Nineteen classes of bacterial species, most common in the samples analyzed, were selected for this analysis (Supplementary Figure 8 and Table 3). The results show that 10 classes of bacteria could be potential carriers of 35 and 7 RGs/MGEs in biofilms collected UPS and DWS, respectively. Among these 10 classes, Nitrospira and Bacteroidetes incertae sedis probably carry the greatest number of diverse genes (27 and 14, respectively) in the UPS-collected biofilms. Nitrospira was correlated with almost all groups of RGs analyzed: chloramphenicol RGs (cat and catB3), aminoglycoside RGs [aadB, aph(3′)-III, AAC(3)-Ib], bacitracin RG (bacA_1), beta-lactam RGs (blaPer-1, blaIMP, blaTEM, blaKPC, blaDHA and cblA), macrolide RGs (ermC, ermY, ermX), one multidrug RG (mdtL), MRGs (cusF, copD), RG of quaternary ammonium compounds (QACs) (qacA), methicillin RG (mecA), trimethoprim RG (dfrF), polymyxin RG (arnA), sulfonamide RGs (sul1, sulA), transposase genes (IS6100, ISS1N) and streptogramin RG [vat(A)]. Bacteroidetes incertae sedis were associated with beta-lactam RGs (blaKPC, blaDHA and cblA), chloramphenicol RGs (cat and catB3), methicillin RG (mecA), trimethoprim RG (dfrF), polymyxin RG (arnA), transposase genes (IS6100, ISS1N), class 1 integron integrase (Intl1), ERY RG (ermC), sulfonamide RG (sul1) and streptogramin RG [vat(A)]. Chloroplast was a possible host of two beta-lactam RGs (blaCTX–M, blaIMP), multidrug RGs (mdtL, tolC), one macrolide RG (mefA_10), cooper RG (cusF) and QACs RG (qacA). Finally, Verrucomicrobiae was found to be correlated with class 3 integron integrase (Intl3), vancomycin RG (vanA) and aminoglycoside RG [aac(6′)-IIa], and Acidobacteria_Gp6 with TTV. Other works have already published evidence on the close association between resistome and microbiome in different compartments (wastewater, soil, and human intestinal tract) (Zhao et al., 2018). However, one cannot be certain that the measured abundance of resistance genes is solely due to the taxonomic composition of biofilms. The correlation between certain bacterial species and RGs could mean that the presence of certain genes is due to and/or favored by the presence of particular strains; or that there is co-selection by MGEs carrying genes encoding resistance to more than one type of antimicrobial agent (e.g., disinfectants and antibiotics); or that it is simply a co-occurrence.

TABLE 3

Summary table of correlation between resistome and microbial communities.

Microbial communities		Genes Biofilms
Phylum	Class	UPS	DWS
Acidobacteria	Acidobacteria_Gp6	TTV	aac(6′)-IIa
Actinobacteria	Actinobacteria	–	–
Bacteroidetes	Bacteroidetes_ incertae_sedis	cat, catB3, cblA, blaDHA, blaKPC, ermC, arnA, sul1, mecA, dfrF, vat(A), ISS1N, IS6100, Intl1	–
	Flavobacteriia	–	–
	Sphingobacteriia	–	–
Gemmatimonadetes	Gemmatimonadetes	ermC	–
Proteobacteria	Alphaproteobacteria	–	–
	Betaproteobacteria	–	–
	Deltaproteobacteria	–	–
	Gammaproteobacteria	–	blaVIM, acrA, mdfF, cusF, dfrA27
Verrucomicrobia	Opitutae	–	ISS1N
	Subdivision3	–	–
	Verrucomicrobiae	aac(6′)-IIa, van(A), Intl3	–
Candidatus Saccharibacteria	Candidatus Saccharibacteria	–	aac(6′)-IIa
Chloroplast	Chloroplast	blaCTX–M, blaIMP, mefA_10, tolC, mdtL, cusF, qacA	–
Cyanobacteria	Cyanobacteria	–	arcA, cusF
Firmicutes	Bacilli	–	–
	Clostridia	–	–
Nitrospirae	Nitrospira	cat, catB3, aph(3′)-III, aadB, AAC(3)-Ib, bacA_1, cblA, blaTEM, blaIMP, blaDHA, blaKPC, blaPer-1, ermC, ermY, ermX, mdtL, cusF, copD, qacA, arnA, sul1, sulA, mecA, dfrF, vat(A), ISS1N, IS6100	–

As for the samples collected DWS, only Gammaproteobacteria showed the correlation pattern with several resistance genes such as beta-lactam (blaVIM), cooper RG (cusF), multidrug RGs (arcA and mdfF) and TMP RG (dfrA27). Other bacterial phyla had links to only one or two RGs. The lack of significant (non-random) links between bacteria and antibiotic resistome DWS from the WWTP once again indicates that the measured abundance of RGs is likely to be unrelated to the presence of ATBs.

Bottery et al. (2020) examined in detail AMR in bacterial communities and argued that the evolution of bacterial resistance is not only the result of antibiotic exposure but also their combined effects along with bacterial interactions within the communities. But in this study, a comparison of the resistome and microbiome of river biofilms proved difficult to identify specific bacterial communities associated (hosts) with ARGs. Bottery et al. (2020) claimed that bacteria living in biofilms have more favorable conditions compared to free-living species (planktonic lifestyle) because the sessile lifestyle increases the probability of compensatory mutations, which allows stepwise generation of low-cost resistance mutations. The role of compensatory mutations in the acquisition of a low-cost resistance mutation has also been supported by some other researchers (Andersson and Hughes, 2011; Dunai et al., 2019). Moreover, Gullberg et al. (2011) argued that the selection of low-cost and no-cost resistance mutations is favored by sub-MIC antibiotic concentrations.

On the other hand, it is important to remember that in this study ARGs, MGEs, and MRGs abundance had been evaluated during point measurements Various parameters and processes, such as physical, chemical (sorption/desorption, degradation), or biological (reproduction, aging, biodegradation), can influence and change the number of genes. Some ATBs are known to be able to persist for a very long time in river sediments. Several works show that antibiotics (Hektoen et al., 1995) and ARGs (Muziasari et al., 2014; Devarajan et al., 2015) can persist in deep layers of sediments. Likewise, many factors suggest the long-term persistence of ARGs in vivo, since, in the absence of fitness costs, resistant mutants will persist even when antibiotic use is discontinued. Andersson and Hughes (2011) cited several examples of studies proving that antibiotic resistance persists for many years before returning to its original state prior to ATB exposure. However, results obtained by Dunai et al. (2019), who tested the effectiveness of compensatory evolution in multidrug-resistant E. coli under antibiotic-free conditions, showed that such conditions can lead to reduced resistance to certain but not all ATBs. Nevertheless, a similar test on biofilm communities in addition to single-cell cultures would have been necessary to confirm these observations, since this lifestyle is more common in most natural ecosystems.

Conclusion

It appears that the presence of resistance markers in environmental bacteria is a complex process that cannot be explained by considering ATB concentration as the main determinant. The results obtained in this study raises many rather fundamental questions about why the antibiotic dose/resistance selection mechanism in biofilms is ineffective. The complex matrix of environmental biofilms likely affects the bioavailability and bioaccessibility of ATBs, which may limit their impact on microbial communities. Specific studies are needed to better characterize the interaction of ATB with biofilm components, especially at the scale of the microenvironment within the biofilm.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Author Contributions

OM participated in sampling campaigns, performed analytical analyses and sample preparation for metagenomic studies and resistome analysis, performed statistical analyses, and drafted the manuscript. MG participated in the metagenomic studies. JL and LM conceived, funded and coordinated the study, organized and participated in the sampling campaigns, and reviewed and edited the manuscript. CD and CG conceived and coordinated the study, reviewed and edited the manuscript. All authors contributed to the article and approved the submitted version.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

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.