Metabolic network analysis reveals microbial community interactions in anammox granules

1Microbial communities mediating anaerobic ammonium oxidation (anammox) represent one of the most energy-efficient environmental biotechnologies for nitrogen removal from wastewater. The process is typically used to treat high strength ammonium wastewaters and offers significant cost savings compared to conventional nitrogen removal processes that require energy-intensive aeration for nitrification and also consume large quantities of organic carbon during denitrification. In practice, anammox-based wastewater treatment systems are combined with a nitritation step in either a single-stage or two-stage bioreactor system. In these engineered ecosystems, a fraction of the ammonium is first oxidized to nitrite by aerobic ammonia oxidizing bacteria. Subsequently, anammox bacteria anaerobically oxidize the remaining ammonium directly to nitrogen gas using the produced nitrite as a terminal electron acceptor.

2Five genera of anammox bacteria have been discovered to date, including Kuenenia, Brocadia, Anammoxoglobus and Jettenia commonly found in activated sludge, and Scalindua commonly found in marine environments. These lineages all have ‘Candidatus' status as they do not exist in pure culture and must be grown in laboratory enrichments. Because anammox bacteria have a slow growth rate, either biofilm reactors that use carrier media or granular sludge reactors are used to retain sufficient biomass in the system. Such reactor configurations support the formation of dense microbial communities that can be readily separated from the liquid wastewater and enriched in the bioreactor. Anammox granules consist of a mixture of cell aggregates and abiotic particles embedded within a matrix of organic extracellular polymeric substances (EPS). EPS present in anammox granules has been found to contain high amounts of protein and polysaccharides, where increased hydrophobic amino acid content has been observed to be a main factor determining granule aggregation ability. The high aggregate density also serves to limit oxygen diffusion into the granule interior, allowing for the proliferation of anammox and other anaerobic bacteria.

3While most studies have focused on understanding the ecophysiology of anammox bacteria and their interactions with autotrophic nitrifying bacteria, little is known about the activity of heterotrophic bacteria in anammox bioreactors. Previous studies based on 16S ribosomal RNA (rRNA) gene clone libraries and amplicon sequencing have shown that heterotrophic bacteria affiliated with the phyla Chlorobi, Bacteroidetes, Chloroflexi and Proteobacteria comprise a large fraction of the microbial community in anammox bioreactors. Despite differences in bioreactor influent composition, these heterotrophic bacteria appear to share high phylogenetic similarity across different anammox wastewater treatment systems, suggesting that a more universal interaction exists between them and anammox bacteria.

4The exact role of heterotrophic bacteria in anammox systems has not yet been determined, though a few clues have been uncovered. A recent metagenomic study revealed that most of the heterotrophic organisms in anammox granules encode the ability to respire nitrate via partial denitrification, possibly completing a nitrite loop with anammox and nitrite oxidizing bacteria (NOB) by reducing nitrate back to nitrite. This activity could contribute to the removal of excess nitrate produced from the system during anammox growth or nitrite oxidation by NOB. Hydrolysis of EPS into soluble compounds and/or the secretion of soluble microbial products by anammox bacteria during cell growth is believed to support denitrification and heterotrophic growth, particularly for anammox bioreactors that receive no external organic carbon substrates. However, the specific metabolite exchange reactions promoting interactions between heterotrophic and anammox bacteria remain poorly understood.

5Here, we combine metagenomic and metatranscriptomic data to examine the gene expression of anammox and heterotrophic bacteria in a laboratory-scale anammox bioreactor and to identify their potential interactions. Metagenomic binning was used to recover near-complete population genomes from members of the microbial community inhabiting anammox granules. Resulting population genomes were used to reconstruct each organism's metabolic network and served as reference platforms for the profiling of gene expression at the community scale. We find that Chlorobi-affiliated bacteria may be highly active protein degraders, catabolizing extracellular peptides bound in the EPS matrix, while respiring nitrate, produced during anammox bacterial growth, to nitrite. Other heterotrophic bacteria may also contribute to the scavenging of detritus and peptides produced by anammox bacteria, and potentially use alternative electron donors, such as H2, acetate and formate to fuel their energy metabolism. These findings improve the understanding of major metabolic activities and interactions occurring between anammox and heterotrophic bacteria and offer the first transcriptional insights on ecosystem function in anammox granules.

