Impact of trace metal supplementation on anaerobic biological methanation under hydrogen and carbon dioxide starvation

Earth’s atmosphere reached a concerning carbon dioxide (CO₂) concentration of 417 parts per million in 2022¹, thus it is urgent to study new technologies for efficient CO₂ capture in waste gas streams. Within this framework, there has been an increased focus on carbon capture and utilization (CCU) technologies². Microbial-mediated bioconversions not only have the potential to mitigate greenhouse gas emissions, but they can generate energy carriers, such as methane (CH₄), and other valuable compounds³. Within nature, carbon fixation microbiomes are ubiquitous, spanning across diverse environments ranging from oceanic deep vents⁴, lakes⁵ and marine sediments⁶, peatlands⁷, to hydrocarbon reservoirs⁸. However, studying and characterizing these natural systems pose significant challenges. Therefore, a simpler and more accessible environment is required to model processes such as methanogenesis. Engineered systems, like bioreactors, offer a viable solution by creating what is often termed as an ‘ecosystem in a vessel’, where bioconversions can be effectively studied and utilized. Additionally, these controlled ecosystems offer the advantage of valorizing industrial wastes by leveraging the metabolic potential of microbial consortia^3,9. At present, the industrial-scale adoption of H₂-mediated biological CCU technologies remains limited², often contingent on primary materials and specialized growth media needed by pure cultures¹⁰. Therefore, it is imperative to investigate the physiology of microbiomes to optimize their biosynthetic capability and harness it proficiently. Additionally, green hydrogen (H₂) provision is not constant, as it relies on an energy surplus from solar or wind to fuel the water electrolysis. Thus, microbial communities’ ability to respond to transient starvation events must be further explored.

The gaseous carbon transformation processes rely on the delicate equilibrium involving the structure and function of several anaerobic microbial species, influenced by numerous biotic and abiotic variables. Specifically, the conversion of H₂ and CO₂ into CH₄ is mainly orchestrated by hydrogenotrophic archaea, including genera such as Methanothermobacter and Methanoculleus, which can use H₂ as an electron source¹¹. Conversely, the role of bacteria in CO₂ methanation has been recently proposed but is still not fully elucidated. For instance, selected species can support biofilm formation and syntrophic interactions to overcome auxotrophies were proposed^12,13. These relationships are centered around the canonical Wood–Ljungdahl pathway (WL)¹⁴ or its alternative version, involving the glycine synthase-reductase (GSR) route^15,16. This metabolism couples the methyl branch of the WL pathway with amino acid (AA) metabolism, thereby regulating carbon fluxes within microbial cells. However, metagenomics on its own falls short in providing a sufficiently detailed representation of the diversity and dynamics inherent in methanogenic processes. Hence, a deeper investigation performed through transcriptome-level analysis is needed for acquiring a comprehensive picture of the cellular activity and metabolic behavior¹⁷. Metatranscriptomics is key for the investigation of highly complex phenomena such as the impact of H₂/CO₂ deprivation and micronutrients on biomethanation, which have not been explored up to now. For example, previous demand-oriented biological methanation experiments have delineated the biosynthetic roles of archaeal and bacterial species¹⁸. Similarly, the response to different carbon regimes has been suggested to be critical for soil health, due to its impact on microbiome-mediated biological processes¹⁹. However, evidence from transcriptome expression under carbon and electron source starvation conditions is still lacking, and additional omics (i.e metabolomics, proteomics, etc.) are required to obtain a comprehensive understanding.

Microbial cells require various micronutrients for their growth and activity, with Iron (Fe), Nickel (Ni), and Cobalt (Co) specifically implicated in methanogenesis²⁰. Earlier studies found trace metals to positively impact digester stability long-term, yet explanations remain case-specific due to anaerobic digestion’s operational variability²¹. The micronutrient role is interconnected with their presence in numerous enzymes and cofactors essential to the process. For instance, Ni is an integral component of coenzyme F430²², essential for the activity of methyl coenzyme M reductase (Mcr)²³, while Co is a constituent of the two cobamide factors in tetrahydromethanopterin S-methyltransferase (Mtr)^24,25. In anaerobes, the precorrin-2 precursor, undergoes oxidation catalyzed by precorrin-2 dehydrogenase (CysG). Subsequently, Co²⁺ is inserted via sirohydrochlorin cobalt/nickel chelatase (CbiK) activity, and the resulting Cobalt(II)-sirohydrochlorin undergoes further processing to generate cobalamin. Alternatively, if sirohydrochlorin is directed towards the synthesis of the Ni-containing Coenzyme F430, the same CbiK catalyzes Ni²⁺ insertion, resulting in Nickel(II)-sirohydrochlorin²⁶. Different concentrations of micronutrients have been previously tested in pure cultures of methanogens to identify metal specific optimal values^27,28,29. Conversely, the effect of trace metals augmentation on mixed communities adapted to feeding through H₂ and CO₂ streams, has never been explored before. In such systems, syntrophies play a fundamental role, and the impact of metals must be viewed from a community perspective, where a domino or cascade effect may occur. Furthermore, not only the transcriptional activity but also the genetic makeup of the microbiota can be crucial. Recent strain-level investigations have highlighted how specific stressors can drive the accumulation of single nucleotide variants (SNVs), whose impact can positively affect enzyme function^30,31.

In this study, a combination of omic techniques were employed to elucidate the resilience of the microbiome in recovering after a starvation period during which no H₂ and CO₂ feeding was provided. This intermittent experimental design involved three stepwise periods of starvation lasting three, five, and ten days, respectively. The community was initially starved, the feeding regime was then restored and, once stable biomethanation was achieved, feeding provision was halted again for an extended time. To assess the effect of trace metals on post-starvation recovery, the growth medium was supplemented with either Ni or Co, both of which act as cofactors for enzymes catalyzing key steps in methanogenesis. These metals could potentially enhance the restoration of methanogenic activity following starvation. H₂ supplementation has been shown to significantly upregulate the methanogenic pathway in both mesophilic and thermophilic hydrogenotrophic methanogens within anaerobic digesters^32,33. This indicates that, following prolonged deprivation of H₂/CO₂, hydrogenotrophic methanogens likely enhance the biosynthesis of key enzymes, necessitating essential cofactors such as trace metals. To analyze microbial cellular responses, the current study employed a methodology that integrates metagenomic and transcriptomic data into genome-scale metabolic models (GEMs) to perform flux balance analysis (FBA). These computational models effectively elucidated the bioconversion capabilities, growth requirements, metabolite exchanges occurring, and susceptibility to inhibitors of the microbial activity^34,35,36. Moreover, a strain-resolved approach encompassing SNVs profiling of the species allowed to reveal fine scale evolutionary mechanisms. This work offers a comprehensive insight into the genetic heterogeneity of carbon fixating microbiomes and explores the effect of trace metals supplementation on community metabolism.

Temporal tracking of biochemical dynamics in CO₂ biomethanation

