Distribution and Inﬂuence on the Microbial Ecological Relationship of Antibiotic Resistance Genes in Soil at a Watershed Scale

: Antibiotic resistance genes (ARGs) are ubiquitous in the environment, with previous studies mainly focusing on the terrestrial ecosystem, which is prone to higher antibiotic application. However, the characteristics, distribution pattern, and driving factors of soil ARGs at the macro scale are still unclear. In this study, the soil ARGs, antibiotics, mobile genetic elements (MGEs), soil properties, toxic metals, polycyclic aromatic hydrocarbons (PAHs), and bacterial community in the Taipu River Basin were analyzed to investigate the distribution and dissemination of ARGs at a watershed scale. The results revealed that ARGs were widespread in the soils along the Taipu River, and that ARG proﬁles varied greatly with different types of land use, but showed regional similarities. The characteristics were mainly determined by antibiotic input and the ARG transmission mediated by MGEs. The order of the contribution of environmental factors to ARG distribution was toxic metals > PAHs > soil properties. Toxic metal pollution was coupled with ARGs through MGE mediation, while PAHs and soil properties were most likely to affect the ARG distribution by shifting the bacterial community. The microbial–ecological relationship changed signiﬁcantly with the enrichment of ARGs, and its impact may extend to the watershed scale. Transposon IS1247 can be used as an indicator of the ARGs impact on the microbial ecological relationship in the soils of the Taipu River Basin.


Introduction
Antibiotic resistance genes (ARGs) are emerging pollutants [1]. China is facing serious ARG pollution [2]. Traces of ARGs have been found in the main water bodies of China, including the Yangtze River, the Yellow River, the Haihe River, and the Pearl River [3][4][5][6], indicating the widespread ARG pollution caused by intensive human activities along these rivers at a macro scale. However, while most of the macroscopic studies on the distribution of ARGs were focused on water bodies and sediments [7,8], it is still unknown whether there is a macro scale migration and transformation of ARGs in the riverside soil.
The Taipu River (E 120°29′~121°4.5′, N 31°0′~31°1.67′), named after the connection between Taihu Lake and Huangpu River, has a total length of 57.6 km and passes through 15 districts in Jiangsu Province, Zhejiang Province, and Shanghai City (Figure 1). It is the largest artificial river in the Taihu Lake Basin. The Taipu River Basin belongs to the plain river network, along which there are textile printing and dyeing areas, industrial areas, farmland areas, aquafarms, and residential areas. A total of 15 sets of sampling sites along the Taipu River at 3 km intervals were sampled, with a total of 45 samples collected, with three replicates for each sample. The land use of sites includes two typical antibiotic pollution resources, three land uses that may be contaminated by antibiotics, and four land uses that have not appeared in the study of ARGs, but have a great impact on the local ecological environment. The sampling sites were designated as S1, S2…S15. Detailed information about the sampling sites was given in Supplementary Table S1.  Supplementary  Table S1).
Soil samples were collected from 0-20 cm of topsoil in July 2019, and immediately returned to the laboratory using ice packs. The soil samples were subdivided into two parts, where one of which was used for DNA extraction and stored at −80 °C, and the other was air-dried (soil particle size < 2 mm) and stored at 4 °C for the determination of physical and chemical properties.

Determination of Environmental Variables
The environmental variables include pH, nitrate nitrogen (NO3 -), ammonium nitrogen (NH4 + ), total carbon (TC), total nitrogen (TN), toxic metals, and polycyclic aromatic hydrocarbons (PAHs) of soil. The soil pH was analyzed in 1:2.5 soil/water after 1 h of shaking at 220 rpm. The determination of nitrate nitrogen and ammonium nitrogen was Soil samples were collected from 0-20 cm of topsoil in July 2019, and immediately returned to the laboratory using ice packs. The soil samples were subdivided into two parts, where one of which was used for DNA extraction and stored at −80 • C, and the other was air-dried (soil particle size < 2 mm) and stored at 4 • C for the determination of physical and chemical properties.

