Enhanced long-term nitrogen removal by organotrophic anammox bacteria under different C/N ratio constraints: quantitative molecular mechanism and microbial community dynamics

Abstract

The anaerobic ammonium oxidation (anammox) process has mainly been applied to NH₄⁺–N-rich wastewater with very low levels of organic carbon (<0.5 g COD per g N). Some anammox bacteria species have the capacity to oxidize organic carbon with nitrate as the electron acceptor. However, little is known about the organotrophic anammox nature of “Ca. Brocadia sinica”. To elucidate the metabolic versatility and microbial succession of “Ca. Brocadia sinica” under TOC/TN stress conditions, the influence of TOC/TN ratios on the nitrogen transformation pathway and the dynamics of microbial communities were investigated. It was found that an appropriate TOC/TN ratio (<0.2) could promote the anammox activity over the short-term. Meanwhile, “Ca. Brocadia sinica” had higher tolerance to higher TOC/TN (>0.4) abiotic stresses. Mass balance indicated that organotrophic anammox could outcompete denitrifiers under a TOC/TN ratio of 0.1–0.2. Quantitative response relationships and pathway analysis revealed that (AOA amoA + AOB amoA + anammox + nrfA)/bacteria, nrfA/(narG + napA), and nrfA were the key functional gene groups determining the organotrophic anammox contribution. Additionally, MiSeq sequencing showed that Planctomycetes, Proteobacteria, Chloroflexi, and Chlorobi were the most abundant phyla in the organotrophic anammox system. Furthermore, higher TOC/TN ratios (>0.40) could result in the community succession of anammox species, in which “Ca. Jettenia caeni” and “Ca. Kuenenia stuttgartiensis” were the dominant organotrophic anammox bacteria species. Overall, combined analyses revealed that the coupling of anammox, DNRA (organotrophic anammox), and denitrification comprised the primary pathway that accounted for TOC and nitrogen removal.

---

1. Introduction

Anaerobic ammonium-oxidizing (anammox) bacteria, which were discovered in a denitrifying bioreactor and mediated by Plantomycete bacteria,^1,2 have the unique ability to combine ammonium and nitrite to form N₂. Understanding of the microbial nitrogen cycle has therefore been fundamentally altered by the discovery of anammox bacteria.³ The anammox-related process, which is a lithoautotrophic nitrogen removal process, has been successfully applied in the treatment of ammonia-rich wastewater with low COD (chemical oxygen demand) : N (nitrogen) ratios due to its cost-effective and energy-efficient qualities.^4,5

Nevertheless, the slow growth rate, low cell yield, and sensitivity to environmental conditions which characterize anammox bacteria have presented major obstacles to the broader application of anammox-related processes. For instance, in a mainstream and side-stream nitrogen removal process, anammox bacteria were unable to avoid the influences of volatile fatty acids (VFA) (expressed as total organic carbon), which existed in large volume in municipal and industrial WWTPs.⁶ Although the addition of organic matter had significant effects on anammox bacteria in many studies, there has been no consensus on which TOC (total organic carbon) to TN (total nitrogen) ratios inhibit/or affect the anammox bacteria.^7–9 In addition, investigations have reported that the most critical point in the competition between autotrophic anammox, organotrophic anammox, and heterotrophic denitrification for TOC is the C : N ratio in the influent.^10,11 Moreover, detailed evidence describing the contribution of organotrophic or mixotrophic anammox processes to nitrogen removal in the presence of TOC has yet to be presented.

To date, six anammox bacteria genera have been detected and proposed using 16S and 23S rRNA gene sequencing: “Ca. Brocadia sinca”, “Ca. Anammoxoglobus propionicus”, “Ca. Jettenia caeni”, “Ca. Kuenenia stuttgartiensis”, “Ca. Scalindua sp.”, and “Ca. Brocadia anammoxidans”.^12,13 A few recent studies^10,14,15 have reported that some anammox bacteria species, including Ca. Jettenia asiatica, Ca. Ananimoxoglobus propionicus, Ca. Brocadia fulgida, and Ca. Kuenenia stuttgartiensis have the capacity to oxidize acetate and propionate. In addition, phylogenetic classification of anammox 16S rRNA¹² revealed that “Ca. Brocadia sinica” was closely related to the group “Ca. Brocadia Scalindua sp.” and “Ca. Brocaida anammoxidans”, which were organotrophic. Thus, it could be speculated that TOC might have been utilized as electron donors by “Ca. Brocadia sinica” and this phylum could success to other organotrophic anammox species. However, whether the organotrophic nature of “Ca. Brocadia sinica” participated in TOC oxidation remains unverified. This inspired us to explore the adaptability and organotrophic nature of “Ca. Brocadia sinica” in the presence of TOC.

In addition, anammox sludge in wastewater treatment plants is a highly complex system of eukaryotes (protozoa, fungi, and microalgae),^16–19 bacteria, archaea, and viruses, in which bacteria are dominant. Molecular biological methods which have been applied to explore the microbial structures in anammox-related system include clone library of 16S rRNA genes, denaturing gradient gel electrophoreses (DGGE) analysis and fluorescence in situ hybridization (FISH). With the recent rapid development of the next-generation sequencing, high-throughput sequencing has been received more attention. 454 pyrosequencing²⁰ and Illumina high-throughput sequencing²¹ have been applied to the investigation of microbial communities in lab-scale and pilot-scale anammox-related systems. However, little is known about the dynamics of microbial communities and functional genes under TOC/TN constraints. Knowledge of microbial community structures and the links to different TOC stresses is therefore essential for understanding the profile of the organotrophic anammox process.

Furthermore, nitrification, denitrification, and anammox could co-exist in anammox-related systems when the TOC was present.²² These nitrogen removal processes involve several functional genes which have played key roles in microbial nitrogen cycling, including anammox 16S rRNA, archaea ammonia monoooxygenase (AOA-amoA), ammonia monoooxygenase (AOB-amoA), nitrite oxidoreductase (nxrA), periplasmic nitrate reductase (napA), membrane-bound nitrate reductase (narG), dissimilatory nitrate reductase (nrfA), copper-containing nitrite reductase (nirK), nitrite reductase (nirS), and nitrous oxide reductase (nosZ).^23,24 Nevertheless, the quantitative response relationship between nitrogen transformation rates and functional genes are unknown in the anammox system. Furthermore, the corresponding dynamics of microbial community structures and functional gene groups on the quantitative molecular level in the organotrophic anammox process are still unclear.

The present study is the first to systematically investigate the microbial community structure dynamics and quantitative molecular mechanisms of nitrogen transformation in anammox systems under different TOC stress constraints. Based on these arguments, this study was performed for the following purposes: (1) to assess the impacts of different TOC/TN ratios on organotrophic anammox growth rates and activity using batch experiments; (2) to evaluate the long-term adaptation of “Ca. Brocadia sinica” and contribution of organotrophic anammox bacteria species under different TOC/TN constraints; (3) to explore the quantitative response relationships between nitrogen transformation rates and functional gene groups, and (4) to investigate the dynamics of nitrogen-related microbial communities under TOC/TN stress conditions.

2. Methods

2.1. Kinetic evaluation and treatment performance of an anammox bioreactor under different TOC/TN ratio constraints

