Decoupling the Dynamics of Bacterial Taxonomy and Antibiotic Resistance Function in a Subtropical Urban Reservoir as Revealed by High-Frequency Sampling

Aquatic environments serve as important reservoirs of antibiotic resistance genes (ARGs), but the information on the high-resolution temporal pattern of ARGs in waterbodies is extremely limited. In this study, the weekly dynamics of ARGs and their relationships with microbial taxonomic communities and environmental variables were analyzed in a subtropical urban reservoir over the period of 1 year using high-throughput approaches. In total, 197 ARGs and 10 mobile genetic elements (MGEs) were detected. The results showed that the bacterial community had a seasonal pattern, while ARGs composition did not exhibit seasonality, thereby indicating the asynchrony or decoupling of temporal patterns of microbial taxonomy and function. More importantly, bacterial abundance and community diversity were more strongly correlated with 17 measured environmental variables than ARGs (36 significant correlations for OTUs, 11 for ARGs). However, stochastic processes appeared to have a minor role in the structuring of the ARG profiles, but a more important role in the structuring of bacterial taxonomic communities. Furthermore, we found that precipitation and turbidity were significantly correlated with the richness and diversity of ARGs, suggesting that multiple environmental factors influence the composition and dynamics of ARGs in complex ways. MGEs were abundant and showed significant positive correlations with ARGs, indicating a plausible influence of MGEs on the variation of ARGs. This is the first study which provides an overview of high-resolution dynamics of ARGs in a subtropical waterbody. Our results improve the understanding of microbial processes and mechanisms of ARGs at fine temporal scale, and offer empirical data of use in the monitoring, assessment and management of the urban water environments.


INTRODUCTION
In aquatic ecosystems many studies have been carried out on nitrogen, phosphorus, and microbial community (Yang et al., 2012;Zhang et al., 2019), however various "emerging" pollutants including antibiotics and antibiotic resistance genes (ARGs) have been rapidly drawing increased attention from the public and government (Pruden et al., 2006;Zhu et al., 2017;Yang et al., 2018). Antibiotics, first discovered by Sir Alexander Fleming in 1928, were regarded as panacea for any microbial infections (World Health Organization [WHO], 2014). However, with the emerging problem of antibiotic resistance, antibiotic resistant bacteria (ARB) and ARGs have been identified as a significant issue, posing a serious threat to the health of humans and the environment (Pruden et al., 2006;Berendonk et al., 2015). Antibiotic resistance is ancient (D'Costa et al., 2011), but with the continuous increase in discharge of overused antibiotics, the development, spread, and enrichment of ARGs has accelerated globally (Laxminarayan, 2014), and heavy metal has been reported to play an important role in the enrichment of antibiotic resistance by co-selection mechanisms (Stepanauskas et al., 2006). More importantly, ARGs can be exchanged and transferred among environmental bacteria, pathogens and nonpathogens via horizontal gene transfer (HGT) carried out by mobile genetic elements (MGEs) (Kruse and Sørum, 1994;Martínez et al., 2015). Aquatic environments are considered to be significant reservoirs of both ARB and ARGs, providing ideal settings for the occurrence and dissemination of ARGs (Marti et al., 2014). Due to long-term misuse and overuse of antibiotics by human, urban water often contains diverse ARGs at high levels (Xu et al., 2016), which may pose high risk to human health during the interaction with receivers. Therefore, ARGs occurrence and their temporal dynamics in the aquatic environment have become an important issue in environmental science .
Microbial communities can vary between environments, and play an important role in drive global biogeochemical cycling, but it is not well known how variation in taxonomic composition relates to function (Louca et al., 2016). Previous studies reported that bacterial community composition closely correlates with ARG profiles and controls the transfer of ARGs in antibiotic-rich environments (e.g., soil, urban wastewater, and sewage sludge) (Forsberg et al., 2014;Su et al., 2015), but the contribution of microbial community shifts to antibiotic resistance function variations in natural waterbodies with low and medium antibiotic concentrations remains largely unknown, especially at a fine temporal scale. Furthermore, microbial communities have been reported to exhibit seasonal variability in the environment (Treusch et al., 2009), and high-frequency sampling analysis could reveal and disentangle detailed temporal dynamics of bacterial community (Lindh et al., 2015). Previous studies have assessed the temporal variation of ARGs and supported the occurrence of seasonal variations of ARGs based on only two or four sampling time points in 1 year Zheng et al., 2018). However, research on the dynamics of ARGs in urban water environment, and their relationship to microbial taxonomic community and environmental factors based on highfrequency time series is still lacking.
Revealing the mechanisms that drive the variation of microbial communities and ARG profiles are major challenges. It has been shown that microbial communities are simultaneously shaped by stochastic (neutral) and deterministic processes (Sloan et al., 2006;Logares et al., 2013). The neutral model is a stochastic process that develops community through births, death and immigration, can successfully predict microbial community assembly dynamics (Hubbell, 2001). However, the neutral model cannot completely descript the community assembly without deterministic process such as environmental factors and competition (Sloan et al., 2006;Ofiţeru et al., 2010). Both stochastic (neutral) and deterministic processes need to be considered when investigating the assembly of microbial communities and ARG profiles, but there is still a lack of studies on the contribution of both processes in structuring ARG profiles Guo et al., 2018).
The urban reservoir in this study is located in a rapidly urbanizing area of Xiamen, it receives numerous pollutants from surrounding housing estates by point and non-point source pollutions. However, this urban reservoir is the most important landscape waterbody in the Jimei district, so the ARB and ARGs in the water may pose ecological risk in the interaction with inhabitants and tourists. In this study, we simultaneously investigated the dynamics of bacterial OTUs and ARGs in the urban reservoir using high-throughput qPCR and 16S rRNA gene sequencing on a weekly basis over a 1-year period. We aimed to (1) characterize the temporal dynamics of both bacterial taxonomic community and antibiotic resistance gene profiles; (2) reveal the relationships between bacterial taxonomic composition and antibiotic resistance function at different levels of taxonomic resolution; and (3) explore the community assembly processes and mechanisms shaping the ARGs composition dynamics. Additionally, we hypothesize that the temporal patterns of microbial OTUs and ARGs are asynchronous, and their temporal patterns are influenced by different processes and mechanisms. To our best knowledge, this is the first study to analyze a high-resolution time series of bacterial OTUs and ARGs communities in urban water environments. Therefore, this study may provide a fundamental data for evaluating the ecological risks of ARGs in urban water environment and offer useful information for government and policy maker.

