甘肃农业大学：食品科学与工程学院（文献讲义）Discrepant gene functional potential and cross-feedings of anammox bacteria Ca. Jettenia caeni and Ca. Brocadia sinica in response to acetate.pdf

Discrepant gene functional potential and cross-feedings of anammox bacteria Ca. Jettenia caeni and Ca. Brocadia sinica in response to acetate Ying Feng a, b, 1 , Yunpeng Zhao a, b, 1 , Bo Jiang a, b , Huazhang Zhao a, b , Qilin Wang c , Sitong Liu a, b, * a Department of Environmental Engineering, Peking University, Beijing, 100871, China b Key Laboratory of Water and Sediment Sciences, Ministry of Education of China, Beijing, 100871, China c Centre for Technology in Water and Wastewater, School of Civil and Environmental Engineering, University of Technology Sydney, Ultimo, NSW, 2007, Australia article info Article history: Received 19 March 2019 Received in revised form 9 August 2019 Accepted 10 August 2019 Available online 10 August 2019 Keywords: Anammox Acetate Cross-feedings Metagenomics Metatranscriptomics abstract Although the enhancement of anammox performance for wastewater treatment due to the addition of small amount of acetate has been reported, discrepant metabolic responses of different anammox species have not been experimentally evaluated. Based on metagenomics and metatranscriptomic data, we investigated the competitiveness between two typical anammox species, Candidatus Jettenia caeni (J. caeni) and Candidatus Brocadia sinica (B. sinica), in anammox consortia under mixotrophic condition, where complex metabolic interactions among anammox bacteria and heterotrophs also changed with acetate addition. Contrary to J. caeni, the dissimilatory nitrate reduction to ammonium pathway of B. sinica was markedly stimulated for improving nitrogen removal. More acetate metabolic pathways and up-regulated AMP-acs expression for acetyl-CoA synthesis in B. sinica contributed to its superiority in acetate utilization. Interestingly, cross-feedings, including the nitrogen cycle, amino acid cross-feeding and B-vitamin metabolic exchange between B. sinica and other heterotrophs seemed to be enhanced with acetate addition, contributing to a reduction in metabolic energy cost to the whole community. Our work not only clarified the mechanism underlying discrepant responses of different anammox species to acetate, but also suggests a possible strategy for obtaining higher nitrogen removal rates in wastewater treatment under low C/N ratio. © 2019 Elsevier Ltd. All rights reserved. 1. Introduction Depending on nutritional conditions, microorganisms may grow autotrophically, heterotrophically or mixotrophically (Smith et al., 1980). The growth rate and activity, as well as metabolic function, of microorganisms may vary with nutrient conditions. Generally, heterotrophs and mixotrophs grow at higher rates than autotrophs, as they are able to obtain more energy from organic substrates due to differences in metabolism (Kim et al., 2013). More importantly, intracellular metabolic pathways, as well as activities and gene functional potential, could be dissimilar even among different species within the same genus. Analysis based on comparative genomics have indicated that closely related strains may exhibit metabolic divergence due to genomic discrepancies (Bombar et al., 2014). Wastewater is a complex mixture consisting of a wide range of organic matter, which inevitably affects the metabolic state of microorganisms in sludge during the treatment process (Le and Stuckey, 2016). Anaerobic ammonia oxidation (anammox), which is regarded as an energy-efficient process for wastewater treatment, has drawn a lot of attention recently (Kartal et al., 2010). Anammox bacteria are able to use, not only inorganic carbon sources, but also organic matter such as acetate and propionate for reducing inorganic nitrates/nitrate ions (NO3 ) to ammonium cations (NH4 þ) via the dissimilatory nitrate reduction to ammonium (DNRA) pathway (Guven et al., 2005; Kartal et al., 2007a, 2007b; Van De Vossenberg et al., 2008). According to previous studies, organics may have an effect on the performance of the anammox reactor (Tang et al., 2014) as well as on the functional gene expression profile (Shu et al., 2015) and the microbial community * Corresponding author. College of Environmental Science and Engineering, Peking University, Yiheyuan Road, No.5, Haidian District, Beijing, 100871, China. E-mail address: liusitong@pku.edu.cn (S. Liu). 1 These authors contributed equally to this work. Contents lists available at ScienceDirect Water Research journal homepage: www.elsevier.com/locate/watres https://doi.org/10.1016/j.watres.2019.114974 0043-1354/© 2019 Elsevier Ltd. All rights reserved. Water Research 165 (2019) 114974

