Metagenomic and viromic data mining reveals viral threats in biologically treated domestic wastewater.

Activated sludge (AS), a common biological secondary treatment process in wastewater treatment plants (WWTPs), is known to remove a large spectrum of microorganisms. Yet little is known about its effect on the entire viral community. After compiling 3 Tbp of next-generation sequencing (NGS) metagenomic/viromic datasets consisted of 119 sub-datasets of influent, effluent, and AS samples from 27 WWTPs, viral removal efficacy is evaluated through data mining. The normalized abundance of viruses suggests effluents exhibit the highest viral prevalence (3.21 ± 3.26%, n = 13) followed by the AS (0.48 ± 0.25%, n = 57) and influents (0.23 ± 0.17%, n = 17). In contrast, plasmids, representing genetic element of bacteria, show higher average prevalence (0.73 ± 0.82%, n = 17) in influents than those of the AS (0.63 ± 0.26%, n = 57) and effluents (0.35 ± 0.42%, n = 13). Furthermore, the abundance-occupancy analysis identifies 142 core phage viruses and 17 non-phages core viruses, including several pathogenic viruses in the AS virome. The persistent occurrence of pathogenic viruses, coupled with non-favorable virus removal by the AS treatment, reveals the hidden virus threats in biologically treated domestic wastewater. The mechanisms for why viruses persist and the possibility that WWTPs are potential hotspots for viral survival deserve attention.

Introduction

Activated sludge (AS), the most widely employed biological treatment by wastewater treatment plants (WWTPs), was first introduced by Arden and Lockett in Manchester, UK [1] and has played critical roles in the removal of pollutants including inorganics, bacteria, and protozoa [2]. Usually, the core design contains three parts, including an aeration tank, a settling tank, and return AS equipment. In brief, the aeration tank provides detention time and ensures that the wastewater influent and AS are well mixed. After a serial biochemical reaction, the pollutants can be removed and the influent can be converted into the effluent. Although a growing body of evidence has shown that it is capable of lowering the concentrations of bacteria [e.g., fecal indicator bacteria, (FIB)] [3,4], much less is known about viral removal efficacy despite a consensus that WWTPs are viral reservoirs, some of which are pathogenic [5]. Human pathogenic viruses identified in WWTPs include adenovirus, enterovirus, rotavirus, norovirus, and polyomavirus, etc. [6]. The majority of the DNA viruses are bacteriophages that affect WWTP performances through shaping the microbial community [7,8]. To date, efforts have been directed to the reduction of FIB and phages per water quality guidelines [9,10].

Because the aforementioned bacterial indicators do not always correlate with waterborne viruses [11], viral removal requires separate assessment with studies reporting findings that are encouraging but also concerning. Ottoson et al. found that no indicators predicted enterovirus genome occurrence in treated wastewater, but E. coli, enterococci, and C. perfringens removal were correlated to the enterovirus genome removal [12]. A one-year survey in Japan (monthly sampling) showed that the removal of noroviruses varied greatly between winter and summer [13]. Ito et al. demonstrated viral-type dependent removal during secondary treatment (physical settling and biological process) using qPCR to target noroviruses (genogroup I and II) [14]. Encouragingly, Qiu et al. found that AS treatment greatly reduced eight types of human pathogenic viruses in Canada [15]. In part due to the choice of low-throughput qPCR technologies targeting FIB, certain bacteriophages, and a few well-studied pathogenic viruses, the effectiveness of AS treatment to reduce viruses at the entire community level is unknown.

Although the bacterial proportion of the WWTP microbiome has been widely investigated by high-throughput sequencing approaches to reveal treatment-relevant community assembly patterns [[16], [17], [18]], the virome of WWTPs remains largely uncharted. It is only recently that protocols developed during the Global Ocean Viromes exploration (GOVs) [19] have been introduced to investigate the viral diversity in WWTPs. Viromic datasets can now be generated with the aid of the particular virus concentration, purification, and extraction steps [20]. However, only a handful of studies have used such sophisticated viral quantification methods to investigate virus diversities in WWTPs. Wang et al. provided the first long-term assessment of viral community dynamics of activated sludge virome [21]. They identified maximum AS viral abundance in the winter season of Hong Kong, inspiring us to profile viral prevalence patterns during biological secondary wastewater treatment of WWTPs globally. To the best of our knowledge, no large-scale metagenomic effort has been put forth to do so.