A time-course investigation was conducted to determine the impact of trace element supplementation (i.e. Ni and Co) on the structure, activity and resilience of the microbiome following substrate-deprivation stress. The experimental design consisted in one start-up period (ST), three progressively longer starvations (SP1, SP2 and SP3), each one followed by feeding restoration (OP1, OP2 and OP3). The process parameters were systematically monitored, with emphasis on microbial growth rate, outlet gas composition, daily CH₄ production, micronutrient uptake, and VFAs profiles (Supplementary Table 1–5, Fig. 1). Through the tracking of these indicators, the microbiome behavior in response to both trace metal supplementation and intermittent H₂/CO₂ feeding was thoroughly characterized at cellular level. The growth curves registered during the initial exponential phase overlap in all the conditions (Fig. 1b). Notably, after each starvation, the Ni cultures were capable of switching faster from lag to growth phase. A sustained growth (OD ~ 0.6) was reached faster in Ni cultures, when compared to CTRL and Co (Fig. 1b), especially in OP2 and OP3. Despite this kinetic difference, the same OD values were reached at the end of the last operation period (OP3), showing an overall similar profile for the microbiomes which, eventually, converged towards the same direction. Daily measures of gas volume and composition were used to estimate the biomethanation performance in terms of CH₄ yield. After the end of SP1 and SP2 the communities supplemented with Ni and Co showed higher CH₄ content, reaching a production 2.5 and 4.5 times higher than control (Fig. 1c). The third starvation (SP3) was more severe, and the microbiome required more days to recover a stable production of CH₄ (above 100 mg/L), which was reached after 5 days (T46). Nonetheless, the Ni-supplemented community proved again to be the one with the best recovery when compared to CTRL and Co conditions. The reactors supplemented with Co showed a more fluctuating behavior (Fig. 1c), but overall, their performance was superior to the CTRL condition across the entire experimental period.

a Study design. Shotgun metagenomic and metatranscriptomic sequencing was performed on the highlighted time point. The sequential experimental periods are highlighted as follows: startup (SP), H₂/CO₂ starvation (SP1, SP2, SP3), and operation (OP1, OP2, OP3). b Measures of OD, (c) average CH₄ production (ml/day) of the microbiome during intermittent operation. The communities investigated are highlighted by a different color: control (Ctrl – blue), Ni-enriched (Ni – orange), Co-enriched (Co – green) and blank (Blank – red). d Co and Ni concentration in mg/L throughout the experiment in the three sets of reactors. The dotted line refers to the average concentration of Ni and Co in the blank reactors. Starvation periods are highlighted by a shaded gray area along the x-axis. Days with the highest differences in methanogenic activity between conditions are marked with “*”.

Full size image

Finally, to validate the uptake and use of Ni and Co by the microorganisms, concentration of trace metals dissolved in the medium was monitored. An initial decrease was measured during the AP and prolonged up to the end of OP2 (day 31), with Co and Ni reaching the minimum concentration level of 0.40 mg/L and 0.215 mg/L (Fig. 1c), respectively. Afterwards, a progressive increase in metals was observed during OP3, with Co and Ni concentrations reaching 0.129 mg/L and 0.296 mg/L at day 45. VFA profiles were also analyzed throughout the experiment. In all conditions, a notable rise in short acid concentration was observed during the first week, predominantly marked by the biosynthesis and accumulation of acetate and propionate (Supplementary Table 5). Subsequently, under intermittent operation, acetate concentration declined until OP3, while propionate levels steadily increased across all reactors (~50 mg/L). Notably, in Ni supplemented culture, acetate stabilized within the range of 50–85 mg/L, while under Co and CTRL conditions, its concentrations were recorded at 7 and 10 mg/L by day 22 (Supplementary Table 5). The VFA accumulation observed occurs upon carbon deprivation, implying that the cellular metabolism switches from methanogenesis to the GCS pathway, with concurrent short acid production.

Microbiome profiling during intermittent feeding

Genome-centric metagenomics was used to explore the microbiome composition of the three growing conditions (CTRL, Ni, and Co) and track its evolution under an intermittent feeding provision regime. The resulting assembly of the shotgun reads covered a total of 5.45 Gb, with an average alignment rate of 95%; the community was nearly entirely represented, a result mainly determined by its low complexity, which allowed an accurate metabolic reconstruction. A low number of species was observed, due to the combined use of a simplification adaptation period, synthetic medium and simple carbon source for growth. The medium-to-high quality Metagenome Assembled Genomes (MAGs) recovered after metagenomic binning comprised 21 species (Fig. 2a, Supplementary Table 6), covering >90% of the overall relative abundance (RA). Notably, Methanothermobacter thermautotrophicus MN_MX1 was dominant across all the samples (average RA 87%), being the sole detected methanogen in both the inoculum and the experimental period. This species is recognized as a hydrogenotrophic archaeon capable of exploiting the bioconversion of CO₂ and/or CO with H₂ to CH₄ in thermophilic conditions³⁷. The elevated RA combined with the known preferences for the provided carbon source suggests M. thermautotrophicus MN_MX1 plays a major role within the community. However, despite the favorable experimental condition, the RA of this archaeon decreased during starvation periods, with values ranging between 80-75% upon feeding reestablishment (Supplementary Table 7). This highlights its sensitivity toward resource limitation in comparison to other components of the community. To delve deeper, the predominant fraction of the microbiome was investigated, focusing the attention on the four bacterial species with RAs exceeding 1% (Fig. 2b) in at least one time point. Caldanaerobacter subterraneus MN_VB17 stood out as the most abundant among the bacterial population. Isolation studies identify this species as a thermophilic fermentor capable of growth on carbohydrates, producing acetate, L-alanine, H₂, and CO₂³⁸. It is important to emphasize that, given the experimental setup, potential syntrophies are not based on CO₂ and H₂ flows, which are externally provided in large excess, but rather on the exchange of other compounds. Conversely, potential competition for CO₂ is likely negligible, considering the ample feeding quantity provided.

a Beta diversity of reconstructed MAGs using the robust Aitchison metric. The communities investigated are highlighted by a different color: control (Ctrl – blue), Ni-enriched (Ni – orange) and Co-enriched (Co – green). b Phylogenetic tree of the reconstructed MAGs; completeness and contamination are reported in the barplot to the right. c Average RA and RPKM values of MAGs across the sampled data points. The full name of the MAGs can be found in Supplementary Table 6. Each microbial community is labeled as follows: starting inoculum (Inc), control (Ctrl), Ni-enriched (Ni) and Co-enriched (Co).

Full size image

Globally, the microbiome exhibited fluctuating dynamics in RA during the CO₂ methanation process (Fig. 2b). Nonetheless, the total number of species remained constant, as indicated by the Chao1 index, while alpha diversity revealed a consistent positive trend at the later time points (Supplementary Table 8). Rather than an increase in species richness, these results imply a higher evenness in RA of species occurring during the later stages (i.e., C. subterraneus MN_VB17, Bacillota sp. MN_19). A noteworthy distinction arises in the Co condition, where the overall diversity value was markedly lower compared to those found in CTRL and Ni (Supplementary Table 8). This result is associated with a lower increase in RA of the rare species observed in the Co reactors. Analysis of beta diversity revealed distinct separation among the three microbial communities (Fig. 2a). The Ni-amended microbiome formed a separate cluster, while the CTRL and Co conditions exhibited greater similarity. Moreover, time points corresponding to recovery periods (T21, T31, T48) clustered together, distant from post-starvation samples (T12, T27), indicating a shift in community structure induced by imposed C-deprivation.

Assessing species metabolic potential and putative interspecies interaction through gene-level survey

