Metabolic network analysis reveals microbial community interactions in anammox granules

Abstract

Microbial communities mediating anaerobic ammonium oxidation (anammox) represent one of the most energy-efficient environmental biotechnologies for nitrogen removal from wastewater. However, little is known about the functional role heterotrophic bacteria play in anammox granules. Here, we use genome-centric metagenomics to recover 17 draft genomes of anammox and heterotrophic bacteria from a laboratory-scale anammox bioreactor. We combine metabolic network reconstruction with metatranscriptomics to examine the gene expression of anammox and heterotrophic bacteria and to identify their potential interactions. We find that Chlorobi-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 anammox bacteria, and potentially use alternative electron donors, such as H₂ , acetate and formate. Our findings improve the understanding of metabolic activities and interactions between anammox and heterotrophic bacteria and offer the first transcriptional insights on ecosystem function in anammox granules.

ARTICLE

Received 21 Sep 2016 | Accepted 23 Mar 2017 | Published 31 May 2017

Christopher E. Lawson1, Sha Wu2, Ananda S. Bhattacharjee2, Joshua J. Hamilton3, Katherine D. McMahon1,3, Ramesh Goel2 & Daniel R. Noguera1,4

Microbial communities mediating anaerobic ammonium oxidation (anammox) represent one of the most energy-efcient environmental biotechnologies for nitrogen removal from wastewater. However, little is known about the functional role heterotrophic bacteria play in anammox granules. Here, we use genome-centric metagenomics to recover 17 draft genomes of anammox and heterotrophic bacteria from a laboratory-scale anammox bioreactor. We combine metabolic network reconstruction with metatranscriptomics to examine the gene expression of anammox and heterotrophic bacteria and to identify their potential interactions. We nd that Chlorobi-afliated 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 anammox bacteria, and potentially use alternative electron donors, such as H2, acetate and formate. Our ndings improve the understanding of metabolic activities and interactions between anammox and heterotrophic bacteria and offer the rst transcriptional insights on ecosystem function in anammox granules.

DOI: 10.1038/ncomms15416 OPEN

Metabolic network analysis reveals microbial community interactions in anammox granules

1 Department of Civil and Environmental Engineering, University of WisconsinMadison, Madison, Wisconsin 53706, USA. 2 Department of Civil and Environmental Engineering, University of Utah, Salt Lake City, Utah 84112, USA. 3 Department of Bacteriology, University of WisconsinMadison, Madison, Wisconsin 53706, USA. 4 Great Lakes Bioenergy Research Center, Wisconsin Energy Institute, University of WisconsinMadison, Madison, Wisconsin 53726, USA. Correspondence and requests for materials should be addressed to R.G. (email: mailto:[email protected]

Web End [email protected] ) or to D.R.N. (email: mailto:[email protected]

Web End [email protected] ).

NATURE COMMUNICATIONS | 8:15416 | DOI: 10.1038/ncomms15416 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 1

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms15416

Microbial communities mediating anaerobic ammonium oxidation (anammox) represent one of the most energy-efcient environmental biotechnologies for

nitrogen removal from wastewater. The process is typically used to treat high strength ammonium wastewaters and offers signicant cost savings compared to conventional nitrogen removal processes that require energy-intensive aeration for nitrication and also consume large quantities of organic carbon during denitrication1. In practice, anammox-based wastewater treatment systems are combined with a nitritation step in either a single-stage2 or two-stage bioreactor system3. In these engineered ecosystems, a fraction of the ammonium is rst 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 acceptor4,5.

Five 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 environments6. 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 rate7, either biolm reactors that use carrier media or granular sludge reactors are used to retain sufcient biomass in the system8. Such reactor congurations 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)911. 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 ability11. The high aggregate density also serves to limit oxygen diffusion into the granule interior, allowing for the proliferation of anammox and other anaerobic bacteria.

While most studies have focused on understanding the ecophysiology of anammox bacteria4,6,12 and their interactions with autotrophic nitrifying bacteria13, 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 afliated with the phyla Chlorobi, Bacteroidetes, Chloroexi and Proteobacteria comprise a large fraction of the microbial community in anammox bioreactors10,14,15. Despite differences in bioreactor inuent composition, these heterotrophic bacteria appear to share high phylogenetic similarity across different anammox wastewater treatment systems16, suggesting that a more universal interaction exists between them and anammox bacteria.

The 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 denitrication, possibly completing a nitrite loop with anammox and nitrite oxidizing bacteria (NOB) by reducing nitrate back to nitrite16. 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 denitrication and heterotrophic growth, particularly for anammox bioreactors that receive no external organic carbon substrates17,18. However, the specic metabolite exchange reactions promoting interactions between heterotrophic and anammox bacteria remain poorly understood.

Here, 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 organisms metabolic network and served as reference platforms for the proling of gene expression at the community scale. We nd that Chlorobi-afliated 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 ndings improve the understanding of major metabolic activities and interactions occurring between anammox and heterotrophic bacteria and offer the rst transcriptional insights on ecosystem function in anammox granules.

ResultsMetagenomic sequencing and binning. Sequencing of whole community DNA extracted from two separate biomass samples yielded a total of 41,063,036 reads after quality ltering (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 ltered DNA reads obtained from the second sampling event (Supplementary Table 2). Approximately 74% of the quality ltered 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) afliated with the phyla Planctomycetes, Chlorobi, Bacteroidetes, Chloroexi, Proteobacteria and the candidate phylum Microgenomates (OP11; Table 1). Together, these genomes accounted for B82% and 59% of the total quality ltered 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 anammox systems16. Surprisingly, the Chlorobi CHB1 genome was virtually identical to the Chlorobi OLB4 genome recovered by Speth et al.16, sharing an average nucleotide identity of 99.8% despite the signicant differences in reactor operation (two-stage versus single-stage), wastewater inuent composition (potato-processing wastewater versus anaerobic digester ltrate) and geographical location (USA versus Netherlands; Supplementary Data 1).