6Sequencing of whole community DNA extracted from two separate biomass samples yielded a total of 41,063,036 reads after quality filtering (Supplementary Table 1). Co-assembly of the resulting reads using CLC Genomics Workbench generated a total of 83,770 contigs with an N50 of 2,376 bp and an N20 of 20,848 bp, accounting for 96% of the quality filtered DNA reads obtained from the second sampling event (Supplementary Table 2). Approximately 74% of the quality filtered mRNA reads could be mapped to this assembly, indicating that it also captured a large proportion of the metatranscriptome. Contigs were subsequently binned into population genomes based on tetranucleotide frequency and differential coverage using MetaBAT. This resulted in the recovery of 17 draft metagenome-assembled genomes (MAGs) affiliated with the phyla Planctomycetes, Chlorobi, Bacteroidetes, Chloroflexi, Proteobacteria and the candidate phylum Microgenomates (OP11; Table 1). Together, these genomes accounted for ∼82% and 59% of the total quality filtered DNA and mRNA reads obtained from the biomass samples, respectively, and therefore represented a major fraction of the microbial community present in the granular sludge reactor (Supplementary Table 2).

7Figure 1 shows a phylogenetic tree of the recovered population genomes based on the protein sequences of 37 conserved bacterial marker genes. Many of the recovered genomes, particularly CHB1 and CFX5, had high similarity to genomes recently recovered from a single-stage reactor that were shown to share high 16S rRNA gene similarity to organisms detected in other anammox systems. Surprisingly, the Chlorobi CHB1 genome was virtually identical to the Chlorobi OLB4 genome recovered by Speth et al., sharing an average nucleotide identity of 99.8% despite the significant differences in reactor operation (two-stage versus single-stage), wastewater influent composition (potato-processing wastewater versus anaerobic digester filtrate) and geographical location (USA versus Netherlands; Supplementary Data 1).

8We used the reads per kilobase per million mapped reads (RPKM) values for metagenomic reads and transcripts that mapped to each MAG as proxies for relative abundance and gene expression, respectively (Fig. 2). A summary of the metatranscriptomic sequencing statistics and read mapping can be found in Supplementary Table 2. Overall, gene expression corresponded with abundance in the anammox granules. Genomes affiliated with Brocadia (AMX1) and Chlorobi (CHB1) dominated the abundance and gene expression of the microbial community in the anammox granules (Fig. 2). AMX1 had a relative abundance and gene expression of ∼62 and 54%, respectively, whereas CHB1 had a relative abundance and gene expression of ∼21 and 26%, respectively. Other organisms that displayed moderate abundance and gene expression in the anammox ecosystem were affiliated with the phyla Chlorobi, Bacteroidetes, Chloroflexi and Proteobacteria (Fig. 2). Interestingly, several low abundance genomes were also observed to be particularly active based on gene expression (that is, PRO1, PRO2 and PLA1), as has been previously observed in other activated sludge ecosystems.

9To examine the functional potential and gene expression of the anammox community, open reading frames (ORFs) were predicted and annotated across each MAG. Subsequently, mRNA transcripts were mapped against all ORFs to identify functions that were highly expressed by the community during steady-state bioreactor operation (Supplementary Data 2). Metabolic pathways were then reconstructed for each of eight near-complete MAGs that displayed high abundance and/or gene expression in the community (>0.5%) using MetaPathways 2.5 (ref. ). A complete list of the inferred pathways and associated enzymes of each MAG can be found in Supplementary Data 3.

10As expected, genes involved in anammox metabolism from AMX1 were among the highest genes expressed in the community (Supplementary Data 2). Most notably, UTAMX1_1243, UTAMX1_1246 and UTAMX1_1249 annotated as hydrazine dehydrogenase (hdh), hydroxylamine oxidoreductase (hao) and hydrazine synthase subunit A (hzsA) were highly expressed, maintaining 10–20-fold coverage above median gene expression levels in the AMX1 genome (Supplementary Fig. 1). Six other hao-like genes were also identified in the AMX1 genome and displayed above median gene expression levels. Reciprocal best-BLAST searches confirmed that these genes were orthologous to hao-like genes from publically available anammox genomes (Supplementary Fig. 1). Among the Hao-like proteins, UTAMX1_1246 had the highest gene expression, suggesting it played an important role in anammox metabolism. This gene is orthologous to previously characterized Hao proteins from ‘Candidatus Kueneia stuttgartiensis' (kustc1061), ‘Candidatus Brocadia anammoxidans' and ‘Candidatus Jettenia caeni', which have been shown to oxidize hydroxylamine to nitric oxide and were also observed to be among the most abundant proteins expressed in these organisms. Phylogenetic analysis showed that AMX1 was closely related to both ‘Candidatus Brocadia sinica' and ‘Candidatus Brocadia fulgida' (Fig. 1). Similar to ‘Ca. Brocadia sinica', neither nirK nor nirS were found in the AMX1 genome. This is consistent with the recently proposed hydroxylamine-dependent anammox mechanism in ‘Ca. Brocadia sinica' that first reduces nitrite to hydroxylamine (instead of nitric oxide), and subsequently converts hydroxylamine and ammonium to hydrazine (Supplementary Fig. 1). It has been proposed that Hao-like proteins lacking a crosslinking tyrosine in the c-terminus may be involved in the initial nitrite reduction reaction to either nitric oxide or hydroxylamine, for which an enzyme has not yet been identified. The hao-like genes UTAMX1_1996, UTAMX1_1759 and UTAMX1_1192 expressed in the AMX1 genome that lack a crosslinking tyrosine might fulfil this role (Supplementary Fig. 1).

