High density sampling reveals the spatiotemporal characteristics of microbial communities in a full-scale municipal wastewater treatment plant

Highlights

	● High density spatiotemporal sampling was adopted to investigate a municipal WWTP.
	● Spatially independent corridors showed high microbial community similarities.
	● Three distinct stable states of the microbial community were observed over a year.
	● Conserved function over microbial community succession was observed in the WWTP.

Abstract

Insights into the microbial communities in municipal wastewater treatment plants (WWTPs) are critical for the optimization of biological nutrient removal process. However, our understanding about the spatiotemporal characteristics of the microbial communities in WWTPs remains limited. In the present study, 264 samples were collected biweekly from four spatially independent corridors in a typical municipal WWTP. The annual compositional and metagenomic characteristics were investigated based on multiple ecological indicators using statistical tests. The results revealed that the microbial community compositions from the four corridors showed significantly high similarities, as revealed by the statistical analysis at the operational taxonomic unit (OTU) level. Consistent with the OTU level results, the functionality of the microbial communities in the four independent corridors also showed significant similarity. In comparison, the dynamics of the microbial community over the year showed two successional peaks of the microbial communities with the spatial similarity, and this resulted in three alternative stable states of the microbial communities in a calendar year. The microbial communities only drifted in July and November, suggesting an uneven community succession pattern driven by seasonal variation in environmental conditions. The functional characteristics were found to be relatively conservative compared to the microbial community succession, which revealed the decoupling between the composition and functionality of the microbial community in the municipal WWTP. The present study provides an in-depth overview of the microbial communities in a municipal WWTP and will be useful for the establishment of the connection between ecological characteristics and the operational stability of WWTPs.

Graphical abstract

Keywords

Activated sludge / Microbial community / Spatial similarity / Successional pattern / Metagenomics

1 Introduction

In typical municipal wastewater treatment plants (WWTPs), carbon and nutrient removal and resource recovery are primarily achieved by complex microbial communities (Xia et al., 2018; Woo et al., 2022). Currently, efficient pollutant removal of WWTPs is primarily mainly achieved by optimizing the microbial community based on operational parameter management. Nevertheless, an understanding of the microbial communities in WWTPs have been overlooked. From the perspective of microbial ecology, WWTP is a semi-artificial ecosystem that accommodates thousands of bacterial genera with high biomass concentrations ranging from 2 to 10 g/L that mediate specific engineered functions (Zhang et al., 2012; Tong et al., 2019; Chen et al., 2023). Recent studies have suggested that the operational stability and performance of WWTPs are closely related with the characteristics of the microbial communities (Cao et al., 2020; Wang et al., 2021; Dottorini et al., 2023). These findings highlight the importance of an in-depth understanding of the microbial communities in WWTPs.

Ecological surveys in WWTPs at spatial and temporal scales have deepened our understanding of the microbial communities within WWTPs (Woo et al., 2022). Numerous studies have verified the existence of shared bacteria (so called ‘core bacteria’) in different WWTPs. Most of these members are involved in functions that are highly related with the WWTP operational efficiency (Zhang et al., 2023). In addition, the distinct spatial heterogeneity of microbial communities in WWTPs has been revealed by global crossover research, although these WWTPs had similar engineered functions (Zhang et al., 2012; Wu et al., 2019). Similar to the spatial differences, the distinct diversity of the microbial communities over time has also been widely reported. Seasonal succession and incomplete annual regression of microbial communities have been verified based on multiple temporal sampling intervals ranging from daily to trimonthly (Jiang et al., 2018b; Wang et al., 2021).

Recent studies have revealed significant differences in the microbial community in WWTPs with similar influent, processes and located in similar regions (Wu et al., 2019; Yang et al., 2020). WWTPs normally include more than one independent corridor to reduce the loading rate of a single reactor. Physicochemical parameters are commonly used to maintain consistent operation and performance of these independent corridors (Cui et al., 2023). However, the ecological features of these corridors during long-term operation have received less attention. Whether the microbial communities in these corridors maintain identical composition during perennial independent operation over time remains unclear. The relationship among the microbial communities fundamentally determines whether these spatially separated corridors should be independently operated as different systems. Beyond the spatial features, the characteristics of the microbial community under the stable operation period of WWTPs remains unclear. Considering that the most recent studies of WWTPs primarily focused on either the long-term sampling of a single WWTP or a single sampling site for each WWTP with different types of processes (Cao et al., 2020; Chen et al., 2022), a high-resolution investigation at both the spatial and temporal scales will contribute to a more comprehensive understanding of microbial communities in municipal WWTPs.