The anammox biomass used in this study was taken from a laboratory-scale sequencing batch reactor (SBR) with a working volume of 2.6 L. The anammox–SBR system has been operated for over 18 months under mesophilic conditions (35 ± 2 °C) with hydraulic residence time (HRT), influent NH₄⁺–N, and NO₂⁻–N concentrations of 4 h, 200 mg L⁻¹ and 220 mg L⁻¹, respectively. According to Shu et al.,²¹ the “Ca. Brocadia sinica” phylum comprised the dominant anammox bacteria in this SBR system. Prior to the addition of TOC, the stock solution of organic carbon sources was prepared by mixing it with acetate and propionate at a volume ratio of 1 : 1. In the case of batch experiments, 10 mL of anammox biomass were dispensed into 60 mL of liquid in 100 mL serum vials with different TOC/TN ratios of 0, 0.05, 0.10, 0.20, 0.41, 0.61, and 0.82 (Table 1). Finally, the concentration values of mixed liquor volatile suspended solids (MLVSS), NH₄⁺–N and NO₂⁻–N in each vial were 2850 mg L⁻¹, 80 mg L⁻¹, and 96 mg L⁻¹, respectively. The above experimental procedures were performed in an anaerobic glove box. All the experimental vials were thereafter incubated at 32 °C and shaken at a speed of 120 rpm in the dark. The water samples for further kinetic analysis were taken from the vials hourly over 8 hours.

Table 1 Short- and long-term experiment conditions

	NH4^+–N (mg L^−1)	NO2^−–N (mg L^−1)	TOC (mg L^−1)	TOC : TN
Batch experiments
Batch test 1	80	96	—	—
Batch test 2	80	96	8.99	0.05
Batch test 3	80	96	17.98	0.10
Batch test 4	80	96	35.96	0.20
Batch test 5	80	96	71.91	0.41
Batch test 6	80	96	107.87	0.61
Batch test 7	80	96	143.82	0.82
Long-term experiments
Phase I (1–62 days)	190	220	—	—
Phase II (63–79 days)	190	220	41.20	0.10
Phase III (80–95 days)	190	220	82.40	0.20
Phase IV (96–106 days)	190	220	164.79	0.40
Phase V (107–120 days)	190	220	329.59	0.80
Recovery (121–140 days)	70	84	—	—

---

To investigate the effects of different TOC/TN ratios on “Ca. Brocadia sinica”, the specific anammox activity (SAA) and specific anammox growth rates (μ_AN) were measured according to the following models.²⁵

where K_S is the half saturation constant; S is the TOC concentration; K_I is the inhibition constant; and SAA_max and μ_{AN max} are the maximum specific anammox activity and specific anammox growth rates, respectively.

For long-term experiments, ∼900 mL of seeding sludge were derived from the initial SBR system and then cultivated in a new SBR system with a working volume of 2.6 L. The new SBR system was operated under the same mesophilic conditions with mineral medium and trace element solution.¹ The new anammox–SBR system was run in a 12 hour-cycle and fed with 190 mg L⁻¹ NH₄⁺–N and 220 mg L⁻¹ NO₂⁻–N. Because the anammox bacteria tolerate to lower DO concentrations (<0.5 mg L⁻¹) well and the dissolved oxygen (DO) in anammox–SBR system was consumed by some of the aerobic bacteria, the dissolved oxygen (DO) was not controlled in the anammox–SBR system. More specifically, the dissolved oxygen (DO) was not controlled in any of the phases, and it measured 1.0–1.3 mg L⁻¹ in the influent and 0.08–0.3 mg L⁻¹ in the effluent. After the adaptation stage, the experimental batches with different ratios of TOC/TN had values of 7.2 mL, 14.4 mL, 28.8 mL, and 57.6 mL (7.24 g TOC per L) mixed solution of sodium acetate and propionate. These were added into the reactor automatically at the end of each feeding period to maintain influent TOC/TN ratios of 0.10, 0.20, 0.40, and 0.80 (Table 1).

2.2. DNA extraction, PCR amplification and Illumina MiSeq sequencing

After a start-up period, ∼0.5 g of anammox biomass was collected from the SBR system at the end of each phase. Then genomic DNA was extracted by using the FastDNA® SPIN Kit for Soil (Mp Biomedicals, Illkirch, France) according to the manufacturer's protocols. DNA purity was checked in agarose gel (1.2%) and its concentration was determined with Nanodrop Spectrophotometer ND-1000 (Thermo Fisher Scientific, USA).

Before sequencing, genomic DNA was amplified by PCR using primer sets 338F (5′-barcode-ACTCCTACGGGAGGCAGCAG-3′) and 806R (5′-barcode-GGACTACHVGGGTWTCTAAT-3′) for the V3–V4 regions of the bacterial 16S rRNA. The amplification reactions were performed in triplicate using the previously described primers and protocols.²⁶ After amplification, the PCR products were purified with the AxyPrep DNA Gel Extraction Kit (Axgen, USA) and quantified with a QuantiFluor™-ST (Promega, USA) according to the instructions. Then the three individual PCR products were pooled in equimolar ratios and amplicon libraries were constructed before sequencing. Finally, the amplicon libraries were prepared and then run on the MiSeq Illumina platform (Shanghai Personal Biotechnology Co., Ltd, Shanghai, China). All the raw sequences have been deposited into the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) database (Accession number: SRR2962328).

2.3. Sequence processing and bioinformatics analysis