This study assesses the virus removal efficacy through profiling the core viruses pre-, during, and post-AS treatment to shed light on removal preferences at the whole community level. To this end, 119 WWTP metagenomic/viromic datasets from 27 WWTPs using conventional or modified AS treatment were first obtained from the National Center for Biotechnology Information (NCBI) and Metagenomic Rapid Annotations using Subsystem Technology (MG-RAST) public servers (See SI and Supplementary file 1 for details). The viral prevalence of sewage influents, activated sludge, and effluents were quantified using an improved metagenomics analysis method to assess virus removal efficacy. New insights gained on viral community characteristics and implications for risk management during biological secondary treatment of WWTPs are discussed.

Materials and methods

Metagenomic and viromic data

In this study, metagenomic and viromic NGS sequencing data were downloaded from NCBI and MG-RAST, consisting of 26 influent/sewage, 59 activated sludge and 34 effluent samples from 27 WWTPs located in 12 different countries encompassing Uganda, USA, 5 Asian (China, Korea, Malaysia, Singapore, Eastern category), and 5 European (Denmark, Italy, Slovenia, Switzerland, UK, Western category) countries (see Table S1, Table S2 and Supplementary file 1). The wastewater types, configurations, and loading capacity characteristics of all surveyed WWTPs were described in supporting information (SI).

Virus identifications and annotations

Except for the directly retrieved assembled contigs, we applied metaSPAdes [22] and MEGAHIT [23] (>34 Gb raw reads) to assemble post-quality-controlled reads (fastp [24] addressed) with default settings. Only contigs/scaffolds longer than 1 kb were subject to the subsequent virus identification step. Thereafter, we counted the assembly numbers and analyzed valid assemblies using a recently published metagenomic virus identification tool named metaviralSPAdes [25]. In brief, two steps named viralVerify and viralComplete were performed. The predicted genes by Prodigal [26] were first searched using hmmsearch (HMMER 3.1b2) against the Pfam-A database v30.0 [27]. Here, we only selected the “virus” category for further viralComplete processing step in which the “similarity” of a given viral contig was calculated using Naive Bayesian Classifier to each known virus from the RefSeq database. The assembly and viral contig numbers were reported in Supplementary file 1, and the viral identifications from all 119 samples were summarized in Supplementary file 2. The predicted genes by Prodigal were aligned against Swiss-Prot database using BLASTP [28] to annotate viral protein functions (Fig. S3).

Specifically, to search for the viral antibiotic resistance genes (ARGs), all viral assemblies were analyzed by ResFinder 3.2 (DNA alignment based method) using 70% identity and 60% minimum length thresholds [29]. In addition, we also used Hidden Markov Model (HMM) to identify the ARGs. Briefly, hmmsearch was used to annotate ARGs (Prodigal predicted) against SARGfam database with a gathering cutoff (--cut_ga) [30]. Thereafter, dN/dS ratio was used to evaluate the selective pressure acting on the virus-carrying ARGs using KaKs_calculator 2.0 with GMYN approximate method [31]. Briefly, MUSCLE alignment between nucleotide sequences of virus-carrying ARGs and their closest hits (against ResFinder database) were used as input for KaKs_calculator 2.0 [32].

Relative abundance calculations

Since the directly downloaded assemblies (n = 18) lacked the contig coverage values, only 101 raw read datasets were used for virus quantification calculations. Therefore, in this study for a given metagenomic/viromic dataset, the relative abundance of the virus was calculated by coverage-weighted contig-abundance using the formula below:

C is the coverage of a contig; L is the length of the corresponding contig; m is the total number of viral contigs; R is the average read length of a dataset; N is the total number of reads of a dataset.

We also calculated the relative viral abundances among WWTP locations (see Table S4) and WWTP treatment processing, respectively. Due to the specific viral concentrating and handling methods, the Slovenia virus relative abundances were much higher than the other surveyed datasets; therefore, we excluded them when comparing virus prevalence among influent, AS, and effluent (as well as performing principal coordinate analysis). We also applied the same strategies to calculate plasmid relative abundances (only targeting >1 kb fragments).

Statistical analyses

