Investigation of the Prevalence of Antibiotic Resistance Genes According to the Wastewater Treatment Scale Using Metagenomic Analysis

Although extensive efforts have been made to investigate the dynamics of the occurrence and abundance of antibiotic resistance genes (ARGs) in wastewater treatment plants (WWTPs), understanding the acquisition of antibiotic resistance based on the WWTP scale and the potential effects on WWTPs is of relatively less interest. In this study, metagenomic analysis was carried out to investigate whether the WWTP scale could be affected by the prevalence and persistence of ARGs and mobile genetic elements (MGEs). As a result, 152 ARG subtypes were identified in small-scale WWTP samples, while 234 ARG subtypes were identified in large-scale WWTP samples. Among the detectable ARGs, multidrug, MLS (macrolide–lincosamide–streptogramin), sulfonamide, and tetracycline resistance genes had the highest abundance, and large and small WWTPs had similar composition characteristics of ARGs. In MGE analysis, plasmids and integrons were 1.5–2.0-fold more abundant in large-scale WWTPs than in small-scale WWTPs. The profile of bacteria at the phylum level showed that Proteobacteria and Actinobacteria were the most dominant bacteria, representing approximately 70% across large- and small-scale WWTPs. Overall, the results of this study elucidate the different abundances and dissemination of ARGs between large- and small-scale WWTPs, which facilitates the development of next-generation engineered wastewater treatment systems.


