The raw water from the selected drinking water treatment plant (DWTP) (Guangzhou, China) originates in the North Rivers of Pearl River, and its phased treatment includes flocculation, sedimentation, filtration, and chlorination (Figure 1). Water samples were collected along the treatment processes of the DWTP and distribution pipeline, including samples for raw water (A), flocculation tank effluent (B), sedimentation tank effluent (C), quartz sand (D), sand filter effluent (E), chlorine-disinfected water (F), and tap water (G). They were respectively collected from a water pumping station (5 L), flocculation tank (18 L), sedimentation tank (18 L), sand filter (36 L), clear water tank (54 L), and tap water (about 2,000 L) (Figure 1). Additionally, quartz sand samples (D) were obtained from the upper part of the quartz-sand layer in the sand filter (50 g). To minimize the temporal variation, three replicate water samples were simultaneously collected from each location (21 samples in total) in July 2019. All of the samples were immediately transported to a laboratory for the experiments. The water sampling and pre-treatment methods were performed as previously described²⁷. DNA for each sample was extracted using a FastDNA SPIN Kit for Soil (MP Biomedicals, CA, United States), and the DNA concentration and purity were measured using a microspectrophotometry (Nano DropND-2000, NanoDrop Technologies, Wilmington, DE, United States).

For each sampling site, temperature, pH, and dissolved oxygen (DO) level were measured in situ using an HQ40D portable multi-parameter water quality analyzer (HACH, Loveland, CO, United States). Turbidity was measured using a 2100Q portable turbidimeter (HACH, Loveland, CO, United States). Ammonia-nitrogen, chemical oxygen demand (COD), and total organic carbon (TOC) were measured using a DR1900 portable spectrophotometer (HACH, Loveland, CO, United States). Adenosine triphosphate (ATP) level was measured with an ATP assay kit (Beyotime, Shanghai, China). All determinations were performed in triplicate.

Bacteria were isolated from the water samples of the different sampling sites. Briefly, 0.2 ml of water samples or proper dilution of the water samples was spread on R2A agar plates with an L-type coating rod. Then, all the plates were incubated at 22°C for 7 days. Subsequently, total heterotrophic bacterial counts (HPCs) were determined by counting the number of colonies on the R2A agar plates. Next, colonies with different morphologies were picked and purified using the streaking method (Massa et al., 1998). The process was repeated until a pure bacterial strain was obtained. To identify the obtained bacteria, all pure isolates were first incubated in R2A broth for 24 h, and then, genomic DNA for each pure isolate was extracted using a DNA extraction kit (Magen Biotech, Guangzhou, China). Subsequently, the 16S rRNA gene was amplified in a 25-μl volume using universal primers 27 F and1492 R (Mas-Ud et al., 2020), and amplification was performed as previously described (Gu et al., 2018). Then, PCR products were sequenced by BGI Tech Solutions Co., Ltd (Guangzhou, China), and 16S rRNA gene sequences were submitted to the NCBI public library of bioinformatics for comparison and identification.¹

Antibiotic resistance of the representative bacterial isolates from the DWTS was determined using the standard Kirby–Bauer disc diffusion method according to Clinical and Laboratory Standards Institute Guidelines (CLSI, 2020). The antibiotics and concentrations of the disks (Oxoid, Basingstoke, United Kingdom) used in this study were as follows: piperacillin-tazobactam (TZP, 100/10 μg), norfloxacin (NOR, 10 μg), clindamycin (DA, 2 μg), cefepime (FEP, 3 μg), ceftriaxone (CRO, 30 μg), ciprofloxacin (CIP, 5 μg), penicillin (P, 10 μg), cefotaxime (CTX, 30 μg), chloramphenicol (C, 30 μg), ampicillin (AMP, 10 μg), trimethoprim-sulfamethoxazole (SXT, 1.25/23.75 μg), erythromycin (E, 15 μg), aztreonam (ATM, 30 μg), meropenem (MEM, 10 μg), cefoxitin (FOX, 30 μg), amikacin (AK, 30 μg), kanamycin (K, 30 μg), tobramycin (TOB, 10 μg), amoxicillin-clavulanate (AMC, 30 μg), tetracycline (TE, 30 μg), gentamicin (CN, 10 μg), ceftazidime (CAZ, 30 μg), and levofloxacia (LEV, 5 μg). Additionally, multiple antibiotic resistance index (MARI) values were calculated for all the isolates to evaluate whether they originated from high-risk antibiotic pollution sources (Labella et al., 2021).

