High-Throughput Screening of Antimicrobial Resistance Genes and Their Association With Class 1 Integrons in Urban Rivers in Japan

Antimicrobial resistance (AMR) is a serious public health concern. Many countries have implemented AMR surveillance programs for humans and animals, but a scheme for monitoring AMR in the environment has not yet been established. Class 1 integrons, which can acquire antimicrobial resistance genes (ARGs) to gene cassettes, were proposed as a candidate to evaluate the anthropogenic impacts on AMR. However, the association between class 1 integrons and ARGs in aquatic environments is less studied and requires further elucidation. This study used high-throughput quantitative polymerase chain reaction (HT-qPCR) to characterize the pollution profiles of ARGs and mobile gene elements (MGEs) in 24 urban rivers in Tokyo and its surrounding area. The abundance of class 1 integron-integrase gene (intI1) and the array of class 1 integron gene cassettes were also determined. In total, 9–53 target genes were detected per sample, and their abundances increased following effluent discharge from wastewater treatment plants. The river and wastewater samples were categorized based on their HT-qPCR profiles, indicating that this method was useful for characterizing the pollution status in aquatic environments. The prevalence of intI1 in the rivers was observed. Some ARGs and MGEs were positively correlated with intI1, indicating that intI1 could be used as a proxy for monitoring these ARGs and MGEs in urban rivers. Long-read sequencing of class 1 integron gene cassettes revealed that one to three ARGs were present in the gene cassettes. Regardless of the sample type, bla GES-24, aadA2, and qacH were dominant in the gene cassettes. The source and spread of class 1 integrons carrying these ARGs in aquatic environments should be further monitored.


INTRODUCTION
The health burden of antimicrobial resistance (AMR) is a crucial issue across the world (O'Neill, 2016). While the overuse or abuse of antimicrobial agents for humans and animals has led to the emergence of antimicrobial resistant bacteria (ARB) and antimicrobial resistance genes (ARGs), the recipient environments could serve as their reservoirs (Nnadozie and Odume, 2019). Therefore, One Health, which is a comprehensive and multisectoral approach to address AMR issues in humans, animals, and the environment, serves as an essential initiative to mitigate the spread of AMR in the society. Many countries have implemented surveillance programs for pathogenic ARB in humans and animals (WHO, 2021a). On the other hand, the dimension of environmental AMR remains unknown (Larsson et al., 2018;Samreen et al., 2021;Zhuang et al., 2021), and the interaction between pathogens and environmental bacteria carrying ARGs could be facilitated by horizontal gene transfer (Martínez, 2019). WHO recently launched the Tricycle project, which is aimed at One Health surveillance by focusing on extended-spectrum beta-lactamase (ESBL)-producing Escherichia coli (Pruden et al., 2021;WHO, 2021b). However, monitoring targets and goals to control AMR in the environment is still challenging owing to the lack of basic information such as ARG abundance and diversity.
Metagenomic analysis and high-throughput quantitative polymerase chain reaction (HT-qPCR) are promising tools for the comprehensive surveillance of ARGs and mobile gene elements (MGEs) in the environment. Metagenomic analysis is a non-target screening method in which no preliminary information of genes is known (Chen et al., 2019b;Hendriksen et al., 2019;Liang et al., 2020;Lira et al., 2020). HT-qPCR, which enables the simultaneous quantification of hundreds of target genes, is generally more sensitive than metagenomic analysis for ARGs and MGEs surveillance (Waseem et al., 2019). Although HT-qPCR provides semiquantitative data, it is still rapid and inexpensive for screening complex AMR profiles in the environment. Therefore, HT-qPCR has been employed to evaluate ARGs and MGEs in various aquatic environments, including rivers (Khan et al., 2019;Lai et al., 2021;Yu et al., 2021), water sources (Han et al., 2020), drinking water (Xu et al., 2016), urban sewage Pärnänen et al., 2019;Quintela-Baluja et al., 2019), and aquaculture systems (Muziasari et al., 2016). HT-qPCR provides information to help determine factors that shape environmental resistome, such as bacterial community, antibiotic concentration (Han et al., 2020), and ARG sources in the environment (Khan et al., 2019).
In urban rivers, wastewater treatment plants (WWTPs) are hotspots that release ARB and ARGs as well as residual antibiotics (Michael et al., 2013;Mao et al., 2015;Guo et al., 2017;Amarasiri et al., 2019). In a previous study, HT-qPCR based on 384 primer sets was used for the comprehensive surveillance of ARGs and MGEs in WWTP influent and effluent in seven European countries, and 289 primer sets showed positive results . HT-qPCR revealed that the total abundances of ARGs in recipient surface water bodies were higher than those at upstream sites, suggesting that WWTP effluent was a major source of ARGs in urban aquatic environments Quintela-Baluja et al., 2019;Lai et al., 2021).
Representative indicators are useful for the efficient and routine monitoring of various ARGs and MGEs. Class 1 integrons have been proposed as an anthropogenic pollution marker for AMR (Amos et al., 2015;Gillings et al., 2015;Zheng et al., 2019;Li et al., 2020). Integrons are bacterial genetic elements that can incorporate multiple exogenous genes, including ARGs, into gene cassettes (GCs) by the site-specific recombination function of integrase (Gillings, 2014). Because the incorporated genes can be expressed through the integronassociated promoter (Collis and Hall, 1995), integron GCs containing ARGs can spread multidrug resistance (Gillings, 2014). Five classes of mobile integrons (class 1-5) are involved in the spread of ARGs as they are frequently associated with transposons and conjugative plasmids (Gillings, 2014). The class 1 integron-integrase gene (intI1) was found to be prevalent in wastewater and river water (Ma et al., 2017). A strong correlation was observed between intI1 and ARGs, such as aminoglycoside resistance genes and sulfonamide resistance genes, in aquatic environments (Gillings, 2014;Ma et al., 2017;Dong et al., 2019;Zheng et al., 2019;Agramont et al., 2020;Nguyen et al., 2021). Previous studies further investigated class 1 integron GCs in wastewater using clone library analysis or next-generation sequencing (Gatica et al., 2016;Ma et al., 2017;An et al., 2018). In wastewater, many ARGs conferring resistance to aminoglycoside, beta-lactam, and trimethoprim were often detected in class 1 integron GCs (Ma et al., 2017;An et al., 2018). As there is a large diversity in the types of GCs (Moura et al., 2009), they can be regarded as fingerprints of AMR in the environment.
In Japan, several studies have reported the presence of antimicrobial-resistant E. coli in rivers (Ham et al., 2012;Urase and Sato, 2016;Gomi et al., 2017;Yamashita et al., 2017;Suzuki et al., 2019;Tsutsui and Urase, 2019). Although some studies investigated specific ARGs in rivers (Nguyen et al., 2019;Liu et al., 2020), to the best of our knowledge, no comprehensive profiles of ARGs and MGEs in urban rivers have been reported in Japan. In the present study, HT-qPCR was performed to determine the prevalence of 67 ARGs and MGEs in 24 urban rivers in Tokyo and surrounding prefectures. Based on the results of HT-qPCR, the applicability of intI1 as a surrogate marker was assessed. The arrays of class 1 integron GCs were determined by using nanopore long-read sequencing to further evaluate the association between class 1 integrons and ARGs in urban rivers.

