Metabolic network analysis reveals microbial community interactions in anammox granules — TEI-tagged view

Christopher E. LawsonpersName, Sha WupersName, Ananda S. BhattacharjeepersName, Joshua J. HamiltonpersName, Katherine D. McMahonpersName, Ramesh GoelpersName, Daniel R. NoguerapersName (2017date). “Metabolic network analysis reveals microbial community interactions in anammox granules.” Nature Communicationstitle·j.

Object of study — the metabolic network and interactions of anammoxterm·key and heterotrophic bacteriaterm·key in anammoxterm·key granules

Research question — what role heterotrophic bacteriaterm·key play in anammoxterm·key granules

Contribution — 17 draft genomes with metabolic reconstruction and metatranscriptomicsterm·key, proposing Chlorobiterm·key as protein degraders

Full text — tagged

The reference list is not part of the text we hold for this paper. Our copy was extracted from the publisher's page and the bibliography did not survive that extraction, so the article below ends before its references. Nothing has been removed by us.

Metabolic network analysis reveals microbial community interactions in anammoxterm·key granules Abstract Microbial communities mediating anaerobic ammonium oxidation (anammoxterm·key) represent one of the most energy-efficient environmental biotechnologies for nitrogen removal from wastewater. However, little is known about the functional role heterotrophic bacteriaterm·key play in anammoxterm·key granules. Here, we use genome-centric metagenomicsterm·key to recover 17 draft genomes of anammoxterm·key and heterotrophic bacteriaterm·key from a laboratory-scale anammoxterm·key bioreactor. We combine metabolic network reconstruction with metatranscriptomicsterm·key to examine the gene expressionterm·var of anammoxterm·key and heterotrophic bacteriaterm·key and to identify their potential interactions. We find that Chlorobiterm·key - affiliated bacteria may be highly active protein degraders, catabolizing extracellular peptides while recycling nitrate to nitrite. Other heterotrophs may also contribute to scavenging of detritus and peptides produced by anammoxterm·key bacteria, and potentially use alternative electron donors, such as H 2 , acetate and formate. Our findings improve the understanding of metabolic activities and interactions between anammoxterm·key and heterotrophic bacteriaterm·key and offer the first transcriptional insights on ecosystem function in anammoxterm·key granules.

Microbial communities mediating anaerobic ammonium oxidation (anammoxterm·key) 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 denitrificationterm·key 1 . In practice, anammox-based wastewater treatmentterm·key systems are combined with a nitritation step in either a single-stage 2 or two-stage bioreactor system 3 . In these engineered ecosystems, a fraction of the ammonium is first oxidized to nitrite by aerobic ammonia oxidizing bacteria. Subsequently, anammoxterm·key bacteria anaerobically oxidize the remaining ammonium directly to nitrogen gas using the produced nitrite as a terminal electron acceptor 4 5 .

Five genera of anammoxterm·key bacteria have been discovered to date, including Kuenenia , Brocadiaterm·key , Anammoxoglobus and Jettenia commonly found in activated sludge, and Scalindua commonly found in marine environments 6 . These lineages all have ‘Candidatus' status as they do not exist in pure culture and must be grown in laboratory enrichments. Because anammoxterm·key bacteria have a slow growth rate 7 , either biofilm reactors that use carrier media or granular sludge reactors are used to retain sufficient biomass in the system 8 . Such reactor configurations support the formation of dense microbial communities that can be readily separated from the liquid wastewater and enriched in the bioreactor. Anammoxterm·key granules consist of a mixture of cell aggregates and abiotic particles embedded within a matrix of organic extracellular polymeric substancesterm·dom (EPSterm·dom) 9 10 11 . EPSterm·dom present in anammoxterm·key 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 11 . The high aggregate density also serves to limit oxygen diffusion into the granule interior, allowing for the proliferation of anammoxterm·key and other anaerobic bacteria.

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

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

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

Results Metagenomicterm·key sequencing and binning Sequencing 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 genomesterm·mth (MAGsterm·mth) affiliated with the phyla Planctomycetes, Chlorobiterm·key, 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 ).

Figure 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 anammoxterm·key systems 16 . Surprisingly, the Chlorobiterm·key CHB1 genome was virtually identical to the Chlorobiterm·key OLB4 genome recovered by Speth et al . 16 , 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 ).

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

Metabolic reconstruction of microbial community To examine the functional potential and gene expressionterm·var of the anammoxterm·key community, open reading frames (ORFs) were predicted and annotated across each MAGterm·mth. 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 MAGsterm·mth that displayed high abundance and/or gene expressionterm·var in the community (>0.5%) using MetaPathways 2.5 (ref. 20 ). A complete list of the inferred pathways and associated enzymes of each MAGterm·mth can be found in Supplementary Data 3 .

Anammoxterm·key metatranscriptomicterm·key insights As expected, genes involved in anammoxterm·key 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 expressionterm·var 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 expressionterm·var levels. Reciprocal best-BLAST searches confirmed that these genes were orthologous to hao -like genes from publically available anammoxterm·key genomes ( Supplementary Fig. 1 ). Among the Hao-like proteins, UTAMX1_1246 had the highest gene expressionterm·var, suggesting it played an important role in anammoxterm·key metabolism. This gene is orthologous to previously characterized Hao proteins from ‘ Candidatus Kueneia stuttgartiensis' (kustc1061) 5 21 , ‘ Candidatus Brocadiaterm·key anammoxidans' 22 and ‘ Candidatus Jettenia caeni' 23 , 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 Brocadiaterm·key sinica' 24 and ‘ Candidatus Brocadiaterm·key fulgida' 25 ( Fig. 1 ). Similar to ‘ Ca. Brocadiaterm·key sinica', neither nirK nor nirS were found in the AMX1 genome. This is consistent with the recently proposed hydroxylamine-dependent anammoxterm·key mechanism in ‘ Ca. Brocadiaterm·key sinica' that first reduces nitrite to hydroxylamine (instead of nitric oxide), and subsequently converts hydroxylamine and ammonium to hydrazine 26 ( 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 5 or hydroxylamine 26 , 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 ).

Denitrificationterm·key gene expressionterm·var in anammoxterm·key granules All heterotrophic organisms in the anammoxterm·key granules encoded capabilities for partial or full denitrificationterm·key ( Fig. 3 ). This was similar to a recent report from a single-stage reactor by Speth et al . 16 Genes involved in denitrificationterm·key were highly expressed in the genomes, consistent with oxidative phosphorylationterm·dom coupled to nitrate respiration being a dominant form of heterotrophic energy metabolism in the anammoxterm·key 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 Chlorobiterm·key 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, Chlorobiterm·key (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 loopterm·key with anammoxterm·key bacteria or support complete denitrificationterm·key, thus enhancing overall nitrogen removal performance in the bioreactor.

Carbon and energy metabolism of heterotrophic bacteriaterm·key Figure 4 shows the gene expressionterm·var profiles of major carbon and energy metabolic pathways across the anammoxterm·key community. Only two low abundance heterotrophic organisms, BCD1 and PRO2, encoded known pathways for CO 2 fixation that would permit mixotrophic growth ( Fig. 4 ; Supplementary Data 3 ). BCD1 displayed above median gene expressionterm·var of PEP carboxylase (UTBCD1_1894), which allows CO 2 fixation into oxaloacetate ( Supplementary Data 2 ); PRO2 displayed above median gene expressionterm·var of Ribulose-1,5-bisphosphate carboxylase (RuBisCO, UTPRO2_2197), which would permit CO 2 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 27 ), it is likely that these genes are not used for CO 2 fixation by heterotrophs in the community.

All 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 Chlorobiterm·key 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 Chlorobiterm·key bacterium OLB4 and OLB6 genomes 16 , although they were present in the genomes of Ignavibacterium album 28 and Melioribacter roseus 29 . An alternative route for NADPH generation in CHB1 and CHB2 is likely through NAD(P) + transhydrogenase 30 , which transfers electrons from NADH to NADP + and had above median gene expressionterm·var in both genomes ( Supplementary Data 2 ).

Genes encoding pathways for pyruvate fermentationterm·dom to acetate and CO 2 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 31 . 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.

The 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 32 . This energy conservation mode was also identified in the closely related I. album 28 and M. roseus 29 genomes, however no orthologs were detected in the CHB1 genome.

The ability to use H 2 , 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 H 2 oxidation. CFX2 also expressed genes encoding a formate hydrogenlyase complex (UTCFX2_2421-2429) that directly links formate oxidation to H 2 production 33 . 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 CO 2 , 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 34 . 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 CO 2 reduction to formate by AMX1, which is the first step in CO 2 fixation via the reductive acetyl-CoA pathway.

Unexpectedly, 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 35 , 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.

Amino acid and carbohydrate catabolism The breakdown of EPSterm·dom produced by anammoxterm·key bacteria has been proposed to be a major organic carbon source supporting heterotrophic growth in anammoxterm·key granules 11 17 . 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 36 37 . 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 38 revealed that all organisms in the anammoxterm·key 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.

Microbial 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.

Vitamin and amino acid auxotrophy Metabolite exchange of amino acids and vitamins is known to shape microbial community assembly 39 . We observed that several of the abundant heterotrophic bacteriaterm·key, 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.

Several of the heterotrophic bacteriaterm·key 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 anammoxterm·key 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.

Discussion The integration of metagenomicterm·key and metatranscriptomicterm·key sequencing allowed us to examine gene expressionterm·var and microbial interactions in a lab-scale anammoxterm·key community at the ecosystem-scale. Our results revealed that anammoxterm·key bacteria affiliated with Brocadiaterm·key (AMX1) and heterotrophic bacteriaterm·key affiliated with Chlorobiterm·key (CHB1) dominated the abundance and gene expressionterm·var of the anammoxterm·key granule community. Other less abundant heterotrophic bacteriaterm·key affiliated with Chloroflexi, Bacteroidetes and Proteobacteria also displayed significant gene expressionterm·var that may contribute to ecosystem function. The high similarity between these genomes and the genomes recently recovered from a single-stage anammoxterm·key reactor 16 suggests that a core microbiometerm·dom exists in anammox-based wastewater treatmentterm·key systems, despite differences in reactor operation and influent wastewater composition.

Remarkably, the CHB1 genome was nearly identical to the Chlorobiterm·key bacterium OLB4 genome recovered by Speth et al 16 (99.8% average nucleotide identity) that achieved the second highest abundance in a full-scale partial-nitritation-anammox reactor from the Netherlands 16 40 , suggesting that this organism interacts consistently with Brocadiaterm·key sp. in anammox-based wastewater treatmentterm·key systems. Figure 7 summarizes the major metabolic interactions proposed between AMX1 and CHB1 in anammoxterm·key granules, based on the metabolic reconstruction and gene expressionterm·var results obtained in this study. We posit that the ecological role of CHB1 in anammoxterm·key granules is to degrade and catabolize extracellular peptides bound in the EPSterm·dom matrix, while respiring nitrate produced during anammoxterm·key bacterial growth to nitrite. This is supported by the high gene expressionterm·var 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 bacteriaterm·key in anammoxterm·key granules. Such preferential degradation of amino acids has also been observed in anoxic water columns with sinking particulate organic carbon 41 .

Because 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 Brocadiaterm·key sp. , which derives carbon via CO 2 fixation and expressed pathways for the synthesis of nearly all amino acids ( Fig. 6 ). This is consistent with 14 C-bicarbonate tracer experiments that showed heterotrophic bacteriaterm·key in anammoxterm·key bioreactors fed with no external organic carbon compounds degraded and utilized cellular components produced by anammoxterm·key bacteria 42 . 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 anammoxterm·key granules, possibly making amino acids available to other microbial community members. The number and high gene expressionterm·var 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 anammoxterm·key granules, similar to particle degrading taxa in marine ecosystems 43 .

It is also possible that Brocadiaterm·key 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 44 . The specific role vitamin exchange plays in shaping anammoxterm·key community interactions currently remains unknown, and should therefore be investigated further in future studies.