Microbial community abundance and gene expression. We 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

2 NATURE COMMUNICATIONS | 8:15416 | DOI: 10.1038/ncomms15416 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms15416 ARTICLE

Table 1 | Genome statistics of MAGs recovered from the anammox community.

Bin ID Locus tag

Taxonomy Completeness(%)

Contamination (%)

Genome size (bp)

# Scaffolds

N50 (scaffolds)

GC (%)

Predicted genes

AMX1 UTAMX1 Bacteria; Planctomycetes;Planctomycetia; Planctomycetales;

Planctomycetaceae; Brocadia

96 1 3,142,066 36 137,788 42.3 2,724

AMX2 UTAMX2 Bacteria; Planctomycetes;Planctomycetia; Planctomycetales;

Planctomycetaceae; Brocadia

96 2 3,437,337 148 32,519 45.1 3,064

PLA1 UTPLA1 Bacteria; Planctomycetes 82 2 4,435,946 642 8,518 60.8 4,062 BCD1 UTBCD1 Bacteria; Bacteroidetes;

Sphingobacteriia; Sphingobacteriales

96 1 3,757,364 58 91,406 41.0 3,180

CFX1 UTCFX1 Bacteria; Chloroexi; Anaerolineae 91 2 2,777,515 34 125,597 52.6 2,576 CFX2 UTCFX2 Bacteria; Chloroexi; Anaerolineae 81 5 2,924,761 111 48,865 57.0 2,622 CFX3 UTCFX3 Bacteria; Chloroexi; Anaerolineae 82 1 2,655,948 281 12,640 61.3 2,676 CFX4 UTCFX4 Bacteria; Chloroexi 81 5 4,654,268 672 8,164 54.8 4,526 CFX5 UTCFX5 Bacteria; Chloroexi; Anaerolineae 84 2 3,776,973 364 12,671 62.8 3,432 CHB1 UTCHB1 Bacteria; Chlorobi; Ignavibacteria;

Ignavibacteriales; Melioribacteraceae

95 0 2,424,598 19 211,700 37.6 2,083

CHB2 UTCHB2 Bacteria; Chlorobi; Ignavibacteria;

Ignavibacteriales; Ignavibacteriaceae

95 0 4,062,427 225 36,474 33.0 3,556

CHB3 UTCHB3 Bacteria; Chlorobi; Ignavibacteria;

Ignavibacteriales

83 2 2,943,821 349 10,790 42.1 2,526

CHB4* UTCHB4 Bacteria; Chlorobi; Ignavibacteria;

Ignavibacteriales

43 0 1,189,945 296 4,153 34.7 1,311

CHB5* UTCHB5 Bacteria; Chlorobi; Ignavibacteria;

Ignavibacteriales

6 2 1,684,758 362 4,608 34.4 1,621

PRO1 UTPRO1 Bacteria; Proteobacteria;

Deltaproteobacteria; Myxococcales;

Polyangiaceae

74 2 3,818,019 45 114,904 69.8 3,321

PRO2 UTPRO2 Bacteria; Proteobacteria;

Betaproteobacteria; Rhodocyclales;

Rhodocyclaceae; Sulfuritalea

68 1 2,355,954 380 7,446 66.9 2,730

CPR1 UTCPR1 Bacteria; Microgenomates 68 0 892,243 91 13,172 37.7 959

*indicates incomplete MAGs not submitted to GenBank.

GenBank accession numbers for each MAG can be found in Supplementary Table 3. GC, guanine-cytosine content.

Supplementary Table 2. Overall, gene expression corresponded with abundance in the anammox granules. Genomes afliated 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 B62 and 54%, respectively, whereas CHB1 had a relative abundance and gene expression of

B21 and 26%, respectively. Other organisms that displayed moderate abundance and gene expression in the anammox ecosystem were afliated with the phyla Chlorobi, Bacteroidetes, Chloroexi 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 ecosystems19.

Metabolic reconstruction of microbial community. To 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 (40.5%) using MetaPathways 2.5 (ref. 20). A complete list of the inferred pathways and associated enzymes of each MAG can be found in Supplementary Data 3.