Sampling
From September 2019 to February 2020, river water samples (n = 30) from 24 rivers in Tokyo and its surrounding prefectures (Kanagawa, Chiba, Saitama, and Ibaraki prefectures) were collected. According to Ministry of Land, Infrastructure, Transport and Tourism in Japan, sewage coverage rates in Tokyo Metropolis, Kanagawa, Chiba, Saitama, and Ibaraki prefectures in 2019 were 99.6%, 96.9%, 75.5%, 81.9%, and 63.0%, respectively. The sampling sites are summarized in Figure 1 and Table 1. They were categorized into three groups based on the rough estimation of the percentages of treated effluent from WWTPs to river flow rates (annual average) at the sampling sites (Table 1). Group A included sampling sites with no WWTPs located upstream. Group B included sampling sites in which the percentages of effluent were estimated to be <10% (1-6%, average: 2%). Group C included sampling sites in which the percentages of effluent could be >10% (11-100%, average: 43%). Additional pollution sources such as the livestock industry and decentralized wastewater treatment facilities were not considered in the grouping. In Tamagawa River, Iruma River, and Arakawa River, water samples were collected from upstream (TM1, TM2, IR1, and AR1) and downstream sites (TM3, TM4, TM5, IR2, and AR2) to evaluate the impact of treated effluent discharged between two sites. Influent (Group INF) and treated effluent (Group EFF) samples (n = 8) were collected from four municipal WWTPs (WWTP A-D) in the region. All plants use a conventional activated sludge process followed by chlorine disinfection to treat domestic wastewater. The treatment capacities of WWTP A-D were 20,000 m 3 /day, 140,000 m 3 /day, 290,000 m 3 /day, and 450,000 m 3 /day, respectively.

Water Quality Analysis
The water temperature, pH, and conductivity were measured onsite using a portable Combo meter (Hanna Instruments, RI, United States). Total coliforms and E. coli were cultured on Chromocult coliform agar (Merck Millipore, MA, United States) at 37°C for 24 h. The water samples were filtered through 0.2-µm mixed cellulose ester membranes (DISMIC-25 AS, Advantec, Japan) and subjected to ammonium analysis. Ammonium concentrations were measured using a salicylate method with a spectrophotometer (TNT830A, DR2800-01B, Hach, CO, United States). Total cell counts (TCC) were determined using a flow cytometer (Accuri C6, BD, NJ, United States) by staining samples with SYBR Green I (Thermo Fisher Scientific, MA, United States) .