All 119 samples' viral identifications were grouped for unique and shared elements. The shared and unique viruses were illustrated with the UpSet plot and Venn Diagram (UpSetR and VennDiagram packages in R). Principal coordinate analysis (PCoA) was performed based on the Bray-Curtis distances (square root transformed) of the major viral relative abundances using “procomp” function in VEGAN package in R. To ensure confidence in the statistical significance of PCoA, the identified viruses were excluded if 1) existing in less than half of the activated sludge, influent or effluent samples, respectively; and 2) exhibiting average relative abundance less than 0.001% across all the samples. Ellipses of different groups were drawn at a confidence level of 0.68, which showed a non-overlapped clustering pattern among categories compared with the maximum number of samples included. To determine whether there were statistically significant differences among WWTP treatment stages, ANOSIM was performed using “anosim” function with 9999 permutations. To determine the viral populations that can best discriminate between WWTP virome treatment stages, we used a supervised learning technique - Breiman's Random Forest algorithm (randomForest package). Datasets from Slovenia and Malaysia were excluded for PCoA analysis among WWTP treatment stages as they showed a clear batch effect, which may be caused by the special viral enrichment procedure in Slovenia datasets and the arbitrarily-controlled sequencing-depth of Malaysia datasets (Table S1 and supplemental file1).

To determine the viruses that were cosmopolitan in the activated sludge, each viral population with an average relative abundance of more than 0.001% was drawn in the occupancy-abundance plot (Fig. 1b) for ecological partitioning to facilitate the identification of core and special viruses. Virus populations that were widely distributed in both at least 60% of Eastern (5 Asian countries) and 60% of Western (USA and European countries) activated sludge [33] were defined as ecological AS generalists.

Fig. 1

(A) boxplot shows the comparison of virus abundance (a1) and plasmid abundance (a2) for influent, activated sludge and effluent samples obtained from 19 WWTPs. Diamond represents the mean value. Asterisk represents group that showed significant (P-value < 0.05) variation. (b) Occupancy-abundance plot of viral populations in 37 Eastern (x-axis) and 24 Western (y-axis) activated sludge (AS) samples. Only viral populations showed an average relative abundance of more than 0.0010% are shown. The size of circles is proportional to the average relative abundance of the corresponding viral population. AS virus generalists (159 viral populations) are identified in at least 60% Eastern (5 Asian countries) and 60% Western (USA and European countries) AS.

Results and discussion

Relative viral abundance increases from influent to activated sludge and effluent

A total of 1779 unique virus populations (212,310 viral-contigs) were identified in 119 WWTP samples (supplementary file 1 and 2). As expected, 72.2% of the viruses were bacterial phages or prophages, confirming that phages and prophages were dominant in WWTPs regardless of locations or treatment configurations (Figs. S1 and S2) [7]. The viral relative abundance averaged 0.84 ± 1.59% (n = 87), increasing from 0.23 ± 0.17% in influent (n = 17), to 0.48 ± 0.25% (n = 57) in activated sludge and finally to 3.21 ± 3.26% (n = 13) in effluent (Fig. 1a).

For quality assurance of the metagenomic analysis, the relative abundance of plasmids which like viruses, were low in abundance, were enclosed for comparison (Fig. 1a). The influents exhibited the highest plasmids prevalence (0.73 ± 0.82%, n = 17), followed by activated sludge (0.63 ± 0.26%, n = 57) and effluents (0.35 ± 0.42%, n = 13).

Activated sludge viral generalists and occurrence of non-bacterial host viruses

Across all 57 AS samples, 159 AS major viral populations (relative abundance > 0.001%) were identified as cosmopolitan generalists widely distributed in WWTPs (Fig. 1b). A vast majority (142 out of 159) of these AS core viruses were phages, infecting 54 bacterial genera (see supplementary file 3). Those genera were commonly observed in human and animal intestinal microbiota, highlighting a tight connection between enteric microbiota and WWTP virome. Generally, the positive abundance-occupancy relationship had been recognized as a widespread macroecological distribution pattern at many scales [34]. However, such a pattern was not observed in this study, which 24 generalist viral populations were not abundant at all (mean abundance lower than the third quantile of major viruses abundance distribution, Fig. 1b), reflecting the wide variety of bacterial generalist (the major host populations) in AS [16].

Several large-sized DNA viruses affecting unicellular protozoa were also identified in AS viral generalist (the full list of AS viral generalists is in Supplementary file 3), these viruses included: 1) four giant DNA viruses infecting amoeba in various soil/water habitats (mean abundance ranged from 0.00063% to 0.0029%); 2) six algae-infecting large DNA viruses (mean abundance ranged from 0.00027% to 0.0019%); 3) one halovirus strain HRTV-5 (mean abundance 0.00040% ± 0.00053%) isolated from Archaeal genus of Halorubrum [35].