Determination of Environmental Variables
The environmental variables include pH, nitrate nitrogen (NO 3 -), ammonium nitrogen (NH 4 + ), total carbon (TC), total nitrogen (TN), toxic metals, and polycyclic aromatic hydrocarbons (PAHs) of soil. The soil pH was analyzed in 1:2.5 soil/water after 1 h of shaking at 220 rpm. The determination of nitrate nitrogen and ammonium nitrogen was performed using a continuous flow analyzer (SEAL, AA3, Germany). The TC and TN were measured by a C/N Analyzer (Vario MAX C/N, Germany) [32]. The soil sample was prepared by solid phase extraction (SPE) after being extracted by the extract (EDTAsodium phosphate buffer-acetonitrile, Mg(NO 3 ) 2 -NH 3 -2 O), and the soil antibiotics were measured through a high performance liquid chromatography-mass spectrometer (HPLC-MS) (3200 QTRAP™ LC/MS/MS system, SCIEX, USA). The details of this are in the Supplementary Information.

DNA Extraction and High-Throughput Quantitative PCR (HT-qPCR)
The soil DNA was extracted using 0.25 g of soil following the instructions of the DNeasy PowerSoil Kit (MOBIO, USA). To ensure the quality, the extracted soil DNA was run on 1% agarose gel electrophoresis, and the DNA concentration was further quantified using Nanodrop2000 (Nano-Drop, USA). The DNA was then diluted with sterile water to 10 ng/uL, and was kept at −20 • C for subsequent experiments. The HT-qPCR was performed through the SmartChip real-time PCR system (WaferGen, USA), which has the capacity to process 5184 nanopore reactions per run. A total of 384 pairs of primers were selected to target 283 common ARGs and 12 mobile genetic elements (MGEs). The MGEs included 8 transposons, 4 common integrons (intI 1, intI 1LC, intI 2, intI 3) and a 16S rRNA gene. The PCR system consisted of nuclease-free PCR grade water (Promega, USA), 1 × LightCycler 480 SYBR Green I Master, 20 mg/mL of bovine serum albumin (BSA), 500 nM of forward and reverse primer, and 20 ng/µL of DNA templates. The enzyme was first activated at 95 • C for 10 min, denatured at 95 • C for 30 s, annealed and extended at 60 • C for 30 s, and cycled 40 times. The melting curve was automatically generated using WaferGen software, and the qPCR results were analyzed by SmartChip qPCR Software. Data including multiple melting peaks and amplification efficiency exceeding the range (1.8-2.2) were removed. Each sample contained three replicates with a cycle threshold (CT) of 31 as the detection limit [33]. The three replicates were only regarded as positive when they were successfully amplified. The relative copy numbers of ARGs and MGEs were calculated on the WaferGen platform according to Equation (1), and normalized to absolute gene copy numbers following Equation (2), so as to minimize the error caused by the 16S rRNA gene abundance difference between the samples. The standardized copy number of ARGs for each bacterial cell was used and calculated as follows: Relative gene copy number = 10(31−CT)/(10/3) Normalized copy number of ARG gene = (Relative ARG gene copy number / Relative 16S rRNA gene copy number) × 4.1 (2) where CT is the threshold, and 4.1 is the average number of 16S rRNA genes in each bacterium (estimated from the ribosome operon copy number database).

Quantitative PCR (qPCR)
The absolute abundance of ARGs was determined using a Roche 480 fluorescence quantitative PCR instrument. For the absolute quantification of ARGs by an external standard method, the standard plasmid with the original concentration of 1.38 × 108 copy/µL was sequentially diluted 10-fold to construct a standard curve. The reaction system (20 uL) included 10 µL of LightCycler 480 SYBR Green I Master Mix, 0.8 µL of 16S rRNA gene primer (10 µM), 2 µL of the DNA template, 0.5 µL of the BSA solution (20 ug/µL), and nuclease-free water. The program of the PCR reaction included pre-denaturation at 95 • C for 5 min, entering the cycle, denaturation at 95 • C for 15 s, annealing at 60 • C for 1 min, extension at 72 • C for 5 s, and circulation 40 times. The dissolution curve was automatically added by the instrument, and nucleic acid-free water was used as a negative control. Each sample of this study was in triplicate.