DNA Extraction
The water samples (100 ml) were filtered through 0.22-µm Isopore polycarbonate membrane filters (Merk Millipore, MA, United States) to harvest bacteria. The filters were dissolved in phenol:chloroform:isoamyl alcohol solution (25:24:1) (Nippon Gene, Japan) and treated by bead beating using a FastPrep 24 Instrument (MP Biomedicals, CA, United States). DNA extraction was performed by using the FastDNA SPIN Kit for Soil (MP Biomedicals) . DNA concentrations were measured with a Nanodrop ND-1000 spectrophotometer (Thermo Fisher Scientific).

Quantitative Polymerase Chain Reaction
Bacterial 16S rRNA genes, intI1, sul1 (sulfonamide resistance gene), and tetA (tetracycline resistance gene) were quantified using LightCycler 480 SYBR Green I Master (Roche, Switzerland) to compare the results of HT-qPCR, and human-specific crossassembly phage (crAssphage) (Dutilh et al., 2014;Stachler et al., 2019) was quantified using LightCycler 480 Probe Master (Roche). The primers and probes used in the analysis are listed in Supplementary Table S1. PCR was performed in triplicate using a LightCycler 480 II (Roche). The PCR conditions for 16S rRNA genes, intI1, sul1, and tetA comprised of 95°C for 5 min, followed by 45 cycles of 95°C for 10 s, 55°C for 20 s, and 72°C for 20 s (detection). Melting curve analysis was performed by increasing the temperature from 65 to 95°C to check for nonspecific amplification. The PCR conditions for crAssphage comprised of 95°C for 5 min, followed by 45 cycles of 95°C for 10 s, 60°C for 50 s, and 72°C for 1 s (detection). A 10fold dilution series (5.0 × 10 1 to 5.0 × 10 6 gene copies/µl) was prepared for standard curves using an artificially synthesized plasmid containing the target genes. The average PCR amplification efficiencies of 16S rRNA genes, intI1, sul1, tetA, and crAssphage were 99.3%, 93.7%, 98.5%, 92.5%, and 97.8%, respectively.

High-Throughput-Quantitative Polymerase Chain Reaction
The DNA extracts of 38 samples were sent to Resistomap Oy (Helsinki, Finland) for HT-qPCR analysis using a SmartChip Real-time PCR system (TaKaRa Bio, Japan) (Stedtfeld et al., 2018). The target genes, including ARGs, MGEs, and 16S rRNA genes, were analyzed using 68 primer sets validated by Primer Set 2.0 (Stedtfeld et al., 2018) Table S2). The PCR reaction mixture (100 nL) was prepared using 1× SmartChip TB Green Gene Expression Master Mix (TaKara Bio, Japan), nuclease-free PCR-grade water, 300 nM of each primer, and 2 ng/μL DNA template. After initial denaturation at 95°C for 10 min, PCR comprised 40 cycles of 95°C for 30 s and 60°C for 30 s, followed by melting curve analysis for each primer set (Wang et al., 2014). The threshold cycle (C T ) of 27 was selected as the detection limit (Muziasari et al., 2016;Muziasari et al., 2017). Amplicons with nonspecific melting curves and multiple peaks were excluded. The mean C T of three technical replicates in each reaction was used to calculate the ΔC T values (C T of detected gene-C T of 16S rRNA gene). The relative abundances of the detected gene to 16S rRNA gene were estimated using the ΔC T method (Schmittgen and Livak, 2008).