The extracted DNA from each sample was sequenced using the Illumina HiSeq 2500 platform (Illumina Inc., San Diego, CA, United States). Over 10 Gb of raw data was generated for each sample (Supplementary Table 1). Quality control and assembly were conducted as described previously (Bai et al., 2021). To obtain low-abundance species, the obtained clean data were mapped to its own Scaftigs to obtain un-using reads, and then all the un-using reads were collated for hybrid packaging. The open reading frames (ORFs) within the contigs of each sample were subsequently predicted using Prodigal.² The gene prediction results were used for gene clustering by Linclust.³ The longest sequence in each cluster was taken as the representative sequence to obtain the non-redundant gene catalogue (Unigenes). Clean reads were mapped to the gene catalogue by BBmap,⁴ and the relative abundance of each gene in each sample was calculated using the number of reads mapped and gene length with the following formula:

where k stands for any identified gene, G_k is the relative abundance of each identified gene, r_k is the number of mapped metagenomic reads for each identified gene, and L_k is the length of each identified gene. ARGs were detected against a structured ARG reference database (SARG version 2.0) with the following parameters: under a cutoff of ≥80% similarity and ≥70% query coverage (Yin et al., 2018). Furthermore, an online analysis pipeline (ARGs-OAP) based on the SARG database was used for the ARG-like sequence annotation from the metagenomic sequencing reads with default settings (Yin et al., 2018). Subsequently, the abundance of ARGs was normalized and expressed as the number of copies of ARGs per cell (copy/cell) (Ma et al., 2017). A contig was retrieved as an ARG-carrying contig (ACC) if it contained ≥ one ARG-ORF (Zhang et al., 2019). The contigs that harbored ARGs were selected, and the protein sequences on the contigs were used to make a comparison against the local NCBI NR database using BLASTP with a cutoff of E-value ≤ 10^–5 (Mao et al., 2014). Subsequently, the MEtaGenome ANalyzer (MEGAN, version 6) was used to align the profiles. If ≥ 50% of ORFs on a contig were attributed to the same taxon, these ARG-carrying contigs (ACCs) were assigned to this taxon (Ishii et al., 2013). ARG-carrying pathogens were determined according to a previous study (Li et al., 2015).