16S rRNA Gene Amplification and Illumina High-Throughput Sequencing
DNA samples were diluted to 20 ng/µL, and a two-step PCR amplification protocol was applied to amplify the V4-V5 region fragment of the bacterial 16S rRNA gene (515F-907R) in each sample. The first step was amplified in triplicate in a reaction system (50 µL), which included 25 µL of 2 × premix Takara Ex Taq (Takara, Japan), 1 µL and 10 µM of forward and reverse primers, PCR-grade water, and 20 mg/mL of BSA. PCR was performed according to the procedure of the bacterial primer pair (95 • C/5 min, 94 • C/30 s, 58 • C/30 s, 72 • C/30 s, 72 • C/10 min, cycled 35 times). Then, 5 µL of each reaction system was loaded on 1% agarose gel, and the PCR quality was controlled by verifying no band was in the negative control. The next step was the purification of the reaction system: the bacterial amplicon was loaded onto 1.5% agarose gel and run at 120 V for 40 min. Next, the DNA purification kit DP214 (Beijing, China) was used to cut the band into about 400 bp fragments, which were then purified. The DNA concentration was again determined by fluorescence, and 300 ng of DNA from each barcode amplicon was pooled in a gene bank of each microbiome. Then, sequencing was performed on the Illumina HiSeq 2500 platform.
Sequences were pre-processed in USEARCH v10.0 [25] and VSEARCH v2.12.0, then denoised using UNOISE3 into zero-radius Operational Taxonomic Units (ZOTUs). The amplicon readings were then combined at double ends (with at least a 50 bps overlap). The primers and barcodes were trimmed, and the chimeras and filtered mass were removed using USEARCH v10.0 (maximum expected error threshold was 1.0). The average length of the remaining high-quality amplicon reads was 392 bps (bacterial). The VSEARCH v2.12.0 and the Silva/UNITE database were used to classify ZOTUs. According to the bioinformatics perspective, ZOTUs with lower mass numbers (<10 sequence compositions) may represent PCR or sequencing errors, and therefore have been removed from the sample. For the observed species (ZOTU number), the Chao1 index was used to describe the α diversity of each sample, and a box plot was generated to compare the bacterial OTU diversity.

Data Analysis
Significant differences among the samples were performed using SPSS 20.0 at p < 0.05. The principal coordinate analysis (PCoA), distance-based redundancy analysis (db-RDA), and variation partitioning analysis (VPA) were performed using the R software 'vegan' package, and the heat map was drawn using the 'pheatmap' package. The Pearson correlation coefficient was calculated using the 'psych' package. The network analysis was visualized using Cytoscape software, and the drawing of other graphics was made using Origin 2018.