Anammox metatranscriptomic insights. As 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 1020-fold coverage above median gene expression levels in the AMX1 genome (Supplementary Fig. 1). Six other hao-like genes were also identied in the AMX1 genome and displayed above median gene expression levels. Reciprocal best-BLAST searches conrmed 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)5,21, Candidatus Brocadia anammoxidans22 and Candidatus Jettenia caeni23, 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 sinica24 and Candidatus Brocadia fulgida25 (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 rst reduces nitrite to hydroxylamine (instead of nitric oxide), and subsequently converts

NATURE COMMUNICATIONS | 8:15416 | DOI: 10.1038/ncomms15416 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 3

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms15416

CHB2 (MWSV01)

76 I. album (CP003418)

BCD1(MWTC01)

100

CHB4

96

100

CHB5

CHB3 (MWSU01)

Bacteroidetes OLB9 (LNBW01)

PLA1 (MWTD01)

100

M. roseus (CP003557)

98

CHB1 (MWSW01) Chlorobi OLB4 (JYPE01)

Chlorobi

Bacteroidetes

Planctomycetes

Proteobacteria

Chloroflexi

Candidate phylum

100

85

Chlorobi OLB5 (LLZO01)

C. tepidum (AE006470)

Bacteroidetes OLB12 (LNFR01)

Ca. Kueneia stuttgartiensis (CT030148) Candidatus. Scalindua brodae (JRYO01)

100

C. thalassium (CP001100)

52

Chlorobi OLB6 (LLZP01)

100

100

Chlorobi OLB7 (LMYZ01)

Azoarcus sp. PAO1 (LARU01)

100

100 68

N. koreensis (CP003178)

C. pinensis (CP001699)

100

Bacteroidetes OLB11 (LNFQ01)

95

82

100

H. hydrossis (CP002691)

Bacteroidetes OLB10 (LNBX01)

Bacteroidetes OLB8 (JZQW01)

100

100

100

S. canadensis (CP003349)

79

100

Ca. Brocadia fulgida (LAQJ01)

100 Ca. Jettenia caeni (BAFH01)

Candidatus Accumulibacter phosphatis clade IIA (CP001715)

CPR1 (MWSR01)

AMX2 (MWTE01)

100 AMX1 (MWTF01)

Chloroflexi OLB15 (LMZS01)

100

100

Ca. Brocadia sinica JPN1 (BAFN01)

Ca. Brocadia sinica OLB1 (JZEK01)

80

100

100

PRO2 (MWSS01)

95

93 S. hydrogenivorans (AP012547)

A. thermophila (AP012029)

K. racemifer (ADVG01)

100

100

100

100

100

T. denitrificans (CP000116)

55

N. europaea (AL954747)

PRO1 (MWST01)

100 N. soli (CP007035)

35

98

Omnitrophica OLB16 (LMZT01)

CFX1 (MWTB01)

100 Chloroflexi OLB14 (LMZR01)

97

100 CFX3 (MWSZ01)

100

CFX2 (MWTA01)

CFX4 (MWSY01)

100

CFX5 (MWSX01)

47

100

Chloroflexi OLB13 (JZRA01)

100

C. aerophila (AP012337)

100

C. aurantiacus (CP000909)

100

Dehalogenimonas sp. WBC-2 (CP011392)

0.3

Figure 1 | Phylogenetic tree of all recovered draft genomes from the anammox bioreactor. Tree includes MAGs recovered from this study (red) and closely related genomes downloaded from the NCBI genome repository. GenBank accession numbers for each genome are provided in parentheses. Branch node numbers represent bootstrap support values. The tree was constructed using RAxML based on a set of 37 concatenated universal single-copy marker genes.

hydroxylamine and ammonium to hydrazine26 (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 oxide5 or hydroxylamine26, for which an enzyme has not yet been identied. The hao-like genes UTAMX1_1996, UTAMX1_1759 and UTAMX1_1192 expressed in the AMX1 genome that lack a crosslinking tyrosine might full this role (Supplementary Fig. 1).

Denitrication gene expression in anammox granules. All heterotrophic organisms in the anammox granules encoded capabilities for partial or full denitrication (Fig. 3). This was similar to a recent report from a single-stage reactor by Speth et al.16 Genes involved in denitrication 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 Chloroexi genome CFX2 (Fig. 3; Supplementary Data 4). Chloroexi (CFX1 and CFX2)

4 NATURE COMMUNICATIONS | 8:15416 | DOI: 10.1038/ncomms15416 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms15416 ARTICLE

60%

Abundance

Expression

AMX1

CHB1

CHB2

BCD1

CFX1

CFX2

PRO1

PRO2

C1 metabolism

Central metabolism

Alcohol deg.

Carbohydrate deg.

Carboxylate deg.

Nitrogen

Other metabolism

Respiration

Sulfur

Fermentation

0 1 2 3 Relative expression

CBB cycle

50%

Formaldehyde assimilation

Formaldehyde oxidation

Formate oxidation to CO2

Wood-Ljungdahl pathway

Gluconeogenesis

Glycogen biosynthesis

40%

Glycolysis

PPP (non-oxidative branch)

PPP (oxidative branch)

Relative abundance

TCA cycle

Ethanol degradation

30%

Glycerol degradation

Acetoin degradation

Chitin degradation

Galactose degradation

Glycogen degradation

20%

Lactose degradation

Acetate to acetyl-CoA

Citrate degradation

Anammox

Nitrate reduction (denitrification)

10%

DNRA

Pyruvate fermentation to acetate

Glyoxylate cycle

H2 production/oxidation

Superoxide radicals degradation

2-Oxoglutarate decarboxylation

0%

2-Oxoisovalerate decarboxylation

AMX1

CHB1

CHB4

CHB5

CHB2

BCD1

CFX1

CFX2

PRO1

PRO2

AMX1

CFX3

CFX5

PLA1

CPR1

CFX4

CHB3

Aerobic respiration (cytochrome c)

Sulfate reduction (dissimilatory)

Sulfide oxidation

Sulfite oxidation

Figure 2 | Abundance and gene expression of organisms represented by the MAGs recovered from the anammox bioreactor. Abundance and gene expression estimates were based on RPKM values of metagenomic reads and transcripts that mapped to each MAG, respectively. See Supplementary Table 2 for mapping details.

nrfHA

Figure 4 | Relative gene expression of major carbon and energy metabolic pathways encoded by each MAG. Colour intensity represents gene expression, based on mapping of metatranscriptomic reads to the metagenomic assembly. Gene expression was relativized by median RPKM values calculated across each ORF in a given MAG (Methods section).

A value of 1 equals median expression in a given genome. Yellow box indicates pathway absence. Anammox metabolic pathway in the AMX1 genome had a relative expression value of 15. See Supplementary Data 3 for a detailed summary of all reconstructed metabolic pathways.

Oxidation state

+5 +3 +2 +1 0 3

narGHIJ

NO3

nirK norBC nosZ norZ

N2

NO2 N2O

NO

NH3

0 2 4 5+

1 3

CHB1

CHB2

CFX1

CFX2

BCD1

PRO1

PRO2

napAB nirS

RPKM (log2)

Figure 3 | Presence and expression of denitrication genes across the recovered heterotrophic genomes. Purple arrows indicate gene presence. Colour intensity represents gene expression (log2 RPKM), based on mapping of metatranscriptomic reads to the metagenomic assembly. An outlined white arrow indicates one or more enzyme subunits had no detectable gene expression. A summary of genes involved in denitrication across the recovered genomes can be found in Supplementary Data 4.

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 denitrication, thus enhancing overall nitrogen removal performance in the bioreactor.

Carbon and energy metabolism of heterotrophic bacteria. Figure 4 shows the gene expression proles 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 xation that would permit mixotrophic growth (Fig. 4; Supplementary Data 3). BCD1 displayed above median gene expression of PEP carboxylase (UTBCD1_1894), which allows CO2 xation into oxaloacetate (Supplementary Data 2); PRO2 displayed above median gene expression of Ribulose-1,5-bisphosphate carboxylase (RuBisCO, UTPRO2_2197), which would permit CO2 xation via the Calvin-Benson-Bassham cycle (Fig. 4; Supplementary Data 3). While all heterotrophs, except CHB1, expressed genes encoding pyruvate:avodoxin/ferredoxin oxidoreductase (PFOR) and

NATURE COMMUNICATIONS | 8:15416 | DOI: 10.1038/ncomms15416 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 5

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms15416

a

MEROPS

AMX1

CHB1

M01

M06

M12B

M23A

M23B

S01C

S08A

S09C

CHB2

BCD1

CFX1

CFX2

PRO1

PRO2

AMX1

CHB1

CHB2

BCD1

CFX1

CFX2

PRO1

PRO2

AMX1

CHB1

CHB2

BCD1

CFX1

CFX2

PRO1

PRO2

b

Legend

TCDB

Cationic amino acid transporter (CAT)

1

2

3

4

5

6

7

0

1

2

3

4

5

Archaeal/Bacterial transporter (ABT)

Rel. expression

L-type amino acid transporter (LAT)

Count

Proton-dependent oligopeptide transporter (POT)

Oligopeptide transporter (OPT)

Polar amino acid uptake transporter (PAAT)

Hydrophobic amino acid uptake transporter (HAAT)

Peptide/Opine/Nickel uptake transporter (PepT)

AMX1

CHB1

CHB2

BCD1

CFX1

CFX2

PRO1

PRO2

Figure 5 | Predicted peptidases and amino acid transporters recovered from the MAGs. (a) Number (bubble diameter) and relative gene expression (bubble colour intensity) of selected peptidases possibly involved in EPS matrix protein degradation. Peptidases were annotated against the MEROPS database66. The subcellular location (extracellular, outer membrane or peroplasm) of each peptidase was predicted using the subcellular localization predictor (CELLO)69. A summary of all predicted peptidases can be found in Supplementary Data 5. (b) Number (bubble diameter) and relative gene expression (bubble colour intensity) of amino acid and peptide transporters predicted across the recovered genomes. Transporters were annotated against the transporter classication database (TCDB)70 and can be found in Supplementary Data 2.

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 cycle27), it is likely that these genes are not used for CO2 xation 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 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 conrmed to be missing from the closely related Chlorobi bacterium OLB4 and OLB6 genomes16, although they were present in the genomes of Ignavibacterium album28 and Melioribacter roseus29. An alternative route for NADPH generation in CHB1 and CHB2 is likely through NAD(P) transhydrogenase30, which transfers electrons from

NADH to NADP and had above median gene expression in both genomes (Supplementary Data 2).

Genes 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 afnity pathway for acetate assimilation compared to the high-afnity acetyl-CoA synthetase (acs) pathway31. 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 oxidore

ductase (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 PFOR32. This energy conservation mode was also identied in the closely related I. album28 and M. roseus29 genomes, however no orthologs were detected in the CHB1 genome.

The ability to use H2, sulte, 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 production33. Organisms encoding the

6 NATURE COMMUNICATIONS | 8:15416 | DOI: 10.1038/ncomms15416 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms15416 ARTICLE

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 generation34. 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 rst step in CO2 xation 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-afnity aa3-type terminal cytochrome c oxidase; CHB2 expressed a high-afnity 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 oxygen35, 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 EPS produced by anammox bacteria has been proposed to be a major organic carbon source supporting heterotrophic growth in anammox granules11,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 substrates36,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 cost38 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.

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

Several 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

AMX1

CHB1

CHB2

BCD1

CFX1

CFX2

PRO1

PRO2

AMX1

CHB1

CHB2

BCD1

CFX1

CFX2

PRO1

PRO2

Aspartate (17)

Glutamate (15)

Arginine (9)

Histidine (4)

Lysine (7)

Asparagine (16)

Glutamine (14)

Serine (20)

Threonine (13)

Alanine (18)

Isoleucine (6)

Leucine (8)

Methionine (5)

Phenylalanine (2)

Tryptophan (1)

Tyrosine (3)

Valine (11)

Cysteine (10)

Glycine (19)

Proline (12)

Biosynthesis Degradation

0 1 2 3

Relative expression

Figure 6 | Relative gene expression of amino acid biosynthetic and degradation pathways encoded by each MAG. Red colour intensity represents gene expression, which was relativized by median RPKM values calculated across each ORF in a given MAG. A value of 1 equals median expression. Top, middle and bottom panels separate amino acids by hydrophilic, hydrophobic and special structured side chains, respectively. Bracketed numbers rank the metabolic cost of amino acid biosynthesis based on values reported by Akashi and Gojobori38, with 1 being the most costly.

NATURE COMMUNICATIONS | 8:15416 | DOI: 10.1038/ncomms15416 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 7

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms15416

N2 NH4+ NO2 N O

N

MQH

H+ H+ + +

Electron transport

CO2

Brocadia(AMX1)

Chlorobi (CHB1)

Figure 7 | Proposed metabolic interactions between Brocadia (AMX1) and Chlorobi (CHB1) in anammox granules. AMX1 xes CO2 and synthesizes amino acids and EPS. CHB1 degrades proteins bound in EPS using extracellular peptidases and subsequently transports and catabolizes short peptides (circles) to central carbon intermediates. Nitrite oxidation and reduction by AMX1 and CHB1, respectively, results in a distributed nitrite loop. Purple arrows indicate nitrogen cycling; orange arrows indicate carbon cycling; light blue arrows indicate vitamin B metabolite exchange. Hatched ovals indicate peptidases. BCAA, branched-chain amino acids. Met(?), methionine synthase present in genome but other steps involved in Met biosynthesis not identied. The presence of a periplasm has been ignored for clarity of the schematic.

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.

DiscussionThe 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 afliated with Brocadia (AMX1) and heterotrophic bacteria afliated with Chlorobi (CHB1) dominated the abundance and gene expression of the anammox granule community. Other less abundant heterotrophic bacteria afliated with Chloroexi, Bacteroidetes and Proteobacteria also displayed signicant gene expression that may contribute to ecosystem function. The high similarity between these genomes and the genomes recently recovered from a single-stage anammox reactor16 suggests that a core microbiome exists in anammox-based wastewater treatment systems, despite differences in reactor operation and inuent wastewater composition.

Remarkably, the CHB1 genome was nearly identical to the Chlorobi bacterium OLB4 genome recovered by Speth et al16(99.8% average nucleotide identity) that achieved the second highest abundance in a full-scale partial-nitritation-anammox reactor from the Netherlands16,40, 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 carbon41.

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 Brocadia sp., which derives carbon via CO2 xation 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

8 NATURE COMMUNICATIONS | 8:15416 | DOI: 10.1038/ncomms15416 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms15416 ARTICLE

Table 2 | Process performance data from the two-stage PNA bioreactor at steady-state operation.

pH COD (mg l 1) Ammonia (mgN l 1) Nitrite (mgN l 1) Nitrate (mgN l 1) VSS (mg l 1) Activity(kg NH4-N kg 1 VSS-day)

Filtrate 8.1 53510.8 716.0 69.6 20.8 NA NA PN-efuent 6.8 40512.2 462.0 424.8 46.2 NA NA Anammox-efuent 6.9 35014.7 56.0 28.0 48.2 NA NA Anammox reactor 6.8-8.0 NA NA NA NA 802 0.3090.0094 PN reactor 6.5-8.5 NA NA NA NA 643 0.4190.0265

COD, chemical oxygen demand; PN, partial nitritation; VSS, volatile suspended solids.

carbon compounds degraded and utilized cellular components produced by anammox bacteria42. The high transcript abundance of CHB1 in the metatranscriptome (25% relative community expression) also suggests that proteolysis by this organism is a signicant component of carbon ux 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 signicant contributor to community assembly in anammox granules, similar to particle degrading taxa in marine ecosystems43.

It 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 gut44. The specic role vitamin exchange plays in shaping anammox 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 loop with Brocadia 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 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 granules45. 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.

In addition to a nitrite loop, heterotrophic bacteria in anammox ecosystems may facilitate complete denitrication, either alone or based on metabolite exchange of denitrication intermediates with other organisms. Indeed, the segregation of denitrication intermediates across different organisms has been shown to reduce the accumulation of growth-inhibiting metabolites (for example, nitrite) and potentially accelerate denitrication rates46. The close spatial association of organisms in anammox granules may encourage such a distributed metabolic network47,48. While complete denitrication 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. Specically, 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.

Aside 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 metabolism49. 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 signicance of H2 and acetate metabolism in anammox granules still remains to be determined.

In 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, Chloroexi and Proteobacteria) and also sheds light on their potential ecological roles and interactions. Our results implicate members of the Chlorobi afliated 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.

While further studies are still required to determine the specic 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.

NATURE COMMUNICATIONS | 8:15416 | DOI: 10.1038/ncomms15416 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 9

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms15416

Methods

Bioreactor operation. A 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 performance50. 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 ltrate obtained from the belt lter press of a local wastewater treatment plant (CVWRF, Salt Lake City, UT, USA). Before entering the anammox reactor, the ltrate was rst passed through a nitritation reactor, where approximately half of the inuent ammonia was oxidized to nitrite. The nitritation reactor had a working volume of 2 l with a hydraulic retention time of 1 day50. Dissolved oxygen in the reactor was maintained at 0.50.5 mg l 1 and bicarbonate present in the inuent ltrate 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 dioxide50. Table 2 summarizes the process performance data of the nitritation and anammox bioreactors during steady-state operation.

DNA and RNA sequencing. Two independent biomass samples were collected from the anammox bioreactor at different time points during periods of high total nitrogen removal efciency (480%). 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 ash 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 Scientic, 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 specic 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 rst strand cDNA synthesis reaction, as well as in the nal library construction steps using manufacturers protocols. DNA from the rst 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.

Metagenomic assembly and binning. Raw paired-end reads from the MiSeq and HiSeq platforms were initially ltered 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/

Web End =https://sourceforge.net/projects/bbmap/) with the parameters minid 0.95 and

ambig random. Resulting mapping les were subsequently used by MetaBAT

v0.26.3 to bin metagenomic contigs into draft genomes, based on the sensitive parameters53. 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 genes55. 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 genes56. The taxonomic afliation of each draft genome was determined using the Phylosift all command. Marker genes were also identied 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). Briey, 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 annotation58. Contig nucleotide sequences were also queried against theSILVA 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 database59. 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 pathway60. All inferred pathways were then manually curated to verify predictions made by Pathologic.