Study Area and Sample Collection
All samples were collected once a week from September 2016 to August 2017 in station G (24 • 36 N, 118 • 04 E) of Xinglinwan Reservoir, which is situated in the Houxi River watershed, Xiamen city, Fujian province, southeast China (Supplementary Figure S1). In total, 51 samples were collected in this study. Several stormwater and wastewater inlets were found near the station G along the reservoir shoreline within a distance of 1000 m. The Xinglinwan Reservoir is a subtropical eutrophic urban reservoir, and is characterized by high turbidity and low water transparency. It has multiple functions in flood control, aquaculture and tourism, exerting a great impact on the surrounding residents.
Surface waters were collected in a sterile container from the top 0.5 m water, immediately transported to laboratory and kept in the dark until pre-treatment within 0.5 h. In order to remove the influence of large plankton and particles, samples were prefiltered through a 200 µm mesh, and then filtered through 0.22 µm polycarbonate membranes (47 mm diameter, Millipore, Billerica, MA, United States) using a vacuum filtration system. To ensure the sufficient microorganisms, the filtrated water volume ranged from 200 to 500 mL for all samples, because we kept a similar filtrating time (about 50 min) for each of membranes. The membranes were placed in sterilized tubes and stored at -80 • C until DNA extraction.

Environmental Parameters
A total of 17 environmental variables were determined in this study (Supplementary Figure S2). The comprehensive trophic state index (TSIc) were calculated according to Yang et al. (2012). Water temperature (WT), pH, dissolved oxygen (DO), turbidity, electrical conductivity (EC), salinity, and oxidation reduction potential (ORP) were measured in situ with a Hydrolab DS5 multiparameter water quality analyzer (Hach Company, Loveland, CO, United States). Chlorophyll a concentrations were quantified by PHYTO-PAM Phytoplankton Analyzer (Heinz Walz GmbH, Eichenring, Germany). Total carbon (TC), total organic carbon (TOC), total nitrogen (TN), ammonium nitrogen (NH 4 -N), nitrate nitrogen (NO 3 -N), nitrite nitrogen (NO 2 -N), total phosphorus (TP), and phosphate phosphorus (PO 4 -P) were measured according to the standard methods described in our previous study (Liu et al., 2013). The precipitation data referred to the cumulative rainfall during the last week before the sampling day, and the data were collected from the Xiamen Meteorological Bureau.

