Unraveling Nitrogen, Sulfur, and Carbon Metabolic Pathways and Microbial Community Transcriptional Responses to Substrate Deprivation and Toxicity Stresses in a Bioreactor Mimicking Anoxic Brackish Coastal Sediment Conditions

Microbial communities are key drivers of carbon, sulfur, and nitrogen cycling in coastal ecosystems, where they are subjected to dynamic shifts in substrate availability and exposure to toxic compounds. However, how these shifts affect microbial interactions and function is poorly understood. Unraveling such microbial community responses is key to understand their environmental distribution and resilience under current and future disturbances. Here, we used metagenomics and metatranscriptomics to investigate microbial community structure and transcriptional responses to prolonged ammonium deprivation, and sulfide and nitric oxide toxicity stresses in a controlled bioreactor system mimicking coastal sediment conditions. Ca. Nitrobium versatile, identified in this study as a sulfide-oxidizing denitrifier, became a rare community member upon ammonium removal. The ANaerobic Methanotroph (ANME) Ca. Methanoperedens nitroreducens showed remarkable resilience to both experimental conditions, dominating transcriptional activity of dissimilatory nitrate reduction to ammonium (DNRA). During the ammonium removal experiment, increased DNRA was unable to sustain anaerobic ammonium oxidation (anammox) activity. After ammonium was reintroduced, a novel anaerobic bacterial methanotroph species that we have named Ca. Methylomirabilis tolerans outcompeted Ca. Methylomirabilis lanthanidiphila, while the anammox Ca. Kuenenia stuttgartiensis outcompeted Ca. Scalindua rubra. At the end of the sulfide and nitric oxide experiment, a gammaproteobacterium affiliated to the family Thiohalobacteraceae was enriched and dominated transcriptional activity of sulfide:quinone oxidoreductases. Our results indicate that some community members could be more resilient to the tested experimental conditions than others, and that some community functions such as methane and sulfur oxidation coupled to denitrification can remain stable despite large shifts in microbial community structure. Further studies on complex bioreactor enrichments are required to elucidate coastal ecosystem responses to future disturbances.


INTRODUCTION
Microorganisms drive and link the biogeochemical cycles of carbon, nitrogen, and sulfur by a variety of redox reactions (Madsen, 2011). Anthropogenic nutrient inputs from land into the ocean constitute a major impact on marine ecosystems, altering seawater and sediment biogeochemistry, and leading to increased primary production that can result in toxic algal blooms and oxygen depletion (Wallenius et al., 2021). Such impacts, combined with ocean warming and consequent seawater stratification and deoxygenation, can further stimulate the production of methane and nitrous oxide, potent greenhouse gases, as well as sulfide and nitric oxide (NO), toxic products of sulfate reduction and denitrification (Van Helmond et al., 2020;Malone and Newton, 2020;Wells et al., 2020;Wallenius et al., 2021). In coastal sediments, ammonium and nitrate can be introduced via agricultural runoff, while sulfide, nitrogen oxides, methane, and ammonium are generated in situ via sulfate reduction, partial denitrification, methanogenesis, and organic matter decomposition, respectively (Egger et al., 2018). Therefore, characterizing microbial communities, interactions, and reactions performed by microorganisms that couple methane, nitrogen, and sulfur cycling is fundamental for understanding biogeochemical cycling and linked greenhouse gas emissions in dynamic coastal ecosystems impacted by anthropogenic activity.
Modeling efforts suggest links between future environmental changes, biogeochemical cycles, and ecosystem functions (Treseder et al., 2012;Nazaries et al., 2013). However, given that most microorganisms are widespread and functionally redundant, they are frequently treated as a "black box" in models-preventing the effective modeling of their reactions and responses. Additionally, recent efforts with laboratory cultures and engineered systems have provided insights about impacts of substrate availability changes on environmentally and economically relevant microbial communities (Chen et al., 2017;Saad et al., 2017;Caffrey et al., 2019;Delgado Vela et al., 2021;Deng et al., 2021;Nie et al., 2021). Yet, few studies have examined microbial community responses to prolonged periods of substrate scarcity or environmental stresses in controlled systems (Bürgmann et al., 2011;Shade et al., 2012). However, these are highly relevant disturbances in coastal ecosystems, where, for instance, nitrogen limitation is a major control on eutrophication (Howarth and Marino, 2006), and sulfide toxicity can lead to mortality of marine life (Grieshaber and Völkel, 1998). Such studies are needed to unravel key microbial interdependencies, competitive interactions, and functional shifts, as well as to comprehend their environmental distribution and resilience under current and future disturbances.
Carbon-, nitrogen-, and sulfur-cycling microbial communities harbor potential for biotechnological applications such as the improvement of wastewater treatment systems. For instance, Deng et al. (2021) proposed that sulfide-driven partial denitrification could be coupled to anaerobic ammonium oxidation (anammox) in future applications, given that rapid oxidation of sulfide to elemental sulfur can prevent toxicity and inhibition of anammox activity. On the other hand, sulfide addition in a controlled aerated bioreactor setting promoted undesirable production of nitrous oxide and nitrite via partial denitrification and dissimilatory nitrite reduction to ammonium (DNRA), respectively (Delgado Vela et al., 2021). These observations indicate that further studies are necessary to understand complex microbial community interactions and to evaluate how they can be best employed in biotechnological applications.
In this research project, we investigated transcriptional stress responses of a complex microbial community enriched in an anoxic bioreactor mimicking dynamic and brackish sediment conditions, where periodic ammonium deprivation, and sulfide and NO toxicity stresses, the chosen stressors in this study, might occur. The culture performed sulfide, ammonium, and methane oxidation at the expense of nitrate via sulfideoxidizing denitrifiers, anammox bacteria, and nitrite/nitratedependent anaerobic methane oxidizers (Arshad et al., 2017). The aims of this study were (i) to understand the effect of periodic ammonium removal on DNRA as a source of ammonium for anammox activity, as well as on general microbial community structure and transcriptional activity, (ii) to enrich sulfide-oxidizing NO reducers (potentially represented by Ca. Nitrobium versatile) while characterizing microbial community structure and transcriptional responses to prolonged sulfide and NO toxicity stresses, and (iii) to unravel potential metabolic reactions utilized by key community members while pursuing the first two aims.

