Microbial activity response to hydrogen injection in thermophilic anaerobic digesters revealed by genome-centric metatranscriptomics

Background The expansion of renewable energy produced by windmills and photovoltaic panels has generated a considerable electricity surplus, which can be utilized in water electrolysis systems for hydrogen production. The resulting hydrogen can then be funneled to anaerobic digesters for biogas upgrading (biomethanation) purposes (power-to-methane) or to produce high value-added compounds such as short-chain fatty acids (power-to-chemicals). Genome-centric metagenomics and metatranscriptomic analyses were performed to better understand the metabolic dynamics associated with H2 injection in two different configurations of anaerobic digesters treating acidic wastes, specifically cheese manufacturing byproducts. These approaches revealed the key-genes involved in methanation and carbon fixation pathways at species level. Results The biogas upgrading process in the single-stage configuration increased the CH4 content by 7%. The dominant methanogenic species responsible for the upregulation of the hydrogenotrophic pathway in this reactor was Methanothermobacter wolfeii UC0008. In the two-stage configuration, H2 injection induced an upregulation of CO2 fixation pathways producing short-chain fatty acids, mainly acetate and butyrate. In this configuration, the abundant species Anaerobaculum hydrogeniformans UC0046 and Defluviitoga tunisiensis UC0050 primarily upregulated genes related to electron transport chains, suggesting putative syntrophisms with hydrogen scavenger microbes. Interestingly, Tepidanaerobacter acetatoxydans UC0018 did not act as an acetate-oxidizer in either reactor configurations, and instead regulated pathways involved in acetate production and uptake. A putative syntrophic association between Coprothermobacter proteolyticus UC0011 and M. wolfeii UC0008 was proposed in the two-stage reactor. In order to support the transcriptomic findings regarding the hydrogen utilization routes, an advanced bioconversion model was adapted for the simulation of the single- and two-stage reactor setups. Conclusions This is the first study investigating biogas reactor metatranscriptome dynamics following hydrogen injection for biomethanation and carbon fixation to short-chain fatty acids purposes. The same microbes showed different patterns of metabolic regulation in the two reactor configurations. It was observed an effect of the specialized acidogenic reactor on the overall microbial consortium composition and activity in the two-stage digester. There were also suggested the main species responsible for methanation, short-chain fatty acids production, and electron transport chain mechanisms, in both reactor configurations. Electronic supplementary material The online version of this article (10.1186/s40168-018-0583-4) contains supplementary material, which is available to authorized users.


Background
The increasing demand for renewable energy sources (RES) promoted the expanded use of windmills and photovoltaic panels. The installed global wind capacity increased by 10.8% in 2017, with China and the USA as the major producers of electricity from wind [1]. This expansion has led to the production of a considerable electricity surplus, which cannot be easily stored in batteries due to high cost or injected into the national grid, since it could cause electrical instabilities. This electricity surplus could be used for H 2 production via water electrolysis [2]. Nevertheless, hydrogen is highly volatile and therefore difficult to store and transport, and is associated with some environmental risks. Alternatively, this H 2 could be used for biogas upgrading purposes (power-to-methane) or for the production of high value-added compounds such as fatty acids (mainly short-chain carboxylates) and alcohols (power-to-chemicals), via anaerobic digestion (AD), generating an energy gain in the form of methane along with organic waste valorization [3][4][5][6].
The process of anaerobic digestion relies on a complex microbial syntrophic chain, which degrades the organic matter into various byproducts such as short-chain fatty acids (SCFAs) and eventually to biogas (i.e., methane and carbon dioxide). Different operating conditions can direct the process toward higher CH 4 or SCFAs yields. Specifically, pH, temperature, and hydraulic retention time (HRT) represent the most important factors in SCFAs accumulation [7]. SCFAs, also known as volatile fatty acids (VFAs), consist of six or fewer carbon atoms and can be used in the synthesis of a wide range of compounds, including biosurfactants, bioflocculants, and bioplastics (polyhydroxyalkanoates (PHAs)) [8,9].
SCFAs are among the main products of cheese whey AD and, despite the low alkalinity typical of this waste, Fontana et al. have recently demonstrated the higher efficiency of a two-stage continuous stirred tank reactor (CSTR) than a single-stage configuration for treating whey permeate and hard-cheese solid waste to produce biogas [10]. However, the different reactor configurations evaluated in the study also indicated the potential to use cheese wastes AD for carbon dioxide fixation to SCFAs, such as butyrate and acetate.
The conversion of these chemical compounds is performed by an intricate set of microbes where different species cooperate or compete to generate the final product. Functional properties of individual species can be explored by reconstructing their genomes from the metagenome. This metagenomic approach allows the identification of the so-called metagenome assembled genomes (MAGs), which can be successfully reconstructed by binning the scaffolds that were previously assembled from the shotgun reads. Among the existing binning strategies [11,12], Campanaro et al. developed a specific methodology that was applied to characterize microbial communities in biogas reactors [13]. This genome-centric approach, combined with the metatranscriptomics, provides a powerful method to uncover the phylogenetic and metabolic properties in anaerobic digestion without relying on culture-dependent techniques. Particularly, the possibility to align the RNA-seq reads to the MAGs allows to investigate the activity of specific microbes and to monitor the modifications occurring during the alteration of the environmental conditions.
The current study aimed to investigate the effect of exogenous H 2 injection on microbial activity in two different thermophilic reactor configurations (single and two-stage continuous stirred tank reactors) treating cheese wastes. Particularly, improvements in biogas methanation and/or CO 2 fixation to SCFAs by CO 2 reduction pathways have been evaluated. Genome-centric metagenomics and metatranscriptomics allowed an in-depth analysis of differential gene expression following H 2 addition by the most abundant MAGs. The outcomes of the present study also reveal putative syntrophic associations between key microbial groups and methanogenic archaea.