11All heterotrophic organisms in the anammox granules encoded capabilities for partial or full denitrification (Fig. 3). This was similar to a recent report from a single-stage reactor by Speth et al. Genes involved in denitrification were highly expressed in the genomes, consistent with oxidative phosphorylation coupled to nitrate respiration being a dominant form of heterotrophic energy metabolism in the anammox granules (Supplementary Data 2 and 4). All heterotrophs expressed respiratory nitrate reductase genes (narGHIJ) that reduce nitrate to nitrite and require the transport of nitrate into the cytoplasm (Fig. 3; Supplementary Data 4). However, it should be noted that the heme-containing membrane anchor subunit (narI) gene had no detectable expression in CHB2 or CFX2, suggesting that this enzyme had low activity in these organisms. A periplasmic nitrate reductase (napABCGH) that does not contain a coupling site for proton motive force (PMF) generation was also expressed by PRO2. Pathways for dissimilatory nitrate reduction to ammonia via pentaheme nitrite reductase genes (nrfHA) were expressed by the Chlorobi genome CHB2 and the Chloroflexi genome CFX2 (Fig. 3; Supplementary Data 4). Chloroflexi (CFX1 and CFX2) and Proteobacteria (PRO1 and PRO2) genomes also expressed genes for the reduction of nitrite to nitric oxide, either through a copper containing nitrite reductase (nirK) or cytochrome cd1 nitrite reductase (nirS; Fig. 3; Supplementary Data 4). The reduction of nitric oxide to nitrous oxide, either through the cytochrome c-dependent nitric oxide reductase (norBC) or the quinol-dependent nitric oxide reductase (norZ), was expressed by the BCD1, CFX2 and PRO1 genomes. Finally, Chlorobi (CHB1 and CHB2), Proteobacteria (PRO1 and PRO2), and the Bacteroidetes (BCD1) genomes expressed genes capable of reducing nitrous oxide to nitrogen gas via nitrous oxide reductase (nosZ). Together, these functions could facilitate a nitrite loop with anammox bacteria or support complete denitrification, thus enhancing overall nitrogen removal performance in the bioreactor.

12Figure 4 shows the gene expression profiles of major carbon and energy metabolic pathways across the anammox community. Only two low abundance heterotrophic organisms, BCD1 and PRO2, encoded known pathways for CO2 fixation that would permit mixotrophic growth (Fig. 4; Supplementary Data 3). BCD1 displayed above median gene expression of PEP carboxylase (UTBCD1_1894), which allows CO2 fixation into oxaloacetate (Supplementary Data 2); PRO2 displayed above median gene expression of Ribulose-1,5-bisphosphate carboxylase (RuBisCO, UTPRO2_2197), which would permit CO2 fixation via the Calvin-Benson-Bassham cycle (Fig. 4; Supplementary Data 3). While all heterotrophs, except CHB1, expressed genes encoding pyruvate:flavodoxin/ferredoxin oxidoreductase (PFOR) and 2-oxoglutarate ferredoxin oxidoreductase (KFOR) that are commonly implicated in the reductive TCA cycle, other key genes of this pathway, such as ATP-citrate lyase and fumarate reductase, were missing from the genomes. As PFOR and KFOR are known to participate in other metabolic pathways (for example, KFOR may participate in the oxidative TCA cycle), it is likely that these genes are not used for CO2 fixation by heterotrophs in the community.

13All organisms expressed genes involved in central carbon metabolic pathways, including glycolysis/gluconeogenesis, glycogen synthesis/degradation, the TCA cycle and the pentose phosphate pathway (PPP; Fig. 4; Supplementary Data 3). Surprisingly, the oxidative branch of the PPP was completely missing in both Chlorobi bacterium genomes (CHB1 and CHB2), suggesting that other enzymes are important for NADPH generation in these organisms. Genes for the oxidative PPP were also confirmed to be missing from the closely related Chlorobi bacterium OLB4 and OLB6 genomes, although they were present in the genomes of Ignavibacterium album and Melioribacter roseus. An alternative route for NADPH generation in CHB1 and CHB2 is likely through NAD(P)+ transhydrogenase, which transfers electrons from NADH to NADP+ and had above median gene expression in both genomes (Supplementary Data 2).