Introduction
The WHO has announced that the rapid development of the spread of antibiotic resistance is one of the top 10 global public health threats facing humanity [1]. Due to the overuse and misuse of antibiotics in human activities such as agriculture, livestock, and medical treatment, antibiotic resistance has rapidly accelerated in environmental ecosystems [2,3]. Aquatic environments are regarded as a major source for the aggregation and dissemination of antibiotic resistance genes (ARGs) because they are easily influenced by human activities [4,5]. In particular, the emergence of ARGs could significantly increase the spread of antibiotic resistance, especially via mobile genetic elements (MGEs) including plasmids, integrons and insertion sequences [4][5][6]. There is growing evidence that wastewater treatment plants (WWTPs) have received special attention in playing a major role in the proliferation of antibiotic resistance into other environmental systems, such as rivers, sediment, and oceans [6]. Because WWTPs are mainly designed for the removal of physical and chemical pollutants, ARGs and antibiotic-resistant bacteria (ARB) are usually not considered for the effective removal of antibiotics in WWTPs. Therefore, WWTPs can be a reservoir and vehicle for the transmission of biological contaminants such as pathogenic bacteria and antibiotic resistance mechanisms to other environmental media [7].
Recent studies have consistently reported relatively high levels of ARGs and MGEs and even heavy metal resistance genes in WWTPs and in each process in wastewater treatment systems [6,8,9]. However, most previous studies preliminarily focused on the identification of ARGs and ARB during the wastewater treatment process, such as activated sludge and anaerobic digestion processes, not considering the scale and treatment capacity of WWTPs [6,8]. According to a previous study, the capacity of a WWTP can affect the bacterial community via various operation practices, such as influent characteristics, sludge retention time (SRT), hydraulic retention time (HRT), reactor volumes, and physical constraints [10]. Recent studies have suggested that the bacterial community plays a key role in ARG prevalence and enrichment during the WWTP process [11]. In addition, coselection of the bacterial community and ARGs is favored during the wastewater treatment process due to the massive amount of nutrient substances and certain operation conditions [12].
Although extensive efforts have been made to investigate the dynamics of the occurrences and abundances of ARGs in the activated sludge process using omics approaches, such as metagenomics, metatranscriptomics, metaproteomics, and real-time qPCR approaches [11,[13][14][15], the understanding of antibiotic resistance and MGE contamination according to WWTP capacity is understudied. In addition, it is still unclear to what extent water-borne ARGs promote the acquisition of antibiotic resistance among bacterial communities via MGEs according to the WWTP scale. Therefore, this study aimed to investigate whether the WWTP scale could be affected by the prevalence and persistence of ARGs and MGEs. For this, shotgun metagenomics was applied because it is a good approach to gain deeper insights into the entire genomic pool from the mixed environmental microbiome [16][17][18].

Study Area and Sampling
Twelve samples (influent, activated sludge, effluent) in this study were collected and two duplicate biological replicate samples were considered for each sampling point from 4 full-scale WWTPs in Busan city, South Korea, in April 2020. The designed A and B WWTP scales were 9.5 × 10 4 m 3 /d and 4.0 × 10 4 m 3 /d, respectively, while the C and D WWTP scales were designed to be 34.0 × 10 4 m 3 /d and 45.2 × 10 4 m 3 /d, respectively. Detailed information on the investigated WWTPs is provided in Table 1. Briefly, all WWTPs were activated sludge processes equipped with a series of treatment units, including coagulation, sedimentation, sand filtration, and disinfection. The major difference between the small and large WWTPs is that small WWTPs use the conventional activated sludge process, which includes the following treatment units: primary settling tanks, modified Ludzack-Ettinger (MLE) tanks and secondary settling tanks. However, large WWTPs consist of an activated sludge process with the same treatment method as small WWTPs and a membrane bioreactor (MBR) system using microfiltration (MF) membrane units. The MLE process was designed for the treatment of wastewater mainly from municipal sources with a total design capacity of 27.5-34.5 × 10 4 m 3 /d and MF membranes of 6.5-10.0 × 10 4 m 3 /d. All samples were stored at −80 • C in a refrigerator until DNA extraction.

DNA Extraction and High-Throughput Shot-Gun Sequencing
Total genomic DNA was extracted using the FastDNA Spin Kit for Soil (MP Biomedicals, CA, USA) according to the manufacturer's suggested protocol. Extracted DNA concentrations and purities were measured using a 1% gel and a Nano300 Micro Spectrophotometer (Allsheng, HangZhou, China). Approximately, 5 µg of each purified DNA sample was used for shotgun paired-end library construction. Briefly, DNA fragments were end-polished, A-tailed, and ligated with the full-length adaptor for Illumina sequencing with further PCR amplification and purification. The insert sizes of libraries were analyzed using an Agilent 2100 bioanalyzer (Agilent Technologies, Palo Alto, CA, USA). Shotgun metagenome sequencing was performed by Macrogen (Seoul, Korea) using the TruSeq DNA PCR-Free Kit (Illumina Inc., San Diego, CA, USA) according to the library protocol, where index coding had been added. Each sample was barcoded and analyzed using a 2 × 101 bp paired-end protocol with an Illumina HiSeq 2000 plat- form. The total output of the 12 metagenome samples was approximately 120 Gb. All metagenomes generated in this study are publicly available via MG-RAST under sample IDs mgs821935, mgs821941, mgs821947, mgs821953, mgs821959, mgs821965, mgs821971, mgs821977, mgs821938, mgs821989, mgs821995, and mgs822001.

Bioinformatics Analysis
Raw sequences from the Illumina HiSeq were trimmed with respect to adapters and low-quality nucleotide stretches using Trimmomatic (version 0.39) [19] with default settings. Trimmed and filtered sequences were then used as input for all further analyses. Clean sequences were analyzed using the ARG-OAP v2.0 pipeline [20] to explore the diversity and abundance of different types and subtypes of ARGs from nine metagenomic datasets in this study. ARG types and subtypes were annotated using default criteria [21]. Metagenomics rapid annotation subsystem technology (MG-RAST) was used to identify bacterial taxonomic classification [22], and a column chart comparing the relative abundance of each class was generated. All clean data were also searched for MGEs against an offline database, including integrons, gene cassettes, insertion sequences and plasmids, downloaded from the INTEGRALL and NCBI RefSeq databases [6,23]. If the nucleotide sequence identity of the best BLASTn hit was ≥90% over an alignment of ≥50 bp with an e-value ≤ 1 × 10 −5 , it was annotated as an integron, insertion sequence or gene cassette. For plasmid annotation, the sequence identity of the BLASTn hit was ≥ 95% with an alignment length ≥ 90 bp [24,25]. Relative abundances of MGEs were calculated using plasmid-and integron-like reads per total metagenomic sequencing read. All statistical analyses were conducted with R software. The correlation coefficient was calculated using Spearman's rank correlation coefficient method, applying a statistical test with r > 0.6 and p < 0.05.

Diversity and Occurrence of ARGs from Metagenomic Analysis
Metagenomics analysis revealed that 152 ARG subtypes belonging to 18 ARG types were identified in the small-scale WWTP samples, while 234 ARG subtypes belonging to 21 ARG types were identified in the large-scale WWTP samples. This could be also due to the different number of inhabitants and equivalent of the scale ( Table 1). The total abundance of ARGs ranged from 0.38 to 133.07 ppm in small-scale WWTPs and from 2.94 to 145.57 ppm in large-scale WWTPs. As shown in Figure 1, similar types of ARGs, such as aminoglycosides, bacitracin, beta-lactams, multidrug resistance, MLS (macrolide-lincosamide-streptogramin), sulfonamide, and tetracycline resistance genes, were detected between small and large WWTPs, despite differences in their treatment capacities. These abundant ARGs are usually associated with antibiotics used extensively in human or veterinary medicine, including as growth promoters [1,26]. Genes for resistance to MLS, multidrugs, and tetracycline were represented by 70% of all wastewater treatment processes (influent, activated sludge, effluent). Mostly, the abundance and diversity of the detectable ARGs in influent (1.59-145.57 ppm) were significantly higher than those in activated sludge (4.81-38.79 ppm) and effluent (0.38-26.39 ppm), indicating that the biological wastewater treatment process could decrease the ARGs effectively. However, the abundances of multidrug resistance genes were still high, with a range of 21.47 to 26.39 ppm after the treatment process. Interestingly, some resistance genes, such as tetracycline and sulfonamide resistance genes, usually decreased in the effluent during the treatment process in small-scale WWTPs (A and B) but increased in large-scale WWTPs (C and D). Generally, the removal efficiencies for all WWTP samples were approximately 90% after the treatment process, which may have led to a significant reduction in ARG detection in the effluent samples. The removal of antibiotics is attributed to biodegradation and biosorption onto activated sludge [27][28][29]. In addition, the MLE anoxic/aerobic tanks in the activated sludge process may not provide good conditions for the proliferation of ARG-carrying bacteria; hence, the replication and dissemination of ARGs were not favored [6,30]. However, it is of concern that tetracycline and sulfonamide resistance genes were still high in the effluent because these resistance genes can persist under harsh Antibiotics 2021, 10, 188 5 of 13 environmental conditions, which may potentially transfer antibiotic resistance to bacteria via HGT (horizontal gene transfer) [21,31].
Considering the relative abundance from heatmap analysis, tet genes, OXA genes, mex genes, mdt genes, and sul genes were the dominant ARGs in the activated sludge samples ( Figure 2). However, aad genes, ere genes, erm genes, mdt genes, sul genes, and tet genes were relatively dominant subtypes in the effluent samples. Tetracycline antibiotics are adsorbed by activated sludge flocs or biofilms and are largely concentrated in sludge at levels as high as 0.1-1 mg/L [32][33][34]. The tetM gene is commonly localized on a conjugative transposon, Tn916, with a broad host range [35], possibly explaining the observed diversity of the tetM-carrying taxa in effluents. TetQ is normally associated with conjugative elements [36], whereas tetW has been found to be associated with a conjugative transposon potentially transferred into other genera [37]. MLS also tends to be hydrolyzed or sorbed. Because of the strong sorption of tetracycline and MLS antibiotics, their mobility in the environment may be facilitated by transport with wastewater [38]. Notably, adsorption is one of the major reasons for the development of ARBs along with HGT [12]. Erm resistance genes can easily be transferred from one host to another [39] since they are usually acquired and associated with mobile elements, such as plasmids [40] and transposons. Aminoglycoside (aadA, aadA1 and aadA2) and beta-lactam (blaVEM and blaOXA) resistance genes are frequently detected ARGs in WWTPs and are closely related to antibiotics important for human infection treatment or veterinary usage [21]. These characteristics possibly explain why tetracycline, MLS, sulfonamide and multidrug resistance genes were mostly high in abundance in the 12 metagenomes.
As such, the NMDS (nonmetric multidimensional scaling) results showed that the ARG subtypes were primarily influenced by the wastewater treatment process, not the WWTP scale (Figure 3). Contrary to expectation, sample types and origins were generally clustered together. Although physiochemical parameters play an important role in ARG prevalence and dissemination, it remains unknown whether ARGs and physiochemical parameters are significantly correlated between large-and small-scale WWTPs. In the present study, correlation analysis was carried out to understand specific relationships according to WWTP scales. However, there were no significant differences between largeand small-scale WWTPs. Physiochemical parameters such as BOD, COD and SS were positively correlated (p < 0.05) with the majority of ARGs in WWTP influents, while negative correlations were observed in effluents. Previous studies reported that pH, COD, and TN/P were positively correlated with ARG abundance in both influent and effluent from municipal WWTPs, and they may be important factors affecting the bacterial community during the treatment process [41,42]. Because raw influent wastewater contains massive amounts of organics and excess nutrients compared to effluent, potentially providing selective pressure for ARB and ARGs [6,12], physiochemical parameters may be positively correlated with major ARG types. Although physiochemical parameters were different according to large and small scale (Table 1), there was no statistically significant difference (p > 0.05) between them, and the overall treatment performance (removal rate from influent to effluent) was not quite different. These results suggested that ARG abundance may not be affected by the WWTP scale, but could be affected by influent characteristics in WWTPs. Some previous studies have shown that selected physiochemical parameters are often weakly correlated with ARG abundance in WWTPs, and abundance of ARGs does not rely on operation conditions such as temperature, HRT, SRT and MLSS [41][42][43][44]. Therefore, more studies are required to resolve and determine the underlying physicochemical, biological, and operating drivers for the dissemination of ARG patterns between large and small WWTPs.
Aminoglycoside (aadA, aadA1 and aadA2) and beta-lactam (blaVEM and blaOXA) resistance genes are frequently detected ARGs in WWTPs and are closely related to antibiotics important for human infection treatment or veterinary usage [21]. These characteristics possibly explain why tetracycline, MLS, sulfonamide and multidrug resistance genes were mostly high in abundance in the 12 metagenomes.

Bacterial Community Difference between the Large and Small WWTPs
A total of 13 bacterial phyla were identified among all WWTP samples. The profile of bacteria at the phylum level showed that Proteobacteria and Actinobacteria were the most dominant bacteria (Figure 4a). Proteobacteria accounted for 10.34% to 20.86% of the bacterial community in the influent, while abundances were significantly higher, ranging from 30.32% to 70.12%, in the activated sludge and effluent from all samples (p < 0.05). At the same time, the proportions of Firmicutes significantly decreased from the influent (28.86-39.64%) to the activated sludge (3.54-4.21%) and effluent (6.84-9.03%). Among the Proteobacteria, β-proteobacteria was the most predominant with a range of 30.9~50.45% among all types of samples between large and small WWTPs, followed by α-proteobacteria for 23.8~43.8%, γ-proteobacteria for 13.1~23.6%, and δ-proteobacteria for 3.7~10.2%. These characteristics were not significantly different between large and small WWTPs. This result is similar to previous results of bacterial communities in most activated sludge studies [6,45], in which Proteobacteria was also the most dominant community. Proteobacteria are generally known to have an important role in the metabolic capacity of degrading organic pollutants in bioreactors, such as nitrogen, phosphorus, and aromatic compounds [45,46]. β-and γ-Proteobacteria are abundant and degrade nutrients such as N, and P in activated sludge of denitrifying reactors [6,45,46]. α-Proteobacteria, nitrite-oxidizing bacteria, and denitrifying capacity were also previously recognized as important members of activated sludge microbial communities [6,45]. Actinobacteria play a role in enhanced biological phosphorus removal systems [47], and high abundances of Actinobacteria were also reported in previous WWTP systems [48,49]. The most dominant genera also showed similar patterns among sample types (Figure 4b).
Bacteroides, Mycobacterium, and Clostridium were represented by 60~70% of all influent samples. Although Mycobacterium (20~60%) still represented the major genera between all activated sludge and effluent samples, shifts in community composition were observed. Burkholderia, Dechloromonas, Nitrospira, and Pseudomonas were the major dominant genera in activated sludge and effluent samples regardless of the WWTP scale. These genera are usually detected with high abundance in activated sludge processes in global WWTPs [45]. In fact, members of Mycobacterium are among the foaming bacteria in activated sludge due to their high hydrophobic cell surface and enrichment in the foam [50]. Therefore, Mycobacterium is known to have good degradation capabilities for aromatic hydrocarbons and nitrogen-containing heterocycles [50]. Nitrospira bacteria are the key nitrite oxidizers in sewage treatment plants [51,52]. The Dechloromonas and Burkholderia genera are known to have an important role in denitrification and the ability to metabolize aromatic compounds [53][54][55]. Among the dominant genera, Mycobacterium and Pseudomonas are opportunistic pathogen-associated genera, and their detectable levels in effluents may pose problems in receiving waters or may enable regrowth in storage tanks.
The results of this study suggest that the microbial communities are generally different according to the capacity of WWTPs, without some exceptions, such as the dominance of Mycrobacterium genera between large-and small-scale WWTPs. This is consistent with previous studies [10,56,57]. Although there are many reasons, similar fitness levels of the microbial community from the influent wastewater community may fit the bioreactor despite different WWTP scales. If specific niches are dominant, such as Mycobacterium, and the dispersal rate is not too high, influent wastewater communities may be assembled in bioreactors. Samples with a relative abundance <1% were classified as "others".

Variation of the Abundance of MGE
As shown Figure 5, the alignment of all the reads to known plasmids showed that large-scale WWTP samples had the highest abundance of all plasmids, which was 1.5-2.0-fold higher than that of the small-scale WWTP samples. In the total numbers of the C and D WWTP samples, known plasmid copies encompassed a range of 0.58 ± 0.16. In contrast, the relative abundance of known plasmids in the A and B WWTP samples was observed to decrease to low levels of 0.33 ± 0.05. Alignment against known integrons showed approximately 1.5-fold and 3-fold higher relative abundances of integronassociated integrases in the large-scale WWTP samples (0.03-0.04%) than in the smallscale WWTP samples (0.008-0.02%) ( Figure 5). The class 1 integronase gene intI1 was predominant in all samples. More types of integrase genes were detected as intI, intI3, and unknown integrase genes, and their encoding gene cassettes were also detected in this study (data not shown). The ARG diversity and abundance of both large-and small-scale WWTPs showed a strong correlation with those of the plasmids, showing R 2 values of 0.66 to 0.84, respectively (Table 2). Considering the different ARG profiles of each sample, the similarity in diversity indices for the parameters of antibiotic resistance suggests that plasmids can aid the persistence of ARGs under different environmental conditions. Plasmids conferring antibiotic resistance are known to be stable in the environment even without antibiotic selection pressures [12,17,58]. In addition, a strong correlation between ARG and integron abundance and diversity, showing an R 2 value of 0.72 to 0.90, was also detected in this study (Table 2). Integrons are suggested to contribute to the exchange and incorporation of ARGs, resulting in the proliferation of bacterial antibiotic resistance in WWTPs [12,17,58]. HGT via MGEs easily occurs in niches with high biomass and extracellular stress, such as sediments [59], activated sludge [60], and anaerobic digestion [61]. Previous studies reported that tet and sul resistance genes have a strong association with MGEs to move from one species to another [21]. tet resistance genes are usually highly abundant and known to be transferred among bacteria into the receiving environment via HGT [62]. Considering the transfer mechanisms of sulfamethoxazole, sul1 is carried on intI, but sul1 has been detected on a broad range of host plasmids; this could lead to widespread detection of sul1 in the wastewater environment [63,64]. The intI1 gene had intensive connections with multiple ARGs, which can potentially confer resistance to multiple classes of antibiotics. This result was also found in previous studies [8,31,58], thereby suggesting that integrons have widespread coexistence of the same or very similar ARG clusters in various environments because integrons can capture multiple gene cassettes to facilitate the coexistence of ARGs [12,58].

Conclusions
This study provides a baseline for the different types of antibiotic resistomes and MGEs between small-and large-scale WWTPs. As expected, large-capacity WWTPs had a relatively higher abundance and occurrence of ARGs and MGEs than small-scale WWTPs due to the high input of organic compounds and biomass, which provide favorable environments for the proliferation of ARGs. However, similar taxa structures between large-and small-scale WWTPs assume that specific WWTP niches correspond to regional environmental conditions. Therefore, more studies should be carried out to understand bacterial taxa according to the WWTP scale and their occurrence trends to control the potential relationship between ARG emergence and removal in WWTPs in the future. In addition, the actual relationship between the bacterial community and HGT according to the WWTP scale should be determined in future studies.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.