ARG Profile and Microbial Spatial Pattern
The ARGs and MGEs were detected in all sampling sites. The absolute abundance of ARGs and MGEs in the collected samples were ranged from 5.00 × 10ˆ4 to 1.04 × 10ˆ8 copies/g, in which the highest abundance of ARGs were detected in S5, followed by S4, S15, and S6 ( Figure 2a). The absolute abundance of ARGs in petroleum storage, the sludge treatment plant, and the pharmaceutical plant was much higher than that in other land use ( Figure 2b). Aminoglycosides, multidrugs, beta-lactams, and macrolides were the most widely distributed classes of ARGs (Figure 2c). The highest number of ARGs were detected in S4 (99), while the lowest were detected in S5 (30). The number of ARGs detected in some sites was similar to that of nearby sites, but differed greatly from sites of the same land-use type. The numbers of ARGs detected in the aquaculture area were in the order of S4 > S12; in a residential area, the number of detected ARGs in S9 were much higher than that of S7 and S11; and in farmland, the number of detected ARGs in S10 were significantly higher than that of S13 and S14, but the number of detected ARGs in S9 was similar to that in S10. Unlike absolute abundance, the number of ARGs detected in petroleum storage, the sludge treatment plant, and the pharmaceutical plant was not higher than that in other land uses (Figure 2d).
There was a strong positive correlation (Pearson's r = 0.82, p = 0.05) between the absolute abundance of MGEs and ARGs (Figure 2e). The normalized abundance represented the degree of ARG enrichment in soil bacteria (Figure 2f). Multidrug and aminoglycoside ARGs accumulated the most in microbial cells. The normalized abundance of aminoglycoside ARGs were the highest in S6, and the normalized abundance of multidrug and tetracycline ARGs were the highest in S5. similar to that in S10. Unlike absolute abundance, the number of ARGs detected in petroleum storage, the sludge treatment plant, and the pharmaceutical plant was not higher than that in other land uses (Figure 2d).  A total of 100 high-abundance detected ARG subtypes and MGEs were profiled for further study (Figure 3). The sampling site was generally enriched with 32 ARG subtypes and MGEs, including 8 aminoglycoside ARGs, 7 multidrug ARGs, 4 macrolide ARGs, several other ARGs subtypes (tetracyclines, sulfonamides, beta-lactams, fluoroquinolone), and MGEs, of which aac-Via was the most abundant ARG subtype. Among the detected MGEs, IS1247, ISEcp1, intI 1, and trb-C were the most abundant. IS1247 was enriched in all of the samples, with the highest abundance in S4 and the lowest abundance in S8. The abundance of mef (B) and sulAfolP in the petroleum storage area (S5) was higher than that in other samples, but it was much lower for aadA7, aph6ic, spcN, aph4ib, tetG&F, QnrB4, and mdtg. The abundance of ISEcp1 and intI 1 in S5 was higher than other sites, but trb-C was the opposite. Sustainability 2021, 13, x FOR PEER REVIEW 8 of 18  The microbial diversity (chao1 index) in the aquafarm and the petroleum storage depot was lower than in other land uses (Figure 4a). Overall, the alpha diversity of the bacterial community was not significantly different whether it was in different land uses or regions (Figure 4b), which contributed 26.7% and 8.34% of the total variance (Figure 4c,d). Proteobacteria and Acidobacteria were the dominant species in the soil bacterial community. Proteobacteria occupied a higher proportion of the bacterial community in the petroleum storage depot than in other land uses (Figure 4c), while Acidobacteria occupied a smaller proportion. Compared with Proteobacteria, the microbial composition of other phyla had more differences in different land uses or regions (Figure 4d). The microbial diversity (chao1 index) in the aquafarm and the petroleum storage depot was lower than in other land uses (Figure 4a). Overall, the alpha diversity of the bacterial community was not significantly different whether it was in different land uses or regions (Figure 4b), which contributed 26.7% and 8.34% of the total variance (Figure 4c,d). Proteobacteria and Acidobacteria were the dominant species in the soil bacterial community. Proteobacteria occupied a higher proportion of the bacterial community in the petroleum storage depot than in other land uses (Figure 4c), while Acidobacteria occupied a smaller proportion. Compared with Proteobacteria, the microbial composition of other phyla had more differences in different land uses or regions (Figure 4d).

The Contribution of Environmental Factors on ARGs and Microbial Community
Redundancy analysis (RDA) and variance partitioning canonical correspondence analysis (VPA) was applied to show the contribution of soil chemical properties (pH, TC, TN, ammonia, nitrate), toxic metals (Be, Cr, Co, Ni, Cu, Zn, As, Pb, Cd, Sb, Ba, Bi), PAHs on ARGs, and bacterial community distribution with different land use ( Figure 5). The vectors representing toxic metals and PAHs had similar directions with the first (interpretation: 42.05%) and second axes (interpretation: 13.29%), which affected the distribution of ARGs, respectively. Soil chemical properties, toxic metals, and PAHs accounted for 36%, 50%, and 36% of the total variance. The common interpretation of these factors was 27%. TN, ammonia, nitrate), toxic metals (Be, Cr, Co, Ni, Cu, Zn, As, Pb, Cd, Sb, Ba, Bi), PAHs on ARGs, and bacterial community distribution with different land use ( Figure 5). The vectors representing toxic metals and PAHs had similar directions with the first (interpretation: 42.05%) and second axes (interpretation: 13.29%), which affected the distribution of ARGs, respectively. Soil chemical properties, toxic metals, and PAHs accounted for 36%, 50%, and 36% of the total variance. The common interpretation of these factors was 27%. The contribution of soil chemical properties, toxic metals, and PAHs to the bacterial community was far less than the ARG distribution. The soil chemical properties independently explained 7% of the total variance, 3% of the toxic metals, and 8% of the PAHs. The toxic metals and PAHs commonly contributed 6% of the total variance. Antibiotics and ARGs interpreted 3% and 9% of the bacterial community variance in different land use.