structure (Leal et al., 2016). Responses of different anammox species to organics seem to be distinct. For example, it was reported that Candidatus Jettenia asiatica (J. aiatica) showed no superiority in growth under mixotrophic conditions compared to autotrophs (Huang et al., 2014), while other study found that the biomass of Candidatus Brocadia fulgida (B. fulgida) showed an increase under certain C/N ratios (Jenni et al., 2014). However, the process by which, organics affect the gene functional potential of different anammox species still remains unclear (Shu et al., 2015). As anammox bacteria grow extremely slowly, it may be meaningful to investigate the mixotrophic metabolism of anammox bacteria further, as they may be acquiring extra energy from organic matter in order to grow faster and thus shorten the reactor start-up period (Kartal et al., 2012). Nutrient sources also affect metabolic interactions of the microbiome. Metabolic interactions are ubiquitous in microbial communities, especially in microscale cell aggregates, which play an important role in the functioning of microbial communities (Cordero and Datta, 2016). From an ecological aspect, metabolic interactions may help to maintain a stable coexistence between bacteria as a strategy to decrease the energy consumption of the community (Guo et al., 2018; Pande et al., 2014). The emergence and maintenance of metabolic interactions depends on many factors, such as nutrient sources (Benomar et al., 2015) and spatial organization (Jiang et al., 2018). Variation in nutrient conditions may influence gene transcription and thereby impact metabolic interaction (Steffen et al., 2014). With the rapid development of meta-omics technology, the subject of metabolic interactions in microbial communities has drawn wide attention and turned into an important topic (Ponomarova and Patil, 2015). Pure anammox culture is extremely hard to obtain (Kuenen, 2008), and many heterotrophs, such as Chloroflexi and Chlorobi, are abundant in these communities (Speth et al., 2016). Recently, metabolic interactions which could perform energy-efficient nitrogen removal from wastewater, such as degradation of extracellular peptide substrates of anammox bacteria by heterotrophs and nitrogen and metabolite cross-feeding between anammox bacteria and heterotrophs, have been found in anammox consortia (Lawson et al., 2017). Cross-feeding is a kind of microbial interaction, in which metabolites could be shared by both the producer and the receiver, thus they can benefit from this process (Zhao et al., 2018). An investigative study of the mechanism underlying metabolic interactions in microbial communities may broaden our insight in regard to the composition and assembly of these communities (Zengler and Zaramela, 2018). However, the process by which organics influence interactions between anammox bacteria and heterotrophs in the consortia still remain unresolved. We investigated competition between two typical anammox species, J. caeni and B. sinica, under autotrophic condition and mixotrophic condition with acetate addition based on batch tests. This phenotype was mapped to the underlying microbiome and further determined by sampling and analyzing autotrophic and mixotrophic anammox consortia. We characterized the gene functional potential, as well as the metabolic network, using levels of anammox species in the anammox community, in order to explore the hypothesis that different anammox species have discrepant responses to acetate addition. As well, potential mechanisms associated with individual anammox species and their metabolic interactions in the consortia were analyzed. Our study provides a novel, detailed insight into mixotrophic metabolism of anammox bacteria, and suggests the possibility of predicting an increase in anammox performance under low C/N ratio. 2. Materials and methods 2.1. Sample collection The anammox consortia used in this study was collected from a 3 L lab-scale sequencing batch reactor (SBR), which had been in operation at 37 C for 280 days (Tang et al., 2018a). It was fed with a synthetic medium solution (Van de Graaf et al., 1995), and the concentration of NH4 þ-N and NO2 -N in the influent were 300 mg L1 . The hydraulic retention time (HRT) was 0.75 d. The pH was maintained at 6.8e7.5, and dissolved oxygen was removed by sparging with N2-CO2 (95/5%) gas. The batch tests were performed in six 250 mL serum bottles with the effective volume of 200 mL, each containing 0.315 g volatile suspended solids (VSS)/L anammox consortia inoculum mentioned above. By referring to the previous study about heterotrophic metabolism of anammox bacteria (Güven et al., 2005), the initial concentrations of NH4 þ-N and NO2 -N in synthetic medium solution were set at a ratio of 1:1 with the concentrations of 50 mg L1 . No sodium acetate was added (COD/ TN ¼ 0) to the three serum bottles of the control group (autotrophic group), while acetate was added to the three serum bottles of the experiment group (mixotrophic group), to maintain a final COD/TN ratio of 0.3 as per a previous study (Feng et al., 2018). The pH of the medium was adjusted to 7.2 by adding 0.1 M NaOH solution (Carvajal-Arroyo et al., 2014). Although the pH was not controlled, it was within a constant range of 7.2e7.5 during the experiment. After being sprayed with a gas mixture of N2-CO2 (95/5%) in order to maintain strict anaerobic conditions, all serum bottles were incubated at 37 C and agitated at 150 rpm in the dark. 1 mL supernatant sample from each bottle were collected using syringes to determine the concentrations of NH4 þ-N, NO2 -N, NO3 -N, and COD, which was done for five times during the experiment. At the point where the nitrate concentration between control and experiment group had statistical difference (p < 0.05 by t-test) and the difference of average nitrate concentration of control and experiment group was more than 10% (Kartal et al., 2007a), both autotrophic and mixotrophic anammox consortia were collected from each serum bottle and transferred into RNase-free tubes. After being rapidly frozen in liquid nitrogen, consortia samples were stored at 80 C for subsequent metagenomic, metatranscriptomic, and metabolomic analyses, which could be applied to explore microbial gene functional potential and cross-feedings (Bahram et al., 2018; Lawson et al., 2017). 2.2. Metagenome and metatranscriptome sequencing Autotrophic anammox consortia samples and mixotrophic anammox consortia samples collected in triplicate from each bottle in batch tests were used for total DNA extraction with the FastDNA Spin Kit for Soil (MP Biotechnology, CA, U.S.). DNA concentration and purity was determined using TBS-380 and NanoDrop 2000, respectively. 1% agarose gels electrophoresis system was used to examine DNA quality. Respective triplicate autotrophic and mixotrophic anammox consortia DNA samples were mixed thoroughly for metagenome sequencing (Jia et al., 2018). Total RNA was extracted from triplicate autotrophic and mixotrophic anammox consortia samples from batch tests using the E.Z.N.A® Soil RNA Midi Kit (Omega BioTek, Norcross, GA, U.S.) according to manufacturer's protocols. RNA quality was assessed with a RNA6000 Nano chip (total RNA) in an Agilent 2100 Bioanalyzer and was determined by the RNA integrity number (RIN). Respective triplicate autotrophic and mixotrophic anammox consortia RNA samples were used for metatranscriptome sequencing. To construct paired-end library, DNA was fragmented to an average size of ~300 bp with Covaris M220 (Gene Company 2 Y. Feng et al. / Water Research 165 (2019) 114974