In this study, the microbial community from four independent corridors of a typical municipal WWTP are characterized biweekly as the biotic feature of the WWTP in a calendar year. Since previous studies have revealed distinct microbial communities across WWTPs over space and time, the primary goal of the present study is to characterize the spatiotemporal properties of the microbial communities from four independent corridors in a typical municipal WWTP during long-term operation. To address this goal, the composition, structure, and successional dynamics of the microbial communities in the WWTP are analyzed. To elucidate the spatiotemporal functionality of the microbial communities, representative samples are selected for metagenomics sequencing. The implication of microbial ecology for guiding WWTP operational strategies is also discussed. The results of the present study provide further understanding of the potential ecological parameters for optimizing biochemical processes in WWTPs.

2 Materials and methods

2.1 Characteristics of the high-capacity WWTP and activated sludge sampling

2.1.1 Overview of the WWTP

The activated sludge was collected biweekly from four physically separated corridors in a full-scale WWTP from January 2022 to January 2023 (Fig.1). The primary portion of the influent of the WWTP was domestic wastewater, while industrial and chemical wastewater only accounted for a small proportion. The influent contained 46.4 mg/L of NH₄⁺−N, 4.6 mg/L of PO₄³⁻−P and 305 mg/L of chemical oxygen demand (COD). The mean influent quantity of the WWTP was 1000000 m³/d. The preliminary treatment units included grilles, lift pump houses, aerated grit chambers, and primary sedimentation tanks. The anaerobic-anoxic-oxic (AAO) process was adopted for simultaneous biological nitrogen and phosphorous removal. The hydraulic retention time (HRT) of the process was maintained at 8–10 h based on the pollutant removal efficiency of the WWTP. The sludge retention time (SRT) was controlled at 12–15 d based on the suspended solids concentration of the biochemical tank. And the average concentrations of NH₄⁺−N, PO₄³⁻−P and COD in the effluent were 2.4, 0.6, and 45 mg/L, respectively (Guo et al., 2017). The WWTP included four independent corridors. Each corridor had the same influent quality, volume, operational conditions. In addition, the four corridors had independent sedimentation tanks and sludge return systems. Overall, the four corridors were identical, except for the physical separation.

Fig.1 Schematic diagram of the spatiotemporal sampling strategy in the wastewater treatment plant (WWTP). The activated sludge from the four corridors was sampled from January 2022 to January 2023 in triplicate, which resulted in 264 samples. All samples were sequenced based on the V3–V4 region of 16S rRNA gene. The samples colored with dark gray were selected for the metagenomic sequencing.

Full size|PPT slide

2.1.2 Sampling strategy and procedures

The activated sludge of the four corridors were sampled in triplicate at the end of the biochemical tanks at one meter below the surface. The retrieved activated sludge samples were concentrated using gravity sedimentation to remove the excess moisture. The sampling procedure at 22 time points from the four corridors in triplicate resulted in a total of 264 samples being sequenced. Eight representative sampling points were selected for the metagenomics analysis to investigate the functional characteristics of the microbial communities during different seasons and in different corridors. Before DNA was extracted, the activated sludge samples were freeze-dried by a freeze dryer (Free Zone, Labconco, USA).

2.2 DNA extraction, amplicon and metagenomic gene sequencing

DNA was extracted from the freeze-dried activated sludge samples by FastDNA® SPIN Kit for Soil (MP Biomedicals, USA) following the standard protocol. The purity and contamination of the extracted DNA were checked by a spectrophotometer (ND-1000, NanoDrop® Technologies, USA). DNA libraries was constructed prior to high-throughput sequencing by bridged polymerase chain reaction (PCR) amplification. The bacterial 16S rRNA gene was captured by the universal primers, 338F (5′-ACTCC TRCGG GAGGC AGCAG-3′) and 806R (5′-GGACT ACCAG GGTAT CTAAT-3′). The PCR reaction system included 0.4 × 10⁻⁶ mol/L of paired primers, 1 × 10⁻⁶ mol/L of four kinds of deoxy-ribonucleoside triphosphate, 1 × 10⁻⁶ mol/L of FastPfu polymerase, 4 × 10⁻⁶ L of 5 times PCR buffer, and 10 ng of template DNA. The bridged PCR reaction was conducted for 30 cycles with an initial denaturation step and a final extension step. Paired-end sequencing of the amplicon library were performed on a Miseq PE 300 platform (250 base pair for both sides, Illumina, USA). For the metagenomic sequencing, the extracted DNA from activated sludge samples was fragmented to an average insert size of 400 bp by Covaris M220 (Gene Company Limited, China) for the subsequent library construction. Then the paired-end library was constructed using NEXTFLEX Rapid DNA-Seq (Bioo Scientific, USA). The high throughput sequencing targeted on the environmental samples was performed on the Illumina Novaseq 6000 (Illumina, USA) following the standard protocol with NovaSeq Reagent Kits. The raw metagenomic and 16S rRNA gene sequencing data sets were uploaded to the Sequence Read Archive Database of the National Center for Biotechnology Information under the bioproject number of PRJNA1029465.

2.3 Amplicon and metagenomic sequences processing