Metabolic functional reconstruction, aimed at profiling the microbiome potential, targeted specific pathways of interest, including hydrogenotrophic methanogenesis, carbon utilization routes, and micronutrient cellular uptake systems. The completeness of essential pathways in microbial metabolism (KEGG modules) was evaluated in the subset of five selected MAGs. M. thermautotrophicus MN_MX1 was categorized as a hydrogenotrophic methanogen (M00617) in accordance with its taxonomy, while all other bacteria were identified as homoacetogens (M00618). As expected, functional analysis of the archaeon revealed a complete hydrogenotrophic pathway (M00567), supporting autotrophic CO₂ assimilation for carbon fixation (Fig. 3a). The pathways for the biosynthesis of coenzymes involved in methanogenesis were present, encompassing methanofuran (M00935), coenzyme M (M00358), coenzyme F420 (M00378), and coenzyme F430 (M00836). Additionally, M. thermautotrophicus MN_MX1 possesses the complete pathway for cobalamin biosynthesis via sirohydrochlorin (M00924), from which coenzyme F430 is derived.

a Metabolic reconstruction of the main pathways in the bacterial population (left, pink cell) and the methanogenic archaeon (right, blue cell). The degree of completeness of C-metabolism KEGG modules is depicted in the heatmap on the left, wherein fully annotated modules are highlighted in orange and incomplete in gray. b Transcriptional profiles of key genes were analyzed before and after the first starvation period (9–12 days). Changes in gene expression levels between conditions (CTRL, Ni, Co) are represented as log₂FC, calculated from three replicates. The log₂FC (left panel) and average FPKM (right panel) of annotated genes involved in selected pathways are reported: fatty acid biosynthesis, glycolysis, pyruvate oxidation, citrate cycle, propionate oxidation and hydrogenotrophic methanogenesis. Additionally, genes encoding for RNA polymerase and transcription factor (Bacteria/Archaea Transcription), as well as ribosomal proteins (Bacteria/Archaea Translation), were considered as basal cell activity indicators.

Full size image

The investigation performed on the bacterial microbiome focused on the capacity to metabolize VFAs and other short-chain molecules, to highlight potential metabolite exchanges with the archaeon and thereby infer putative syntrophic interaction. Notably, the module associated with the canonical WL pathway was incomplete across all MAGs (M00377). Nevertheless, the GCS pathway was prevalent among bacteria, allowing for the possibility of syntrophic acetate oxidation (M00621) and was reconstructed by manually checking the presence of specific genes (Fig. 3a). Specifically, the formate-tetrahydrofolate ligase (fhs), a molecular marker for syntrophic acetate oxidizing metabolism, was identified among the annotated genes of C. subterraneus MN_VB17. This species lacks formate dehydrogenase, but it has a formate transporter (fdhC), suggesting potential metabolite exchange with M. thermautotrophicus MN_MX1. Following uptake by the archaeon, formate undergoes conversion into CO₂ through formate dehydrogenase (NAD⁺), serving as a substrate for hydrogenotrophic methanogenesis. Other bacteria also share the GCS route (Fig. 3b), along with one or more formate dehydrogenases (Fdh, HylABC). The only exception was identified in Bacilli MN_MA21, which converts acetate into acetyl-CoA via acetyl-CoA synthetase (ACSS).

Exploration of transporters for metals, crucial for their cellular uptake, identified all genes encoding subunits of the cobalt/nickel transport system (cbiO, cbiM, cbiN, cbiQ) as present in M. thermautotrophicus MN_MX1, as well as C. subterraneus MN_VB17 and Bacillota MN_MA19. Other transport-related genes were included in the search because, although substrate-specific, their selectivity does not preclude affinity for other compounds. Specifically, transporters for iron (ABC.FEV.S, ABC.FEV.P, ABC.FEV.A), tungstate (TupABC), molybdate (MobABCF), phosphate (PstABCS), zinc (ZnuABC), and biotin (BioY, EcfT, EcfA1, EcfA2) were also present at the genomic level.

Transcriptome activity highlighted the potential syntrophic lifestyle of selected microbes

Investigation at transcriptomic level provides a mechanistic view of the metabolic process and helps to assess the impact of intermittent feeding on methanogenesis at molecular level. Additionally, transcriptomics enabled a verification of the metals’ effect on the post-starvation recovery of specific taxa. The nutritional stress induced a decrease in OD and lower growth rates as predicted by the peak-to-trough ratio (PTR) of most MAGs (Supplementary Table 9). Across all conditions, ten species exhibited a negative log₂PTR at the end of SP1, whereas three showed a slightly positive value. This latter result contrasts with the metabolic rewiring observed, which appeared to favor Tissierellaceae sp. MN_MX6 and Bacillota MN_MA19. These species demonstrated increased expression of ribosomal proteins, possibly reflecting a higher rate of protein biosynthesis (Fig. 3b). Nevertheless, previous studies have shown that when Escherichia coli experiences AA deprivation, it upregulates AA biosynthesis pathways to adapt to the nutrient downshift³⁹. Therefore, under stress, bacteria can activate specific stress-response mechanisms without necessarily shutting down transcription and translation machinery⁴⁰. Conversely, the translational and transcriptional activity of M. thermautotrophicus MN_MX1 was lower after SP1, indicative of decreased metabolic rate resulting from the starvation period. This stall was evident from the decline observed in the CH₄ production curve (Fig. 1b) and corroborated by an overall downregulation of genes associated with the hydrogenotrophic metabolism (Fig. 3b). Specifically, the expression levels of genes encoding methyl coenzyme M reductase (mcrA, mcrB, mcrG) were either downregulated in CTRL (p-value < 0.05) and Co reactors or exhibited no transcriptional changes in the Ni condition, with average log₂FC values of -1.05, -1.24, and 0.34, respectively. Notably, the activity level of substrate specific transported was altered, with the cobalt/nickel transport system (CbiMNOQ) being moderately upregulated in all conditions in the archaeon and downregulated in the bacteria (Supplementary Table 10). Additionally, the expression of zntA, encoding a transmembrane protein that facilitates transport of Zn²⁺ and Co²⁺, was upregulated in multiple bacteria upon starvation, but statistically significant only in Tissierellaceae sp. MN_MX6 (p-value = 0.0005). Conversely, the permease and substrate-binding subunits of the ABC-type cobalt importer (CbtJKL), were showing more than one-fold change increase in expression in Tissierellaceae sp. MN_MX6 and C. subterraneus MN_VB17 (Supplementary Table 10).