Among the 17 non-bacterial host viruses, animal and human-specific viruses were identified (Table S3): 1) three avian poxviruses affecting wild birds (with Aquatic bird bornavirus 1 showed high prevalence > 0.02%, twice in Hong Kong Shatin WWTP, Fig. S4). The mean abundance of this virus was up to 0.0027% ± 0.0046% in activated sludge samples studied; 2) one ranavirus strain (European sheatfish virus of 0.0013% ± 0.0022%) causing high mortality rates in sheatfish industry; 3) Raccoonpox virus of 0.00068% ± 0.00086%, which showed natural tumor tropism without human pathogenicity, therefore could be applied as an oncolytic virus for clinical tumor treatment; 4) one viral strain(Mocis latipes granulovirus of 0.00023% ± 0.00017%) infecting the geographically widespread moth. Such common presences of animal (avian and fish virus) and human (Raccoonpox virus) specific viruses in the biological treatment process indicated potential viral threats as host-specific waterborne viruses may cause a wide range of diseases in humans and/or other animals [36].

Additionally, we detected sporadic plausible signal of human pathogenic viruses in AS and effluent samples, including partial hits to SARS virus human coronavirus HKU1 (HCoV-HKU1) of 0.000066% in AS of Hong Kong Shatin wastewater treatment plant collected in March 2009 and of 0.00001.5%–0.000029% in the effluent of Korea Gwangju wastewater treatment plant collected prior 2018, implying the potential to use wastewater virome to monitor the spread of the virus in a community [37].

Distinct viral community in influent, AS and effluent

Despite three outliers of effluent samples based on RNA sequencing [38], the viral profiles of influents, AS and effluent showed significant separations (p < 10−4, ANOSIM) in PCoA (Fig. 2), suggesting different viral compositions. The effluents shared a higher number of viruses with the AS samples (977 viruses) than with the influents (872 viruses), indicating that activated sludge treatment could alter the virus community and could serve as a distinct viral source [39]. Based on random forest analysis, we identified nine phage populations (supplementary file 4) as the main drivers for the separations (Fig. 2) among AS, influent and effluent. Briefly, Bacillus phage, Cronobacter phage, Mannheimia phage, Mycobacterium phage, and Rhizobium phage were more prevalent in WWTP effluent, while the remaining Lactobacillus phage and Synechococcus phage were outstanding for activated sludge, and Pseudomonas phage and Roseobacter phage protruded influent samples [7].

Fig. 2

Principal coordinate analysis (PCoA) plot depicts Bray-Curtis distance between the viral compositions of samples from different WWTP treatment stages. Ellipses were drawn at a confidence level of 0.68. Virome from influent, AS and effluent cluster separately (P < 10−4, analysis of similarity, ANOSIM). Effluent outliers mentioned in the main text are enclosed in the frame with the dashed line.

The viral ARG threats

Recently, ARGs are widely studied in WWTP systems due to human health concerns [[40], [41], [42]]. The HMM searching found that a total of 1002 ARGs carrying viral contigs, and the antibiotics mainly included fusidic acid (average abundance of 0.070%), polymyxin (0.055%), tetracycline (0.045%), vancomycin (0.044%), and Quinolone (0.015%), etc. The phages Pseudomonas, Vibrio, and Synechococcus served as the main ARG-carrying populations (Fig. 3a). Using a homology-based method, we detected ten viral ARG contigs, and all of which belong to phages. Two viral contigs with multidrug resistant operon were detected (Fig. 3b), in which ARGs were closely related to the integrase activities. The finding of such integrase-enabled mobility of phage-encoded ARGs added to the increasing repertoire of virus-driven mobile ARGs operons [[43], [44], [45]], whose prevalence and persistence can now be investigated in further studies. Additionally, the low dN/dS ratio of virus-carrying ARGs identified in this study (Fig. 3b) supported a negative selection of antibiotic resistance in WWTP virome [46]. In summary, despite the low detection frequency [47,48], the persistence of ARGs on virus genome (indicated by the negative selection) coupling with the ARGs mobility enabled by viral transduction as well as integrase revealed here suggested a non-negligible influence of virus on WWTP resistome.

Fig. 3

Virus-carrying antibiotic resistance genes (ARGs). (a) Relative abundance of ARG-carrying viral populations in different WWTP treatment stages (left) and their corresponding ARGs composition (right). (b) Schematic genetic organization of multidrug resistant gene clusters on viral-contigs identified in WWTP samples. Genes on the viral-contigs are colored according to their functional categories that ARGs are colored in red; mobile gene elements in blue; phage related genes in green and other genes in gray. dN/dS ratio depicts the selective pressure of virus carrying ARGs are shown in the table.