Biogas reactors' configuration
The experiment was carried out in a single and two-stage CSTRs, denoted as R1 and R2-R3, respectively, as described by Fontana et al. [10]. The reactors had a total working volume of 3 L; they were continuously stirred at 150 rpm using magnetic stirrers and equipped with thermal jackets to maintain the operating temperature at 55 ± 1°C. Initially, all reactors were inoculated with thermophilic inoculum obtained from Snertinge biogas plant, Denmark. The experiment was divided in two phases: phase I, when the reactors were fed exclusively with cheese whey permeate and cheese waste powder, and phase II, when exogenous H 2 was added to the reactors (Fig. 1). The H 2 gas was provided in the reactors R1 and R2 by using a peristaltic pump operating at a flow rate of 1.7 mL/min and was injected through Al 2 O 3 ceramic membranes having pore size of 1.2 μm. For the two-stage CSTR system, the output gas and the liquid effluent from the acidogenic reactor (R2) was transferred pneumatically to the methanogenic reactor (R3) via a connecting tube (Fig. 1). In order to maximize the gas-liquid mass transfer, another peristaltic pump was connected with each reactor to recirculate the output gas. The recirculation flow rate was set to 0.7 mL/L min. During both experimental phases, the organic loading rate (OLR) was set at 2.4 g COD/L day, and thus, the hydraulic retention time (HRT) was maintained at 15 days (split in 3 and 12 days for R2 and R3, respectively). The characteristics of the feedstock are described in Fontana et al. [10]. The influent feedstock was automatically provided four times per day using controlled peristaltic pumps. Sodium hydroxide addition in R1 and R3 was applied whenever the pH dropped below 6.5. The H 2 flow rate was defined in relation to the stoichiometry of hydrogenotrophic methanogenesis reaction (4 H 2 /1 CO 2 mol/mol) as described by Bassani et al. [14]. Half H 2 volume (~820 mL/L day) was injected to gradually adapt the system to the new condition. Mass balance calculations were carried out as indicated in the Additional file 1.
The H 2 utilization efficiency was calculated as previously described [15], using the following Eq. (1): The CO 2 conversion efficiency was calculated as follows (2) Sample collection For metatranscriptomic analyses, triplicate samples (~30 mL each) were collected from R1 and R3 at steady-state reactor operation of phase I (i.e., period with stable biogas production with a daily variation lower than 10% for at least 5 days), and after 1 week from the H 2 injection (phase II). Replicates obtained from phase I were indicated as R1_a, R1_b, R1_c, R3_a, R3_b, R3_c, while replicates from phase II were denoted as R1H_a, R1H_b, R1H_c, R3H_a, R3H_b, R3H_c.

