Modulation of carbon-to-nitrogen ratio shapes the microbial ecology in a methanol-fed recirculating marine denitrifying reactor

Acclimation of denitrifying biomass in recirculating reactor

The medium was the commercial seawater Instant Ocean (IO) medium (Aquarium systems; Mentor, OH, USA), the same used in the original denitrification process. The IO medium (30 g/L) was autoclaved, and one mL/L of sterilized trace element solution (FeSO₄•7H₂O 0.9 g/L; CuSO₄•5H₂O 0.03 g/L; MnCl₂ •4H₂O 0.234 g/L; Na₂MoO₄•2H₂O 0.363 g/L) was added.

The denitrifying biomass was provided in the form of biofilm, immobilized on polyethylene supports of the Bioflow 9-type (1.2 mm diameter) that had developed in the denitrification reactor at the Montreal Biodome, but was preserved in 20% (v/v) glycerol (in IO medium) at −20 °C (Payette et al., 2019). Frozen supports were thawed and submerged for 5 days in the IO medium at 4 °C, supplemented with 0.3% (v/v) methanol. The biofilm was also scraped from the supports and inoculated in the same medium. The biomass free or attached to the supports was transferred to the reactor with seventy empty Bioflow 9-type supports. These empty supports were previously washed in a 35% hydrochloric acid bath overnight and then rinsed a minimum of three times with demineralized water.

The reactor consisted of a glass vessel with an airtight 550-mL working volume, installed in a glove box continuously purged with nitrogen gas and equipped with an oximeter to confirm the absence of O₂ (Cucaita, Piochon & Villemur, 2021). A 500-mL volume of culture medium supplemented with 0.3% methanol and 21.4 mM NaNO₃ (300 mg-N/L) (C/N of 3) was added to the reactor. A 4-L tank filled with the same medium was appended to the circuit during the acclimation period. The medium was recirculated with a peristaltic pump at 20 mL/min, and the reactor was run at room temperature (ca. 22 °C). The excessive volume of gas produced due to the denitrifying activities was evacuated in a graduated bottle inverted in water and connected to the reactor. Gas production was recorded by water displacement. All equipment, including culture medium, was sterilized before the start of the experiment; however, the operation of the reactor was not carried out under sterile conditions, although it was operated in closed environment of the glove box. The culture medium was replaced regularly as soon as NO₃⁻ was exhausted. When biofilm was apparent on the empty supports, supports with the original biofilm were removed. When denitrifying activities (removal of NO₃⁻ and NO₂⁻) was stable, the reactor was emptied, and half supports were set aside and stored at 4 °C in IO medium for future used. From this, the C/N assays were carried out with appropriate NO₃⁻ and methanol concentrations.

C/N assays

The reactor was run sequentially under eight different C/N ranging from 1.5 to 7.5 according to Table 1, by adjusting the concentration of either NO₃⁻ or methanol. The maximum addition of Na⁺ concentration with NaNO₃ in the reaction was 42.8 mM, which added less than 10% Na⁺ equivalent in the reactor (30 g/L salt, 513 mM NaCl equivalent). The reactor was operated with medium recirculated at 20 mL/min (not with the 4L tank) in sequential batch mode, according to a cycle of three phases: a filling phase, a reaction phase that lasted until complete NO₃⁻ reduction, and an emptying phase. The eight conditions were applied twice, which we refer as Period 1 (P1) and Period 2 (P2), with the same biofilm supports (no new inoculation). The reactor ran from November 13, 2019, to June 14, 2020 (31 weeks). During Condition 6 (Period 2), the running operation was abruptly interrupted because of the COVID-19 crisis, which lasted 68 days, during which the supports were transferred at 4 °C in IO medium containing NO₃⁻ and methanol. During P1, the operating time for each condition ranged from 96 to 288 h with 1 to 3 medium changes. During P2, the same conditions were reapplied over a shorter interval for Conditions 1 to 5 (48 to 72h) as each condition was running with no medium changes. For Conditions 6 to 8, the reactor was running during much longer time (120 to 240 h) because of the decrease in the denitrification performance probably due to the storing period (COVID crisis).

Liquid samples (one mL) were taken at sampling ports and centrifuged; the supernatant was used to determine the concentrations NO₃⁻, NO₂⁻, ammonium and methanol. The pellets of some of these samples were kept for DNA or protein (or both) extractions. A 10-mL volume of gas was also taken from the head space to measure N₂O and CO₂.

At the end of each condition at the time of the emptying phase, three supports were taken to determine the protein content of the biofilm and to extract the nucleic acids (DNA and RNA). These supports were replaced by three others from the supports that were set aside previously, in order to keep a fixed number of supports. At the emptying phase, the medium was also recovered to collect the suspended biomass and determine its protein content. The biomass collected to extract the RNAs was scraped from supports under the anoxic environment of the glove box and immersed in the extraction buffer (50 mM Tris-HCl, 100 mM EDTA, 150 mM NaCl, pH 8.0) and water-saturated phenol pH 4.3 (v/v). The sample tubes were then immediately frozen in liquid nitrogen and stored at −70 °C until use. For DNA extraction and total protein assays, the biomass from the supports was scraped and then frozen at −20 °C in five mL of sterile culture medium until use. To recover the biomass in suspension, 300 mL of the culture medium was centrifuged at 16,000× g and the pellet was mixed with five mL of sterile medium and stored at −20 °C until use.

Analytic methods

Measurements of the concentrations of NO₃⁻, NO₂⁻ and protein content, and the determination of the specific NO₃⁻ reduction rates and the specific NO_x (NO₃⁻ + NO₂⁻) rates were carried out according to Cucaita, Piochon & Villemur (2021). N₂O and CO₂ concentrations (in ppmv) in the headspace were determined by gas chromatography. Headspace samples (10 mL) were collected using a Pressure Lok gastight glass syringe (VICI Precision Sampling Inc., Baton Rouge, LA, USA) and were injected through the injection port of a gas chromatograph equipped with a thermal conductivity detector and electron-capture detector (7890B seriesGCCustom, SP1 option 7890-0504/0537; Agilent Technologies). The amount of N₂O (in µmole) in the reactor was calculated as described in Mauffrey et al. (2017). The NO₂⁻ and N₂O showed transitory dynamics with appearance that peaked after a certain time, then reduced completely (Figs. 1G and 1H). To quantify these dynamics, we calculated the area under the curve of the amount of NO₂⁻ or N₂O by the time (mM-h) they appeared and then disappeared (Data S1). The CO₂ production rates were calculated from the linear portion of CO₂ produced by the time, which occurred mainly during the first 5–15 h (Data S1). This period corresponded when NO₃⁻ reduction was occurring.