DNA Extraction and High-Throughput Sequencing
Total DNA was directly extracted from each membrane using the FastDNA SPIN Kit (MP Biomedical, Santa Ana, CA, United States) according to the manufacturer's protocol, and subsequently the concentration and quality of the DNA were determined using a NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, United States).
To investigate the bacterial taxonomic community, the V3-V4 hypervariable regions of the 16S rRNA gene were amplified, purified, and quantified according to our previous procedure . The PCR products were pooled and sequenced on the Illumina Miseq platform (Illumina, Inc., San Diego, CA, United States) using a paired-end (2 × 250 bp) sequencing strategy. All the raw sequences were quality-controlled using MOTHUR v1.39.0 (Schloss et al., 2009). The unoise3 pipeline was used to pick operational taxonomy units (OTUs) at 3% dissimilarity level (Edgar, 2010). OTU sequences were taxonomically classified by USEARCH (sintax) against the Greengenes database (DeSantis et al., 2006). All eukaryotic, chloroplast, archaeal, mitochondrial and unknown sequences were removed from the data set. Finally, all sequences were normalized to 50,014 sequences for each of 51 samples, and these sequences were clustered into 9834 OTUs.
The raw data for 16S rRNA gene analysis in this study have been deposited in the NCBI sequence read archive (SRA) database under the BioProject number PRJNA510463 and the accession number SRP173857.

High-Throughput Quantitative PCR (HT-qPCR)
High-throughput quantitative PCR was performed to detect the abundance of ARGs using the Wafergen SmartChip realtime qPCR platform (Wafergen Biosystems, Fremont, CA, United States). A total of 296 primer sets targeting 285 resistance genes, 8 transposase genes, the class I integron-integrase gene (intI), the clinical class 1 integron-integrase gene (cIntI), and 16S rRNA gene were used. These primer sets and PCR reaction protocols were listed in our previous studies Liu et al., 2018). Three technical replicates were performed for each sample, and non-template negative control was also included for each primer set. The absolute abundance of ARGs and MGEs was calculated by normalizing their abundance to the absolute 16S rRNA gene copy number, which was estimated by real-time quantitative PCR  as described below.

Real-Time Quantitative PCR (qPCR)
Real-time quantitative PCR was used to quantify the absolute copy number of 16S rRNA genes in all samples. It was performed on a Lightcycler 480 instrument (Roche, Basel, Switzerland). The 20 µL reaction contained 10 µL 2× LightCycler 480 SYBR Green I Master Mix (Roche, Basel, Switzerland), 7 µL nucleasefree PCR-grade water, 2 µL diluted DNA, and 0.5 µM of each primer. The reactions were performed in triplicate with negative controls. The following thermal cycling conditions were used: initial incubation at 95 • C for 5 min, followed by 40 cycles of 95 • C for 15 s, 60 • C for 1 min, and 72 • C for 15 s.
A six-point calibration curve was generated from 10-fold dilutions for an external standard calculation. The melting curves were used to analyze the specificity of PCR products, and PCR efficiency ranged from 96 to 105% in this study.