DNA and RNA extraction and sequencing
Samples were centrifuged at 10,000 rpm for 10 min and the supernatant was discarded leaving~1 g of pellet. To avoid RNA degradation, 3.5 mL of phenol/chloroform (pH 6.7/8.0) premixed with isoamyl alcohol (25:24:1) (Amresco, Incorporated) was added to the pellet after centrifugation and the samples were immediately processed for RNA extraction, as previously reported [16].  The abundance shifts of the species (MAGs) and the differential gene expression between the two phases were evaluated. MAG metagenome assembled genome, R1 single-stage reactor, R2 acidogenic reactor of the two-stage, R3 methanogenic reactor of the two-stage integrity was determined using agarose gel electrophoresis. Samples were sequenced, using NextSeq 500 sequencing technology and Nextera XT kit (Illumina, San Diego, CA) for library preparation (150 + 150 bp). The quality and the quantity of the extracted DNA and RNA were also determined using both NanoDrop (ThermoFisher Scientific, Waltham, MA, USA) and Qubit fluorometer (Life Technologies, Carlsbad, CA, USA).
Additional liquid samples (4 mL) for metagenomic analyses were collected in three replicates and in three time points (days 41, 52, and 61) at steady-state condition from R1 and R3. DNeasy PowerSoil Kit (QIAGEN, Germany) was used for genomic DNA extraction with minor modifications (purification with 1 mL of Phe:Chl:IAA pH 8, Sigma-Aldrich, DK). Microbial community composition was determined using the V3-V4 hypervariable regions of 16S rRNA gene using universal primers Pro341F and Pro805R [18]. Amplicon preparation and sequencing were performed at BMR Genomics S.r.l. (Padua, Italy) using Illumina MiSeq platform. Raw sequencing data were submitted to the sequence read archive database (SRA) of NCBI under the BioProject PRJNA490620 and the accessions SAMN10054307-SAMN10054315 for R1 samples and SAMN10054316-SAMN10054324 for R3 samples. Data analysis was performed using CLC Workbench software (V.8.0.2) with microbial genomics module plug in (QIAGEN Bioinformatics, Germany) as previously described [5] with "quality limit" parameter set at 0.01. In brief, chimeras filtering, operative taxonomical units clustering, taxonomic assignment (with Greengenes v13_5 database), and alpha and beta diversity (Unweighted UniFrac) calculation were done using standard parameters. Additionally, relative abundances of the OTUs were used to determine the beta diversity value (Whittaker method) and all the possible comparisons between replicates, time points, and reactors were performed.

Reads alignment, gene expression calculation, and statistical analysis
Gene expression analysis was performed considering as reference the global metagenome assembly. Gene finding and annotation were reported in a previous study [10] and in brief were performed as follows: gene prediction on the entire assembly was determined with Prodigal (v2.6.2) run in metagenomic mode [19]. Protein-encoding genes were annotated using reverse-position specific BLAST algorithm and using Clusters of orthologous groups COG and Pfam database [20,21]; only results with e value lower than 1e −5 were retained. Genes were also annotated according to Kyoto Encyclopedia of Genes and Genomes (KEGG) using GhostKOALA [22] and to EggNOG 4.5.1 using eggNOG-mapper [23]. Filtered reads were aligned on reference metagenome assembly using bowtie2 (v2.2.4) [24] and the number of reads mapped per each gene was determined from the "sam" file using HTSeq (v0.6.1) [25] with the options "-count" and "intersection-non empty." Each gene was previously assigned to the correspondent MAG with a binning strategy previously described [10]. MAGs abundance was considered directly related to scaffold coverage, which was determined by aligning the shotgun reads on the assembly as described by Campanaro and co-workers [26]. Coverage values determined for MAGs were visualized with MeV [27]. To evaluate the changes in abundance of the main MAGs after H 2 injection (phase II), a comparison with the coverage values of the previously described MAGs, before H 2 injection (phase I) [10], has been carried out. The statistical analysis was performed independently for each MAG using edgeR software [28] and the differentially expressed genes were filtered considering the p value (pVal.Tgw < 0.05) and the coverage ratio (> 2-fold change). KEGG pathway maps were obtained with the "KEGG Mapper Search&Color Pathway" tool [29]. To identify the COG and KEGG functional classes statistically enriched of differentially expressed genes, the procedure described by Treu et al. [30] was followed. Briefly, 10,000 random samplings of n genes (where n is the number of genes for each COG or KEGG class) were performed on the entire data set of genes expressed (data set S3) using a Perl script implementing the "rand()" function. Assuming differentially expressed (DE) as the number of differentially expressed genes in a group of n genes, the fraction of randomly selected samples having differentially expressed genes equal to or higher than DE was calculated. If this fraction was lower than the significance level (0.05), the enrichment of genes differentially expressed in the n genes was considered significant. Additional statistical analyses have been performed on a selection of MAGs; Fisher's exact test was applied to define the KEGG pathways including a significant number of differentially expressed genes. Finally, multiple test correction was performed for calculating false discovery rate (data set S4).
Canonical correspondence analysis (CCA) was performed using the R functions implemented in VEGAN v2.4-4 [31], while correspondence analysis (CA) based on Pearson calculation and non-metric multidimensional scaling (NMDS) were performed with R following the procedure described by Torondel et al. [32].