Bioreactor Operation Under Regular and Experimental Conditions
In this study, we used a previously described complex coculture (Arshad et al., 2017). Briefly, this coculture was created by combining biomass from a marine enrichment culture containing sulfide-oxidizing and anammox bacteria (Russ et al., 2014) together with a freshwater nitrate/nitrite-dependent anaerobic oxidation of methane (N-AOM) enrichment culture dominated by Ca. Methanoperedens archaea and Ca. Methylomirabilis bacteria (Ettwig et al., 2016). After approximately 1 year, the stable coculture was characterized for substrate consumption and community composition (Arshad et al., 2017). This same stable coculture was subjected to a 10-week substrate starvation (ammonium removal) experiment (Figure 1, box 1), and later, after regular operation conditions were reestablished, to a 7week sulfide and nitric oxide stress experiment (Figure 1, box 5). Briefly, under regular operation, the reactor culture was kept under anoxic conditions and fed daily with 0.67 mmol of sulfide, 1.4 mmol of ammonium, 3 mmol of nitrate, and 488 mmol of methane. The medium contained a solution of salts, trace elements, minerals, and vitamins as previously described (Arshad et al., 2017), and the pH was kept at 7.1 (Figure 1, box 0 and 4). Under experimental operation, the following conditions were modified. During the ammonium removal experiment, the ammonium concentration in the medium was gradually decreased from 7 to 0 mM over 1 month (Figure 1, box 1). Two more weeks passed without ammonium addition to the reactor in order to ensure that no residual ammonium was present, and the first biomass sample for metatranscriptomic sequencing (T1) was collected and stored at −80 • C. After 10 more weeks without ammonium, another biomass sample for metatranscriptomic sequencing was collected, (T2) and ammonium was fully reintroduced into the medium. For exact dates, see Supplementary Table 1.
In preparation for the next experiment, the reactor was moved into a laminar flow hood. Methane flow was stopped at the time of reactor moving and reintroduced after 15 days (Figure 1, box 3). Once regular operation conditions and substrate consumption were reestablished, the sulfide and NO toxicity stress experiment began. Initially, the NO concentration in the reactor headspace was gradually increased to 1-5% for 10 weeks so that the community could adapt (Figure 1, box 4). Then, NO concentrations were increased to ∼10-13%, and all substrates except sulfide were removed from the reactor, resulting in only sulfide, NO, and medium were provided to the reactor for the following 7 weeks (Figure 1, box 5). A biomass sample for metatranscriptomic sequencing was collected before any NO was added to the reactor (T4) and after the 17 weeks (T5). Originally, the experiment was planned to be conducted in 20 weeks, as during the 10 last weeks the reactor would be maintained under only sulfide and NO. However, the COVID-19 lockdown in the Netherlands resulted in restricted access to the laboratory, and the timing had to be changed from 20 to 17 weeks. Additional metatranscriptomics samples in between experiments were included in this study: T0, approximately 2 months before the ammonium removal experiment, and T3, 6 months after the ammonium removal experiment and 11 months before the sulfide and NO stress experiment (Figure 1 and Supplementary Table 1).
During regular and experimental conditions, the reactor was checked daily for general parameters including pH with an Applisens electrode (Applikon, Delft, Netherlands) and nitrate and nitrite concentrations with MQuantTM colorimetric test strips (Merck, Darmstadt, Germany). During experimental conditions, nitrate and nitrite were also measured with the Griess assay as previously described (Arshad et al., 2017), and ammonium was measured fluorometrically after reaction with 10% orthophthaldialdehyde as previously described (Taylor et al., 1974). NO in headspace samples was measured with a Sievers Nitric Oxide analyzer (NOA280i; GE Power Water & Process Technologies, Boulder, CO, United States), and sulfide in liquid samples was measured with the methylene blue assay (Moest, 1975) (HACH, Loveland, CO, United States). Sulfate was determined via the barium precipitation method using the Sulfate Assay Kit following the instructions of the manufacturer (Sigma-Aldrich, St. Louis, MI, United States). Sample measurements were carried out in technical triplicates.