Statistics
The Bray-Curtis similarity was calculated using PRIMER v7.0 (Clarke and Gorley, 2015), and non-linear regression of time lag with community similarity was analyzed to reveal the temporal patterns of bacterial OTUs and ARGs communities. The nonlinear regression was the fit to the sinusoidal curve to explore the seasonality in SigmaPlot v12.0 (Systat Software Inc., Chicago, IL, United States). The diversity indexes were calculated using the vegan package in R v3.3.1. The SIMPER analysis was performed to identify the contribution of bacterial OTUs and ARGs to the temporal variation of communities or profiles (Clarke and Gorley, 2015).
Spearman's and Pearson's correlations were calculated using SPSS v20.0 (IBM Corp, Armonk, NY, United States), and the Mantel test was performed in an R environment with the vegan package. Further, a Procrustes test was conducted with the vegan package in R environment to explore the synchronicity of bacterial OTUs communities and ARG profiles (Forsberg et al., 2014). Network analysis was used to reveal Spearman correlations between ARGs and bacterial taxa with the picante package in R v3.3.1, and the correlation coefficients |ρ| > 0.6 between microbial taxa and ARGs, and |ρ| > 0.4 between MGEs and ARGs, were considered at P < 0.01. Network information was generated by the Gephi v8.2 and Cytoscape v3.6.1.
In addition, we used redundancy analysis (RDA) to explore the relationship between environmental factors and the compositions of bacterial OTUs and ARGs subtypes. We further used partial redundancy analysis (pRDA) to evaluate the relative contribution of bacterial taxonomic community, environment variables and MGEs to the ARG profiles of all samples by R v3.3.1 with the vegan package. To reveal the assembly processes underlying both bacterial OTUs and ARGs compositions, we used the neutral community model, which can evaluate the importance of stochastic processes in explaining the assembly of bacterial taxonomic and functional communities (Sloan et al., 2006). This model predicts the relationship between the frequency at which taxa occur in a set of local communities (in this case, bacterial OTUs and ARGs at one sampling point) and their abundance across the wider metacommunity (bacterial OTUs and ARGs sampled over a year). The neutral models were performed in R v3.3.1 using the Hmisc, stats4, and minpack.lm packages.

Dynamics of Environmental Variables, Bacterial OTUs, and ARGs Compositions
The major physicochemical variables exhibited a great change over time (Supplementary Figure S3). Comprehensive TSIc varied from middle eutrophic level to hypereutrophic level, and nutrient concentrations were high throughout the study period, with higher value in 2017 than 2016. Water temperature ranged from 15 to 35 • C during the 1-year sampling. Higher precipitation values were recorded in September 2016 and June 2017.
Time series analysis indicated relative stability at bacterial phylum and ARG type level, but very rapid or abrupt fluctuations among both bacterial OTUs and ARG subtypes (Figure 1). The relative abundance of ARG subtypes exhibited distinct variations across the year at weekly intervals ( Figure 1D). For example, the relative abundance of qacE 1-02 (multidrug resistance gene) varied from 0 to 0.5 and its relative abundance changed almost every week. Time-lag regression analysis showed that the bacterial community displayed a good fit to a sinusoidal curve ( Figure 1E), but ARG profiles failed to fit such a curve ( Figure 1F); indicating that bacterial community composition exhibited seasonality while ARGs composition did not. The general variation of ARG similarity was narrower compared to the bacterial taxonomic community (Supplementary Figure S6), suggesting that the compositions of resistomes were more stable. Additionally, SIMPER analysis indicated that Proteobacteria contributed most to the variation within the microbial community (31.86%), while beta-lactamase resistance gene contributed most to the variation of ARGs (23.41%) ( Supplementary  Table S1). Furthermore, Procrustes test revealed asynchrony or decoupling between bacterial ARGs and OTUs communities (r = 0.347, P < 0.01, Figure 2A). However, the Pearson's correlation indicated a significant but weak relationship between bacterial ARGs and OTUs communities ( Figure 2B).

Co-occurrence Patterns Between Bacterial Taxonomic Community and Antibiotic Resistome
Network analysis was performed to reveal the significant and robust correlations between ARG subtypes and bacterial OTUs. The whole network consisted of 1,547 nodes and 2,139 edges, and showed a modularity index of 0.762. Aminoglycoside resistance genes had the most edges (489) with bacterial taxa, followed by tetracycline (474), MLSB (454), and beta-lactamase (318) resistance genes (Figures 3A,B). Meanwhile, Proteobacteria, Bacteroidetes, Firmicutes, and TM7 were significantly correlated with most ARGs. The co-occurrence network showed that the hubs of all modules were ARGs (Supplementary Figure S7A), indicating the important roles of ARGs in the network.
We also used Mantel tests to evaluate the relationships between bacterial phyla and ARG types. Similar to the results of the network analysis, most ARGs were significantly and strongly correlated with Bacteroidetes, Firmicutes, and Proteobacteria (Supplementary Table S2). However, sulfonamide resistance genes did not show any significant correlation with the bacterial community, and MGEs only significantly correlated with Chloroflexi, Firmicutes, and TM7.

