Metagenomic Analysis Reveals the Fate of Antibiotic Resistance Genes in a Full-Scale Wastewater Treatment Plant in Egypt

Wastewater treatment plants (WWTPs) are recognized as hotspots for the dissemination of antibiotic resistance genes (ARGs) and antibiotic-resistant bacteria (ARBs) in the environment. Our study utilized a high-throughput sequencing-based metagenomic analysis approach to compare the ARG abundance profiles of the raw sewage, treated effluent and activated sludge samples from a full-scale WWTP in Egypt. In addition, the difference in microbial community composition due to the treatment process was assessed. As a result, 578 ARG subtypes (resistance genes) belonging to 18 ARG types (antibiotic resistance classes) were identified. ARGs encoding for resistance against multidrug, aminoglycoside, bacitracin, beta-lactam, sulfonamide, and tetracycline antibiotics were the most abundant types. The total removal efficiency percentage of ARGs in the WWTP was found to be 98% however, the ARG persistence results indicated that around 68% of the ARGs in the influent could be found in the treated effluent. This finding suggests that the treated wastewater poses a potential risk for the ARG dissemination in bacterial communities of the receiving water bodies via horizontal gene transfer (HGT). The community composition at phylum level showed that Proteobacteria, Bacteroidetes, Firmicutes, and Actinobacteria were the most abundant phyla in all datasets. Although the relative abundance of several pathogenic bacteria in the influent declined to less than 1% in the effluent, the taxonomic assignments at species level for the effluent and sludge metagenomes demonstrated that clinically important pathogens such as Escherichia coli, Klebsiella pneumonia, and Aeromonas caviae were present. Overall, the results of this study would hopefully enhance our knowledge about the abundance profiles of ARGs and their fate in different wastewater treatment compartments that have never been examined before.