Modeling of the experimental setup
Software simulations of the single and two-stage reactors were carried out using the complex anaerobic digestion simulation suite (BioModel) developed by Angelidaki and co-workers [33,34], which was later extended to account for biogas upgrading experiments with hydrogen injection [35]. For improved simulation accuracy, previously optimized model parameters were taken from Kovalovszki et al. [36]. The tool was further enabled to simulate two-stage reactor configurations, through the inclusion of intermediate compounds and gases produced in the acidogenic reactor (R2) in the feed of the methanogenic reactor (R3). In addition, syntrophic acetate oxidizing microbial groups were also considered in the model, making it possible to simulate the dynamic interactions between acetoclastic and hydrogenotrophic methanogenic groups.

Results and discussion
Metagenomic and metatranscriptomic investigations were performed at two time points; the first point referred to the reactors' steady-state performance before H 2 injection (phase I) and the second occurred 1 week after H 2 injection (phase II). To verify the stability of the microbial community during the reactors' stable operation, an additional set of metagenomic samples was collected from R1 and R3 at multiple time points and was analyzed using 16S rRNA gene amplicon sequencing. The overview of the sequencing depth obtained with the different NGS data type showed that the microbial community under consideration was well captured (Additional file 1: Table S1). OTUs taxonomy showed that the biological process was adequately captured in terms of microbial composition; results from beta diversity demonstrated that the microbial community was stable during time, and thus, the selected point chosen for in depth analysis was representative of the steady-state period (Additional file 1: Table S2 and Figure S1). In particular, the overall OTUs' taxonomic distribution in both R1 and R3 was in agreement with the profile of the reconstructed MAGs (Additional file 1: Table S2). PCoA results and OTUs relative abundances (i.e., 1.5 average fold change) revealed negligible variations among the different time points in R1 (Additional file 1: Table S2 and Figure S2). Regarding the reactor R3, the dominant OTUs abundances were coherent with MAGs coverages (Fig. 2a and Additional file 1: Table S2). The observed differences in the PCoA results were mainly attributed to the microbial diversity of a minor subset of OTUs in the middle sampling point. Considering the results from the biochemical parameters, the reactor operation was stable, indicating that this OTUs subset was not primarily involved in the methanation process.
The reconstructed MAGs identified in the microbial community represented more than 60% of the entire microbiome. Therefore, the results from the current work covered successfully the majority of the transcriptional changes occurring in the reactors excluding only a minor fraction of the information present in the shotgun reads. Moreover, the total number of protein-encoding genes identified in the assembly was slightly higher than 196,000, out of which 80% had at least 1 read in 1 of the samples examined and 27% had 10 or more reads. Consequently, the outcomes of the identified genes confirm that the transcriptional study was representing the expression of a considerable fraction of the total genes in the microbiome.
In order to acquire a global overview, analysis of the total expressed genes (not assigned to MAGs) was carried out considering COG classification (Additional file 1: Figure S3, Additional file 2: Dataset S3). In both reactor configurations, the most differentially expressed categories (excluding R and S categories, representing the general and unknown functions, respectively) belonged to the carbohydrate and amino acid transports and metabolisms. However, a high fraction of genes within the C category (energy production and conversion) was also differentially expressed in both single and two-stage reactors. Analysis of the expressed genes was subsequently performed in a genome-centric perspective to decipher the roles of the individual MAGs. The investigation was focused on the most abundant and active species, having more than 1000 expressed genes after H 2 injection. However, the analysis was exceptionally expanded to two MAGs (Methanothermobacter wolfeii UC0008 and Tepidanaerobacter acetatoxydans UC0018) that were considered of particular interest, despite the fact that they showed less than 1000 expressed genes.