Metatranscriptomic analysis. Metatranscriptomic reads were quality ltered and merged as described above. Subsequently, rRNA sequences were ltered from the metatranscriptomic data set using SortMeRNA v2.0, based on multiple rRNA databases for bacterial, archaeal and eukaryotic sequences61. Resulting non-rRNA reads were mapped back to all assembled metagenomic contigs using BBMap v35.92 with the parameters minid 0.95, which species a minimum alignment

identify of 95% corresponding to the well-established criteria for identifying microbial species using average nucleotide identify62 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 parameter63 and normalized by sequencing depth and ORF length and expressed as RPKM values64. This allowed for the proling 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.

The 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).

Carbohydrate hydrolase and peptidase identication. Carbohydrate hydrolases and peptidases potentially involved in the breakdown of EPS were identied 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 conrmedbased on queries against the NCBIs Conserved Domain Database67 and Pfam database68 using MOTIF (http://www.genome.jp/tools/motif/

Web End =http://www.genome.jp/tools/motif/). The subcellular location of identied proteins was predicted using CELLO v2.5, which classies proteins using support vector machines trained by multiple feature vectorsbased on n-peptide composition69.

Data availability statement. All 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 les for the whole genome sequences (MAGs) described in this study can be found under accession numbers listed in Supplementary Table 3.

References

1. Kartal, B., Kuenen, J. G. & van Loosdrecht, M. C. M. Sewage treatment with anammox. Science 328, 702703 (2010).

2. Sliekers, A. O. et al. Completely autotrophic nitrogen removal over nitrite in one single reactor. Water Res. 36, 24752482 (2002).

3. Van Dongen, U., Jetten, M. S. M. & Van Loosdrecht, M. C. M. The SHARON-Anammox process for treatment of ammonium rich wastewater. Water Sci. Technol. 44, 153160 (2001).

4. Strous, M. et al. Deciphering the evolution and metabolism of an anammox bacterium from a community genome. Nature 440, 790794 (2006).

5. Kartal, B. et al. Molecular mechanism of anaerobic ammonium oxidation. Nature 479, 127130 (2011).

6. Kartal, B. et al. How to make a living from anaerobic ammonium oxidation. FEMS Microbiol. Rev. 37, 428461 (2013).

7. Strous, M. et al. Missing lithotroph identied as new planctomycete. Nature 400, 446449 (1999).

10 NATURE COMMUNICATIONS | 8:15416 | DOI: 10.1038/ncomms15416 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms15416 ARTICLE

8. Abma, W. R. et al. Full-scale granular sludge Anammox process. Water Sci. Technol. 55, 2733 (2007).

9. Ni, B. J. et al. Microbial and physicochemical characteristics of compact anaerobic ammonium-oxidizing granules in an upow anaerobic sludge blanket reactor. Appl. Environ. Microbiol. 76, 26522656 (2010).

10. Gonzalez-Gil, G., Sougrat, R., Behzad, A. R., Lens, P. N. L. & Saikaly, P. E. Microbial community composition and ultrastructure of granules froma full-scale anammox reactor. Microb. Ecol. 70, 118131 (2014).

11. Hou, X., Liu, S. & Zhang, Z. Role of extracellular polymeric substance in determining the high aggregation ability of anammox sludge. Water Res. 75, 5162 (2015).

12. Ali, M. et al. Physiological characterization of anaerobic ammonium oxidizing bacterium Candidatus Jettenia caeni. Environ. Microbiol. 17, 21722189 (2015).

13. Ma, Y., Sundar, S., Park, H. & Chandran, K. The effect of inorganic carbon on microbial interactions in a biolm nitritation-anammox process. Water Res. 70, 246254 (2015).

14. Li, X. R. et al. The bacterial diversity in an anaerobic ammonium-oxidizing (anammox) reactor community. Syst. Appl. Microbiol. 32, 278289

2009:

15. Gonzalez-Martinez, A. et al. Microbial community analysis of a full-scale DEMON bioreactor. Bioprocess Biosyst. Eng. 38, 499508 (2015).

16. Speth, D. R., in t Zandt, M. H., Guerrero-Cruz, S., Dutilh, B. E. & Jetten, M. S. M. Genome-based microbial ecology of anammox granules in a full-scale wastewater treatment system. Nat. Commun. 7, 11172 (2016).

17. Ni, B. J., Ruscalleda, M. & Smets, B. F. Evaluation on the microbial interactions of anaerobic ammonium oxidizers and heterotrophs in Anammox biolm. Water Res. 46, 46454652 (2012).

18. Liu, Y. et al. Assessment of heterotrophic growth supported by soluble microbial products in anammox biolm using multidimensional modeling. Sci. Rep. 6, 27576 (2016).

19. Lawson, C. E. et al. Rare taxa have potential to make metabolic contributions in enhanced biological phosphorus removal ecosystems. Environ. Microbiol. 17, 49794993 (2015).

20. Konwar, K. M. et al. MetaPathways v2.5: quantitative functional, taxonomic and usability improvements. Bioinformatics 31, 33453347

2015:

21. Maalcke, W. J. et al. Structural basis of biological no generation by octaheme oxidoreductases. J. Biol. Chem. 289, 12281242 (2014).

22. Schalk, J., De Vries, S., Kuenen, J. G. & Jetten, M. S. M. Involvement of a novel hydroxylamine oxidoreductase in anaerobic ammonium oxidation. Biochemistry 39, 54055412 (2000).

23. Shimamura, M. et al. Another multiheme protein, hydroxylamine oxidoreductase, abundantly produced in an anammox bacterium besides the hydrazine-oxidizing enzyme. J. Biosci. Bioeng. 105, 243248

2008:

24. Oshiki, M., Shinyako-Hata, K., Satoh, H. & Okabe, S. Draft genome sequence of an anaerobic ammonium-oxidizing bacterium,Candidatus Brocadia sinica. Genome Announc. 3, 34 (2015).

25. Ferousi, C. et al. Identication of the type II cytochrome c maturation pathway in anammox bacteria by comparative genomics. BMC Microbiol. 13, 265 (2013).

26. Oshiki, M., Ali, M., Shinyako-Hata, K., Satoh, H. & Okabe, S. Hydroxylamine-dependent anaerobic ammonium oxidation (anammox) by Candidatus Brocadia sinica. Environ. Microbiol. 18, 31333143 (2016).

27. Zhang, S. & Bryant, D. The tricarboxylic acid cycle in cyanobacteria. Science 7002, 15511554 (2011).

28. Liu, Z. et al. Complete genome of Ignavibacterium album, a metabolically versatile, agellated, facultative anaerobe from the phylum Chlorobi. Front. Microbiol. 3, 115 (2012).

29. Kadnikov, V. V. et al. Genomic analysis of Melioribacter roseus, facultatively anaerobic organotrophic bacterium representing a novel deep lineage within Bacteriodetes/Chlorobi group. PLoS ONE 8, e53047 (2013).

30. Spaans, S. K., Weusthuis, R. A., van der Oost, J. & Kengen, S. W. M. NADPH-generating systems in bacteria and archaea. Front. Microbiol. 6, 127 (2015).

31. Wolfe, A. J. The acetate switch. Microbiol. Mol. Biol. Rev. 69, 1250 (2005).32. Buckel, W. & Thauer, R. K. Energy conservation via electron bifurcating ferredoxin reduction and proton/Na translocating ferredoxin oxidation.

Biochim. Biophys. Acta 1827, 94113 (2013).33. McDowall, J. S. et al. Bacterial formate hydrogenlyase complex. Proc. Natl Acad. Sci. USA 111, E3948E3956 (2014).

34. Sawers, G. The hydrogenases and formate dehydrogenases of Escherichia coli. Antonie van Leeuwenhoek 66, 5788 (1994).

35. Richardson, D. J. Bacterial respiration: a exibe process for a changing environment. Microbiology 146, 551571 (2000).

36. Hase, C. C. & Finkelstein, R. A. Bacterial extracellular zinc-containing metalloproteases. Microbiol. Rev. 57, 823837 (1993).

37. Siezen, R. J. & Leunissen, J. A. Subtilases: the superfamily of subtilisin-like serine proteases. Protein Sci. 6, 501523 (1997).

38. Akashi, H. & Gojobori, T. Metabolic efciency and amino acid composition in the proteomes of Escherichia coli and Bacillus subtilis. Proc. Natl Acad. Sci. USA. 99, 36953700 (2002).

39. Mee, M. T., Collins, J. J., Church, G. M. & Wang, H. H. Syntrophic exchange in synthetic microbial communities. Proc. Natl Acad. Sci. USA 111, E2149E2156 (2014).

40. Abma, W. R., Driessen, W., Haarhuis, R. & Van Loosdrecht, M. C. M. Upgrading of sewage treatment plant by sustainable and cost-effective separate treatment of industrial wastewater. Water Sci. Technol. 61, 17151722 (2010).

41. Van Mooy, B. A. S., Keil, R. G. & Devol, A. H. Impact of suboxia on sinking particulate organic carbon: Enhanced carbon ux and preferential degradation of amino acids via denitrication. Geochim. Cosmochim. Acta 66, 457465 (2002).

42. Kindaichi, T., Yuri, S., Ozaki, N. & Ohashi, A. Ecophysiological role and function of uncultured Chloroexi in an anammox reactor. Water Sci. Technol. 66, 25562561 (2012).

43. Datta, M. S., Sliwerska, E., Gore, J., Polz, M. & Cordero, O. X. Microbial interactions lead to rapid micro-scale successions on model marine particles. Nat. Commun. 7, 11965 (2016).

44. Degnan, P. H., Taga, M. E. & Goodman, A. L. Vitamin B12 as a modulator of gut microbial ecology. Cell Metab. 20, 769778 (2014).

45. Ali, M. et al. Source identication of nitrous oxide emission pathways from a single-stage nitritation-anammox granular reactor. Water Res. 102, 147157 (2016).

46. Lilja, E. E. & Johnson, D. R. Segregating metabolic processes into different microbial cells accelerates the consumption of inhibitory substrates. ISME J. 10, 15681578 (2016).

47. Kindaichi, T. et al. In situ activity and spatial organization of anaerobic ammonium-oxidizing (anammox) bacteria in biolms. Appl. Environ. Microbiol. 73, 49314939 (2007).

48. Vlaeminck, S. E. et al. Aggregate size and architecture determine microbial activity balance for one-stage partial nitritation and anammox. Appl. Environ. Microbiol. 76, 900909 (2010).

49. Kartal, B. et al. Candidatus Brocadia fulgida: An autouorescent anaerobic ammonium oxidizing bacterium. FEMS Microbiol. Ecol. 63, 4655 (2008).

50. Kotay, S. M., Mansell, B. L., Hogsett, M., Pei, H. & Goel, R. Anaerobic ammonia oxidation (ANAMMOX) for side-stream treatment of anaerobic digester ltrate process performance and microbiology. Biotechnol. Bioeng. 110, 11801192 (2013).

51. Joshi, N. A. & Fass, J. N. Sickle: A sliding-window, adaptive, quality-based trimming tool for FastQ les Available at https://github.com/najoshi/sickle (2011).

52. Salzberg, S. L. FLASH: fast length adjustment of short reads to improve genome assemblies tanja mago. Bioinformatics 27, 29572963

2011:

53. Kang, D. D., Froula, J., Egan, R. & Wang, Z. MetaBAT, an efcient tool for accurately reconstructing single genomes from complex microbial communities. PeerJ 3, e1165 (2015).

54. Parks, D. H., Imelfort, M., Skennerton, C. T., Hugenholtz, P. & Tyson, G. W. CheckM: assessing the quality of microbial genomes recovered fromisolates, single cells, and metagenomes. Genome Res. 25, 10431055

2015:

55. Dupont, C. L. et al. Genomic insights to SAR86, an abundant and uncultivated marine bacterial lineage. ISME J. 6, 11861199 (2012).

56. Darling, A. E. et al. PhyloSift: phylogenetic analysis of genomes and metagenomes. PeerJ 2, e243 (2014).

57. Hyatt, D. et al. Prodigal: prokaryotic gene recognition and translation initiation site identication. BMC Bioinformatics 11, 119 (2010).

58. Konwar, K., Hanson, N., Page, A. & Hallam, S. J. MetaPathways: a modular pipeline for constructing Pathway/Genome Databases from environmental sequence information. BMC Bioinformatics 14, 202 (2013).

59. Karp, P. D. et al. Pathway Tools version 19.0 update: software for pathway/ genome informatics and systems biology. Brief. Bioinform. 17, 877890 (2015).

60. Karp, P. D., Latendresse, M. & Caspi, R. The pathway tools pathway prediction algorithm. Stand. Genomic Sci. 5, 424429 (2011).

61. Kopylova, E., No, L. & Touzet, H. SortMeRNA: Fast and accurate ltering of ribosomal RNAs in metatranscriptomic data. Bioinformatics 28, 32113217 (2012).

62. Konstantinidis, K. T. & Tiedje, J. M. Genomic insights that advance the species denition for prokaryotes. Proc. Natl. Acad. Sci. USA 102, 25672572 (2005).

63. Anders, S., Pyl, P. T. & Huber, W. HTSeq-A Python framework to work with high-throughput sequencing data. Bioinformatics 31, 166169 (2015).

64. Mortazavi, A., Williams, B. A., McCue, K., Schaeffer, L. & Wold, B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat. Methods 5, 621628 (2008).

NATURE COMMUNICATIONS | 8:15416 | DOI: 10.1038/ncomms15416 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 11

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms15416

65. Lombard, V., Golaconda Ramulu, H., Drula, E., Coutinho, P. M. & Henrissat, B. The carbohydrate-active enzymes database (CAZy) in 2013. Nucleic Acids Res. 42, 490495 (2014).

66. Rawlings, N. D., Barrett, A. J. & Finn, R. Twenty years of the MEROPS database of proteolytic enzymes, their substrates and inhibitors. Nucleic Acids Res. 44, D343D350 (2015).

67. Marchler-Bauer, A. et al. CDD: NCBIs conserved domain database. Nucleic Acids Res. 43, D222D226 (2015).

68. Finn, R. D. et al. The Pfam protein families database: Towards a more sustainable future. Nucleic Acids Res. 44, D279D285 (2016).

69. Yu, C.-S., Lin, C.-J. & Hwang, J.-K. Predicting subcellular localization of proteins for Gram-negative bacteria by support vector machines based on n-peptide compositions. Protein Sci. 13, 14021406 (2004).

70. Saier, M. H. et al. The transporter classication database (TCDB): Recent advances. Nucleic Acids Res. 44, D372D379 (2016).

Acknowledgements

This work was supported by funding from the Central Valley Water Reclamation Facility (CVWRF) and the North Davis Sewer District awarded to R.G. and S.W.; Partial support was also provided by funding from the National Science Foundation (CBET-1435661 and MCB-1518130) awarded to D.R.N. and K.D.M., and a Postgraduate Scholarship-Doctoral (PGS-D) awarded by the National Sciences and Engineering Research Council of Canada (NSERC) to C.E.L.; The authors acknowledge Masaru K. Nobu for helpful discussion with the manuscript.

Author contributions

C.E.L., K.D.M., R.G. and D.R.N. designed the study. S.W. and A.S.B. operated the bioreactor and performed the sampling and sequencing. C.E.L. analysed the data

and wrote the manuscript. J.J.H. and A.S.B. contributed to the data analysis. J.J.H., A.S.B., K.D.M., R.G. and D.R.N. reviewed and provided valuable edits to the manuscript.

Additional information

Supplementary Information accompanies this paper at http://www.nature.com/naturecommunications

Web End =http://www.nature.com/ http://www.nature.com/naturecommunications

Web End =naturecommunications

Competing interests: The authors declare no competing nancial interests.

Reprints and permission information is available online at http://npg.nature.com/reprintsandpermissions/

Web End =http://npg.nature.com/ http://npg.nature.com/reprintsandpermissions/

Web End =reprintsandpermissions/

How 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).

Publishers note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional afliations.

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the articles Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

Web End =http://creativecommons.org/licenses/by/4.0/

r The Author(s) 2017

12 NATURE COMMUNICATIONS | 8:15416 | DOI: 10.1038/ncomms15416 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications