Urbanization drives riverine bacterial antibiotic resistome more than taxonomic community at watershed scale

Although the occurrence and distribution of antibiotic resistance genes (ARGs) in various aquatic ecosystems are well explored, understanding of the ecological processes and mechanisms governing the composition and dynamics of bacterial ARGs still remains limited across space and time. Here, we used high-throughput approaches to detect spatial patterns of bacterial ARGs and operational taxonomic units (OTUs) in an urbanizing subtropical watershed, Xiamen, southeast China over a five-year period. At watershed scale, the OTU profiles were undergoing a directional change, but the ARG profiles showed a high stability or stochastic change over time. Compared with the upstream and midstream, the richness, absolute abundance, normalized abundance and diversity of ARGs were significantly higher in the downstream waters. Our results revealed a clear rural-urban disparity in ARG and OTU profiles which were mainly governed by deterministic and stochastic assembly processes, respectively. With the increase of urban building area along the river, the ecological processes of ARG profiles shifted from stochastic to deterministic. In downstream waters, the bacterial ARG profiles were much more stable than bacterial OTUs. Further, our results indicated that both human-dominated environment (e.g., land use) and mobile genetic elements (MGEs) played an important role in shaping the ARG profiles and dynamics. Overall, this was a response to spatially extensive human-landscape interactions that included urban development in the river downstream region, which were common across subtropical coastal cities of China and can alter the ARG profile dynamics along rural-urban gradient. Therefore, watershed management actions aiming at reducing threats posed by ARGs in urbanizing watershed should first consider the surrounding urbanization level and the mode and intensity of human activity. Our findings also imply that due to the decoupling of bacterial function and taxonomy, both aspects should be studied separately.


Introduction
Bacterial resistance to antibiotics is one of the most serious threats to human health, endangering life-saving antibiotic therapies worldwide (Bengtsson-Palme and Larsson, 2015;Neu, 1992;Sommer et al., 2009). Currently, the widespread use of antibiotics causes increasing occurrence, enrichment and spread of antibiotic resistance genes (ARGs) in various ecosystems of the world (Alonso et al., 2001;Smillie et al., 2011;Su et al., 2017;Zhu et al., 2013). From the time when ARGs were identified as emerging environmental contaminants (Pruden et al., 2006;Rysz and Alvarez, 2004), they have received increasing attention from researchers, public and governments (Berendonk et al., 2015; physical gradient changes from the head water to river mouth. The complexity of riverine ecosystems may also result from exchanges between upstream and downstream waters, and between surficial and underground environments (Wohl, 2017). Despite this, most studies on riverine ARGs focused mainly on lotic systems associated with humanpolluted water bodies, such as municipal sewage, hospital sewage and livestock wastewater (Luo et al., 2010;Rodriguez-Mozaz et al., 2015), iterating only lateral connectivity. However, there exists limited knowledge on the properties of the microbial antibiotic resistome in the context of the above four dimensions of river connectivity at watershed scale (Fig. S1).
As it is widely known, human activities are responsible for the distribution and dissemination of ARGs in various environmental systems Pruden et al., 2012). As the urban river environment frequently suffers relatively harsh ecological perturbation caused by human activities, the precise quantification of the risks related to ARGs is extremely difficult . Urbanization arising from human activities alters both biotic and abiotic ecosystem properties (Grimm et al., 2008). Land use change characterizing urbanization is one of the most important ways by which humans change the surrounding physical and chemical environmental factors (e.g. temperature and nutrient load) and collectively cause a disturbance to riverine ecosystem known as the urban stream syndrome via both point and non-point source pollutions (Dodds et al., 2015;Liao et al., 2018). Land use change can substantially degrade water quality in a complex way (Foley et al., 2005). Around the world, rates of land change will increase greatly as urbanization continue to expand and human populations grow (Foley et al., 2005;Grimm et al., 2008). A rapid and large scale urbanization took place in China on an unprecedented scale during recent decades (Zhu et al., 2011), but our understanding about the effect of urbanization, especially land use, on distribution of ARGs in rivers draining China's urbanizing watersheds is still limited. Thus, it is important to conduct research on the ARG dynamics in riverine systems that drain urbanizing watersheds.
Moreover, it is not well known how changes in microbial taxonomic composition relate to ecosystem function at watershed scale Louca et al., 2016). In plant communities, variation in species composition are often coupled to variation in functional trait composition (Robroek et al., 2017). Whereas microbial communities can display opposite results that variable taxonomic composition often coincide stable metabolic functional community . To date, we still know little about how variation in antibiotic resistance function relates to bacterial taxonomic composition at a watershed scale. Network theory is developing in its application across ecology, social science, and economics (Newman, 2006), and network analysis can be used to identify nodes that are of particular importance within a given network by virtue of their central locations and high connectedness at community level (Liu et al., 2019a).
In microbial ecology, understanding the processes driving natural community assembly is a major challenge (Amos et al., 2015;Liu et al., 2018). The study of community assembly processes has established much about the mechanisms which mainly determine the community structure in general, but still remains a major challenge in ecology of ARGs Fang et al., 2019). Basically, there are two main processes, namely stochastic process and/or deterministic process that simultaneously shape microbial communities Mo et al., 2018;Sloan et al., 2006). The understanding of the ecological process of ARGs is crucial for prevention and control of the enrichment and spread of ARGs but which process (stochastic or deterministic) is the main driver remains largely unknown in rivers at watershed scale.
This study aimed at understanding the occurrence, spatial and temporal distribution of ARGs in an urbanizing river at watershed scale. Here, we employed high-throughput approaches to (1) characterize the spatiotemporal variations of antibiotic resistome and bacterial taxonomic communities in the context of river connectivity; (2) discern the relative importance of stochastic and deterministic processes structuring bacterial ARGs and bacterial operational taxonomic units (OTUs); and (3) gain insights into the impact of rural-urban gradient, particularly land use, on the bacterial antibiotic resistome over five years. We hypothesized that the mechanisms underlying the composition of ARGs may be shifted with lateral river connectivity as the surrounding environment is under urbanizing condition; and bacterial resistance functional and taxonomic communities might be shaped by the different assembly mechanisms. F. Peng, et al. Environment International 137 (2020) 105524 2. Materials and methods