Single-stage reactor: power-to-methane
The single-stage reactor (R1) exhibited a pH trend ranging between 6.3 and 7.3 during phase I (Fig. 2b). Total VFAs were highly concentrated (9.7 ± 1.1 g/L) and composed mainly of acetate (6.1 ± 1.0 g/L) (Fig. 2b). These conditions inhibited the activity of methanogenic archaea, resulting in a CH 4 yield equal to 31% of the theoretical value, which is 350 mL CH 4 /g COD (Fig. 2a and Table 1). High VFAs concentrations lower the pH of the reactor, and thus, concomitantly lead to alteration of the microbial activities [37]. This effect is especially evident during the anaerobic digestion of acidic substrates characterized by poor buffering capacity [10]. In particular, methanogens are the most sensitive species to over acidification events, since their optimal growth rate ranges between the pH values of 6.5 and 8.5 [38]. Total VFAs increased by~1 g/L in phase II, mainly due to higher butyrate concentration (Fig. 2b). This increase could be caused by the high acetate levels present in this reactor (6.7 ± 0.8 g/L), which may have hampered syntrophic butyrate oxidation [39]. Despite the further over acidification during phase II, the CH 4 yield in R1 increased by 10% compared to the previous experimental phase (Fig. 2a and Table 1), indicating a positive effect of H 2 injection on the methanogenic consortia.
Three MAGs were identified as dominant (77% of the microbial community) in R1 during phase I, specifically Coprothermobacter proteolyticus UC0011, Anaerobaculum hydrogeniformans UC0046, and Defluviitoga tunisiensis UC0050 (Fig. 3, Additional file 1: Figure S4, Additional file 2: Dataset S1). This microbial core reached 85% of relative abundance after H 2 injection, with C. proteolyticus UC0011 as the dominant species (61% relative abundance) (Fig. 3, Additional file 1: Figure S4, Additional file 2: Dataset S1). These results highlight the strong microbial selection The orange and green arrows highlight the DNA/RNA sampling points for the single and the two-stage configuration, respectively operated by both the feed characteristics such as the acidic pH and low buffer capacity, and the increased H 2 partial pressure inside the reactor. A significative correlation between C. proteolyticus UC0011 in phase II and H 2 content in the reactor was also highlighted by statistical analysis (Additional file 1: Figure S5). Transcriptional data showed that C. proteolyticus UC0011 responded to H 2 addition by differentially expressing genes related to carbon metabolism, specifically the pyruvate metabolic pathway ( Table 2 and Fig. 4). Genes associated with the pyruvate dehydrogenase complex and pyruvate-formate lyase (Aco, Ace, and Pfl), both involved in acetyl-CoA production, increased their expression by~3-fold in C. proteolyticus UC0011 (Fig. 5 and Additional file 2: Dataset S2). This upregulation suggests that C. proteolyticus UC0011 is involved in the acetate accumulation observed in R1 (Fig. 4 and Fig. 5). In contrast, expression of the ATP-dependent protease Clp was inhibited by~4-fold in C. proteolyticus UC0011 (Fig. 4 and Additional file 2: Dataset S2), indicating a specific repression of the proteolytic activity of this enzyme, which causes H 2 release [40,41].
Analysis of A. hydrogeniformans UC0046 revealed the differential expression of genes encoding ABC transporters related to amino acid translocation across the plasma membrane (Table 2 and Fig. 4). Expression of acetyl-CoA   [42]. Moreover, the gene coding for Buk enzyme, which catalyzes the final step for butyrate formation, is frequently used as biomarker for the identification of butyrate-producing communities [43]. Therefore, these results indicate that A. hydrogeniformans UC0046 contributes to the increased butyrate concentration found in R1 after H 2 injection. Only 11 genes of D. tunisiensis UC0050 were differentially expressed after H 2 addition, and 3 were involved in quorum sensing activities (Table 2). RpoD (sigma 70) (Additional file 2: Dataset S2) is the main bacterial sigma factor responsible for housekeeping gene transcription [44], and showed decreased expression. This regulation pattern suggests an inhibition of basal gene expression in this species during phase II.  The known syntrophic acetate-oxidizing bacterium (SAOB) T. acetatoxydans UC0018 [41] slightly decreased in abundance after H 2 addition (Fig. 3, Additional file 1: Figure S4, Additional file 2: Dataset S1). Transcriptomic data showed that T. acetatoxydans UC0018 differentially expressed genes encoding ABC transporters and enzymes involved in amino acid and sugar metabolism ( Table 2 and Fig. 4). Sugar intake was decreased by the downregulation of specific ABC transporters (Mgl permease) as well as genes related to glucose and galactose metabolism (nag sugar kinase, fruK, fba) by 4-and 8-fold, respectively (Additional file 2: Dataset S2). This regulation suggests the existence of a "feedback mechanism" to limit excessive acetate production via sugar catabolism.
Regarding the methanogenic consortia, there was a clear dominance of one archaeal species, the hydrogenotrophic M. wolfeii UC0008, which was reduced in abundance by half after H 2 injection (Fig. 3, Additional file 1: Figure S4, Additional file 2: Dataset S1). A significative reduction was also observed for the less abundant hydrogenotrophic Methanothermobacter thermautotrophicus UC0010 (Fig. 3, Additional file 1: Figure S4, Additional file 2: Dataset S1). The growth inhibition of Archaea may be partly due to the low alkalinity intrinsic to cheese whey permeate, along with the increased acidification of the system following H 2 addition (Fig. 2b) [38]. Despite this inhibition, H 2 injection induced a significant upregulation of the hydrogenotrophic methanogenesis pathway in M. wolfeii UC0008 (Additional file 1: Figure S6). Such transcriptional behavior led to increased CH 4 content in the biogas, with a~54% CO 2 conversion efficiency (Table 1).
The high accumulation of acetate measured in this reactor could also be related to the lack of aceticlastic methanogens, whose growth was probably not favored by the conditions established in the single-stage configuration. In addition to the difficulty in maintaining the pH in a proper range for methanogenic growth, toxicity related to the accumulation of cations (i.e., potassium) or lipids has been previously hypothesized [10].