14Genes encoding pathways for pyruvate fermentation to acetate and CO2 via PFOR, phosphate acetyltransferase (pta) and acetate kinase (ack) were expressed in the CHB2 and PRO2 genomes (Fig. 4; Supplementary Data 3). While these genes could facilitate ATP synthesis via substrate level phosphorylation, it is also possible that the pta-ack reactions encoded by CHB2 and PRO2 are used for acetate consumption rather than production. These reactions form a lower affinity pathway for acetate assimilation compared to the high-affinity acetyl-CoA synthetase (acs) pathway. Under this scenario, however, only PRO2 could use acetate as a sole carbon source because it is the only genome that encodes the glyoxylate cycle (that is, malate synthase and isocitrate lyase). The AMP-forming acs gene was expressed in all genomes, whereas the ADP-forming acs gene was expressed in CHB1, CHB2, CFX1, CFX2 and PRO2 (Supplementary Data 2 and 3). Since these reactions are also reversible, with acetate formation resulting in ATP conservation, it is unclear whether the acs and/or pta-ack acetate pathways expressed by the microbial community are involved in acetate formation, acetate assimilation or both.

15The CHB2 genome also expressed an eight-gene cluster (UTCHB2_2599 to UTCHB2_2606) ferredoxin:NAD+ oxidoreductase (rnfCDGEAB) complex and two putative electron-bifurcating hydrogenase gene clusters (hydABC; UTCHB2_709-711 and UTCHB2_1049-1051). The Rnf complex could couple the formation of a proton gradient with ATP synthesis by oxidizing reduced ferredoxin produced from the electron-bifurcating hydrogenase and/or PFOR. This energy conservation mode was also identified in the closely related I. album and M. roseus genomes, however no orthologs were detected in the CHB1 genome.

16The ability to use H2, sulfite, formate, ethanol and other simple organic compounds as electron donors could also support the energy metabolism of the community (Fig. 4; Supplementary Data 3). CHB2, CFX2 and PRO2 expressed genes encoding a NAD-dependent hydrogenase that could be involved in H2 oxidation. CFX2 also expressed genes encoding a formate hydrogenlyase complex (UTCFX2_2421-2429) that directly links formate oxidation to H2 production. Organisms encoding the ability to use formate for energy conservation included CFX2 and PRO2 (Fig. 4; Supplementary Data 3). These organisms expressed genes encoding NAD-dependent formate dehydrogenase that catalyse the oxidation of formate to CO2, donating electrons to NAD+ that could be used to generate a PMF. PRO2 also expressed genes encoding formate dehydrogenase O (FDH-O; UTPRO2_1096 – UTPRO2_1099) that would donate electrons to the quinone pool and could also be used for PMF generation. Because no organism encoded the ability to produce formate via pyruvate formate lyase, it is possible that formate is made available to the community via CO2 reduction to formate by AMX1, which is the first step in CO2 fixation via the reductive acetyl-CoA pathway.

17Unexpectedly, genes encoding aerobic respiration machinery were expressed by some community members, despite the lack of oxygen in the bioreactor. CHB1, CFX2 and PRO2 expressed genes encoding a low-affinity aa3-type terminal cytochrome c oxidase; CHB2 expressed a high-affinity cytochrome d ubiquinol oxidase gene (Fig. 4; Supplementary Data 3). While this could allow these organisms to respire oxygen under aerobic conditions, it is also possible that these oxidases have evolved for protection against oxygen, perform some unknown function, or are post-transcriptionally repressed. All organisms were also found to express genes encoding superoxide dismutase (Fig. 4; Supplementary Data 3), which may additionally function to protect against oxidative stress.

18The breakdown of EPS produced by anammox bacteria has been proposed to be a major organic carbon source supporting heterotrophic growth in anammox granules. In agreement with this, genes encoding a wide range of peptidases were expressed in the heterotrophic genomes (Fig. 5; Supplementary Data 5). In particular, CHB1 and CHB2 encoded many extracellular subtilisin-like serine peptidases and metallopeptidases (Fig. 5; Supplementary Data 5). These genes were among the highest expressed ORFs in the genomes and may be involved with growth on proteinaceous substrates. Consistent with this function, the CHB1 and CHB2 genomes expressed genes for the transport and catabolism of peptides and amino acids to central carbon intermediates (Figs 5 and 6; Supplementary Data 2 and 3). This would allow CHB1 and CHB2 to use amino acids as a carbon and energy source, in addition to direct assimilation into protein biosynthesis. Interestingly, ranking of the amino acids according to their biosynthetic cost revealed that all organisms in the anammox granules lacked catabolic pathways for many amino acids with a high biosynthetic cost (Fig. 6; Supplementary Data 3). This may suggest that heterotrophic community members have been selected to preferentially degrade amino acids that are less costly to synthesize, while directly incorporating more costly amino acids into protein synthesis. Other heterotrophic genomes (BCD1, CFX1, CFX2, PRO2 and PRO1) also expressed genes involved in the transport and degradation of amino acids (Figs 5 and 6; Supplementary Data 2). These organisms could potentially take advantage of substrates made available by CHB1 and CHB2, or contribute to extracellular protein degradation as well.