Limited, China). Paired-end library was prepared with TruSeqTM DNA Sample Prep Kit (Illumina, San Diego, CA, USA). Adapters containing the full complement of sequencing primer hybridization sites were ligated to the blunt-end fragments. Paired-end sequencing was accomplished on Illumina HiSeq 4000 platform (Illumina Inc., San Diego, CA, USA) at Majorbio Bio-Pharm Technology Co., Ltd. (Shanghai, China) with HiSeq 3000/4000 PE Cluster Kit and HiSeq 3000/4000 SBS Kits. Total RNAs of consortia samples was subjected to an rRNA removal procedure with the Ribo-zero Magnetic kit according to the manufacturer's instruction (Epicentre, an Illumina® company). cDNA libraries were constructed with TruSeq™ RNA sample prep kit (Illumina). The barcoded libraries were paired end sequenced on the Illumina Hiseq 3000 platform at Majorbio Bio-Pharm Technology Co., Ltd. (Shanghai, China) with HiSeq 4000 PE Cluster Kit and HiSeq 4000 SBS Kits according to the manufacturer's instructions (www.illumina.com). Finally, raw metagenomics and metatranscriptomics datasets were deposited in the NCBI Sequence Read Archive database, and accession numbers were listed in Table S1. 2.3. Metagenomic assembly and binning First, raw metagenomic reads were trimmed by stripping adaptor sequences and ambiguous nucleotides with SeqPrep version 1.1 based on default parameters (Prensner et al., 2011). The trimmed sequences were quality filtered with Sickle version 1.33 based on a minimum quality score of 20 and a minimum sequence length of 50 bp (Wu et al., 2016). Next, contigs and scaffolds were assembled individually for each sample via IDBA-UD with default parameters (Peng et al., 2012). Generated scaffolds were binned into draft genomes based on abundance and the tetranucleotide frequency using the MetaBAT version 0.32.5 with the sensitive model and a minimum contig size of 1500 bp (Kang et al., 2015). CheckM version 1.0.7 was employed to assess the completeness and the quality of recovered draft genomes through single-copy marker genes, which are specific within a phylogenetic lineage (Parks et al., 2015). Draft genomes were deposited in GenBank, and accession numbers were listed in Table S2. 2.4. Phylogenetic analysis of recovered draft genomes Phylosift version 1.0.1 was used to construct the phylogenetic tree (Darling et al., 2014). Marker genes of 37 reference genomes were selected. The concatenated protein alignments of the recovered draft genomes and reference genomes were aligned with MAFFT version 7.310 (Katoh et al., 2002) and concatenated again, using homemade scripts. Finally, a maximum likelihood phylogenetic tree was conducted using the RAxML version 8.2.11 (Stamatakis et al., 2008). Molecular Evolutionary Genetics Analysis (MEGA) software version 7.0 were used to visualize the phylogenetic tree (Kumar et al., 2016). 2.5. Metagenomic and metatranscriptomic analysis First, raw metagenome reads of each sample were mapped to all contigs which were assembled into high-quality draft genomes (completeness 70% and contamination 10%) with the bbmap version 37.75 (minid ¼ 0.95 and ambig ¼ random). Then, the coverage of each bin was calculated by adding the reads which were mapped to the contigs of each bin and normalized by genome size. The open reading frames (ORFs) of each recovered draft genome were annotated through Prodigal version 2.6.3 with the ‘meta’ option for metagenomes, where the minimum nucleotide length was set as 60 (Hyatt et al., 2010). ORFs were then queried against the carbohydrate-active enzymes database (CAZy, accessed on July 2016), the eggNOG database (accessed on December 2016) and Kyoto Encyclopedia of Genes and Genomes pathway database (KEGG, accessed on August 2017) with DIAMOND software, where the blast e value was set as 1e-5 (Buchfink et al., 2014). OrthoANI software was employed to calculate average nucleotide identity (Lee et al., 2016). Genome comparison was visualized via BLAST Ring Image Generator (BRIG) (Alikhan et al., 2011) and Easyfig software (Sullivan et al., 2011). The raw metatranscriptomic reads were trimmed and quality filtered with the SeqPrep and Sickle software respectively. rRNA reads were removed through SortMeRNA version 2.0 (Kopylova et al., 2012) aligning to the SILVA 128 version database (accessed on February 2017). Non-rRNA reads were mapped to all contigs which were assembled into high-quality draft genomes with bbmap software. Read counts of each gene were calculated using htseq-count v0.9.1 with the ‘intersection strict’ parameter (Anders et al., 2015) and normalized as transcripts per million (TPM) values (Wagner et al., 2012). The relative gene expression and pathway expression were calculated by relativizing the median TPM across the draft genome, and by the median TPM value of each reaction in the pathway, respectively. 2.6. LC-MS-based metabolomic profiling and quantitation analysis Metabolic products were extracted from the anammox consortia samples in autotrophic and mixotrophic groups. The protocols used to extract metabolic products were conducted according to those described previously (Guo et al., 2017; Tang et al., 2018b). In brief, anammox consortia suspension was collected by centrifugation, washed using Phosphate buffer saline (PBS), and sonicated using a sonicator. Then, previously cooled methanol was added to the mixture, and the precipitated protein was removed by centrifugation. The supernatants were dried in a nitrogen gas stream. The residues were used to determine metabolite contents via LC-MSbased metabolomic analysis. All tests were conducted in quadruplicate. 3. Results 3.1. Higher nitrogen removal rate with acetate addition In order to investigate the effect of acetate on the nitrogen removal rate of anammox consortia, batch tests were conducted. Based on nitrogen and chemical oxygen demand (COD) consumption (Fig. 1), 155 min was selected as the sludge sample collection point, for NO3 -N concentrations in mixotrophic groups were significantly (p < 0.05 by t-test) and 40% lower than autotrophic groups at this point, which enabled the determination of discrepant expression genes and metabolic pathways of autotrophic and mixotrophic anammox consortia. At 155 min, average NH4 þ-N removal rates were 311.5 ± 15.7 and 284.0 ± 7.2 mg N/(Ld), and the average NO2 -N removal rates were 384.0 ± 41.9 and 359.2 ± 14.2 mg N/(Ld) in autotrophic group and mixotrophic group, respectively. Importantly, the average NO3 -N accumulation rate in the mixotrophic group was 31.0 ± 19.0 mg N/ (Ld), which was significantly lower than that in the autotrophic group (70.9 ± 8.0 mg N/(Ld)) (p < 0.05 by t-test). The DNO3 - -N/ DNH4 þ-N ratio in the mixotrophic group (0.11 ± 0.06) was signifi- cantly lower compared to that of the autotrophic group (0.23 ± 0.03) (p < 0.05). In mixotrophic group, a 59.0 ± 3.7% COD was degraded by anammox consortia with a rate of 245.9 ± 12.4 mg/(Ld) until this point (Fig. 1(d)). Y. Feng et al. / Water Research 165 (2019) 114974 3