Methanol was measured by gas chromatography coupled with a flame ionization detector (GC-FID). Nitrogen gas was used as a carrier gas. The culture medium samples were supplemented with 5% isopropanol as an internal standard. The standard range included solutions of 0.01 to 0.2% methanol. The total amount of methanol consumed was determined at the end of the reaction phase. The methanol consumption rates followed a first order dynamics, where a k factor was calculated as [ln(C_t/C₀)]/time with C_t methanol concentration at the sampling time and C₀ the initial methanol concentration (Data S1). Ammonium concentration was measured by colorimetric method as described by Mulvaney (1996).

Polymerase chain reaction-denaturing gradient gel electrophoresis and 16S rRNA amplicon sequencing

Supports were taken at the end of each condition and period for a total of 16 samples. Biofilm from the original biofilm was also taken for DNA extraction. DNA extractions of biomass on supports were performed as described (Geoffroy et al., 2018). The PCR amplifications (one PCR for each sample) of the V3 region of the 16S rRNA genes for denaturing gradient gel electrophoresis (DGGE) experiments were performed as described (Lafortune et al., 2009) with the 341f and 534r primers.

The bacterial composition of the biofilm was determined by amplicon sequencing of the 16S rRNA V3–V4 regions with the S-D-Bact-0341-b-S-17 (5′TCGTCGGCAGCGTCAGAT GTGTATAAGAGACAG-CCTACGGGNGGCWGCAG 3′) and S-D-Bact-0785-a-A-21 (GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG-GACTACHVGGGTATCTAAT CC) primers (Klindworth et al., 2013) linked with Illumina sequences (underlined sequences). The amplifications were carried out in 25-µL volume containing the commercial Taq polymerase buffer, 0.25 mM MgSO₄, 0.2 µg BSA, 200 nM of each primer, 0.5-unit AccuPrime™ Taq DNA Polymerase (Invitrogen, Carlsbad, CA, USA) and 50 ng DNA extract. The reactions were run at 94 °C 5 min, 30 cycles at 94 °C 30 s, 55 °C 30 s, and 68 °C 30 s, and finally at 68 °C 10 min. Residual primers and primer dimers were removed with the AMPure XP magnetic beads kit according to the manufacturer (Beckman Coulter, USA). The amplicons were then amplified a second time in which indexes were added. The amplifications were carried out in a 25-µL volume containing the Taq polymerase buffer, 0.25 mM MgSO₄, 0.2 µg BSA, 400 nM index primers, 0.5-unit AccuPrime™ Taq DNA Polymerase (Invitrogen, Carlsbad, CA, USA) and 5 µL amplicon. The reactions were run at 94 °C 3 min, 8 cycles at 94 °C 30 s, 55 °C 30 s, and 68 °C 30 s, and finally at 68 °C 5 min. Residual primers and primer dimers were once again removed with AMPure XP magnetic beads kit, and the purified amplicons were quantified with the Quant-iT PicoGreen dsDNA Assay kit (Invitrogen). The amplicons were sent for sequencing by Illumina MiSeq PE250 (250 bp paired-end sequencing reactions) at the Centre d’expertise et de services Génome Québec (Montréal, QC, Canada). Sequence quality per base was >32. Sequencing reads were processed by the dada2 pipeline at the Galaxy server (https://usegalaxy.eu/) to generate ASVs (amplicon sequence variants). The parameters were: (1) FilterAndTrim of fastqR1 and R2 with Truncate reads at quality threshold set to 2, Trim start of each read set to 15, and Remove short reads set to 20; (2) learn errors with Magnitude of number of bases to use for learning set to 8, Error function set to loess to estimate error rates from transition counts, and Maximum number of times to step through the selfconsistency loop set to 10; (3) dada: Remove sequencing errors, (4) mergePairs with Minimum length of the overlap set to 12; (5) removeBimeraDenovo, with homology set to 80%. The taxonomic assignment was carried out at the Silva web site (https://www.arb-silva.de/) with sequence identity cut off set to 80%. Number of reads ranged from 27 268 to 33 031 for the recirculating reactor samples, and 423 208 reads from the original biofilm sample. Rarefaction analyses were performed with the Analytic Rarefaction v1.3 software downloaded from the UGA Stratigraphy Lab web site (http://stratigrafia.org/) with increments of 10. All samples showed rarefaction saturation. ASVs affiliated to the same lineage were grouped together, in order to establish the proportions of each taxon in the samples. Raw sequencing data were deposited in Sequence Read Archive (SRA) at the National Center for Biotechnology Information (NCBI: https://www.ncbi.nlm.nih.gov/) (Data S2).

Metatranscriptomes

The RNA extraction was previously described by Mauffrey, Martineau & Villemur (2015). RNA quality was verified by agarose gel electrophoresis and by spectrophotometry (Nanodrop) with ratio OD ₂₆₀/OD ₂₈₀ > 1.8. The RNA samples were sent for sequencing using the Illumina Method (NovaSeq 6000 S4 PE100). Library preparation and sequencing were performed by the Centre d’expertise et de services Génome Québec (Montréal, QC, Canada). Ribosomal RNA were depleted using the Ribo-ZeroTM rRNA Removal Kits (Meta-Bacteria; Epicentre, Madison, WI, USA). RNAseq reads were deposited in SRA (Data S2). Raw sequencing reads were trimmed using fastq_quality_filter (1.0.2) to remove low quality reads (score < 20) (Gordon & Hannon, 2010) , and the paired reads were then merged.

Statistical analysis

Pearson correlation analyses were performed by GraphPad Prism v10.4.2. Relative gene expression profiles of strain GP59 were analyzed by Principal coordinate analysis (PCoA) (Canoco version 5.15; Ter Braak & Šmilauer, 2018), with Log (1 * X + 1) calculating matrix of distances, using percentage difference (Bray–Curtis distance). PERMANOVA (Vegan package version 2.6.2, command Adonis2) were performed for significance with 999 permutations (Anderson, 2017).

Acclimation of the denitrifying biomass

The Biodome denitrification system was dismantle in 2006, from which we took all the biofilm supports and preserved them in glycerol at −20 °C. Previous works by Laurin et al. (2006) demonstrated no differences in the recovery of the denitrifying activities of the biofilm when it was stored at −20 °C or −70 °C for 17 months. Thirteen year later, the frozen supports were used to inoculate the recirculating reactor for the colonization of new supports. This strategy has been successfully used in our previous studies (Payette et al., 2019; Villemur et al., 2019). We refer to this biofilm as the Original Biofilm to distinguish it from the one that developed in the recirculating reactor. The reactor was operated under sequential-batch conditions with the IO medium, the same as that used by the Biodome, and supplemented with 0.3% methanol and 21.4 mM NO₃⁻ (0.891 g-C/L and 0.2996 g-N/L: C/N of 3). These concentrations were based on a previous report by Payette et al. (2019) that succeeded culturing the biofilm under these conditions. After 37 days, stable denitrifying activities (reduction of NO₃⁻ and NO₂⁻) occurred in the reactor, and fresh biofilm was visually observed on the supports (Fig. S1). This biofilm is referred to as the RR biofilm (RR for recirculating reactor). The C/N assays were then carried out. Day 0 (D0) was set after the 37-day acclimation period.

C/N assays

The reactor with supports fully colonized by denitrifying biomass was operated sequentially under eight different C/N (1.5 to 7.5) (Table 1). These conditions were performed twice (P1 and P2) with the same biomass to assess the reproducibility of the results. The reactor was operated for 31 weeks during these two periods, where several physico-chemical parameters were measured (Figs. 1A to 1H). The pH was stable (7.6–8.0) and the protein content on the supports (35), representing the total microbial biomass, ranged from 10 to 25 mg. No ammonium was detected under any conditions.

NO₃⁻ reduction was occurring with transient NO₂⁻ and N₂O accumulation (Figs. 1C and 1D and 1H; Data S1). On average, the NO_x (NO₃⁻ + NO₂⁻) reduction rates (Fig. 1A) were 77% of the NO₃⁻ reduction rates (Fig. 1B), due to the transient NO₂⁻ accumulation (Fig. 1C). Although we operated the reactor twice under all eight C/N ratios, we encountered a significant obstacle, the COVID-19 crisis, where the reactor had to be abruptly shut down during P2-Condition 6, and stored in cold room for several weeks, which negatively impacted the denitrification performance (Fig. 1). If we assume the NO_x reduction rates represent the denitrification rates (NO and N₂O reduction have minimal influence on these rates), the performance of the reactor ranged from 0.166 to 0.328 mM h⁻¹ if we exclude post-COVID data. These rates averaged 0.119 mM h⁻¹ in P2-Conditions 6 to 8. We estimated the residuals of NO_x, NO₃⁻ and NO₂⁻ after four operating days (Figs. 1A–1C). We observed these residuals in P1 only under Conditions 7 and 8, with 28.5% residual NO_x in Condition 8. In P2, Conditions 6 to 8 showed higher levels of residuals due to the lower reactor performance with 50.6% residual NO_x in Condition 8.

The peak of N₂O concentrations appeared mainly in the first 10-20 h of operating conditions and ranged from 0.04 mM to 0.4 mM (Fig. 1G and 1H). No residual N₂O was found at the end of each condition. The amount of methanol consumed ranged from 0.031% to 0.147% (v/v), which corresponded to between 9% to 53% of the initial methanol concentrations (i-Methanol) (Fig. 1E). The methanol consumption rates (Fig. 1E) followed first-order dynamics, while the CO₂ production rates (Fig. 1F) were at their maximum during the first 5-15 h (Data S1).

Correlation analyses were first performed between (i) the i-Methanol, (ii) the initial NO₃⁻ concentrations (i-NO₃⁻), (iii) the C/N applied to the reactor, and the physico-chemical parameters measured during the operating conditions (Table 2). These analyses revealed that changes in C/N did not correlate with the denitrification dynamics (NO₃⁻ reduction rates, the NO_x reduction rates, the NO₂⁻ and N₂O dynamics; p > 0.05). However, these C/N changes did correlate with the methanol consumption rates, and the CO₂ production rates (0.05 < p < 0.001; Table 2). Correlations (p < 0.05) occurred between the i-NO₃⁻ and (i) the NO₃⁻ reduction rates, (ii) the NO_x reduction rates and (iii) the NO₂⁻dynamics (Table 2), meaning higher i-NO₃⁻ resulted in higher denitrification dynamics. The i-Methanol correlated with the methanol consumption rates and the CO₂ production rates (0.05 < p <0.001; Table 2), meaning higher i-Methanol resulted in higher carbon dynamics. Correlation analyses were also performed between the respective measured parameters (Table 2). As expected, strong correlation (p < 0.0001) was observed between the NO₃⁻ reduction rates and the NO_x reduction rates (Table 2). The N₂O dynamics did not correlate with any conditions.

Bacterial community

The evolution of the bacterial community of the RR biofilm during the different operating conditions and the Original Biofilm was assessed by performing PCR amplification of the 16S rRNA genes and by analysing the composition of the amplicons by DGGE and amplicon sequencing. We first compared the DGGE migration profiles of the Original Biofilm between a fresh sample collected in 2004, when the Biodome system was operational, and the frozen stock sample used to inoculate the reactor. The results showed that the two main DNA bands of M. nitratireducenticrescens and H. nitrativorans were present, although the M. nitratireducenticrescens band is more intense (higher proportion) in the 2004 sample. These results suggest that no major changes in the bacterial community occurred in the frozen stock.

The polymerase chain reaction-denaturing gradient gel electrophoresis (PCR-DGGE) migration profiles showed striking differences between the Original Biofilm sample and the samples collected from the reactor (Fig. 2A). The same migration profiles were obtained with samples collected during the second period, suggesting that the interruption of the reactor during the COVID crisis did not deeply affect the composition of the bacterial community. The migration profile of the D0 sample revealed some similarities with those derived from samples collected during the operating conditions. The PCR-DGGE migration profiles of the later samples did not reveal substantial changes between them, suggesting that the bacterial community remained stable during the operating conditions.

We performed 16S amplicon sequencing on total DNA extracted from the RR biofilm samples. The first sample was collected at D0 (C/N of 3.0; 21.4 mM NO₃⁻, 0.3% methanol) to provide the bacterial diversity profile at day 0. Two samples came from the same conditions (C-5, 21.4 mM NO₃⁻ and 0.45% methanol: C/N of 4.5) but one in P1 (C5-1; Day 30) and the second in P2 (C5-2; Day 120). Finally, the sample from condition 6 (21.4 mM NO₃⁻ and 0.75% methanol; C/N of 7.5) during P2 (C6-2; Day 200) was chosen because it came from the reactor operated after the COVID break. For comparison, we added data that we obtained from the Original Biofilm of the Biodome denitrification system (estimated C/N of 0.86) (Labbé, Parent & Villemur, 2003).

We identified forty-two taxa among the five biofilm samples (Data S3) and determined their proportions in the biofilm samples. Figure 2B illustrates the proportions of the most abundant taxa observed in the samples. Principal coordinate analysis (PCoA; Data S3) was performed on the bacterial population profiles of these samples. Due to the uniqueness of the samples (no replicate), we could not perform statistical analysis such as PERMANOVA to determine the significance of variations between the five profiles. However, we believe it is reasonable to conclude that the bacterial population profiles of the C5-1, C5-2 and C6-2 biofilm samples are similar, and those of the D0 and Original Biofilm samples differ from each other, and from the other ones (Fig. 2C). These results concur with the PCR-DGGE results. Therefore, for the analysis of the composition of the bacterial community in the biofilm samples, we considered the C5-1, C5-2 and C6-2 bacterial profiles as a whole, distinct from those of the Original Biofilm sample and the D0 biofilm sample.

Thus, higher proportions of Hyphomicrobiales, Roseovarius spp. and Oceanibaculaceae, and a much lower proportion of Methylophaga spp. were found in the Original Biofilm sample compared to the four RR biofilm samples. In the D0 biofilm sample, the highest relative abundance of Methylophaga spp. (69%) was found, which is consistent with the DGGE migration profile of the D0 biofilm sample, where an intense DNA fragment associated with M. nitratireducenticrescens was observed. For the C5-1, C5-2 and C6-2 biofilm samples, increases in the relative abundance of Marinicella spp., Stappiaceae spp., Muricauda (Flagellimonas) spp. and Oceanibaculaceae were observed compared to the D0 biofilm sample with an approximately 50% decrease in relative abundance of Methylophaga spp. These three biofilm samples had higher alpha diversity indexes compared to the Original Biofilm (Table 3; Data S3). Due to the high proportion of Methylophaga spp. in the D0 biofilm sample, the alpha diversity indexes were lower than those of the four other biofilm samples. Regarding beta diversity (Table 3), examination of common taxa between samples revealed that taxa found in the RR biofilm samples shared the highest number of taxa. All these results showed that the original bacterial community evolved through diversification in the recirculating reactor from Day 0, and then remained relatively stable during the operating conditions.

Metatranscriptome analysis

We derived metatranscriptomes of the microbial community from four RR biofilm samples to determine (1) the relative transcript profiles of M. nitratireducenticrescens strain GP59, the main denitrifier in the RR biofilm; and (2) the potential metabolic pathways active in the RR biofilm.

The metatranscriptomes were derived during P1 from: Condition 1 (C1-1; C/N of 1.5; 10.7 mM NO₃⁻, 0.075% methanol; Day 4), Condition 4 (C4-1; C/N of 3.0; 21.4 mM NO₃⁻, 0.3% methanol; Day 24). Condition 6 (C6-1; C/N of 7.5; 21.4 mM NO₃⁻, 0.75% methanol; Day 36), and Condition 8 (C8-1; C/N of 3; 42.8 mM NO₃⁻, 0.6% methanol; Day 89). Condition 1 and Condition 6 were chosen as the lowest (1.5) and highest (7.5) C/N applied to the reactor, respectively. Although Condition 4 and Condition 8 had the same C/N (3.0), they were adjusted with distinct levels of NO₃⁻ and methanol. In addition, we derived the metatranscriptome of the Original Biofilm for comparison.

Relative transcript profiles of M. nitratireducenticrescens strain GP59

The proportions of transcript reads associated with strain GP59 in the metatranscriptomes were similar across the four RR biofilm RNA samples (average 7.15%), four times higher than in the Original Biofilm RNA sample (Data S4). The relative gene expression profiles of strain GP59 were then derived from the five RNA samples. Due to the uniqueness of the sample (no replicate), we were unable to perform statistical analysis such as PERMANOVA to determine the significance in the variations between the five expression profiles. However, we took this opportunity to compare these profiles with those derived from our previous studies. These include:

1-The Original Biofilm cultured in vials under batch-static conditions in artificial seawater (ASW, 2.75% NaCl), with two different NO₃⁻ concentrations and temperatures: < <300N23C> > (21.4 mM NO₃⁻, 0.15% methanol, C/N of 1.5, 23 °C); and < <900N30C> > (64.2 mM NO₃⁻, 0.45% methanol, C/N of 1.5, 30 °C) (Payette et al., 2019; Villemur et al., 2019).

2- The Original Biofilm cultured in vials under batch-static conditions in low salt concentration ASW: < <0%NaCl> > (ASW with 0% NaCl, 21.4 mM NO₃⁻, 0.15% methanol, C/N of 1.5, 23 °C) (Payette et al., 2019; Villemur et al., 2019).

3-Planktonic cultures (triplicate) of strain GP59 cultured in a Methylophaga-specific culture medium (Methylophaga-1403, 0.3% methanol; C/N of 3) under oxic conditions with NO₃⁻ (21.4 mM, PlkON) or without NO₃⁻ (PlkO), and under anoxic conditions with 21.4 mM NO₃⁻ (PlkAN) (Lestin & Villemur, 2024).

All these cultures differed in culture media (IO, ASW, low salt ASW, Methylophaga-1403), C/N, operating mode (batch-static, continuous, recirculating sequential-batch), and physiology (biofilm, planktonic). The objective of this analysis was to reveal metabolism pathways of strain GP59 affected by the different culture conditions. PCoA clearly showed distinct clustering of each category of cultures, explained by 68.5% variations (Fig. 3; Data S5). PERMANOVA were performed with the planktonic cultures (triplicate: Methylophaga-1403 medium, anoxic or oxic, presence or absence of NO₃⁻, planktonic), and the four RR biofilm samples assuming that their relative transcript profiles of the RR biofilm RNA samples were similar enough to be considered as replicates (IO medium, anoxic, biofilm, presence of NO₃⁻). The results showed, as perceived in Fig. 3, that the medium/physiology (IO (biofilm) vs M1403 (planktonic)), the operating conditions (oxic vs anoxic), and the methanol and NO₃⁻ concentrations had substantial effects (p < 0.01) on the relative transcript profiles of strain GP59 (Data S5).

We examined the relative transcript levels of genes involved in the N cycle. Thirty-three genes or gene clusters were selected for their involvement in denitrification, N assimilation, regulation and NO-responses (Data S5). PCoA was performed to derive the preferential expression profiles in the different culture types. Figure 4 illustrates the PCoA that explained 78% of the variations. These genes and cultures fell into three specific quadrants: (i) the lower left quadrant with the planktonic oxic cultures (PlkO, PlkON); (ii) the lower right quadrant with the 0%NaCl biofilm (non-marine conditions) and the Original Biofilm; and (iii) upper quadrants with the planktonic anoxic cultures (PlkAN) and the RR biofilm samples. PERMANOVA showed this time that only the media/physiology (IO (biofilm) vs M1403 (planktonic)), and the operating conditions (oxic vs anoxic) had substantial effects (p < 0.01) on the relative transcript profiles of these genes.

Higher relative transcript levels of several genes in the N assimilation pathways (ammonium transporters, nasAnirDB, glnA, glnK, glnLG) and the denitrification genes nor2 and narK12f were found in the Original Biofilm and the 0%NaCl biofilm (Fig. 4). The regulatory genes, fnr, ntrXY and narXL, and the denitrification genes nos, nirK and nar2 showed higher relative transcript levels in the PlkO and PlkON (oxic) cultures. The PCoA pattern suggests that the relative expression profiles are similar in the PlkAN (anoxic planktonic) cultures and the RR biofilm cultures. For instances, the denitrification genes nar1, nor1 and its associated regulatory genes norRE had higher relation transcript levels in these cultures, as well as nnrS1, nsrR1 and dnrN, which are adjacent genes in strain GP59 genome, and hmp1 and nnrS2 that encode putative regulators or proteins with NO-response function.

We then examined the relative expression profiles of genes involved in the C cycle. Among the 112 genes or gene clusters identified, we selected 17 genes or gene clusters for PCoA, in which their relative transcript levels varied more than five times (significant differences inferred by EdgeR analysis with FDR <0.05) between at least two conditions (Data S5). We also included in the PCoA the relative transcript profiles of the eleven riboswitches (Ribo1 to Ribo11) identified in the genome of strain GP59 (Table 4; Data S7). Ribo1 to Ribo7 are c-di-GMP riboswitches; Ribo9 and Ribo10, cobalamin riboswitches adjacent to each other and flanked by genes involved in cobalamin synthesis; Ribo8 and Ribo11, thiamine pyrophosphate and S-adenosyl methionine riboswitches, respectively. PERMANOVA showed that the media/physiology (IO (biofilm) vs M1403 (planktonic)), the operating conditions (oxic vs anoxic), and the methanol concentrations had substantial effects (p < 0.01) on the relative transcript profiles of these genes.

Table 4:

List of riboswitches and their flanked gene.

Riboswitch number and type	Location	Gene product
Riboswitch1: c-di-GMP	82867	Thrombospondin
Riboswitch2: c-di-GMP	822774	Pyrroloquinoline quinone precursor peptide PqqA
Riboswitch3: c-di-GMP	986642	TonB-dependent receptor
Riboswitch4 c-di-GMP	1390371	VPLPA-CTERM sorting domain
Riboswitch5 c-di-GMP	1405201	PEP-CTERM domain
Riboswitch6 c-di-GMP	1406231	PEP-CTERM domain
Riboswitch7 c-di-GMP	1407172	PEP-CTERM domain
Riboswitch8: Thiamine pyrophosphate	2153463	Phosphomethylpyrimidine synthase ThiC
Riboswitch9 and 10: Cobalamin	2296147	Cobalamin synthesis gene cluster
Riboswitch11: s-adenosylmethionine	3207836	Methionine adenosyltransferase

DOI: 10.7717/peerj.20129/table-4

Notes:

Riboswitch9 and 10 are adjacent to each other. Location: from genome annotation: GenBank accession number CP021973.1. See Data S7 for complete detail.

The relative transcript levels of the four xoxF genes (encoding putative lanthanide [Ln³⁺]-dependent methanol dehydrogenases) were higher in all four RR biofilm samples (Fig. 5), as well as two genes encoding the pyruvate carboxylase (pycAB; HCO₃⁻ to oxaloacetate) and the carbonic anhydrases (cynT; dissolved CO₂ to HCO₃⁻), involved in CO₂ assimilation. The relative transcript levels of the methanol dehydrogenase gene cluster (mxa) and one of the formate dehydrogenase gene clusters (fdh2) were higher in the PlkO, PlkON (oxic) cultures. Interestingly, the relative transcript levels of 10 out of 11 riboswitches were higher in the RR biofilm samples (Fig. 5). The other riboswitch (Ribo5; upstream of an ORF encoding PEP-CTERM domain, Table 4) had higher relative transcript levels in the PlkAN (planktonic) cultures.

Functional diversity

In addition to the microbial diversity provided by the metatranscriptome data, we can infer the potential metabolic pathways expressed in the microbial taxa active in the biofilm samples. Therefore, contigs with ORF that clustered according to their taxonomic affiliation provided an estimate of the metabolic activities occurring in the associated taxa in the biofilm. Among the ninety-three identified taxa, twenty-five taxa, in which their relative transcript levels of contigs were >2% TPM-RSEM in at least one RNA sample, were examined for putative metabolic pathways expressed in these taxa (Table 5; Data S6). Genes involved in denitrification (the four reductases) and in the carbon cycle (Fig. 6) were further examined as these pathways could explain the potential of these taxa to harvest energy and to assimilate the C1 carbon for their growth and maintenance. The bacterial pathways chosen were as follows.

1-Methanol oxidation to formaldehyde by the methanol dehydrogenase (Mdh), which provides gain of reductive equivalents by reducing cytochromes for the respiratory chain (Anthony, 1982).

2-Formaldehyde dissimilation pathways to formate by formaldehyde dehydrogenases (Fad; glutathione-dependent or independent) or by the tetrahydromethanopterin (THMPT) pathway, and formate to CO₂ by formate dehydrogenases (Fdh). Gain of reductive equivalents (NADH) was obtained in these dissimilation pathways (Anthony, 1982; Vorholt, 2002).

3-Formaldehyde assimilation pathways through the ribulose monophosphate (RuMP). Associated with RuMP is the Entner–Doudoroff pathway, with carbon assimilation to pyruvate and glyceraldehyde3P, or carbon dissimilation by regenerating Ru5P with a lost of CO₂ but gain of energy (NADPH) (Anthony, 1982; Vorholt, 2002).

4-Formaldehyde assimilation pathway through the methyl transfer (spontaneous) to tetrahydrofolate (THF) generating methylene-THF (meTHF), which is in turn incorporated in the biomass by the serine pathway or putatively by the reverse glycine cleavage system (RevGCS) (Bar-Even et al., 2013). Formate can be transformed to meTHF by the THF pathway (Anthony, 1982).

5-The Calvin cycle and the ethylmalonyl-CoA pathway (Alber, 2011) were also examined as associated with the pentose pathway or the Serine pathway. These pathways can provide additional carbon in the biomass through CO₂ incorporation (Anthony, 1982).

Nineteen taxa contain transcripts associated with denitrification genes (Table 6), among which genes encoding all four reductases were found in the genera Hyphomicrobium, Marinobacter, Methylophaga, Sedimenticola and Stappia. Besides Methylophaga spp. and Hyphomicrobium spp., Paracoccus spp. was the only other taxon with transcripts associated with methanol dehydrogenase (Table 6). The formaldehyde dissimilation pathways (Fad, THMPT, Fdh; Fig. 6) were found in eighteen taxa, among which three have THMPT and Fdh, and eight have Fad and Fdh (none have the three; Table 6). Only Methylophaga spp. use the RuMP pathway (Fig. 6). The THF pathway and the serine pathway (Fig. 6) were found in 11 and 10 taxa, respectively (6 taxa contain both; Table 6), while GCS was found in all taxa. The Calvin cycle and ethylmalonyl-CoA pathway (Fig. 6) were found in six and four taxa, respectively (Table 6).

Table 6:

Transcripts associated with genes involved in denitrification, and in formaldehyde assimilation and dissimilation pathways.

			Dissimilation			Assimilation
Taxon	Denitrification	Mdh	Fad	THMPT	Fdh	THF	RuMP	Calvin	Ser	EM	ED
Aequorivita*	nir, nor, nos	-	-	-	-	+	-	-	-	-	-
Bradymonas*	nap, nir	-	+	-	+	-	-	-	+	-	-
Cand. Promineofilum	nir, nos	-	-	-	+	-	-	-	-	-	-
Cohaesibacter	nir, nor, nos	-	+	-	+	-	-	-	-	-	-
Cycloclasticus*		-	+	-	+	-	-	-	+	-	-
Dechloromarinus*	nar, nor, nos	-	+	-	+	-	-	+	+	-	-
Geofilum	nor	-	-	-	-	+	-	-	-	-	-
Hyphomicrobium	nap, nir, nor, nos	+	-	+	+	+	-	-	+	+	-
Idiomarina*	nar, nir, nor	-	+	-	-	-	-	-	-	-	-
Lentimicrobium		-	-	-	+	-	-	-	-	-	-
Leptolyngbya*		-	+	-	-	-	-	+	-	-	-
Marinicella	nir, nor	-	+	-	-	-	-	-	-	-	-
Mariniphaga		-	-	-	-	-	-	-	-	-	-
Marinobacter*	nar, nir, nor, nos	-	+	-	+	-	-	-	+	-	+
Methylophaga	nar, nir, nor, nos	+	-	+	+	+	+	+	-	-	+
Muricauda (Flagellimonas)	nor, nos	-	-	-	-	+	-	-	-	-	-
Oceanibaculum	nar	-	+	-	+	+	-	-	+	-	-
Paracoccus	nar, nir, nor, nos	+	+	-	+	+	-	+	+	+	+
Pseudodesulfovibrio*	nar	-	-	-	+	-	-	+	-	-	-
Roseovarius*	nir, nor, nos	-	+	-	+	+	-	-	+	+	-
Sedimenticola*	nar, nir, nor, nos	-	-	+	+	+	-	+	+	-	-
Spirochaeta*		-	-	-	-	-	-	-	-	-	-
Stappia*	nap, nir, nor, nos	-	+	+	+	-	-	+	+	+
Sunxiuqinia	nir, nor	-	-	-	-	+	-	-	-	-	-
Tenericutes		-	-	-	-	-	-	-	-	-	-

DOI: 10.7717/peerj.20129/table-6

Notes:

*Pathways were identified by the KEGG mapper reconstruct analysis with data from the metatranscriptomes, and by the KEGG pathways analysis with a strain representative of the taxon when available.

+: presence; -: absence.

See Fig. 6 for pathway description.

Mdh, Methanol dehydrogenase (mxa mainly); Fad, Formaldehyde to formate by the formaldehyde dehydrogenase (either glutathione-dependent or -independent); THMPT, Process of formaldehyde to formate through the Tetrahydromethanopterin (THMPT) pathway; Fdh, Formate dehydrogenase; THF, Process of formate to Methylene tetrahydrofolate (meTHF) through the THF pathway; RuMP, Formaldehyde assimilation by the Ribulose-5P pathway, Calvin, Calvin cycle here is defined as taxon containing transcripts encoding the phosphoribulokinase (Prk) and the ribulose-1,5-bisphosphate carboxylase/oxygenase (rubisco) enzymes; Ser, Formaldehyde assimilation by the serine pathway; EM, Ethylmalonyl-CoA pathway; All taxa have complete pathways of glycine cleavage system (GCS), tricarboxylic acid cycle (TCA) and glycolysis; ED, Entner–Doudoroff pathway.

C/N and the denitrification performance

Our results showed that different C/N did not have a significant impact on the reactor denitrification performance (NO_x/NO₃⁻ reduction rates, NO₂⁻ and N₂O dynamics), but did have an impact on the methanol consumption rates and the CO₂ production rates. This impact may be linked to the initial concentrations of methanol (i-Methanol) applied to the reactor, as these concentrations also correlated with the methanol consumption rates and the CO₂ production rates. Interestingly, the initial concentrations of NO₃⁻ (i-NO₃⁻) did positively correlate with the denitrification performance, but such correlation was not observed with the i-Methanol. Our results suggest that opposite effects occur between the i-NO₃⁻ and the i-Methanol applied to the reactor, resulting in a null effect by C/N on the denitrification performance, but not on the carbon metabolism.

In methanol-fed denitrification processes, methylotrophic bacteria use methanol as source of carbon and energy, where it is first converted to formaldehyde by the methanol dehydrogenase with transfer of electrons to cytochromes and then to the respiratory chain. Formaldehyde can either be assimilated into the biomass, or dissimilated in formate and CO₂ by enzymatic reactions generating NADH (Anthony, 1982). Under denitrifying conditions (anoxic), NO₃⁻ and N oxides replace O₂ as terminal electron acceptors for respiration. NO₃⁻ acts as an electron acceptor with the NO₃⁻ reductase system (Nar system) contributing to the proton motive force, while NO₂⁻, NO and N₂O accept electrons from electron pools generated from the respiratory chain (Simon, 2011). Our results suggest that the level of NO₃⁻ reached a saturation point where the denitrification performance could not increase. This was not the case with methanol, where it can contribute to both energy and cell growth and maintenance. Therefore, as the methanol concentration increased, higher methanol consumption (input) and CO₂ production (output) occurred without the need for greater denitrifying activities.

The reactor was able to reduce most, if not all, all i-NO₃⁻ applied to the reactor after four operating days, although the higher concentrations (32.1 and 42.8 mM) generated higher levels of residual NO₂⁻, which would be eventually reduced over time. The reactor did not generate ammonium, and negligible amount of N₂O was found at the end of the reaction phase. However, high proportions of residual methanol (47 to 91%) were found at the end of the reaction phase, which is not a desired outcome for a denitrification system. The release of methanol in the effluent of a water treatment can be harmful to the ecosystem. For example, Pungrasmi, Playchoom & Powtongsook (2013) showed that applying a C/N of 4 with methanol in their denitrification system resulted in the best denitrification performance. However, due to the concentration of residual methanol exceeded toxic level for fish development, a C/N of 3.3 had to be applied in their system when treating the effluent from a recirculation (tilapia) aquaculture system (RAS).

We acknowledge that analysing a single reactor limited the significance of our findings. Running multiple reactors in parallel is challenging to maintain synchronization of conditions and can be a burden on human resources, as the reactors operate continuously and require daily attention. The initial strategy was to operate all eight conditions three times (three periods) with the same biofilm, but due to the COVID crisis, only two could be performed. We noticed however that, as the biofilm ages, the microbial community evolved from the early stages, and the biofilm thickens with active and dead zones, increasing the variability of the operating conditions. Ideally, each condition should have been conducted with fresh acclimated biofilm in triplicate (i.e., eight reactors with fresh acclimated biofilm, in triplicate).

Evolution of the bacterial community

Both approaches (16S amplicon sequencing and metatranscriptomes) to monitor the evolution of the bacterial diversity profiles in the biofilm during the operating conditions showed that the bacterial community established in the reactor diverged substantially from the Original Biofilm collected in the Biodome denitrification system in the proportion of the common taxa. These changes can be explained by the operating conditions: laboratory-scale sequential batch-mode recirculating reactor versus full-scale continuous operating mode in the Biodome denitrification system. Although both processes used the same medium (IO), the input water in the Biodome system contained additional matter from animal and plant waste, or non-ingested food (Parent, Labbé & Villemur, 2002), and additional microorganisms from the aquatic ecosystem (e.g., protozoa) (Laurin et al., 2008) that might have influenced its microbial community.

The results obtained with 16S amplicon sequencing revealed that a high proportion of Methylophaga spp. established in the recirculating reactor during the acclimation phase (D0), and that subsequently the biofilm populations increased in diversity at the expense of Methylophaga spp. (although still constituting a high proportion of the bacterial community), in favor of heterotrophic bacteria such as the genera Stappia, Marinobacter, Marinicella, Oceanibaculum and Paracoccus. Methylophaga spp., Hyphomicrobium spp., Paracoccus spp., Marinobacter spp. and Marinicella spp. have been found in other marine denitrification systems (Güven, 2009; Lu, Chandran & Stensel, 2014; Furukawa et al., 2016; Wang et al., 2023). Some members of these taxa have denitrification capacity such as Stappia spp., Paracoccus spp. and Marinobacter spp. Contrary to what was stated by Duc et al. (2018) and echoed by other reports about Shoji et al. (2014), the sulfur oxidizing bacteria Marinicella were not shown to carry NO₃⁻ reduction to N₂ in this report. Results from metatranscriptomes, however, revealed that nirK and nor genes were found to be affiliated to Marinicella spp. None of the Marinicella sp. sequenced genomes in GenBank contain the full set of denitrification reductases. Stappia spp. are marine chemoorganotrophic bacteria that can oxidize CO, and some species contain a gene for the large subunit of ribulose-1,5-bisphosphate carboxylase/oxygenase (RuBisCO) (cbbL) that may suggest the ability to couple CO utilization to CO₂ fixation (Weber & King, 2007).

Potential metabolic pathways enabling the transformation of formaldehyde for carbon assimilation or dissimilation (Anthony, 1982; Vorholt, 2002; Yurimoto, Kato & Sakai, 2005) could explain the persistence and growth of heterotrophs in the reactor, as these pathways were found in the observed taxa. For instances, the oxidation of formaldehyde and formate generates NADH by their respective dehydrogenases. Formaldehyde is produced by the action of the methanol dehydrogenase, carried by Methylophaga spp., Paracoccus spp. and Hyphomicrobium spp., at the periplasm (Anthony, 2004). Due to high proportion of Methylophaga spp. and high level of methanol, it is possible that excess of formaldehyde diffused into the biofilm and was available to the other organisms because of their promiscuity (Chistoserdova & Kalyuzhnaya, 2018). The emergence of these heterotrophic bacteria could also be the result of subsequent colonization of the biofilm after the early colonization of the supports by the methylotrophs such as Methylophaga spp., providing complex carbon to the heterotrophs such as carbohydrates in the extracellular polymeric substances of the biofilm matrix or organic waste from dead cells as the biofilm aged.

Relative expression profile of strain GP59

M. nitratireducenticrescens strain GP59 has adapted to the recirculating reactor resulting in colonization of supports in high proportion. Following the early stage of colonization, heterotrophs increased in proportion during the operating conditions. Regarding the relative transcript levels of strain GP59, these were influenced by the physiology state (biofilm vs planktonic), the culture media, and the operating conditions (oxic vs anoxic).

High proportions (22.5 to 68.9%) of 16S rRNA gene sequences affiliated with Methylophaga spp. were found in all RR biofilm samples, six to 18 times higher than in the Original Biofilm (3.7%). However, despite this low level in the Original Biofilm, transcript reads associated with strain GP59 and other Methylophaga species in the Original Biofilm together accounted for approximately 50% of the overall relative transcript profile of the Original Biofilm, suggesting that Methylophaga spp. were more active in the Original Biofilm than the other taxa. For example, transcript reads associated with Hyphomicrobium spp. in the Original Biofilm RNA sample represented <1% relative transcript levels, whereas it represented >40% in 16S rRNA gene sequences in the DNA sample. Although relative transcript levels of strain GP59 had 4-fold increases on average in the RR biofilm RNA samples, these levels associated with the other Methylophaga species were almost negligible, suggesting that strain GP59 has better adapted to the anoxic recirculating conditions than the other Methylophaga species. The denitrification system at the Montreal Biodome, although containing a deoxygenation reactor, was not entirely airtight and might have sufficient O₂ (0.3–0.8 mg O₂/L; Labbé et al., 2003) for the development of obligate aerobic Methylophaga sp. such as M. frappieri strain JAM7 that we isolated from this system (Auclair et al., 2010).

The relative gene expression profiles of strain GP59 in the RR biofilm were clearly different from those observed in (i) the Original Biofilm, (ii) the ASW-medium, batch-static biofilm cultures, and (iii) the planktonic cultures with the specific Methylophaga medium (Methylophaga-1403). These results strongly suggest that the medium, the physiological conditions (planktonic, biofilm) and the operating conditions (anoxic, oxic) have a significant impact on the overall gene expression profiles in strain GP59 (Di Capua et al., 2022).

The relative expression profiles of genes involved in N-cycle in strain GP59 had a similar pattern in the four RR biofilm samples and in the PlkAN cultures, suggesting that the operating conditions (planktonic cultures versus recirculating biofilm reactor) did not modify the expression profiles of strain GP59 regarding the N-cycle despite changes in C/N in the recirculating reactor. The relative expression patterns of the gene clusters of the nar1 and nor1 systems and genes involved in NO response appear to have been more engaged under these conditions, unlike the nar2 and nor2 systems, which were more involved under the oxic conditions or in the Original Biofilm. We demonstrated in a previous study (Lestin & Villemur, 2024) that nar1 (Nar: NO₃⁻ reductase) had higher expression levels in strain GP59 planktonic cultures than those of nar2, suggesting higher denitrifying activities in the RR biofilm and the PlkAN cultures. These higher activities would generate higher levels of NO, explaining the higher relative transcript levels of genes involved in NO response.

All four xoxF genes had higher relative transcript levels in the four RR biofilm samples. The XoxF-type lanthanide (Ln³⁺)-dependent methanol dehydrogenases are widespread among bacteria. Ln³⁺ is an important co-factor that regulates the expression of the xoxF and the mxa gene cluster (encoding Ca²⁺-dependent methanol dehydrogenase) (Chistoserdova & Kalyuzhnaya, 2018). We do not know whether Ln³⁺ was included in the IO commercial formulation, or whether it exists naturally as trace element. However, the fact that the xoxF genes had higher relative transcript levels in the RR biofilm samples in contrast to those of mxa which were higher in the oxic planktonic cultures, suggests important roles of XoxF in the carbon dynamics in biofilm under anoxic conditions as alternative methanol dehydrogenases (Chistoserdova & Kalyuzhnaya, 2018). Mustakhimov et al. (2013) proposed a connection between XoxF in providing electrons from methanol oxidation and the denitrification pathway in the methylotrophic Methylotenera mobilis JLW8.

Along with the xoxF expression patterns, similar results were observed with 10 of the 11 riboswitches with higher relative transcript levels in the RR biofilm samples. Examination of the genes flanking the eleven riboswitches showed diverse functions. Four of them (Ribo4, 5, 6, 7; c-di-GMP type) are located in close proximity (18 kb apart), and each flanked a gene encoding a protein with PEP-CTERM domain, which appears to be involved in exopolysaccharide-associated protein sorting systems (Haft et al., 2006; Villemur et al., 2019). Furthermore, genes surrounding these coding sequences encode proteins involved in exopolysaccharide biosynthesis and export. In Gram-negative bacteria, many synthase-dependent exopolysaccharide secretion systems are post-translationally regulated by an inner-membrane c-di-GMP receptor (Whitney & Howell, 2013). Interestingly, Ribo5 had higher relative transcript levels in planktonic anoxic cultures (PlkAN), while the three others had higher levels in the RR biofilm samples. These results suggest that specific controls of these pathways affecting the bacterial surface are important for the fate of the bacterial physiology (biofilm vs planktonic). It is known that c-di-GMP is a key factor in bacteria to regulate several genes involved in the planktonic-biofilm transition. One of the mechanisms of biofilm development is the surface attachment, which requires specific adherent organelles such flagella, fimbriae or pili (Martínez & Vadyvaloo, 2014).