Two-stage configuration, acidogenic reactor: power-tochemicals
The hypothesis for applying the H 2 only into the acidogenic reactor (R2) of the two-stage configuration was that it could better withstand a potential pH increment that could be caused by the transformation of CO 2 into methane. This reactor indeed maintained a stable pH (~4) throughout the process. The main VFA produced in R2 was butyrate (3.9 ± 0.7 g/L), which increased by~1 g/L after H 2 injection (Fig. 2c). Bifidobacterium crudilactis UC0001 was the dominant species inhabiting this reactor, and it showed a change in abundance after H 2 addition, decreasing from 82 to 52% of the total microbiome (Fig. 3, Additional file 1: Figure S4, Additional file 2: Dataset S1). In contrast, the heterofermentative lactic acid bacteria Leuconostoc pseudomesenteroides UC0016 strongly increased in abundance (~3-fold). This variation could be related to the higher butyrate concentration present in R2 during phase II (Fig. 2c), since the lactose fermentation to lactate by L. pseudomesenteroides UC0016 can enhance the cross-feeding of Clostridiales species involved in the conversion of lactic acid to butyrate [45]. It was indeed observed a~4-fold increase of Clostridiaceae sp. UC0025, Clostridiaceae sp. UC0028, and Clostridium sp. UC0030 during phase II (Fig. 3, Additional file 1: Figure S4, Additional file 2: Dataset S1). Moreover, butyrate production by clostridial-type fermentation is also known to be favored under high H 2 partial pressures [46][47][48][49][50], as can occur during exogenous H 2 injection in R2. The significative effect of butyrate increase in phase II on microbial distribution was also evidenced by statistical analysis (Additional file 1: Figure S5).