3.2. Draft genomes obtained by metagenomic binning Sequencing of DNA extracted from autotrophic and mixotrophic anammox consortia yielded a total of 35,270,644 and 47,527,352 raw reads and 30,542,924 and 44,096,156 clean reads, respectively, after quality control. Next, the clean reads were assembled, and 124,673 and 152,645 contigs were generated with N50 of 2215 bp and 2391 bp. More than 85% of the quality filtered mRNA reads could be mapped to the assembly. After binning, 14 high-quality draft genomes (completeness 70%; contamination 10%) were obtained (Table 1) (Parks et al., 2015), making up 64e69% of the original sequencing reads. A phylogenetic tree of recovered draft genomes is shown (Fig. 2). Bacteria in the anammox community mainly belonged to the phyla Planctomycetes, Chloroflexi, Proteobacteria and Cyanobacteria, while an unclassified bacteria class (CPR1) was also detected. Further details regarding metagenomic and metatranscriptomic read mapping statistics are presented in Supplementary Data. Phylogenetic analysis showed that MAGs AMX1 and AMX2 were closely related to J. caeni and B. sinica, which shared average nucleotide identity values of 99.96% and 99.89%, respectively. Circular maps of two draft genomes (AMX1 and AMX2) aligned to each other are shown in Fig. S1, and the Fig. 1. Concentrations of (a) NH4 þ-N, (b) NO2ˉ-N, (c) NO3ˉ-N in autotrophic and mixotrophic groups in batch tests, and (d) concentration of COD in mixotrophic group. Error bars are defined as s.e.m. (n ¼ 3, biological replicates). * represents p < 0.05 by two-tail t-test. Table 1 Genome statistics of 14 draft metagenome-assembled genomes recovered from the anammox community. Bin ID Taxonomy Completeness (%) Contamination (%) Genome size (bp) Number of scaffolds N50 length GC (%) Predicted genes AMX1 Bacteria; Planctomycetes; Planctomycetia; Planctomycetales; Planctomycetaceae; Jettenia 93 3 3,667,307 88 64,658 40.0 3120 AMX2 Bacteria; Planctomycetes; Planctomycetia; Planctomycetales; Planctomycetaceae; Brocadia 98 2 3,912,452 92 73,161 42.3 3488 PLA1 Bacteria; Planctomycetes 94 7 4,191,635 132 55,477 63.5 3592 PLA2 Bacteria; Planctomycetes 98 2 4,212,226 64 97,350 65.6 3367 CFX1 Bacteria; Chloroflexi; Anaerolineae 98 5 10,306,408 689 28,785 55.1 9340 CFX2 Bacteria; Chloroflexi 89 3 3,596,698 250 24,529 56.9 3451 CFX3 Bacteria; Chloroflexi; Anaerolineae 94 1 3,418,579 52 106,613 60.2 3228 CFX4 Bacteria; Chloroflexi 97 5 4,405,196 173 48,484 56.5 4195 CFX5 Bacteria; Chloroflexi 99 2 5,283,747 22 384,718 51.2 4504 PRO1 Bacteria; Proteobacteria; Betaproteobacteria 95 1 3,183,021 96 79,901 66.5 3346 PRO2 Bacteria; Proteobacteria; Alphaproteobacteria 78 3 2,465,431 337 9252 66.3 2696 PRO3 Bacteria; Proteobacteria; Gammaproteobacteria 84 1 2,290,772 67 42,710 68.3 2187 CYA1 Bacteria; Cyanobacteria 79 0 2,542,560 13 241,538 68.7 2184 CPR1 Bacteria; unclassified bacteria 74 1 796,921 115 9410 44.3 887 4 Y. Feng et al. / Water Research 165 (2019) 114974