Following sequencing, all the raw reads were initially merged using FLASH (Version 1.2.11, http://ccb.jhu.edu/software/FLASH/) software, and then Trimmomatic (Version 0.33, http://www.usadellab.org/cms/?page=trimmomatic) was used to remove barcodes, primers, and low quality reads. After filtration, the remaining set of effective sequences were clustered into operational taxonomic units (OTUs) using Usearch (Version 8.1, http://drive5.com/uparse/) at a 97% similarity setting. Then, the taxonomy was assigned using an RDP classifier (Version 2.2, http://sourceforge.net/projects/rdp-classifier/) via Silva SSU database (Release 119, http://www.arb-silva.de) with a set confidence threshold of 70%. Furthermore, appropriate subsample depth analyzing was conducted to avoid unequal sampling depth biases during comparison of microbial diversity and to ensure adequate sample depth while retaining the lowest sequences.²⁷ Subsequently, alpha diversity statistics including Chao 1 estimator, ACE estimator, Shannon index, Simpson index, Good's coverage, and rarefaction curves at a distance of 0.03 (equivalent to 97% similarity), were calculated for five samples using the Mothur program (Version 1.30.1, http://www.mothur.org/wiki/Main_Page).²⁸ On the basis of the UniFrac metric, beta diversity statistics including cluster analysis (CA), (un)weighted UniFrac distance metrics, and principal co-ordinates analysis (PCoA) were also conducted to evaluate the similarities and differences of five samples.

2.4. Quantitative real-time PCR

Compared with RNA-based techniques, DNA-based methods are more robust alternative for measuring the total microbial growth rates and activity for the given environmental conditions because they account for all active, dormant, and deceased cells.²⁹ Thus, quantitative real-time PCR, a DNA-based technology, was used to explore the “key players” in the nitrogen removal process and its quantitative molecular mechanism. The absolute abundance of bacteria 16S rRNA, anammox 16S rRNA, AOB amoA, AOA amoA, nosZ, nirS, nirK, narG, napA, and nrfA were quantified in triplicate by Mastercycler ep realplex (Eppendorf, Hamburg, Germany) based on the SYBR Green II method using the described primers and protocols.²⁶ The plasmids containing bacteria 16S rRNAs, anammox bacteria 16S rRNAs, Nitrobacter 16S rRNAs, Nitrospira 16S rRNAs and other functional genes (i.e. amoA, nosZ, nirS, nirK, narG, napA, nrfA, mcrA, dsrA) were manufactured by Sangon Biological Technology Company (Sangon Biological Engineering Co., China).

To generate the standard curve, plasmid DNAs were diluted to yield a series of concentrations, each with a 10-fold difference. The efficiencies of the qPCR assays ranged from 102% to 110%, and the R² value for each standard curve line exceeded 0.98. The C_t value (threshold cycle) was determined to quantify the copy number of all of the above mentioned genes. The qPCR amplification was conducted in 10 μL reaction mixtures, consisting of 5 μL SYBR® Premix Ex Taq™ II (Takara, Japan), 0.25 μL of each primer, 1 μL of genomic DNA and 3.5 μL dd H₂O. The primers and qPCR protocols followed Shu et al.²¹ Each qPCR reaction was run in triplicate.

2.5. Statistical analysis

Analytical measurements of NH₄⁺–N, NO₂⁻–N, NO₃⁻–N, TN, and TOC in influent and effluent wastewater in the anammox bioreactor were performed according to the standard methods.³⁰ The specific anammox activity and growth rates were measured, and the kinetic parameters were fitted using the secant method embedded in AQUASIM 2.1d.³¹

The ecological association between rates of nitrogen transformation and functional gene groups were evaluated via stepwise regression analyses (SPSS 20, USA). Furthermore, the direct and indirect contribution of different functional gene groups on nitrogen transformation rates were determined using path analysis according to the methods described by Pang et al.³² Values of direct effects (path coefficients)³³ were derived by the simultaneous solution of the normal equations for multiple linear regression using SPSS Statistics 20 (IBM, USA). Indirect effects were obtained from simple correlation coefficients between functional genes using SPSS Statistics 20 (IBM, USA) (http://www-01.ibm.com/software/analytics/spss/).³²

3. Results and discussion

3.1. Short-term experiments and kinetics evaluation

To explore the specific anammox growth activity (SAA) and anammox growth rates (μ_AN) in eight batch tests at different TOC/TN ratios, the rates of NH₄⁺–N consumption were measured and the kinetics were fitted according to the AQUASIM 2.1d.³¹ As shown in Fig. 1, the kinetics matched well with the corresponding experimental measurements. Generally speaking, with increased ratios of TOC/TN from 0.05 to 0.20, the evaluated μ_AN increased from 0.1093 d⁻¹ to 0.1236 d⁻¹ (Fig. 1a–c). The results revealed that appropriate TOC/TN ratios could improve the anammox activity. However, the μ_AN declined from 0.1183 d⁻¹ to 0.0926 d ⁻¹ when the TOC/TN ratios were increased from 0.41 to 0.82 (Fig. 1d–f). Compared with batch test 1 (0.1102 d⁻¹), the results indicated that the maximum μ_AN (μ_{AN max}) was 0.1236 d⁻¹ in the presence of a TOC/TN ratio of 0.20, which was 10.8% higher than that observed in batch test 1. In addition, the μ_AN in batch test 5 and batch test 6 were also higher than batch test 1. Therefore, these results indicated that “Ca. Brocadia sinica” have tolerance to higher TOC stress conditions (TOC/TN > 0.20).

Fig. 1 (a)–(f) The kinetic fitted and measured NH₄⁺–N consumption profiles in six 8 h batch tests under different TOC/TN ratios; (g) the actually observed and model-fitted relationships between TOC/TN ratios and specific anammox activity using substrate inhibition kinetics; (h) relationships between TOC/TN ratios and specific anammox growth rates.

Additionally, the dependence of SAA and μ_AN on the different TOC/TN stresses were well described by the substrate inhibition kinetics model. The results from Fig. 1g–h showed that 0.0524 kg N per kg VSS per d and 0.1356 d⁻¹ were the maximum SAA (SAA_max) and μ_{AN max}, respectively. The inhibition constants of nitrite and TOC were 4.8651 mmol L⁻¹ (122 mg N per L) and 2.4762 mmol L⁻¹ (196 mg L⁻¹ TOC), respectively. Meanwhile, the 95% confidence interval showed in Fig. 1g and h further indicate that the specific anammox growth rate under TOC/TN abiotic stresses could be described by eqn (1).

Kartal et al. found “Ca. Kuenenia stuttgartiensis” have the unique ability to use TOC as electron donors to reduce the nitrate and nitrite to ammonium.³⁴ Winkler et al. also found that both “Ca. Brocadia fulgida” and “Ca. Anammoxoglobus propionicus” were capable of oxidizing VFAs.¹⁰ Moreover, it was found that the appropriate influent C/N ratio for this process is 0.5 g COD per g NH₄⁺–N. Huang et al. found that “Ca. Jettenia asiatica” have an organotrophic anammox nature and that they could obtain high rates of NH₄⁺–N conversion at low COD : N ratios (<1.5).¹⁴ In addition, Du et al. found that anammox bacteria could outcompete heterotrophic denitrifiers. In this study, “Ca. Brocadia sinica” have higher rates of specific anammox growth and activity at TOC/TN ratios ranging from 0.05 to 0.61. There are two possible explanations for this result. One is that “Ca. Brocadia sinica” oxidized the TOC to CO₂ with nitrate and/or nitrite as the electron acceptor; the other is that the surplus TOC were consumed by heterotrophic microorganisms. These two possible explanations demanded exploration and either rejection or confirmation based on detailed evidence from long-term treatment performance and mass balance.

3.2. Treatment performance of organotrophic anammox and mass balance at different TOC/TN ratios

In order to evaluate the treatment profiles of the organotrophic anammox process for nitrogen and TOC removal, a long-term experiment was carried out over 140 days. As shown in Fig. 2, experiment was performed without the addition of TOC in phase I (1–62 days). During the phase I, the respective NH₄⁺–N, NO₂⁻–N, and total nitrogen removal (TN) efficiencies were 96.19 ± 2.52%, 99.85 ± 0.46%, and 89.62 ± 1.94%, respectively. Correspondingly, NH₄⁺–N, NO₂⁻–N, and TN nitrogen transformation rates reached at 0.36 ± 0.01, 0.45 ± 0.01, and 0.74 ± 0.02 kg N per m³ per d. The stoichiometric ratio of NH₄⁺–N, NO₂⁻–N, and NO₃⁻–N were 1 : (1.24 ± 0.04) : (0.20 ± 0.03), which was accordant with the theoretical values for the anammox process.¹ During phase II (63–79 days), the NH₄⁺–N removal efficiencies slightly declined to 93.28 ± 3.17% (Fig. 2). However, NO₂⁻–N and TN efficiencies increased to 100% and 92.01 ± 2.81%. The NH₄⁺–N, NO₂⁻–N, and TN transformation rates were 0.35 ± 0.01, 0.45 ± 0.02, and 0.76 ± 0.03 kg N per m³ per d, respectively (Fig. 2). Nevertheless, the stoichiometric ratios of NH₄⁺–N, NO₂⁻–N, and NO₃⁻–N were changed to 1 : (1.31 ± 0.07) : (0.12 ± 0.05). During phase III (80–95 days), the NO₂⁻–N and TN removal efficiencies did not vary significantly. However, the NH₄⁺–N removal efficiency decreased to 88.17 ± 2.78%. In addition, the NH₄⁺–N transformation rates decreased to 0.33 ± 0.01 kg N per m³ per d. During phase IV (96–106 days), the NH₄⁺–N and TN removal efficiencies dropped to 70.75 ± 10.98% and 86.85 ± 5.02%, respectively. The NH₄⁺–N transformation rates further decreased to (0.26 ± 0.04) kg N per m³ per d. During phase V (96–106 days), the effluent NH₄⁺–N levels increased remarkably and the NH₄⁺–N and TN removal efficiencies decreased sharply to 10.23 ± 6.74% and 60.23 ± 3.34%, respectively. In the recovery phase (120–140 days), the operational strategy without TOC was performed. Because the anammox was seriously suppressed under high TOC/TN stress, the influent NH₄⁺–N and NO₂⁻–N levels were set at 70 and 84 mg L⁻¹, respectively. Meanwhile, 10 mL high-loaded anammox sludge from another anammox reactor was added at the beginning of day 125.³⁵ The results from Fig. 2a indicate that this type of high-loaded sludge could be beneficial for the rapid recovery of anammox activity.

Fig. 2 Long-term performance of anammox–SBR system under different TOC/TN ratio stresses. (a) Nitrogen concentration; (b) nitrogen removal efficiency; (c) nitrogen transformation rates; (d) nitrogen load.

In general, the results from the long-term treatment performance of assessment of the organotrophic anammox system under different TOC/TN ratio conditions indicated that the autotrophic anammox process was the primary process accounting for nitrogen removal in phases I–III, which was consistent with the results from the short-term experiment in this study as well as the work of Tang et al.⁹ In addition, due to the weakening of the anammox process in phase IV, the decreased production of NO₃⁻–N indicated that organotrophic anammox and heterotrophic denitrifiers could coexist for purposes of nitrogen removal. There are two reasons which may explain these results: (1) organotrophic anammox bacteria could consume the acetate and propionate to reduce the nitrate and nitrite to ammonia;^11,34 (2) heterotrophic denitrification and denitritation could contribute to the TN removal.

The severely inhibited anammox in phase V along with increased NH₄⁺–N concentration and further decreased NO₃⁻–N effluent production indicate that autotrophic anammox could not outcompete heterotrophic microorganisms under high TOC/TN ratios. Moreover, it was found that lower ratios of TOC/TN (<0.2) were useful for simultaneous removal of nitrogen and TOC but anammox growth rates were suppressed by higher TOC/TN ratios (>0.40).

Based on the above results, it could be seen that the average NO₃⁻–N production was 35.37 ± 6.18 mg L⁻¹ without the addition of acetate and propionate in phase I, which was very close to the theoretical production. However, during phases II–V, the average nitrate production was far less than that in phase I. For instance, the average nitrate production was 20.32 mg L⁻¹ in phase II. Ahn et al.³⁶ have reported that 1 g nitrate would be reduced when 3.31 g chemical oxygen demand (COD) was consumed. Thus, 15.05 mg (35.37 − 20.32 = 15.05) nitrate would be reduced in another pathway when 49.82 mg COD (15.05 × 3.31 = 49.82) was consumed. In phase II, the TOC/TN ratio was 0.19 (equal to 110 mg L⁻¹ COD) and the TOC removal efficiency was 36.65 ± 4.56%, revealing that 40.32 mg L⁻¹ was consumed in phase II. It is obvious that nitrate was in too short a supply to foster denitrification in phases III–V. It can be concluded that organotrophic anammox bacteria which consumed TOC and thus reduced the levels of nitrate existed in phases II–V.³⁷

Some studies have reported a total of four possible nitrogen and TOC consumption processes: the autotrophic anammox process eqn (2); the organotrophic anammox process eqn (3); the denitritation process eqn (4); and the denitrification process eqn (5).^5,12,38 To further understanding the nitrogen pathway in the organotrophic anammox process, the removal of NH₄⁺–N, NO₂⁻–N, NO₃⁻–N, and TOC were calculated based on nitrogen, TOC mass balance and eqn (2)–(5).

Autotrophic anammox:

NH₄⁺ + 1.32NO₂⁻ + 0.066HCO₃⁻ + 0.13H⁺ → 1.02N₂ + 0.26NO₃⁻ + 0.066CH₂O_0.5N_0.15 + 2.03H₂O (2)

Organotrophic anammox:

NH₄⁺ + NO₃⁻ + HCOO⁻ + H⁺ → N₂ + CO₂ + 3H₂O, (ΔG′_o = −522.4 kJ mol⁻¹) (3)

Denitritation:

NO₂⁻ + 0.19CH₃CH₃CH₂COOH + H₂CO₃ → 0.037C₅H₇NO₂ + 0.481N₂ + 1.14H₂O + HCO₃⁻ + CO₂ (4)

Denitrification:

NO₃⁻ + 0.29CH₃CH₃CH₂COOH + H₂CO₃ → 0.034C₅H₇NO₂ + 0.483N₂ + 1.54H₂O + HCO₃⁻ + 0.986CO₂ (5)

As shown in Table 2, during phases II–V of the organotrophic anammox, only a small contribution to the overall nitrogen removal was made. In general, as the TOC/TN ratio increased from 0.10 to 0.80, the average nitrogen removal pathway percentage observed in the autotrophic anammox process decreased from 76.03% (phase II) to 2.42% (phase V). Correspondingly, the heterotrophic denitritation contribution increased from 7.05% (phase II) to 80.81% (phase V). These results indicate that autotrophic anammox and organotrophic anammox could not outcompete heterotrophic denitrification at high TOC/TN ratios. This was because the growth rate of the heterotrophic denitrification process enjoys have much shorter doubling times (2–16 hours), which renders the overall growth rate higher than that of anammox bacteria (8–12 days).¹²

Table 2 Contribution of different pathways to nitrogen removal on the basis of COD (Chemical oxygen demand) and nitrogen mass balance at different TOC/TN ratios

Consumption	Removal route	Phase II	Phase III	Phase IV	Phase V
NH4^+–N removal (mg L^−1)	Anammox (autotrophic)	150.19	135.96	116.91	7.84
	Anammox (organotrophic)	22.87	30.04	15.53	11.46
NO2^−–N removal (mg L^−1)	Anammox (autotrophic)	198.25	179.47	154.32	10.35
	Denitritation	27.73	46.51	71.66	215.63
NO3^−–N removal (mg L^−1)	Anammox (organotrophic)	22.87	30.04	15.53	11.46
	Denitrification	—	—	—	—
Biomass increased (mg N)	Anammox (autotrophic)	9.91	8.97	7.72	0.52
	Anammox (organotrophic)	0.75	0.99	0.51	0.38
	Denitritation	1.03	1.72	2.65	7.98
Total biomass increased (mg N)	—	11.69	11.69	10.88	8.87
Total nitrogen removal (mg L^−1)	—	378.72	400.27	367.5	256.95
Average percentage of nitrogen removal routes (%)	Anammox (autotrophic)	76.03	69.06	67.48	2.42
	Anammox (organotrophic)	11.88	14.76	8.31	8.77
	Denitritation	7.05	11.19	18.78	80.81
	Denitrification	—	—	—	—
	Other pathways	5.04	4.99	5.43	7.99

---

Although the TOC/TN ratios were further increased, the organotrophic anammox contribution rates of 11.88%, 14.76%, 8.31%, and 8.77%, respectively, did not vary significantly. Interestingly, the maximum organotrophic anammox contribution was achieved at the TOC/TN ratio of 0.20. These results indicated that the organotrophic anammox activity had been maximized in this TOC/TN ratio. Meanwhile, the organotrophic anammox could outcompete heterotrophic denitrification in phases II–III. The reason for these favorable results is that organotrophic anammox could utilize nitrate and nitrite as electron acceptors to oxidized acetate and propionate under these conditions. Other nitrogen pathways such as fermentation and sulfur-based autotrophic denitrification varied from 5.04% to 7.99%.⁷

3.3. Quantification of nitrogen-related 16S rRNA and functional genes

In order to gain insight into the variation of nitrogen-related functional genes under TOC/TN stress conditions, samples of anammox biomass were taken from the end of each phase and the copy numbers of all above mentioned 16S rRNA and functional genes were quantified. As shown in Fig. 3a–f, the gene copy numbers of bacterial 16S rRNAs were in the same order of magnitude in phases I–V. The gene copy numbers of anammox 16S rRNAs were also in the same order of magnitude ranging from 2.50 × 10⁸ to 4.07 × 10⁸ copies per g wet sludge even though the anammox were severely suppressed, indicating that the abundance of anammox had no remarkable reduction in the system, which was consistent with the finding of Du et al.²³ It was evident that higher TOC/TN ratios (>0.40) could result in the community succession of anammox species in this study.

Fig. 3 Absolute abundance of microbial communities and functional genes in the anammox–SBR system. (a) Total bacterial and anammox bacterial 16S rRNA; (b) AOA amoA, AOB amoA, and nxrA; (c) napA, narG, and nrfA; (d) nirK, nirS, and nosZ; (e) mcrA and dsrA; (f) ratios of different functional gene groups. Error bars represent standard deviation calculated from three independent experiments.

The copy numbers of three nitrification groups, including AOA amoA, AOB amoA and nxrA genes, were summarized in Fig. 3b. With the increase in TOC/TN ratios from 0–0.40 during phases I–IV, the absolute abundance of AOB amoA and nxrA gradually increased from 1.66 × 10⁶ to 1.18 × 10⁷ copies per g wet sludge, 4.87 × 10⁴ to 1.04 × 10⁵ copies per g wet sludge, respectively. However, higher TOC/TN ratios in phase V had adverse effects on these two genes, suggesting that the activity of AOB and NOB (nitrite oxidizing bacteria) have tolerate to TOC/TN ratios of 0.1–0.4 well. In addition, the absolute abundance of AOA amoA increased from 1.11 × 10² to 1.21 × 10² copies per g wet sludge in phases III–V, suggesting that AOA had positive relationship with TOC. The results indicated that AOA might have an organotrophic nature, which was consistent with the work of Mußmann et al.³⁹

As illustrated in Fig. 3c, the dissimilatory nitrogen reduction gene nrfA increased from 2.02 × 10⁴ to 6.95 × 10⁵ copies per g wet sludge during phases I–V as the TOC/TN ratios increased. Kartal et al. found that some anammox bacteria species could engage in “disguised denitrification” in the form of the reduction of nitrate and nitrite (termed DNRA).³⁴ As depicted in Fig. 3f, the variation of nrfA/anammox and anammox/bacteria had a high degree of consistency, indicating that DNRA or organotrophic anammox had greatly contributed to nitrogen removal in the organotrophic anammox process. In addition, Fig. 3c indicated that the absolute gene abundance of napA increased in phases I–IV, while it decreased in phase V. As compared to the napA gene, the absolute gene abundance of the narG gene in phases I–III was nearly one to two orders of magnitude higher than in phases IV–V. The conversion of NO₃⁻ → NO₂⁻ is catalyzed by the napA and narG genes. Therefore, the results indicated that high TOC/TN ratios in phase IV–V had no significant influence on the rates of NO₃⁻–N reduction.

As shown in Fig. 3d, the absolute abundance of nirS in phase IV was nearly one order of magnitude greater than other four phases. The absolute abundance of the nirK gene in phases I–IV was nearly the same order of magnitude, while it slightly decreased in phase V. The absolute abundance of the nosZ gene increased from 2.28 × 10⁴ to 6.59 × 10⁵ copies per g wet sludge during phases I–V as the TOC/TN ratios increased. These results showed that the absolute abundance of denitrification genes did not vary dramatically, suggesting that denitrification was not dominant process accounting for nitrogen removal under the different TOC/TN stress conditions; this was accordant with the results of Table 2.

Furthermore, as shown in Fig. 3e, the absolute abundance of the mcrA gene increased gradually during phases I–IV, while it decreased in phase V. In terms of the dsrA gene, copy numbers increased from 4.14 × 10³ to 2.98 × 10⁵ copies per g wet sludge. Notably, as illustrated in phases IV–V, methanogens could compete with sulfate reduction bacteria for the TOC.

Taken together, it is plausible that the co-existence of autotrophic anammox, organotrophic anammox (or DNRA), and denitritation were the primary pathway that accounted for nitrogen removal under TOC/TN stress conditions.

3.4. Functional gene groups determining nitrogen transformation rates

In this study, single functional genes and functional gene group ratios were used as input variables and introduced into stepwise regression analyses. Results showed nitrogen transformation in phases I–II differed from that in phases III–V. Therefore, the key functional groups in phases I–II and in phases III–V were separately investigated to characterize the nitrogen transformation pathway under different TOC/TN stress conditions.

As shown in Table 3, (AOA amoA + AOB amoA + anammox + nrfA)/bacteria, nrfA/(narG + napA), and nrfA were the key factors for NH₄⁺–N transformation in phases I–V. The variable, (AOA amoA + AOB amoA + anammox + nrfA)/bacteria showed a positive relationship with NH₄⁺–N transformation in phases I–II. The AOA amoA and AOB amoA genes are often regarded as NH₄⁺–N to NO₂⁻–N oxidation markers, and nrfA is often regarded as a NO₃⁻–N to NH₄⁺–N reduction makers.²⁴ Thus, results showed that the coupling of partial nitrification, anammox, and DNRA was involved in NH₄⁺–N conversion. The nrfA/(narG + napA) group was the key functional group for NO₃⁻–N consumption in phases III–V. Both narG and napA were regarded as marker genes for NO₃⁻–N to NO₂⁻–N in the first denitrification step. Hence, nrfA/(narG + napA) and nrfA denoted NO₃⁻–N reduction, showing a positive relationship and a negative relationship with NH₄⁺–N conversion, respectively. This result suggested that dissimilatory nitrate reduction to ammonium affected NH₄⁺–N conversion in phases III–V.

Table 3 Quantitative response relationships between nitrogen transformation rates (mg L⁻¹ d⁻¹) and functional genes abundance (copies per g sludge) in long-term experiments. Anammox represent the absolute abundance of anammox bacterial 16S rRNA. Bacteria represent the absolute abundance of total bacterial 16S rRNA

Stepwise regression equations	R^2	p value
(TOC/TN = 0 and TOC/TN = 0.10)
NH4^+–N = 385.638(AOA amoA + AOB amoA + anammox + nrfA)/bacteria + 313.386	1.000	0.013
NO2^−–N = 73.546anammox/bacteria + 458.162	0.995	0.024
NO3^−–N = −8679.981nosZ/nirS + nirK + 761.040	0.983	0.008
(TOC/TN = 0.20, TOC/TN = 0.40, TOC/TN = 0.80)
NH4^+–N = −2.709nrfA/(narG + napA) − 0.002nrfA + 1143.337	0.997	0.008
NO2^−–N = 8.067 × 10^−6nirK − 246.61AOB amoA/(nxrA + anammox + nirK + nirS) + 475.431	1.000	0.017
NO3^−–N = 0.103(nirK + nirS)/AOB amoA − 0.002nrfA + 1578.045	0.982	0.005

---

The NO₂⁻–N transformation rates were jointly determined by anammox/bacteria, nirK, and AOB amoA/(nxrA + anammox + nirK + nirS) in phases I–V. The variable anammox/bacteria in the NO₂⁻–N equation which was denoted NO₂⁻–N consumption showed a positive relationship with NO₂⁻–N transformation in phases I–II. The results suggested that anammox were the primary NO₂⁻–N removal pathway under low TOC/TN ratios (≤0.10). The variable nirK genes were regarded as marker genes for NO₂⁻–N to NO in the denitritation step, showing a positive correlation with NO₂⁻–N conversion. The variable AOB amoA/(nxrA + anammox + nirK + nirS) group which was denoted NO₂⁻–N accumulation showed a negative relationship with NO₂⁻–N conversion. This result indicated that the denitritation step was the primary factor which accounted for the NO₂⁻–N transformation at high TOC/TN ratios (≥0.20), while nitrification was inhibited in phases III–V.

The NO₃⁻–N transformation rates were collectively determined by nosZ/(nirS + nirK), (nirK + nirS)/AOB amoA, and nrfA in phases I–V. nirS and nirK genes are often regarded as NO₂⁻–N to NO genes. The nosZ(/nirS + nirK) in phases I–II was negatively associated with NO₃⁻–N transformation rates, suggesting that denitrification was not significant under low TOC/TN ratios (≤0.10). The (nirK + nirS)/AOB amoA in phases III–V was denoted NO₂⁻–N consumption showed a positive relationship with NO₃⁻–N conversion. The more NO₂⁻–N that was consumed, the more NO₃⁻–N was transformed; this was because lower NO₂⁻–N concentrations can reduce its toxic effects on denitrifiers.³²

Based on these results, the simultaneous presence of anammox, denitrification, and organotrophic anammox (DNRA) were confirmed in the anammox–SBR system under different TOC/TN constrains. This coupling of anammox, denitrification and DNRA can assist in the simultaneous removal of nitrogen and organic carbon in a single system, rather than over a sequential chain of treatments.⁴⁰

3.5. Molecular mechanism of nitrogen transformation pathways

Path analysis was further performed to explore the contribution of functional gene groups to the process of nitrogen removal in phases I–V. Path analysis revealed that (AOA amoA + AOB amoA + anammox + nrfA)/bacteria has a direct positive effect on the NH₄⁺–N transformation rate (0.216) under TOC/TN ratios ≤ 0.10 (Fig. 4 a). However, nrfA/(narG + napA) and nrfA also has a lower direct negative effect on the NH₄⁺–N transformation rate (−1.101 and −0.143, respectively) under TOC/TN ≥ 0.20. The indirect positive effect of (AOA amoA + AOB amoA + anammox + nrfA)/bacteria via nrfA on the NH₄⁺–N transformation rate was 0.107. The indirect positive effect of nrfA via (AOA amoA + AOB amoA + anammox + nrfA)/bacteria on the NH₄⁺–N transformation rate was 0.058, which was lower than the aforementioned indirect effect. These results suggested that (AOA amoA + AOB amoA + anammox + nrfA)/bacteria was the better predictive variable and a key functional group for the NH₄⁺–N transformation rate. The results also indicated that coupling nitrification and anammox was the primary factor behind NH₄⁺–N removal under TOC/TN ratios ≤ 0.10, while DNRA was the pathway contributing to NH₄⁺–N accumulation under TOC/TN ≥ 0.20.

Fig. 4 Path diagrams estimating the effects of the functional gene groups on (a) NH₄⁺–N transformation rate, (b) NO₂⁻–N transformation rate, and (c) NO₃⁻–N transformation rate under different TOC/TN ratios constrains.

The direct effect of anammox/bacteria and nirK on NO₂⁻–N transformation rates were positive and values were 0.590 and 0.841, respectively. The direct effect of AOB amoA/(nxrA + anammox + nirK + nirS) on the NO₂⁻–N transformation rate was negative (−0.869) and was lower in strength than the positive direct effects. The indirect effect of anammox/bacteria via AOB amoA/(nxrA + anammox + nirK + nirS) on the NO₂⁻–N transformation rate was positive (0.172), and the indirect effect of nirK via AOB amoA/(nxrA + anammox + nirK + nirS) on the NO₂⁻–N transformation rate was also positive (0.108) (Fig. 4b). However, the indirect effect of AOB amoA/(nxrA + anammox + nirK + nirS) via anammox/bacteria and nirK on the NO₂⁻–N transformation rates were negative. The results showed that anammox/bacteria was the key factor for NO₂⁻–N removal under TOC/TN ratios ≤ 0.10, and nirK was the better predictive variable for NO₂⁻–N removal under TOC/TN ratios ≥ 0.20. The results also indicated that anammox was the primary process responsible for NO₂⁻–N removal in phases I–II, while denitritation was the primary pathway responsible for NO₂⁻–N conversion in phases III–V, which was congruent with the discussion in Section 3.4.

It was found that nosZ/(nirS + nirK) and nrfA had direct negative effects on the NO₃⁻–N transformation rates (−0.758 and −1.012, respectively) in phases I–V. However, (nirK + nirS)/AOB amoA had positive effects on the NO₃⁻–N transformation rate (0.149) (Fig. 4c). The indirect positive effects of (nirK + nirS)/AOB amoA and nrfA via nosZ/(nirK + nirS) were 0.032 and 0.133, respectively. The indirect negative effects of (nirK + nirS)/AOB amoA and nosZ/(nirK + nirS) via nrfA were −0.234 and −0.117, respectively. The results indicated that (nirK + nirS)/AOB amoA and nrfA were the best predictive variable and a major contributing factor to the determination of the NO₃⁻–N transformation rate. The results also suggested that denitrification was not primary factors in phases I–V, while nrfA and denitritation was the primary process responsible for NO₃⁻–N and NO₂⁻–N removal under at TOC/TN ratios ≥ 0.20, supporting the results presented in Section 3.4 and in Table 3.

Compared with functional gene copy numbers, the results have evidenced that key functional gene groups can serve as integrative variables to characterize nitrogen transformation rates in the organotrophic anammox process.⁴¹ On the quantitative molecular level, these analyses clearly revealed that the co-existence of autotrophic anammox, denitrification, and DNRA (organotrophic anammox) could useful for simultaneous nitrogen and TOC removal within the organotrophic niche.

3.6. Dynamics of bacterial communities and functional generalists

In this study, MiSeq sequencing was applied to investigate the dynamics of microbial communities and functional generalists in the organotrophic anammox process. After the filtration of raw sequences, 14 505–21 586 effective reads were obtained from five samples. As shown in Table S2,† OTUs were in the range of 1112–1750, and increased in phases II–V as the TOC/TN ratio increased. Two estimators, which were Good's coverage and Simpson, showed unremarkable variation. However, Shannon, Chao 1, and ACE estimators varied significantly. These results indicated that TOC/TN stresses did not significantly effect the richness of the bacteria in the system, while they could significantly promote the diversity of related bacteria in the organotrophic anammox process. The rarefaction curves of five samples were also drawn in Fig. S1.† Results showed that the rarefaction curves did not reach a plateau in all phases, indicating that rare species could continue to emerge when the sequence depth exceed 14 000. As displayed in Fig. S2,† five samples belonged to four groups, indicating that lower TOC/TN ratios (<0.2) did not significantly influence microbial diversity.

In this study, effective sequences were assigned to phyla, class, order, family, and genera using RDP classifier via Silva SSU database. As shown in Fig. 5, a total of 15 bacterial phyla were identified. Planctomycetes was the most dominant phylum in all phases, accounting for 23.40–37.83% (averaging 31.25%). The other dominant phyla were Proteobacteria (11.45–36.92%, averaging 22.22%), Chloroflexi (1.95–28.91%, averaging 17.93%), and Chlorobi (3.54–14.01%, averaging 11.78%). Previous studies have reported that Proteobacteria, Chloroflexi and Planctomycetes were the significant phyla in the nitritation–anammox system.⁴² Fig. 5 clearly showed that Planctomycetes was more abundant than Proteobacteria, Chloroflexi, and Chlorobi, suggesting that Planctomycetes played the dominant role for nitrogen removal. In addition, as TOC/TN ratios increased from 0 to 0.20, the percentage of Planctomycetes increased from 28.48% to 37.83%. Then, it decreased to 29.59% as the TOC/TN ratio increased to 1.5%. The results indicated that Planctomycetes had higher tolerance for TOC/TN stresses, and that the appropriate addition of acetate and propionate was more favorable to enriching organotrophic anammox bacteria species in the phyla of Planctomycetes; these findings were consistent with the results reported in Section 3.1.

Fig. 5 Distribution of phyla in the different phase based on the taxonomy annotation from SILVA SSU database using RDP classifier. The thickness of each ribbon represents the abundance of each taxon. The absolute tick above the inner segment and relative tick above the outer segment stand for the reads abundances and relative abundance of each taxon, respectively. Others refer to those unclassified reads. The data were visualized using Circos (Version 0.67, http://circos.ca/).

Notably, among Proteobacteria, β-Proteobacteria was the most dominant in all phases, followed by α-Proteobacteria, γ-Proteobacteria, δ-Proteobacteria, and ε-Proteobacteria (Fig. S3a†). Beside the class, the results from Fig. S3b and c† showed that the following major orders, namely Ignavibacteriales, Clostridiales, Phycisphaerales, Brocadiales, Rhodobacterales, and Rhodocyclales, and their corresponding families were shared by five phases, suggesting that these dominant populations played pivotal roles for nitrogen removal in the system.

As shown in Fig. S3d,† a total of 144 genera were assigned and 27 of them were dominant genera. Meanwhile, these genera were identified as belonging to 8 functional groups (Fig. 6). Among them, the anammox group involved “Ca. Brocadia sinica”, “Ca. Jettenia asiatica”, and “Ca. Kuenenia stuttgartiensis”. The relative abundance of “Ca. Brocadia sinica” were 4.32%, 4.43%, 4.73%, 5.3%, and 3.44%, respectively. Compared with “Ca. Brocadia sinica”, the relative abundance of “Ca. Jettenia asiatica” and “Ca. Kuenenia stuttgartiensis” were not detected in the phase I. These results indicated that “Ca. Brocadia sinica” was the dominant anammox bacteria species in the autotrophic system. With the addition of volatile fatty acids, the quantity of “Ca. Brocadia sinica” decreased. This was due to the ability of “Ca. Brocadia sinica” to oxidize volatile fatty acids.³⁴ During phases II–III, the percentage of “Ca. Jettenia asiatica” increased from 4.02% to 4.12%, while it decreased significantly in phase V. “Ca. Kuenenia stuttgartiensis” was detected at the lower TOC/TN ratio and decreased remarkably in phases IV and V. However, the quantity of “Ca. Kuenenia stuttgartiensis” increased from 3.01% to 3.67% with the further increased in the TOC/TN ratio from 0.4 to 0.8. This was due to the K-strategy survival ability of “Ca. Kuenenia stuttgartiensis”,⁴³ which prefers this strategy to low TOC conditions. However, “Ca. Kuenenia stuttgartiensis” was able to adapt to a higher TOC niche during the long-term acclimation. The appearance of “Ca. Jettenia asiatica” and “Ca. Kuenenia stuttgartiensis” in phases II–V resulted in a lack of significant variation in the anammox copy numbers despite suppression of the nitrogen removal activity in phase V.

Fig. 6 The relative abundance of total 9 functional genera in five samples.

Based on the results from the nitrogen treatment performance and MiSeq sequencing, it can be concluded that organotrophic anammox bacteria species have the capacity to oxidize acetate and propionate. In addition, “Ca. Brocadia sinica”, “Ca. Jettenia asiatica”, and “Ca. Kuenenia stuttgartiensis” have an organotrophic nature under the appropriate TOC/TN stress conditions. Furthermore, it was evident that higher TOC/TN ratios (>0.40) could result in the community succession of anammox species and alter the character of the microbial communities observed in this study.

4. Conclusion

Short- and long-term experiments revealed that the appropriate TOC/TN ratios enabled the maximum growth rate of the anammox to increase to 0.1356 d⁻¹. TOC biomass balance revealed that organotrophic anammox accounted for 14.76% of the nitrogen loss at a TOC/TN ratio of 0.20. Quantitative molecular analyses and pathway results confirmed that the coupling of anammox, DNRA, and denitrification was pivotal to the nitrogen transformation pathway in the organotrophic anammox process. MiSeq sequencing indicated that Planctomycetes, Proteobacteria, Chloroflexi, and Chlorobi were the most abundant phyla in the system. Furthermore, “Ca. Brocadia sinica” had tolerance to higher TOC stress conditions, and single environmental factors such as TOC/TN ratios could lead to microbial succession between “Ca. Brocadia sinica”, “Ca. Jettenia asiatica”, and “Ca. Kuenenia stuttgartiensis”. However, the quantitative molecular mechanism of the organotrophic anammox process in the anammox–SBR system needs further study using ¹³C-labeled DNA/RNA-SIP techniques. Moreover, the molecular mechanism for microbial succession between “Ca. Kuenenia stuttgartiensis”, “Ca. Jettenia asiatica”, and “Ca. Kuenenia stuttgartiensis” should also be further explored using metagenomic and metatranscriptomic approaches.

Acknowledgements

This study was financially supported by the National Natural Science Foundation of China (51308453) and the Science & Technology Innovation Program of Shaanxi Province (2011KTZB03-03-01).

References

1. A. A. Van de Graaf, P. de Bruijn, L. A. Robertson, M. S. Jetten and J. G. Kuenen, Microbiology, 1996, 142, 2187–2196.
2. M. Strous, J. G. Kuenen and M. S. Jetten, Appl. Environ. Microbiol., 1999, 65, 3248–3250.
3. J. Guo, Y. Peng, L. Fan, L. Zhang, B. J. Ni, B. Kartal, X. Feng, M. S. Jetten and Z. Yuan, Environ. Microbiol., 2015, 16, 1–15.
4. T. Lotti, R. Kleerebezem, C. van Erp Taalman Kip, T. L. Hendrickx, J. Kruit, M. Hoekstra and M. C. van Loosdrecht, Environ. Sci. Technol., 2014, 48, 7874–7880.
5. B. Kartal, N. M. Almeida, W. J. Maalcke, H. J. Camp, M. S. Jetten and J. T. Keltjens, FEMS Microbiol. Rev., 2013, 1–34.
6. S. Jenni, S. E. Vlaeminck, E. Morgenroth and K. M. Udert, Water Res., 2014, 49, 316–326.
7. C. Chen, X. Huang, C. Lei, T. C. Zhang and W. Wu, Bioresour. Technol., 2013, 148, 172–179.
8. S.-Q. Ni, J.-Y. Ni, D.-L. Hu and S. Sung, Bioresour. Technol., 2012, 110, 701–705.
9. C.-j. Tang, P. Zheng, C.-h. Wang and Q. Mahmood, Bioresour. Technol., 2010, 101, 1762–1768.
10. M. K. H. Winkler, J. Yang, R. Kleerebezem, E. Plaza, J. Trela, B. Hultman and M. C. M. van Loosdrecht, Bioresour. Technol., 2012, 114, 217–223.
11. B. Kartal, J. Rattray, L. A. van Niftrik, J. van de Vossenberg, M. C. Schmid, R. I. Webb, S. Schouten, J. A. Fuerst, J. S. Damsté and M. S. Jetten, Syst. Appl. Microbiol., 2007, 30, 39–49.
12. J. G. Kuenen, Nat. Rev. Microbiol., 2008, 6, 320–326.
13. B. Kartal, L. van Niftrik, J. T. Keltjens, H. J. Op den Camp and M. S. Jetten, Adv. Microb. Physiol., 2012, 60, 212.
14. X.-L. Huang, D.-W. Gao, Y. Tao and X.-L. Wang, Chem. Eng. J., 2014, 253, 402–407.
15. B. Kartal, L. Van Niftrik, J. Rattray, J. L. Van De Vossenberg, M. C. Schmid, J. S. Damsté, M. S. Jetten and M. Strous, FEMS Microbiol. Ecol., 2008, 63, 46–55.
16. H. Kamyab, M. F. M. Din, S. K. Ghoshal, C. T. Lee, A. Keyvanfar, A. A. Bavafa, S. Rezania and J. S. Lim, Waste Biomass Valorization, 2016, 1–10.
17. H. Kamyab, M. F. M. Din, S. E. Hosseini, S. K. Ghoshal, V. Ashokkumar, A. Keyvanfar, A. Shafaghat, C. T. Lee, A. Asghar Bavafa and M. Z. A. Majid, Clean Technol. Environ. Policy, 2016, 1–11.
18. H. Kamyab, M. F. M. Din, A. Keyvanfar, M. Z. A. Majid, A. Talaiekhozani, A. Shafaghat, C. T. Lee, L. J. Shiun and H. H. Ismail, Energy Procedia, 2015, 75, 2400–2408.
19. T. Zhang, M.-F. Shao and L. Ye, ISME J., 2011, 6, 1137–1147.
20. E. Isanta, T. Bezerra, I. Fernández, M. E. Suárez-Ojeda, J. Pérez and J. Carrera, Bioresour. Technol., 2015, 181, 207–213.
21. D. Shu, Y. He, H. Yue, L. Zhu and Q. Wang, Bioresour. Technol., 2015, 196, 621–633.
22. T. Lotti, R. Kleerebezem, Z. Hu, B. Kartal, M. Jetten and M. van Loosdrecht, Water Res., 2014, 66, 111–121.
23. R. Du, Y. Peng, S. Cao, C. Wu, D. Weng, S. Wang and J. He, Bioresour. Technol., 2014, 162, 316–322.
24. W. Zhi, L. Yuan, G. Ji and C. He, Environ. Sci. Technol., 2015, 49, 4575–4583.
25. C.-J. Tang, P. Zheng, L.-Y. Chai and X.-B. Min, Chem. Eng. J., 2013, 230, 149–157.
26. D. Shu, Y. He, H. Yue and Q. Wang, Chem. Eng. J., 2016, 290, 21–30.
27. C. K. Lee, C. W. Herbold, S. W. Polson, K. E. Wommack, S. J. Williamson, I. R. McDonald and S. C. Cary, PLoS One, 2012, 7, e44244.
28. P. D. Schloss, S. L. Westcott, T. Ryabin, J. R. Hall, M. Hartmann, E. B. Hollister, R. A. Lesniewski, B. B. Oakley, D. H. Parks and C. J. Robinson, Appl. Environ. Microbiol., 2009, 75, 7537–7541.
29. S. J. Blazewicz, R. L. Barnard, R. A. Daly and M. K. Firestone, ISME J., 2013, 7, 2061–2068.
30. E. W. Rice, L. Bridgewater and A. P. H. Association, Standard methods for the examination of water and wastewater, American Public Health Association, Washington, DC, 2012.
31. P. Reichert, Aquasim 2.0-user manual, computer program for the identification and simulation of aquatic systems, Swiss Federal Institute for Environmental Science and Technology (EAWAG), 1998, p. 219.
32. Y. Pang, Y. Zhang, X. Yan and G. Ji, Environ. Sci. Technol., 2015, 49, 13550–13557.
33. D. F. Alwin and R. M. Hauser, Am. Sociol. Rev., 1975, 37–47.
34. B. Kartal, M. M. Kuypers, G. Lavik, J. Schalk, H. J. Op den Camp, M. S. Jetten and M. Strous, Environ. Microbiol., 2007, 9, 635–642.
35. C.-J. Tang, P. Zheng, T.-T. Chen, J.-Q. Zhang, Q. Mahmood, S. Ding, X.-G. Chen, J.-W. Chen and D.-T. Wu, Water Res., 2011, 45, 201–210.
36. Y. Ahn, I. Hwang and K. Min, Water Sci. Technol., 2004, 49, 145–153.
37. Y. Liang, D. Li, X. Zhang, H. Zeng, Y. Yang and J. Zhang, Bioresour. Technol., 2015, 193, 408–414.
38. Y.-H. Ahn, Process Biochem., 2006, 41, 1709–1721.
39. M. Mußmann, I. Brito, A. Pitcher, J. S. S. Damsté, R. Hatzenpichler, A. Richter, J. L. Nielsen, P. H. Nielsen, A. Müller and H. Daims, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 16771–16776.
40. G. D. Song, S. M. Liu, H. Marchant, M. M. M. Kuypers and G. Lavik, Biogeosciences, 2013, 10, 6851–6864.
41. G. Ji, C. He and Y. Tan, Ecol. Eng., 2013, 55, 35–42.
42. Z.-r. Chu, K. Wang, X.-k. Li, M.-t. Zhu, L. Yang and J. Zhang, Chem. Eng. J., 2015, 262, 41–48.
43. W. R. Van der Star, W. R. Abma, D. Blommers, J.-W. Mulder, T. Tokutomi, M. Strous, C. Picioreanu and M. Van Loosdrecht, Water Res., 2007, 41, 4149–4163.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c6ra04114k

‡ These authors contributed equally to this work.

---