Activated sludge treatment does not favor viral removal

The contrasting prevalence of virus and plasmids in influents, activated sludge, and effluent samples (Fig. 1) provided unequivocal evidence for effective bacterial removal by activated sludge treatment but poorer removal of viruses. Assessment of viral removal in WWTPs is challenging and fraught with methodological uncertainties. Targeting indicator viruses by qPCR, studies have shown that viral occurrences are closely related to seasons, and some suggested virus removal was independent of the WWTP capacity [13,49]. The complexity of WWTP virus removal process has also been illustrated by a weak correlation between viral removal efficiency and treatment hydraulic retention time, but the correlation with genome size, capsid protein, and disinfection strategies was stronger [50]. It is well known that physical treatment (e.g., sedimentation process) can remove some viruses from the influents, but questions remain as to how biological treatments remove viruses. Further, the biological mechanisms for viral persistence are unclear. Beyond activated sludge treatment, membrane bioreactor (MBR) is believed to be more efficient in removing phages when dealing with effluent post-AS treatment [11], but this remains controversial since the investigation included only several cultivable phages using enumeration and transmission electron microscopy (TEM) methods. The Bardenpho design, interestingly, seems to exhibit great viral removal ability for 11 frequently occurring pathogenic viruses targeted by qPCR [51], although the mechanism is again unclear. Further research is needed to examine the linkages between a wide range of wastewater treatment technologies and virus reductions, especially those with known health consequences as well as at the community level. Encouragingly, a newly published work has shown that an algal wastewater treatment system appears to effectively reduce viral presence using a metagenomic approach [52]. However, they only compared the effluent viral communities between single algal and conventional WWTP installed in one city, systematic assessment facilitated by large-scale metagenomic data mining presented here might bring more reliable insight on the viral removal efficacy of WWTPs and expand our knowledge about viral persistency during biological wastewater treatment. Taken together, the higher virus prevalence in effluents than activated sludge and influent and the viral compositional changes all imply viruses, and by extension, pathogenic viruses, are difficult to remove using the century-old activated sludge-based biological domestic wastewater treatment technology. Given the rising zoonoses, whether WWTPs act as a hotspot for virus survival deserves attention [53].

Innovations needed to shed light on virome in WWTPs

A combination of improved metagenomic sequencing and more widely used PCR or qPCR should offer a clearer picture of viromes in WWTP, the functions of virus including viral hazardousness. The detection of HAV, norovirus, and polyomavirus, even just in a few metagenomic and viromic datasets in this study, means that higher sensitivity methods such as PCR or qPCR can be applied, followed by other tests such as the serologic assays of viral functionalities [54]. Additionally, prior metagenomic/viromic studies quantified WWTPs viral prevalence directly based on the number of viral-contigs identified [55,56]. Without weighting each viral-contig by its corresponding read coverages within metagenome, it is “semi-quantitative”. While useful in delineating the virus occurrence pattern and the functional clusters of virus-encoded genes, it leads to biases in viral community composition thus cannot be used to draw meaningful conclusions about the WWTP viral prevalence. This study represents an improvement in metagenomic viral quantification by considering the different read coverages of assembled genomic fragments to calculate the relative abundance of viruses based on coverage-weighted contig abundance [21].

Conclusions

3 Tbp of NGS metagenomic/viromic datasets, consisting of 119 sub-datasets of influent, effluent, and AS samples from 27 WWTPs in 12 countries were analyzed to evaluate virome characteristics by activated sludge treatment. Our read coverage based relative abundance analysis showed that the effluents exhibit the highest viral prevalence followed by the AS and influents, which is opposite to the bacteria removal. Meanwhile, the abundance-occupancy analysis identified 142 core phage viruses and 17 non-phages core viruses, including several pathogenic viruses in the AS virome. The persistent pathogenic viruses, coupled with non-favorable biological virus removal, reveal the viral threats in global domestic wastewater. We speculate that the chlorine disinfection may lead to the relative abundance increases, which is supported by our coverage-weighted contig-abundance analysis. Previous research also showed the rotaviruses still exhibit higher persistence and strong infectivity after UV and free chlorine disinfection [57,58], which highlights the significant impact of our results. Further studies are guaranteed to explain 1) why viruses persist and 2) the possibility that WWTPs are potential hotspots for viral survival.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.