Dynamics of integron structures across a wastewater network – Implications to resistance gene transfer

Highlights • Prevalence and types of integrons vary widely across a wastewater network.• Human-impacted class 1 integrons carrying ARGs dominate, but are highest in hospital wastewater.• New qPCR assays are reported that segregate "anthropogenic" class 1 integrons that carry ARGs (aint1) vs empty structures (eaint1).• Recycled activated sludge has triple the “empty structures” per integron relative to WWTP liquid effluents.• Integron dynamics help identify wastewater compartments with elevated ARG transfer.


Introduction
Increasing antibiotic resistance is an evolutionary response to the overuse of antibiotics and environmental pollution, and is a global health concern (Chen and Zhang, 2013;Devarajan et al., 2016;Quintela-Baluja et al., 2015;Xu et al., 2016). Antibiotic use selects for resistance strains in human and animal wastes, which can be released to the environment via wastewater, spreading antibiotic resistance genes (ARGs) and bacteria across nature. Selection is often through versatile gene transfer mechanisms, which allow relatively unconstrained ARG migration across phylogenetic, ecological, and geographic borders. As such, the use of antibiotics in one habitat (e.g., hospitals) can impact resistomes in other locations across the environment. Therefore, understanding the vectors and transfer mechanisms by which bacteria share and acquire ARGs is important, maybe more important than studying the levels of ARGs themselves.
A key group of ARG transmission vectors in bacteria are integrons that are carried on mobile genetic elements (MGEs), including plasmids and transposons (Cambray et al., 2010). Integrons are genetically mobile bacterial recombination systems (i.e., mobile integrons) that allow the acquisition and expression of promoterless protein-coding sequences (gene cassettes; GC). They can confer rapid adaptation and selective advantage to host bacteria (often under specific environmental pressures). A large array of GCs that carry ARGs have been found in integrons associated with MGEs (Cambray et al., 2010;Partridge et al. 2009), which are very diverse and associated with ARG transmission (Cambray et al., 2011).
Class 1 integrons that carry ARGs are found in both pathogenic and non-pathogenic strains, although not all class 1 integrons are the same. A sub-group of class 1 integrons have been related with "clinical" settings (hereafter called clinical class 1 integrons, "clintI1") (Gillings et al., 2008). clintI1 differs from other class 1 integrons in that its integron-integrase gene is conserved (5'CS end), but it carries an additional segment at its 3'CS end, which is composed of qacEΔ1 (low-level resistance to quaternary ammonium compounds), sul1 (sulfonamide resistance), and orf5 (unknown function) (Gillings et al., 2008). Like other class 1 integrons, clintI1 is co-regulated by SOS response systems, suggesting reordering and acquiring existing/exogenous cassettes may be activated under stressed conditions when genetic adaptation is needed (Barraud and Ploy, 2015;Cambray et al., 2011). However, integron regulation is not fully understood -in particular how it influences the dynamics of ARG transmission and antibiotic resistance spread. Understanding the role of integrons within a system is further confused by the fact that many integrons, including clintI1, can be found "empty"; e.g., carrying no ARGs (Chakraborty et al. 2013;Chang et al., 2007;Dillon et al., 2005;Moura et al., 2012;Sá et al., 2010;Toval et al., 2015). The incorporation of ARGs into clintI1 also occurs in the absence of antibiotics (Pérez-Valdespino et al., 2016), suggesting integron and ARG ecology are not solely driven by antibiotic use.
Here we examine integron dynamics by quantifying different integrons and microbiota across a wastewater network that includes community and hospital sources, sewer lines, a wastewater treatment plant (WWTP), and upstream and downstream receiving waters (Fig. 1). Integron types, states, and possible bacterial hosts (called here "potential bacteria carrying mobile integrons" or PBCMI) were compared across compartments. The goal was to delineate the spatial dynamics of integrons, especially human-impacted class 1 integrons, to help explain where ARG transmission really occurs within wastewater systems. However, to do this, we needed to design new probe-primer sets that better segregate clintI1 from intI1, and probes for quantifying class 1 structures that do not carry ARGs.

Study area and sampling procedure
Sampling was performed during the summer in a city in northern Spain, which has limited heavy industry or agriculture. A schematic of the network is shown in Fig. 1. Samples were collected over three weeks, including effluents from two hospitals (HP_A and HP_B); the parallel community sewer line (CM); the WWTP influent (INF) and effluent (EFF); recycled activated sludge (RAS); and the river water column and sediments 100 m upstream (RU and SRU) and downstream (RD and SRD) of the WWTP discharge. Further details of the wastewater network, sampling regime, and sample handling and processing methods have been reported previously (Quintela-Baluja et al., 2019). Temperature, dissolved oxygen (DO), pH, and specific conductivity were measured in situ using hand-held probes (Mettler Toledo™, FG3 FiveGo™, and Jenway Model 350 pH Meter).

DNA extraction
Biomass was harvested for DNA extraction via vacuum filtration using sterile 0.22-μm membrane disc filters (Millipore, Billerica, MA, USA) or by centrifugation at 12,000 rpm for 30 min. Extraction was performed using the Fast DNA Spin Kit for Soils (MP Biomedicals, USA), according to the manufacturer's instructions. DNA was stored at − 20 • C prior to subsequent analysis.

Design of new primers and probes to quantify alternate clinical class 1 integron structures
The terminology of Gillings et al. (2008) has been used for many years to define class 1 integron recombination platforms. Pre-clinical class 1 integrons are elements that contain a class 1 integrase, which includes a complete tni transposition module (pclintI1). Conversely, "clinical" class 1 integrons (clintI1) contain qacEΔ, sul1, and orf5 in their structure. A random subset of recently published plasmid sequences was examined, and it was determined the old primer set for clintI1 poorly discriminated from pclintI1 structures ( Figure S1, Tables S1-S3). Therefore, we designed a new set of primers that are more specific to the qacEΔ, sul1, and orf5 structure, which we call the aint1 assay.
Our new name represents "anthropogenically impacted class 1 integrons" and is an appropriate name for class 1 integrons that have been impacted by human activity and not just clinical factors. In parallel we also designed a new assay for anthropogenic class 1 integrons that do not carry ARGs, which we call empty structures (eaint1). The design process of the new assays is provided under Results.

Quantitative PCR Assays and related data analysis
IntI1, intI2, and intI3 as well as aint1 and eaint1 were quantified using the qPCR TaqMan method (Table 1). TaqMan qPCR reactions were conducted using SsoAdvanced™ Universal Probes Supermix (BioRad), employing the following thermocycler program: 3 min of initial denaturation at 95 • C, and 40 cycles of 5 sec denaturation at 95 • C and 30 sec annealing/extension (see temperatures in Table 1). In addition, total bacteria and total coliforms were quantified using a SYBR Green based qPCR assay. SYBR Green reactions were conducted using SsoAd-vanced™ Universal SYBR® Green Supermix (BioRad), employing the following thermocycler program: 2 min of initial denaturation at 98 • C, and 40 cycles of 5 sec denaturation at 98 • C and 5 sec annealing/ extension at 60 • C (total bacteria) or 55 • C (total coliforms).
The limit of detection (LOD) and limit of quantification (LOQ) were calculated using a 2-fold dilution series covering the range of 1 to 500 gene copies (gc) from each standard and analysed eight times (R script available at https://github.com/cmerkes/qPCR_LOD_Calc). The LOD was defined as the lowest concentration that can be reliably detected, i. e., at which at least 95% of standard reactions are expected to amplify. The LOQ was defined as the lowest standard that can be quantified with a coefficient of variation less than 0.35 (Klymus et al., 2020) (Table S4).
All assays were carried out in triplicate using the BioRad CFX C1000 System (BioRad, Hercules, CA USA). To minimise inhibition, extracted DNA was diluted to a working concentration of 5 ng/ul and a DNA internal positive control was used in SYBR Green reactions. Briefly, an additional reaction was performed for each sample, containing 10 4 genes copies of the green fluorescent protein gene (gfp). Inhibition was considered to be minimal when the difference in Ct values (ΔCt) of gfp between the spiked sample and the positive control was below 1 Ct. Standard curves for each primer set were constructed using plasmid clones of the target sequences, these were loaded in triplicate in each qPCR run, in parallel with the amplification of test samples. All samples were run in triplicate, using the optimum qPCR conditions determined for each primer set.
Related statistical analyses and data manipulation were performed using the R environment (Team, 2006) with a significant cutoff of α = 0.05. Normality and the variance homoscedasticity were tested by the Shapiro-Wilk test and Levene's test, respectively. Previous conditions were not met for all integron datasets. Therefore, the Krustall-Wallis test was performed to assess statistically significant differences. Conover's-test was performed for pairwise comparison between sampling sites, and adjusted p-values were calculated with the Holm-Bonferroni method.

Target
Primer cut off, checked for chimaeras, and reflective sequences were obtained for each phylotype. Taxonomic classification used QIIME (Caporaso et al., 2010) based on the uCLUST method (Prasad et al., 2015). See Quintela-Baluja et al. (2019) for details. All data analysis and visualisations were conducted using R (Team 2006) through the Rstudio IDE (http://www.rstudio.com/). OTU counts and associated taxonomic assignments were imported and merged into a phyloseq object (McMurdie and Holmes, 2013). All samples were rarefied to ensure the same number of reads per sample (8,577), which corresponded to the sample with the fewest number of sequences, resulting in 5,698 unique OTUs.
Assessment of sample completeness (sample coverage) and alpha diversity were calculated with the R package iNEXT (Hsieh et al., 2016). Alpha-diversity was calculated using Hill numbers (D) with q from 0 to 3 and were expressed in units of "effective" number of species. As q increases, low-abundance OTUs are assigned less weight, while high-abundance OTUs are assigned more weight; this enables the investigation of diversity at different scales. At q=0, D is simply the number of observed OTUs (richness). At q=1, Dexp (Shannon), and at q=2, D is equal to inverse Simpson's index. A cluster dendrogram of community composition dissimilarity (Bray-Curtis, average neighbour clustering) was calculated with the R package vegan (Oksanen et al., 2019).
To determine the potential bacteria carrying mobile integrons (PBCMIs) in each compartment, consensus sequences from class 1 (anthropogenic and non-anthropogenic), 2, and 3 integrons were aligned against the NCBI nr database (database accessed on June 2019) using BLAST® (McGinnis and Madden, 2004). Accession numbers from sequences with 100% identity were downloaded and taxonomy assigned with the R package taxonomizr. Only families with more than one entry in the NCBI nr database were considered as potential integron hosts and used to filter the OTU table from the samples. The PBCMIs (OTU level) distribution across the network was visualized using a two-dimensional hierarchical clustering in conjunction with a heatmap of relative abundances.

Biomarker signature analysis
Linear discriminant analysis effect sizes (LEfSe) was employed to identify PBCMI that are statistically significant across sampling locations (Segata et al., 2011). With the normalised relative abundance matrix, LEfSe uses the Kruskal-Wallis rank-sum test to detect features with significantly different abundances between the assigned taxa and performs Linear Discriminant Analysis to estimate the effect size of each feature. Only taxa with average abundances >1% were considered significant. A significance level alpha of 0.05 and an effect size threshold of 2x were used for all the biomarkers evaluated in this study.

Co-occurrence between integrons and microbial taxa
A correlation matrix was developed by calculating all possible pairwise Spearman's rank correlations among 267 bacterial families, the three integrases (intI1, intI2, and intI3), and the 3'CS conserved region of aint1 present in samples from the study (n = 27). A correlation between two items was considered if the Spearman's correlation coefficient (ρ) was ≥ 0.5 and the p value was ≤ 0.05. To reduce the chances of obtaining false-positive results, the p values were adjusted with a multiple testing correction using the Benjamini-Hochberg method (Benjamini and Hochberg, 1995). The pairwise correlations of the bacterial families and integrases formed their co-occurrence networks. Network analyses were performed in R environment and was further visualised and explored to identify its topological properties (i.e., clustering coefficient, shortest average path length and modularity) in Gephi (Bastian et al., 2009).

New qPCR probe-primer sets for "human-impacted" and "empty" class 1 integrons
Although multiple probe-primers for class 1 integrons and associated gene cassettes already exist Holmes et al., 2003;Stokes et al., 2001), one set has dominated use, i.e., clintI1 (hereafter called the "old set") (Gillings et al., 2015). This old set targets the integron-integrase gene (5'CS end). However, when one compares that target sequence with more recent sequencing data, it becomes clear the region does not clearly discriminate between "clinical" and "pre-clinical" class 1 integrons ( Figure S1) (pre-clinical here meaning the evolutionary precursor of clintI1). However, we identified a conserved 3'CS region (i.e., attc/qacΔE1/sul1) that is unique to clinical (i.e., human-impacted or anthropogenic) class 1 integrons (see Figs. 2 and S2, and Tables S1-S3 for details). Using this region, we designed and validated a new set of probe-primers that better segregate human-impacted class 1 from non-impacted class 1 structures ( Table 1).
The assay design was as follows. Sequences of each individual region were first aligned using MUSCLE in UGENE v.33 to identify optimal target regions for more specific primers ( Figure S2). This revealed selective targets within amplicon-containing sequences at intersecting loci between attc/qacΔE1 and attI/attc (Fig. 2). Using this region, primers and TaqMan probes were designed using OligoArchitect™ (http://www .oligoarchitect.com) with a melting temperature (Tm) of 60 • C. The selected primers and probes were aligned against NCBI using BLAST® ( McGinnis and Madden, 2004) to rule out non-specific binding to other bacterial targets. We call our new probe-primer assay aint1 to avoid confusion with the old set. The new probe-primer set was assessed in terms of historic clinical and non-clinical integron structures (Tables S2  and S3), and shows aint1 better targets the conserved 3'CS segment of the clinical structure (attc/qacΔE1/sul1) and does not overlap with pre-clinical structures. Henceforth we only refer to aint1 hereafter.
Using a similar approach, a new probe-primer assay was also designed to quantify empty anthropogenic class I integron structures (which we call eaint1) ( Table 1). The need to segregate between empty structures and integrons carrying GCs was evident as all of the empty class I integron structures in our plasmid sequence dataset were associated with anthropogenic class 1 structures. The old assay did not uniquely detect anthropogenic class 1 integron structures nor indicate whether they were carrying ARGs in the GC. Combining aint1 and eaint1 allows one to quantify anthropogenic class 1 integron dynamics more exactly, segregating GCs that do not contain ARGs.

Integron pool and distribution across our urban wastewater network
Quantification of integron structures across our wastewater network showed that class 1 integrons were present in all compartments (Table 2 and S5; Fig. 3). In contrast, class 2 integron levels were low everywhere and below detection in river samples collected upstream of the wastewater effluent discharge point or in the downstream water column. Class 3 integrons were also low in upstream samples and lower across the wastewater network, except in hospital wastewater samples. Following quantification of integron gene levels per bacterial cell across the network (i.e., as genome equivalent), hospital effluents contained significantly higher levels of integron genes per cell (p < 0.05), much greater than any other compartment (Fig. 3, Table 2). Metadata on chemical conditions in each compartment is shown in Table S6.
Within the WWTP itself, integron abundances decreased between 93 to 98% from influent to effluent (log removals of 1.18 ± 0.28 [mean ± SD] for intI1, log 1.66 ± 0.26 for intI2, and log 1.21 ± 0.16 for intI3). These removals directly paralleled bacterial removals (Table 2 and S5), implying the removal of bacterial hosts and MGEs occur in tandem across the treatment steps. This is confirmed by intI1, aint1, and intI3 copies per bacterial cell being roughly constant across the system (Table  S7). Despite large reductions in integrons across the WWTP, absolute class 1 and 3 integron gene abundances increased in the river downstream of the WWTP, suggesting integrons were still reaching the environment.
The new probe-primer design showed that the ratio of aint1 to intI1 did vary across "raw" wastewater compartments (Fig. 4A). aint1 was the dominant integron structure downstream of the biological wastewater treatment unit (p <0.05) (Fig. 4B), with relatively low levels of eaint1. In contrast, the percentage of eaint1 in RAS was consistently higher, double that in WWTP influents and hospital samples, and much higher than any other compartment (Fig. 4C-D), except downstream river locations.

Identifying potential bacteria carrying mobile integrons (PBCMI)
Integrase gene abundance across this wastewater network correlates with the type and abundance of bacterial populations, which was previously observed for ARGs within this network (Quintela-Baluja et al., 2019). Therefore, data here lends itself well to statistical analysis on most probable populations carrying key integrons, especially intI1 vs aint1 vs eaint1.
Screening of the NCBI nr database found 28 potential bacterial families known to host intI1: 14 for aint1, 6 for intI2, and 4 for intI3 (Table   S8). Only intI1 has been previously found in Gram-positive bacteria, whereas the 3'CS conserved region of aint1 has only been found in Gram-negative bacteria. Of the identified families with mobile integrons (from the database), 21 were found in one or more samples in this study (Table S8).
Beta-diversity analysis was used to compare microbial community diversities among wastewater network locations. For this analysis, the dataset was re-sampled to ensure the same number of reads per sample (8,577), which corresponded to the sample with the fewest number of sequences, resulting in 5,698 OTUs. The OTUs from PBCMI families were filtered and re-sampled (from the raw reads) to the sample with the fewest number of sequences (1,266), resulting in 392 unique OTUs. The percentage of total estimated species richness in each sample resulted in 95.4 -99.0% for the total microbial community, and 94.4 -98.1% for the PBCMI, which suggests the overall microbial communities were adequately sampled in both cases (Table S9).
Hill number diversity indices were calculated with the parameter q from 0 to 2. Note that as q increases, rare species are given less weight and contribute less toward the "effective number of OTUs" in a sample. Hill numbers with q equal to 0 and 1 showed that raw wastewater samples had lower diversity compared with upstream river samples (both water column and sediment), and from the downstream river -target for the class 1 integron (clintI1, Barraud et al., 2010): -new target for the empty-class 1 integron structure (eaint1; this study); and -new target location for the clinical class 1 integron at the attc/qacΔE1 interface (aint1; this study). sediments (Fig. S3) Specifically, Hill numbers (q equal to 0 and 1) indicated that bacterial diversity was greater in samples upstream of the WWTP (RU, SRU) because of more rare taxa, which was not apparent when only PBCMI taxa were considered (Fig. S4). Conversely, the Bray-Curtis dissimilarity dendrogram showed that total microbial community compositions differed among network compartments, which included four main clusters with distinct ecosystems (cut off = 0.72) (Fig. S5). The first cluster contained samples associated with raw wastewater (community, hospital, and WWTP influent). The second cluster contained RAS, WWTP effluent, and the downstream river water and sediment samples. The third and fourth clusters included the river water column and sediment samples upstream of the WWTP, respectively.
The PBCMI sub-community followed similar patterns, but the cut off value of the cluster structure was significantly lower (cut off = 0.53). This observation, in conjunction with evident two-dimensional hierarchical clustering seen in the heatmap of PBCMI relative abundances (Fig.  S6), suggests that a PBCMI core sub-community exists in all network compartments. However, Fig. S6 shows that WWTP releases clearly impact river microbial communities, increasing PBCMIs in the downstream environment (Figs. 5 and S6). The family Comamonadaceae is ubiquitous, which is often associated with wastewater (e.g., Aeromonadaceae and Enterobacteriaceae), and those more associated with the environment (e.g., Rodobacteriaceae or Sphingomonadaceae).

Microbial communities and their connectivity across the network
Characterising microbial communities across a wastewater network (in terms of diversity, evenness, and taxonomic composition) is key to identifying linkages among compartments and determining microbial contributions from internal versus external sources. LEfSE analysis showed the relative presence of the family Erysipelotrichaceae, Lachnospiraceae, Bacteroidaceae, and Ruminococcaceae tend to best define community wastewater. In contrast, hospital wastewater was characterised by higher abundances of Prevotellaceae, Enterobacteriaceae, and Desulfovibrionaceae. Weeksellaceae and Moraxellaceae were associated with WWTP influent, whereas RAS was best defined by Rhodobacteraceae, Frankiaceae, Dietziaceae, Gordoniaceae, and Intrasporangiaceae (Fig. S7).
Each compartment of the network had characteristic, reflective signature families, suggesting local enrichment and-or selection. To examine this, we combined the biomarker discovery tool LefSE and cooccurrence network analysis to assess the relative transmission of ARGs associated with mobile integrons across compartments. The cocorrelation network consists of 259 nodes (255 taxa and 4 mobile integrons) and 4.465 edges with an average degree or node connectivity of 34.479. The average network distance between all pairs of nodes (average path length) was 2.94 edges with a network diameter of 7 edges.
Bacterial families, comprising PBCMI, that most explain differences between network compartments (LefSE analysis) are highlighted in Fig.  S7. Network analysis produced four different modules (Fig. S8), equivalent to the microbial composition dissimilarity clusters. The network showed possible ARG dissemination pathways via mobile integrons from human-associated taxa (Module I) to bacteria in human-impacted environments (Module II) and non-human environmental bacteria (Modules III and IV). LEfSE analysis revealed that Module I include bacterial families defining raw wastewater (CM, HP, and INF) and all mobile integrons (Fig. S9A). Conversely, Module II mostly contains RAS and river sediment taxa with limited PBCMI contributions (RAS, RD, SRD, and EFF), emphasising the connectivity of RAS bacteria with river sediments (Fig. S8B). Module III contains taxa defining the river water column (Fig. S8C). Finally, Module IV contains mostly bacterial families in upstream river sediments (SRU). However, some bacteria families from the downstream river sediments also are present in Module IV (Fig.  S8D).

Discussion
Here we show that intI1 and aint1 dominated all compartments in our wastewater network, with intI2 and intI3 being much less abundant at all sampling points. All integron markers, except for empty structures (eaint1), were highest in hospital wastewater where antibiotic use is more intensive (Escudero-Oñate et al., 2017;Rodriguez-Mozaz et al., 2015). This is consistent with previous data that showed the relative number of ARGs per bacterial cell were at least one order of magnitude higher in hospital wastewater compared with the other compartments (Quintela-Baluja et al., 2019). Interestingly, while intI2 levels per bacterial cell declined along the wastewater process stream, the relative levels of intI1, aint1, and intI3 per cell did not change amongst Fig. 3. Comparison of the relative abundance of A) class 1; B) class 2; and C) class 3 integron genes per bacterial cell across the network. Krustall-Wallis rank sum test and the Wilcoxon-Mann-Whitney test was used to determine if each integron gene was enriched or reduced at each sample site compared with the population mean (dashed line). Significance: *0.05 and 'ns' not significant. Sample sites: RU, upstream water column; SRU, upstream sediments; HP_A and HP_B, hospital wastewater; CM, community wastewater; INF, WWTP influent; RAS, recycled activated sludge; EFF, WWTP effluent; RD, downstream river water column; and SRD, downstream sediment.
Although the scale of differences may be specific to our network, we suggest hospital wastewaters tend to select and retain a greater array of "defence" genes in their aint1 due to greater stressors in their environment (e.g., therapeutic drugs) and-or richer organic carbon levels (wastes dominated by faeces from patients). Similar uniqueness of hospital wastewater was also seen in Sweden (Kraupner et al., 2021). Conversely, bacteria in RAS grow under nutrient-limited conditions and under longer residence times (e.g., than a gut environment); thus, RAS bacteria are under persistent nutrient stress compared to the gut. Most activated sludge WWTPs recycle RAS for days rather than hours under comparatively dilute but fluctuating biotic and abiotic conditions. Such conditions might promote rearrangement in aint1 and loss of ARGs, but apparently not reacquisition, possibly because selective agents (e.g., antibiotics) are more dilute or not bioavailable.
The high prevalence of empty cassettes in RAS and their apparent low presence in the WWTP effluents has important implications to pathways of integron-mediated ARG spread from WWTPs to the environment. Until now, only two studies have used high-throughput sequencing to study intI1 cassette diversity in environmental samples Relative abundance of potential bacteria carrying integron genes in wastewater network compartments. The right side of the table shows the bacterial families which host integron genes (dark blue). The numbers show the relative abundance of the potential bacteria hosting integrons in each compartment. RU, upstream river water column; SRU, upstream river sediments; HP, hospital wastewater (HP_A and HP_B); CM, community wastewater; INF, WWTP influent; RAS, recycled activated sludge; EFF, WWTP effluent; RD, downstream river water column; and SRD, downstream river sediment. (Gatica et al., 2016;Stalder et al., 2014). While Stalder et al. (2014) sequenced clone libraries of intI1 cassettes from wastewater compartments, Gatica et al. (2016) targeted recombination sites that flank the gene-cassettes of intI1 genes in WWTP effluents. Both studies found a high diversity of gene cassettes; however, many of their detected sequences had unknown functions. Stalder et al. (2014) found their WWTP tended to reduce the diversity of gene cassette arrays relative to raw wastewater and had elevated levels of empty cassettes. Whereas Gatica et al. (2016) showed a marginal persistence of β-lactam resistance gene cassettes in soil microbiomes impacted by WWTP effluents.
Work here shows the nature and content of integron structures, and associated microbial communities differ widely across compartments, further confirming "local" evolutionary ecosystems drive selection in sub-ecosystems embedded within larger networks (Baquero et al., 2008;Quintela-Baluja et al., 2019). Local ecosystems and their communities are spatially separate, and each is driven by different habitat and other selective factors. This includes specific community selection (Munck et al., 2015) and ARG carriage (Quintela-Baluja et al., 2019). Parallel differences exist among prevalent mobile integrons, and their state relative to ARG carriage (empty or not). However, even though different wastewater compartments have different communities, continual hydraulic interconnection allows microbial dispersal and gene flow, which explains linkages among local resistomes, including local mobile integrons.
Co-occurrence analysis ( Fig.s S8 and S9) highlights the association of PBCMIs within microbial communities in raw wastewater, especially compared with communities less impacted by raw wastewater. The importance of this observation is most evident in river sediments, where resident bacteria and their mobile integrons reflect unsettled floc in WWTP effluents. As such, sediments may be place of environmental ARG transfer, although there would need to be selective drivers for gene acquisition to occur, such as other pollutants in the sediments. Our previous work supported this hypothesis where downstream sediments were largely comprised of RAS bacteria (Quintela-Baluja et al., 2019). A recent study by Di Cesare et al. (2020) supported this hypothesis by showing that the vertical distribution of intI1 in lake sediments mirrored the deposition of bacteria from the water column and to sediments below.
Consistent with previous work, sediments downstream of our WWTP were enriched in ARGs, probably due to outputs from the WWTP (Di Cesare et al., 2020;Quintela-Baluja et al., 2019), which is consistent with greater abundances of aint1 compared with upstream sediments. Although it is speculation, we suspect the co-release of potentially selective agents in WWTP effluents might further promote such persistence. Specifically, the dynamic process of gene cassette acquisition in aint1 and the presence of a qac gene in its structure (resistance to biocides) potentially confers advantageous phenotypes that may be behind the fixation and spread of this integron structure .
A noteworthy new observation here, which supports the potential for ARG transfer in sediments downstream of WWTPs, is elevated levels of eaint1. The two compartments with the highest relative levels of empty clinical class 1 integron cassettes are the RAS and downstream sediments, which implies they are hydraulically connected and may be locations of reacquisition of ARGs due to higher levels of empty mobile integrons. A comparison of microbial communities between RAS and downstream sediments (Quintela-Baluja et al., 2019) furthers shows this connection and suggests downstream sediments could be locations of elevated horizontal gene transfer.

Conclusions
Here we show that integrons have their own ecology within a wastewater network and can be used to speculate on the dynamics of ARGs within such systems. However, integron dynamics appear more complex than previously reported, especially the spatial distribution of ARGs, and relationships with human-impacted and empty class 1 integron cassettes. This recognition was only possible by developing a more exact probe-primers for anthropogenic class 1 integrons that better segregate class 1 integron types and determine whether they carry ARGs or not. We are clearly getting closer to understanding what really happens in wastewater networks relative to AMR spread, but new evidence here suggests ARG transmission in wastewater networks may differ than previously beleived. Regardless, we hope this work will guide new and improved WWTPs and their operations in the future.

Data availability
The 16S rRNA sequencing data were submitted to EBI with accession number PRJEB46784 http://www.ebi.ac.uk/ena/data/view/PRJEB3216.

Declaration of Competing Interest
NA.