The raw amplicon sequencing data generated by Miseq PE 300 platform were pre-treated using Usearch (version 11) in Linux system (Edgar, 2013). In short, the barcode and primers of each amplicon was cut after sequencing. The clean amplicons with low sequencing qualities (error rates > 1%) were filtered out. The paired-end reads with high quality were merged to obtain the full-length amplicons. The reads with low abundance were filtered out after the quality control. The high-quality reads were clustered into unique operational taxonomic units (OTUs) using Usearch. Chimeric sequences were identified and deleted from the amplicon data set according to the Ribosomal Database Project (RDP) reference database (Cole et al., 2014). Each qualified OTU was further aligned to the SILVA rRNA database to obtain reliable taxonomies (Quast et al., 2013). The data processing resulted in a total of 5891409 high quality reads from 264 samples. The retrieved reads were clustered into 821 OTUs. All samples were rarefied to 10000 reads (based on the lowest read depth in all samples) prior to further analysis. The rarefaction curves of the richness index for each sample were checked to ensure that the sequencing depth of the high-throughput sequencing was technically adequate.

For the metagenomic sequencing data, the adaptors of paired-end reads were cut after sequencing. Fastp (version 0.20.0) was used to identify and remove the low-quality reads or those with lengths less than 50 base pair (Chen et al., 2018). MEGAHIT (version 1.1.2) was used to assemble the high-quality reads into contigs after quality control (Li et al., 2016). The contigs with lengths less than 300 base pair were deleted from the assemblage data set. Prodigal was used to predict the open reading frames (ORF) based on each assembled contig (Hyatt et al., 2010). Then the Cluster Database at High Identity with Tolerance (CD-HIT, version 4.6.1) was used to construct a non-redundant gene catalog with a sequence identity and coverage rate of 90% (Fu et al., 2012). The abundance of genes retrieved from the high-quality contigs were counted with 95% similarity threshold using the Short Oligonucleotide Analysis Package (SOAPaligner, version 2.21) (Li et al., 2008). The gene annotation procedure was conducted by Diamond (version 0.8.35) (Buchfink et al., 2021) against the Kyoto Encyclopedia of Genes and Genomes (KEGG) database with an E-value threshold of 10⁻⁵ (Kanehisa and Goto, 2000).

2.4 Bioinformatic analysis

To minimize the bias of a single estimator on the microbial community diversity, multiple alpha diversity estimators including the OTU richness, ACE, Chao1, and Shannon indices, were calculated using the R package, “Vegan”. The Spearman correlation coefficients of the alpha diversity indices were calculated to investigate the spatiotemporal characteristics of the microbial communities from different activated sludge samples. A principal coordinate analysis (PCoA) was visualized using the “ggplot2” and “vegan” packages in R software. The Spearman correlation coefficients of the bray-Curtis dissimilarity among different corridors were compared using a Mantel test. The retrieved results of Mantel test was characterized as network diagram on Gephi 0.9.2 (Bastian and Jacomy, 2009). The Spearman correlation coefficients of the bacterial species with relative abundance greater than 0.1% were also calculated to characterize the co-occurrence patterns of the bacteria at different time and/or space.

2.5 Mathematical statistics

Shapiro-Wilk test was performed to examine whether the ecological parameters follows a normal distribution. The homoscedasticity of the data sets was checked using either Bartlett’s test or the Fligner-Killeen test. The selection of test method depends on whether the data followed a normal distribution. Depending on the normality of the data set, pairwise significant differences in the retrieved values were analyzed using either a Student’s t-test or a Wilcoxon rank sum test. The significance of the differences in the variance of the alpha diversity indices were examined by either analysis of variance (ANOVA) or the Kruskal–Wallis rank sum test, depending on the normality and homoscedasticity of the data set. Post hoc comparisons were conducted using either Tukey’s honest significant differences test or pairwise Wilcoxon rank sum test.

3 Results

3.1 The spatial similarity of the microbial communities in the four independent corridors

A total of 821 bacterial OTUs were detected in all corridors with the absence of OTU763 in corridor 3 that revealed high similarity of the microbial communities among the different corridors at the OTU level (Fig.2). The microbial community richness indices among the four corridors were distributed over a large range of 355–681 during annual operation (Fig.2). There was no significant difference of the richness indices among the four corridors, as did other alpha diversity indices (Fig.2; Wilcoxon rank sum test, p < 0.1; the other indices are shown in Fig. S1). The pairwise Pearson correlation coefficients of the OTU abundance ranged from 0.661 to 0.848 (Mantel test, p < 0.01), indicating significantly high similarities of the microbial community compositions of the four corridors at the OTU level (Fig.2). The richness index of the microbial communities in the different corridors showed high correlations over the year (0.557–0.816, p < 0.01; Fig.2). The correlation coefficients of the ACE and Chao1 indices between the different corridors ranged from 0.432 to 0.701 (p < 0.01). In addition, the coefficients of the composited indices also ranged from 0.681 to 0.800 (p < 0.01). In addition, the Shannon index showed high correlations of 0.370–0.520 between the four individual corridors (p < 0.01, Fig.2), indicating high alpha diversity similarities among the four independent corridors. The pairwise Pearson correlation coefficients of the Bray-Curtis dissimilarity metrics ranged from 0.634 to 0.789 (Fig.2, Bray Curtis dissimilarity with Mantel test, p < 0.01), indicating significantly high similarities of the microbial communities among the four corridors. These results suggested that the composition and structures of the microbial communities in the four individual corridors showed high similarity over time.