Introduction
The antibiotic resistance phenomenon was addressed by the World Health Organization (WHO) and the European Center for Disease Prevention and Control (ECDC) as one of the major public health problems worldwide in the second decade of the current century [1]. In addition, the review report on antimicrobial resistance published by the UK government in collaboration with Wellcome Trust in May 2016 indicated that the rates of infectious diseases increased in low and middle-income countries with inefficient sanitation systems and safe water sources [2]. Furthermore, the extensive use of antibiotics to treat infections led indirectly to the development of microbial drug resistance within these populations. Therefore, the effective management of the utilization of antibiotics in disease treatment as well as reducing their disposal into the environment would play an important role in the mitigation actions against the spread of antibiotic resistance nationwide [2]. Sewage represents a model for an anthropogenically-impacted ecosystem. It involves human-associated and environmental bacteria coupled with trace concentrations of pharmaceutical residues, biocides, and heavy metals in one medium. The elevated toxin levels in WWTPs can drive microorganisms to acquire ARGs via horizontal gene transfer (HGT) as part of their defense mechanism for survival [3,4]. As a result, ARGs could be encoded by new species leading to the spread of antibiotic resistance phenotype in the environment [5].
Wastewater treatment technologies are classified into two main categories: conventional and non-conventional (advanced) [6]. Conventional wastewater treatment is relatively economical and highly efficient in the removal of organic substances and a wide range of pathogenic bacteria. However, it fails to adequately eliminate ARBs and ARGs from wastewater [7]. The conventional wastewater treatment consists of primary, secondary and sometimes tertiary phases [8]. Primary stage is mainly mechanical that aims to reduce solid particles such as grease, grit, and sand. However, the secondary treatment is a biological phase aims at reducing organic matter content using aerobic or anaerobic microbial degradation [9]. On the other hand, non-conventional wastewater treatment uses advanced technologies such as membrane bioreactors (MBR), moving bed biofilm reactor (MBBR), fixed bed bioreactors (FBR), and new class of nanomaterials. Although, the non-conventional methods generate a better-quality treated wastewater, they are expensive and not widely spread worldwide [9].
Since 99% of the bacteria in different ecosystems are unculturable [10], the culturebased methods have a limited capacity to reflect a realistic view for the ARG occurrence in environmental samples. Molecular-based approaches such as high throughput DNA sequencing and metagenomic analysis became the best alternative for a comprehensive ARG profiling in diverse environments [11]. Multiple reports validated the power of metagenomic approach in determining the impact of wastewater treatment on the abundance and diversity of ARGs and ARBs in a variety of sites in WWTPs [12][13][14]. For example, Yang et al., [15] investigated the fate of ARGs in a full-scale WWTP using metagenomic analysis. As a result, 271 ARGs belonged to 18 different resistance types were detected along with various treatment units in the examined station. Another metagenomic analysis study of ARG profiles in aerobic activated sludge (AAS) and anaerobic digested sludge (ADS) by Guo et al., [16] concluded that some environmental bacteria such as Nitrosomonas and Clostridium were found to be potential hosts for ARGs contributing in the spread of antibiotic resistance in the downstream environments. In a recent shotgun metagenomic analysis study by Yoo et al., [17], authors investigated two activated sludge and two anaerobic digestion sludge samples from WWTPs. As a result, a total of 181 ARGs encoding for resistance against 18 antibiotic classes were detected.
Available reports for the effects of wastewater treatment on ARGs sometimes showed conflicting results due to different sampling methods and experimental designs as well as various bioinformatic analysis tools [18]. While some reports indicated that ARG abundance was reduced as a result of the treatment process [15,[19][20][21], other studies claimed that the wastewater treatment caused the elevation of ARG concentration in the environment [22][23][24][25][26].
Despite the extensive research conducted on ARGs and ARBs in WWTPs in the last two decades, there are no implementations for reference guidelines or legislations determining the minimum permissible concentration of ARGs in the environment [27].
Moreover, the evidenced prevalence of ARGs in WWTPs indicates that the problem has not been tackled effectively [5]. Therefore, continuous research efforts may add a valuable knowledge towards reaching a sustainable solution for a potential health-threatening problem. The present work fills a research gap by providing a novel source of a publicly available metagenomic datasets from a formerly uninvestigated geographic area.
This study aims to identify the abundance profiles of ARGs in the raw influent, treated effluent and activated sludge of a full-scale WWTP using metagenomic analysis. In addition, assessing the impact of treatment process on ARG elimination by calculating the ARG persistence and removal efficiency percentages in the studied WWTP. Moreover, comparing the estimated ARG abundance profiles to their counterparts of other selected environments. Finally, identifying the microbial community composition at multiple taxonomic ranks for the collected metagenomic datasets.

Sample Collection and Processing
Tezmant wastewater treatment plant (Tz-WWTP) is located in Beni Suef governorate, Egypt ( Figure 1). It serves a population of about 296,537 with a mean daily inflow rate 55,000 m 3 [28]. To avoid seasonal variations in sewage content, six sampling events were made on monthly basis in winter (28 December 2016, 29 January and 28 February 2017) and summer (29 June, 30 July, and 29 August 2017). Winter influent samples were collected each in 2-L volumes in a sterile glass bottle and labeled as INF_W1, INF_W2, and INF_W3. Those for effluent were collected each in 10 L sterile polypropylene tanks and labeled as EFF_W1, EFF_W2, and EFF_W3. Similarly, summer samples were labeled as INF_S1, INF_S2, INF_S3, EFF_S1, EFF_S2, and EFF_S3. Due to the stability of microbial composition of the activated sludge as previously reported [15], only one winter and one summer samples were collected during January and July 2017 sampling events respectively. Samples were collected in sterile 50 mL falcon tubes and labeled as SL_W and SL_S Totally, fourteen manual grab samples were collected as six samples of the raw influent (INF), six samples of the treated effluent (EFF) and two samples of activated sludge (SL) (Figure 1). All samples were fixed at the sampling site by adding one part of sterile TE buffer (concentration?) to five parts of each sample in volume/volume ratio [29]. Then, samples were transported to the laboratory in an icebox for processing within four hours from sampling.
In the laboratory, pellets from INF and SL samples were collected by centrifugation at 9000 rpm for 15 min on refrigerated centrifuge (Eppendorf ® , Model 5810R, Hamburg, Germany). On the other side, EFF samples were filtered through MF-Millipor Membrane Filter, 0.1 µm pore size (catalog number VCWP14250, Merck-Millipore, Burlington, MA, USA) using Millipore Masterflex I/P peristaltic pump (Cole-Parmer Co., Vernon Hills, IL, USA) with 142 mm stainless steel filter holder. Pellets and membranes were stored at −80 • C until DNA extraction.

DNA Extraction
Genomic DNA extraction from INF and SL samples were conducted using the Power-Soil DNA isolation kit (catalog number 12888-50, MO BIO Laboratories Inc., Carlsbad, CA, USA) according to the manufacturer's protocol with some modifications [29]. For each EFF sample, the filter membrane, was cut into small pieces in a 50 mL sterile tube and covered by the modified TE buffer solution (10 Mm tris pH 8.0 10 mM EDTA) before vortexing for a short time to mix using vortex mixer. The remaining extraction steps for cell lysis and DNA extraction were performed as done for INF and SL samples.
To avoid any variation in the DNA extraction for samples from the same location and due to the limited sequencing funds, DNA from the three-winter influent samples (INF_W1,  INF_W2  ARG-like sequences were extracted and classified using ARGs-OAP v.2.0 software pipeline with Structured Antibiotic Resistance Genes reference database version 2 (SARG 2) [30]. SARG v2.0 represents one collection of 12,307 manually curated non-redundant ARG sequences derived from the antibiotic resistance genes database ARDB [31], the comprehensive antibiotic resistance database CARD [32], and National Center for Biotechnology Information Non-Redundant Protein Sequence Database NCBI-NR version on 21 July 2016. These sequences are hierarchically organized into 24 ARG classes and 1208 sub-classes respectively. The detection of ARG-like sequences in all metagenomic datasets started by a local pre-screening phase to extract the potential ARG-like sequences and 16s rRNA sequences using Ublast [33] search against SARG2 and Greengenes [34] reference databases respectively. Then, the extracted ARG-like sequences were annotated and classified using the online galaxy instance of ARGs-OAP v2.0 software (http://smile.hku.hk/SARGs, accessed 19 April 2019). The software parameters for ARGs annotation were adjusted to E-value cutoff ≤ 10 −5 , minimum sequence identity 90% and the minimum alignment length 25 amino acids. The abundance ratios of ARG types and subtypes across all metage-nomic datasets were expressed in three different normalization methods including the number of ARG copies to the number of 16S rRNA gene copies extracted from the metagenomic dataset (ARGs copies/16s rRNA gene copies), ARG copies to the number of the universal single-copy marker genes (USCMG) detected in the studied samples (ARGs copies/cell number copy) and the number of ARGs sequences in one million sequences of the metagenomic dataset (ppm). For simplicity, ARG abundance was reported using ARG copies/16s rRNA gene copies format throughout the study. Principal component analysis (PCA) was conducted to compare the ARG abundance profiles of the current study with those of other selected environments using R ver. 3.5.2 (20 December 2018).
The removal efficiency percentage of ARG types and subtypes from influent to effluent at Tz-WWTP was calculated using the following formula [13]: (1) where: (Data for VSS were provided by Tz-WWTP operation authority). Abundance1 is the abundance ratio for ARG type or subtype/16s rRNA genes of the influent samples.
Abundance2 is the abundance ratio for ARG type or subtype/16s rRNA genes of the effluent samples.
Taxonomic profiling was determined using MetaPhlAn 3.0 software [35]. MetaPhlAn3 accepts shotgun metagenomic reads as input and generates a list of the detected microbes classified into different taxonomic clades to the species rank associated with their relative abundances. The Graphical Phylogenetic Analysis software (GraPhlAn) [36] was used for the phylogenetic tree visualization of the examined metagenomes. All metagenomic datasets have been deposited to the Rapid Annotation using Subsystems Technology for Metagenomes server, MG-RAST [37], and accession numbers for all metagenome were listed in Table 1. Heatmaps of ARGs abundance profiles, taxonomic profiles, and Venn diagram were created using complexheatmaps, ggbiplot and vegan packages in RStudio package [38] for R version 3.5.2 (20 December 2018).

Metagenomic Sequencing
The total number of the high-quality sequences ready for downstream bioinformatics analysis was around 661 million reads with an approximate base pair size of 99 Gb. The number of reads per sample ranged from 81-134 million (Mean:110 ± 27 million reads) with a size range of 12-20 Gb (Average: 16.5 ± 4 Gb) Table 1.
The rarefaction curve created by MG-RAST for the read counts against the number of species retrieved from each sample showed a saturation tendency indicating that the sequencing depth is sufficient to reflect a satisfactory picture for the microbial composition of the examined environments ( Figure S1).

Profiling of ARGs
The mean percentages of the total number of ARG-like sequences derived from influent, effluent, and activated sludge datasets were 0.53%, 0.12%, and 0.06% respectively (Table S1). Influent showed the highest total average abundance of ARGs types among all datasets with a value of 3.39 ARGs/16S rRNA (Table 2). This number was reduced by more than 2 folds in the effluent to 1.55 ARGs/16S rRNA and by almost 3 folds in the activated sludge to 1.17 ARGs/16S rRNA ( Figure 3A).
In total, 578 non-redundant ARG subtypes belonging to 18 ARGs types were detected in all metagenomic samples (Table S3). The percentages of the most prevalent ARG subtypes with a minimum average abundance score of 0.01 ARGs/16s rRNA genes in any sample constituted 86%, 86%, and 89% of the total ARG subtype average abundance in the influent, effluent, and activated sludge samples respectively and covered 15/18 of the total number of the identified ARG types (Figure 4). (ARG subtypes with minimum abundance ratio/16s rRNA were selected for the heatmap analysis and scaled to Zscore values).
As part of the analysis results of all metagenomic datasets using ARGs-OAP V2.0 pipeline, the abundances of the detected ARG subtypes normalized to ARG copies/16s rRNA gene copies were compared to 51 other metagenomic datasets representing different environmental sources (drinking water, livestock, sediment, sewage treatment plants, and ocean) (Table S4). PCA analysis showed that INF, EFF, and SL samples were generally grouped closer to livestock samples more than sewage treatment plant along the PC1 axis ( Figure 5). Besides, the two INF samples were clustered closer to each other than the other four EFF and SL samples which were grouped together closely.

ARGs Persistence in Influent
A total of 371/549 (67.50%) ARG in the influent belonging to 17 ARG types were detected in the treated effluent and defined as the persistent ARGs ( Figure 6). Of them, 269 persistent ARGs were shared by activated sludge samples as well.
Resistance against aminoglycoside, beta-lactam, tetracycline, bacitracin, sulfonamide, chloramphenicol, and multidrug represented more than 90% of the total abundance of the persistent ARG types from influent to effluent in all metagenomes ( Figure S2 and Table 3). The total abundance percentage of the persistent ARGs per 16s rRNA gene copies accounted for 99.56% in the influent, 99.63% in the effluent, and 99.75% in the activated sludge (Table S5). The top ten most abundant persistent ARG subtypes in the influent were represented by three aminoglycoside-resistance genes including aadA encoding for aminoglycoside nucleotidyltransferase enzyme, aph(3")-I and aph(6)-I both encode for aminoglycoside phosphotransferases, two tetracycline resistance genes represented by tetA for tetracycline efflux pump and tetQ for tetracycline-resistant ribosomal protection protein, one sulfonamide resistant dihydropteroate synthase sul1, one bacitracin resistant undecaprenyl pyrophosphate recycling bacA, one chloramphenicol exporter gene cmlA, a multidrug inner membrane transporter mexF and unclassified gene related to cyclic AMP regulatory protein cAMP.
Moreover, several variants of clinically important ARGs were detected among the persistent subtypes in the effluent and/or sludge samples including beta-lactamase genes (CTX-M, KPC, VIM, NDM, IMP and SHV), and quinolone resistance genes (qnrS and qnrB).
Furthermore, 280 genes in the influent represented the non-persistent ARG subtypes where they were eliminated or could not be detected in the effluent and activated sludge metagenomes (Table S6). The top ten non-persistent ARG subtypes in the influent included two genes that were removed from the effluent represented by tetH encoding for tetracycline efflux protein and FOX-4 associated with beta-lactamase enzyme. Additionally, three beta-lactamases (SHV-39, CfxA3, OXA-211), a gene for membrane fusion protein associated with multidrug resistance (emrK), trimethoprim dihydrofolate reductase dfrA10, methyltransferase enzyme associated with MLS resistance ermC, and a DNA-binding transcriptional regulator gene gadX were removed from the sludge. While an integron-encoded dihydrofolate reductase (dfrB3) belonging to the trimethoprim ARG type was among the top 10 ARG subtypes in the influent that could not be found in both EFF and SL samples.

Removal Efficiency of ARG Types and Subtypes from Influent to Effluent
Since the persistent ARGs are discharged into the environment through the treated wastewater and sludge, their abundance profiles were of more concern than those of non-persistent counterparts. The overall removal efficiency percentage of the persistent ARG types from influent to effluent was 97% of the total detected ARGs per 16s rRNA gene copies. The highest removal percentage (99.20%) was associated with ARGs type conferring resistance to fosmidomycin antibiotic followed by kasugamycin ARGs type (99.16%). Multidrug resistance was the least removed type (84%) followed by trimethoprim, sulfonamide, and beta-lactam ARG types with removal efficiency percentages ranging from 96.48% to 96.88% (Table 3).
For the persistent ARG subtypes, Fox7 gene associated with beta-lactam resistance had the highest removal efficiency percentage of 99.92% (Table S7). Furthermore, Fox5, Fox8, and class B beta-lactamase genes conferring resistance to beta-lactam antibiotics in addition to 23S ribosomal RNA methyltransferase gene namely ermC encoding for resistance against MLS antibiotics had the second-highest removal percentage ranging from 99.81% to 99.90%. The removal efficiency of the rest of persistent ARG subtypes from influent to effluent ranged from 19.31% in multidrug amrB to 99.79% in beta-lactam per-5 genes. Interestingly, there were three genes smeE, smeD, and smeF conferring resistance against multidrug in the influent with low abundance scores accounting for 5.66 × 10 −5 , 2.71 × 10 −5, and 1.26 × 10 −5 ARGs/16s rRNA respectively that were enriched from influent to effluent by percentages ranging from 15% to more than 200%.

Microbial Composition of Metagenomic Datasets
The taxonomic analysis showed that Proteobacteria, Bacteroidetes, Firmicutes, and Actinobacteria were the most abundant phyla representing around 99% of the total microbial composition of the six metagenomic samples (Figure 7). The average abundance of Proteobacteria phylum in the influent, effluent and activated sludge samples were 61%, 48%, and 61% respectively (Table 4). Bacteroidetes was the second most abundant phylum with average abundances of 22%, 34%, and 27% in the influent, effluent and activated sludge respectively. At the genus level, the picture was more divergent where 200 genera were detected in all samples. The total abundance of the top 10 genera ranged from 72% to 92% in the six samples (Table S8). Aeromonas, Prevotella, and Pseudomonas from the Proteobacteria phylum and Bacteroides from the Bacteroidetes phylum were among the most prevalent genera detected in all metagenomes. The relative abundance of some genera increased from the influent to the effluent or, the activated sludge samples along with the process of wastewater treatment ( Figure 8A). For example, the abundance percentage of Pseudomonas accounted for 4% in INF_W while it jumped by 10 folds into 41% in SL_W representing the most dominant genus in the winter sludge sample. Furthermore, Archobacter and Tolumonas accounted for 6% and 0.2% respectively in INF_W whereas they had percentages of 12% and 1.7% in EFF_W. Similarly, Bifidobacterium, Rivicola and Alcaligenes genera in INF_S with abundance percentages of 0.5%, 0.6% and 0.01% respectively were enriched into 4%, 6% and 3.5% respectively in EFF_S. For the species level, a total of 542 species were detected in all samples. INF_W and INF_S had the highest number of species as 413 and 361 respectively. While these numbers declined into 162, 209, 107, and 77 for EFF_W, EFF_S, SL_W, and SL_S respectively. Only a total of 27 species were assigned with a minimum abundance of 1% in any sample ( Figure S3). The most abundant species belonged to Gammaproteobacteria class (13 species) followed by four species in Bacteroidia, three species in Betaproteobacteria, two species in each of Clostridia and Bacteroidetetes_unclassified in addition to one species in each of Epsilonproteobacteria, Negativicutes and Actinobacteria classes respectively. This group of bacteria represented 76%, 80%, 84%, 79%, 92% and 89% of the total abundance for INF_W, INF_S, EFF_W, EFF_S, SL_W and SL_S metagenomes respectively. Aeromonas caviae and Aeromonas media from Proteobacteria phylum were the two microorganisms among the top ten most abundant species that were shared by all samples with variable percentages ranging from 2% to 18% for the first microorganism and from 2% to 8% for the second one Table S9. Furthermore, it was observed that a gut bacterium Prevotella copri from the Bacteroidetes phylum dominated both winter and summer effluent samples with percentages of 17% and 11% in EFF_W and EFF_S respectively. The distribution pattern of bacterial species throughout the wastewater treatment process from influent to effluent and sludge in Tz-WWTP were shown in Figure 8B.

Discussion
This study is the first determining a comprehensive ARGs catalog in a full-scale WWTP in Egypt using metagenomic analysis. A comprehensive profile of ARG subtypes (578 genes) encoding for resistance to multiple antibiotic classes was recovered from the investigated metagenomic samples. Similar results were obtained from a four-year study by Yang et al., [39] for the activated sludge samples from WWTP using metagenomic analysis where a wide range of ARGs conferring resistance to aminoglycoside, tetracycline, sulfonamide, multidrug, and chloramphenicol were the most prevalent classes. Although WWTPs were not originally designed to eliminate ARGs from the raw sewage, the reduction in ARG abundance from influent to effluent was observed. However, many of ARG subtypes in the influent could be detected in the effluent and activated sludge indicating that conventional wastewater treatment was inefficient in the complete elimination of around 68% of the detected ARGs from the influent in Tz-WWTP. Previous studies provided contradictory results with respect to the impact of the treatment process on the abundance of ARGs in WWTPs. Lira et al., [18] demonstrated that the abundances of ARGs, ARBs, and plasmids decreased significantly from raw wastewater to UV-treated effluent samples as a result of wastewater treatment. Yuan, Guo & Yang [40] found that 40% of the examined ARGs conferring resistance to erythromycin and 80% of those associated with tetracycline resistance were not removed from a series of chlorinated wastewater effluent samples. However, Mao et al., [41] who tracked the flow of 30 ARGs through treatment units of WWTPs demonstrated that twelve ARGs including tet(A), tet(B), tet(E), tet(G), tet(H), tet(S), tet(T), tet(X), sul1, sul2, qnrB, and erm(C) were enriched from influent to effluent during the treatment process. Divergence in ARG abundances among studies could be referred to differences in sewage sources and treatment process of WWTPs in addition to sequencing depth, size of ARGs reference database, and alignment identity cutoffs [42].
The dominance of multidrug resistance type in the influent samples was comparable to that obtained from clinical settings using susceptibility testing techniques in Egypt [43]. They indicated that 93% of the tested bacterial isolates (n = 126) from hospitalized patients exhibited multidrug resistance phenotype. Zhang, Zhang & Fang, [44] suggested that WWTPs could be one potential source for the increase of multidrug bacteria in water environments. Moreover, the prevalence of multidrug resistance in the aerobic treatment of wastewater was previously reported [39]. This was explained by the presence of a high concentration of several pollutants like heavy metals and several toxins in wastewater which might stimulate bacteria to select for resistance mechanisms against multiple drug types to maintain their survival in the contaminated environment [25].
Concerning the increase in sulfonamide abundance from influent to effluent and sludge, it could be interpreted as the sewage treatment system had a little impact on the removal of the genes belonging to this drug class [42]. Furthermore, Wang et al., [45] indicated that the high persistence of sul1 and sul2 genes in WWTPs was due to the chemical stability of sulfonamide antibiotics and their increased solubility in wastewater. In another research at the same WWTP of the current study, Younes et al., [28] demonstrated that the residual concentration of sulfamethoxazole antibiotic (sulfonamide class) represented the highest concentration in the influent and effluent samples respectively. This observation could be correlated with our results of the elevated abundance of sulfonamide ARG type in Tz-WWTP. Further investigations are needed to study the correlation between antibiotic concentration and ARG abundance in the studied WWTP.
The abundance ratios of some ARG subtypes within the influent samples showed seasonal variations between summer and winter which in turn affected the corresponding scores in the treated effluent and sludge samples (Table S3). For example, the abundance ratio of tetQ gene encoding for resistance against tetracycline antibiotic was 5.26 × 10 −2 ARGs/16s rRNA in the summer influent sample however, the value for the same gene was almost doubled in the winter influent to be 1.18 × 10 −1 ARGs/16s rRNA. Furthermore, the multidrug resistance genes in the influent mexF and mexT had higher abundance ratios in summer (1.47 × 10 −1 and 8.83 × 10 −2 ARGs/16s rRNA) more than in winter (4.15 × 10 −2 and 3.53 × 10 −2 ARGs/16s rRNA). More research is required to investigate the impacts of seasonal variations on the ARGs abundance profiles in Tz-WWTP. However, some previous reports claimed that the abundance of specific ARGs in the winter was higher than that recorded in the summer [46,47]. Harnisz et al., [48] explored the impact of seasonal variations on the prevalence of selected ARGs in the effluent of 17 WWTPs and 7 points at river catchments using qPCR. Authors indicated that blaTEM, tet(A), ermF, sul1, and aac (6 )-Ib-cr genes were more abundant in winter and spring than in summer and autumn.
PCA analysis of ARG subtypes abundances of the current study compared to publicly available data of other environmental samples indicated that metagenomes of Tz-WWTP were clustered closer to livestock metagenomes. This could be interpreted as some of the sewage sources of Tz-WWTP were coming from nearby rural areas in the Beni Suef governorate where the wastes of the household or farm animals were possibly discharged into the sewer system affecting the microbial composition and ARGs content of the raw wastewater inflow. Furthermore, the total abundance ratios for ARG subtypes in all metagenomic samples of our study which ranged from 1.09 to 3.59 ARGs copies/16s rRNA gene copy was comparable to the range obtained in a study by Li et al., 2015 for livestock wastewater samples which ranged from 5.4 × 10 −1 to 3.1 ARGs copies/16s rRNA gene copy.
The persistent ARGs in the influent represented more than 99% of the total abundance of all the detected genes in the effluent and sludge which suggested that the raw influent was the main source of ARG subtypes detected across these two sampling locations. This finding was in agreement with the results obtained previously by Yang et al., 2014. The existence of persistent ARGs associated with high clinical significance as part of the discharged effluent or activated sludge in Tz-WWTP; even if they were found in low abundances; raised a question mark about the potential emergence of these resistance phenotypes in the downstream environment and urged for additional surveillance plans.
The microbial community composition of sewage is complex and reflects not only microbes of human and animal origins but also other environmental bacteria from sewer systems [33]. A study by Zhou et al., [49] demonstrated that the abundance and diversity of the studied ARGs were driven by the shifts occurring in the microbial composition of the samples collected from highly polluted areas. Additionally, Hultman et al., [19] indicated that although the treatment process succeeded in minimizing the abundance of bacterial hosts for ARGs, some bacteria in the effluent harbored ARGs that were not carried by these groups in the influent indicating that shifts in host bacteria occurred as a result of wastewater treatment. In our study, the taxonomic composition at the phylum level of the samples collected from different treatment units was in concordance to those obtained from WWTPs in previous reports [11,15,50]. These results validated the presence of a common taxonomic profile within the dominant phyla including Proteobacteria, Bacteroidetes, Firmicutes, and Actinobacteria in various compartments of WWTPs regardless of their geographical locations [17,51]. Members of the Proteobacteria phylum were considered to be a major player in the removal of organic pollutants such as phosphates and nitrates from wastewater during the conventional treatment process [52].
The divergent results of the microbial community composition at genus level among WWTPs of different geographic locations could be referred to variations in the pore size of filters used to collect bacteria from wastewater, DNA extraction method applied, sequencing technology, and other environmental factors such as pH and temperature [53]. In the current study, the existence of some bacterial genera with known pathogenic members such as Aeromonas, Pseudomonas, Eubacterium, Laribacter, Escherichia, Klebsiella, Archobacter, and Acinetobacter among the prevalent genera in at least one sample calls for long term surveillance for the potential risks of ARGs dissemination in the pathogenic members of these groups (for example Aeromonas caviae, Aeromonas veronii, Escherichia coli, Klebsiella pneumonia, and Eubacterium rectale) at the receiving water bodies of the treatment station.
Available studies showed that members of Aeromonans and Pseudomonas genera conferred resistance against various beta-lactamases of classes A, B, C, and D, quinolones in addition to multidrug resistance [54][55][56]. In another study by Hu et al., 2016 Pseudomonas spp were found to be potential hosts for ARGs encoding for resistance-nodulation-division (RND) family transporters which are part of bacterial efflux pumps. The high abundance of Pseudomonas in the activated sludge samples of the current study was in agreement with a previous study by Moura et al., 2010 indicating that bacteria of Pseudomonas genus were known to be important players in the microbial degradation of polyphosphates in activated sludge of both municipal and industrial WWTPs.
Our results for genera associated with fecal bacteria such as Bacteroides and Prevotella were in concordance with a recently published study by Quintela-Baluja et al., [57] who compared microbial community composition and ARGs across different treatment compartments in both domestic and hospital WWTPs. Authors indicated that the levels of coliform bacteria were relatively higher in the effluent than in the activated sludge suggesting that liquid wastewater had a more serious impact on the downstream environment regarding fecal contamination and possible risks to human health. The high species abundance of a human gut microbe Prevotella copri in the raw influent at Tz-WWTP could be explained by its prevalence in the stool of the inhabitants in the Beni Suef governorate, Egypt. Likewise, Shankar et al., [58] showed that Prevotella was the most prevalent genus in the studied gut microbial communities of Egyptian teenagers while Bacteroides was the dominant genus in the gut microbiome of the children from the USA. Authors indicated that the shifts that occur in the microbial communities were related to the differences in diets of Mediterranean and Westernized populations while fibers and polysaccharides were the major components in the first group, proteins and fats represented the main part in the second one. Moreover, other scientific reports indicated that the abundance of Prevotella copri in gut microbiota was associated with some human diseases as rheumatoid arthritis (RA) [59] ankylosing spondylitis [60], and insulin resistance [61].
The capacity of ARGs to mobilize among community members through HGT is driven by the co-existence of resistance genes with a mobile genetic determinant (plasmid, transposon, or integron) harbored in the same genetic context and hosted by pathogenic bacteria [62]. Therefore, the current profiling study provides limited information concerning the ARGs mobility potential in Tz-WWTP suggesting that further study is required to investigate the genomic context of the detected ARGs in Tz-WWTP and their possible hosts.

Conclusions
The research results of the current study confirmed the power of metagenomic analysis to reveal a comprehensive profile of ARGs across different treatment stages in a full-scale WWTP. A total of 578 ARG subtypes belonging to 18 different ARG types were detected in INF, EFF, and SL samples of Tz-WWTP. Although the total average ARG subtype abundance ratios declined as a result of the wastewater treatment process in Tz-WWTP, 371 ARG subtypes persisted in the chlorinated effluent and sludge samples posing potential risks of ARGs spread in the receiving environment. The abundance ratios of some ARG subtypes within the influent samples showed seasonal variations between summer and winter.
Some ARGs of clinical relevance such as beta-lactamase genes (CTX-M, KPC, VIM, NDM, IMP, SHV), and quinolone resistance genes (qnrS, qnrB) were among the persistent ARGs in the effluent and sludge metagenomes at Tz-WWTP.
The microbial community composition of all examined samples in Tz-WWTP had a common core of abundant phyla including Proteobacteria, Bacteroidetes, Firmicutes, and Actinobacteria. At the class level, Gammaproteobacteria, Bacteroidia, Betaproteobacteria, Clostridia, Bacteroidetes-unclassified, Epsilonproteobacteria, Negativicutes, and Actinobacteria represented the top prevalent taxa in Tz-WWTP. Concerningly, several pathogenic bacteria were among the assigned species in the effluent and sludge samples despite the treatment process. These bacteria included Aeromonas caviae, Aeromonas veronii, Escherichia coli, Klebsiella pneumonia, and Eubacterium rectale raising a research question for the possible contribution of these pathogens to the spread of ARGs in the receiving water bodies.
In fact, some difficulties were met during the computational analysis due to the inadequate hardware resources suitable for the processing of large size datasets. Furthermore, the limited number of samples was challenging in statistical analysis. Moreover, the analysis of short reads (150 bp each) has a limited potential to retrieve complete ARG sequences with their genomic context. Therefore, further studies are encouraged to extend this work by analyzing the generated datasets using metagenomic assembly. Consequently, the ARG mobility potential can be examined to estimate the risk level of the detected ARGs to human health.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/su132011131/s1, Figure S1: Rarefaction curve for the species count per metagenmic sample against the number of sequenced reads as estimated by MG-RAST analysis., Figure S2: Average abundance of the persistent ARG types from influent to effluent and activated sludge metagenomes, Figure S3: Taxonomic cladogram representing the top most abundant species present in metagenomic samples, Table S1: Summary of the total abundance ratio of the extracted ARG sequences from all metagenomic datasets using ARGs-OAP v2.0, Table S2: Comparison of ARGs diversity across ARG types in the Influent, Effluent, and Activated sludge samples, Table S3: ARGs subtype abundance per 16s rRNA gene copies in the influent, effluent and activated sludge samples, Table S4: PCA matrix for ARG subtype abundance in Tz-WWTP compared to other environmental samples as indicated by ARGs-OAP V2.0 analysis, Table S5: Persistent ARG subtypes from influent (INF) to effluent (EFF) and sludge (SL)in Tz-WWTP, Table S6: Non-Persistent ARG subtypes from influent (INF) to effluent (EFF) and sludge (SL) in Tz-WWTP, Table S7: Average abundance of the top 100 persistent ARG subtypes with their removal effeciency percentages from Influent to Effluent metagenomes, Table  S8: The relative abundance of taxonomic clades at genus level, Table S9: The relative abundance of taxonomic clades at species level with minimum abuancance of 1% in any metagenomic sample in Tz-WWTP.