Amplicon Sequencing of Class 1 Integron Gene Cassettes
Class 1 integron GCs were amplified from all DNA extracts of rivers and WWTPs samples using 5′CS (5′-GGCATCCAAGCA GCAAG-3′) and 3′CS (5′-AAGCAGACTTGACCTGA-3′), which are specific to the conserved segments of both ends of class 1 integron GCs (Levesque et al., 1995;Ma et al., 2017;. The thermal conditions of the first PCR were as follows: 25 cycles of 94°C for 30 s, 55°C for 1 min, and 72°C for 2.5 min, with a final extension at 72°C for 10 min. TaKaRa EX Taq Hot Start version (TaKaRa Bio) was used for PCR. The first PCR product was purified using the MinElute PCR Purification Kit (Qiagen, Germany). The second PCR was performed using the same primers with nanopore sequencing adapters (underlined): 5′CS-adp (5′-TTTCTGTTGGTGCTGATATTGCGGCATC CAAGCAGCAAG-3′) and 3′CS-adp (5′-ACTTGCCTGTCG CTCTATCTTCAAGCAGACTTGACCTGA-3′). The thermal conditions of the second PCR were as follows: 15 cycles of 94°C for 30 s, 55°C for 1 min, and 72°C for 2.5 min. A final extension at 72°C for 10 min was added. The second PCR products were purified using the MinElute PCR Purification Kit and checked by electrophoresis on 1.0% agarose gel at 100 V for 20 min.
For multiplex nanopore sequencing, barcoding adapters were attached to the second PCR products using PCR Barcoding Expansion Pack 1-96 (Oxford Nanopore Technologies, United Kingdom) and LongAmp Taq 2× Master Mix (New England BioLabs, MA, United States). The barcoding PCR involved the following steps: initial denaturation at 95°C for 3 min, followed by 15 cycles of 95°C for 15 s, 62°C for 15 s, and 65°C for 5 min. A final extension step at 65°C for 5 min was added. The products were purified using AMPure XP (Beckman Coulter, CA, United States). Finally, equal amounts of the barcoded PCR products were pooled and mixed with DNA CS (Ligation Sequencing Kit 1D, Oxford Nanopore Technologies), NEBNext FFPE DNA Repair Buffer (New England BioLabs), NEBNext FFPE DNA Repair Mix (New England BioLabs), Ultra II End-prep reaction buffer (New England BioLabs), and Ultra II End-prep enzyme mix (New England BioLabs). They were incubated at 20°C for 5 min and at 65°C for 5 min. After purification, adapter ligation was performed using the Ligation Sequencing Kit 1D (SQK-LSK109) (Oxford Nanopore Technologies). The prepared library was loaded onto an FLO-MIN106D flow cell (R9.4.1) (Oxford Nanopore Technologies) and sequenced on a MinION (Oxford Nanopore Technologies). Base-calling and debarcoding were then performed using Guppy (version 5.0.16) (Oxford Nanopore Technologies) with the superaccuracy mode. Reads with shorter than 500 bp of sequence length and lower than Q10 of mean quality were excluded using Filtlong (version 0.2.0) (https://github.com/rrwick/Filtlong) from further analysis. Error correction of the filtered reads was performed using Canu (version 2.1.1) (Koren et al., 2017) with default parameters. ARGs were detected using Staramr (version 0. 7.2) (https://github.com/phac-nml/staramr) with the setting of identity ≥90% and overlap ≥60%. The nucleotide sequence data are available at the DDBJ Sequence Read Archive under the accession number DRA013066.

Microbial Community Analysis
Microbial community structures were analyzed for all samples from rivers and WWTPs by amplicon sequencing of V4 regions of 16S rRNA genes. A primer set of 515F (5′-GTGCCAGCMGCCGCGGTAA-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′) with the adapter sequence was used (Caporaso et al., 2011). Paired-end sequencing analysis was performed on the Miseq platform (Illumina, CA, United States) using MiSeq Reagent Kit v3 kit (2 × 300 bp) at Bioengineering Lab (Japan). Quality filtering was conducted using the FASTX-Toolkit (version 0.0.14) (http:// hannonlab.cshl.edu/fastx_toolkit/) to extract reads, which showed a perfect match with the primer sequences. Chimeric sequences and noise were removed in DADA2 in QIIME 2.0 pipeline (version 2021.4) (Bolyen et al., 2019). Phylogenetic analysis was performed using the q2 feature-classifier plugin of QIIME 2.0 with reference sequences in Greengene (version 13_8) (DeSantis et al., 2006). Operational taxonomic units were defined by a sequence similarity threshold of 97%. The nucleotide sequence data are available at the DDBJ Sequence Read Archive under the accession number DRA013028.

Statistical Analysis
Statistical tests were performed using BellCurve for Excel (version 3.21) (Social Survey Research Information Co., Ltd., Japan). Cluster analysis based on the Ward method and principal coordinate analysis based on the Bray-Curtis dissimilarity index was performed using R (version 4.0.5) with the vegan package. Network analysis was performed based on Spearman's rank correlation coefficient between the relative abundances of the target genes and the genus-level abundances of the microbial community. Spearman's rank correlation coefficient was calculated using the psych package in R. Associations with a correlation coefficient of >0.600 (p < 0.05) were visualized using Gephi (version 0.9.2) (Bastian et al., 2009).

Comparison Between Quantitative Polymerase Chain Reaction and High-Throughput-Quantitative Polymerase Chain Reaction
The ratio of target genes against 16S rRNA gene such as intI1/16S rRNA gene, sul1/16S rRNA gene, and tetA/16S rRNA gene assessed using conventional qPCR was compared with the results of HT-qPCR to evaluate the quantitative performance of HT-qPCR (Supplementary Figure S1). Both data showed significantly positive correlations (Pearson's correlation coefficient r = 0.876 for intI1, r = 0.827 for sul1, and r = 0.613 for tetA; p < 0.05). However, the conventional qPCR results were 6.2, 2.6, and 5.0 times higher than those of HT-qPCR for intI1, sul1, and tetA, respectively.

Profiles of River Water Resistome
The relative abundances of the ARGs and MGEs in the samples are shown in Figure 2. There was a positive correlation between total ARGs/16S rRNA genes and total MGEs/16S rRNA genes (Pearson's correlation coefficient r = 0.840, p < 0.05) (Supplemetary Figure S2). The total relative abundances of the ARGs and MGEs were not significantly different among Groups A, B, and C (Steel-Dwass test, p > 0.05). ARGs conferring resistance to aminoglycoside, multidrug, and sulfonamide were dominant among the analyzed ARGs in the river water samples. No general relationship was noted between the relative gene abundances and group category. The total relative gene abundances of ARGs and MGEs at the downstream sites of Tamagawa River, Iruma River, and Arakawa River were 2.4-4.7 times higher than those at the upstream sites. In particular, the downstream site of Arakawa River (AR2), which was located immediately after the effluent discharge of a large WWTP in Arakawa River, showed relative gene abundances equivalent to those of the influent samples. The relative gene abundances of the influent samples were generally higher than those of the other groups, whereas wastewater treatment reduced the relative gene abundances by 30-62%. The abundances of MLSB and tetracycline resistance genes in the influent samples decreased following treatment. In contrast, the relative abundances of multidrug resistance genes as well as aminoglycoside and sulfonamide resistance genes increased or remained constant in the effluent samples. For MGEs, the samples with higher abundances of ARGs/16S rRNA genes also showed higher abundances of MGEs/16S rRNA genes. Extraordinary higher abundances of MGEs/16S rRNA genes were observed at the sampling site in Yaji River (YJ in Group A) than those at the sites in the other tributaries of Tamagawa River.
The river and wastewater samples were grouped based on the relative abundances of ARGs and MGEs by Cluster analysis, as shown in Figure 3. They were categorized into three major clusters: Clusters 1-3. Cluster 1 included wastewater samples. Sub-clusters 1a and 1b contained effluent and influent samples, respectively. The sampling sites of Sakai River (SI) and downstream of Arakawa River (AR2) in Group C, which were considerably affected by effluent, were also included in subcluster 1a. The river samples in Groups A-C were mixed in Cluster 2. Cluster 3 was separated into sub-clusters 3a, 3b, and 3c. Sub-clusters 3a and 3b included the upstream sites of Tamagawa River (TM1 and TM2 in Group A) and its tributaries, and the downstream sites of Tamagawa River (TM3-TM5 in Group C) and several other rivers were included in sub-cluster 3c.
As compositional changes in ARGs and MGEs were observed between upstream and downstream sites in Tamagawa River, the impact of effluent was also observed in cluster analysis for Iruma River (IR1: Cluster 3a, IR2: Cluster 2a) and Arakawa River (AR1: Cluster 3c, AR2: Cluster 1a). Furthermore, after the entry of WWTP effluent into the river, the emergence or increase of ARGs such as ESBL genes (bla VEB , cfxA, and bla GES ), MLSB genes (ermF and mphA), and tetracycline resistance genes (tetQ, tetX, and tetE), as well as MGEs (Tp614 and Tn3), were observed.

Microbial Community Structures
The microbial community structures were analyzed by amplicon sequencing of 16S rRNA genes. Comamonadaceae, Flectobacillus, and Flavobacterium were dominant in the river water samples. Arcobacter, Acinetobacter, and Bacteroides were dominant in the influent samples (Supplementary Figure S3). The community structures were compared using principal coordinate analysis ( Figure 4). Axis 1 distinguished river water and wastewater samples, whereas axis 2 distinguished influent and effluent samples. The microbial community structures in Tamagawa River (upstream and downstream) and its tributaries were similar. However, the community structures in Ohori River (OH) and Hokota River (HO) in Group A were different from those in the other sites in Group A. The community structures in Sakai River (SI), Tsurumi River (TR), Sumida River (SM), and downstream of Aarakawa River (AR2) in Group C were more related to those in effluent or influent features. Although a clear transition of community structure driven by effluent discharge was observed in Iruma River (IR1 and IR2) and Arakawa River (AR1 and AR2), the difference was not clear in Tamagawa River (TM1-2 and TM3-5).

Characterization of Class 1 Integron GCs in River Samples
Class 1 integron GCs were analyzed by amplicon sequencing by a MinION nanopore sequencer. Amplicon sequencing depths of 12.7-29.1 Mb with average raw read lengths of 512-1,080 bp were obtained (Supplementary Table S4). The lengths of the most GCs were found to range from <500 bp to 2000 bp after quality filtering (Supplementary Figure S6). Then, 571-1,228 contigs ranging from 500 to 4,580 bp detected in each sample were analyzed (Supplementary Table S5). The median contig size ranged from 593 to 1,489 bp (Supplementary Table S5). While 35% of the contigs did not contain ARGs, 57%, 8%, and 1% of the contigs contained one, two, three ARGs, respectively. Overall, 148 GCs carrying ARGs were detected, including 65 GCs with one ARG, 76 with two ARGs, and 7 with three ARGs (Supplementary Table S6). The percentages of GCs without ARGs were >65% in TM1, TM2, AK, KA, MA (upstream of Tamagawa River and its tributaries), and IR1 (upstream of Iruma River) in Group A (Supplementary Figure S7). Conversely, >90% of the contigs in some rivers in Groups A-C and Group INF carried ARGs (Supplementary Figure S7). In the upstream and downstream sites of Tamagawa River and Iruma River, the percentages of GCs with ARGs increased from 35% (TM2) to 63% (TM3) and from 22% (IR1) to 68% (IR2), respectively (Supplementary Figure S7). Figure 7 shows the relative abundances of contigs of class 1 integron GCs containing ARGs. The information of representative class 1 integron GCs containing ARGs is summarized in Table 3. While ARGs encoding resistance to aminoglycoside, beta-lactam, and multidrug were dominant, specific features were also observed in different samples. bla GES-24 (beta-lactam resistance), aadA2, and qacH (multidrug resistance) were prevalent in the GCs in most of the samples. aac(6′)-31 and aadA1 (both aminoglycoside resistance) were more abundant in the samples other than some rivers in Group A. ere(A) (MLSB resistance) was frequently detected in the river samples. At the same time, it was not dominant in Groups INF and EFF. The GC containing two ARGs (aadA2-qacH) was found in rivers and effluent, but it was rare in the influent samples. Tandem array of aac(6′)-IIabla OXA-21 -catB3 (aminoglycoside resistance, beta-lactam resistance, and phenicol resistance) was only detected in influent samples of WWTP A and B.

DISCUSSION
AMR surveillance in wastewater and aquatic environments is required to fill the gap of One Health. The present study was the first to employ HT-qPCR to determine the prevalence of ARGs and MGEs in urban rivers in Japan. Among the genes detected from >70% of the river water samples, sul1 and qacEdelta1 were Frontiers in Environmental Science | www.frontiersin.org February 2022 | Volume 10 | Article 825372 8 also reported as "core wastewater ARGs and MGEs", which were present in all influent and effluent samples of 12 WWTPs in Europe . Moreover, ARGs and MGEs, such as aadA, strB, ermF, intI1, tnpA, and ISPps, were also categorized as "persistent ARGs and MGEs", which remained in >90% of the effluent samples . As these genes were detected in most of the influent and effluent samples in the present study, the core/persistent ARGs and MGEs associated with wastewater could be prevalent in aquatic environments in Japan. Because intI1 and ISPps were even detected at upstream sampling sites of Tamagawa River (TM1) and Iruma River (IR1) with lower human activity impact, these genes could possibly serve as sensitive markers of anthropogenic pollution. Further study is necessary to identify the prevalence and sources of these genes in upstream area.
Effluent from WWTPs affect the resistome in recipient rivers (Rodriguez-Mozaz et al., 2015;Cacace et al., 2019;Khan et al., 2019;Pärnänen et al., 2019;Lai et al., 2021). The significantly positive correlation between intI1 and crAssphage in the studied rivers suggests that class 1 integrons and the associated ARGs and MGEs could be originated from human feces (Chen et al., 2019a;Karkman et al., 2019;Agramont et al., 2020;Nguyen et al., 2021). Clear shifts in the relative abundances and profiles were observed after the entry of treated effluent in Tamagawa River, Iruma River, and Arakawa River, which is consistent with the qPCR results in Tamagawa River . Drastic changes in the microbial community were also observed from upstream to downstream in Iruma River and Arakawa River, suggesting that the microbial community in the effluent could mostly determine the resistome in these rivers. The resistomes in downstream of Arakawa River (AR2) and Sakai River (SI), which had higher percentages of effluent in river flow (AR2: 34%, SI: 52%), were clustered with the effluent resistome. The similarity in resistome profiles between effluent and recipient rivers in urban areas was also demonstrated using HT-qPCR Khan et al., 2019). Conversely, the classification of sampling sites based on the ratio of effluent to river flow rates (Groups A-C) was not always associated with the relative abundance and profiles of  ARGs and MGEs. Even in Group A, exceptionally higher abundances of ARGs and MGEs were observed in Yaji River (YJ), Hokota River (HO), Sakura River (SK), and Ohori River (OH). As categorization of the sites in this study was simple, the actual contribution of effluent at the sampling occasion could be different from the estimated percentages. The performance of WWTPs in removing ARB and ARGs can fluctuate (Harnisz et al., 2020), and additional pollution sources such as livestock industry and decentralized treatment facilities in rural areas should be considered, especially in rivers in rural area. More intensive sampling considering watershed characteristics is necessary to demonstrate the specific resistome profiles in these rivers. It is important to know the hosts of ARGs and MGEs to determine the health risks. While single-cell sorting (Chijiiwa et al., 2020;Wang et al., 2020); emulsion, paired isolation, and concatenation PCR (Hultman et al., 2018); and high-throughput chromosome conformation capture (Stalder et al., 2019) can more directly identify the hosts of ARGs and MGEs, cooccurrence of specific taxa and these genes is also informative to explore potential hosts (Quintela-Baluja et al., 2019;Han et al., 2020;Yu et al., 2021). Interestingly, popular taxa related to some ESKAPE (Enterobacteriaceae, Acinetobacter, Aeromonadaceae, and Pseudomonas) were screened by network analysis. The associations observed in Enterobacteriaceae and Aeromonadaceae (tetA, aadA, and catB3), Acinetobacter (aadA1, aadA2, and strB), and Pseudomonas (aadA6) were endorsed by the comprehensive antibiotic resistance database (Alcock et al., 2020). Although network analysis showed that class 1 integrons and the related genes (intI1, sul1, and qacEdelta1) were associated with miscellaneous taxa, the whole genome database revealed that class 1 integrons are mostly carried by FIGURE 5 | Network analysis of co-occurrence pattern among microbial taxa, ARGs, and MGE based on Spearman's rank correlation coefficient. Associations with Spearman's rank correlation coefficient of >0.600 (p < 0.05) were selected. Different modules are shown in different colors. The size of nodes represents the number of associations, and the width of lines is proportional to Spearman's rank correlation coefficient. three families: Enterobacteriaceae, Pseudomonadaceae, and Moraxellaceae in Gammaproteobacteria . Spearman's rank correlation coefficients between class 1 integrons and these taxa, such as Enterobacteriaceae, Pseudomonas, and Acinetobacter, were only 0.108-0.363, which suggests that other methods should be used to validate the result of network analysis.
Many studies reported that class 1 integrons are a promising indicator of anthropogenic pollution of ARGs Pärnänen et al., 2019;Zheng et al., 2019). The high prevalence and correlation with other ARGs and MGEs in the river samples suggests that intI1 is a representative target in rivers. Some ARGs and MGEs that showed stronger correlation  with intI1 could be genetically associated with class 1 integron GCs. For instance, sul1 and qacEdetlta1 are typically fused genes in the 3′ conserved segment of class 1 integron GCs (Gillings, 2014). The other selected ARGs, such as dfrA27 (Wei et al., 2008), dfrA1 (Zhao et al., 2020), strB (Le-Vo et al., 2019), ereA (Malek et al., 2015), aadA6 (Mirahsani et al., 2016), aadA2 (Ahmed and Shimamoto, 2004), and bla GES (Maehana et al., 2021), were also detected from class 1 integron GCs of Gram-negative bacteria. For MGEs, class 1 integron-dfrA5-IS26 element was found in E. coli (Dawes et al., 2010), and transposition genes such as tnpA were associated with class 1 integrons (Ghaly et al., 2017). These reports were consistent with the results of HT-qPCR, demonstrating that HT-qPCR can dissect the relationship between class 1 integrons and other ARGs/MGEs in aquatic environments.
ARGs conferring resistance to aminoglycoside, beta-lactam, multidrug, MLSB, phenicol, and trimethoprim were frequently acquired in class 1 integron GCs in urban rivers and wastewater samples analyzed in the present study, which is consistent with previous reports (Ma et al., 2017;An et al., 2018;Gatica et al., 2019). The acquisition of specific ARG types by class 1 integrons was demonstrated by whole-genome database analysis . Although various GC types were detected in the present study, major GCs were not found in other studies analyzing class 1 integron GCs in wastewater by amplicon sequencing using the same primer set (Ma et al., 2017;An et al., 2018). As the composition of class 1 GCs was found to be different in river water, sewage, feces, and livestock wastewater (Ma et al., 2017), the diversity of class 1 integron GCs in aquatic environments could likely depend on geographical and socioeconomic settings.
The percentage of GCs containing ARGs was lower in upstream rivers in Group A than that in the other groups. Moreover, the composition of GCs containing ARGs was not necessarily identical among different rivers. Thus, GC profiles in aquatic environments may indicate AMR fingerprints in each watershed. Quantitative monitoring of class1 integrons and qualitative features of its GCs can be integrated for efficient resistome monitoring in aquatic environments. For example, aac (6′)-31 and aadA1 in GCs were more abundant in the samples other than upstream rivers, indicating the impact of wastewater effluent. aadA2-qacH was present in rivers and effluent but not in influent, suggesting that class 1 integron GCs containing aadA2-qacH could be enriched in wastewater treatment. While ere(A) was frequently detected in rivers and wastewater samples, ere(A) acquired by class 1 integrons was more abundant in rivers than influent and effluent samples. This gap suggests that the genetic context of ere(A) could be different in rivers and wastewater.
Common GC types containing bla GES-24 , aadA2, or qacH were observed in river (Groups A-C) and wastewater (Groups INF and EFF) samples. HT-qPCR also revealed a significantly stronger correlation between class 1 integrons and ARGs such as bla GES , aadA2, and qacF/H. Class 1 integron GCs containing aadA2 or qacH were previously reported in a riverine system (Amos et al., 2018). The prevalence of bla GES in class 1 integron GCs in wastewater effluent in Cyprus and Israel was demonstrated, while bla OXA associated with class 1 integron GCs was dominant in effluent samples in other European countries (Gatica et al., 2016). Although the GCs containing aadA2 or qacH were present in the database of integron GCs (INTEGRAL) (Moura et al., 2009), the GC containing bla  has not yet been registered in the database. bla GES-24 encodes a variant of GES, which is class A beta-lactamase. GES-4, -5, -6, -14, which are characterized by a substitution of Gly170Ser, show carbapenem hydrolysis activity (Bontron et al., 2016). As GES-24 has the same substitution, it has potential to hydrolyze carbapenem. bla  was carried by bacteria, such as Aeromonas hydrophila, Citrobacter freundii, Enterobacter cloacae, Klebsiella pneumoniae, and Pseudomonas aeruginosa (Alcock et al., 2020). Four tandem copies of bla GES-24 were detected from class 1 integron GCs on the plasmids of A. hydrophila, which was isolated from clinical wastewater in Japan (Maehana et al., 2021). Although Aeromonas spp. in aquatic environments could be the key host of bla GES-24 , it is also possible that a wide range of bacteria can carry class 1 integron GCs containing bla GES-24 . This ARG was prevalent in rivers in urban areas and wastewater as well as in upstream rivers. More consideration should be given to the dissemination and evolution of bla GES variants associated with class 1 integrons in aquatic environments.
As revealed using conventional qPCR, the relative abundances of the tested genes were 2.six to six times higher than those revealed using HT-qPCR. As the same thermal conditions were employed for all target genes in HT-qPCR, the amplification efficiency may not always be optimized for each gene (Waseem et al., 2019). Although the results of both methods were highly correlated, the quantified values should be carefully interpreted when they are compared with other studies. Moreover, a hydrolysis probe-based HT-qPCR protocol should be compared for more specific quantification (Khan et al., 2019). The technical limitations of class 1 integron GCs analysis include primer coverage of class 1 integron GCs and the accuracy of longread sequencing. The primer set for class 1 integron GCs (5′CS and 3′CS) that were used in this study was also used in other studies on amplicon sequencing of class 1 integron GCs in aquatic environments (Ma et al., 2017;An et al., 2018). However, the coverage of this primer set was 23.6% of 2,153 integrons in the database , which suggests a greater diversity of class 1 integron GCs in the environment. Novel ARGs have been discovered from class 1 integron GCs such as sul4 (sulfonamide resistance) (Razavi et al., 2017) and gar (garosamine-specific aminoglycoside resistance) (Bohm et al., 2020). Therefore, a comprehensive approach such as metagenomic analysis and the amplicon sequencing approach are necessary to reveal the whole picture of class 1 integron GCs in aquatic environments. As co-occurrence of class 2 and 3 integrons with specific ARGs have been reported (Lai et al., 2021), the different integrons could contribute to the spread of specific ARGs in aquatic environments (Gillings, 2014;Deng et al., 2015;An et al., 2018). Regarding longread sequencing, nanopore technology can circumvent the assembly errors of short reads, while the sequencing error rates of long-read sequencing are generally higher (Weirather et al., 2017). Although error correction of raw reads was applied and ARGs with high identity and overlap values in the polished Frontiers in Environmental Science | www.frontiersin.org February 2022 | Volume 10 | Article 825372 reads were explored in this study, the validation by short-read sequencing with higher accuracy could compensate for the limitations of long-read sequencing.

DATA AVAILABILITY STATEMENT
The nucleotide sequence data are available at the DDBJ Sequence Read Archive under the accession numbers DRA013066 and DRA013028.

AUTHOR CONTRIBUTIONS
IK and KN conducted sampling. IK, KN, and MS provided experimental data. IK, KN, and MS analyzed data. IK and KN wrote the main text. All authors corrected and approved the final text.

FUNDING
This work was supported by JSPS KAKENHI (Grant numbers 19K21984 and 21K18742).