Relationships Between Environmental Factors, Bacterial OTUs, and ARGs
Mantel tests showed aminoglycoside, beta-lactamase, MLSB, tetracycline, and vancomycin resistance genes were significantly correlated with precipitation (Supplementary Table S3). RDA further indicated that water temperature, chlorophyll a, salinity, total carbon, TOC, total nitrogen, nitrate nitrogen, and total phosphorus were significantly correlated with the β-diversity of the bacterial community, while salinity, total carbon, nitrate nitrogen, total phosphorus, and precipitation were strongly correlated with ARGs compositions, suggesting that nutrients were key factors that influence the dynamics of both ARGs and bacterial communities (Figure 4). We found that the bacterial community (R 2 = 0.365) showed a much better fitted to the neutral model (Figure 4) than ARGs (R 2 = 0.079), indicating that stochastic process played an important role in microbial community assembly but not in the resistome profiles.
The absolute abundance of ARGs was positively correlated with total carbon and TOC, while ARGs richness was significantly and positively correlated with turbidity and precipitation ( Figure 5A). Additionally, turbidity and precipitation also exhibited a positive correlation with richness and α-diversity of bacterial community, but were not significantly correlated with 16S rRNA gene abundance ( Figure 5A). Partial redundancy analysis showed that the bacterial community, MGEs and environmental factors separately explained 19.23, 19.87, and 9.4% of the variation of ARG profiles, respectively. Interaction between FIGURE 3 | Significant correlation among bacterial OTUs, ARGs, and MGEs. Classification of nodes within the network depicting the co-occurrence patterns of microbial taxa and ARGs (A), MGEs and ARGs (C). Other, other bacterial phyla; others, other ARGs. The Spearman's absolute rho was calculated at OTU and ARG subtype levels, and the number of edges was summarized at phylum and type levels, respectively. The edges correspond to a strong (Spearman's |ρ| > 0.6 for microbial taxa and ARGs, Spearman's |ρ| > 0.4 for MGEs and ARGs) and significant (P < 0.01) correlation between nodes. Frequency distributions of correlation coefficients between ARGs and microbial taxa (B) and between ARGs and MGEs (D) at P < 0.01. all variables explained 4.63% of the variation, and 42.49% of the variation was unexplained.

Correlations Between ARGs and MGEs
The absolute abundance of ARGs was positively correlated with MGEs (Supplementary Figure S8), and the normalized abundance of MGEs was significantly correlated with the richness and normalized abundance of ARGs (Supplementary Table S4). The network analysis further indicated that aminoglycoside, beta-lactamase, MLSB, and tetracycline resistance genes were significantly correlated with MGEs (Figures 3C,D), and the co-occurrence pattern revealed that Tp614 and tnpA-03 play the most important roles in the network (Supplementary Figure S7B). However, the integrase gene intI-1 only showed a significant correlation with multidrug and vancomycin resistance genes (Supplementary Table S5), but the absolute abundance of transpose gene tnpA-05 showed a strong correlation to all ARG types except vancomycin resistance genes ( Supplementary Table S6). Additionally, the transposase genes IS613, tnpA-03, tnpA-04, and tnpA-07 also showed a significant correlation with some types of ARGs. Partial redundancy analysis showed that 19.87% of the variation of ARG profiles was explained by MGEs (Figure 5B), indicating that MGEs were an important contributor to the change of ARG profiles during our study period. These results suggested that MGEs played a crucial role in the change, accumulation and dissemination of ARGs in this urban landscape reservoir.