nucleotide identity between them was 75.93%. 3.3. Prominent anammox species B. sinica and the gene functional potential with acetate addition The bacterial relative abundance and gene expression levels were calculated using the transcripts per million (TPM) value (Moitinho-Silva et al., 2017). Anammox bacteria were highly enriched in the consortia, comprising approximately 55% of the whole community (Fig. S2). The six abundant species were J. caeni (AMX1 ~51%), Rhodocyclaceae (PRO1 ~20%), Cyanobacteria (CYA1 ~8%), Anaerolineae (CFX3 ~6%), B. sinica (AMX2 ~5%) and Anaerolineae (CFX2 ~2%). Gene expression abundance of J. caeni and B. sinica in the autotrophic group were 78.9 ± 2.7% and 19.3 ± 2.6% respectively, while these were 73.9 ± 4.9% and 23.3 ± 4.5%, respectively, in the mixotrophic groups. Interestingly, PRO1 was the only species whose gene expression abundance was significantly higher in the mixotrophic group compared to that in the autotrophic group (p < 0.05 by t-test). In order to investigate the overall influences of acetate on J. caeni and B. sinica, genes with a significantly changed expression were selected and classified according to cluster orthologous gene (COG) function (Fig. S3). In J. caeni, expressions of a total of 284 genes were significantly down-regulated upon acetate addition, while 36 upregulated genes were expressed. The main decreasing COG function observed in J. caeni was that of [C] energy production and conversion, containing 31 down-regulated genes. On the contrary, acetate greatly facilitated gene expression of B. sinica and a total of Fig. 2. Phylogenetic tree of all recovered draft genomes from the anammox consortia. Draft metagenome-assembled genomes recovered from this study are shown in red, and closely related genomes downloaded from the NCBI are shown in black. GenBank accession numbers of each genome are also presented. Bootstrap support values are represented at branch nodes. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) Y. Feng et al. / Water Research 165 (2019) 114974 5