Study area and sampling
The Houxi River watershed (24°34′-24°46′N, 117°55′-118°7′E), covering approximately 205 km 2 (Yu et al., 2014), is situated in Xiamen, Fujian province, southeast China (Fig. 1). It has a subtropical monsoon climate and is subject to seasonal changes in hydrology and aquatic environmental conditions (Liu et al., 2013). Houxi River is a typical urbanizing river and has undergone obvious effects of human activity in the downstream region due to the rapid urban expansion of the Xiamen city during the past two decades. Specifically, there has been an increasing urban buildings area in the lower reaches of the Houxi River watershed. It therefore constitutes an ideal system for an in-depth understanding the impact of urbanization on the functional profile of ARGs at watershed scale.
A total of 100 surface water samples were collected along the Houxi River from ten sites in January and July from 2013 to 2017, therefore each site had ten replicate samples ( Fig. 1 and Fig. S2). Site 1 was located at the upstream area with very minor anthropogenic interference/disturbance within the forest catchment (Fig. S2). Sites 2-5 were located at the riverine and lacustrine zones of Shidou and Bantou reservoirs, as midstream area, and sites 6-10 were located in the downstream in a typical urbanizing area and referred to as downstream area in this study (Fig. S2). A volume of 500 mL water samples from each site were collected, and immediately transported to the laboratory. Before filtering, we mixed the water first to homogenize the bacterioplankton community for each sample. Samples were screened with a 200 μm mesh to remove larger particles, and then filtered through a 0.22 μm pore size polycarbonate filters (47 mm diameter, Millipore, Billerica, MA, USA). The filters with bacteria were stored at −80°C until DNA extraction.
The watershed boundaries and the river system were delineated using Digital Elevation Model in ArcGIS 10.4 software (ESRI, Redlands, CA, USA). The watershed land use data in dry seasons of 2013, 2015 and 2017 were extracted from Landsat 8 images which were downloaded from the data sharing infrastructure of earth system science in China (http://www.geodata.cn/). Four main types of land use were used: (1) natural and semi-natural areas, including forest and grassland; (2) agricultural areas, including dry land and paddy fields (although some small scale agriculture also takes place within the forest); (3) built areas, including urban areas and rural settlements; and (4) other areas (mainly barren area), which do not fall in any of the above categories. The percentages of each land use category were calculated in ArcGIS 10.4 software.

DNA extraction and sequencing
The filters were cut into small pieces using flame-disinfected scissors in a sterilized biological safety cabinet (Airtech, Huntington Beach, CA, USA). Then, FastDNA SPIN Kit (MP Biomedicals, Santa Ana, CA, USA) was used to extract the DNA from the filters according to the manufacturer's instructions (Liu et al., 2019b). Spectrophotometric analysis with NanoDrop ND-1000 (Thermo Fisher Scientific, Waltham, MA, USA) was used to determine the concentration and quality of the extracted DNA. Total DNA was stored at −20°C until further use.
Universal primers (341F: 5′-CCTAYGGGRBGCASCAG-3′ and 806R: 5′-GGACTACNNGGGTATCTAAT-3′) were used to amplify V3 and V4 hypervariable regions of the bacterial 16S rRNA gene. PCR reactions were performed in 30 μL in triplicate containing 15 μL of Phusion® High-Fidelity PCR Master Mix (New England Biolabs, Ipswich, MA, USA); 0.2 μL of forward and reverse primers, and 10 ng of template DNA. PCR conditions were comprised of an initial denaturation at 98°C for 1 min, 30 cycles of denaturation at 98°C for 10 s, annealing at 50°C for 30 s, and an elongation phase at 72°C for 60 s followed by a final extension at 72°C for 5 min. After confirming the length of PCR products (about 470 bp) using a 2% agarose gel electrophoresis, the products of the three replicates were mixed in equimoral amounts and purified with GeneJET Gel Extraction Kit (Thermo Scientific, Waltham, MA, USA). Sequencing libraries were generated using NEB Next® Ultra™ DNA Library Prep Kit for Illumina (NEB, Ipswich, MA, USA) following the manufacturer's instructions, and index codes were added. The library quality was assessed using the Qubit@ 2.0 Fluorometer (Thermo Fisher Scientific, Waltham, MA, USA) and Agilent Bioanalyzer 2100 system (Agilent Technologies, Palo Alto, CA, USA). Bar-coded DNA fragments for each sample were sequenced using a paired-end method on Illumina HiSeq platform (Illumina, Inc., San Diego, CA, USA).
These raw sequences were merged using FLASH (Magoč and Salzberg, 2011). Quality control, including barcode and primer sequence removal, was done in QIIME (Caporaso et al., 2011) and the resulting sequence data for each sample were made publically available at the Sequence Read Archive (SRA) database of the NCBI under project number PRJNA383082 and the accession number SRP104354. These sequence data were further processed by assembling high quality sequences (i.e., sequences with maximum number of consecutive lowquality base = 3; minimum of continuous high-quality base = 75% of total read length; maximum number of ambiguous bases = 0) (Caporaso et al., 2011;Liu et al., 2019b). Chimeric sequences were discarded prior to further analysis (Edgar et al., 2011). The UPARSE pipeline was used to pick operational taxonomic units (OTUs) at 97% similarity level (Edgar, 2013). The taxonomy of bacterial OTUs were assigned using the RDP classifier at 80% confidence threshold mapped against the Greengenes database (DeSantis et al., 2006).

High-throughput quantitative PCR (HT-qPCR) and real-time quantitative PCR (qPCR)
Before conducting quantitative PCR, the DNA was diluted to 30 ng/ μL for HT-qPCR and 3.0 ng/μL for qPCR, respectively . The SmartChip Real-time PCR system (Wafergen Biosystems, Fremont, CA, USA) was used to perform HT-qPCR to quantify the abundance of ARGs in each sample. A total of 296 primer pairs targeting 285 ARGs conferring resistance to most major classes of antibiotics, eight transposase genes, two integrase genes and one 16S rRNA gene, and the reactions were used following our previous works Liu et al., 2018). All qPCR reactions were performed in triplicate and each chip included a non-template negative control. Wafergen software automatically generated the melting process. The results of HT-qPCR were analyzed as described previously Looft et al., 2012).
To quantify the absolute abundance of 16S rRNA gene, a SYBR® Green approach was implemented on Lightcycler 480 instrument (Roche, Basel, Switzerland). Each 96-well plate included both standard curve and negative control (Schmittgen and Livak, 2008). Each PCR amplification was set up in 20 μL with triplicate under the program as described in our previous study (Fang et al., 2019).

