Metagenomic insights into viral community restructuring and microbial interactions in waste stabilization ponds

Waste stabilization ponds (WSPs), also referred to as wastewater lagoons, are widely implemented as nature-based, cost-effective technologies for municipal wastewater treatment^1,2. By exploiting long hydraulic retention times, solar irradiation, algal–bacterial coupling, and food-web interactions, WSPs achieve contaminant removal with minimal energy input, making them particularly suitable for regions characterized by abundant land and high solar exposure³. Configured sequentially as anaerobic, facultative, and maturation ponds, these systems progressively reduce organic matter, nutrients, and pathogens through interacting physicochemical and biological processes, including sedimentation, solar-driven inactivation, elevated pH associated with algal photosynthesis, and biological regulation^4,5.

WSP technology remains globally prevalent, with extensive implementation across North America, Europe, and Australia^6,7,8. Their resilience and operational simplicity have supported widespread adoption, particularly in developing and peri-urban settings³. However, while bacterial and protozoan reduction in WSPs has been extensively characterized, viral ecology across the treatment gradient, particularly in maturation ponds, remains comparatively underexplored⁹.

Viruses constitute one of the most abundant biological entities in aquatic systems, often occurring at concentrations 10–1000-fold higher than other microorganisms^10,11. Beyond their public health relevance, viruses may influence microbial community structure and nutrient cycling, lysogeny, and horizontal gene transfer, thereby shaping microbial community structure and nutrient cycling. In engineered anaerobic digesters, viral dynamics have been shown to explain up to 41% of microbial community variance, exceeding the contribution of abiotic factors¹². Advances in shotgun metagenomics have expanded the known virosphere, revealing extensive viral diversity across environmental systems^12,13, including the discovery of large bacteriophages and giant eukaryotic viruses^14,15. These developments underscore the importance of virome-resolved analyses for understanding ecosystem-level processes.

In wastewater treatment plants (WWTPs), viral communities are typically dominated by bacteriophages, with comparatively fewer detected human-pathogenic viruses¹⁰. Phage–host associations can influence microbial stability, organic matter turnover, and treatment efficiency¹⁶. However, WSPs differ fundamentally from mechanized systems in that they impose extended retention times, intense solar exposure, oxygen stratification, diel pH fluctuations, and enhanced sedimentation¹⁷. These features likely create strong ecological filtering pressures on viral populations, particularly in downstream units.

Among the sequential compartments of WSP systems, the maturation pond represents the terminal polishing stage, where effluent quality and pathogen attenuation are finalized^4,17. Viral reduction in WSPs results from synergistic mechanisms, including ultraviolet irradiation, reactive oxygen species generation, elevated temperature and pH, adsorption to particulates, sedimentation, and biological interactions^9,17. Maturation ponds, in particular, integrate these stressors with intensified food-web interactions and particle aggregation, potentially functioning not only as viral attenuation zones but also as ecological reactors where viral communities are reorganized^4,17.

Microbial trophic dynamics may further influence viral persistence. Grazers such as ciliates and nanoflagellates can suppress microbial populations and ingest virus-like particles¹⁸. Previous network-based analyses in WSPs have reported negative associations between protozoa and viral or bacterial taxa, suggesting top-down regulatory effects⁴. Nevertheless, virus removal efficiency varies considerably between systems¹⁹, and reliance on limited qPCR targets or fecal indicators may obscure broader virome restructuring^20,21. A community-resolved perspective is therefore needed to evaluate apparent viral reduction and ecological reassembly.

Ecological network approaches provide a complementary framework for investigating intra- and inter-domain associations among viruses, bacteria, and microeukaryotes^19,22,23. While correlation-based networks are inherently hypothesis-generating^24,25, they can reveal modular organization, co-variation patterns, and candidate trophic linkages. Positive associations often reflect shared environmental filtering or host–phage coupling, whereas negative associations may indicate niche divergence, antagonism, or grazing-linked dynamics^23,26. Integrating virome profiling with cross-domain network analysis can therefore help explore how viral relative abundance changes intersect with microbiome restructuring along treatment gradients.