19Microbial community members also expressed diverse genes involved in the hydrolysis of carbohydrate bonds (Supplementary Data 6). In particular, CHB2 and BCD1 expressed many ORFs encoding glycoside hydrolases, some of which were predicted to be extracellular. The high number of these genes in the CHB2 genome appeared to be a major function differentiating it from the CHB1 genome, which had the lowest number of predicted glycoside hydrolyases in the community (Supplementary Data 6). CHB2 expressed several genes that may act on glycosidic bonds found in polysaccharides, whereas BCD1 expressed a range of glycoside hydrolase genes that may use oligosaccharides, mucins and glycolipids as substrates and play an important role in carbohydrate breakdown for the community.

20Metabolite exchange of amino acids and vitamins is known to shape microbial community assembly. We observed that several of the abundant heterotrophic bacteria, namely CHB1, CHB2, CFX1 and CFX2, were missing pathways for the synthesis of many hydrophobic amino acids (Fig. 6; Supplementary Data 3). This was in contrast to the AMX1 genomes and other heterotrophic genomes (BCD1, PRO1 and PRO2) that expressed pathways for the synthesis of most amino acids.

21Several of the heterotrophic bacteria were also missing key genes involved in B-vitamin biosynthesis. The CHB1, CHB2 and CFX2 genomes were missing key genes involved in thiamin (vitamin B1) biosynthesis (thiamine-phosphate synthase and thiamine-monophosphate kinase), biotin (vitamin B7) biosynthesis (adenosylmethionine-8-amino-7-oxononanoate aminotransferase and bioin synthase) and adenosylcobalamin (vitamin B12) biosynthesis (cobalamin synthase and adenosylcobinamide-phosphate synthase; Supplementary Data 3). CFX1, BCD1, PRO2 and PRO1 also lacked key genes involved in adenosylcobalamin de novo synthesis. These genes were expressed in the AMX1 genome, suggesting that anammox bacteria may support B-vitamin requirements for the community. Genes encoding thiamin ABC-type transporters were expressed in the PRO1 and PRO2 genomes, whereas the CFX1, CFX2, PRO1 and PRO2 genomes expressed ABC-type transporters for vitamin B12 (btuFCD; Supplementary Data 2). The AMX1 (UTAMX1_2357), CHB1 (UTCHB1_731), CHB2 (UTCHB2_2713) and BCD1 (UTBCD1_3076) genomes expressed genes encoding the TonB-dependent vitamin B12 transporter (btuB), which facilitates vitamin B12 translocation across the outer membrane and may also be involved in uptake.

22The integration of metagenomic and metatranscriptomic sequencing allowed us to examine gene expression and microbial interactions in a lab-scale anammox community at the ecosystem-scale. Our results revealed that anammox bacteria affiliated with Brocadia (AMX1) and heterotrophic bacteria affiliated with Chlorobi (CHB1) dominated the abundance and gene expression of the anammox granule community. Other less abundant heterotrophic bacteria affiliated with Chloroflexi, Bacteroidetes and Proteobacteria also displayed significant gene expression that may contribute to ecosystem function. The high similarity between these genomes and the genomes recently recovered from a single-stage anammox reactor suggests that a core microbiome exists in anammox-based wastewater treatment systems, despite differences in reactor operation and influent wastewater composition.

23Remarkably, the CHB1 genome was nearly identical to the Chlorobi bacterium OLB4 genome recovered by Speth et al (99.8% average nucleotide identity) that achieved the second highest abundance in a full-scale partial-nitritation-anammox reactor from the Netherlands, suggesting that this organism interacts consistently with Brocadia sp. in anammox-based wastewater treatment systems. Figure 7 summarizes the major metabolic interactions proposed between AMX1 and CHB1 in anammox granules, based on the metabolic reconstruction and gene expression results obtained in this study. We posit that the ecological role of CHB1 in anammox granules is to degrade and catabolize extracellular peptides bound in the EPS matrix, while respiring nitrate produced during anammox bacterial growth to nitrite. This is supported by the high gene expression of respiratory nitrate reductase and extracellular peptidases in the CHB1 genome, as well as transporters and genes involved in amino acid catabolism. Interestingly, the CHB1 genome lacked catabolic pathways for amino acids that carry a high biosynthetic cost, suggesting that less costly amino acids (for example, glutamine) may be preferentially used as carbon and energy sources by heterotrophic bacteria in anammox granules. Such preferential degradation of amino acids has also been observed in anoxic water columns with sinking particulate organic carbon.