Three main metabolic routes were chosen to evaluate the species activity at the end of the experimental period: (a) the hydrogenotrophic methanogenesis (b) the cobalamin biosynthesis, and (c) Ni transport. M. thermautotrophicus MN_MX1 exhibited the highest activity in the Ni reactors, where the mcr gene displayed an FPK value 1.5 times higher than in the CTRL and Co reactors (Supplementary Table 11). The expression of mcr demonstrated a significant correlation with both CH₄ production rate and yield, indicating that monitoring this gene accurately reflects the performance of the biomethanation system. Remarkably, genes encoding for the NAD⁺ formate dehydrogenase (fdh) were exclusively active in the CTRL and Co-supplemented conditions. This finding is supporting the idea that Ni-supplied cultures were functioning optimally with the provided carbon source, while in other conditions, the archaeon had to rely on additional carbon compounds besides CO₂ for growth. The impact of Co on the transcriptional activity was assessed by exploring the expression of genes encoding for Cobalamin- and Corrinoid-dependent enzymes. The expression levels of cobalamin biosynthesis genes, exclusively observed in M. thermautotrophicus MN_MX1, were higher in the CTRL condition compared to Co, with an average log₂FC of +0.85 FPK (Supplementary Table 11). Conversely, in Tissierellaceae sp. MN_MX6, enzymes catalyzing carbon skeleton rearrangements dependent on cobamide or cobalamin exhibited increased activity in presence of Co. Specifically, genes encoding two aminomutases⁴¹, D-ornithine 4,5-aminomutase (oraE) and beta-lysine 5,6-aminomutase (kamD), were expressed 10 and 65 times more in the Co condition compared to CRTL, respectively. In the same way, the activity of Ni-related genes, including transporters and hydrogenases using Ni as cofactor, was investigated (Supplementary Table 12). The expression of peptide/nickel transport system (ABC.PE.P, ABC.PE.P1, ABC.PE.S, ddpE, ddpF) genes was increased in the Ni-supplemented condition for Tissierellaceae sp. MN_MX6 and C. subterraneus MN_VB17. The two MAGs displayed an average FC increase of +2.86 and +0.75 FPK, respectively. Notably, the gene nikA, responsible for the Ni²⁺ binding component of the ABC-type Ni²⁺ transporter⁴², reached nearly 15,000 FPK compared to the 1000 FPK observed in the CTRL condition (Supplementary Table 12).

Reconstruction of microbial interaction through genome-scale metabolic modeling

Insights into metabolic interactions driving the CO₂ methanation were obtained by reconstructing community-level GEMs throughout time and experimental conditions. The reconstruction resulted in a comprehensive community-specific model comprising on average 14 species, 21,378 reactions and 18,754 metabolites. Subsequently, FBA simulations maximizing the community growth were used to discern the exchanges occurring among taxa and with the growth medium. The space of solutions of the simulation were constrained using selected biochemical measures (i.e. consumption and production rates of VFA, H₂, CO₂, and micronutrients) and gene expression data to generate condition-specific simulations (Supplementary Table 13). GEMs performances were validated assessing the accuracy of CH₄ prediction, the key product of CO₂-biomethanation. Results revealed high positive correlation between the measured and predicted metabolic fluxes (Pearson correlation, r = 0.74, p-value = 4.96e⁻⁴), thus demonstrating the reliability of GEMs in reproducing the fluctuating dynamics characterizing the starvation and growth phases (Fig. 4a). Additionally, differences between the simulations of stress starvation (i.e. T12 and T27) and those under continuous feeding emerged when considering the overall spectrum of exchanged metabolites with the medium (Fig. 4b). As expected, the primary differences stem from the consumption of gasses introduced into the system, although a positive uptake of ammonia (NH₃) mainly related to protein biosynthesis was also observed during the operational periods, particularly in Ni-supplemented samples at the last time point. Indeed, microbiomes supplemented with Ni exhibited the highest CH₄ production predictions, with a positive efflux of ~12 mmol*DW⁻¹*h⁻¹ at time point T9 and T31, suggesting a greater methanogenic activity. At a species-level resolution, the central product of each microbiome is represented by CH₄, only produced by M. thermautotrophicus MN_MX1 (Supplementary Fig. 2). This species plays a pivotal role in the microbiome not only for CH₄ production, but also because it is capable of incorporating NH₄⁺ into AAs essential for microbial cell proliferation (Supplementary Fig. 3). Analysis performed at T21 evidenced that M. thermautotrophicus MN_MX1 can use NH₃ to synthesize AAs with high biosynthetic metabolic costs (eg. Arg and Thr)⁴³. NH₄⁺ is imported through energy-dependent Amt transporters and enters the high-affinity glutamine synthetase-glutamate synthase (GS-GOGAT) pathway. Part of these overall products are released by MN_MX1 and uptaken by the dominant bacteria (Fig. 4c, d, Supplementary Fig. 4), such as C. subterraneus in all conditions and Bacilli MN_MA21 in control and Co-supplemented cultures. Specifically, they appear to be involved in energy-production processes associated with the degradation of energetically rich AAs, releasing inorganic NH₃ and AAs with lower biosynthetic cost (e.g. Gly and Ala). The AAs oxidation is a basic bioenergetic pathway of the bacteria kingdom⁴⁴.

a Correlation between measured and predicted CH₄ production. Pearson correlation coefficient and its statistical significance are reported. The communities investigated are highlighted by a different color: control (Ctrl – blue), Ni-enriched (Ni – orange) and Co-enriched (Co – green). Each experimental time point was represented by a specific symbol: circle (T9), oblique cross (T12), square (T21), vertical cross (T27), diamond (T31) and star (T48). b PCA of the overall predicted metabolic fluxes per condition and across all the time points. Arrows illustrate the contribution of each variable, in this case the most relevant exchanges of the community with the environment, to the principal component. Illustration of the intricate network of interactions of the control (c) and Ni-supplemented (d) microbiome at time point T21. Focus is given to variations in exchanges of AAs driving the inter-species interactions in the two conditions. The size of arrows and dots are proportional to the predicted exchange and the species average RA, respectively. To each compound exchanged a specific color was assigned: ammonia (NH₃ – dark brown), alanine (Ala – light blue), arginine (Arg – magenta), aspartate (Asp – dark orange), glycine (Gly – dark green), glutamate (Glu – purple), serine (Ser – slate gray), threonine (Thr – yellow) and tyrosine (Try – gray).

Full size image

Variant analysis revealed stress-induced genomic heterogeneity of microbial species

Merely assessing species abundance and activity is insufficient for explaining the genetic dynamics within microbiota. While RA may remain constant, it does not ensure genetic stability, as significant genomic alterations can still occur. Here, variants were investigated at the species level, focusing on SNVs within dominant MAGs. A total of 8452 unique SNVs were identified across seven samples, categorized as synonymous (23.5%), nonsynonymous (67.1%), and intergenic (9.4%). The time-course trajectory of microbial strains associated with dominant species was scrutinized by analyzing the shift in SNV frequency across experimental stages. Results unveiled the accumulation of consistent genetic differences in the MAGs, potentially impacting the phenotype. For example, Tissierellaceae sp. MN_MX6 harbored a higher number of SNVs in the Ni condition; on the other hand, M. thermautotrophicus MN_MX1 genome was more conserved, with a number of variants ranging between 10 and 20, depending on the time point and condition. Nonsynonimous SNVs (nsSNVs) were further classified based on their specificity for control, Ni- and Co-enriched conditions, totaling 2184, 3487, and 2971, respectively. After the removal of variants shared among different conditions (Methods 2.9), 2592 and 2118 unique nsSNVs were specifically linked to the Ni and Co conditions. This implies that, in terms of the overall number of variants, ~25–30% of them are not favored by the selective parameters imposed. Theoretically, a SNVs can be fixed over time when it provides a phenotypic advantage to the microbe, thus, to identify these variants, only SNVs with a frequency exceeding 0.5 in the last two time points were considered.