Relationship among ARGs, MGEs, and Microbial Community
Network analysis was performed to analyze the relationships between MGEs and ARG subtypes (Figure 6), which included 6 soil microbial phyla, 100 high-abundance ARG subtypes, and 20 MGEs (Pearson's r > 0.6, p < 0.01). Three modules were set up in the network. In module I, aminoglycoside ARG subtypes were the main part and MGEs were at the edge, while MGEs played a different role in Module II and III as hubs. Module II consisted of multiple ARGs, with ISEfm1, IS256, and IS26 closely linked to ARGs and other MGEs as hubs. The scale of Module III was small, but the contributing ARGs were widely distributed in the environment. Among them, intI 1 was significantly correlated with the higher abundance of ARG subtypes (qnrD, mepA, ttgA, and ttgB). Different from the coupling relationship between ARGs and MGEs, the main phyla of bacterial communities formed a relatively independent module in the co-occurrence network, and had no obvious association with ARGs and MGEs ( Figure S1). The co-occurrence network of ARGs, MGEs, and bacterial communities in distinct land use showed that various ARG enrichment degrees had different effects on the bacterial community structure (Figure 7). The average degree of the co-occurrence network of ARGs, MGEs, and bacteria had a positive correlation with the absolute abundance of ARGs (R = 0.88, p = 0.02) (Figure 8). The co-occurrence network of ARGs, MGEs, and bacterial communities in distinct land use showed that various ARG enrichment degrees had different effects on the bacterial community structure (Figure 7). The average degree of the co-occurrence network of ARGs, MGEs, and bacteria had a positive correlation with the absolute abundance of ARGs (R = 0.88, p = 0.02) (Figure 8).

The Characteristics of ARG Distribution
The results showed that ARGs and MGEs were widely enriched in aquafarms, pharmaceutical factories, and sewage treatment with a massive application of antibiotics [34][35][36]. In this study, the abundance of ARGs and MGEs within these land uses were higher than in other areas (except petroleum storage areas) (Figure 2). In addition, the distribution characteristics of ARGs in petroleum storage areas were extremely different from other sites, which further showed the important impact of land uses on ARG distribution.
Although the land-use types had a certain influence on the distribution of ARGs, the regional features were more obvious within the scope of the watershed scale. There was little difference in the number of ARGs detected in different land uses (Figure 2d). In the case of the same land use type, the distribution of ARGs and MGEs at various locations were quite dissimilar. Moreover, no apparent regional variations were found in the microbial distribution at the order level ( Figure S2). The differences in the changes of ARGs and bacterial communities under the spatial geography gradient were compared. The spatial geography difference had a higher contribution to the total variation, and was more obvious in the ARG distribution ( Figure S3). To some extent, the regional correlation represented the ARG transmission at the macro scale, which was also supported by the strong positive correlation between the abundance of ARGs and MGEs in the soil of the Taipu River Basin.

Driving Factors of ARG Distribution
Antibiotics and MGEs were the drivers of ARG distribution, as documented in previous studies [37,38]. The correlation analysis showed the role of two driving factors in the spatial pattern of ARGs. A typical aminoglycoside ARG subtype (aac3Via) was identified with the highest absolute abundance, while no significant correlation was observed with MGEs. The same was true for aph4ib (the co-occurrence hub of ARGs and MGEs), which implied that MGEs were not the main driving force for the enrichment of aminoglycoside ARGs in the environment. In addition, the correlation analysis showed a positive correlation between MGEs and ARGs, except aminoglycoside ( Figure S4). Aminoglycoside ARGs revealed a significant positive correlation with antibiotics, but not MGEs, indicating that the enrichment of aminoglycoside ARGs were related to the antibiotic application. In contrast, multidrug ARGs and other ARGs were linked with MGEs. Some ARG subtypes were even specifically related to MGEs. For example, tet44 was only significantly associated with IS1133, while IS3 was only related to erm (36). That is to say, the distribution of aminoglycoside ARGs was related to the input of antibiotics, and the transmission of ARGs mediated by MGE was the main factor affecting the distribution pattern of other ARGs.

The Effects of Environmental Factors on ARG Distribution
Toxic metals were the key factors affecting the distribution of ARGs and MGEs [39]. The abundance of multidrug ARGs was significantly related to Cr, Co, and As ( Figure S5). A significant positive correlation was found between macrolide ARGs and Cr, Co, Pb, and Bi. In addition, tetracycline, sulfonamide, and rifamycin ARGs were positively correlated with one or several toxic metals. CzcA, a typical HRG, was confirmed to have a positive coupling relationship with ARGs and a class 1 integron [22,40]. CzcA was highly enriched in the soil along the Taipu River and was strongly related to intI 1 at the watershed scale, indicating the compound pollution risk of toxic metals and ARGs in the Taipu River Basin. On the contrary, the influence of soil properties and PAHs on ARG distribution was generally not observed at the watershed scale. However, PAHs can provide selective pressure for ARG enrichment by selecting microorganisms [41]. The mechanism of soil properties (such as pH) providing selective pressure was similar to that of PAHs [42]. The investigation of the Taipu River Basin revealed soil properties and PAHs contributed more to bacterial communities and less to differences in MGEs than toxic metals (Table S2). Therefore, selective pressure may be the main mechanism for the soil properties and PAHs affecting ARG distribution.

The Effect of ARG Enrichment on Microbial Ecological Relationship
The results suggested that the effect of ARG enrichment on the bacterial community structure did not extend to the watershed scale, but did exist. For land use with a low abundance of ARGs, such as residential areas and farmland, the bacterial community was relatively independent in the co-occurrence network and had little connection with ARGs. However, in land use with a high abundance of ARGs, the bacterial community was not only closely related to ARGs and MGEs, but also divided into multiple modules in a co-occurrence network. To some extent, the differences in networks characteristics in various land uses reflected the changes in the influence of ARG enrichment on the bacterial community. In other words, the structural differences of the network indicated the feasibility of ARG enrichment as a pollutant affecting the microbial ecological relationship. The reason why pollution had not risen to the watershed scale was that the ARGs shown in the co-occurrence network were not highly enriched in the area. However, due to the significant correlation between ARGs and MGEs, the ecological risk caused by ARG migration may appear within the watershed scale.

The Indicator of ARG Transmission and Pollution Risk
The Intl1 has been proven to be a common integron closely related to multidrug resistance [43][44][45][46]. The results and previous studies also showed that the intI 1 played a key role in the distribution of high-abundance ARGs and HRGs [47]. However, the ARGs that were significantly related to intI 1 revealed a non-significant impact on the ecological relationship of the soil bacterial community in the Taipu River Basin. Actually, the dominant MGE will be taken over as selective pressure changes or new genes enter the pool [48]. IS1247 was a common transposase in the co-occurrence network of ARGs and bacterial communities under different land uses, and can spread between different bacteria and pathogens [49,50]. It has been found to have a close correlation with ARG subtypes in water and the guts of insects, fish, and humans [51][52][53]. IS insertions can move regions adjacent to alter the ARG expression, enhance the niche adaptability of bacteria, and affect the structure of the bacterial communities [54][55][56]. Jose, etc., suggested that the co-occurrence network of integrons and microorganisms reflects the specific hosts of integrons, the dominant bacteria, and their relationship under resistance pressure [57]. This study showed that IS1247 was widespread in the Taipu River Basin, and was strongly related to the co-occurrence network properties (average degree, density, average clustering coefficient) of ARGs, MGEs, and the bacterial community ( Figure 8), but not environmental factors, which indicated that IS1247 was a common and stable element that can be used to measure the impact of ARGs on the properties of bacterial co-occurrence networks. In other words, intI 1 can be regarded as an indicator of MGE-mediated ARG transmission under the influence of toxic metals, while IS1247 was more suitable as an indicator of the influence of ARGs on the bacterial community structure.

Conclusions
Our study provides comprehensive insights into the migration, pollutant coupling, and microbial ecological risks of soil ARGs at the macro scale. The ARGs were widespread in the soils of the Taipu River Basin. Due to the selection induced by the antibiotic application and MGE-mediated ARG transmission, the distribution of ARGs formed the spatial characteristics of land use type differences and regional similarities at the watershed scale. MGE-mediated ARG transmission had a greater impact on the complex spatial pattern caused by urbanization at the macro scale. In this process, toxic metals formed a coupling pollution with ARGs through co-selection mediated by MGEs. In contrast, PAHs and soil properties affected the ARG distribution mainly through the shift in bacterial community. In general, the contribution of toxic metal-induced co-selection to the ARG distribution was greater than the contribution of PAHs and soil properties in the Taipu River Basin. The enrichment of ARGs significantly changed the microbial ecological relationship in certain areas of the Taipu River Basin, and these changes might be able to reach to the basin scale. Intl1 and IS1247 could be used as molecular indicators of the potential transmission capacity of ARGs and the influence on the soil bacterial community in the Taipu River Basin, respectively. This study provides new insights into the migration of ARGs in riverside soils and their impact on the ecological relationship of bacteria at a macro scale.