Fig.2 Spatial similarities between the microbial communities in the different corridors of the WWTP. A: the global distribution of retrieved OTUs in the different corridors; B: Annual distribution of Richness index of the microbial communities from each corridor. Wilcoxon rank sum test was adopted to verify the significance of pairwise Richness index differences among different corridors; C: The pairwise correlations of OTU distribution in each corridor based on Spearman’s rank correlation coefficients. OTU table was rarefied to 10000 per sample. OTUs with rarefied abundance of higher than 10 in any sample were kept in the data sets. The abundances of OTUs were log-transformed; D: The pairwise correlations of richness (ACE, Chao1 and Richness indices) between the microbial communities in different corridors base on Spearman’s rank correlation coefficients; E: The pairwise correlations of evenness (Shannon index) between the microbial communities in different corridors base on Spearman's rank correlation coefficients; F: The correlations of successional pattern of the microbial communities in different corridors. The distance matrixes of the microbial communities from each corridor were calculated based on Bray-Curtis dissimilarity. Mantel test was adopted to show the Pairwise correlations of the distance matrixes between each corridor.

Full size|PPT slide

3.2 The successional pattern of the microbial communities in a typical year

A microbial community dynamics analysis was performed to characterize the successional pattern of the four independent corridors. The average richness index of the microbial communities was stable at 589.4 during the first 29 d and then decreased from 577.3 to 551.8 on day 47 (GBD04, Student’s t-test, p < 0.05; Fig.3). The value then gradually increased from 558.5 to 639.4 in 106 d (GBD06–GBD10) and remained stable for 129 d (GBD10–GBD17; Student’s t-test, p > 0.1). At the end of the year, the richness index dramatically decreased from 635.3 to 523.5 in 11 days and remained stable for another 45 d (GBD19–GBD22; Student’s t-test, p > 0.1). Among the 369 d of the sampling period, the richness index of the microbial communities significantly changed over 135 d (36.6%, related to six sampling points), while remaining relatively stable for 234 d (63.4%, related to 16 sampling time points). Other estimators showed similar trends as the richness index (Fig. S2). The dynamics of the alpha diversity of the microbial communities in each corridor showed high consistency with the others, while the outliers only had the minor portion, and this may have been caused by the bias during the sampling and data processing procedures (Fig.3 and S3).

Fig.3 The temporal variation of the alpha diversity (Richness index) of the microbial communities from different corridors. A: The annual distribution of the alpha diversity (Richness index) of the microbial communities from different corridors with pairwise Student’s t-test. The two highest and lowest values of each sampling time point were removed; B: The annual variation of alpha diversity (Richness index) in the four independent corridors. Each curve represents one biological replicate. The diversity curve is smoothed by interpolation.

Full size|PPT slide

The Bray-Curtis dissimilarities were statistically compared to further investigate the successional dynamics of the microbial community over the year. According to the hierarchical clustering (Bray-Curtis dissimilarity, average linkage method, Fig. S4), the microbial communities were divided into three groups by sampling time, except for a few outliers. A clear temporal succession of the bacterial communities was reflected in the PCoA (Fig.4). Two distinct successions of the microbial communities occurred at the transitions of spring and summer (GBD07 – GBD17) and fall and winter (GBD17–GBD22). The two microbial community successions divided the community into three alternative stable states that showed significantly different alpha diversities (Fig.3, Student’s t-test for richness index, p < 0.05) and t distributions in the first two principal coordinates of the PCoA (Fig.4, Bray-Curtis dissimilarity, p < 0.05). Despite the seasonal succession over the year, the microbial communities from the different corridors significantly clustered together at each sampling time point and showed no significant differences. The microbial communities from the different corridors showed similar annual successional patterns (t distribution, p < 0.05).

Fig.4 Principal coordinates analysis (PCoA) of the annual samples from the four different corridors at OTU level. The density curve and boxplot show that the distributions of the three groups were significantly different. The two axes represent a total of 48.3% interpretation. The confidence ellipse based on the t distribution has a confidence level of 0.05.

Full size|PPT slide