Additionally, the Grantham distance (Gd) metrics was used as a proxy for the biochemical dissimilarity between two AAs based on three key side chain chemical properties: composition, polarity, and molecular volume⁴⁵. This distance highlights the impact of missense substitutions, with a higher distance probably determining a higher phenotypic effect. According to these restrictions, only 26 and 30 nsSNVs were identified exclusively in Ni and Co conditions, respectively. In the Ni-enriched cultures, Tissierellaceae sp. MN_MX6 had a variant in the marker gene encoding for syntrophic acetate oxidizing metabolism, the formate-tetrahydrofolate ligase (fhs). This SNV changed the AA from asparagine to serine (N221S), a modification with a relatively low impact (Gd: 46) (Fig. 5). Conversely, under Co-enriched conditions, C. subterraneus MN_VB17 accumulated SNVs on key genes of the GCS pathway, specifically encoding the phosphate acetyltransferase (pta) and the component E1 of the pyruvate dehydrogenase complex (pdhA). The two SNVs involved a change from valine to proline and to isoleucine respectively (V54F and V396I), both with a mild potential impact (Gd < 50). A third SNVs was identified in Bacillota MN_MA19 (Q355K) on the gene encoding beta-lysine 5,6-aminomutase (kamD), an enzyme with pyridoxal 5’-phosphate and adenosylcobalamin coenzymes (Fig. 5). The observed shift from glutamine to lysine (Gd: 53), suggests a potentially moderate effect on the protein functionality. Additionally, during OP1, genes responsible for high-affinity transport systems exhibited an accumulation of nsSNVs. However, it is noteworthy that these variants reverted to the baseline frequency levels in subsequent experimental stages (OP2 and OP3). The ABC-type Ni²⁺ transporter (nikA) and the permease of the peptide/nickel transport system (ABC.PE.P, ddpF) evidenced the same peculiar pattern in Tissierellaceae sp. MN_MX6.

Both the original peptide (light blue) and the one with the variants (light brown) were modeled using AlphaFold and aligned. The AA changes are highlighted in orange and reported by concatenating the “old AA”, the “position” and the “new AA”. The predicted structures were superimposed and visualized with ChimeraX. The surface of the modified protein is reported in red for the modified AA and in green for the original one. The 3D structure was rotated from left to right to provide three side-views. a N221S in fhs of Tissierellaceae sp. MN_MX6. b R131L in cobC of C. subterraneus MN_VB17. c Q355K in kamD of Bacillota MN_MA19.

Full size image

To ascertain whether the primary evolutionary force influencing the microbiome was the intermittent feeding regimen rather than the supplementation of additional micronutrients, an investigation was conducted into the SNVs shared among conditions. A total of 765 nsSNVs were identified, of which 127 had a frequency equal or higher than 0.5 in the last two time points. Genes under selective pressure involved the R131L replacement identified on the alpha-ribazole phosphatase (cobC) of C. subterraneus MN_VB17 (Fig. 5). This enzyme is responsible for the catalysis of the last step in the anaerobic pathway of adenosylcobalamin biosynthesis⁴⁶, as well as participating in a salvage pathway that recycles cobinamide into adenosylcobalamin⁴⁷. The change from arginine to leucine entails moving from a charged side-chain AA to one with non-polar sidechains (Gd: 102). Additionally, in Bacillota MN_MA19 a fixed SNV was mapped on por, whose encoded pyruvate:ferredoxin oxidoreductase is involved in catalyzing the decarboxylation of pyruvate during the alternative WL pathway. More into detail, the variant results in a switch from cysteine to tyrosine with differences in the R groups, the first being hydrophilic and the second relatively hydrophobic due to its aromatic side chain (Gd: 194). The spatial localization of the mutated AA was also confirmed through alignment of the protein sequence against the Swiss-Prot database, juxtaposing the MAG-predicted peptide with annotations from other organisms (Supplementary Figs. 5–8).

Ni and Co are important metal ions for methanogenesis often supplemented at industrial scale. However, their effect on the process performance and on CH₄ productivity at metabolic level is still unclear. More specifically, it is important to determine under which circumstances the addition of these compounds can be beneficial to the carbon bioconversion. Our experimental setup highlighted that Ni was favoring growth after a starvation period, especially when the lack of carbon compounds was extended up to 10 days. Ni-enriched conditions reached consistently higher microbial growth, suggesting that this is the most valuable treatment for increased cell productivity and stress resistance. Indeed, the positive effect of Ni on methanogenesis has already been demonstrated in anaerobic digestion systems, including mono-digestion of maize silage⁴⁸ and wheat stillage⁴⁹, as well as in pure cultures of Methanosarcina barkeri^20,50,51. However, this study is the first instance of successful metal amendment for CO₂ biomethanation implementing insight at molecular level on a mixed hydrogenotrophic-centric microbiome. Conversely, supply of Co did not affect the microbial growth rate but provided only a marginal benefit to the methanogenic activity. This dissociation between microbiome growth and moderately increased CH₄ production suggests a lack of direct correlation among these two metabolisms. In fact, methanogenesis, associated with energy generation, is influenced by diverse physiological processes, especially at high H₂ concentrations⁵². Microbial cells, beyond biomass generation, allocate energy to repair mechanisms, homeostasis maintenance, or compound synthesis under stressful conditions⁵³.

The simplified microcosm, dominated by the hydrogenotrophic M. thermautotrophicus, did not experience any major shift following the absence of feeding provision. The system demonstrated solid CH₄ production (>80% purity) throughout the whole operational period, comparable with similar previous setups for CO₂ biomethanation⁵⁴. Notably, following each starvation period, the microbiome demonstrated the ability to restore its activity within 2-3 days, underscoring resilience to sequential feeding interruptions of increasing duration. Recovery from starvation within 1 week of operation has been previously demonstrated in fed-batch reactors supplemented with acetate or glucose and H₂⁵⁵. Specifically, the dominant methanogen, Methanoculleus sp., was replaced by Methanobacterium sp., in diverse biomethanation systems⁵⁵. In a continuous stirred tank reactor testing multiple “on-off” H₂-feeding regimes, a similar shift in the microbiome population occurred, evidencing an increased RA in the acetoclastic archaeon (i.e., Methanothrix sp.), showcasing high adaptive potential¹⁸. In contrast to previous literature, our study revealed a persistence of the same dominant archaeon, M. thermautotrophicus, despite the imposed starvation periods. The results underline the fast growth kinetics of this microbe, as well as its proficiency in recovering from feeding stresses when compared to potential competing archaeal species. However, caution should be exercised when attributing superior competitive advantages to M. thermautotrophicus, since biochemical and operational conditions (e.g., experimental setup, temperature), as well as the complexity of the tested inoculum, vary considerably across studies. To gather more information regarding the microbiome resilience and decipher the role of Ni and Co in the recovery process, FBA simulations were applied and disentangled the dynamics governing the species interactions. Notably, in the exchanges between M. thermautotrophicus MN_MX1 and the bacteria, AAs emerged as prominent components. Specifically, serine and threonine are predicted to be released by the archaeon, whereas glycine, alanine and NH₃ produced by bacterial metabolism are uptaken by the methanogen. The key role of AA in the community may be attributed to the maximum H₂ concentration tolerated by each degradation pathway⁵⁶. Previous studies have demonstrated that serine and threonine possess highly tolerant degradation pathways, making them effective energetic substrates in selective environmental conditions. Conversely, glycine and alanine degradation exhibit low tolerance to H₂ (<10 Pa H₂), thereby serving as effective building blocks for the archaeon. With respect to C. subterraneus, the most abundant bacterium, the predicted interaction is supported by the extensive repertoire of hydrolase genes, including glycosidases, esterases and proteases, encoded in its genome³⁸. This enzymatic versatility offers a broad spectrum of carbon and energy substrates to the species. Moreover, alanine has been previously recognized as a byproduct of C. subterraneus metabolism during carbohydrates utilization³⁸. Based on the predicted interspecies dynamics, these features have potentially contributed to the fitness of C. subterraneus MN_VB17 and led to the hypothesis of its syntropic role with the archaeon. While GEMs offer a detailed platform for investigating species’ metabolic potential, they possess inherent limitations. Defining gene-protein-reaction rules, particularly for transporters, remains a challenge and can significantly influence the solution space. Moreover, GEMs primarily describe trophic interactions through molecule uptake, conversion, and secretion, overlooking other crucial macroecological processes. For instance, interspecies interactions driven by senescence, necromass formation, and biomass recycling through heterotrophic growth are not captured by GEMs. These macroecological dynamics, in conjunction with molecular exchanges, can play a crucial role in shaping microbial consortia equilibrium. For instance, in co-cultures of Clostridium phytofermentans and E. coli, necromass recycling has been shown to substantially enhance community biomass productivity⁵⁷. Similarly, in thermophilic iron-oxide mats, viral predation of primary autotrophs, alongside metabolite exchange, supports heterotrophic growth⁵⁸. Such combined processes are likely influencing the microbial community in the present study as well. Here, M. thermautotrophicus MN_MX1, the main consumer of CO₂ and H₂, serves as the primary resource population. Continuous cell death due to aging likely drives material and energy flow to potential secondary-resource specialists, such as C. subterraneus MN_VB17 and Acetivibrionales sp. MN_VB16⁵⁹, thereby increasing the complexity of the trophic network.