The 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 loopterm·key with Brocadiaterm·key sp ., as proposed by Speth et al . 16 Amino acid catabolism by CHB1 would also result in the release of ammonia, allowing it to be removed from the system by Brocadiaterm·key sp ., thus further enhancing overall nitrogen removal. It is also possible that nitrate reduction to nitrite by CHB1 and other heterotrophic organisms helps stimulate anammoxterm·key metabolism. Recent microsensor studies have shown that nitrite concentrations can become limiting in the interior of anammoxterm·key granules 45 . Therefore, the ability of CHB1 to rapidly recycle nitrate back to nitrite may be critical for supporting anammoxterm·key growth, which is required for their own proliferation.

In addition to a nitrite loopterm·key, heterotrophic bacteriaterm·key in anammoxterm·key ecosystems may facilitate complete denitrificationterm·key, either alone or based on metabolite exchange of denitrificationterm·key intermediates with other organisms. Indeed, the segregationterm·dom of denitrificationterm·key intermediates across different organisms has been shown to reduce the accumulation of growth-inhibiting metabolitesterm·dom (for example, nitrite) and potentially accelerate denitrificationterm·key rates 46 . The close spatial association of organisms in anammoxterm·key granules may encourage such a distributed metabolic network 47 48 . While complete denitrificationterm·key would also serve to remove nitrogen from the ecosystem, it may potentially be detrimental to anammoxterm·key bacterial growth because it would consume nitrite, rather than make it available for Brocadiaterm·key 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 Brocadiaterm·key sp. , similar to NOBterm·dom in single-stage partial nitritation-anammox bioreactors.

Aside from proteinaceous substrates, many organisms in the anammoxterm·key community were capable of using other electron donors for energy conservation, including acetate, H 2 , formate or simple organic compounds. It is possible that these substrates are produced during fermentationterm·dom reactions and/or the breakdown of polysaccharide compounds present in the EPSterm·dom. Chlorobiterm·key organisms related to CHB2 may be the biggest contributor to this process, based on expression of PFOR, hydrogenases and ADP-forming ACS involved in H 2 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 anammoxterm·key metabolism 49 . Nevertheless, oxidative phosphorylationterm·dom coupled to nitrate respiration, rather than fermentationterm·dom, 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 H 2 and acetate metabolism in anammoxterm·key granules still remains to be determined.

In conclusion, our combined approach of genome-centric metagenomicsterm·key with metatranscriptomicsterm·key allowed us to obtain a much deeper understanding of ecosystem function in anammoxterm·key granules that has broad implications for anammox-based wastewater treatmentterm·key. This study provides a robust analysis of the metabolic networks underlying several poorly characterized taxa frequently detected in anammoxterm·key granules (Planctomycetes, Chlorobiterm·key, Bacteroidetes, Chloroflexi and Proteobacteria) and also sheds light on their potential ecological roles and interactions. Our results implicate members of the Chlorobiterm·key affiliated with CHB1 that are broadly found in anammoxterm·key granules as highly active protein degraders, possibly liberating amino acids bound in the EPSterm·dom matrix for themselves and other organisms in the community. These amino acids likely originate from anammoxterm·key 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 anammoxterm·key bacteria, and may potentially use alternative electron donors, such as H 2 , acetate and formate. These substrates may be made available through the hydrolysis and fermentationterm·dom of EPSterm·dom and detritus.

While 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 anammoxterm·key and heterotrophic bacteriaterm·key in anammoxterm·key 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 EPSterm·dom biosynthesis, highlights the need for continued work on these novel and industrially important ecosystems.

Methods Bioreactor operation A 5 l continuously fed sequencing batch reactor (anammoxterm·key bioreactor) has been operating in the GoelpersName Laboratory at the University of Utah for the past 6 years, achieving stable anammoxterm·key performance 50 . The anammoxterm·key bioreactor was originally inoculated with anammoxterm·key 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 treatmentterm·key plant (CVWRF, Salt Lake City, UT, USA). Before entering the anammoxterm·key 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 50 . 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 anammoxterm·key 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 anammoxterm·key 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 anammoxterm·key bioreactor were maintained by continuously purging a mixture of 95% nitrogen gas and 5% carbon dioxide 50 . Table 2 summarizes the process performance data of the nitritation and anammoxterm·key bioreactors during steady-state operation.