Proteobacteria, Actinobacteria, Chloroflexi, Bacteroidetes, Nitrospirae, Firmicutes, Acidobacteria, and Saccharibacteria were the most abundant eight phyla in the WWTP (91.8%–100% in the different samples). All the phylum showed significant variations in their relative abundance over the year (Fig.5). Proteobacteria, Actinobacteria, Chloroflexi and Bacteroidetes are widely reported to constitute a large proportion of the microbial communities in activated sludge systems (Wu et al., 2019; Yang et al., 2020). Although the relative abundances of Proteobacteria and Bacteroidetes significantly fluctuated over the year (13.1%–58.3% and 1.92%–18%, respectively), these two phyla only had slightly significant difference between corridors 2 and 3 (Wilcoxon rank sum test, p > 0.1, Fig.5). Actinobacteria was the second most abundant phylum in the communities, and the relative abundance ranged from 11.7% to 36.6% in the four corridors during the year. The distributions of Actinobacteria in corridors 1 and 2 were significantly different from corridors 3 and 4 (Wilcoxon rank sum test, p > 0.05, Fig.5). Since little attention has been paid to Actinobacteria in the activated sludge systems, some of its members are recognized as filamentous or bulking bacteria and tend to dominate those reactors with an elevated sludge volume index (Knoop and Kunst, 1998; Anderson et al., 2012). However, it is worth noting that no sludge bulking was observed during the experimental period. Chloroflexi ranged from 3.28% to 24.6% in the microbial communities and remained relatively stable in the four corridors. Some members of Chloroflexi were involved in the degradation of complex organics such as polysaccharides and proteins (Kindaichi et al., 2012; Luo et al., 2022).

Fig.5 Taxonomic similarities between the microbial communities in different corridors. A: Overview of the annual dynamics of the top seven phyla in the WWTP. The x-axis represents the 22 sampling time points during the experimental period; B: Spatial similarities of top eight phyla among the four corridors based on Wilcoxon rank sum test. Post hoc comparisons were also conducted by pairwise Wilcoxon rank sum tests.

Full size|PPT slide

3.3 The functional characteristics of the microbial communities over space and time

The Shannon index of the three levels of the KEGG metabolic pathways and Orthology (KO) of each corridor were calculated. The results revealed that the functional composition of the microbial communities among the different corridors showed significant similarities, except there was a mild difference of KO between corridor 2 and corridor 3 (Fig.6 and Fig.6B). Metabolism was used as an example, and the energy metabolism, global and overview maps, and amino acid metabolism were the top three secondary pathways. The relative abundances of the three pathways accounted for 91.1%–91.7% of metabolism and 74.5%–74.8% of the total genes. The nitrogen metabolism pathway, which is highly related to the biological nitrogen removal process, was the primary energy metabolism pathway and showed little difference among the four corridors (Fig.6). The PCoA based on three KEGG metabolic pathways and the gene level further confirmed the spatial stability of the functional composition in the different corridors. The functional compositions of the community in the four corridors (green triangle) clustered into one group and were significantly different from other time points (Fig.7).

Fig.6 The alpha diversity indices of the microbial communities in the four corridors over the year. A and C: The pairwise difference of the Shannon indices among different corridors or different sampling time (corridor 2) at multiple KEGG metabolic pathway analyzed by t-test; B and D: The relative abundance of the KEGG metabolic pathways at different samples.

Full size|PPT slide

Fig.7 Principal coordinates analysis (PCoA) of the annual samples from corridor 2 at different functional levels. In particular, all corridors were sampled at batch GBD11, which resulted in 12 green triangle points. A: The functional characteristics of the microbial communities at KEGG level 2; B: The spatial characteristics of the microbial communities at KEGG level 3; C: The spatial characteristics of the microbial communities at KO level; D: The spatial characteristics of the microbial communities at gene level.

Full size|PPT slide

Compared with the seasonal succession of the microbial community, the temporal variation in functions was relatively conservative. After a slight decrease, the Shannon index of the metabolic pathways recovered to the initial states (GBD1) at GBD17 but dramatically decreased to the minimum at GBD22 and did not increase to the initial states within the year (Fig.6 and Fig.6D). In addition, the Shannon index of the KO continuously rose throughout the year and was also significantly different from the previous year. The functional composition of the microbial communities was stable in the summer (GBD07–GBD17) but showed significant variation between fall and winter (GBD17–GBD22), which was similar as the variation in the microbial communities.

4 Discussion

4.1 Temporal diversity of microbial communities over the spatial stability in the WWTP

Typical WWTPs normally include more than one relative independent corridor with similar influent and capacities. A global-scale research of WWTPs reported that the microbial communities manifested large spatial differences (Wu et al., 2019). However, the spatial similarities of the independent corridors in a single WWTP that is supposed to be an important piece of information for WWTP operation optimization, has been rarely reported. This study provided multiple pieces of evidence for the spatial similarity of the microbial communities in the different corridors with biweekly-level intensive longitudinal sampling. The ecological analysis indicated that the microbial communities remained spatially stable over seasonal dynamics (Fig.1 and Fig.4). The combined ecological characteristics including the OTU composition, alpha diversity index, and Bray-Curtis dissimilarity statistically proved the similarity of the microbial communities in the different corridors (Fig.1). Therefore, the four independent corridors of the WWTP could be treated as the same engineering ecosystem, and operated under similar parameters. In addition, the microbial communities in the four corridors showed similar successional patterns (Fig.3 and Fig.4). Since the current microbial ecological theories have concluded that microbial community assembly is driven by deterministic and/or stochastic processes (Jiang et al., 2018a; Chen et al., 2022), the result of the present study provided strong evidence for microbial community assemblages during long-term operation driven by deterministic processes. The operational parameters and seasonal condition played more important roles than stochastic operational disturbances in the microbial community assemblages.