Statistics
Differences in bacterial functional (i.e. ARGs) and taxonomic (i.e. OTUs) community compositions were determined by Bray-Curtis matrices, which is one of robust approaches in community ecology. Before calculating the Bray-Curtis matrices and Euclidean distance for environmental factors, all data variables were log-transformed except pH to improve data homoscedasticity and normality. Non-metric multidimensional scaling (NMDS) ordination was employed to investigate the distribution of ARG and OTU communities along the levels of urbanization with the PRIMER 7.0 (PRIMER-E, Plymouth, United Kingdom) (Clarke and Gorley, 2015). We investigated differences between ARG communities across longitudinal river sections, five sampling years and two seasons (dry and wet). Significant difference (P) and the degree of separation (Global R) among the above groups were tested by the analysis of similarities (ANOSIM). Global R ranges between 0 and 1, with Global R = 0 indicating no separation and Global R = 1 indicating complete separation.
To evaluate the relative role of stochastic processes on bacterial functional (i.e. ARGs) profile and taxonomic (i.e. OTUs) community assembly, both the Sloan neutral community model (NCM) Sloan et al., 2006) and null model approach based on Raup-Crick coefficients (Chase et al., 2011;Raup and Crick, 1979) were employed. The NCM can analyze 16S rRNA gene and functional gene, especially samples recovered from environment (Sloan et al., 2006). The parameter Nm is an estimate of the dispersal between communities, N is the metacommunity size and m is the immigration rate. The parameter R 2 predicts the overall fit to the model, as positive R 2 indicates the fit to the neutral model and negative R 2 indicates no fit (Sloan et al., 2006). The null model approach based on Raup-Crick dissimilarity coefficients was developed by Chase et al. (2011). This approach creates a re-scaled probability metric ranging from −1 to 1, indicating whether local communities are more dissimilar (approaching 1), as dissimilar (approaching 0), or less dissimilar (approaching −1), than expected by random chance (Raup and Crick, 1979).
To test the effect of dispersal, the linear regression of geographical distance with Bray-Curtis similarity was used (Louca et al., 2016). To assess the temporal pattern of community dynamics, the linear regression of time lags (square root-transformed) with community similarity based on Bray-Curtis similarity was also performed .
To examine the linkages among determinants of urbanization (especially land use types), bacterial taxonomic community and antibiotic resistome, we conducted partial least squares-path modeling (PLS-PM) using the R package "plspm" (Wetzels et al., 2009). To do this, we first ran a permutation test with 999 simulations on all the ecological variables. These variables that significantly influenced the antibiotic resistome were classified into the following four groups: land use group (including built area and other area), physicochemical variable The barcharts showing the richness, abundance and diversity of ARGs (c-e) and OTUs (h-j), respectively. Green, blue and red colors indicate upstream, midstream, and downstream, respectively. Up represents upstream; Mid represents midstream, and Down represents downstream. Significant differences (P < 0.01) between groups are indicated by different letters of the alphabet. Statistical analysis is Mann-Whitney U test. Data are mean ± standard error. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) F. Peng, et al. Environment International 137 (2020) 105524 group (covering dissolved oxygen, turbidity, chlorophyll-a, salinity, TOC, NH 4 -N, NO 3 -N, NO 2 -N, PO 4 -P and TN:TP), taxonomic community (containing different bacterial phyla), and mobile genetic elements (MGEs). These groups were used in the PLS clusters to determine the direct and indirect effects on antibiotic resistome. Subsequently, the overall prediction performance of the model was validated using the goodness-of fit index (GoF) which is defined as the geometric mean of the average R 2 . The higher value indicates better prediction performance.
We used network analysis to identify co-occurrence patterns among ARGs, MGEs and bacterial taxa. Correlation matrices were created through all possible pairwise Spearman's rank correlations. If a Spearman's correlation coefficient (|ρ|) was ≥0.6 and the P-value was < 0.01, the correlation was considered statistically robust (Ma et al., 2017b). To reduce the false-positive results, the P-value was adjusted with a multiple testing correction using the false discovery rate (FDR) method. Network analyses were performed in R environment using the "psych" and "qvalue" packages (Liu et al., 2019a) and visualizations were conducted on the interactive platform of Cytoscape v3.6.1 (available at http://www.cytoscape.org/) and Gephi 0.9.2 (available at https://gephi.org/). We also identified network modules following the Louvain community detection algorithm (Blondel et al., 2008).
Statistical significance tests between groups were the nonparametric Kruskal-Wallis H test and Mann-Whitney U test, and both were performed in SPSS v22.0 (IBM Corp., Armonk, NY, USA).

Changes of bacterial antibiotic resistome and taxonomic communities
We found a significant shift in bacterial ARG profiles and OTU community compositions along the rural-urban gradient, which was well corresponded to different land use types ( Fig. 2; Fig. S3 and Fig.  S4). Rural-urban variability among the ARG profiles greatly overwhelmed their seasonal and inter-annual differences (Table 1 and Fig.  S5). The ANOSIM showed a strong and significant difference (Global R = 0.765, P < 0.01) between the upstream, midstream and downstream waters. The difference of ARGs between January and July was not statistically significant (Global R = 0.015, P = 0.12) (Table 1). Similarly, NMDS results revealed a significant and moderate separation in bacterial taxonomic community along the rural-urban gradient (Global R = 0.507, P < 0.01) (Fig. 2f).
A total of 248 ARGs and 10 MGEs were detected in this study. Further, we identified 106,164 and 238 ARGs,and 16,193,20,525,and 28,614 OTUs in the upstream, midstream and downstream of Houxi River, respectively (Fig. 2b, g). The abundance of both ARGs and OTUs exhibited the highest values in downstream samples (Fig. 2d, i; Table  S2). However, the bacterial functional and taxonomic diversities showed different patterns. The highest ARGs diversity, based on both ARG richness and Shannon-Wiener index, was observed in downstream (Fig. 2c, e). Contrary to bacterial ARGs, the highest OTU richness and diversity were observed in upstream waters (Fig. 2h, j). For ARG richness, there was no significant difference between upstream and midstream waters (Fig. 2c, d and e). But the OTU richness and Shannon-Wiener diversity along the river continuum exhibited a significant difference among each other (Mann-Whitney U test, P < 0.01; Fig. 2h, j).

Shared and unique ARGs among upstream, midstream and downstream
A total of 89 out of 248 ARGs, conferring resistance to all nine types of antibiotics were shared between upstream, midstream and downstream waters, accounting for 84.0%, 54.3% and 37.4% of the total number of ARGs, respectively ( Fig. 2b; Fig. S6a and Table S3). The shared ARGs made up 94.7%, 77.2% and 82.1% of the total ARG absolute abundance for upstream, midstream and downstream, respectively. The majority of shared ARGs dominated the downstream waters, and only one tetracycline resistance gene and three other resistance genes (β-lactam, MLS and vancomycin resistance genes) were abundant in the upstream and midstream, respectively (Fig. S6b).
The unique ARG richness in upstream, midstream and downstream was 1, 9 and 67, respectively. Only 16 ARGs were shared between upstream and downstream sections, while 66 ARGs were shared between midstream and downstream ( Fig. 2b and Fig. S6a).

Fit to the neutral model of ecological processes
The frequency of ARG occurrence along the Houxi River showed no fit to the neutral model (R 2 = −0.204; Fig. 3a; note that negative R 2 values only occur when there is no fit to the model), indicating that neutral processes were not an important mechanism in structuring the ARG profile. In fact, the occurrence frequency of most abundant ARGs that had high relative abundance tended to be more frequent than expected, thus contributing most to the deviation of ARG profile from the neutral model. By contrast, the frequency of the OTU occurrence in the river fitted well to the neutral model (R 2 = 0.789; Fig. 3a), suggesting that stochastic processes could explain a high proportion of the bacterial taxonomic community. Further, the re-scaled probability metric created by the null model approached −1, indicating that the ARG profiles in the downstream waters were less dissimilar than expected by random chance (Fig. 3b). Interestingly, in the upstream and midstream, where human activities were less intense, the re-scaled probability metrics approached 0, indicating that the ARG profiles were as dissimilar as expected by random chance (Fig. 3b). Therefore, the antibiotic resistome in downstream might have greatly contributed to the above observed deviation of the total community from the neutral model (Fig. 3).

Correlation between antibiotic resistome and bacterial taxonomic community
We found distinct spatial patterns between antibiotic resistance function and bacterial taxonomic communities along the Houxi River over five years. The composition of ARGs in the downstream section showed the lowest β-diversity, whereas the taxonomic composition was different and exhibited the highest β-diversity in downstream waters in comparison with those in upstream and midstream waters (Fig. 2a, f;  Fig. S5a). For ARGs, there was no significant relationship through the time-lag analysis (P = 0.949), suggesting that the component of ARGs was relatively stochastic change at the time scale, but for OTUs, the component of community was unstable and undergoing a directional change over time (R = -0.398, P < 0.001; Fig. S7).
The ARG profiles exhibited a higher variability than OTU profiles along the rural-urban gradient ( Fig. S5 and Fig. S8). Further, the ARG profile and OTU community showed a decoupled behavior, however, we did find that there was a weak but significant positive relationship between ARG profile turnover and taxonomic community turnover (r = 0.225, P < 0.001; Fig. S8c). Both antibiotic resistance and taxonomic community compositions were significantly and negatively Higher R value indicates stronger compositional difference between groups (n = 100). ** P < 0.01.
F. Peng, et al. Environment International 137 (2020) 105524 related to geographical distance and environmental factors, except for ARGs in the upstream (P = 0.951) and ARGs with geographical distance in the midstream (P = 0.344), and ARGs exhibited a slightly higher correlation coefficient than OTUs at the watershed scale (Fig.  S7). The network between ARGs and bacterial taxa showed that Acinetobacter belonging to Gamaproteobacteria was the most densely connected node, followed by Megamonas (Firmicutes) and Comamonas (Betaproteobacteria) (Fig. S8b), indicating these bacterial species might be the host candidates of ARG, or these bacterial species and ARGs probably originated from a similar source.

Relationship between environmental variables and antibiotic resistome
Along the rural-urban gradient, most environmental factors, such as electrical conductivity (EC), salinity and nutrients exhibited a significant change (Table S1). Our path modeling analysis indicated that the abundance variation of the 10 MGEs detected had the highest direct effect on antibiotic resistome variation ( Fig. 4 and Table S4). Land use change showed a relatively low direct negative effect on ARGs, but it had the highest indirect impact on the ARGs (Fig. 4b). Meanwhile, bacterial community had a very weak effect on ARGs (Fig. 4), suggesting a decoupling function and taxonomy in the bacterial community.

Co-occurrence patterns among ARGs and MGEs
There were 135 nodes and 565 significant edges among the co-occurrence of ARGs and MGEs in all 100 samples (Table S5). Among them, aminoglycoside resistance genes had the most connections with other ARGs (Fig. S9). ARGs had a low co-occurrence with the 10 MGEs detected, suggesting that these ARGs might have a low potential of horizontal gene transfer rate. However, tetracycline resistance genes were closely connected with the 10 MGEs (Fig. S9), suggesting that in the presence of MGEs, tetracycline resistance genes might be easily transferred across bacterial species.
The numbers of nodes (significant edges) in co-occurrence network of ARGs and MGEs were 114 (652), 98 (186) and 119 (188) in upstream, midstream and downstream, respectively (Fig. S10). Most network properties, including clustering coefficient, network centralization, average degree and network density, were largest in the upstream section indicating that the upstream had the most complex and unique network ( Fig. S10; Table S5). In addition, the modularity indices of the networks from the upstream, midstream and downstream waters were 0.698, 0.649 and 0.656, respectively, suggesting that these networks had a well modular structure property (Table S5).   Peng, et al. Environment International 137 (2020) 105524 4. Discussion In contrast to chemical contaminants, ARGs are capable of selfpropagating and persisting for longer periods of time and are less likely to be rendered harmless by dilution (Pruden et al., 2006). Most of the existing studies, however, describe spatial distributions with limited or short-term monitoring efforts, especially in lotic ecosystems subjected to human activities (Luo et al., 2010;Berendonk et al., 2015;Liu et al., 2018). Data presented here were collected for five years across both the dry (January) and wet (July) seasons, and showed similar spatial patterns. It's generally accepted that the human-dominated environment plays an important role in the enrichment and transmission of ARGs in the urban waterbodies . However, there is a critical knowledge gap about the mechanisms and drivers involved in the environment at watershed scale under anthropogenic stresses (Larsson et al., 2018). Our previous study showed that spatial differences in ARG richness and abundance overwhelmed their seasonal variations in the Houxi River, and only a few abundant ARGs could contribute to largest of the total ARGs abundance in four seasons . Here we found that ARGs exhibited different profiles along the rural-urban gradient in the Houxi River watershed across five years. Therefore, we went one step further elucidating the mechanism and the main driving factors affecting ARG profile at watershed scale to provide useful information for controlling and reducing the risk of ARGs.

Increased ARG richness and abundance along the urbanization level
Along the rural-urban gradient, the richness and abundance of ARGs significantly rose as the increasing anthropogenic stresses (Fig. 2c, d;  Fig. S3; Fig. S4 and Table S2). About one third of detected ARGs (89 out of 248) were shared between the three sections (i.e. upstream, midstream and downstream) of the Houxi River, with almost all ARGs enriched in the downstream waters ( Fig. 2b and Table S2). The downstream was a rapidly urbanizing region, characterized by much more building type of land use. This result is in line with the previous observations that the concentrations of ARGs would be distinctly high under human-dominated urban conditions (Ouyang et al., 2015;Yang et al., 2017b;Zheng et al., 2018;Zhu et al., 2013). Anthropogenic activities have been demonstrated to have a strong effect on the occurrence and enrichment of ARGs Luo et al., 2010;Ma et al., 2017a;Zhang et al., 2015). High abundance and diversity of ARGs tend to be linked with a higher probability of dissemination and the high abundance can also reflect an antibiotic selective pressure (Bengtsson-Palme and Larsson, 2015). Human activities, including antibiotic therapy, agricultural runoffs and urban wastewater treated or untreated discharging, can lead the release of antibiotics which have promoted the occurrence, evolution and spread of ARGs in the aquatic environment Zheng et al., 2018;Zhu et al., 2017b). Therefore, the increase of urbanization level (represented by urban land percentage) from nature-dominated upstream to human-dominated downstream may result in the significant increase of richness and abundance of ARGs in the downstream waters.
Noticeably, the normalized abundance of ARGs (ARGs/16S rRNA gene) was much higher in the downstream than upstream and midstream stations ( Fig. 2; Fig. S4), indicating that ARGs not only originated from the external environment but also could self-proliferate within the bacterial cell in downstream waters. Also ARGs can be released into the aquatic environment accompanied by their bacterial hosts from surrounding environment (Fang et al., 2019). This is supported by our previous study which showed that the normalized abundance of ARGs was higher in the south/central Chinese lakes than north ones largely due to stronger human activities .
Among the 248 ARGs detected in this study, multidrug resistance genes could be detected in all samples (Fig. S3). This may imply that multidrug-resistance may facilitate broad bacterial distribution in diverse and changing environments. One possible reason would be that in natural environment, because of the wide use of different antibiotics in daily life and animal production, there was a wide diversity of antibiotics Vaz-Moreira et al., 2014). Multidrug resistance genes are capable of conferring resistance to several different types of antibiotic compounds. In addition, multidrug resistance genes can be driven not only by different antibiotics but also by other toxic compounds including biocides, heavy metals, organic solvents and detergents (Alonso et al., 2001;Martinez, 2008).
Ten MGEs, including 2 integrons and 8 transposons, exhibited the highest abundance in downstream waters ( Fig. S3; Fig S4 and Table S2). Some studies reported that MGEs can play a major role in antibiotic resistance development and transferring of ARGs among microorganisms in waste water, even in the clinical environment (Karkman et al., 2018;Martinez et al., 2015). So the emergence and spread of the 10 MGEs detected might be a major challenge for controlling ARGs widespread transmission. Perhaps due to the increase of human activities, the abundance of MGEs in the downstream sites was significantly higher than that in upstream and midstream (Fig. S4), as a previous study showed that one of MGEs, intI1, can be used as a proxy for anthropogenic pollution (Gillings et al., 2015).