DNA and RNA sequencing Two independent biomass samples were collected from the anammoxterm·key 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 MAGsterm·mth and total RNA was extracted from the second biomass sample to examine steady-state community gene expressionterm·var. 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.

Metagenomicterm·key assembly and binning Raw paired-end reads from the MiSeq and HiSeq platforms were initially filtered using Sickle v1.33 (ref. 51 ) 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. 52 ) 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 metagenomicterm·key contigs into draft genomes, based on the ‘sensitive' parameters 53 . CheckM v1.0.3 (ref. 54 ) was used to estimate the contamination and completeness of each draft genome based on 111 essential single-copy marker genes 55 . Draft genomes have been deposited to GenBank under the accession numbers provided in Supplementary Table 3 .

Phylogenetic analysis of recovered draft genomes Phylogenetic analysis of the recovered draft genomes was accomplished using Phylosift v1.0.1, based on a set of 37 universal single-copy marker genes 56 . 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.

Genome annotation and metabolic reconstruction Metabolic reconstruction of the recovered draft genomes was performed using MetaPathways v2.5 (ref. 20 ). Briefly, ORFs were predicted using Prodigal v2.0 (ref. 57 ), 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 58 . 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 59 . 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 60 . All inferred pathways were then manually curated to verify predictions made by Pathologic.

Metatranscriptomicterm·key analysis Metatranscriptomicterm·key reads were quality filtered and merged as described above. Subsequently, rRNA sequences were filtered from the metatranscriptomicterm·key data set using SortMeRNA v2.0, based on multiple rRNA databases for bacterial, archaeal and eukaryotic sequences 61 . Resulting non-rRNA reads were mapped back to all assembled metagenomicterm·key 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 62 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 63 and normalized by sequencing depth and ORF length and expressed as RPKM values 64 . This allowed for the profiling of gene expressionterm·var across the recovered draft genomes. Relative gene expressionterm·var 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.

The relative abundance and gene expressionterm·var 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).

Carbohydrate hydrolase and peptidase identification Carbohydrate hydrolases and peptidases potentially involved in the breakdown of EPSterm·dom were identified in each genome based on BLASTP searches against the CAZy (accessed on September 2014) 65 and MEROPS release 10.0 (ref. 66 ) 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 67 and Pfam database 68 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 69 .

Data availability statement All data associated with this project can be found at the NCBI under BioProject PRJNA343219. Illumina HiSeq metagenomicterm·key data can be found under BioSample SAMN05785375, Illumina MiSeq metagenomicterm·key data can be found under BioSample SAMN05785373, and Illumina HiSeq metatranscriptomicterm·key data can be found under BioSample SAMN05785376. Annotated GenBank files for the whole genome sequences (MAGsterm·mth) described in this study can be found under accession numbers listed in Supplementary Table 3 .

Additional information How to cite this article: LawsonpersName, C. E. et al . Metabolic network analysis reveals microbial community interactions in anammoxterm·key granules. Nat. Commun. 8, 15416 doi: 10.1038/ncomms15416 (2017).

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

Supplementary Material Supplementary Figure and Supplementary Tables Genome comparison of the CHB1 and OLB4 genomes recovered from the lab-scale reactor employed in this study and the Olburgen PNA reactor (Speth et al., 2016), respectively.

Predicted open reading frame annotations and gene expressionterm·var for each MAGterm·mth recovered from the anammoxterm·key community.

MetaCyc pathway/genome database for each MAGterm·mth recovered from the anammoxterm·key community.

Summary of genes and corresponding gene expressionterm·var values (RPKM) involved in denitrificationterm·key across MAGsterm·mth.

Summary of predicted peptidases across MAGsterm·mth based on annotations against the MEROPs database. * indicates peptidases identified using the CDD database or MOTIF.

Summary of genes predicted to be involved in carbohydrate metabolism based on annotations against the CAZy database.