24Because CHB1 lacks many pathways for amino acid biosynthesis, it is likely that peptide and amino acid substrates used as carbon and energy sources originate from Brocadia sp., which derives carbon via CO2 fixation and expressed pathways for the synthesis of nearly all amino acids (Fig. 6). This is consistent with 14C-bicarbonate tracer experiments that showed heterotrophic bacteria in anammox bioreactors fed with no external organic carbon compounds degraded and utilized cellular components produced by anammox bacteria. The high transcript abundance of CHB1 in the metatranscriptome (25% relative community expression) also suggests that proteolysis by this organism is a significant component of carbon flux in anammox granules, possibly making amino acids available to other microbial community members. The number and high gene expression of amino acid transporters across the other heterotrophic genomes supports this interaction, suggesting that protein degradation by CHB1 may be a significant contributor to community assembly in anammox granules, similar to particle degrading taxa in marine ecosystems.

25It is also possible that Brocadia sp. supply B-vitamins to members of the microbial community, based on the expression of key genes involved in their biosynthesis and auxotrophies present across the heterotrophic genomes. The exchange of B-vitamins has been postulated to shape the structure and function of some microbial communities, such as those harboured in the human gut. The specific role vitamin exchange plays in shaping anammox community interactions currently remains unknown, and should therefore be investigated further in future studies.

26The oxidization of amino acids coupled to the reduction of nitrate to nitrite by CHB1 would serve to enhance overall nitrogen removal in the ecosystem through the completion of a nitrite loop with Brocadia sp., as proposed by Speth et al. Amino acid catabolism by CHB1 would also result in the release of ammonia, allowing it to be removed from the system by Brocadia sp., thus further enhancing overall nitrogen removal. It is also possible that nitrate reduction to nitrite by CHB1 and other heterotrophic organisms helps stimulate anammox metabolism. Recent microsensor studies have shown that nitrite concentrations can become limiting in the interior of anammox granules. Therefore, the ability of CHB1 to rapidly recycle nitrate back to nitrite may be critical for supporting anammox growth, which is required for their own proliferation.

27In addition to a nitrite loop, heterotrophic bacteria in anammox ecosystems may facilitate complete denitrification, either alone or based on metabolite exchange of denitrification intermediates with other organisms. Indeed, the segregation of denitrification intermediates across different organisms has been shown to reduce the accumulation of growth-inhibiting metabolites (for example, nitrite) and potentially accelerate denitrification rates. The close spatial association of organisms in anammox granules may encourage such a distributed metabolic network. While complete denitrification would also serve to remove nitrogen from the ecosystem, it may potentially be detrimental to anammox bacterial growth because it would consume nitrite, rather than make it available for Brocadia sp. Specifically, the heterotrophic organisms capable of using nitrite as an electron acceptor (CHB2, CFX1, CFX2, PRO1 and PRO2) via nitrite reductase (NirK/S or NrfHA) may compete with Brocadia sp., similar to NOB in single-stage partial nitritation-anammox bioreactors.

28Aside from proteinaceous substrates, many organisms in the anammox community were capable of using other electron donors for energy conservation, including acetate, H2, formate or simple organic compounds. It is possible that these substrates are produced during fermentation reactions and/or the breakdown of polysaccharide compounds present in the EPS. Chlorobi organisms related to CHB2 may be the biggest contributor to this process, based on expression of PFOR, hydrogenases and ADP-forming ACS involved in H2 and acetate production. The availability of such substrates may contribute to niche differentiation by heterotrophs in the community, while also providing additional reducing equivalents to fuel anammox metabolism. Nevertheless, oxidative phosphorylation coupled to nitrate respiration, rather than fermentation, was found to be the main energy metabolism expressed across the heterotrophic genomes, consistent with the high availability of nitrate and nitrite in the ecosystem. As such, the significance of H2 and acetate metabolism in anammox granules still remains to be determined.

29In conclusion, our combined approach of genome-centric metagenomics with metatranscriptomics allowed us to obtain a much deeper understanding of ecosystem function in anammox granules that has broad implications for anammox-based wastewater treatment. This study provides a robust analysis of the metabolic networks underlying several poorly characterized taxa frequently detected in anammox granules (Planctomycetes, Chlorobi, Bacteroidetes, Chloroflexi and Proteobacteria) and also sheds light on their potential ecological roles and interactions. Our results implicate members of the Chlorobi affiliated with CHB1 that are broadly found in anammox granules as highly active protein degraders, possibly liberating amino acids bound in the EPS matrix for themselves and other organisms in the community. These amino acids likely originate from anammox bacteria, which also provide nitrate and potentially essential B-vitamins to the community, supporting anaerobic respiration and growth. Other heterotrophic community members also appear to contribute to the scavenging of detritus and peptides produced by anammox bacteria, and may potentially use alternative electron donors, such as H2, acetate and formate. These substrates may be made available through the hydrolysis and fermentation of EPS and detritus.