The present study also provided an ecological understanding of the biological wastewater treatment system, which is a typical complex semi-artificial ecological system. In previous studies, constant operational parameters and performance with stable composition and function of the microbial community was expected in a WWTP (Zhu and Anderson, 2017). Interestingly, this study elucidated two peaks of community succession in a natural year that resulted in three different alternative stable states of the microbial community with a legal effluent of the WWTP. The microbial communities significantly drifted at GBD9 (July 30th) and GBD18 (November 21st), while they remained relatively stable for most of the year (Fig.3 and Fig.4). The alternative stable states of the microbial communities over seasonal factors has been widely observed in many WWTPs (Ju and Zhang, 2015; Wang et al., 2021). In WWTPs, the microbial communities are affected by complex uncontrollable factors and controllable parameters (Yu et al., 2023). Since the operational parameters are relatively stable during annual operation, some uncontrollable seasonal conditions might significantly impact the microbial community (Yu et al., 2023). The seasonal variation in the environmental condition under natural environmental conditions might have been responsible for microbial community succession among the three alternative stable states. Although it is still difficult to conclude the precise succession pattern of the microbial community, the results obtained in this study provided an overview that the potential microbial community dynamics under the stable operation of municipal WWTPs have been overlooked. Due to the inhomogeneous succession of the microbial community over time, sampling intervals should be determined in combination with the variation of the environmental conditions and operating parameters when developing an ecological survey in WWTPs. According to the results of our study, the biweekly sampling strategy captured most of the features of the microbial communities. However, weekly, even daily sampling are recommended during the intense successional period of microbial communities (e.g., GBD9 and GBD18 in this study). In general, more comprehensive studies based on well-designed sampling plans are required to clarify the peak period of the microbial community succession during environmental fluctuations or parameter adjustments.

4.2 Spatiotemporal conservation of function over microbial community succession

How the microbial community structure relates to the functional performance is a continuous topic in microbial ecology. From the perspective of engineering, stable microbial community under constant operational parameters is expected in WWTPs. However, the results of the present study revealed that the microbial communities drifted among three alternative stable states. In addition, the diverse microbial communities could lead to the uncertainty in WWTP operations. Furthermore, the metagenomics analysis showed that the functional composition of the microbial communities remained relatively conserved under the three alternative states (Fig.3 and Fig.4). This result suggested that the stability of the WWTP could be related to the functionality of a microbial community rather than the taxonomic relationships. In other words, stable nutrient removals were achieved by different microbial communities with nutrient metabolic capacities under seasonal variation in the environmental conditions. In most cases, microbial communities tend to drift into multiple alternative stable states during environmental variation (Shade et al., 2012; Philippot et al., 2021). However, different microbial communities might contain similar functional compositions. Such phenomena have been reported in multiple habitats including wastewater treatment reactors, the ocean, and with plant tissue (Fernández et al., 1999; Louca et al., 2016; 2018). The decoupling of function and taxonomy can be fundamentally explained by the Insurance Hypothesis (Yachi and Loreau, 1999). When the abundance of a certain member decreases, other members in the same ecological niche can complement the required functions, resulting in the functional stability of the ecosystem during community succession. Considering that specific functions related to nutrient removal are required in WWTPs, more attention should be paid to the microorganisms that possess key functions instead of focusing on their taxonomy.

5 Conclusions

A successful application of microbial ecology in wastewater treatment engineering was achieved in the present study. The spatial stability of the microbial communities in the independent corridors was verified. Distinct temporal differences in the microbial communities over spatial stability were observed. Two seasonal community successions appeared in a natural year, which divided the microbial communities into three alternative stable states. The functional composition of the microbial community was conservative compared to their dynamics, resulting in a decoupling between composition and function of the microbial communities in the WWTP. For typical municipal WWTPs with multiple parallel corridors, a single corridor could represent the microbial community characteristics of the WWTP. Future studies should focus on the relationships between the compositional and functional characteristics of microbial communities in WWTPs.

References

[1] Anderson I, Abt B, Lykidis A, Klenk H P, Kyrpides N, Ivanova N. (2012). Genomics of aerobic cellulose utilization systems in Actinobacteria. PLoS One, 7(6): e39331

[2] BastianM H S, JacomyM (2009). Gephi: an open source software for exploring and manipulating networks. In: Proceedings of International AAAI Conference on Weblogs and Social Media. Menlo Park, California: Association for the Advancement of Artificial Intelligence