Building on our previous characterization of DNA viral communities in the same WSP system, which focused primarily on the inlet and facultative pond⁹, the present study extends this framework to the maturation pond and integrates matched bacterial and microeukaryotic datasets to investigate viral and cross-domain association networks across the treatment gradient. In this study, we applied shotgun metagenomics to characterize DNA viral communities across the full treatment train—including the maturation pond—of a full-scale WSP system in Giza, Egypt. By integrating virome profiles with matched bacterial (16S rRNA gene) and microeukaryotic (18S rRNA gene) datasets, we aimed to (i) characterize stage-resolved restructuring of viral communities from influent to maturation pond, (ii) evaluate taxa-specific patterns of relative abundance decrease and downstream enrichment, and (iii) assess intra- and cross-domain co-occurrence network organization across the treatment gradient. Through this framework, we provide a community-level perspective on how natural wastewater treatment systems function as both selective viral barriers and zones of ecological restructuring.

WSP description and sample collection

This study was conducted at a full-scale WSP system located in Atfih City, Giza Governorate, Egypt (29°27′32.6″ N, 31°14′37.0″ E). The plant is designed for a nominal capacity of 45,000 m³ d⁻¹, with typical operating flows ranging between 5,000 and 7,000 m³ d⁻¹. The treatment train comprises anaerobic, facultative, and maturation ponds arranged in series, providing an overall hydraulic retention time of approximately 35–38 days. Treated effluent is reused for tree irrigation in accordance with the Egyptian code for safe agricultural reuse of treated wastewater⁴.

To resolve stage-dependent viral restructuring, water samples were collected monthly over one year from four treatment stages: inlet, anaerobic pond, facultative pond, and maturation pond outlet^4,9. Samples were collected at approximately 30 cm below the water surface using sterile containers, immediately placed on ice, and transported to the National Research Centre laboratories within 5 h. For each sampling event, 1 L of water was allocated for molecular analyses and 3 L for physicochemical characterization. For integrative multi-domain analysis, a subset of 12 representative samples was selected for shotgun metagenomic sequencing of DNA viral communities (n = 4 per stage for inlet, facultative, and maturation; one per season), no technical sequencing duplicates were included. Importantly, the corresponding DNA extracts from these same samples were used for 16S and 18S rRNA gene amplicon sequencing, enabling direct cross-domain comparisons of viruses, bacteria, and microeukaryotes across identical spatial and temporal conditions. This design enabled explicit comparison of upstream compartments with the maturation pond as the terminal ecological unit.

Physicochemical analysis

In situ measurements of temperature, dissolved oxygen (DO), and DO saturation were performed using an AD360 DO meter (Adwa Instruments, Europe). pH, electrical conductivity (EC), and total dissolved solids (TDS) were measured using Jenway meters (models 3510 and 4510). Laboratory analyses were conducted according to standard methods^27,28, to quantify biochemical oxygen demand (BOD), chemical oxygen demand (COD), ammonium nitrogen (NH₄–N), nitrate nitrogen (NO₃–N), nitrite nitrogen (NO₂–N), total Kjeldahl nitrogen (TKN), total nitrogen (TN), and total phosphorus (TP). The detailed physicochemical dataset for this full-scale WSP system has been reported in our previous publication²⁹; therefore, these data are not re-presented in the current manuscript.

DNA extraction and amplicon sequencing

Water samples were filtered through sterile membrane filters, which were aseptically sectioned and subjected to DNA extraction using the DNeasy PowerLyzer PowerSoil Kit (Qiagen, Germany), optimized for inhibitor-rich environmental matrices. Mechanical and chemical lysis were achieved via bead-beating with 0.1 mm glass beads and kit-provided buffers. Following inhibitor removal and column purification, DNA was eluted in 100 µL buffer. Concentration and purity were assessed by NanoDrop spectrophotometry and agarose gel electrophoresis. To contextualize viral dynamics within the broader microbiome, bacterial and microeukaryotic communities were characterized from matched DNA extracts. The V4–V5 region of the 16S rRNA gene was amplified using primers 515 F and 907R³⁰, while the V4 region of the 18S rRNA gene was amplified using primers A-528 F and B-706R³¹. PCR reactions (25 µL) contained 10× FastPfu buffer, dNTPs, 0.4 µM of each primer, bovine serum albumin, FastPfu DNA polymerase (TransGen Biotech, China), and 20 ng template DNA. Thermocycling conditions consisted of 95 °C for 5 min; 25 cycles of 95 °C for 30 s, 55 °C for 30 s, and 72 °C for 90 s; followed by a final extension at 72 °C for 10 min. Amplicons were purified and sequenced on the Illumina HiSeq 4000 platform (2 × 250 bp).