Bacterial antibiotic resistance functional and taxonomic communities assembling via different mechanisms
Both the neutral and null models were used to determine whether microbial community assembly could be explained by the stochastic process (Chase et al., 2011;. In our study, taxonomic community fitted well with the neutral model at the watershed scale, revealing that stochastic processes were the main drivers of taxonomic community variability (Fig. 3). This finding is consistent with previous studies in which the neutral model could well explain the patterns of community assembly in the prokaryotic and microeukaryotic plankton Sloan et al., 2006). However, the ARG profile did not fit to the neutral model, indicating that stochastic processes were not an important mechanism driving the ecological process of ARGs. This result is inconsistent with our previous study on ARGs in Xidong Reservoir . As we known, there are several factors (disturbance, habitat connectivity and size, productivity, predation and resource availability) influencing the relative importance of stochastic/ deterministic processes in local community assembly (Zhou et al., 2014). This is probably due to the less human-driven environmental changes of the reservoir in previous study, and the environmental condition was more homogeneous than along the urbanizing river of this study. As a previous study mentioned, the more homogeneous physicochemical conditions appear to favor neutral processes (Logares et al., 2013;Sloan et al., 2006). Another reason might be the scale effect (e.g. the habitat size), that at small space scale like a reservoir, processes related to random chance may be stronger; and at the mesoscale with strong environmental gradients, environmental variables play a dominant role in determining species or gene distributions (Karst et al., 2005). Indeed, along the urbanization level, we found different assembly mechanisms structuring ARG profile (Fig. 3, Fig. S7), supporting the findings on the scale dependency of the assembly processes. Along the whole rural-urban gradient, the effect of anthropologic activities (i.e. environmental effect) exceeded progressively the random effect, which is supported by the result that strong environmental gradients might favor environmental filtering over neutral processes (Logares et al., 2013). As a result, the main assembly mechanism of ARG profile shifted from stochastic processes in upstream and midstream to deterministic processes in downstream waters (Fig. 3b). In downstream regions with high urbanization level, there were more and strong disturbances caused by human. Along this urbanization gradient or level, the ARG profile became more stable as illustrated in the NMDS ordination ( Fig. 2a and Fig. S5), and the similar result was reported from an urban reservoir across one year (Fang et al., 2019). Indeed, the taxonomic F. Peng, et al. Environment International 137 (2020) 105524 community showed a high variability in the downstream largely due to more productive environments which normally result in a higher biodiversity (Chase, 2010). And in downstream waters, high taxonomic variability could often coincide with stable functional structure because the microbial functional structure and taxonomic composition may be shaped by separate processes Vanwonterghem et al., 2014).