30While further studies are still required to determine the specific niche space occupied by members of the community, the metabolic characterization provided here advances our insight on the ecological roles played by both anammox and heterotrophic bacteria in anammox granules. The high expression of enigmatic functions by the community, including hypothetical proteins and HAO-like proteins, together with our incomplete understanding of the biochemical pathways underlying EPS biosynthesis, highlights the need for continued work on these novel and industrially important ecosystems.

31A 5 l continuously fed sequencing batch reactor (anammox bioreactor) has been operating in the Goel Laboratory at the University of Utah for the past 6 years, achieving stable anammox performance. The anammox bioreactor was originally inoculated with anammox biomass from the City College of New York (Civil Engineering Department) that was enriched from activated sludge. The bioreactor was fed with anaerobic digester filtrate obtained from the belt filter press of a local wastewater treatment plant (CVWRF, Salt Lake City, UT, USA). Before entering the anammox reactor, the filtrate was first passed through a nitritation reactor, where approximately half of the influent ammonia was oxidized to nitrite. The nitritation reactor had a working volume of 2 l with a hydraulic retention time of 1 day. Dissolved oxygen in the reactor was maintained at 0.5±0.5 mg l−1 and bicarbonate present in the influent filtrate kept the pH buffered between 6.5 and 8.5. The anammox bioreactor was operated at room temperature with a hydraulic retention time of 2 days and no biomass wasting, which encouraged the formation of dense granules. The pH of the anammox bioreactor was maintained at between 6.8 and 8.0 and the reactor contents were completely mixed using a stir bar and gas sparging. Anaerobic conditions in the anammox bioreactor were maintained by continuously purging a mixture of 95% nitrogen gas and 5% carbon dioxide. Table 2 summarizes the process performance data of the nitritation and anammox bioreactors during steady-state operation.

32Two independent biomass samples were collected from the anammox bioreactor at different time points during periods of high total nitrogen removal efficiency (>80%). DNA was extracted separately from each biomass sample to improve the recovery of MAGs and total RNA was extracted from the second biomass sample to examine steady-state community gene expression. 500 mg of biomass was collected from the reactor, centrifuged at 4500, r.p.m. for 5 min at 4 °C to remove the supernatant, and instantly flash frozen in liquid nitrogen. Total genomic DNA and RNA was extracted from biomass pellets using the PowerMax Soil DNA Isolation Kit (MoBio Laboratories, USA) and the PureLink RNA mini kit (Life Technology, NY, USA), respectively, according to the manufacturers' protocols. Genomic DNA was quality-checked using agarose gel electrophoresis and a Nanodrop ND-2000c (Thermo Fisher Scientific, USA). Following RNA extraction, residual genomic DNA was removed from total RNA using an on-column PureLink DNase set (Life Technologies, NY, USA). Total RNA quality and quantity was subsequently checked using the Bioanalyzer RNA 6000 Nano Assay (Agilent, Santa Clara, CA, USA) to ensure only high-quality nucleic acids were used for downstream analysis. Total RNA of 100 ng was used to construct strand specific RNA-Seq libraries with the Encore Complete Prokaryotic RNA-Seq DR Multiplex System (NuGEN, San Carlos, CA, USA). Non-rRNA in RNA-Seq libraries were enriched by selective priming during the first strand cDNA synthesis reaction, as well as in the final library construction steps using manufacturer's protocols. DNA from the first sample was sequenced on the Illumina MiSeq platform (Illumina, CA, USA) to generate 300 bp paired-end reads (550 bp mean insert size). DNA and RNA from the second sample were sequenced on the Illumina HiSeq 2000 platform to generate 125 bp paired-end reads (180 bp mean insert size) at greater sequencing depth. All DNA and RNA sequencing was performed at the Huntsman Cancer Institute (HCI), University of Utah. Raw DNA and RNA sequences can be found on the National Center for Biotechnology Information (NCBI) website under BioProject PRJNA343219.