Shotgun metagenomic sequencing of DNA viral communities

For virome characterization, approximately 1 µg of total DNA per selected sample was used for library preparation and shotgun metagenomic sequencing. Libraries were constructed and sequenced on the Illumina NovaSeq 6000 platform (paired-end 2 × 150 bp) at Novogene (Tianjin, China), generating approximately 20 Gb raw data per sample. This approach enabled stage-resolved profiling of DNA viral assemblages without prior viral enrichment, thereby capturing viral DNA signals derived from both extracellular and intracellular sources within the broader microbial community. Sequencing coverage estimates for the analyzed DNA viral communities were close to 100% across all samples, indicating adequate representation of the detectable viral diversity under the applied metagenomic workflow.

Bioinformatic processing

Amplicon sequencing data (16S and 18S) were processed using the LotuS2 pipeline v2.22^32,33. Raw reads were quality-filtered and denoised using DADA2 to infer exact amplicon sequence variants (ASVs)³⁴. Taxonomic assignments were performed using the RDP classifier with an 80% bootstrap confidence threshold against the SILVA v138 database (prokaryotes) and PR2 v4.12.0 database (microeukaryotes)^35,36. Sequencing coverage was assessed for both amplicon datasets to evaluate whether the sequencing depth adequately represented the sampled communities. Sequencing coverage was assessed using Good’s coverage estimator, as implemented in the microeco R package, with coverage ranging from 98.66% to 99.69% (mean, 99.29%) for the 16 S rRNA dataset and from 99.49% to 99.97% (mean, 99.82%) for the 18 S rRNA dataset, indicating high coverage of the detectable bacterial and microeukaryotic communities across the analyzed samples.

Shotgun metagenomic reads were quality-controlled using KneadData v0.6.1 to remove low-quality sequences and potential host contamination³⁷. Clean reads were taxonomically classified using the Microbial Community Profiling (MoP) pipeline^38,39, aligning reads against the NCBI Virus RefSeq database (release 218) under default parameters⁴⁰. The proportion of reads assigned to NCBI Virus RefSeq was calculated for each metagenomic sample. Virus-classified reads represented a small fraction of total reads, ranging from 0.0065% to 1.44% across samples.

Network analysis

To account for the compositional nature of sequencing-derived abundance data, association networks were inferred using the SPIEC-EASI (Sparse InversE Covariance estimation for Ecological Association Inference) framework. Separate networks were constructed for viral, bacterial, microeukaryotic, and integrated cross-domain datasets. Prior to network inference, low-abundance taxa were filtered using a relative abundance threshold of 0.001 (0.1%) to reduce sparsity. Network inference was performed in R using the microeco package with the SPIEC-EASI approach^41,42. Network topology was characterized based on node degree, positive and negative associations, and module assignment, with highly connected taxa identified according to node degree. Network visualization was performed using Gephi v0.9.2⁴³.

Statistical analysis