The plasmids, integrons, and insert sequences were identified based on aclame, integration, and ISfinder databases with the following parameters: under a cutoff of identity ≥80% and query cover ≥70%. Metagenome binning, bin checking, and taxonomic annotation were performed as previously described (Poghosyan et al., 2020). The metagenome bins that were >50% complete with <10% contamination were retained for building the phylogenomic tree using GToTree (https://github.com/AstrobioMike/GToTree/wiki). The abundance [reported as reads per kilobase per million (RPKM)] of genes was conducted as previously described (Potgieter et al., 2020). The occurrence of ARGs in chromosomes or plasmids was also determined as previously described (Camacho et al., 2009), and MAFFT (Katoh and Standley, 2013) was used to align the amino acid sequences of IntI genes recovered from the metagenomes together with the reference sequences of class 1 (IntI1, AAQ16665.1) integrases. Raw data have been submitted to the NCBI Sequence Read Archive (SRA) under accession number SRR17063149-17063169.

The results of the physicochemical and microbiological parameters of water samples are illustrated in Table 1. The water temperature was 27–28°C, and the water pH was 7.63–8.16. The turbidity (NTU) decreased significantly from 69.47 to 2.62 NTU after the settling process. Water DO maintained a steady state during the water treatment process. However, the ATP and HPC gradually decreased. After chlorination, no cultivable bacteria were detected. However, in the tap water, 10 CFU/ml HPC were detected. COD_Mn decreased gradually with the treatment process and increased slightly after entering the pipe network. TOC was always kept at a very low level. Other parameters, such as TN, also demonstrated a downward trend. They remained at a low concentration in the chlorine-disinfected water, even in the tap water. All the physicochemical qualities of the chlorine-disinfected water met the drinking water standards.

The dominant bacterial phyla and genera were identified based on the metagenomic data in every sampling site (Figure 2). A total of 95 bacterial phyla were identified. The relative abundances of the top nine bacteria at the phyla level in all samples are illustrated in Figure 2A. Proteobacteria and Actinobacteria were the dominant phyla in the DWTP, and they shifted largely during the water treatment process. After flocculation, the relative abundance of Proteobacteria was greatly improved and then remained at a certain ratio until chlorination, whereas the proportion of Actinobacteria decreased sharply after flocculation, followed by a slight increase after the settling process. However, the fate of the Actinobacteria was the same as the Proteobacteria; it decreased greatly after chlorination but maintained a certain ratio.

Previous studies have indicated that δ-Proteobacteria could be used as a biomarker to distinguish source water and tap water (Han et al., 2020). In other words, δ-Proteobacteria exhibit relatively low chlorine resistance, but α-Proteobacteria, β-Proteobacteria, γ-Proteobacteria, and Actinobacteria demonstrate strong resistance to the water treatment processes, including chlorination. Similarly, another author found that α-Proteobacteria and Actinobacteria were dominant in all the chlorine-disinfection samples (Li et al., 2020). Thus, in this study, we observed that the phyla of Proteobacteria and Actinobacteria were still predominant in the finished water. Moreover, the microbial community structure in the sand biofilm was significantly different from the water samples. The relative abundance of Nitrospirae, Planctomycetes, Chloroflexi, and Acidobacteria in the sand biofilm was higher than in the water samples. We also observed that the microbial community in the sand biofilm shaped the sand filter effluent. This finding is consistent with our previous study (Hou et al., 2018) and others (Oh et al., 2018; Zhou et al., 2021). Eustigmatophyceae, Firmicutes, and Planctomycetes in the finished water increased when they passed through the pipe network into the tap water, which might be involved in the regrowth and secondary contamination of bacteria in the pipe network (Douterelo et al., 2016). Some previous studies have demonstrated that the bacterial community is strongly correlated with the resistance profile (Jia et al., 2015; Zhang L. et al., 2018; Chen et al., 2019). In particular, Proteobacteria and Acidobacteria are important hosts of ARGs (Ma et al., 2019), and some studies have detected significant positive correlations between them and ARGs (Su et al., 2015; Xu et al., 2020). Thus, the changes in the relative abundance of Proteobacteria and Acidobacteria might explain the variations in the abundances of ARGs during the water treatment processes.

As evident in Figure 2B, at the genus level, the dominant genera were Polynucleobacter, Acidovorax, Limnohabitans, and Nitrospira, of which Acidovorax has been reported to be the host of ARGs in the DWTP (Jia et al., 2019). After chlorine disinfection, the relative abundance of Stenotrophomonas greatly increased and became the predominant bacteria. Stenotrophomonas has been reported to be typical chlorine-resistant bacteria (Luo et al., 2021), and it harbors many ARGs (Gil-Gil et al., 2020). Some studies reported that chlorine disinfection leads to an increase in the abundance of antibiotic-resistant bacteria (Khan et al., 2016), and chlorine disinfection affects the horizontal transfer of ARG through various pathways, including increasing the relative abundance of MGEs (Zhang et al., 2019). Moreover, some researchers have demonstrated that bacteria that survived chlorine disinfection might have higher levels of antibiotic resistance than other bacteria (Jia et al., 2015). Furthermore, the increased ARG abundance of chlorine-resistant bacteria might increase antibiotic resistance in pathogenic bacteria. Therefore, chlorine-resistant bacteria deserve attention.

Metagenomic data based on the structured SARG database were used to profile the antibiotic resistome in the DWTS. All the identified ARG-like sequences were classified into 13 ARG types (e.g., aminoglycoside and bacitracin resistance gene, as illustrated in Supplementary Material 1) according to the set script. As illustrated in Figure 3A, the abundance of 13 ARG types detected in all samples was normalized as the number of copies of ARGs per cell (copy/cell). Gene-encoded bacitracin (0.3 × 10^–2 to 0.40 copies/cell, mean value = 0.23 ± 0.11 copies/cell), multidrug (0.02 to 0.23 copies/cell, mean value = 0.13 ± 0.05 copies/cell), and sulfonamide (0.2 × 10^–3 to 0.7 × 10^–1 copies/cell, mean value = 0.05 ± 0.02 copies/cell) ARGs were detected at the highest levels of abundance in the A–E samples. Previous studies have also reported that bacitracin, multidrug, and sulfonamide ARGs were predominant in drinking water samples (Ma et al., 2019; Zhang et al., 2020d).

However, aminoglycoside, beta-lactam, and rifamycin ARGs were notably enriched after chlorination. Moreover, in tap water, beta-lactam and multidrug ARGs retained a high abundance when they passed through the drinking water distribution system, which might threaten human health. From another perspective, the bacitracin ARG was mainly distributed in the B and C samples, which account for about 36 and 30%, respectively (Supplementary Figure 1 and Supplementary Material 2). The multidrug ARGs were mainly distributed in the B and F samples, which account for 20 and 26%, respectively (Supplementary Figure 1), while the sulfonamide ARGs were distributed in the A, B, C, and E samples. Similarly, as evident from Figure 3B, bacitracin and sulfonamide ARGs can be significantly removed by sand filters and chlorination. In other words, the DWTS effectively removed these ARGs, while the other ARGs cannot be eliminated very well. Some ARG types increased through the DWTS, such as multidrug, beta-lactam, aminoglycoside, and rifamycin (Figure 3B).

In total, 464 ARG subtypes were annotated (e.g., bacA, sul1 and multidrug_transporter, illustrated in Supplementary Material 1). Figure 3C illustrates the profile of the top 20 ARG subtypes during the drinking water treatment processes. Bacitracin_bacA (0.1 to 0.5 copies/cell) had the highest levels of abundance in the A–E samples, followed by sulfonamide_sul1 (0.2 × 10^–1 to 0.9 × 10^–1 copies/cell) and multidrug_mdtB (0.2 × 10^–1 to 0.5 × 10^–1 copies/cell). After chlorination, gene-encoded multidrug (acrB, amrB, and ceoB), bacitracin (bcrA), sulfonamide (sul1 and sul2), chloramphenicol (floR), tetracycline (tetA and tetG), polymyxin (arnA), macrolide-lincosamide-streptogramin (ereA), and beta-lactam (OXA-21) decreased sharply. Previous studies have reported that tetA could be removed using free chlorine, and this is similar to our results (Amarasiri et al., 2020). However, some researchers have demonstrated that the abundance of sul1 and sul2 increased after chlorination (Liu et al., 2018; Zhu et al., 2021). Subsequently, the shift of these ARGs remains controversial when chlorine is used.

Additionally, multidrug (multidrug_transporter, mexE, smeB, smeC, smeD, smeE, and smeF), aminoglycoside [aadA, aph (6)-1, aph (3’)-IIb], beta-lactam (metallo-beta-lactamase, TEM-1, TEM-117, TEM-157, TEM-187, TEM-205, TEM-75, TEM-91, and TEM-118), and rifamycin (rifampin monooxygenase) were enriched. With respect to the phenomenon that some ARGs increased after chlorine disinfection, some researchers have pointed out that free chlorine could promote the transformation of ARGs through a series of cell responses, for example, raising levels of reactive oxygen species (ROS), ROS-mediated DNA damage, and bacterial membrane damage (Zhang et al., 2021). Our results indicated that chlorination reduced the diversity of the ARG types. Additionally, the Venn diagrams for the different drinking water treatment processes display the number of unique and shared ARG subtypes at each sampling site (Figure 3D). The shared 33 ARG subtypes, which constituted 10% of the total ARG subtypes, were universally present from the raw water to tap water (Supplementary Material 3). The shared ARG subtypes were distributed in aminoglycoside (n = 1), aminoglycoside (n = 2), bacitracin (n = 2), beta-lactam (n = 2), fosmidomycin (n = 2), macrolide-lincosamide-streptogramin (n = 2), multidrug (n = 18), rifamycin (n = 1), tetracycline (n = 1), vancomycin (n = 1), and unclassified (n = 1), respectively. Thus, the largest variety of shared ARG subtypes belonged to the multidrug ARGs. These shared ARG subtypes persistent in the DWTS will lead to continuous pollution. It is also worth noting that the unique ARG subtypes in process F (n = 14) that escaped from chlorine disinfection were assigned to beta-lactam (n = 11), followed by macrolide-lincosamide-streptogramin (n = 2) and chloramphenicol (n = 1), respectively. Furthermore, the unique ARG subtypes in process G (n = 24) were assigned to aminoglycoside (n = 1), beta-lactam (n = 18), bleomycin (n = 1), macrolide-lincosamide-streptogramin (n = 2), multidrug (n = 1), and vancomycin (n = 1), respectively. Therefore, it is evident that the unique ARG subtypes in tap water rebounded when passed through the distribution system, which requires further attention. Shared ARG subtypes or unique ARG subtypes in the processes F and G pose a risk to human health when people drink the water. Additionally, these ARG subtypes might be transmitted to human living areas; for example, when people use water to wash food, such as fruit and vegetables, the ARGs will contaminate the fruit and vegetable surfaces and affect health through the food chain. Recently, scholars have discovered that ARGs in vegetables or fruit pose risks to ecological environment health (Sun et al., 2021). Furthermore, when people use the water to take a bath, the ARG subtypes might also remain on the skin (Nielsen and Jiang, 2019).

As evident from Figure 4A, most of the ACCs were annotated as fragments of Proteobacteria (46.6%). Additionally, these contig-harbored ARGs were mainly resistant to bacitracin (26.08%), aminoglycoside (4.74%), multidrug (4.31%), and beta-lactam (3.45%). Similarly, the Proteobacteria phylum was identified as the most abundant bacterial host of ARGs (25.3–94.6%) in DWTPs (Zhao et al., 2022). Proteobacteria phylum was also regarded as a major pool of ARGs in sewage sludge (Xu et al., 2018), industrial mariculture systems (Wang et al., 2018), and the human gut (Hu et al., 2016). This is supported by the findings of our research. As illustrated in Figure 4B, before chlorine disinfection (processes A–E), the main putative ARG hosts were Acidovorax and Polynucleobacter, which were also the dominant bacteria in the raw water (Figure 2B and Supplementary Table 2). In the flocculation tank effluent, we found that Acidovorax and Stenotrophomonas were the main putative ARG hosts, while, in the sedimentation effluent, the main putative ARG carriers were Acidovorax, Polynucleobacter, and Stenotrophomonas. The main ARG carrier in the quartz sand sample was Pseudomonas, whereas Nitrospira, the predominant bacteria in the quartz sand biofilm, was not the dominant ARG host. Therefore, sand filters might be low risk for the transmission of ARGs. In the sand filter effluent, the main ARG carrier was Limnohabitans. Therefore, in general, we can see that Acidovorax and Stenotrophomonas were the most abundant ARG carriers of bacteria before chlorine disinfection.

However, after chlorine disinfection (processes F and G), the main ARG carriers were Stenotrophomonas and Escherichia. Interestingly, after the chlorine disinfection, the relative abundance of Stenotrophomonas also greatly increased and became the predominant bacteria (Figure 2B). Stenotrophomonas has been reported to be a highly chlorine-resistant Gram-negative bacterium that can regrow post disinfection (Shekhawat et al., 2020, 2021), and it is regarded as a global emerging Gram-negative multiple drug-resistant organism (Brooke, 2012). Similarly, Escherichia was widely reported to be an indicator of drug-resistant genes (Zhang et al., 2020e) and easily escaped the oxidation by disinfectants (Mounaouer and Abdennaceur, 2016). We assumed that these bacteria might regrow in the drinking water distribution system, or there was secondary pollution in the pipe network. Therefore, these bacteria are of high concern for public health safety in drinking water environments.

From all the water samples, as illustrated in Figure 4C, the dominant putative ARG hosts include Acidovorax (6.0%), Polynucleobacter (4.3%), Pseudomonas (3.4%), Escherichia (1.7%), Klebsiella (1.5%), Polaromonas (1.1%), Mycobacterium (0.9%), and Limnohabitans (0.6%). Acidovorax was the most abundant bacterium that carried bacitracin ARGs as Polynucleobacter, Polaromonas, and Limnohabitans. A previous study also reported that Acidovorax was the main carrier of bacitracin_bacA (Asante and Sekyere, 2019), and Polynucleobacter was an abundant ARG host in drinking water sources and frequently carried bacitracin ARGs (Bai et al., 2019). It was reported that the putative ARG host distribution pattern might be most likely responsible for the overall ARG profile (Bai et al., 2019). Therefore, these bacteria were most likely to drive the shift of bacitracin ARGs. Mycobacterium drove the shift of rifamycin ARGs, and Pseudomonas, Escherichia, and Klebsiella carried the largest variety of ARGs. According to the further analysis of ARG-carrying contigs (ACCs), we found that 49 ACCs were annotated as pathogens (Supplementary Material 4). They were mainly classified as Stenotrophomonas maltophilia (n = 15), P. aeruginosa (n = 12), E. coli (n = 7), Klebsiella pneumoniae (n = 6), Acinetobacter baumannii (n = 3), Rhodococcus erythropolis (n = 2), Pseudomonas putida (n = 2), Mycobacterium smegmatis (n = 1), and Enterobacter aerogenes (n = 1). Thus, we can see that the most common pathogens were S. maltophilia and P. aeruginosa in the DWTS. Moreover, Escherichia spp. mainly carried tetracycline and sulfonamide ARGs (Figure 4C), while Stenotrophomonas spp. mainly carried aminoglycoside, beta-lactam, and multidrug ARGs (Figure 4C). Additionally, pathogens were distributed in processes A (n = 41), B (n = 40), C (n = 40), D (n = 12), E (n = 41), F (n = 26), and G (n = 20). In particular, S. maltophilia was persistent in every process (Figure 4D and Supplementary Material 5), and the sustained survival of these bacteria in the DWTS could lead to the horizontal propagation of ARGs and even increase risks to human health.

As illustrated in Figure 5A, ARGs present in the DWTS were more prevalent in bacterial chromosomes than in plasmids. Specifically, ARGs for bacitracin, rifamycin, vancomycin, quinolone, polymyxin, kasugamycin, and fosmidomycin were just prevalent in chromosomes. These ARGs might not contribute to horizontal gene transfer. While ARGs for multidrug, beta-lactam, aminoglycoside, tetracycline, sulfonamide, trimethoprim, and chloramphenicol were prevalent mainly in plasmids. Significantly, chloramphenicol was only presented in plasmids. Therefore, we must focus on these kinds of ARGs with a higher risk of transmission. Additionally, the phylum of Proteobacteria harbored the highest abundance of plasmids (Figure 5B), followed by Firmicutes. Meanwhile, the abundance of plasmids was higher than that of insert sequence (IS) and integron (Figure 5C). The most abundant MGE plasmids appeared in the sand filter effluent, followed by the flocculation tank effluent.

However, in the quartz sand biofilm, the abundance of MGE plasmids was the lowest, and the distribution trend of MGE ISfinder in drinking water treatment processes was similar to the MGE plasmids. Interestingly, ARG abundance in the quartz sand biofilm was also the lowest (Supplementary Figure 2). At the same time, the Shannon index indicated that microbial diversity was the most abundant in the quartz sand biofilm. These results demonstrate that the probability of horizontal transmission of ARGs was the least in the quartz sand filters and the highest in the flocculation tank effluent. Integron class1 is a major player in disseminating antibiotic resistance genes across pathogens and commensals. Therefore, in this study, we also looked for integron class1 and found a higher prevalence in the DWTS; the abundance of integron class1 in the flocculation tank effluent water was the highest than the other samples (Supplementary Figure 3). Additionally, it decreased after chlorine disinfection and significantly decreased in the process of drinking water transmission, as well as MGE plasmids and MGE ISfinder. Plasmids and integrons are prone to recombination (Mazel, 2006); thus, chlorine disinfection and residual chlorine in the distribution system could prevent the dissemination of ARGs.

A network analysis based on the metagenomic approach was used to discern the co-occurrence patterns between ARG subtypes and demonstrated a good advantage (Ju et al., 2016; Zhang et al., 2020a). As illustrated in Figure 6 and Supplementary Table 3, a total of 90 pairs of ARG subtypes under a cutoff of p < 0.01 and ρ > 0.8 were discovered from 29 ARG subtypes assigned to 12 ARG types. ARG subtypes assigned to the same ARG type (i.e., intra-types) tend to be co-occurring. Overall, ten out of the 90 pairs were intra-type co-occurrence correlations. For instance, emrE and qacEdelta1 for multidrug resistance correlated (r = 0.85), as well as the sul1 and sul2 (r = 0.93) for sulfonamide resistance and tetA and tetG (r = 0.88; Supplementary Table 3) for tetracycline resistance. Co-occurrence patterns were also observed between different ARGs subtypes (i.e., inter-types). For example, aac (3)-I for aminoglycoside resistance strongly correlated with ksgA (r = 1) for kasugamycin resistance, bacA for bacitracin resistance correlated with sul1 for sulfonamide resistance (r = 0.96), floR for chloramphenicol resistance correlated with tetG for tetracycline resistance (0.95), and qacEdelta1 for multidrug resistance correlated with sul1 for sulfonamide resistance (0.93). Recent research also addressed the significant positive correlations among these ARG types and subtypes (Zhang S. et al., 2018; Zhang et al., 2020a).

A network analysis was performed to explore the co-occurrence patterns among ARG subtypes, MGEs, and bacterial genera. Co-occurrence patterns between ARG subtypes and bacterial genera could be a good way to track the potential hosts of ARGs in their environment. We hypothesized that the ARGs and bacterial taxa had a strong and significant positive correlation under a cutoff of p < 0.01 and ρ > 0.8. As illustrated in Supplementary Figure 4 and Supplementary Material 6, the co-occurrence results indicated that 42 genera might be potential hosts for 27 ARG subtypes. Among the 42 bacterial genera, Acidovorax, Curvibacter, Variovorax, Comamonas, and Hydrogenophaga carried more diverse ARG subtypes (14, 14, 13, 12, and 12, respectively) than other genera. For example, Acidovorax was identified as the possible host of the multidrug-resistance gene (emrE and qacEdeltal), sulfonamide-resistance gene (e.g., sul1 and sul2), tetracycline-resistance gene (tetA and tetG), trimethoprim-resistance gene (dfrA1), beta-lactam-resistance gene (OXA-21 and class-A-beta-lactamase), chloramphenicol-resistance gene (floR), bacitracin-resistance gene (bacA), macrolide-lincosamide-streptogramin-resistance gene (ereA), and aminoglycoside-resistance gene [aac(6_)-I] (Supplementary Material 6). Likewise, Curvibacter was the possible host of the bacitracin-resistance gene (bacA), beta-lactam-resistance gene (OXA-21 and class-A-beta-lactamase), multidrug-resistance gene (emrE and qacEdeltal), sulfonamide-resistance genes (e.g., sul1 and sul2), tetracycline-resistance gene (tetA and tetG), and trimethoprim-resistance gene (dfrA1).

Additionally, a total of 57 pairs with p < 0.01 and ρ > 0.8 were observed between the MGEs and bacterial genera. For instance, intI1 had a significant and positive correlation with ten bacterial genera (Supplementary Material 6). IS21 correlated positively and significantly with 15 bacterial genera, and IS3 had a strong positive collection with 10 bacterial genera. Furthermore, a total of 31 pairs with p < 0.01 and ρ > 0.8 were observed between the MGE and ARG subtypes. IntI1 had a significantly positive correlation with eight ARG subtypes, and IS21 had a significantly positive correlation with nine ARG subtypes. These results implied the important role of IntI1 and IS21 in antibiotic resistance dissemination. It has been reported that MGE-associated ARGs tend to persist in the treatment at DWTPs and might pose a higher risk of resistance dispersal (Zhao et al., 2022). Similarly, recent literature has also reported the increase in the proportion of MGE-associated ARGs in wastewater treatment plants (Che et al., 2019).

In this study, after metagenomic binning, 33 high-quality metagenome-assembled genomes (MAGs) were obtained, which mainly included Gammaproteobacteria (n = 13), Alphaproteobacteria (n = 5), Betaproteobacteria (n = 2), Planctomycetota (n = 2), Bacteroidota (n = 1), Actinobacteriota (n = 2), and Nitrospirota (n = 1) (Figure 7). The general characteristics of these 33 MAGs are illustrated in Supplementary Table 4. Of the 33 MAGs, the phylum Proteobacteria (n = 20) dominated the microbial community in the DWTS, followed by Planctomycetota (n = 2) and Actinobacteriota (n = 2). The MAG identified as S. maltophilia-like (bin.195) harbored the greatest number of ARG subtypes (n = 8), including multidrug (n = 6; smeD, semE, multidrug_transporter, mexE, semB, and smeC), beta-lactam (n = 1; metallo-beta-lactamase), and aminoglycoside [n = 1; aph(3’)-IIb] (Supplementary Table 5). Meanwhile, Burkholderiaceae-like (bin.379) also contained bacitracin (n = 1; bacA). Additionally, the MAG (bin.84) identified as R. erythropolis-like contained aminoglycoside [n = 2; ant (3”)-Ih-aac (6’)-Iid and ant (2”)-I] and rifamycin (n = 1; rifampin monooxygenase). Interestingly, based on our genome data of strain F2 identified to the genus of Stenotrophomonas sp., which also harbored smeD, semE, multidrug_transporter, semB, and smeC (Supplementary Material 7 and Supplementary Figure 5), these results can be verified against each other.

A total of 93 isolates were recovered from the DWTS using traditional culture technology (Supplementary Tables 6, 7). The isolates belong to 30 genera based on a 16S rRNA gene sequencing analysis, and they were found to belong to five phyla: Proteobacteria, Firmicutes, Actinobacteria, Deinococcus-Thermus, and Bacterioidetes. Proteobacteria (42%) and Firmicutes (38%) were the dominant phyla. As illustrated in Supplementary Table 7, at the genus level, Bacillus was the most abundant genus (20%), followed by Pseudomonas (11%) and Acidovorax (9%), respectively. Moreover, 47 isolates within 20 genera were recovered from the sand filter, representing the maximum microbial diversity in all the sampling sites. Bacteria isolated from this process accounted for 51% of all the isolates, meaning that bacteria in the sand filter were the most abundant. These results were consistent with the microbial diversity analysis based on the metagenomic data in this study.

It is worth noting that some pathogens and opportunistic pathogens such as P. aeruginosa, Bacillus cereus, Staphylococcus hominis, E. coli, and Aeromonas enteropelogenes were detected in water treatment processes. Fortunately, only one bacterium was detected after chlorination, and it was not a pathogen or opportunistic pathogen, which indicated that chlorination effectively removed high-risk bacteria in the DWST. However, it does not demonstrate that the tap water was completely microbial safe. In our study, five isolates were obtained from the tap water, indicating secondary pollution in the drinking water pipe network. Previous literature has reported that disinfectant could induce culturable bacteria and various pathogenic bacteria into a viable but non-culturable (VBNC) state (Lin et al., 2017) in the DWTS, which would result in these bacteria not being detected using traditional culture methods. Therefore, a comprehensive microbiological safety assessment of the drinking water could combine molecular methods with culture-based methods.

The antimicrobial susceptibility testing results (Table 2) demonstrated that selected isolates had high resistance to FOX (56%), 48% of the selected isolates were resistant to P, and 28% were resistant to AMC. However, the resistance to TOB, DA, and E was lower (4, 8, and 8%, respectively). In the Proteobacteria, 88% of the isolates were resistant to AMC and 63% were resistant to AMP. However, the incidence of resistance to TE, TOB, CIP, CAZ, and MEM was very low. In the Firmicutes, 71% of the isolates were resistant to FOX, 59% to P, and only 12% to TE, DA, and E. Bacteria were regarded as multidrug resistant (MDR) if they were resistant to three or more groups of antibiotics. In this study, 28% of the selected isolates belong to MDR, and the occurrence of MDR phenotypes was mainly related to the bacterial species. The MARI was applied to verify the origin of the antibiotic resistance in the isolates. The MARI value of the selected isolates recovered from the DWTS was 0.09–0.50. Compared to other environments, the MARI value of the isolates that recovered from the DWTS was lower. For example, the MARI value of Salmonella enterica that recovered from broiler chickens and retail shops was 0.2–0.6 (Elkenany et al., 2019), and E. coli that isolated from clinical specimens, including eye and blood cells, had higher MARI values of 0.82 and 0.74, respectively (Masoud et al., 2021). Additionally, according to the MARI value, the B3, B11, D38, D20, and F2 isolates might be from a source where antibiotics were used frequently.