Progressive urbanization driven by human activities promoted the stability of ARG profile
Crucially, deterministic process played a relatively important role in explaining the assembly of ARG profile at watershed level, especially in the downstream waters. Further, PLS-PM results showed that the abundance of the 10 MGEs contributed most to ARGs variation ( Fig. 4 and Table S4). Therefore, MGEs might be the main factors directly affecting ARGs distribution in the natural environment through horizontal gene transfer which is a major mechanism for ARG spread (Hu et al., 2016;Ma et al., 2017a;Ochman et al., 2000). As MGEs were identified as one of carriers of ARGs, the genetic linkage and co-selection among them may facilitate ARGs dissemination (He et al., 2016;Zheng et al., 2018).
The origin or enrichment of detected antibiotic resistance genes and the absence or presence of a selective pressure have major implications for managing risks (Bengtsson-Palme and Larsson, 2015). Change of land use in watersheds is a key representative of human activities which has been shown to make a strong contribution to ARG distribution and dynamics Qiao et al., 2018). The composition of ARGs showed a significant difference among different types of land use, indicating that the later can affect the variation of ARGs in a complex way. We found that land use change showed an indirect influence on ARGs by significantly and strongly affecting local environmental conditions, supporting the fact that surrounding landscapes and domestic activity could play a key role in selecting different chemical composition impacting the river (Amos et al., 2015). The inorganic nutrient content (as well as heavy metal), which was as a pollution indicator from municipal wastewater and agriculture, might have positive relationship with the prevalence of some resistance gene levels (Czekalski et al., 2015;Jia et al., 2017). Due to the extension of built-up areas the population density has increased, as have the antibiotics used for medical treatment and animal production . Meanwhile, via urban wastewater and sewerage drainage, a diverse mixture of antibiotics with other pollutants and their metabolites and antibiotic resistant bacteria from surrounding environments can reach the natural rivers, thereby altering the distribution and magnitude of ARGs in the receiving waters, especially in the high-density residential areas (Amos et al., 2015;Marti et al., 2014;Suzuki et al., 2017). In addition, the wastewater discharge as well as the storm waters conduced to the occurrence and spread of ARGs in the receiving water (Jia et al., 2017). Distribution patterns of ARGs are therefore influenced by the multiple dimensions of river connectivity (Ickes et al., 2005;Allan and Castillo, 2007;Wohl, 2017). The present data suggest a more pronounced influence of land use changes, rather than temporal (season) variation ( Fig. 2 and Table 1). Land use changes associated with urbanization drive micro-climate change which is considered to be one of the factors that potentially influence ARGs through the transmission of ARGs and could thereby contribute to the emergence of new resistant strains (Calero-Cáceres et al., 2017;Grimm et al., 2008).
We found that the bacterial taxonomic community showed a weak but significant relationship with ARGs, and played a minor role in explaining the variation of ARGs (Fig. 4 and Fig. S8c). Probably for the sampling effect existing, our sampling sites were in a river ecosystem with low antibiotic selective pressure, and the acquisition of the resistance gene may decrease the fitness of bacteria, so the bacteria tended to loss antibiotic resistance under low antibiotic condition in upstream and midstream waters (Xu et al., 2016). Under this condition, the antibiotic resistant bacteria are rare and the overall carriers of ARGs are weakly relevant under low antibiotic conditions. But as the ARG hosts, we did find that there were co-occurrence relationships between ARGs and some bacterial taxa. For example, we found higher degree of co-occurrence between ARGs and members of Firmicutes and Proteobacteria. Of these bacterial taxa, the genera Acinetobacter, Megamonas, and Comamonas had the highest number of significant correlations with ARG subtypes (Fig. S8b). In aquatic environment, single bacterial taxon significantly correlated with multiple ARG subtypes, suggesting that a bacterial host can harbour multiple ARGs. Recently, Berendonk and his colleagues found that Gammaproteobacteria and Firmicutes are the most frequent carriers of acquired ARGs in diverse environments (Berendonk et al., 2015).
There are some potential limitations that merit further discussion. First, our analyses were at bacterial community level, we did not know the exact host bacteria for each of ARGs at species level. Network results on the co-occurrence patterns between ARGs and bacterial taxa indicated the possible host information of ARGs, but further studies are needed to do and finally verify the ARG bacterial hosts at species level. In this study, however, we focused on different patterns and processes of ARGs and bacteria ecology at community level along rural-urban gradient. The evaluation of exact bacterial host for each ARG seems to go beyond the scope of this study. Second, our findings were based on one watershed investigation, and a key priority for future study should either concentrate on a larger spatial scale investigation, and/or focus on additional field controlled experiments. Even so, our results could contribute to filling the knowledge gap about the distribution pattern and community assembly mechanism of bacterial ARGs, and the relationship between bacterial antibiotic resistance function and taxonomic composition in human-impacted rivers at a watershed scale. We propose that past changes in urban land use are a strong contributing factor or indicator in ongoing ARGs enrichment and dynamic in other regions that underwent similar urbanization changes, therefore watershed management policies aiming at reducing threats posed by ARGs should first consider local urbanization level (i.e. the mode and intensity of human activities).

Conclusions
By characterizing the spatiotemporal distribution of ARGs in a rapidly urbanizing watershed over the period of five years, we found that the richness, abundance and normalized abundance of ARGs exhibited the highest value in the downstream with highest urbanization level. The ARGs profiles showed a stronger deterministic process, although their hosts (bacteria communities) were strongly shaped by stochastic processes at watershed scale. Interestingly, the dominant driver of assembly of ARG profiles shifted from stochastic process in upstream and midstream waters to deterministic process in downstream waters. Further, land use, especially urbanized area percentage, mainly had an indirect effect on ARG profiles distribution. Taken together, our results illustrated a repeated spatial distribution pattern of ARGs in five years across dry and wet seasons, and different community assembly mechanisms of ARGs at different river sections. These findings indicate that anthropologic activities acted as an important ecological driver of ARG dynamics in the Houxi River. Our results also highlight that both bacterial antibiotic resistance functional and taxonomic community compositions should be considered and distinguished in the future studies, due to the decoupling of function and taxonomy at community level.

Declaration of Competing Interest
The authors declare no conflicts of interest. F. Peng, et al. Environment International 137 (2020)