Differences in viral community composition across treatment stages and seasons were evaluated using permutational multivariate analysis of variance (PERMANOVA) and analysis of similarities (ANOSIM) based on Bray–Curtis dissimilarities^44,45. Permutational analysis of multivariate dispersions (PERMDISP) was additionally performed using Bray–Curtis dissimilarities to assess the homogeneity of multivariate dispersion among groups⁴⁶. Principal coordinates analysis (PCoA) was used to visualize beta-diversity patterns. Permutation-based statistical tests were conducted using 999 permutations, with statistical significance defined at p = 0.05⁴⁷. All statistical analyses were conducted using R version 4.1.0 (https://www.r-project.org/)⁴⁸ and PRIMER v7.0.21 (Quest Research, New Zealand)⁴⁷.

Stage-driven restructuring of the DNA virome between the inlet and maturation pond

PCoA revealed pronounced stage-dependent structuring of DNA viral communities between the inlet and maturation pond samples (Fig. 1). The first ordination axis (PCo1), explaining 59.2% of total variance, separated upstream compartment (i.e., inlet) from the maturation pond samples. PERMANOVA confirmed significant stage effects (R² = 0.53, p = 0.033), and ANOSIM further supported strong compositional separation (R = 0.75, p = 0.027). However, PERMDISP detected significant heterogeneity in multivariate dispersion between the inlet and maturation pond communities (F = 6.185, p = 0.001), indicating that the significant PERMANOVA result may reflect both differences in group location and within-stage dispersion. In contrast, seasonal effects were negligible (PERMANOVA R² = 0.23, p = 0.888; ANOSIM R = − 0.27, p = 0.727), indicating that treatment stage exerted substantially stronger influence on viral community structure than temporal variability. Notably, maturation pond samples formed a discrete cluster distinct from inlet, highlighting its role as a terminal ecological compartment within the system. These patterns indicate progressive virome restructuring along the treatment gradient, culminating in a compositionally differentiated maturation pond assemblage.

Taxonomic reorganization and maturation pond-specific differentiation

Across all stages, the virome was dominated by tailed double-stranded DNA bacteriophages assigned to class Caudoviricetes. However, their proportional contribution shifted along the treatment train, increasing from 73.86% at the inlet to 79.21% in the maturation pond. Human-associated viral families, including Adenoviridae and Poxviridae, were detected at the inlet but were not detected in the maturation pond. Inoviridae declined markedly (2.23% to 0.13%), while Microviridae exhibited moderate reduction (0.27% to 0.16%). In contrast, several viral groups emerged or were enriched in the maturation pond. Iridoviridae (8.1379%) and signals assigned to Herpesvirales (0.0168%) were detected exclusively at the maturation stage, indicating stage-specific ecological differentiation (Fig. 2A). These maturation-associated signals are consistent with potential co-variation with locally established microbial or microeukaryotic communities. Seasonally, Caudoviricetes remained dominant during autumn, spring, and summer (> 84%) but declined markedly in winter (45.86%) (Fig. 2B), reflecting moderate seasonal modulation without overriding stage-driven structuring. Viral alpha-diversity indices differed significantly among WSP stages. Shannon diversity showed a significant stage effect (Kruskal-Wallis; p = 0.015), and similar patterns were observed for Simpson index (p = 0.030), inverse Simpson index (p = 0.030), and Pielou’s evenness (p = 0.048). Diversity and evenness were generally highest in inlet samples and lower in the facultative and maturation ponds, suggesting a downstream shift toward a less diverse and less even detectable viral DNA assemblage. Collectively, these results indicate that the maturation pond was not simply a reduced-diversity endpoint, but a compositionally distinct downstream compartment characterized by selective changes in viral DNA profiles and stage-specific taxonomic signatures.

Taxa- and species-specific relative abundance changes

Relative abundance changes between the inlet and maturation pond were highly taxon dependent (Fig. 3). At the family/class level, Poxviridae and Adenoviridae showed a 100% decrease in relative abundance, indicating that these taxa were detected at the inlet but not detected in the maturation pond under the sequencing and analytical conditions used. Inoviridae showed a substantial decrease in relative abundance (94%), whereas Microviridae showed a more moderate decrease (41.7%). In contrast, the Caudoviricetes signal increased in relative abundance downstream (− 182.9%), indicating relative enrichment rather than attenuation.

At the species level, several phage-associated signals showed a 100% decrease in relative abundance, including Cronobacter phage ENT47670, Lactococcus phages bIL286, bIL309, and ul36, Mycobacterium phage Brocalys, and others. Conversely, increased relative abundance was observed for Pseudomonas phage MD8 (− 331.8%), Burkholderia virus Bcep22 (− 1107.4%), and Planktothrix phage PaV-LD (− 68,848.3%). These negative values indicate downstream enrichment in relative abundance and should not be interpreted as “negative removal,” absolute viral increase, or direct evidence of infective replication. The very large enrichment value for Planktothrix phage PaV-LD likely reflects its very low inlet relative abundance combined with higher detection in the maturation pond. Overall, these results indicate that viral DNA profiles changed selectively across the WSP treatment train, with some taxa decreasing to below detection in the maturation pond and others showing downstream relative enrichment as part of broader community compositional restructuring.

Cross-domain ecological organization across matched viral, bacterial, and microeukaryotic datasets

The integrated SPIEC-EASI-based cross-domain network comprised 329 nodes and 4544 inferred conditional associations (Fig. S1), including 161 bacterial ASVs, 104 microeukaryotic ASVs, and 64 viral taxa. Cross-domain associations represented the majority of the network (3097 edges; 68.2%), with the highest number of links observed between bacteria and microeukaryotes (1218 edges), followed by viruses and bacteria (1105 edges), and viruses and microeukaryotes (774 edges). Positive associations predominated overall (2757 edges; 60.7%), particularly in the virus–bacteria and virus–microeukaryote compartments, whereas bacteria–microeukaryote links showed a more balanced distribution of positive and negative associations. The network resolved into four modules, with the largest modules containing taxa from all three domains, indicating an integrated rather than domain-isolated association structure (Fig. S1). Thus, M1 and M2 were comparatively enriched in bacterial and viral nodes, whereas M3 was dominated by microeukaryotic taxa. Highly connected nodes included bacterial ASVs affiliated with Actinobacteriota, Chloroflexi, and Cyanobacteria and viral taxa affiliated with Iridoviridae and Caudoviricetes. The highest-degree nodes were Actinobacteriota-associated ASV101 (degree = 86) and ASV260 (degree = 85), followed by ASV387 and ASV498 (degree = 78 each). Among viral taxa, Invertebrate iridovirus 25 (degree = 77), Invertebrate iridovirus 22 (degree = 72), and Pseudomonas phage phiPSA1 (degree = 71) (Table 1). Overall, the network suggests broad cross-domain covariation across the WSP treatment gradient.

Intra-domain viral network organization

The virus-only co-occurrence network comprised 69 viral taxa (nodes) connected by 109 significant associations (edges) (Fig. 4). The network consisted of 68 positive associations (62.4%) and 41 negative associations (37.6%), indicating structured patterns of viral co-variation across treatment stages within the WSP system. Taxonomically, the network was dominated by Caudoviricetes, with comparatively smaller contributions from Iridoviridae, unclassified viral groups, and other viral taxa. Community organization resolved into several interconnected modules of varying sizes, reflecting stage-associated restructuring of viral assemblages throughout the treatment process. Positive associations likely reflect shared ecological responses or similar distribution patterns across treatment stages, whereas negative associations may indicate contrasting abundance trajectories or differing ecological preferences among viral taxa (Table 2; Fig. 4). For example, positive associations were detected among Salicola phage CGphi29, Enterobacteria phage If1, Mycobacterium phage Milly, and Enterobacteria phage HK629, suggesting partially overlapping co-abundance patterns among these bacteriophages within the WSP system. In contrast, negative associations involving Enterobacteria phage HK629 were observed with Escherichia phage pro483 and Escherichia virus P2, indicating contrasting abundance trajectories among specific Caudoviricetes-associated phages across treatment stages.

Bacterial community co-occurrence network

The bacterial association network inferred using the SPIEC-EASI compositional data-aware framework comprised 161 ASVs connected by 2805 positive conditional associations (Fig. S2). The network showed a relatively high level of connectivity, with an average degree of 34.85 and a density of 0.218, indicating that bacterial ASVs were broadly interconnected across the WSP system. The short average path length of 1.24 and network diameter of 2 further indicate a compact topology, where most taxa were separated by only one or two network steps (Table 2). This suggests that the bacterial community formed a highly integrated association structure rather than a fragmented network.

The clustering coefficient of 0.253 indicates moderate local cohesiveness, suggesting that some bacterial ASVs formed locally connected neighborhoods. However, the relatively low modularity value of 0.145 shows weak compartmentalization, indicating that the network was not strongly divided into isolated modules. Instead, the bacterial community appeared to be organized as a broadly connected network with limited modular separation. The network centralization (0.245) and heterogeneity (0.390) further indicated variation in connectivity among bacterial ASVs. Node-level analysis identified highly connected ASVs distributed primarily across modules M1–M3 (Table 3). ASV16, affiliated with Actinobacteriota (Actinobacteria; PeM15), showed the highest degree (74) and eigenvector centrality (1.000), whereas ASV15 (Actinobacteriota) exhibited the highest betweenness centrality (345) among the ten most highly connected ASVs.

Zi–Pi analysis further differentiated the structural roles of these highly connected nodes. ASV157, assigned to the genus Propioniciclava, was classified as a module hub (Zi = 2.549, Pi = 0.517). ASV34, ASV238, ASV576 (Planktothrix NIVA-CYA 15), ASV50, and ASV18 (Cyanobium PCC-6307) were classified as connectors, whereas ASV16, ASV15, ASV54, and ASV13 (CL500-29 marine group) were peripheral nodes despite their relatively high degree values (Table 3). Collectively, these results demonstrate that high node connectivity did not necessarily correspond to a hub or connector role and provide a quantitative characterization of bacterial network organization beyond visual inspection of the complete network.

Taxonomically, the highly connected nodes were predominantly affiliated with Actinobacteriota, with additional representation from cyanobacterial taxa, including Planktothrix and Cyanobium (Table 3). All associations involving the ten most highly connected ASVs were positive in the inferred network. The network also contained genera with pathogenic potential (e.g., Klebsiella, Mycobacterium, Acinetobacter, and Streptococcus), although only Klebsiella pneumoniae was resolved to species level; therefore, confirmation of clinical relevance would require species- or strain-resolving approaches.

Microeukaryotic community co-occurrence network

The microeukaryotic association network comprised 104 eukaryotic ASVs connected by 306 inferred associations. Positive associations predominated, accounting for 207 edges (67.6%), whereas negative associations represented 99 edges (32.4%). The node table showed that the network was dominated by Opisthokonta (53 ASVs) and Alveolata (25 ASVs), followed by Archaeplastida (11 ASVs), Stramenopiles (7 ASVs), Rhizaria (3 ASVs), Amoebozoa (1 ASV), and four unclassified eukaryotic ASVs. The network was mainly structured by Opisthokonta-associated links, particularly Opisthokonta–Opisthokonta associations (95 edges) and Alveolata–Opisthokonta associations (78 edges). These were followed by Archaeplastida–Opisthokonta (29 edges) and Opisthokonta–Stramenopiles associations (28 edges), indicating that fungal, rotifer-related, ciliate, algal, and stramenopile taxa were tightly embedded within the same association structure (Fig. 5).

Several highly connected nodes belonged to Opisthokonta, including Cryptomycota, Filinia longiseta, and other rotifer- or fungal-related ASVs. The strong representation of rotifers, fungi, ciliates, algae, and stramenopiles suggests that the WSP microeukaryotic network may reflect multiple ecological processes, including organic-matter decomposition, and shared responses to pond conditions. Notably, Amoebozoa was represented by Arcella vulgaris, which formed four associations: three positive and one negative. Although Amoebozoa was not numerically dominant, its connection to Opisthokonta, Alveolata, and unclassified eukaryotic taxa suggests that testate amoebae may participate in the broader microeukaryotic association structure of the WSP system (Fig. 5).

Stage-driven restructuring of the DNA virome

Our results demonstrate clear stage-driven restructuring of the DNA virome across sequential WSP compartments. Viral community composition was strongly dictated by treatment stage rather than seasonal variation, consistent with previous observations from the same WSP system^4,9. Although the inlet and facultative-pond metagenomic datasets overlap with those analyzed by Rizk et al. (2026)⁹, the present study provides an independent analytical contribution by incorporating the maturation pond and examining virus–virus and virus–bacteria–microeukaryote association networks, analyses that were not undertaken in the previous study. The influent was characterized by enteric DNA viruses and bacteriophages associated with human fecal waste, whereas downstream facultative and maturation ponds exhibited turnover toward viruses linked to the developing pond ecosystem. The dominance of tailed bacteriophages within Caudoviricetes in early-stage ponds aligns with global reports of phage-rich raw sewage viromes^49,50,51. As wastewater progressed through the system, viruses associated with human gut bacteria and enteric hosts declined, while viral groups more typical of environmental aquatic systems became proportionally more prominent. By the maturation pond, the virome included a greater representation of viruses potentially associated with algae, protozoa, and environmental bacteria^52,53,54, reflecting community compositional transition from fecal-dominated influent to a stabilized, phototrophic pond environment.

Such successional patterns are consistent with the shifting microbiological milieu across treatment stages. The inlet is characterized by high organic loading and fecal microbial inputs, whereas maturation ponds resemble enriched aquatic environments with active phytoplankton and zooplankton communities. Similar stage-resolved viral differentiation between influent, secondary treatment, and final effluent has been reported in other wastewater systems^49,50,51. Notably, the stage effect in our system was more pronounced than seasonal variability, reinforcing the dominant role of ecological filtering imposed by sequential anaerobic, facultative, and maturation compartments⁴. The prolonged retention and relative hydraulic stability of WSPs may buffer seasonal fluctuations, leading to reproducible stage-specific viral profiles.

Taxa-specific removal patterns of viruses

Viral attenuation within the WSP system was clearly family dependent. Large dsDNA families such as Poxviridae and Adenoviridae, detected primarily in the inlet, were not detected in the maturation pond. Reported adenovirus removal in maturation ponds often ranges around 1–2 log units^4,55, and differences between studies may partly reflect methodological variation, including targeted qPCR versus shotgun metagenomics approaches.

Virus reduction in algal-based systems is attributed to multiple interacting mechanisms, including solar irradiation, reactive oxygen species production, elevated dissolved oxygen, high pH, sedimentation, and adsorption to biomass^56,57. In shallow maturation ponds, these stressors operate concurrently. Attachment of viruses to settling particulates or biomass can facilitate removal via sedimentation^58,59, while wastewater-derived natural organic matter may enhance reactive oxygen species formation relative to natural aquatic systems^59,60. Collectively, these mechanisms are likely to contribute to the selective attenuation observed for certain viral taxa.

Conversely, large eukaryote-associated viral families, including Mimiviridae and Iridoviridae, were detected in the maturation pond but not at the inlet under the sequencing and analytical conditions used. This downstream occurrence is best interpreted as a stage-associated relative enrichment of viral DNA signals, potentially linked to local community restructuring, environmental filtering, particle-associated DNA retention, or differential detectability. Mimiviridae are known to infect amoebae and have been isolated from natural pond environments⁶¹, while Iridoviridae infect a range of aquatic organisms including invertebrates and vertebrates⁵⁴. The maturation pond, characterized by increased protozoan and planktonic activity, may provide ecological niches supporting persistence of detectable viral DNA signals. Additionally, viruses may accumulate within aggregates, biofilms, or settling particles in later ponds^54,59, facilitating detectable persistence without necessarily indicating active viral proliferation.

These observations underscore that viral removal in WSPs is not uniform across taxa. Viral fate is influenced by intrinsic properties such as particle size, genome type, envelope structure, and host dependence²¹. Consequently, reliance on a single viral surrogate may not adequately represent system-wide performance. Prior recommendations advocating multi-target viral monitoring in wastewater treatment are therefore well justified^21,59,62.

Cross-domain interactions in the WSP system

The cross-domain network revealed that viruses, bacteria, and microeukaryotes were embedded within a shared conditional association structure across the WSP treatment gradient. Positive associations predominated, suggesting that many taxa exhibited similar abundance trajectories across treatment stages. However, given the strong physicochemical gradients in WSP systems, including changes in pH, dissolved oxygen, organic loading, solar exposure, and UV irradiation, these associations should be interpreted cautiously. Positive edges may reflect shared environmental filtering, stage-associated community restructuring, host-community context, or particle-associated DNA dynamics, whereas negative edges may indicate contrasting abundance trajectories or differential responses to pond conditions^9,23.

The positive association between Gordonia phage Bantam and the ciliate taxa Tetrahymena spp. and Vorticella may reflect shared abundance trends or shared environmental filtering, rather than direct trophic interactions. Although ciliates can ingest bacteria and can internalize viral particles, the present network cannot determine whether these processes underlie the observed association⁶³. Conversely, the negative association between Pseudomonas virus F116 and Rotifera may reflect contrasting abundance trajectories or differing habitat preferences. Rotifers and other filter feeders are capable of ingesting viral particles and virus-infected bacteria^64,65. Alternatively, host ecological preferences may drive inverse abundance trends²³.

Intra-domain ecological organization across microbial domains

Comparison of the intra-domain SPIEC-EASI networks revealed domain-specific patterns of microbial organization across the WSP treatment gradient. The bacterial network showed a more integrated association structure, suggesting broad co-variation among bacterial ASVs in response to stage-specific pond conditions. In contrast, the viral network appeared more compartmentalized, indicating that viral DNA signals may be structured by differences in host-community context, environmental persistence, or treatment-stage distribution. The microeukaryotic network showed an intermediate pattern, suggesting partial compartmentalization among major eukaryotic groups and uneven connectivity among taxa. Overall, these differences indicate that viruses, bacteria, and microeukaryotes respond differently to the WSP treatment gradient. The intra-domain networks of viruses, bacteria, and microeukaryotes consistently revealed structured and modular organization across the WSP system, with positive associations predominating in all domains^25,26. This pattern suggests that treatment-stage gradients act as strong environmental filters synchronizing microbial abundance trajectories rather than generating widespread antagonistic restructuring. The viral network was dominated by Caudoviricetes, consistent with the well-documented prevalence of dsDNA bacteriophages in wastewater systems^9,49,50,51. The strong modular structure and high connectivity among phage taxa indicate coordinated host-linked dynamics¹⁷. Limited negative virus–virus associations suggest niche differentiation among host groups, such as segregation between cyanophage-associated signals and phages infecting heterotrophic bacteria, a pattern consistent with ecological partitioning in aquatic systems^66,67.

The bacterial network was similarly modular, dominated by Proteobacteria and Actinobacteriota, which are common and responsive phyla in wastewater environments^4,67. The prevalence of positive associations reflects shared environmental preferences and potential metabolic complementarity, whereas structured cross-phyla negative links indicate competitive or divergent ecological optima under changing physicochemical conditions^9,31.

The microeukaryotic network exhibited overwhelming positive connectivity, particularly among Opisthokonta and Alveolata. In pond ecosystems, protistan communities respond strongly to light, nutrient gradients, and bacterial prey availability⁴. Their cohesive structure suggests synchronized responses to treatment-stage conditions rather than pervasive exclusion. Although protistan grazers can influence bacterial populations and may indirectly affect viral dynamics^22,23, the present SPIEC-EASI network cannot determine whether such trophic processes underlie the observed patterns. Accordingly, these associations should be interpreted as compositionality-aware conditional associations, reflecting potential co-variation, shared environmental filtering, or contrasting abundance trajectories across WSP stages. Collectively, the intra-domain networks suggest stage-associated microbial restructuring across the treatment gradient, but they should not be interpreted as direct evidence of trophic coupling, predation, or competition.

Limitations and future directions

A key limitation of this study is that shotgun metagenomic sequencing was performed on total environmental DNA without prior enrichment of virus-like particles. Consequently, detected viral DNA may originate from free virions, intracellular viral genomes, integrated prophages, cell- or particle-associated viruses, or extracellular and relic DNA adsorbed to suspended solids. Changes in viral relative abundance should therefore be interpreted as shifts in DNA viral community composition rather than evidence of active viral replication, infectivity, or proliferation. In particular, downstream enrichment of some bacteriophage signals may reflect changes in host abundance, particle-associated DNA retention, differential recovery, or low initial abundance rather than confirmed lytic activity. Although sequencing coverage was close to 100% across the analyzed samples, total-community metagenomics may also have limited sensitivity for rare viral taxa; thus, non-detection should not be interpreted as evidence of absolute absence.

The study was additionally restricted to DNA viruses and did not characterize RNA viruses, potentially underrepresenting overall viral diversity and associated virus–microbe relationships within the WSP system. Furthermore, the absence of field and laboratory negative controls prevents complete exclusion of background or reagent-associated contamination, particularly for low-abundance viral signals, which should therefore be interpreted cautiously.

Future studies should integrate total metagenomics with virus-like particle enrichment, DNase pretreatment, prophage-resolved genome reconstruction, virus–host linkage approaches, metatranscriptomics, and targeted qPCR or ddPCR, together with appropriate field and laboratory negative controls. Incorporating both DNA and RNA viromics with time-series sampling and physicochemical measurements would enable more comprehensive characterization of viral diversity and help distinguish free virions, active infections, prophages, and residual viral DNA, thereby providing a more mechanistic understanding of viral fate and virus–microbe relationships in WSP systems.

This study demonstrates that DNA viral communities in a full-scale WSP system are primarily structured by treatment stage rather than season. Viral composition exhibited clear stage-dependent turnover along the treatment train, with the maturation pond forming a distinct assemblage relative to upstream compartments. Viral relative abundance changes were strongly taxon dependent. Human-associated viral families, including Adenoviridae and Poxviridae, were detected in the inlet but were not detected in the maturation pond under the sequencing and analytical conditions applied, whereas other viral groups, such as Microviridae, showed partial persistence. In contrast, several bacteriophage-associated DNA signals increased in relative abundance downstream, indicating that viral community shifts reflect selective reduction and stage-associated compositional restructuring rather than uniform removal or confirmed viral proliferation. Network analyses revealed predominantly positive intra- and inter-domain associations among viruses, bacteria, and microeukaryotes, consistent with coordinated responses to shared environmental gradients across treatment stages. The maturation pond was characterized by structured and modular viral organization, supporting its role as the terminal ecological compartment of the system. Overall, these findings indicate that viral attenuation in WSPs is selective and embedded within broader microbiome reorganization. Stage-resolved virome analysis may therefore contribute to a more comprehensive understanding of microbial dynamics in nature-based wastewater treatment systems.