Nucleic Acid Extractions and Sequencing
Biomass samples for metagenomic sequencing were collected during regular and experimental periods (Figure 1) and stored at −20 • C. All DNA extractions were performed using the FIGURE 1 | Timeline of experiments. Zero indicates when the bioreactor was stably maintained; 1, period of the ammonium removal experiment; 2, after the ammonium removal experiment and restabilization period; 3, reactor moved into a fume hood in preparation for the next experiment; 4, preparatory period for the nitric oxide and sulfide toxicity stresses experiment; 5, nitric oxide and sulfide toxicity stresses experiment. Samples for metagenomics (G0-G5), metatranscriptomics (T0-T5), and RT-qPCR (R0-R1 and R3-R4) were collected as indicated across experimental time points and in between them. The dashed line in the ammonium removal experiment indicates that removal was conducted stepwise (in mM: 7 → 4.5 → 3.5 → 2.5 → 0). Exact dates and experiment durations are provided in Supplementary Table 1. DNAeasy Power Soil Kit (Qiagen, Hilden, Germany) following the instructions of the manufacturer with two modifications: bead beating was performed in a TissueLyser (Qiagen, Hilden, Germany) for 10 min, and autoclaved MilliQ water was used instead of the kit's buffer in the last elution step. RNA extractions were initially performed (for T1 and T2) using the RiboPure Bacteria Kit (Thermo Fisher Scientific, Waltham, MA, United States) following the instructions of the manufacturer. During the stress experiment (T4 and T5), this method no longer resulted in sufficient extracted RNA for sequencing, potentially due to extensive extracellular polymeric substance production, and thus, the method was changed to RNeasy Power Soil Kit (Qiagen, Hilden, Germany), which was also used for two additional samples (T0 and T3). All RNA samples were treated with DNAase I at 37 • C for 30 min from the RiboPure Bacteria Kit (Thermo Fisher Scientific, Waltham, MA, United States). RNA extractions were performed in three biological replicates. Both RNA extraction for RT-qPCR and metatranscriptomics were subjected to the same DNAase treatment. After this treatment, RT-qPCR RNA extractions were tested for DNA contamination. As explained below, RNA samples detected background DNA contamination, equivalent to CT values of RNA extractions of negative controls (water). In addition, we also performed PCR on RNA extractions for RT-qPCR and did not observe any PCR products, thus, validating the effectiveness of the DNAase treatment in trace DNA removal for all RNA extractions employed. DNA and RNA concentrations were measured with a Qubit 2.0 fluorometer using the dsDNA and RNA HS kits (Thermo Fisher Scientific, Waltham, MA, United States). DNA and RNA quality were determined using a NanoDrop Spectrophotometer ND-1000 (Isogen Life Science, Utrecht, Netherlands) and a Bioanalyzer 2100 (Agilent, Santa Clara, CA, United States), respectively.
Metagenomic libraries were prepared using the Nextera XT Kit (Illumina, San Diego, CA, United States) following the instructions of the manufacturer. Enzymatic tagmentation was performed with 1 ng of DNA, followed by the incorporation of indexed adapters and amplification of the library. Amplified DNA libraries were then purified, and their quality and concentration were determined as aforementioned. Libraries were sequenced on an Illumina MiSeq platform (San Diego, CA, United States) using the MiSeq Reagent Kit v3 (San Diego, CA, United States), generating 300-bp paired-end reads. RNA samples of T1 and T2 were rRNA-depleted with MICROBExpress Kit (Thermo Fisher Scientific, Waltham, MA, United States) and MEGAclear kit (Ambion, Life Technologies, Carlsbad, CA, United States). Subsequently, 0.1-4 µg of RNA from T1 to T3 were used to construct strand-specific RNA-Seq libraries. Non-rRNA in RNA-Seq libraries were enriched by selective priming during the first-strand cDNA synthesis reaction, as well as in the final library construction steps using TruSeq Stranded mRNA sample preparation guide (Illumina proprietary catalog RS-122-9004DOC). RNA from these samples was sequenced with an Illumina MiSeq platform (Illumina, San Diego, CA, United States), generating 151-bp singleend reads. Metatranscriptomic libraries of T0, T4, and T5 were constructed using TruSeq stranded mRNA library Kit (Illumina, San Diego, CA, United States). These RNA-Seq libraries were sequenced with an Illumina NovaSeq 6000 platform, generating 151 bp paired-end reads (San Diego, CA, United States). All metatranscriptomes were triplicate RNA extractions and sequencing.
Metatranscriptomic reads were quality trimmed with Sickle v1.33 (Joshi and Fass, 2011) using the sickle se (single end) or pe (paired-end) options for Sanger sequencing (-t sanger). Trimmed transcripts were mapped against the annotated metagenome with Bowtie2 (Langmead and Salzberg, 2012) (bowtie2 -D 10 -R 2 -N 1 -L 22 -i S,0,2.50 -q -a -p 30), allowing only one mismatch. Index stats files were imported into RStudio to calculate transcripts per million (TPM) values according to the formula TPM = [number of reads/(gene length/10 3 )]/10 6 , which were used as unit of gene transcription. Aware of differing extraction and sequencing methods for T1 and T2, these two data points have only been used in separate analyses. TPM values were used for bubble plot generation with the packages ggplot on RStudio. All figures were layout edited using Adobe Illustrator CC v22.1.