[3] Buchfink B, Reuter K, Drost H G. (2021). Sensitive protein alignments at tree-of-life scale using DIAMOND. Nature Methods, 18(4): 366–368

[4] Cao J, Zhang T, Wu Y, Sun Y, Zhang Y, Huang B, Fu B, Yang E, Zhang Q, Luo J. (2020). Correlations of nitrogen removal and core functional genera in full-scale wastewater treatment plants: Influences of different treatment processes and influent characteristics. Bioresource Technology, 297: 122455

[5] Chen G, Lin L, Wang Y, Zhang Z, Cao W, Zhang Y. (2023). Unveiling the interaction mechanisms of key functional microorganisms in the partial denitrification-anammox process induced by COD. Frontiers of Environmental Science & Engineering, 17(8): 103

[6] Chen S, Zhou Y, Chen Y, Gu J. (2018). FASTP: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics, 34(17): i884–i890

[7] Chen W, Wei J, Su Z, Wu L, Liu M, Huang X, Yao P, Wen D. (2022). Deterministic mechanisms drive bacterial communities assembly in industrial wastewater treatment system. Environment International, 168: 107486

[8] Cole J R, Wang Q, Fish J A, Chai B L, Mcgarrell D M, Sun Y N, Brown C T, Porras-Alfaro A, Kuske C R, Tiedje J M. (2014). Ribosomal Database Project: data and tools for high throughput rRNA analysis. Nucleic Acids Research, 42(D1): D633–D642

[9] Cui B, Zhang C, Fu L, Zhou D, Hou M. (2023). Current status of municipal wastewater treatment plants in North-east China: implications for reforming and upgrading. Frontiers of Environmental Science & Engineering, 17(6): 73

[10] Dottorini G, Wagner D S, Stokholm-Bjerregaard M, Kucheryavskiy S, Michaelsen T Y, Nierychlo M, Peces M, Williams R, Nielsen P H, Andersen K S, Nielsen P H. (2023). Full-scale activated sludge transplantation reveals a highly resilient community structure. Water Research, 229: 119454

[11] Edgar R C. (2013). UPARSE: highly accurate OTU sequences from microbial amplicon reads. Nature Methods, 10(10): 996–998

[12] Fernández A, Huang S, Seston S, Xing J, Hickey R, Criddle C, Tiedje J. (1999). How stable is stable? Function versus community composition. Applied and Environmental Microbiology, 65(8): 3697–3704

[13] Fu L, Niu B, Zhu Z, Wu S, Li W. (2012). CD-HIT: accelerated for clustering the next-generation sequencing data. Bioinformatics, 28(23): 3150–3152

[14] Guo J, Li J, Chen H, Bond P L, Yuan Z. (2017). Metagenomic analysis reveals wastewater treatment plants as hotspots of antibiotic resistance genes and mobile genetic elements. Water Research, 123: 468–478

[15] Hyatt D, Chen G L, Locascio P F, Land M L, Larimer F W, Hauser L J. (2010). Prodigal: prokaryotic gene recognition and translation initiation site identification. BMC Bioinformatics, 11(1): 119

[16] Jiang X T, Ye L, Ju F, Li B, Ma L P, Zhang T. (2018a). Temporal dynamics of activated sludge bacterial communities in two diversity variant full-scale sewage treatment plants. Applied Microbiology and Biotechnology, 102(21): 9379–9388

[17] Jiang X T, Ye L, Ju F, Wang Y L, Zhang T. (2018b). Toward an intensive longitudinal understanding of activated sludge bacterial assembly and dynamics. Environmental Science & Technology, 52(15): 8224–8232

[18] Ju F, Zhang T. (2015). Bacterial assembly and temporal dynamics in activated sludge of a full-scale municipal wastewater treatment plant. ISME Journal, 9(3): 683–695

[19] Kanehisa M, Goto S. (2000). KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Research, 28(1): 27–30

[20] Kindaichi T, Yuri S, Ozaki N, Ohashi A. (2012). Ecophysiological role and function of uncultured Chloroflexi in an anammox reactor. Water Science and Technology, 66(12): 2556–2561

[21] Knoop S, Kunst S. (1998). Influence of temperature and sludge loading on activated sludge settling, especially on Microthrix parvicella. Water Science and Technology, 37(4−5): 27–35

[22] Li D H, Luo R B, Liu C M, Leung C M, Ting H F, Sadakane K, Yamashita H, Lam T W. (2016). MEGAHIT v1.0: a fast and scalable metagenome assembler driven by advanced methodologies and community practices. Methods, 102: 3–11

[23] Li R, Li Y, Kristiansen K, Wang J. (2008). SOAP: short oligonucleotide alignment program. Bioinformatics, 24(5): 713–714

[24] Louca S, Parfrey L W, Doebeli M. (2016). Decoupling function and taxonomy in the global ocean microbiome. Science, 353(6305): 1272–1277