The metabolic advantage of partner species was further explored diving into their genetic makeup. Indeed, the tracking of variants over time was already successful in identifying the genomic contribution to species selection actively involved in antibiotic response⁶⁰, dental implant diseases⁶¹, CO₂ fixation process³⁰ and resistance to high NH₃ concentrations³¹. Here, this strategy highlighted the genetic stability of M. thermautotrophicus, even during a prolonged starvation period. The small number of nsSNVs mapped on the archaeon and their coverage depth in different conditions, suggest the presence of only one dominant strain with no genetic susceptibility to either metals supplementation or intermittent feeding. Conversely, the bacterial population showed higher genetic heterogeneity, likely due to the competition arising for the gaseous substrates during the initial growth phase and the starvation periods. The selective pressure on microbial populations may have induced variants providing a selective advantage when changes in environmental conditions are occurring⁶². This diversity is characterized by the presence of several nsSNVs on genes encoding for enzymes involved in key carbon metabolisms. The mutated AAs localization in the protein 3D structure was validated through alignment with reference sequences from similar species, suggesting potential impacts on catalytic activity of the respective SNVs. Under Ni supplementation, the mutation in Tissierellaceae sp. MN_MX6 falls within the formyltetrahydrofolate synthetase domain of Fhs (AAs 19-558)⁶³. Additionally, in Bacillota MN_MA19, the encoded AA change is situated near the binding domains of Por for manganese and thiamine diphosphate (AAs 963 and 962-965), respectively⁶⁴. Altogether, these observations support the hypothesis that Tissierellaceae sp. MN_MX6 and Bacillota MN_MA19 exhibit enhanced carbon utilization capability in response to Ni-amendment.

Transcriptomic profiling allowed the identification of key microbial players, including the putative syntroph Tissierellaceae sp. MN_MX6, rendering an exhaustive picture of the microbiome metabolic activity, as already demonstrated in previous studies^17,65,66,67. Methanobacteriaceae and Tissierellaceae families were related with biogas production in anaerobic digestion systems⁶⁸, thus their co-occurrence in our setup is fostering the possibility of a syntrophic relationship. In parallel, the investigation of the anaerobic cobalamin biosynthetic pathway, exclusively present in M. thermautotrophicus MN_MX1, exhibited lower expression in the Co-supplemented condition compared to the control. This is potentially attributable to the increased availability of the essential trace-metal for the biosynthetic route of the cofactors. Conversely, in Ni-supplemented cultures the higher expression of high-affinity transporters, including the permease ABC.PE and the ABC-type NikABCDE in Tissierellaceae sp. MN_MX6 and C. subterraneus MN_VB17, suggest a physiological response to the elevated metal levels. Finally, through transcriptome analysis the metabolic determinant behind the beneficial effect of Ni on methanogenesis was elucidated. Specifically, M. thermautotrophicus MN_MX1 exhibited a higher expression in Ni-supplied cultures for genes encoding enzymes requiring this metal as a cofactor, including coenzyme F420 hydrogenase (Frh) and H₂:CoB-CoM heterodisulfide ferredoxin reductase. The latter is a complex formed by the heterodisulfide reductase (HrdABC) and the F420-non-reducing hydrogenase (MvhAGD)^69,70 coupling ferredoxin and heterodisulfide reduction via electron bifurcation in hydrogenotrophic archaea⁷¹. The Ni-supplement have significantly enhanced the CH₄ yield potentially by favorably modulating the expression of the [NiFe] hydrogenase MvhAGD, previously identified in complex with HdrABC under Ni-limiting conditions⁶⁹. Thereby the recycling of intermediates between the initial and final steps of methanogenesis will be facilitated.

Research on trace metal supplementation for enhancing biosynthetic potential remains mostly limited to anaerobic species growing in pure culture. In contrast, the dynamics within multi-species cultures are more intricate to monitor and unravel. The current study was designed to involve not only Ni and Co supplementation but also incorporate stepwise increments in starvation periods. Through metabolic modeling and ad hoc SNVs investigation the molecular determinants behind the fluctuations of biochemical indicators and transcriptome analysis were identified. Specifically, the integration of multi-omics analysis and FBA revealed that the CO₂ capture process was almost entirely sustained by M. thermautotrophicus across all tested conditions. Methanogenesis emerged as a central energy-providing metabolism, potentially supporting the growth of not only the archaeon but the entire microbiome. Putative interspecies exchanges suggest that energy, in the form of thermodynamically degradable AAs, primarily flowed from the archaeal to the bacterial populations. While metals amendments did not significantly alter the aforementioned dynamics, Ni supplementation favored the metabolic processes restoration. Increased Ni availability resulted in the upregulation of methanogenesis genes, specifically those encoding for Frh, MvhAGD, and Mcr, facilitating recovery from starvation phases by M. thermautotrophicus MN_MX1. Under this condition, the archaeon exhibited higher biosynthetic activity as predicted by the increase in NH₃/NH₄⁺ uptake. This response enhanced carbon and nitrogen recycling within the microbiome’s trophic network. Despite the intermittent feeding regime, the composition of the microbiome remained stable under Ni and Co treatment. Our model indicates that the most abundant bacteria play a crucial role in mitigating the metabolic deficiencies of the archaeon, which had not demonstrated yet the ability to grow on basic anaerobic media as an isolate. These findings suggest a complex syntrophic interaction between bacteria, such as C. subterraneus MN_VB17, and M. thermautotrophicus MN_MX1. Overall, this study provided a unique opportunity to deepen our understanding of how the microbial population diversity in CO₂ biomethanation systems, particularly under metals amendments, responds to subsequent nutritional stress.