Two-stage configuration, methanogenic reactor: powerto-chemicals
The methanogenic reactor of the two-stage configuration (R3) maintained the pH between 6.7 and 7.5 during phase I, exhibiting lower accumulation of total VFAs (primarily acetate) than the single-stage configuration (3.4 ± 1.3 g/L) (Fig. 2d). These operating conditions resulted in a CH 4 yield equal to 80% of the theoretical value (Fig. 2a, d). However, H 2 addition induced an increase in total VFAs (primarily butyrate), which doubled in concentration to 6.1 ± 0.3 g/L (Fig. 2d). The CH 4 yield was highly reduced under these conditions ( Fig. 2a and Table 1), and the increased butyrate and acetate levels along with the decreased CH 4 content seen in phase II indicate that CO 2 fixation toward SCFAs overtook the methanation pathways. The most abundant MAGs were the same as those found in the single stage R1, specifically C. proteolyticus UC0011, A. hydrogeniformans UC0046, and D. tunisiensis UC0050, which accounted for 47% of the microbiome in the reactor (Fig. 3, Additional file 1: Figure S4, Additional file 2: Dataset S1). This microbial core reached 81% of relative abundance after H 2 injection, and D. tunisiensis UC0050 was the dominant species (54% of the total community) (Fig. 3, Additional file 1: Figure S4, Additional file 2: Dataset S1). Transcriptomic data indicated that D. tunisiensis UC0050 differentially expressed genes involved in carbon metabolism and fixation pathways for energy production ( Table 2 and Fig. 4). It was observed a~4-fold increase in NADH:ubiquinone oxidoreductase expression (the NuoE subunit, forming the NADH dehydrogenase module), which may function as electron acceptor for the also consistently highly expressed flavodoxin FldA (Fig. 5 and Additional file 2: Dataset S2). The NADH:ubiquinone oxidoreductase enzyme is indeed a proton pump (also known as complex I), which couples electron transfer with the translocation of four protons through the membrane [51]. This electron flow can be mediated via a reduced flavodoxin, such as FldA, which acts as intermediate between central carbon metabolism (e.g., TCA cycle) and complex I [52]. Thus, the upregulation of these genes suggests an increased activity of the electron transfer chain via H 2 oxidation [51], and may be involved in syntrophic relationships with hydrogenotrophic species throughout the increased proton extrusion from the cell.
Similarly to D. tunisiensis UC0050, C. proteolyticus UC0011 also differentially expressed genes related to carbon fixation pathways ( Table 2 and Fig. 4). Specifically, C. proteolyticus UC0011 boosted the reductive tricarboxylic acid cycle (rTCA) by a~6-fold increase in the expression of pyruvate:ferredoxin oxidoreductase (PFOR) and phosphoenolpyruvate carboxykinase (Pck) (Fig. 5 and Additional file 2: Dataset S2). Such regulation indicates an uptake of the excess acetate for pyruvate production, from which other central metabolic intermediates can be formed. C. proteolyticus UC0011 also regulated genes involved in amino acids metabolism, including a~4-fold upregulation of the ATP-dependent protease Clp, along with enzymes metabolizing various amino acids (arginine, alanine, glutamate, tryptophan, aspartate) (Fig. 4, Additional file 2: Dataset S2). Since H 2 is one of the main products derived from proteins and amino acids degradation by C. proteolyticus [40,41], it cannot be excluded that this microbial species can form syntrophic association with hydrogen-scavenger microorganisms, such as the hydrogenotrophic methanogen M. wolfeii UC0008. Previous studies pointed out a synergistic effect operated by the co-existence of proteolytic anaerobes and hydrogen-consuming methanogens, revealing an augmented cell growth and protein degradation efficiency [53]. A partnership between C. proteolyticus and archaeal species belonging to the Methanothermobacter genus has also been recently proposed [54,55].
A. hydrogeniformans UC0046 responded to H 2 injection by upregulating a H + -ATPase (NtpB) (Fig. 5 and Additional file 2: Dataset S2) that extrudes protons through ATP hydrolysis, and by downregulating the expression of the Na + /proline symporter (PutP) and Na + /H + antiporters (MnhC, NhaC) (Additional file 2: Dataset S2). These mechanisms are used by various anaerobic bacteria to regulate internal pH and to control the transmembrane electrochemical gradient [56]; however, a syntrophic mechanism also cannot be excluded for this species. Additionally, A. hydrogeniformans UC0046 upregulated the coenzyme F420-reducing hydrogenase (FrhA), the electron transfer flavoprotein (ETF: FixA), and the NADH:ubiquinone oxidoreductase (NuoE), all known to be involved in mechanisms of electron flow and energy production [51,52] (Fig. 5 and Additional file 2: Dataset S2). As for D. tunisiensis UC0050, the upregulation of these genes suggests an involvement of A. hydrogeniformans UC0046 in syntrophic relationships with hydrogen-scavenger microorganisms.
The less abundant SAOB T. acetatoxydans UC0018 increased by almost 8-fold after H 2 injection (Fig. 3, Additional file 1: Figure S4, Additional file 2: Dataset S1). There was an upregulation of glucose metabolism (FruK, Fba), as well as sugar and branched-chain amino acid ABC transporters (Rbs, Mgl, and LivK) in this species (Fig.4, Fig. 5, Additional file 2: Dataset S2). T. acetatoxydans UC0018 also upregulated the rTCA key enzymes pyruvate:ferredoxin oxidoreductase and phosphoenolpyruvate carboxykinase (PorA and PckA) by 8-and 4-fold, respectively ( Fig. 5 and Additional file 2: Dataset S2), indicating an acetate uptake probably aimed to increase the energy store capacity. It is indeed known that the utilization of the TCA cycle in the reductive direction by many autotrophic anaerobes is aimed at producing metabolic intermediates via acetyl-CoA incorporation [57]. The rTCA upregulation seen in T. acetatoxydans UC0018 indicates the different metabolic direction taken by T. acetatoxydans, which did not act as a SAO by upregulating enzymes for acetate oxidation. A significative correlation between T. acetatoxydans UC0018 in phase II and acetate concentration in the reactor was also indicated by statistical analysis (Additional file 1: Figure S5).
The archaeal consortium was composed of three equally abundant methanogens: M. wolfeii UC0008, the generalist Methanosarcina thermophila UC0006 [58], and the methylotrophic Methanomassiliicoccus sp. UC0009 [59]. Only the latter decreased in abundance after H 2 injection (Fig. 3, Additional file 1: Figure S4, Additional file 2: Dataset S1), showing a~4-fold reduction and indicating that this species may be more sensitive to the new condition. The remarkable decrease of Methanomassiliicoccus sp. UC0009 and therefore the CH 4 produced by the methylotrophic pathway could be also one determinant of the lower methanation seen in R3 after H 2 addition. Additionally, although M. thermophila UC0006 and M. wolfeii UC0008 remained quantitatively stable and upregulated the aceticlastic (M. thermophila UC0006) and hydrogenotrophic pathways in phase II (Additional file 1: Figures S7 and S8), the drop in CH 4 content seen in R3 was unchanged (Fig. 2a and Table 1). However, it is worth to highlight that the presence of M. thermophila UC0006 in R3 may have allowed the lower accumulation of acetate compared to the single-stage configuration.