Diverse Microbial Community Members Cycled Methane, Nitrate, Nitrite, Ammonium, and Sulfide
For over 1 year since the combination of the two parent enrichment cultures, the bioreactor ran under stable conditions (Figure 1, box 1; Arshad et al., 2017). During the ammonium removal experiment, the bioreactor received the same daily quantities of nitrate (3 mmol), methane (488 mmol), and sulfide (0.67 mmol) as previously described (Arshad et al., 2017). A comprehensive investigation of daily substrate consumption was carried out in the aforementioned study (corresponding to the zero, Figure 1, box 1), which showed that sulfide was completely consumed with a nitrate removal rate of 2.6 mmol per day. The residual nitrate amount fluctuated around 0.4 mmol/day, while no nitrite could be detected in the bioreactor (Arshad et al., 2017). In the present investigation, ammonium was removed from the medium (Figure 1, box 1), and nitrate, nitrite, and ammonium substrate consumption was monitored for a period of 12 days prior to metatranscriptomics biomass (T2) collection (10 weeks without ammonium). Nitrate had a residual concentration of 3 mM, while ammonium, nitrite, and sulfide concentrations remained below the detection limit. The overall nitrate and nitrite concentrations were consistent with the previously determined concentrations (Arshad et al., 2017).
During the sulfide and NO toxicity stress experiment (Figure 1, box 5), methane, nitrate, and ammonium were not provided to the reactor. Only sulfide and NO were added as substrates in order to enrich sulfide-oxidizing nitric oxide reducers. While sulfide concentrations remained below detection limit for the entire experiment, indicating continuous complete sulfide removal, nitrite accumulated to 870 µM during the second week of the bioreactor receiving only sulfide and NO, and remained between 100 and 400 µM during the remainder of the experiment. This indicated that the community removed residual nitrate (834 µM) via denitrification and DNRA for several weeks. Biomass decay was evident from ammonium concentrations (0-500 µM), the visual increase in reactor turbidity, and to a degree, change in granule color from red to black, until the end of the experiment.
The composition of the bioreactor microbial community was investigated via a time series of metagenomic sequencing (Figure 1, G0-G5). In a previous study with this bioreactor, community members included several proteobacteria, Ca. Kuenenia stuttgartiensis, Ca. Scalindua brodae, Ca. Methanoperedens nitroreducens, Ca. Methylomirabilis species, and a novel bacterium within the Nitrospirota phylum, Ca. Nitrobium versatile, potentially linking sulfur and nitrogen cycling (Arshad et al., 2017;Zecchin et al., 2018;Umezawa et al., 2020Umezawa et al., , 2021. Coassembly of eight samples resulted in 30 high (>90% complete, <5% contaminated) and 29 medium (>50% complete, <10% contaminated) quality metagenome-assembled genomes (MAGs) (Figure 2). We will first introduce the most dominant members, and in the next section, we will discuss their change in abundance over time. The only archaeon detected in FIGURE 2 | Up-to-date bacterial core gene (UBCG) tree of 92 concatenated genes extracted from high-and medium-quality metagenome-assembled genomes (MAGs) in this study. Genome completeness (in yellow) and contamination (in purple) values in percent are indicated to the right. MAGs in bold were selected from detailed characterization due to large shifts in abundance across experiments.
Two of the Planctomycetota MAGs were classified by GTDB-Tk as the anammox bacteria Ca. Kuenenia stuttgartiensis (MAG 32) and Ca. Scalindua rubra (MAG 55). The one Nitrospirota MAG 39 was Ca. Nitrobium versatile, which we further characterized in this study. The two MAGs affiliated to Methylomirabilota were classified by GTDB-Tk as Ca. Methylomirabilis oxyfera (MAG 37) and Ca. Methylomirabilis sp002634395 (MAG 38). However, phylogenetic and average amino acid identity (AAI) analyses revealed that the first MAG represented a species of Ca. Methylomirabilis lanthanidiphila [100% AAI to the type genome (Versantvoort et al., 2018)], and the latter was a novel species sharing only 81% AAI to Ca. Methylomirabilis limnetica (Supplementary Figures 1, 2). The genome-to-genome distance between MAG 38 and Ca. M. limnetica ranged from 0.15 to 0.59, with a probability of DNA-DNA hybridization >70% ranging from 0 to 0.08%, supporting AAI results. We have named this novel species Ca. Methylomirabilis tolerans due to its persistence over the sulfide and NO stress experiment (Figure 3). Genome analyses indicated that Ca. Methylomirabilis tolerans had a similar metabolic potential for NO dismutation coupled to methane oxidation via intracellular oxygen generation as previously described, sharing several other characteristics with Ca. Methylomirabilis lanthanidiphila ( Table 1). Finally, several previously identified and novel putative sulfide oxidizers were identified based on key gene analyses (Figures 4, 5)

Metatranscriptomic Analyses Reveal Active Metabolic Pathways in Key Microbial Community Members Cycling Methane, Nitrogen, and Sulfur
To investigate transcriptional responses of the bioreactor microbial community to experimental conditions and to unravel

Reaction
Ca. M. oxyfera (Versantvoort et al., 2018) Ca. M. lanthanidiphila  Ca. M. tolerans (this study) Ca. M. limnetica (Graf et al., 2018) Carbon metabolism microbial pathways likely utilized for metabolic activity in the reactor, we conducted a time series of metatranscriptomic sequencing. With these data, we calculated transcript per million (TPM) values as units of gene transcription. We had particular interest in microorganisms that achieved highest abundances or persisted across experimental time points, and selected several key species for in-depth metabolic reconstruction: two key putative sulfur-oxidizing microorganisms (MAG 39 Nitrobium versatile and MAG 65 Thiohalobacteraceae), two key methane oxidizers (MAG 36 Methanoperedens nitroreducens and MAG 38 Methylomirabilis tolerans), and two key nitrogen-cycling microorganisms-the ammonium oxidizer Ca. Kuenenia stuttgartiensis (MAG 32) and a putative denitrifying Acidobacterium (MAG 60 Thermoanaerobaculia 1). The Nitrospirota MAG 39 Nitrobium versatile, which was the most abundant microbial community member based on genome coverage before the ammonium removal experiment (Figure 3, G0), was estimated to be a 100% complete genome with 5% contamination (Figure 2). Based on metatranscriptomic analyses, we infer that the most likely metabolism performed by this organism was sulfide oxidation coupled to denitrification (Figure 4). A previously unidentified sulfide:quinone oxidoreductase sqr gene in Ca. Nitrobium versatile was one of the mostly highly transcribed functional genes (TPM = 0.2 ± 0.05) of this organism in a thriving period (T0, Supplementary Table 2), along with cytochrome c-oxidizing nitric oxide reductase genes norBC (respectively, TPM = 1.9 ± 0.7 and 2.8 ± 1.2). Sulfide could, thus, be oxidized to elemental sulfur, which might be the substrate for a dissimilatory sulfite reductase (dsrABC) operating in the oxidative direction, despite the two dsrA copies in the genome being phylogenetically related to reductive dsrA genes (Supplementary Figure 3) and the presence and transcription of a dsrD gene. Interestingly, both dsrABC copies in the genome were transcribed, one approximately three times more than the other (Supplementary Table 2, average TPM 0.3-0.4 vs. 0.1). Sulfite could be oxidized to sulfate via transcribed adenosine-phosphosulfate (APS) reductase (aprAB) and sulfate adenylyltransferase (sat) or via a sulfite:cytochrome c oxidoreductase (sorAB) (Figure 4). Although we did not provide reduced sulfur compounds to the reactor other than sulfide, two genes encoding putative sulfide-generating enzymes were identified and transcribed: a sulfhydrogenase polysulfide reductase (hydB subunit beta only, TPM = 0.1 ± 0.06), and a thiosulfate reductase/polysulfide reductase (chain A only, phsA, present twice in the genome, TPMs = 0.09 ± 0.05 and 0.016 ± 0.018). This organism was severely affected by the ammonium removal experiment, as indicated by genome coverage (Figure 3). TPM values suggest that Ca. Nitrobium versatile became a rare community member but remained transcriptionally active across all time points. When the reactor was primed for the sulfide and NO toxicity experiment (T4), transcripts for 161 genes were detected, including norC and phsA but not sqr or norB, while by the end of the experiment (T5), still, transcripts for 130 genes were found (Supplementary Table 2).  Supplementary Table 4. The gene nrfH was absent in MAG 39 Ca. Nitrobium versatile (in gray), but it was previously identified in the first MAG representing this organism (Arshad et al., 2017). In MAG 32 Ca. Kuenenia stuttgartiensis, the internal compartment represents the anammoxosome.
Interestingly, in MAG 39 Nitrobium versatile, several other functional genes were transcribed at low levels: membranebound and periplasmic nitrate reductase genes (narGHI and napA), ammonium-forming cytochrome c 552 -nitrite reductase (nrfA), and anaerobic dimethyl sulfoxide reductase (dmsABC). Although previously identified (Arshad et al., 2017), napB and nirBD were absent in the genome, and nrfH was present in the genome, but no transcription was detected. All genes encoding subunits of electron transport chain proteins were transcribed (Figure 4). Two genes encoding complexes with homology caa 3 -type low-affinity cytochrome c oxidase proteins in Acidobacteria were identified: (i) ctaCFED, which was followed by a downstream sco assembly protein-encoding gene and a cytochrome c 6 (homologous to petJ, K08906), which we have putatively denominated ctaX, for it could be part of the complex, and (ii) sco followed by a downstream ctaDEFC. All these subunits were transcribed, except for ctaF of the latter complex. A Rieske bc 1 complex was encoded by petBCD, which was preceded by a cytochrome c protein-encoding gene immediately upstream. Two other copies of petBC were present in the genome. All these subunits were transcribed except for the first petB (Supplementary Table 4). Finally, the MAG 39 Nitrobium versatile had two ammonium transporter-encoding genes (amtBtype), of which only the second was transcribed, solely in T0 (Supplementary Table 4).
The gammaproteobacterial MAG 65 Thiohalobacteraceae, representing an organism that seemed enriched over the sulfide and NO toxicity experiment (Figure 3), was estimated to be a 98.1% complete genome with 0.3% contamination (Figure 2). Based on metatranscriptomic analyses, we infer that the most likely metabolism performed by this organism was sulfide oxidation coupled to denitrification (Figure 4). A sulfide:quinone oxidoreductase-encoding sqr gene was among the functional genes with the highest transcription by the end of the sulfide and NO toxicity experiment (T5; TPM = 0.43 ± 0.36; Supplementary Tables 2,4). Also highly transcribed were the genes dsrA (TPM = 1.07 ± 0.78), dsrB (TPM = 0.95 ± 0.5), dsrC (TPM = 4.3 ± 1.1), as well as additional dsrABC copies. Sulfite oxidation to sulfate could proceed via proteins encoded by sorA, soeABC (quinone-sulfite dehydrogenase), as well as aprAB and sat, which were all transcribed. Sulfur oxidation system soxBAZYX genes for thiosulfate oxidation and a thiocyanate dehydrogenase tcdh gene were transcribed, indicating that these substrates could also be oxidized, although they were not provided to the reactor. No soxCD gene was identified in the genome. A cytochrome cd 1 -nitrite reductase nirS gene was highly transcribed (TPM = 1.74 ± 0.78), along with all other genes in the denitrification pathway, including a nirK gene encoding a copper-containing nitrite reductase. Ribulose-bisphosphate carboxylase cbbSL genes were among the most highly transcribed in the genome, and all genes encoding subunits of the electron transport chain were transcribed.
MAG 36 (Ca. Methanoperedens nitroreducens) was 99.4% complete with 4.6% contamination, and had all genes in the reverse methanogenesis pathway (Figure 4 and Supplementary Table 4) as well as potential for carbon fixation via the Wood-Ljungdahl pathway with carbon monoxide dehydrogenase/acetyl-CoA synthase (CODH/ACS)-encoding genes and acetate production or assimilation via an acetyl-CoA synthetase-encoding acs gene. All these genes were transcribed by the end of the sulfide and NO toxicity experiment (T5), as well as electron transport chain-encoding genes (Figure 4). The succinate dehydrogenase sdhC subunit was not present in the genome. Remarkably, the genome had genes encoding two copies of the membrane-bound nitrate reductase (narGHI), an ammonium-forming nitrite reductase (nrfAH), a nonelectrogenic cytochrome c-oxidizing (cNOR) nitric oxide reductase (norBC), and an electrogenic quinone-oxidizing (qNOR) nitric oxide reductase (norZ). Methyl-coenzyme M reductase genes mcrBDGA had TPM values around ∼200-400 before the sulfide and NO toxicity experiment (T4), and ∼20-40 after (T5, Supplementary Table 2), indicating that, although Ca. Methanoperedens nitroreducens was still abundant by the end of the experiment (Figure 3), it suffered a degree of inhibition. Three hypothetical genes had highest transcription at T5 and were upregulated relative to T4 (Supplementary Table 2).
A second key methane oxidizer was represented by MAG 38 (Ca. Methylomirabilis tolerans), 94.5% complete with 2% contamination. This genome had all genes for carbon fixation via the Calvin cycle and for the methane oxidation via particulate methane monooxygenase (pmoABC) coupled to nitric oxide dismutation, with four putative nod genes present and transcribed. Additionally, nitrate, nitrite, and nitric oxide reductase-encoding genes were present and transcribed: nxrABC, ammonium-forming cytoplasmic NADH-nitrite reductase nirBD genes, electrogenic cytochrome c-oxidizing (sNOR) nitric oxide reductase-encoding norSY genes, non-electrogenic quinoneoxidizing (gNOR) nitric oxide reductase-encoding norGH genes, and norZ. All four genes encoding the two subunits of the putative nitric oxide dismutases NOD-1 and NOD-2 had some of the highest transcription (TPM ∼7-60) when the community was being primed for the sulfide and NO toxicity experiment, receiving ∼1-5% external NO in the reactor along with methane, ammonium, and sulfide (T4), as well as pmoABC (TPM ∼6-22) and the lanthanide-dependent methanol dehydrogenase-encoding xoxF gene (TPM = 5.5 ± 3.5) (Supplementary Table 2). Although Ca. Methylomirabilis tolerans was still abundant after the sulfide and NO toxicity experiment (Figure 3), TPM values for all aforementioned functional genes decreased 10-to 100-fold from T4 to T5, indicating that the microorganism suffered a degree of inhibition (Supplementary Table 2).
MAG 60 Thermoanaerobaculia 1, 88.3% complete with 1.5% contamination, had the highest normalized genome coverage when the reactor was being primed for the sulfide and NO toxicity experiment (Figure 3, G4). Based on metatranscriptomic analyses, the most likely metabolism performed by this organism was heterotrophic denitrification (Figure 4). A nitrous oxide reductase nosZ gene was among the functional genes with highest transcription (T4, TPM = 0.21 ± 0.17), followed by narGHI and nrfAH. Interestingly, hydrogenase-2 and low-affinity cytochrome c oxidase-encoding ctaCDEF genes were also transcribed, as well as all electron transport chain-encoding genes, including two RNF complexes (Supplementary Tables 2,4).
Finally, Ca. Kuenenia stuttgartiensis, represented by MAG 32, 95.6% complete with 0.6% contamination, was one of two microorganisms that persisted throughout all regular and experimental conditions (Figure 3). Interestingly, the transcription of hydrazine synthase hzsABC genes decreased 18-27× during the ammonium removal experiment (T1 to T2), and 2-5× during the sulfide and NO toxicity experiment (T4 to T5), suggesting that substrate deprivation was a stronger stressor than toxicity. MAG 32 had all genes encoding substrate oxidation and electron transport chain proteins in the anammox pathway (de Almeida et al., 2016), as well as carbon fixation via the Wood-Ljungdahl pathway, which were transcribed at all time points (Figure 4). The second microorganism that persisted throughout all regular and experimental conditions was represented by MAG 44, 93.1% complete with 3.4% contamination, a divergent Planctomycetes genome. Based on metagenomic and transcriptomic analyses, the most likely metabolism performed by this organism was heterotrophic denitrification via narGHI and nirS (Supplementary Table 2). While normalized genome coverage of this MAG varied only between 16× and 50× across regular and experimental conditions, summed TPM values for all genes transcribed in each time point were highest in the preparatory phase for the sulfide and NO toxicity experiment, as well as before and after the experiment (T3-5; TPM = ∼168-361) in comparison with previous time points (T0-2; TPM = ∼33-43; Supplementary Table 2).

Microbial Community Transcriptional Responses to Substrate Removal and Toxicity Stresses Provide Insights Into Community Dynamics and Resiliency
Gene-centric analyses aimed to elucidate the transcriptional responses of the reactor community to two experimental conditions. The first, ammonium removal from the medium, while methane, sulfide, and nitrate were still provided to the reactor, tested the strength of microbial interactions for the supply of ammonium via DNRA to community members-in particular, to anammox bacteria. The second, sulfide and NO addition, while no other substrates were provided to the reactor, aimed to enrich sulfide-oxidizing nitric oxide reducers and to test the limits of microbial community resiliency to these stresses.
A major response to ammonium removal was a general inhibition of methane, ammonium, and sulfide oxidizers, as indicated by downregulation of their key genes mcrA, pmoA, hzsA, sqr, soxB, sorA, and dsrABC (Figure 5). RT-qPCR of mcrA, pmoA, and hzsA confirmed these trends (Supplementary Table 3). Not surprisingly, anammox bacteria were particularly affected, with an approximately 18-fold decrease in hzsA transcription. On the other hand, nitrogen fixation nifDHK genes, ammonium transporter amtB genes, and the ammonium-forming nitrite reductase nrfA genes were upregulated, likely to compensate for the ammonium limitation ( Figure 5). Ca. Methanoperedens nitroreducens accounted for 97.6% of nrfA transcription and 48.2% of amtB transcription 10 weeks after ammonium was completely removed from the medium (T2). Other taxa that shared a large proportion of amtB transcription in T2 were Ca. Methylomirabilis tolerans (14.5%), Ca. Methylomirabilis lanthanidiphila (9.6%), and Ca. Kuenenia stuttgartiensis (3.2%). The hzsA gene of Ca. Kuenenia stuttgartiensis had at least 20× higher transcription than hzsA of Ca. Scalindua rubra across all time points, although Ca. Scalindua rubra had higher genome coverage before experiments were conducted (Figure 3). Finally, minor changes included the upregulation of the periplasmic nitrate reductase napA, nirB, soeA, and quinone-thiosulfate dehydrogenase doxD genes.

DISCUSSION
In this study, we conducted metagenomic and metatranscriptomic analyses in order to unravel metabolic pathways, microbial dynamics, and responses to stresses in a bioreactor community mimicking anoxic and brackish coastal sediment conditions. Microbial interactions and resiliency were tested by removal and addition of key substrates. Removal of ammonium from the medium had major effects on microbial activity and community structure: Ca. Methanoperedens nitroreducens dominated DNRA activity (Figure 5), and after ammonium was reintroduced in the medium, became the most abundant community member (Figure 3). Ca. Kuenenia stuttgartiensis replaced Ca. Scalindua rubra as the most abundant anammox bacterium, several sulfide oxidizers replaced Ca. Nitrobium versatile, and Ca. Methylomirabilis tolerans dominated over Ca. Methylomirabilis lanthanidiphila (Figures 3, 5). These results suggest strong microbial cooperation between Ca. Methanoperedens nitroreducens and Ca. Methylomirabilis species for the exchange of nitrite under the abundance of methane. Before the ammonium removal experiment, TPM values indicate that Ca. Methanoperedens nitroreducens reduced nitrate to nitrite, which was reduced to nitric oxide mostly by Ca. Methylomirabilis lanthanidiphila. Nitric oxide was then shared among three main organisms: Ca. Methylomirabilis lanthanidiphila, Ca. Methylomirabilis tolerans, and Ca. Nitrobium versatile. However, under ammonium limitation, Ca. Methanoperedens nitroreducens further reduced nitrite to ammonium via DNRA, and this had a cascade effect that changed microbial community structure, resulting in Ca. Methylomirabilis tolerans dominating and Ca. Methylomirabilis lanthanidiphila and Ca. Nitrobium versatile concomitantly decreasing in abundance and activity. These results suggest that Ca. Methylomirabilis tolerans was a more competitive organism in scavenging nitric oxide. However, the mechanisms remain to be elucidated. Both Ca. Methylomirabilis species encoded only a lanthanide-dependent methanol dehydrogenase and could be enriched in the bioreactor due to cerium being provided as part of the medium (Versantvoort et al., 2018Guerrero-Cruz et al., 2019).
Anammox bacteria were more strongly limited by ammonium removal than Ca. Methylomirabilis species, presumably due to the dual competition for both ammonium and nitrite, the key substrates for anammox energy metabolism. Although Ca. Methylomirabilis organisms have a lower affinity for nitrite [K s = 7 µM; Ca. M. oxyfera and Ca. M. lanthanidiphila (Guerrero-Cruz et al., 2019)] than anammox bacteria [K s = 0.2-3 µM; Ca. K. stuttgartiensis and Ca. Scalindua sp. (Oshiki et al., 2016)], additional competition for ammonium could have caused anammox to be less fit. Yet, Ca. K. stuttgartiensis [K s < 5 µM for ammonium (Strous et al., 1999)] seemed to better withstand ammonium starvation than Ca. Scalindua rubra [K s = 3 µM for ammonium in several Ca. Scalindua species (Oshiki et al., 2016)], indicating that the affinity of Ca. K. stuttgartiensis for ammonium could be higher than that of Ca. Scalindua. Ca. Nitrobium versatile, on the other hand, seemed to rely on nitric oxide for sulfide oxidation, decreasing in abundance as Ca. Methanoperedens nitroreducens [Ks = 2.1 ± 0.4 mg N L −1 for nitrate (Lu et al., 2019)], and Ca. Methylomirabilis tolerans increased (Figure 3). All three organisms had potential to generate ammonium (Figure 4), but Ca. Methanoperedens nitroreducens dominated DNRA activity ( Figure 5). Therefore, we conclude that DNRA from Ca. M. nitroreducens alone could not sustain anammox activity under ammonium, nitrite, and organic carbon (except methane) limitation, favoring Ca. Methylomirabilis tolerans, which likely outcompeted Ca. N. versatile and anammox bacteria for nitrite and nitric oxide. These results are in contrast to estuary ecosystems in which anammox activity could be sustained by DNRA, positively correlating to sulfide and sediment organic carbon content (Lisa et al., 2014). High sediment organic carbon to nitrate ratio and ferrous iron availability have been reported to favor DNRA over denitrification in estuary ecosystems (Kessler et al., 2018). Having methane as the only organic carbon source (at saturation) and a high nitrate load (3 mmol/day) could have, thus, modulated DNRA activity to become insufficient to sustain anammox and changed microbial community structure.
We also targeted the enrichment of sulfide-oxidizing NO reducers while characterizing microbial community responses to sulfide and NO toxicity. We found genomic potential and transcriptional evidence for this metabolism in Ca. Nitrobium versatile, the dominant microorganism in the bioreactor community before the ammonium removal experiment (Figure 3). However, Ca. Nitrobium versatile became a rare community member after this experiment and could no longer be enriched. Instead, MAG 65 Thiohalobacteraceae, representing an organism likely performing sulfide oxidation coupled to denitrification, increased in abundance (Figures 3-5). The sulfide and NO toxicity experiment revealed unusual resiliency of Ca. M. nitroreducens and Ca. M. tolerans, which persisted as the most abundant community members (Figure 3) even though methane and nitrate were no longer provided to the reactor for the 7 weeks of the experiment, consistent with the downregulation of mcrA and pmoA (Figure 5 and Supplementary Table 3). Given that the hydraulic retention time of the reactor was 5 days, it is unlikely that our sequencing results and the dominance of Ca. M. nitroreducens and Ca. M. tolerans simply reflect older decaying biomass, which would have been removed.
MAG 65 Thiohalobacteraceae, MAG 44 Planctomycetes 1, MAG 60 Thermoanaerobaculia 1, and Ca. Kuenenia stuttgartiensis followed as the most abundant community members at the end of the experiment. While the first three could make use of sulfide or organic carbon from decaying biomass as electron donors and residual nitrate or nitrite and nitric oxide as electron acceptors, Ca. K. stuttgartiensis could have used ammonium generated from decaying biomass and nitrite from residual nitrate reduction or nitric oxide . Of these four MAGs, only MAG 65 Thiohalobacteraceae increased in abundance and had both nirK and norB transcription increased during the sulfide and NO toxicity experiment (Figure 5 and Supplementary Table 2), while the other three suffered a degree of inhibition inferred from both decreases in genome coverage (Figure 3) and marker gene transcriptional activity (Figure 5 and Supplementary Table 2). Given that sulfide remained below detection limit (0.15 µM) during all time points, we infer that this inhibition is attributed to nitrite (∼100-400 µM) and NO (∼10-13% of headspace). While nitrite and NO were likely toxic for several community members, MAG 65 Thiohalobacteraceae seemed to have taken advantage of these compounds as terminal electron acceptors.
MAG 65 Thiohalobacteraceae had similar metabolic potential as the gammaproteobacterium Thiohalobacter thiocyanaticus, which has sulfur oxidation genes encoding FccAB, DsrABC, AprAB, Sat, SoeABC, and SoxXABYZ, as well as thiocyanate dehydrogenase and carbon fixation via the Calvin cycle (Tsallagov et al., 2019). Additionally, MAG 65 had three sqr copies and two sorA copies. Interestingly, MAG 39 N. versatile also had a sulfur oxidation pathway that included sqr, dsrABC, sorA or aprAB, and sat. Umezawa et al. (2020) described the YTD gene cluster composed of genes yedE-like, tusA, dsrE-like, chp-1, and chp-2, which encoded proteins for sulfur disproportionation in Nitrospirota (Umezawa et al., 2020). As in this previous study, we also detected an incomplete YTD gene cluster in Ca. N. versatile (MAG 39), indicating that it might lack sulfur disproportionation potential via this cluster, in contrast to the Nitrospirota microorganism species 45J (Umezawa et al., 2020), which shares 64% average amino acid identity (AAI) to Ca. N. versatile. The lack of a complete YTD gene cluster is also described in the Nitrospirota microorganism Ca. Sulfobium mesophilum (Zecchin et al., 2018;Umezawa et al., 2020), which has sqr and dsrABCD as well as napAB and nrfAH, and shares 57% average amino acid identity to MAG 39 N. versatile. These AAI values indicate that Ca. N. versatile and species 45J are likely part of the same Ca. Nitrobium genus but distinct species, and that together with Ca. Sulfobium mesophilum, they form a family-level taxonomic group (Luo et al., 2014). While conducting our study, Umezawa et al. (2021) isolated a sulfur-disproportionating microorganism that was named Dissulfurispira thermophila, matching the previously described genus Ca. Nitrobium (Arshad et al., 2017).
Ca. Sulfobium mesophilum and Ca. N. versatile both have dsrA genes that affiliate with bacterial-type reductive sequences, and a dsrD gene, suggested as a potential marker for sulfate reduction, absent in sulfur oxidizers that utilize the reverse Dsr pathway (Rabus et al., 2015;Anantharaman et al., 2018;Dalcin Martins et al., 2018). However, our metatranscriptomic results suggest that Ca. N. versatile was performing sulfide oxidation. Therefore, it seems that Ca. N. versatile, similar to the deltaproteobacterium Desulfurivibrio alkaliphilus, was yet another example of a microorganism disguised as a sulfate reducer (Thorup et al., 2017). Interestingly, Ca. Sulfobium mesophilum and Desulfurivibrio alkaliphilus could couple sulfide oxidation to DNRA, a metabolic potential also present in MAG 39, with transcripts of nrfA detected in our study (T0; TPM = 0.030 ± 0.031). Given that norBC had significantly higher transcriptional activity than nrfA and other denitrification genes (T0, TPM norB = 1.9 ± 0.7, TPM norC = 2.8 ± 1.1), it is more likely that Ca. N. versatile coupled sulfide oxidation to nitric oxide reduction. Given that Ca. N. versatile could not be enriched under sulfide and NO, we hypothesize that MAG 65 Thiohalobacteraceae represented a more competitive microorganism. However, the mechanism and substrate affinities remain to be elucidated.
In our study, methanotrophs could withstand prolonged periods of ammonium, methane, and nitrate deprivation as well as exposure to elevated concentrations of nitrite and nitric oxide. During these disturbances, Ca. Methanoperedens nitroreducens and Ca. Methylomirabilis tolerans did not thrive, as indicated by downregulation of mcrA and pmoA, but tolerated stresses and persisted as abundant community members. These results suggest that methane oxidation could be a relatively stable community function in coastal ecosystems under the stresses investigated in this study and as long as the system stays anoxic. Future investigations are needed to elucidate such dynamics. Interestingly, the Ca. Methanoperedens nitroreducens genome (MAG 36) encoded nitric oxide reductases (Figure 4), which had not been described before in this organism. The qNORencoding gene, present in the same contig as mcr, nar, and nrf genes, had increased transcriptional activity by the end of the sulfide and NO stress experiment (Figure 5). Future studies should further evaluate lateral transfer of nor genes in Methanoperedenaceae, which seem prone to acquire novel metabolic traits via later gene transfer events (Leu et al., 2020), potentially via recently described borgs (Al-Shayeb et al., 2021), as well as metabolic flexibility of Ca. Methanoperedens nitroreducens under methane deprivation. Hypothetical genes in Ca. Methanoperedens nitroreducens identified in this study with high transcription and upregulation in response to stresses could be targets for these future investigations.
Sulfide oxidizers had differential resilience to the stresses investigated in this study, dynamically changing in abundance and transcriptional activity across regular and experimental conditions. However, sulfide was completely removed at all time points, likely due to functional redundancy in the microbial community, indicating that sulfide oxidation could also be a relatively stable community function under the investigated stresses. Finally, in our study, denitrification was a dominant nitrogen-cycling pathway, as previously suggested in coastal sediments of the Bothnian Sea (Rasigraf et al., 2019)particularly, the nitric oxide reduction step, as indicated by norB TPM values (Supplementary Table 2). Anammox activity, as indicated by hzsA TPM values, also had a significant contribution to nitrogen cycling, and dominated at the end of the sulfide and nitric oxide stress experiment (T5). This community function was highly impacted by ammonium deprivation, but was restored when favorable conditions were reestablished. On the other hand, DNRA, as indicated by nrfA TPM values, was a minor nitrogen-cycling pathway, but under ammonium deprivation, it became more significant. Future studies should investigate whether anammox and DNRA similarly oscillate in coastal ecosystems under the stresses investigated in this study.

CONCLUSION
This study contributes to the elucidation of metabolic pathways for carbon, sulfur, and nitrogen cycling in a bioreactor community mimicking anoxic and brackish coastal sediment conditions, as well as shifts in microbial abundance and transcriptional activity in response to prolonged substrate deprivation and exposure to toxic compounds. Together with follow-up studies, these results will help in understanding complex microbial interactions and functions in dynamic coastal ecosystems, and should be considered into future modeling efforts that aim to predict coastal ecosystem responses to environmental change.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

AUTHOR CONTRIBUTIONS
HJMO, MSMJ, CUW, and PDM planned and designed the study. PDM, MJEM, JMK, AA, and HTO executed the research project. PDM and MJEM wrote the manuscript with input from all co-authors. All authors read and approved the final version submitted.