Inoculum

The microbial community utilized in this study underwent a simplification process conducted in the laboratory. The microbiota derived from a thermophilic lab-scale trickle-bed reactor using biogas and digested municipal biowaste as nutrient sources⁷² was obtained from the Technical University of Denmark. Once stabilized, the inoculum was acclimated to a feeding regimen of H₂ and CO₂ in a 4:1 proportion within fed-batch reactors. The simplification process is evident from the observed reduction in alpha diversity, as indicated by a decrease in the Simpson index from 0.22 at generation 0 of the adaptation period³⁰ to 0.029 at the initial time point of the present study (Inc_T0). This adaptation transpired over 30 generations, involving incubation at 55 °C and continuous agitation at 130 rpm, with the provision of basal medium (BA)⁷³, yeast extract (0.2 g/L final concentration), and Wolin’s vitamin solution⁷⁴. Detailed composition of BA and vitamins solutions is described in Supplementary Table 1. After each 21-day generation, cultures were reinoculated in BA medium, using 10% of the final volume as inocula. This iterative process was repeated and maintained until the beginning of the experiment.

Experimental setup and design

The experimental set-up comprised three different conditions conducted in triplicate using 1 L bottles with an initial working volume of 270 ml. Each bottle was inoculated with 27 ml of the fed-batch liquid of a simplified community, containing 5.4 mg of yeast extract, 2.7 ml of Wolin’s vitamin solution at 1X, 2.7 ml of Na₂S·7-9H₂O (25 g/L), and 234.9 ml of BA medium within an anaerobic chamber. Prior to inoculation, all solutions, except for the inoculum, were flushed with N₂ for 10 min to ensure anaerobiosis. For each condition, a bottle without inoculum was included to validate the reliability of growth rate and micronutrient depletion measurements. The fed-batch reactors were incubated at 55 °C on an orbital shaker (NBIOTEK NB-T205L, Korea) with a constant agitation of 130 rpm. Daily, the pressure of each bottle was released, and cultures were replenished with 600 ml of H₂ and 150 ml of CO₂. To investigate the impact of H₂/CO₂ starvation on the microbial community, feeding was interrupted for distinct periods: 3 days (from day 9 to 11), 5 days (from day 22–26), and 10 days (from day 32–41).

Micronutrients supplementation

To delineate distinct experimental conditions during starvation, Ni and Co were introduced as trace metals, given their pivotal roles in various enzymes and cofactors for methanogenesis. Micronutrient solutions were prepared using Nickel (II) Chloride hexahydrate (NiCl₂ · 6H₂O) (Sigma-Aldrich, MO, USA) and Cobalt (II) chloride hexahydrate (CoCl₂ · 6H₂O) (Sigma-Aldrich, MO, USA) salts. The supplementation level was determined based on concentrations reported in previous studies, with maximum concentrations of 10 mg/L for Co and 5 mg/L for Ni, as well as their concentrations in the BA medium (0.023 mg/L for Ni and 0.012 mg/L for Co, respectively)⁷⁵. The aim of the supplementation to the cultures was to attain a final concentration of 1 mg/L for both trace metals.

Analytic measures

The gas composition in the batches was assessed on a daily basis. Two milliliters of gas were collected from the headspace of fed-batch bottles before feeding and analyzed using a gas chromatograph (8860 GC, Agilent Technologies, CA, USA) equipped with a thermal conductivity detector (TCD). For volatile fatty acid (VFA) analysis, one milliliter of the liquid sample was acidified with 40 µl of orthophosphoric acid (85%, Sigma-Aldrich, MO, USA), then centrifuged at 13,000 rpm for 12 min. One milliliter of the supernatant was transferred to a vial, and 100 µl of isocaproic acid (Sigma-Aldrich, MO, USA) was added as an internal standard for a final concentration of 100 mg/L. VFA concentrations were determined using an 8860 GC (Agilent) equipped with a flame ionization detector (FID) and a DB-FFAP fused silica capillary column (30 m, 0.25 mm ID, film thickness 0.25 µm), with helium as the carrier gas. The growth rate was estimated by measuring the optical density (OD) of each fed-batch in triplicate using 125 µl of the liquid phase per replica with a Tecan SPARK microplate reader at a wavelength of 600 nm. For the detection of trace metals dissolved in the microbial medium, measurements were conducted via ICP by MATANlab (BTS-Biogas s.r.l., Affi, VR, Italy). After thawing the samples, approximately 27 g of each were weighed in Falcon tubes, and the water was evaporated overnight in an oven at 105 °C. Subsequently, salts were resuspended in 8 ml of a mixture of hydrochloric acid and nitric acid in a 3:1 ratio, and the volume was adjusted to a final volume of 20 ml with ultrapure water.

DNA/RNA sample collection and sequencing

To extract DNA, liquid samples were collected before and after each starvation phase, resulting in five distinct time points. On the final day (day 48), co-extraction of DNA and RNA was carried out following the manufacturer’s instructions, as summarized below. Prior to DNase treatment, aliquots were withdrawn from each sample to subsequently extract RNA. The remaining samples underwent DNAse treatment in accordance with the manufacturer’s protocol. For each extracted sample, quality and concentration were evaluated using Nanodrop2000 (ThermoFisher Scientific, Waltham, MA, USA) and Qubit fluorometers (ThermoFisher Scientific, Waltham, MA, USA), respectively. Furthermore, RNA quality was assessed using an Agilent 2100 Bioanalyzer with RNA 6000 Nano reagents and RNA Nano Chips (Agilent Technologies, Santa Clara, CA, USA), following the manufacturer’s instructions. Ribosomal RNA was removed using the FastSelect -5S/16S/23S Kit (QIAGEN GmbH, Germany). Library preparation utilized the Nextera DNA Flex Library Prep Kit (Illumina Inc, San Diego, CA), and sequencing was performed on the Illumina NovaSeq platform by the sequencing facility of the Department of Biology (University of Padova, Italy).

Metagenomics and gene-level mining

A genome-centric metagenomic pipeline was employed for the retrieval of microbial genomes. Following sequencing, Illumina reads underwent filtration using Trimmomatic (v0.39-1)⁷⁶ and BBMap (38.90) to remove adapters and low-quality reads. Assembly procedures were conducted with Megahit (v1.2.9)⁷⁷, while Bowtie2 (v2.4.4)⁷⁸ and SAMtools (v1.16)⁷⁹ facilitated the generation of contigs coverage profiles essential for retrieving MAGs. Multiple binning tools were employed, specifically MetaBAT (v1.2.15)⁸⁰, MetaBAT2 (v2.2.15)⁸¹, VAMB (v4.13)⁸², and Maxbin2 (v2.2.7)⁸³. The completeness and contamination of reconstructed MAGs were assessed using checkM2 (v1.0.2)⁸⁴. Subsequently, MAGs underwent aggregation and filtration utilizing dRep (v3.4.0)⁸⁵ to selectively retain medium-to-high quality MAGs in accordance with MIMAG standards. The final set of MAGs was processed using coverM (v0.6.1)⁸⁶ to determine their RA in each sample. GTDB-tk (v2.1.1)⁸⁷ with database r214 was employed for taxonomic classification assignment. Functional classification was achieved through eggNOG-mapper (v2.1.10)⁸⁸ by aligning translated protein sequences, obtained with Prodigal (v2.6.3)⁸⁹, against an internal database of ortholog groups. The Kyoto Encyclopedia of Genes and Genomes (KEGG)⁹⁰ served as the reference database for the metabolic functional reconstruction of the identified MAGs, focusing on specific pathways of interest such as hydrogenotrophic methanogenesis. A KEGG module was considered complete if, at most, only one block (representing a group of reactions) was missing. MAGs that were reconstructed received a customized nomenclature, which comprised the closest taxonomic level achieved, the suffix “MN_“ denoting the experiment, two letters designating the binning software utilized (MX: Maxbin2, MA: MetaBAT, MB: MetaBAT2, VB: Vamb), and an incremental number.