Simulation of hydrogen utilization routes in R1 and R3
The same bacterial species had different regulatory responses in the two reactors. This diverse regulation could be due to the reactor configurations, resulting in different physicochemical conditions, and consequently different H2 utilization capability.
To support the experimental findings based on gene expression results (Fig. 4 and Fig. 5), which showed the main metabolic pathways undertaken by the MAGs, a computational model of the two reactor configurations was also developed. Moreover, mass balance calculations contributed to clarify the processes occurring in the reactors after H 2 injection (Additional file 1). It was found that approximately 40% of the H 2 moles injected in the single-stage configuration (R1) per day were effectively utilized to produce CH 4 . Software simulation results for R1 showed trends similar to those obtained experimentally. In particular, both simulated methane production and total VFA concentration curves agreed with the measured values, although being slightly lower (Additional file 1: Figures S9  and S10). However, the remaining 60% of the added H 2 was not enough to account for the butyrate increase seen in the same reactor (Additional file 1). Therefore, the additional utilization of internal H 2 produced by lactose fermentation to acetate and butyrate should be considered. The most reasonable hypothesis in terms of demand for indigenous H 2 moles suggests that propionate reduction increased butyrate (Additional file 1). Concerning the acidogenic reactor of the two-stage configuration (R2), the fractions of exogenous H 2 moles utilized for the butyrate augmentation were~97%, 60%, and 30%, based on the substrate reduced (CO 2 , acetate, and propionate, respectively). Since the amount of propionate in R2 was negligible and butyrate increase cannot be primarily based on CO 2 reduction (considering the 30% residual hydrogen in the effluent gas from R3), the most probable explanation for butyrate production is via the acetate reduction. This was further confirmed by the slight decrease of acetate concentration in R2 (Fig. 2c). Additionally, the acetate rise seen in the methanogenic reactor of the same configuration (R3) during phase II mostly relied on butyrate oxidation (~60%), also considering the augmented butyrate feeding from R2. Overall, model simulations of the methane production and changes in total VFA concentration were in agreement with the above assumptions (Additional file 1: Figures S15 and S16). Finally, since onlỹ 10% of the acetate rise in R3 could be explained via the acetogenic pathway, it is reasonable that the remaining 30% of the acetate mole increase was likely due to an accumulation effect, which may have inhibited the acetogenic pathway.
Computer-aided simulations, combined with mass balance calculations, indicate the most probable H 2 availabilities in the two reactor configurations and therefore the different metabolic routes for H 2 utilization used by the anaerobic digestion microbiome.

Conclusions
H 2 injection induced different transcriptional regulation responses in the same MAGs inhabiting the two reactors.
Specifically, they favored methanation in the single-stage reactor (power-to-methane), and SCFAs production in the two-stage configuration (power-to-chemicals). The above finding was also confirmed by model simulations. Gene expression results revealed that C. proteolyticus UC0011 and A. hydrogeniformans UC0046 mainly upregulated pathways involved in acetate and butyrate production. However, a 7% increase in CH 4 content in the biogas of the single-stage reactor was observed, mainly due to the dominant hydrogenotrophic M. wolfeii UC0008. In contrast, a doubling of total SCFAs by CO 2 fixation was evidenced in the two-stage configuration, with A. hydrogeniformans UC0046 and D. tunisiensis UC0050 upregulating genes involved in electron transport chains. Interestingly, the SAOB T. acetatoxydans UC0018 did not act as acetate-oxidizer in either reactor configuration, but primarily inhibited sugar metabolism in the single stage and boosted acetate uptake via the reductive TCA cycle in the two-stage configuration. A putative syntrophism between C. proteolyticus UC0011 and M. wolfeii UC0008 was proposed in the serial reactor configuration.

Additional files
Additional file 1: 16S rRNA gene amplicon results (Table S2 and Figures S1-S2), differentially expressed genes not assigned to MAGs ( Figure S3), relative abundance of the 50 MAGs ( Figure S4