3.1 Seasonal variations in the biofilm microbiome
The response of the biofilm microbiome to seasonal variations in a denitrification biofilter of the Maidao wastewater treatment plant was investigated. DNA was extracted from summer and winter biofilm samples, respectively, and the two metagenomes were sequenced using Illumina HiSeq. A total of 9.44 million raw paired-end reads were generated with a length of 150 bps, 9.33 million high-quality clean reads were retained after quality control. 1,064,671 contigs for gene prediction were obtained using the multiple hybrid assembly strategy, and 1,644,199 gene sequences were yielded (Table S1). Bacteria were overwhelmingly dominant among biofilm microorganisms, with over 99.5% of the gene sequences annotated as bacteria and 0.05% of the residual as viruses, eukaryote, and archaea.
For the Bacteria domain, a total of 79 phyla were detected in all samples. As shown in Fig. 1a, the top 5 phyla were Proteobacteria (69.4%), Bacteroidetes (24.5%), Actinobacteria (0.9%), Chloroflexi (0.7%) and Firmicutes (0.6%). The bacteria phyla mentioned above are commonly found in municipal WWTPs worldwide with different biological treatment processes, and Proteobacteria and Bacteroidetes account for the largest proportion of microorganisms (de Celis et al., 2020; Wang et al., 2016). Proteobacteria is a highly diverse phylum that involved in various anaerobic reactions (Greses et al., 2018). A study of nirS-Type Denitrifying Bacterial Community from 18 WWTPs found that Proteobacteria was dominated in all samples (Zhang et al., 2019). Bacteroidetes phylum has a strong metabolic capacity for proteins, lipids and other complex macromolecular organic matter. Chloroflexi is involved in the anaerobic degradation of glucose and complex organic matter (e.g., cellulose) and plays a role in nitrite reduction (Hug et al., 2013).
Although the major phyla were the same in samples from different seasons, their respective abundances were quite different. It is noteworthy that the percentage of Bacteroidetes phylum was 28.3% in the summer sample, 71.7% in the winter sample (Table S2). In contrast, Chloroflexi phylum was accounted for 90.7% in the summer sample and 9.3% in the winter sample. The higher relative abundance of Bacteroidetes in winter samples might be explained by some subdivisions of Bacteroidetes affiliated with psychrophilic and psychrotolerant organisms, which are enriched in cold environments and develop mechanisms of cold resistance (Ducey et al., 2010). Within Chloroflexi, Anaerolineae was the predominant class, followed by Chloroflexia. The relative abundance of Anaerolineae was reported to be decreased with temperature decrease (Damodara Kannan et al., 2020). Moreover, Chloroflexia's representative nutritional mode is photosynthesis (Ward et al., 2018), which may grow on the external layer of the biofilm and is thus strongly influenced by temperature, consistent with the previously reported greater abundance of green algae with increasing temperature (Uribe-Lorio et al., 2019). The abundance of Proteobacteria, Actinobacteria and Firmicutes varied little in the winter and summer samples, showed strong stability to seasonal changes. Previous studies on activated sludge found that the abundance of Proteobacteria in activated sludge samples was much higher in winter than in other seasons, while the abundance of Bacteroidetes was lower than in other seasons (Ju et al., 2014; Liu et al., 2016). A possible explanation for the different results is the structural difference between biofilm and activated sludge.
Protein sequences from metaproteomic were compared with NR database to obtain taxonomic annotations and abundance information. Subsequently, taxonomic annotation results of metagenomics and metaproteomic were compared at the class level for different season samples. Betaproteobacteria class and Alphaproteobacteria class were dominant bacteria (> 55%) in the classification results of both seasons and both annotation methods. For both season samples, Betaproteobacteria were much more abundant in the metaproteomic classification than in the metagenomic, while the abundance of Alphaproteobacteria was lower than that in the metagenomic classification (Fig. 1b & c). Therefore, the majority of functional microorganisms in denitrification biofilter belonging to Betaproteobacteria. Betaproteobacteria are the majority of denitrifying bacteria found in denitrification processes (Wang & Chu, 2016). It is noteworthy that Flavobacteriia, an essential member of biofilm microorganisms, was significantly abundant in the winter sample, in agreement with a previous study that Flavobacteriia increased in abundance under cold water conditions (Kraemer et al., 2020), the enhanced abundance of Flavobacteriia in winter ensured the removal performance of nitrogen (Zhou et al., 2020).
3.2 Proteomic expression profiling and differentially expressed proteins identification
The label-free protein expression quantification method was used to evaluate the effect of seasonal variations on denitrification biofilter microorganisms. The LC-MS/MS analysis generated 99,307 Peptide-Spectrum Matches, 19,309 peptides and 27,305 unique peptides in the database of 375,757 spectra, and a total of 7,767 proteins were identified (Fig. S1a). The gene ontology (GO) analysis grouped all proteins into three major categories: biological process (BP), cellular component (CC), and molecular function (MF) (Fig. S1b). The catalytic activity and binding were the dominant categories in MF ontology. Among the BP ontology, the majority of proteins were attributed to cellular processes and metabolic processes. And cellular anatomical entity in cellular component had large numbers of proteins. Based on the annotation result of KEGG analysis, most proteins were primarily involved in global and overview maps, energy metabolism and carbohydrate metabolism pathways belonging to the metabolism category (Fig. S1c). The Principal Component Analysis (PCA) score plot exhibited that the biofilm samples from summer and winter were significantly separated, indicating that the proteomic expression profiling of the biofilm shifted with the seasons (Fig. 2a).
To further access the dynamic changes of proteins in biofilms during summer and winter, 2835 differentially expressed proteins (DEPs) were identified from the two groups of samples using a Screening Criteria of fold change ≥2.00 or ≤0.50 and P < 0.05. There were 1596 up-regulated and 1239 down-regulated DEPs in summer samples compared with the winter samples (Fig. 2b). The top 15 DEPs ranked according to P-value were listed in Table S3, which may play an important role in seasonal adaptation and may be further investigated as candidate genes in the future.
GO enrichment analysis based on DEPs was conducted to clarify the functional differences of samples between seasons (Fig. 2c). The results showed that in response to seasonal variations, GO terms belonging to BP and MF ontologies had high enrichment factors with strong significance (P<0.05). In BP ontology, cellular process, ATP metabolic process, small molecule metabolic process, organic substance metabolic process and biosynthetic process were significantly enriched. In MF ontology, most of the DEPs were enriched in carbohydrate derivative binding, small molecule binding, heterocyclic compound binding and organic cyclic compound binding. Detailed information about GO enrichment analysis is shown in Table S4. In addition, the most enriched pathways related to the binding of ribonucleotide, nucleotide, and nucleoside phosphate. Nucleotide-binding proteins are essential protein families involved in various critical cellular processes such as cell signaling (Xiao & Wang, 2016). In addition, several GO terms related to transporter activity and biosynthetic and metabolic processes of organic acids were also statistically significantly enriched.
KEGG pathway enrichment analysis was further conducted based on DEPs. DEPs were assigned to 247 KEGG pathways, with most of them (116) related to metabolism (M), 48 pathways to human disease (HD), 28 to organismal systems (OS), 24 cellular processes (CP) pathways, 17 pathways for environmental information processing (EIP), and 14 pathways for genetic information processing (GIP). Among the enriched KEGG pathways (Fig. 2d), seven pathways, including quorum sensing, two−component system, ribosome, benzoate degradation, butanoate metabolism, citrate cycle (TCA cycle), and cysteine and methionine metabolism were significantly enriched (p < 0.01). Except for the two-component system, the other six pathways are mainly enriched by up-regulated DEPs. The ribosome and cysteine and methionine metabolism had the highest significant enrichment levels (p = 0.0007). Quorum sensing plays a vital role in arranging bacterial gene regulation, participating in the regulation of biofilm formation and extracellular polymer production processes, enabling bacteria to adapt to changes in environmental conditions (Acet et al., 2021). Temperature is a significant factor affecting Quorum sensing (Wang et al., 2019), and a study of anaerobic ammonium oxidation biofilm reactor found that a decrease in temperature weakened quorum sensing and deteriorated nitrogen removal performance (Liu et al., 2021). The enrichment of the quorum sensing pathway implies its role in the adaptation of biofilm microbiome to seasonal variations. The two-component system is a fundamental stimulus-response coupling mechanism that enables microorganisms to sense and respond to changes in environmental conditions (Kuang et al., 2018). Following seasonal variations, some proteins involved in the two-component system pathway were differentially regulated. There was significant regulation of proteins involved in the degradation/metabolism of organic acids (such as Benzoate, Butanoate) and amino acids (Cysteine and methionine metabolism). Benzoate conversion into TCA cycle intermediates via benzoyl-CoA (Garrido-Sanz et al., 2018). Butanoate metabolism is directly linked to the TCA cycle via succinate and indirectly linked via pyruvate and acetyl-CoA (Goudarzi et al., 2015). Similarly, GO terms associated with the organic acid metabolic process were significantly enriched. In bacteria, methionine is synthesized from aspartate and cysteine is metabolized to pyruvate through various pathways (Park & Lee, 2010; Samuilov et al., 2018), whereas aspartate and pyruvate are important compounds in the TCA cycle. All these pathways are directly or indirectly linked to the TCA cycle, through which their degradation products generate the reducing agent NADH to be used as an electron donor for denitrification (Emmanuel et al., 2019a). The nitrogen metabolic pathway was also one of the significant pathways with 181 up-regulated and 61 down-regulated DEPs, which validates the enhanced efficiency of denitrification in summer. And the most significant pathway, ribosomes, was also significantly enriched in the up-regulated DEP (p = 0.007, Fig. S2), demonstrating the extensive enhancement of protein expression and synthesis in summer to meet the greater demands of cell growth (Zhong et al., 2019).
3.3 DEPs involved in the nitrogen metabolism
Based on the metagenomic results, the nitrogen metabolism pathway was the most active in the pre-denitrification biofilter, with the denitrification module and DNRA module, the most represented functions; and there were significant differences in denitrification, nitrate/nitrite transport systems, and nitrate assimilation modules (Fig. S3a, b). According to the results of metaproteomics, the expression of proteins for nitrogen metabolism was lower than carbon metabolism and Glyoxylate and dicarboxylate metabolism in summer. In contrast, the expression of proteins for nitrogen metabolism decreased significantly in winter, the sixth in the ranking of the abundance of all pathways (Fig. S3c). The above metagenomic and metaproteomic results showed that the nitrogen metabolism pathway differs significantly between summer and winter.
To further investigate the variation of nitrogen metabolism in different seasons, key genes involved in the nitrogen metabolism pathway were identified by KEGG analysis of DEPs. The KEGG pathway map (Fig. 3a) showed that the nitrogen metabolism pathway was enriched mainly by up-regulated proteins. EC 188.8.131.52 is known as nitrate reductase (alpha subunit, beta subunit and gamma subunit, NarGHI), and EC 184.108.40.206 is the periplasmic nitrate reductase (NapA), both responsible for the conversion of nitrate to nitrite. Nitrite is reduced to nitric oxide by nitrite reductase (NirKS), denoted as EC 220.127.116.11 in the KEGG pathway map. EC 18.104.22.168 marks the conversion of nitric oxide to nitrous oxide by nitric oxide reductase subunit B (NorB). The last step of denitrification is catalyzed by nitrous-oxide reductase (NosZ), denoted as EC 22.214.171.124. The key enzymes mentioned above were commonly expressed in both seasons, with the abundance of up-regulated proteins was much higher than that of down-regulated proteins, and the expression of denitrification related enzymes in summer was at least as twice high as in winter samples. It could explain the decreased denitrification performance of the pre-denitrification biofilter system in winter (Table 1). In addition, nitrate transporter (Nrt) was entirely enriched by up-regulated proteins to uptake extracellular nitrate. The ammonia input-related enzymes such as carbonic anhydrase (Can) and glutamate dehydrogenase (Gdh) were enriched by up-regulated proteins. Glutamine synthetase (Gln) was highly enriched by down-regulated proteins and was expressed more vigorously in winter, and glutamate synthase (Glu) was slightly more expressed in summer than in winter.
Four potential metabolic reactions in nitrogen metabolism were identified, including denitrification, dissimilatory nitrate reduction, assimilatory nitrate reductase and nitrification (Fig. 3b). Genes encoding the alpha and beta subunits of nitrate reductase (narG and narH) were generally actively expressed in summer samples. However, the genes LPB072_22010, Tmz1t_2635 and Tchl_2559 encoded the alpha subunit of nitrate reductase, AzCIB_2182 and LPB072_22005 encoding the beta subunit, as well as the gene Dtpsy_0430 encoding gamma subunit, were significantly better expressed in winter samples. These bacteria belong to the family Comamonadaceae and Zoogloeaceae of the order Betaproteobacteria, core microorganisms in low-temperature bioreactors (Gonzalez-Martinez et al., 2018; Huang et al., 2014), and can therefore accumulate at low temperatures in winter to improve the performance of denitrification. The genes encoding narG and narH can also encode nitrite oxidoreductase. However, under the premise of anaerobic conditions in the denitrification filter, nitrite oxidoreductase cannot be expressed and nitrification reaction does not occur (Kampschreur et al., 2009). In the DNRA pathway, the first step of nitrate reduction is consistent with the denitrification pathway. The presence of genes encoding nirBD and nrfAH was detected by metagenomics (Tian & Wang, 2021), but the expression of related proteins was not detected, thus DNRA pathway excluded in this study. Nitrite reductase was the most actively expressed enzyme in this study, with NirS accounting for 56% of the DEPs involved in nitrogen metabolism. The species composition of the denitrifying bacterial community encoding NirS in summer was significantly different from that in winter. The most abundant bacteria, such as Thauera chlorobenzoica (Tchl_2591), Dechloromonas aromatica (Daro_3323) and Acidovorax sp. JS42 (Ajs_1912), all encoded and expressed genes for nitrite reductase (NirS) and were significantly up regulated in summer. The expression of nirS-encoding genes was much higher in summer than that in winter, and the expression of nirK-encoding genes showed the opposite pattern. Genes encoding NirK showed low expression and was significantly downregulated (p < 0.01) in summer samples, suggesting that the NirK may play an important role in nitrite reduction in winter. The Competition between denitrifying bacteria encoding NirS and NirK could be responsible for the low expression of genes encoding NirK in summer (Shi et al., 2019). Based on the expression of genes encoding nitric oxide reductase, nitrous oxide in denitrification biofilters may be produced mainly by Elizabethkingia, Rubrivivax and Poseidonibacter in winter. Notably, a substantial proportion of the genes encoding NosZ were not significantly differentially expressed in the samples of the two seasons, whose protein expression was ~1.7-fold higher and ~3.0-fold higher than that of the up-regulated and the down-regulated protein. It is consistent with a previous study in which season had less influence on the nosZ community (Ross et al., 2020). Certain bacteria, such as Capnocytophaga (chg: AXF12_04010), Polaromonas (pna: Pnap_3762) and Hydrogenophaga (hyr: BSY239_3898) showed the capability to express assimilatory nitrate reductase in this study. In addition, Polaromonas and Hydrogenophaga, both belonging to Comamonadaceae class, can also produce Nrt to transport extracellular nitrate and nitrite to the intracellular compartment. Even though metagenomics detected multiple nitrogen metabolic pathways, the metaproteomic results validated only two pathways: assimilative nitrate reduction and denitrification. And the genes for denitrification showed the highest protein expression.
3.4 Metabolomic analysis
Microbial cells can adapt to certain shifts in environmental conditions by exhibiting corresponding responses through changes in relevant metabolic processes (Dai et al., 2017). It is crucial to reveal community metabolic profiling and identify the relationship between metabolic activities and seasonal variations. By resolving the end-products and intermediates of cellular metabolism, metabolomics is supposed to be the most sensitive indicator of the phenotype of the whole community and allows to derive critical metabolic processes (Roume et al., 2015). Therefore, to complement the results of metagenomic and metaproteomic analyses, metabolomic analysis was performed to investigate the changes in metabolites of the two groups of biofilm samples.
3.4.1 Biofilm metabolite profile and impact of season variations on community metabolome
Using GC-MS, a total of 248 metabolites were identified in biofilms. PCA was performed among the QC samples and tested samples. The score plot (Fig. 4a) revealed a clear separation of biofilm samples from the winter and summer groups along with the first principal component (PC1), explaining 35.4% of the variance within the data. The QC samples were well grouped, tightly clustered with each other and separated from the tested samples, which confirmed the good analytical performance and data quality of GC-MS and supported that the separation between groups was related to biological variations (Cala et al., 2018). The corresponding loading plot (Fig. 4b) reflected the distribution of metabolites in the biofilm samples, with metabolites further away from the center (such as Lanosterol, Glycerol 3-phosphate and Phosphogluconic Acid) contributing more to the distinction of the samples. Further multivariate analysis using supervised OPLS-DA model was performed to predict variables of interest to discriminate the summer group from the winter group (fa), and the two principal components explained 49.3% of the variance. The corresponding stochastic model of R2 and Q2 values showed that the intercept of the Q2 was less than 0 and the model was not overfitted (Fig. S4b). These results indicated that seasonal variations induced alteration in the biofilm metabolite profile.
The student’s t-test (Fig. S4c) showed that a total of 80 metabolites were significantly different (p < 0.05) between the two groups, together with the screening criterion of OPLS-DA model VIP > 1, a total of 66 differential metabolites were identified. All differential metabolites were assigned to the KEGG database, and 28 of which were matched and classified: Peptides (8), Organic acids (6), Carbohydrates (4), Hormones and transmitters (3), Lipids (1), Nucleic acids (1), Steroids (6) and Vitamins and Cofactors (2) (Fig. S4d). Pathway topology analysis was performed to reveal the perturbed metabolic pathways based on the identified differential metabolites and changes in their concentrations. Eight pathways with pathway impact values above 0.1 were selected as perturbed metabolic pathways (Fig. 4d, Table S5), including pyruvate metabolism, butanoate metabolism, glycerolipid metabolism, glycerophospholipid metabolism, pantothenate and CoA biosynthesis and alanine, aspartate and glutamate metabolism, valine, leucine and isoleucine biosynthesis and citrate cycle (TCA cycle). Further enrichment analysis (Fig. 4e) confirmed that pyruvate metabolism, butanoate metabolism, pantothenate and CoA biosynthesis, alanine, aspartate and glutamate metabolism and TCA cycle were significantly affected (p < 0.05) by seasonal variations. Fold change thresholds were included based on the selected differential metabolites, three metabolites were up-regulated (foldchange > 1.5) in the summer group, and five metabolites were down-regulated (foldchange < 0.67). These eight metabolites were significantly regulated and selected as potential biomarkers (Fig. 4c). The metabolic profiling association between potential biomarkers and perturbed metabolites was further analyzed (Fig. 5b). Changes in potential biomarkers were found to be correlated with changes in most metabolites involved in the perturbation pathway (|r| > 0.4, p < 0.05), suggesting that changes in biomarkers can reflect reprogramming of metabolic pathways during seasonal changes.
Notably, most perturbed pathways were directly related to nitrogen metabolism, such as the TCA cycle and the metabolism of pyruvate, butanoate, alanine, aspartate, and glutamate. Several amino acids and organic acids related to TCA cycle were significantly perturbed by seasonal variation (Fig. 5a). Amino acids are essential to bacteria and components of all proteins, playing a critical role in microbial physiology processes, such as forming bacterial cell walls, synthesizing key enzymes, and being utilized as carbon source and nitrogen source (Aliashkevich et al., 2018; Cava et al., 2011). Aspartate, L-Aspartic acid, L-Tyrosine and L-Phenylalanine were up-regulated, while ornithine and L-Histidine were down-regulated. In addition to the regulation of amino acids, five organic acids (TCA cycle intermediates), including pyruvic acid, isocitric acid, butanedioic acid, fumaric acid and malic acid were significantly elevated in winter. The denitrification pathway could be associated with the degradation of amino acids and organic acids, with their degradation products going through the TCA cycle and producing NADH used as electron donor for denitrification (Emmanuel et al., 2019b). The biosynthesis and metabolism of amino acids and organic acids are closely linked through the TCA cycle (Zhao et al., 2019), therefore the metabolite alteration caused by seasonal variation affect both carbon metabolism and nitrogen metabolism. Moreover, many amino acids can act as compatible solutes due to their kosmotropicity, which allows them to stabilize membranes and proteins against environmental stresses (He et al., 2016; Xu et al., 2020). The accumulation of TCA cycle intermediates, ornithine and L-Histidine in winter may be a defense or coping strategy against the cold temperature (He et al., 2016).
3.4.2 Correlation study of biofilm microbiome and community metabolites
An integrated analysis of metabolomic and metaproteomic was performed for the two groups to explore and define the relationship between microbial communities and metabolic functions. Correlation analysis of the association between microbial communities with differential metabolites provided insight into the composition and function of microbial communities (Fig. 5c). Bacteroidetes with higher abundance in winter was negatively associated (p < 0.05) with up-regulated metabolites such as L-Aspartic acid, D-Fructose and mevalonic acid. While Chloroflexi with higher abundance in summer was negatively associated with down-regulated metabolites such as ornithine, L-Histidine and isocitric acid, suggesting that alterations in metabolite profiling of microbial communities reflect shifts in microbial community dynamics. Among the potential biomarkers, nicotinamide was positively associated with Streptophyta and Chlorophyta belonging to Viridiplantae, as nicotinamide is a metabolite synthesized by all plants and can react with 2-oxoglutamic acid to form L-glutamic acid, which is the largest contributor to the intermediate fluxes of the TCA cycle (Ohata et al., 1993; Takahashi et al., 1999). In addition, nicotinamide was also positively associated with fungi such as Ascomycota and Mucoromycota, and its presence in certain fungi was confirmed by previous studies (Laffont & Arnoux, 2020). Only trace amounts of lanosterol were detected in summer samples, while it was significantly elevated in winter. Lanosterol was significantly associated with Chordata and Rotifera, which belong to metazoan, and it is the initial biosynthetic sterol in the metazoa sterol pathway (Kodner et al., 2008). Other potential biomarkers were also positively or negatively associated with eukaryotes, such as fungi, metazoa and viridiplanta, as eukaryotes may respond to seasonal changes with different resistance and sensitivity (Zhang et al., 2020). Our study revealed that seasonal variations directly impact microbial communities and metabolites and that there is a correlation between them. However, the relationship between metabolites and microbial communities in this study was purely correlational without further controlled experiments, and the control experiment should be conducted to find a detailed relation in the future.