Genome-centric metatranscriptomics

RNA-seq reads underwent quality filtering using Trimmomatic (v0.39-1)⁷⁶ and BBDuk (39.01) before alignment onto MAGs with Bowtie2 (v2.4.4)⁷⁸. The htseq-count script from HTSeq (v2.0.3)⁹¹ was employed for counting the number of RNA fragments per gene in stranded mode during transcriptome data analysis. Fragment counts were normalized separately for each MAG, with size factors calculated using the Estimate Size Factors function. These size factors were applied to raw counts, resulting in normalized counts. Normalization was conducted by gene and scaled by a factor of 10 ^ 3, yielding Fragments Per Kilobase (FPK) values. Differential expression analysis with DESeq2⁹² was defined using the design formula “Reactor + Timepoint”. For each reactor, the log₂ fold change (log₂FC) was calculated by comparing the gene expression on day 12 with that measured on day 9. Analysis was conducted employing a Wald test and a false discovery rate threshold of 0.05, along with a log₂FC threshold of 1.

Metabolic reconstruction and flux-balance analysis

For computational metabolic reconstruction, MAGs that demonstrated a RA greater than 1e⁻⁴% were selected. The construction of GEMs from these MAGs was facilitated by the gapseq tool (v.1.2, sequence DB md5: 748b938), which involved choosing an appropriate database of metabolic reactions and adhering to standard gene-matching criteria to restrict the search within defined taxonomic limits⁹³. During draft GEMs reconstruction, taxonomic data were used to select the optimal biomass reaction according to gapseq parameters. The default bit score threshold was used to determine reactions having sequence support or not. Functional models were generated by gapfilling on the predicted minimal medium using the module gapseq medium. Following this, the micom platform (v.0.21.3) was employed to integrate species-specific GEMs into a community-level model that accounts for the metabolic exchange fluxes among individual species, as well as between the collective microbiome and its environment⁹⁴. The RA of the species was integrated in the reconstruction to reflect the composition of the community. Constraints at the community level were defined, taking into account the measured concentrations of VFA and the consumed feedstock (i.e. H₂ and CO₂). To determine rates, which are denoted in mmol * gDW⁻¹ * h⁻¹, estimations of dry weight and the working volume of the culture were used (Supplementary Table 13). Specifically, the final dry biomass was normalized according to the measured OD to infer the dry weight of the community at each sampling point. To account for variation due to measurement errors the constraints were imposed considering a degree of flexibility equivalent to two standard deviations around the measured value. Standard deviations were estimated as uncertainty propagation of systematic error of 8860 GC measurements. Flux-balance analysis was performed with the goal of optimizing microbiome biomass growth, applying a cooperative trade-off strategy (tradeoff = 0.7). Additionally, RNA data for the available sampling points were utilized to formulate condition-specific community GEMs (CoCo-GEMs), in line with previously established methods⁹⁵. Cytoscape (v.3.10.2) was used to represent the predicted metabolic exchanges⁹⁶. Only interplays between species with a RA above 0.1% and predicted fluxes above 0.09 mmol * gDW⁻¹ * h⁻¹ were represented.

SNVs analysis and filtering

The software InStrain (v1.6.3)⁹⁷ was applied to the reconstructed MAGs to perform variant calling. The InStrain profile module optional parameters were set as follows: -min_mapq 2, –min_read_ani 0.98 and –min_genome_coverage 1. The InStrain outputs underwent additional processing to filter out low-confidence variants. Initially, all SNVs located within 150 base pairs from the 3′ and 5′-end of each scaffold were excluded. This decision was made due to the tendency of coverage to decrease in those genomic regions, leading to less reliable results. Subsequently, variants were discarded if the difference between the SNV coverage and the average scaffold coverage fell outside the range of [–100; +100]. Throughout the subsequent analysis, SNVs shared among the CTRL, Ni, and Co conditions were examined separately. This approach aimed to distinguish between the effects of starvation and those of metal supplementation. Lastly, proteins encoded by genes affected by SNVs were modeled using AlphaFold⁹⁸ to predict the 3D structures of the mutated polypeptides, which were then compared to their original counterparts, as previously described⁹⁹. Protein sequences were aligned using the Clustal Omega tool, based on the ClustalW algorithm¹⁰⁰, against reference structures from Swiss-Prot. This comparison enabled the assessment of the potential impact of SNVs on the target proteins. The results were visualized using ChimeraX (v1.7.1)¹⁰¹, by superimposing the 3D structures of the wild-type and mutated proteins, with regions impacted by amino acid changes highlighted.

Estimation of replication rates

Species replication rates were estimated across 8 million read subsamples using CoPTR (v1.1.4) with default parameters¹⁰². This method calculates the slope of read coverage across genomes or MAGs, based on the assumption that rapidly replicating microbes exhibit a higher number of DNA copies close to the origin of replication. This slope, known as the PTR, quantifies DNA synthesis and generation rates in terms of a dimensionless coverage decay rate, serving as a proxy for microbial replication rates. To ensure reliable PTR estimates, only MAGs with a minimum of 5000 aligned reads were included in this analysis.

Limitations

The current study was designed to involve not only Ni and Co supplementation but also incorporate stepwise increments in starvation periods, which poses challenges in pinpointing the specific selective conditions that drive the evolution of the microbiome. Here an ad hoc SNVs investigation was employed to address the molecular determinants behind the fluctuations observed. Nonetheless, the possibility to directly link the identified SNVs to the strain of origin is still missing since current software works on MAGs and not strain-assembled genomes. Multiple approaches leveraging single-cell technologies are presently under development, founded on the concept that Single-Amplified Genomes exhibiting closely related SNVs can be regarded as belonging to the same strain. While this area of research is still considered nascent, its potential is noteworthy, offering a distinctive opportunity to enhance our comprehension of the interplay between microbial population diversity and the operational conditions of biomethanation systems. Moreover, the computational methods used to deduce the inter-species interactions could be subject to biases difficult to prevent. First, the reconstruction of GEMs for MAGs from extreme ecosystems is prone to gaps due to incomplete genomes or missing gene-protein-reaction annotations. This can result in the underestimation of a species’ true metabolic potential, which impacts community-level simulations by obscuring syntrophic interactions. Second, while FBA is a robust method for estimating inter-species metabolic exchanges, it falls short in modeling ecological dynamics that are potentially present in the analyzed community. Specifically, the role of vertical trophic chains, such as the consumption of the primary producer necromass by heterotrophs, was not explored. This is a key ecological process that may significantly affect metabolite exchange within the community and warrants independent experimental validation to assess its relevance.