The Contrasting Temporal Patterns Between Bacterial OTUs and ARGs
Our high-frequency time series-based study revealed clear differences in the temporal patterns of microbial communities and ARG profiles, along with the decoupling of the temporal dynamics of bacterial taxonomy and antibiotic resistance function in an urban waterbody. The composition of ARG subtypes varied frequently over the year but was relatively stable at the ARG type level. None of the ARG subtypes was present in all 51 samples, and no distinct seasonal patterns in ARG profiles could be determined, but microbial taxonomic community had significant seasonal pattern.
ARGs in this urban reservoir were detected at relatively high abundance and richness. The diversity of ARGs in this study was greater than that reported from Chinese suburban/rural lakes/reservoirs . The high abundance and richness in this waterbody may be driven by the combined human activities of sewage discharge and agricultural/urban runoff, as it is generally known that these activities increase the number of ARGs in waterbodies Ahmed et al., 2018;Liu et al., 2018). Moreover, aminoglycoside and multidrug resistance genes dominated the abundance of ARGs in this study, while aminoglycoside and beta-lactamase resistance genes dominated the types or richness of ARGs. In contrast, Jia et al. (2017) found in swine wastewater and downstream water that tetracycline and aminoglycoside resistance genes were the dominant ARGs, perhaps due to the different antibiotic selection pressures in different environments. However, the primer limitation of this study may result in the bias of detected ARGs, as the PCRbased method only covered limited number of ARG types. Additionally, a majority of ARGs investigations have been carried out in municipal wastewater treatment plants. An et al. (2018) showed significant changes of ARG profiles based on samples just collected in February and August, while another study based on monthly samples for 1 year did not find obvious seasonal variation in methicillin-resistant gene although it varied between FIGURE 5 | Variation of ARG profiles explained by different ecological variables. (A) Spearman's correlation between bacterial OTUs, ARGs, MGEs, and environmental parameters. ORP, oxidation reduction potential. * P < 0.05; * * P < 0.01. (B) Partial redundancy analysis illustrating the effects of bacterial community (relative abundance of ten major phyla), MGEs and environment factors on ARGs distribution pattern. Environmental factors referred to significant correlated factors, including salinity, total carbon, nitrate nitrogen, total phosphorus, and precipitation. months and sampling sites (Börjesson et al., 2009). Yang et al. (2013) suggested that some types of ARGs in activated sludge showed seasonal fluctuations during a 4-year investigation, while most types of ARGs did not show seasonal fluctuations, but these results were based on samples collected only in summer and winter. However, according to our study, ARG subtypes were extremely dynamic and some of them frequently showed rises and falls over short-time periods (i.e., in just one week) but did not show seasonal pattern, so the previous low-frequency based studies may not completely and effectively reflect actual seasonal dynamic patterns of ARGs at fine temporal scale. Overall, our study based on high-frequency sampling over long timescales provide better understanding of the detailed temporal pattern of ARG profiles.
Our data revealed that some bacterial OTUs (e.g., Planktothrix) fluctuated rapidly, while major phyla remained stable at a broad taxonomic level. The microbial community showed an obviously seasonal variation and a trend of resilience during 1 year (12 months). In the microbial community, Proteobacteria and Actinobacteria were the dominant phyla, consistent with previous studies in drinking water system (Jia et al., 2015) and rivers (Jia et al., 2017). Additionally, the high relative percentage of Cyanobacteria was previously reported to have significant impacts on the variation of ARGs . Wang et al. (2018) showed that Proteobacteria, Actinobacteria and Cyanobacteria in the environment can carry some ARGs, confirming the high normalized abundance of ARGs and the vital roles of these bacteria in the dissemination of ARGs. Furthermore, we found a significant seasonality in the microbial taxonomic community, and the Procrustes test indicated the asynchronous temporal patterns between microbial OTUs and ARGs compositions (Figure 2). Previous studies have shown that ARGs strongly correlate with microbial communities in antibiotic-rich environments including human gut  and landfill leachate . Recently, Ma et al. (2017) indicated that ARGs and bacterial communities did not show a strong correlation in drinking water, which is similar to the result of our study under low or medium antibiotic pressure. Altogether, the weak correlation between ARGs and bacterial communities in this study may be due to the relatively low antibiotic and heavy metal selection pressure in the natural environment.