33Raw paired-end reads from the MiSeq and HiSeq platforms were initially filtered using Sickle v1.33 (ref. ) based on a minimum quality score of 20, a minimum sequence length of 100 bp and allowing for no ambiguous bases. Paired-end reads were then merged using FLASH v1.2.11 (ref. ) and co-assembled using the de novo assembler of CLC Genomics Workbench v7.0.3 (CLCbio, Arhus, Denmark) based on default parameters (word size=20 and bubble size=50). The per-base coverage depth across all contigs was calculated by mapping raw reads from each sample against the co-assembled contigs using BBMap v35.92 (https://sourceforge.net/projects/bbmap/) with the parameters ‘minid=0.95' and ‘ambig=random'. Resulting mapping files were subsequently used by MetaBAT v0.26.3 to bin metagenomic contigs into draft genomes, based on the ‘sensitive' parameters. CheckM v1.0.3 (ref. ) was used to estimate the contamination and completeness of each draft genome based on 111 essential single-copy marker genes. Draft genomes have been deposited to GenBank under the accession numbers provided in Supplementary Table 3.

34Phylogenetic analysis of the recovered draft genomes was accomplished using Phylosift v1.0.1, based on a set of 37 universal single-copy marker genes. The taxonomic affiliation of each draft genome was determined using the Phylosift ‘all' command. Marker genes were also identified in 33 publicly available genomes closely related to the recovered draft genomes and used to build a phylogenetic tree. Marker genes were concatenated and aligned using Phylosift and a maximum likelihood tree was generated using RAxML v8.2.4 with the automatic protein model assignment algorithm (PROTGAMMAAUTO) and 100 bootstraps.

35Metabolic reconstruction of the recovered draft genomes was performed using MetaPathways v2.5 (ref. ). Briefly, ORFs were predicted using Prodigal v2.0 (ref. ), based on a minimum nucleotide length of 60, and queried against the SEED subsystems (accessed March 2013), Clusters of Orthologous Groups (COG, accessed December 2013), RefSeq (accessed September 2015) and MetaCyc (accessed October 2011) protein databases using the optimized LAST algorithm (E value, 1E-6) for functional annotation. Contig nucleotide sequences were also queried against the SILVA small subunit (SSU) 123 database to identify the taxonomy of recovered 16S rRNA genes. Annotated genomes were then used to reconstruct the metabolic network of each organism using Pathway Tools and the MetaCyc database. Pathway/Genome Databases were created for each genome, where pathway inference was based on a set of rules used by the Pathway Tools prediction algorithm Pathologic, including the presence of all key reactions and the completeness of the reconstructed pathway. All inferred pathways were then manually curated to verify predictions made by Pathologic.

36Metatranscriptomic reads were quality filtered and merged as described above. Subsequently, rRNA sequences were filtered from the metatranscriptomic data set using SortMeRNA v2.0, based on multiple rRNA databases for bacterial, archaeal and eukaryotic sequences. Resulting non-rRNA reads were mapped back to all assembled metagenomic contigs using BBMap v35.92 with the parameters ‘minid=0.95', which specifies a minimum alignment identify of 95% corresponding to the well-established criteria for identifying microbial species using average nucleotide identify and ‘ambig=random', which ensures reads with multiple top-scoring mapping locations are assigned randomly to a single location. Read counts were calculated for each predicted ORF using htseq-count v0.6.1 with the ‘intersection strict' parameter and normalized by sequencing depth and ORF length and expressed as RPKM values. This allowed for the profiling of gene expression across the recovered draft genomes. Relative gene expression within each genome was calculated by relativizing the expression (RPKM value) of each ORF by the median RPKM value calculated across the genome. Pathway expression levels were calculated based on averaging the RPKM values for each reaction in a given pathway. For reactions with multiple genes, the highest expressed gene between multi-copy genes was selected, whereas the lowest expressed gene between multi-gene enzyme complexes (for example, hydrazine synthase, hzsABC) was selected.

37The relative abundance and gene expression of the recovered genomes was calculated from the total number of DNA or mRNA reads that mapped to the genome, divided by the genome length (read count/genome size).

38Carbohydrate hydrolases and peptidases potentially involved in the breakdown of EPS were identified in each genome based on BLASTP searches against the CAZy (accessed on September 2014) and MEROPS release 10.0 (ref. ) databases, respectively. Manual searches were also performed across all genomes based on gene annotations, and protein domains with hydrolytic activity were subsequently confirmed based on queries against the NCBI's Conserved Domain Database and Pfam database using MOTIF (http://www.genome.jp/tools/motif/). The subcellular location of identified proteins was predicted using CELLO v2.5, which classifies proteins using support vector machines trained by multiple feature vectors based on n-peptide composition.

39All data associated with this project can be found at the NCBI under BioProject PRJNA343219. Illumina HiSeq metagenomic data can be found under BioSample SAMN05785375, Illumina MiSeq metagenomic data can be found under BioSample SAMN05785373, and Illumina HiSeq metatranscriptomic data can be found under BioSample SAMN05785376. Annotated GenBank files for the whole genome sequences (MAGs) described in this study can be found under accession numbers listed in Supplementary Table 3.

40How to cite this article: Lawson, C. E. et al. Metabolic network analysis reveals microbial community interactions in anammox granules. Nat. Commun. 8, 15416 doi: 10.1038/ncomms15416 (2017).

41Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.