[25] Louca S, Polz M F, Mazel F, Albright M B N, Huber J A, O’Connor M I, Ackermann M, Hahn A S, Srivastava D S, Crowe S A. . (2018). Function and functional redundancy in microbial systems. Nature Ecology & Evolution, 2(6): 936–943

[26] Luo J, Fang S, Huang W, Wang F, Zhang L, Fang F, Cao J, Wu Y, Wang D. (2022). New insights into different surfactants’ impacts on sludge fermentation: focusing on the particular metabolic processes and microbial genetic traits. Frontiers of Environmental Science & Engineering, 16(8): 106

[27] Philippot L, Griffiths B S, Langenheder S. (2021). Microbial community resilience across ecosystems and multiple disturbances. Microbiology and Molecular Biology Reviews, 85(2): e00026–20

[28] Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, Peplies J, Glöckner F. (2013). The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Research, 41(D1): D590–D596

[29] Shade A, Peter H, Allison S D, Baho D L, Berga M, Burgmann H, Huber D H, Langenheder S, Lennon J T, Martiny J B. . (2012). Fundamentals of microbial community resistance and resilience. Frontiers in Microbiology, 3: 417

[30] Tong J, Tang A, Wang H, Liu X, Huang Z, Wang Z, Zhang J, Wei Y, Su Y, Zhang Y. (2019). Microbial community evolution and fate of antibiotic resistance genes along six different full-scale municipal wastewater treatment processes. Bioresource Technology, 272: 489–500

[31] Wang Y, Ye J, Ju F, Liu L, Boyd J A, Deng Y, Parks D H, Jiang X, Yin X, Woodcroft B J. . (2021). Successional dynamics and alternative stable states in a saline activated sludge microbial community over 9 years. Microbiome, 9(1): 199

[32] Woo S G, Sewell H L, Criddle C S. (2022). Phylogenetic diversity of NO reductases, new tools for nor monitoring, and insights into N₂O production in natural and engineered environments. Frontiers of Environmental Science & Engineering, 16(10): 127

[33] Wu L, Ning D, Zhang B, Li Y, Zhang P, Shan X, Zhang Q, Brown M R, Li Z, Van Nostrand J D. . (2019). Global diversity and biogeography of bacterial communities in wastewater treatment plants. Nature Microbiology, 4(7): 1183–1195

[34] Xia Y, Wen X, Zhang B, Yang Y. (2018). Diversity and assembly patterns of activated sludge microbial communities: A review. Biotechnology Advances, 36(4): 1038–1047

[35] Yachi S, Loreau M. (1999). Biodiversity and ecosystem productivity in a fluctuating environment: the insurance hypothesis. Proceedings of the National Academy of Sciences of the United States of America, 96(4): 1463–1468

[36] Yang Y, Wang L, Xiang F, Zhao L, Qiao Z. (2020). Activated sludge microbial community and treatment performance of wastewater treatment plants in industrial and municipal zones. International Journal of Environmental Research and Public Health, 17(2): 436

[37] Yu J, Tang S N, Lee P K H. (2023). Universal dynamics of microbial communities in full-scale textile wastewater treatment plants and system prediction by machine learning. Environmental Science & Technology, 57(8): 3345–3356

[38] Zhang B, Sun C, Lin H, Liu W, Qin W, Chen T, Yang T, Wen X. (2023). Differences in distributions, assembly mechanisms, and putative interactions of AOB and NOB at a large spatial scale. Frontiers of Environmental Science & Engineering, 17(10): 122

[39] Zhang T, Shao M F, Ye L. (2012). 454 pyrosequencing reveals bacterial diversity of activated sludge from 14 sewage treatment plants. ISME Journal, 6(6): 1137–1147

[40] Zhu J J, Anderson P R. (2017). Effect of long-term successive storm flows on water reclamation plant resilience. Water Research, 111: 1–9

CRediT Authorship Contribution Statement

Zhaoyang Li: Data curation, methodology, validation, visualization, writing–original draft and editing. Liang Zhang: Conceptualization, funding acquisition, writing–review and editing, supervision, project administration. Jinghan Li: Investigation, data curation, visualization. Da Kang: Writing–review and editing. Jialin Li: Writing–review and editing. Shujun Zhang: resources, data curation, investigation. Xiaoyu Han: resources, data Curation. Bin Ma: Funding Acquisition, writing–review, project administration. Yongzhen Peng: Resource, Supervision, Project administration, Writing–Review.

Acknowledgements

This work was financially supported by the National Natural Science Foundation of China (Grant Nos. 52122005 and U23A20675).

Conflict of Interests

Yongzhen Peng is Editorial Board Members of Frontiers of Environmental Science & Engineering. The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Electronic Supplementary Material

Supplementary material is available in the online version of this article at https://doi.org/10.1007/s11783-025-1943-5 and is accessible for authorized users.

Open Access

This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

RIGHTS & PERMISSIONS

2025 The Author(s) 2025. This article is published with open access at link.springer.com and journal.hep.com.cn