The Co-occurrence Pattern of Bacterial OTUs and ARGs
Previous studies have confirmed that network analysis is a useful and reliable method to decipher the co-occurrence patterns between bacterial taxa and ARG subtypes, and to track potential hosts of ARGs (Barberán et al., 2012;Li et al., 2015;Ma et al., 2017;Feng et al., 2018). The module hubs in the network function as keystone of the module (Faust and Raes, 2012). All module hubs in our co-occurrence network were affiliated to ARGs, implying that ARGs were important in maintaining the structure of the network and that HGT of ARGs was frequent among bacteria in the microbial community. The high frequency of aminoglycoside resistance genes in the network may be due to their high abundance and richness in the community. However, the abundance and richness of tetracycline resistance genes were relatively low in our studied system, but they connected tightly with various bacterial taxa. This suggests that the tetracycline resistance genes were shared among multiple bacterial taxa ). An increasing number of studies have reported the widespread distribution of tetracycline resistance genes (Czekalski et al., 2014;Zhou et al., 2017), as this antibiotic has been widely used in China (Hvistendahl, 2012).
In the co-occurrence network, Proteobacteria, Bacteroidetes and Firmicutes were the most prevalent phyla. These taxa were predicted as possible ARG hosts, and/or they originated from a similar source including upstream wastewater and surface runoff, suggesting the potential dissemination risk of ARGs with these taxa. In fact, Proteobacteria include a wide variety of pathogens and occur frequently as opportunistic pathogens, and have been found in soils that contained almost all major ARG types (Forsberg et al., 2014). Moreover, Jang et al. (2018) found that Bacteroidetes can be potential hosts of the intI-1 gene. In this study, Bacteroidetes were significantly correlated to all major ARG types except sulfonamide resistance genes (Figure 3 and Supplementary Table S2), suggesting the potential transmission of ARGs to human pathogenic bacteria. The Firmicutes were strongly correlated to aminoglycoside resistance genes in sewage sludge , and the ARGs were shown to frequently switch hosts from Firmicutes to other bacterial phyla . In our study, Firmicutes were not the dominant bacteria but may carry diverse ARGs and further cause ARGs pollution to different environments or ecosystems. Nevertheless, the co-occurrence patterns revealed by this study were based on correlation, so the exact ARG carrying host and transmission pathway of ARGs need to be verified by further study.

Assembly Mechanisms Driving the Dynamics of ARGs
Our results indicated that the variation in ARGs could be explained by the bacterial community, MGEs and environmental factors. Environmental factors including salinity, total carbon, nitrate nitrogen, total phosphorus and precipitation influenced the temporal variation of ARGs, while precipitation was the strongest factor. A recent study reported that nutrients explained a large part of ARGs variation and that the nutrients combined with other factors could drive the distribution of ARGs . The seasonal stormwater runoff induced by rainfall can increase the absolute abundance of ARGs and MGEs in an aquatic environment  and could influence the ARGs load through increased input from runoff and soil resuspension (Di Cesare et al., 2017). In accordance with previous studies, we found that precipitation displayed a significant influence on the ARG composition and structure, and was positively correlated with the richness and normalized abundance of ARGs, supporting the view that rainfall contributed to the loading of ARGs in the river and reservoir waters in urban region. In addition, rainfall significantly influenced pathogenic bacterial abundance because stormwater was able to significantly contribute to the occurrence and elevated concentration of pathogenic bacteria in subtropical waters (Ahmed et al., 2018). Our results further showed that turbidity was strongly correlated to richness and diversity of bacteria and ARGs. A turbidity increase is normally associated with rainfall-induced contamination, suspended particulate matter and sediment disturbance. This may further drive changes in bacterial community and ARG profiles. It has been demonstrated that turbidity was significantly correlated to bacterial diversity, reiterating the influence of runoff on the bacterial community (Peter and Sommaruga, 2016). Taken together, our results indicated that precipitation had both direct and indirect effects on bacterial taxonomic communities and resistomes, supporting the idea that stormwater with agriculture/urban runoff, with contaminants from the urban region could increase loadings of bacterial ARGs thereby increasing the risk of ARGs dissemination (Garner et al., 2017). In contrast, we did not find any significant relationship between water temperature and the ARGs composition or abundance, corroborating the fact that the ARGs in this subtropical urban reservoir did not have a significant seasonal pattern. This is inconsistent with a previous finding that water temperature was a potential factor driving the dynamics of ARG profiles (Zheng et al., 2018), and this inconsistency may be attributed to a difference in sampling frequency and climate conditions, as the range of temperature was narrower in this subtropical reservoir compared to other studies in temperate zone. In the Zheng et al. (2018) study, water samples were only collected in March, June, September, and December, and their river showed a wider temperature range than our studied reservoir. In our study, the 42.49% of ARGs variance remained unexplained, may partly be due to the environmental factors that were not considered, such as anthropogenic factors and the co-selection of antibiotics and heavy metals (Stepanauskas et al., 2006;Liu et al., 2018;Chen et al., 2019).
The bacterial community has been identified as one of the key drivers that shape the ARG profiles in antibioticrich environments (Forsberg et al., 2014;Su et al., 2015). However, our results of pRDA, Procrustes test and Spearman's relationship all showed that microbial community was weakly correlated with ARG profiles. More interestingly, our results indicated that the neutral process played a more important role in the assembly of microbial taxonomic communities than in ARG profiles, as the neutral model only explained a minor part (7.9%) of the variation in ARG profiles. Therefore, dispersal and stochastic processes appeared to have more influence on bacterial taxonomic communities compared to ARG profiles. This may be one of the reasons why microbial OTUs and ARGs communities exhibited a decoupling change over 1 year. These results indicated that the vertical transmission (mother-to-child transmission) of ARGs in this environment may not play a leading role.
Our study showed that MGEs could partly explain the variation of ARG profiles in the urban reservoir, and the absolute abundance of MGEs showed a significant and positive correlation with ARGs. Further, MGEs that were always abundant in our samples, could potentially enhance the probability of HGT of ARGs. Previous studies reported that environmental contaminants at low concentration have a potential role in the dissemination of ARGs by promoting HGT, such as antibiotic and non-antibiotic pharmaceuticals (Andersson and Hughes, 2014;Wang et al., 2019). In this study, the rainfall combined with possible contaminants and MGEs may largely facilitate the proliferation and spread of ARGs among different bacterial taxa. These results suggest that MGEs played an important role in the temporal variation of ARGs in the urban waterbody, implying the potential risk that the ARGs may be transferred by MGEs from aquatic environments to human pathogens (Martínez et al., 2015).

CONCLUSION
Our weekly study revealed the asynchronous succession patterns of bacterial taxonomic communities and ARG profiles, indicating a decoupling dynamic between bacterial taxonomy and function on a temporal scale in a subtropical urban reservoir. The bacterial taxonomic community was significantly correlated with water temperature and exhibited seasonality, while ARG profiles was strongly correlated with precipitation and did not show distinct seasonality. Precipitation had both direct and indirect effects on the temporal patterns of bacterial community and ARG profiles, reflecting the effect of stormwater runoff on temporal dynamics of ARG profiles. The diversity and composition of bacterial taxonomic communities appeared to be more responsive to environmental variables than ARGs, and the composition of ARG profiles appeared to be governed less by stochastic processes. Additionally, MGEs were positively correlated with ARGs and had significant contribution to the variation of ARGs community, suggesting the potential risk of HGT. This study may help to better understand the temporal dynamic of ARG profiles in urban water environments at fine temporal scale, and highlight that further high-frequency and long-term time series research should be undertaken to reveal the specific drivers and deep mechanisms of ARGs assembly and dynamics.

DATA AVAILABILITY
The raw data for 16S rRNA gene analysis in this study have been deposited in the NCBI SRA database under the BioProject number PRJNA510463 and the accession number SRP173857.

AUTHOR CONTRIBUTIONS
JY designed the experiments. PF, FP, and XG carried out the sample collection and determined the environmental parameters. PF, FP, and PX conducted the PCR, highthroughput qPCR, high-throughput sequencing, and bioinformatics. PF and JY wrote the first draft of the manuscript. All authors contributed to and have approved the final manuscript.