Hydrogenotrophic methanogens of the mammalian gut: Functionally similar, thermodynamically different—A modelling approach

Summary

Methanogenic archaea occupy a functionally important niche in the gut microbial ecosystem of mammals. Our purpose was to quantitatively characterize the dynamics of methanogenesis by integrating microbiology, thermodynamics and mathematical modelling. For that, in vitro growth experiments were performed with pure cultures of key methanogens from the human and ruminant gut, namely Methanobrevibacter smithii, Methanobrevibacter ruminantium and Methanobacterium formicium. Microcalorimetric experiments were performed to quantify the methanogenesis heat flux. We constructed an energetic-based mathematical model of methanogenesis. Our model captured efficiently the dynamics of methanogenesis with average concordance correlation coefficients of 0.95 for CO₂, 0.98 for H₂ and 0.97 for CH₄. Together, experimental data and model enabled us to quantify metabolism kinetics and energetic patterns that were specific and distinct for each species despite their use of analogous methane-producing pathways. Then, we tested in silico the interactions between these methanogens under an in vivo simulation scenario using a theoretical modelling exercise. In silico simulations suggest that the classical competitive exclusion principle is inapplicable to gut ecosystems and that kinetic information alone cannot explain gut ecological aspects such as microbial coexistence. We suggest that ecological models of gut ecosystems require the integration of microbial kinetics with nonlinear behaviours related to spatial and temporal variations taking place in mammalian guts. Our work provides novel information on the thermodynamics and dynamics of methanogens. This understanding will be useful to construct new gut models with enhanced prediction capabilities and could have practical applications for promoting gut health in mammals and mitigating ruminant methane emissions.

Keywords:

Carbon dioxide – Hydrogen – Mathematical models – Thermodynamics – Ribosomal RNA – Fermentation – Methane – Methanogens

Introduction

Methanogenic archaea inhabit the gastro-intestinal tract of mammals where they have established syntrophic interactions within the microbial community [1–3] playing a critical role in the energy balance of the host [4,5]. In the human gut microbiota, the implication of methanogens in host homeostasis or diseases is poorly studied, but of growing interest [6]. Methanobrevibacter smithii (accounting for 94% of the methanogen population) and Methanosphaera stadtmanae are specifically recognized by the human innate immune system and contribute to the activation of the adaptive immune response [7]. Decreased abundance of M. smithii was reported in inflammatory bowel disease patients [8], and it has been suggested that methanogens may contribute to obesity [9]. In the rumen, the methanogens community is more diverse though still dominated by Methanobrevibacter spp., followed by Methanomicrobium spp., Methanobacterium spp. [10] and Methanomassillicoccus spp [11]. However, the proportion of these taxa could vary largely, with Methanomicrobium mobile and Methanobacterium formicium being reported as major methanogens in grazing cattle [12]. Though methanogens in the rumen are essential for the optimal functioning of the ecosystem (by providing final electron acceptors), the methane they produce is emitted by the host animal, contributing to global greenhouse gas (GHG) emissions. In the gastrointestinal tract of mammals, major rumen methanogens [13] and the dominant human archaeon M. smithii [14], are hydrogenotrophic archaea without cytochrome (membrane-associated electron transfer proteins). Cytochrome-lacking methanogens exhibit lower growth yields than archaea with cytochromes [15]. However, this apparent energetic disadvantage has been counterbalanced by a greater adaptation to the environmental conditions prevailing in the gastrointestinal tract [16], and by the establishment of syntrophic interactions with feed fermenting microbes. This syntrophic cooperation centred on hydrogen allows the anaerobic reactions of substrate conversion to proceed close to the thermodynamic equilibrium [17,18] (that is with Gibbs free energy change close to zero).

To our knowledge, thermodynamic considerations on human gut metabolism have been poorly addressed in existing mathematical models [19–22], although ttheoretical frameworks have been developed in other domains to calculate stoichiometric and energetic balances of microbial growth from the specification of the anabolic and catabolic reactions of microbial metabolism [23,24], and advances have been done to link thermodynamics to kinetics [25–27]. For the rumen, thermodynamic principles have been incorporated already into mathematical research frameworks because of their important role in feed fermentation. Thermodynamic studies have been performed to investigate theoretically (i) the profile of fermentation [28], (ii) alternative routes for hydrogen utilization [29], and (iii) the effect of hydrogen partial pressure on glucose fermentation and methanogenesis [30,31]. However, studies assessing quantitative comparisons between experimental data and model predictions using thermodynamic-based approaches have been of limited success in providing accurate predictions [32,33], probably due to missing controlling factors such as NADH oxidation and the dynamics of hydrogen partial pressure [31]. Another key factor explaining inaccurate predictions in existing rumen models is the lack of a dynamic representation of the microbial methanogens group. In this respect, new knowledge on the extent of methanogenesis and metabolic differences between the microbial members of this group could help to improve existing gut models. To our knowledge, none of the existing gut models integrates thermodynamics aspects to describe microbial growth of methanogens. Accordingly, our objective in this work was to develop a dynamic model with thermodynamic basis to quantify metabolism kinetics and energetic patterns that could inform on metabolic specificities between gut methanogenic archaea. Additionally, this model development aimed at providing tools for the analysis of ecological aspects (e.g., competitive exclusion principle) of the methanogenic community that can be instrumental when designing, for example, nutritional strategies for methane mitigation in ruminants. The model was built upon quantitative data characterizing the in vitro dynamics of hydrogen utilization, methane production, growth and heat flux of three hydrogenotrophic methanogenic species representing major human and ruminant genera: Methanobrevibacter smithii, Methanobrevibacter ruminantium and Methanobacterium formicium.

Material and methods

In vitro growth experiments

Archaeal strains and growth media

Archaeal strains used in the study were M. ruminantium M1 (DSM 1093), M. smithii PS (type strain DSM 861), and M. formicium MF (type strain DSM 1535). All archaeal strains were purchased from DSM previously. The strains were frozen at -80°C, thawed and grown to get actively growing population before the beginning of the study. Balch growth media was prepared as previously described [34] and composition is summarized in S1 Table. Before autoclaving, 6 ml of media were distributed in Balch tubes (26 ml total volume) under CO₂ atmosphere.

Experimental design and measures

Starter cultures were grown until reaching optical density at 660 nm (OD₆₆₀) of 0.400 ± 0.030. Optical density was measured on a Jenway spectrophotometer (Bibby Scientific). Then, 0.6 ml were used to inoculate one experimental tube. Commercially prepared high purity H₂/CO₂ (80%/20%) gas mix was added to inoculated tubes by flushing for 1 min at 2.5 Pa. Mean initial OD₆₆₀ and pressure values are summarized in S2 Table. Growth kinetics for each strain were followed over 72 h. The experiment was repeated twice. Each kinetics study started with 40 tubes inoculated at the same time. At a given time point, two tubes with similar OD₆₆₀ values were sampled. The tubes were used for measuring gas parameters: pressure was measured using a manometer and composition of the gas phase was analysed by gas chromatography on a Micro GC 3000A (Agilent Technologies, France). The GC was equipped with two columns, MS-5A using argon as carrier gas and set to 100°C and PPU using helium and set to 75°C. The GC was calibrated using a certified gas standard mixture (Messer, France) containing methane, oxygen, hydrogen, carbon dioxide, and nitrogen. Approximately 2 ml of the sampled gas was injected in the GC for analysis. After gas sampling, one of the tubes was centrifuged 10 min at 13 000 g. The microbial pellet was weighed and stored at -20°C in 2 ml screw-cap tubes containing 0.4 g of sterile zirconia beads (0.3 g of 1 mm and 0.1 g of 0.5 mm).

DNA extraction and qPCR quantification of 16S rRNA genes

One ml of lysis buffer (50mM NaCl, 50 mM TrisHCl pH 7.6, 50 mM EDTA, 5% SDS) was added directly to the frozen microbial pellet before homogenizing for 2 × 30 s at 5100 tours/min in a Precellys bead-beater (Bertin Instruments). Samples were centrifuged for 3 min at 14 000 g and the liquid phase transferred to a new tube before adding 600 μl of phenol–chloroform–3-methyl-1-butanol (25:24:1) solution. After centrifugation at 14 000 g for 3 min, the aqueous phase was transferred to a fresh tube and 500 μl of chloroform were added. The chloroform-washing step was repeated twice with centrifugation at 14000 g for 3 min between steps. The final volume of the aqueous phase was measured and DNA precipitation was initiated by adding 70% of the volume of isopropanol 100% and 10% of the volume of sodium acetate 3M. Sedimentation at 14 000 g for 30 min was again performed and the resulting DNA pellet was washed with 500 μl of ethanol 70% and dissolved in 50μl of molecular biology quality water. The extraction yield was checked on a Nanodrop 1000 Spectrophotometer (Thermo Fisher Scientific, France) and extracts run on a FlashGel System (Lonza, Rockland, Inc) to check integrity.

Copies of 16S rRNA genes were quantified using a qPCR approach. Primers used are those of Ohene-Adjei et al [35]; reaction assay and temperature cycles were as described previously [36]. Triplicate qPCR quantification was performed on 20 ng of extracted DNA. Amplifications were carried out using SYBR Premix Ex Taq (TaKaRa Bio Inc., Otsu, Japan) on a StepOne system (Applied Biosystems, Courtabeuf, France). Absolute quantification involved the use of standard curves that had been prepared with gDNA of Methanobrevibacter ruminantium DSM 1093. PCR efficiency was of 103%. Results were expressed as copy numbers per ng of extracted DNA per g of microbial pellet. M. smithii and M. ruminantium strains used in this study possess two copies of 16S rRNA genes in their genomes. The number of cells was computed by dividing 16S copy numbers by 2.

Microcalorimetry

Microcalorimetric experiments were performed to determine the heat flux pattern of each methanogen. Metabolic activity and microbial growth were monitored by using isothermal calorimeters of the heat-conduction type (A TAM III, TA Instruments, France) equipped with two multicalorimeters, each holding six independent minicalorimeters, allowed continuous and simultaneous recording as a function of time of the heat flux produced by 12 samples. The bath temperature was set at 39°C; its long-term stability was better than ± 1x10^-4°C over 24h. Each minicalorimeter was electrically calibrated. The specific disposable 4 mL microcalorimetric glass ampoules capped with butyl rubber stoppers and sealed with aluminium crimps were filled with 1.75 mL of Balch growth media and overpressed with 2.5 Pa of H₂/CO₂ 80%/20% gas mixture for 30 s. There was no significant difference in pressure at the beginning of the study. They were sterilized by autoclave and stored at 39°C until the beginning of the microcalorimetric measurements. Actively growing cultures of methanogens (OD₆₆₀ of 0.280±0.030 for M. smithii, 0.271±0.078 for M. ruminantium and 0.142±0.042 for M. formicium) were stored at -20°C in order to diminish microbial activity before inoculation. Cultures were thawed for 30 min at ambient temperature and inoculation was carried out by injecting 0.25 mL of the culture through the septum of the overpressed microcalorimetric ampoules just before inserting them into the minicalorimeters. Samples took about two hours to reach the bath temperature and yield a stable zero baseline. Blank experiments were also carried out by inserting ampoules that were not inoculated and, as expected, no heat flux was observed confirming the medium sterility. Each experiment was repeated thrice.

The heat flux (dQdt), also called thermal power output P, was measured for each methanogen and blank samples with a precision ≥ 0.2 μW. The heat flux data of each sample were collected every 5 minutes during more than 10 days. The total heat Q was obtained by integrating the overall heat flux time curve using the TAM Assistant Software and its integrating function (TA Instruments, France).

Classically, the heat flux-time curve for a growing culture starts like the S-shaped biomass curve (a lag phase followed by an exponential growth phase) but differs beyond the growth phase, the heat flux being then modulated by transition periods [37]. Heat flux data can be used to infer the microbial growth rate constant provided the existence of a correlation between isothermal microcalorimetry data and microbiological data (e.g., cell counts) at early growth [38]. During the exponential growth phase, microbial growth follows a first-order kinetics defined by the specific growth rate constant μ_c (h^-1). Analogously, the heat flux follows an exponential behaviour determined by the parameter μ_c as described by [37,38].

The growth rate constant μ_c can be determined by fitting the exponential part of the heat flux-time curve using the fitting function of the TAM Assistant Software. In our case study, careful selection of the exponential phase of heat flux dynamics was performed to provide a reliable estimation of the maximum growth rate constant from calorimetric data.

Mathematical model development

Modelling in vitro methanogenesis

The process of in vitro methanogenesis is depicted in Fig 1. The H₂/CO₂ mixture in the gas phase diffuses to the liquid phase. The H₂ and CO₂ in the liquid phase are further utilized by the pure culture to produce CH₄. Methane in the liquid phase diffuses to the gas phase.

Model construction was inspired by our previous dynamic models of human gut [19] and rumen in vitro fermentation [39] followed by certain simplifications. The model considers the liquid-gas transfer of carbon dioxide. Due to the low solubility of hydrogen and methane [40], the concentration of these two gases in the liquid phase was not modelled. We assumed that the dynamics of concentrations in the gas phase are determined by kinetic rate of the methanogenesis. To incorporate thermodynamic information, instead of using the Monod equation in the original formulation, we used the kinetic rate function proposed by Desmond-Le Quéméner and Bouchez [26]. The resulting model is described by the following ordinary differential equations

In our model, the stoichiometry of methanogenesis is represented macroscopically by one catabolic reaction (R₁) for methane production and one anabolic reaction (R₂) for microbial formation. We assumed that ammonia is the only nitrogen source for microbial formation. The molecular formula of microbial biomass was assumed to be C₅H₇O₂N [40].

In the model, the stoichiometry of the reactions is taken into account via the parameters Y, YCO2, YCH4, which are the yield factors (mol/mol) of microbial biomass, CO₂ and CH₄. The microbial yield factor Y was extracted from literature. We assumed that M. smithii and M. ruminantium have the same yield (being both Methanobrevibacter). This yield factor was set to 0.006 mol biomass/mol H₂, using the methane-based molar growth yield of 2.8 g biomass/mol CH₄ estimated for M. smithii [41] and the Eqs (9) and (11). Similarly, the yield factor for M. formicium was set to 0.007 mol biomass/mol H₂ using the methane-based molar growth yield of 3.5 g biomass/mol CH₄ reported by Schauer and Ferry [42]. The fraction of H₂ utilized for microbial growth (reaction R₂) is defined by the yield factor Y (mol of microbial biomass/mol of H₂). Now, let f be the fraction of H₂ used for the catabolic reaction R₁. Reaction R₂ tells us that for every 10 moles of H₂ used in R₂, we get 1 mol of microbial biomass. Hence, it follows that

If Y is known, the fraction f can be computed from Eq (9).

The yield factors of CO₂ and CH₄ can be expressed as functions of the the fraction f:

The model has two physicochemical parameters (k_La, KH,CO2) and four biological parameters (μ_max, K_s, Y, k_d). The initial condition for sCO2 is unknown and was also included in the parameter vector for estimation. The Henry’s law coefficients are known values calculated at 39°C using the equations provided by Batstone et al. [40]. An implementation of the model in the Open Source software Scilab (https://www.scilab.org/) is available at https://doi.org/10.5281/zenodo.3271611.

Theoretical model to study interactions among methanogens

The experimental study of microbial interactions requires sophisticated in vitro systems under continuous operation such as the one developed by Haydock et al. [43]. In our work, we explored by means of mathematical modelling how the methanogens can interact under in vivo conditions. For this theoretical study, we elaborated a toy model based on the previous model for in vitro methanogenesis. Let us consider the following simple model for representing the consumption of hydrogen by the methanogenic species i under an in vivo scenario of continuous flow

Parameter identification

Before performing the numerical estimation of the model parameters, we addressed the question of whether it was theoretically possible to determine uniquely the model parameters given the available measurements from the experimental setup. This question is referred to as structural identifiability [44]. Structural identifiability analysis is of particular relevance for model whose parameters are biologically meaningful, since knowing the actual value of the parameter is useful for providing biological insight on the system under study [45]. Moreover, in our case, we are interested in finding accurate estimates that can be further used as priors in an extended model describing the in vivo system.

We used the freely available software DAISY [46] to assess the structural identifiability of our model. Physical parameters (k_La, KH,CO2) were set to be known. The model was found to be structurally globally identifiable. In practice, however, to facilitate the actual identification of parameters and reduce practical identifiability problems such as high correlation between the parameters [47], we fixed some model parameters to values reported in the literature. The transport coefficient k_La, the Henry’s law coefficient KH,CO2, and the dead cell rate constant k_d were set to be known and were extracted from Batstone et al. [40]. Therefore, only the parameters μ_max, K_s and initial condition of sCO2 were set to be estimated. To capitalize on the calorimetric data, we further assumed that μ_max was equal to the specific rate constant μ_c estimated from the heat flux-time curve. By this, only the affinity constant for each strain and the initial condition of sCO2 were left to be estimated.

The parameter identification for each methanogen was performed with the IDEAS Matlab^® toolbox [48] (freely available at http://genome.jouy.inra.fr/logiciels/IDEAS). The parameter identification was performed with the data of the in vitro growth experiments. The measured variables are the number of moles in the gas phase (H₂, CH₄, CO₂). The Lin’s concordance correlation coefficient (CCC) [49] was computed to quantify the agreement between the observations and model predictions.

Results

Methanogens biomass

Archaea-specific primers targeting the 16S rRNA gene were used to enumerate microbial cells in each pure culture. Three hours post inoculation microbial numbers varied from 7.62×10⁷ to 2.81×10⁸ and reached 10⁹ after 72 hours of incubation. S3 Table summarizes microbial numbers at different sampling times.

Calorimetric pattern of methanogens

Fig 2 displays a representative isothermal calorimetric curve for each methanogen. The three measured heat flux dynamics of each methanogen were found to follow similar energetic patterns. M. smithii and M. formicium exhibited a lag phase of a few hours, while M. ruminantium was already metabolically active when introduced into the minicalorimeter though several attempts were made to obtain a lag phase by changing storage conditions and thawing the culture just before inoculating the microcalorimetric ampoules. The pattern of heat flux for all tested methanogens is characterized by one predominant peak which was observed at different times for each methanogen. M. smithii exhibited a second metabolic event occurring at 60 h with an increase of heat flux. The same phenomenon was observed for M. formicium but at a lower intensity that started at 140 h. The process was considered completed when the heat flux ceased marking the end of the metabolic activity. It is noted that M. formicium produced a small peak at 14 h (Fig 2). A similar peak, but of much smaller size, was observed on the other curves obtained with this methanogen. M. smithii also exhibits a small peak (occurrence of 3 out of 3) at 7.4 h shown in the inset of Fig 2.

The total heat (Q_m) produced during the methanogenesis process that took place under the present experimental conditions was, on average, -5.5 ± 0.5 J for the three methanogens (for M. ruminantium, the missing initial part of the heat flux-time curve was approximately estimated by extrapolating the exponential fit). As we shall see below, this experimental value is consistent with the theoretically expected value.

Estimation of thermodynamic properties

We defined two macroscopic reactions to represent the catabolism (R1) and anabolism (R2) of the methanogenesis (see Modelling in vitro methanogenesis section). All thermodynamic properties result from the contribution of both catabolic and anabolic reactions. The calculations of total heat (Q_m), enthalpy (ΔH_m), Gibbs energy (ΔG_m) and entropy (ΔS_m) of the methanogenesis are detailed in S4 Table. The estimated overall heat produced during the methanogenesis process under our experimental conditions was in average Q_m = −5.66J. This heat results from the sum of the heat of the catabolic reaction (Q_c) and the heat of the anabolic reaction (Q_a). From the total heat of the methanogenesis, the anabolic reaction contributes to 7% of the metabolic heat for M. smithii and M. ruminantium. For M. formicium, the contribution of the anabolic reaction to the metabolic heat is 9%. It is also interesting to note that there is a very good agreement between the theoretical value calculated above and the overall heat experimentally determined by microcalorimetry (-5.5 ± 0.5 J).

Table 1 shows the thermodynamic properties per mole of biomass formed during methanogenesis of M. ruminantium, M. smithii and M. formicium on H₂/CO₂. These properties are compared with values found in the literature for other methanogens grown on different substrates.

Dynamic description of in vitro kinetics

The developed mathematical model was calibrated with the experimental data from in vitro growth experiments in Balch tubes. Table 2 shows the parameters of the dynamic kinetic model described in Eqs 2–6. The reported value of μ_max for each methanogen corresponds to the average value obtained from three heat flux-time curves. From Table 2, it is concluded that M. smithii exhibited the highest growth rate constant, followed by M. ruminantium and finally M. formicium. In terms of the affinity constant K_s, while M. smithii and M. ruminantium are of the same order, the affinity constant for M. formicium is lower in one order of magnitude.

Fig 3 displays the dynamics of the compounds in the methanogenesis for the three methanogens. The data of the figure is available at https://doi.org/10.5281/zenodo.3469655. Experimental data are compared against the model responses. Table 3 shows standard statistics for model evaluation. The model captures efficiently the overall dynamics of the methanogenesis. Hydrogen and methane are very well described by the model with average concordance correlation coefficients (CCC) of 0.98 and 0.97 respectively. For carbon dioxide, CCC = 0.95.

Fig 4 displays the dynamics for the methanogens as measured by the 16S rRNA gene, as well as the dynamics of biomass as predicted for the model. As observed, the microbes follow a typical Monod-like trajectory.

Discussion

Our objective in this work was to quantitatively characterize the dynamics of hydrogen utilization, methane production, growth and heat flux of three hydrogenotrophic methanogens by integrating microbiology, thermodynamics and mathematical modelling. Our model developments were instrumental to quantify energetic and kinetic differences between the three methanogens studied, strengthening the potentiality of microcalorimetry as a tool for characterizing the metabolism of microorganisms [52]. This modelling work provides estimated parameters that can be used as prior values for other modelling developments of gut microbiota.

Energetic and kinetic differences between methanogens

Methanogenesis appears as a simple reaction described by a single substrate kinetic rate on H₂. The microcalorimetry approach we applied revealed that this simplicity is only apparent and that hydrogenotrophic methanogens exhibit energetic and kinetic differences. Methanogenesis is indeed a complex process that can be broken down in several stages. The dominant metabolic phase is represented by one peak that occurs at different times. The magnitude of the peak differs between the methanogens and also the slope of the heat flux trajectories. For M. smithii and M. formicium the main peak was preceded by a small increase in heat flux which translates in a metabolic activity that remains to be elucidated. For M. ruminantium, we do not know whether the small peak exists since the initial part of the curve is missing. For M. smithii and M. formicium, it was observed a metabolic event after the main peak around 60 h for M. smithii and 140h M. formicium. This peak could be representing cell lysis process [38]. The return time of the heat flux to the zero baseline was also different. The energetic difference is associated with kinetic differences that translate into specific kinetic parameters, namely affinity constant (K_s) and maximum growth rate constant (μ_max). Previously, energetic differences between methanogens have been ascribed to the presence or absence of cytochromes [15]. These differences are translated into different yield factors, H₂ thresholds, and doubling times. The kinetic differences revealed in this study for three cytochrome-lacking methanogens indicate that factors other than the presence of cytochromes might play a role in the energetics of methanogenesis. Interestingly, calorimetric experiments showed that M. ruminantium was metabolically active faster than the other methanogens, characteristic that could explain the predominance of M. ruminantium in the rumen [53]. Looking at the expression of the affinity constant (Eq (8)), the differences between the affinity constants among the methanogens can be explained by the differences between the by the harvest volume v_harv and the yield factors. Note that in the kinetic function developed by Desmond-Le Quéméner and Bouchez [26], the maximum growth rate did not have any dependency on the energetics of the reaction. Our experimental study revealed that μ_max is species-specific and reflects the dynamics of the heat flux of the reaction at the exponential phase. This finding suggests that a further extension of the kinetic model developed by Desmond-Le Quéméner and Bouchez [26] should include the impact of energetics on μ_max. Since our study is limited to three species, it is important to conduct further research on other methanogens to validate our findings. In this same line, to enhance the evaluation of the predictive capabilities of our model, a further model validation is required with independent data set under different experimental conditions (e.g. continuous mode operation) to those used in this study.

Energetic analysis

Regarding the energetic information for different methanogens summarized in Table 1, it is observed that the thermodynamic behaviour of the three methanogens is analogous to that observed for Methanobacterium thermoautotrophicum [50]. The values reported in Table 1 show indeed that the methanogenesis on H₂/CO₂ is characterized by large heat production. The growth is highly exothermic, with a ΔH_m value that largely exceeds the values found when other energy substrates are used. The enthalpy change ΔH_m, which is more negative than the Gibbs energy change ΔG_m, largely controls the process. Growth on H₂/CO₂ is also characterized by a negative entropic contribution TΔS_m which, at first sight, may look surprising since entropy increases in most cases of anaerobic growth [54]. However, this can be understood if one remembers that TΔS_m corresponds in fact to the balance between the final state and the initial state of the process, that is

Methanogenesis on H₂/CO₂ is particular because the final state of its catabolic reaction (1 mol CH₄ + 2 mol H₂O) involves a smaller number of moles than the initial state (4 mol H₂ + 1 mol CO₂), which results in a significant loss of entropy during the process. For spontaneous growth in such a case, the ΔH_m must not only contribute to the driving force but must also compensate the growth-unfavourable TΔS_m, which means that ΔH_m must be much more negative than ΔG_m [55]. For this reason, methanogenesis on H₂/CO₂, which is accompanied by a considerable decrease of entropy and a large production of heat, has been designed as an entropy-retarded process [50]. More generally, von Stockar and Liu [55] noticed that when the Gibbs energy of the metabolic process is resolved into its enthalpic and entropic contributions, very different thermodynamic behaviours are observed depending on the growth type. These thermodynamic behaviours are: aerobic respiration is clearly enthalpy-driven (ΔH_m ≪ 0 and TΔS_m > 0), whereas fermentative metabolism is mainly entropy-driven (ΔH_m < 0 and TΔS_m ≫ 0). Methanogenesis on H₂/CO₂ is enthalpy-driven but entropy-retarded (ΔH_m ≪ 0 and TΔS_m < 0), whereas methanogenesis on acetate is entropy-driven but enthalpy-retarded (ΔH_m > 0 and TΔS_m ≫ 0). In the present case, the highly exothermic growth of M. ruminantium, M. smithii and M. formicium on H₂/CO₂ is largely due to the considerable decrease of entropy during the process: in fact, 50% of the heat produced here serves only to compensate the loss of entropy. A proportion of 80% was found for M. thermoautotrophicum [50], which results from the fact that their TΔS_m and ΔH_m values are, respectively, 2.7 and 1.7 times larger than ours. This difference might be due to the differences in the temperature of the studies, namely 39°C in our study vs 60°C in the study by Schill et al. [50].

Do our results inform on ecological questions such as species coexistence?

The competitive exclusion principle [56] states that coexistence cannot occur between species that occupy the same niche (the same function). Only the most competitive species will survive. Recently, by using thermodynamic principles, Großkopf & Soyer [27] demonstrated theoretically that species utilizing the same substrate and producing different compounds can coexist by the action of thermodynamic driving forces. Since in our study the three methanogens perform the same metabolic reactions, the thermodynamic framework developed Großkopf & Soyer [27] predicts, as the original exclusion principle [56], the survival of only one species. By incorporating thermodynamic control on microbial growth kinetics, Lynch et al [57] showed theoretically that differentiation of ATP yields can explain ecological differentiation of methanogens over a range of liquid turnover rates. This theoretical work predicts that for a fixed liquid turnover rate, only one species survives. For the continuous culture of microorganisms, it has been demonstrated that at the equilibrium (growth rate equals the dilution rate) with constant dilution rates and substrate input rates, the species that has the lowest limiting substrate concentration wins the competition. From Eq (12), the number of moles of hydrogen of the species ng,H2,i* at the steady state is

Using the model parameters of Table 2, we studied in silico three possible competition scenarios, assuming a constant environment (constant dilution rate D). Two dilution rates were evaluated: D = 0.021 h^-1 (retention time = 48 h) and D = 0.04 h^-1 (retention time = 25 h). A retention time of 48 h corresponds to values measured in small ruminants [58] and to humans as we used in our gut model [19]. For higher retention times, the results obtained for 48 h hold. For D = 0.021 h^-1, we obtained that ng,H2,Ms*=0.32mmol,ng,H2,Mr* = 0.68 mmol, ng,H2,Mf*=0.28mmol, where the subindex Ms, Mr, Mf stand for M. smithii, M. ruminantium and M.formicium. From these results, it appears that under a constant environment, M. formicium will win the competition. Since ng,H2,Ms*<ng,H2,Mr*, M. ruminantium will be extinguished before M. smithii. For D = 0.04 h^-1, we obtained that ng,H2,Ms*=0.49mmol, ng,H2,Mr*=1.42mmol, ng,H2,Mf*=1.57mmol, and thus M. smithii wins the competition. To win the competition, M. ruminantium requires longer retention times than its competitors. Retention times of digesta longer than 48 h are physiologically uncommon, thus the presence of M. ruminantium in the gut ecosystem can be explained, for example, from known adhesion properties (both M. ruminantium and M. smithii genes encode adhesin-like proteins [59,60]. To illustrate these aspects, we built a multiple-species model with the three methanogens using Eqs (12) and (14). The parameter b was set to 0.5 h^-1 and the hydrogen flux production qH2 rate was set to 0.02 mol/min. Fig 5A displays the dynamics of the three methanogens for the first scenario (D = 0.021 h^-1). It is observed that at 50 d only M. formicium survives. This result, however, is not representative of what occurs in the rumen where the three methanogens coexist [5,61]. It is intriguing that in our toy model it is M. formicium that wins the competition, bearing in mind that M. ruminantium and M. smithii are more abundant than M. formicium [5,53]. Fig 5 shows that selective conditions favour the survival of one species. Similar results can be obtained for the human gut by including the effect of pH on microbial growth [22] and setting the gut pH to select one of the species.

On the basis of the competitive exclusion principle, it is thus intriguing that having a very specialized function, methanogens are a diverse group that coexist. Gut ecosystems, therefore, exhibit the paradox of the plankton introduced by Hutchinson (1961) that presents the coexistence of species all competing for the same substrate in a relatively isotropic or unstructured environment [62]. In the case of the rumen, our modelling work suggests that in addition to kinetic and thermodynamic factors, other forces contribute to the ecological shaping of the methanogens community in the rumen favouring the microbial diversity. Indeed, methanogenic diversity in the rumen results from multiple factors that include pH sensitivity, the association with rumen fractions (fluid and particulate material), and the endosymbiosis with rumen protozoa [5,53]. For the human gut, ecological factors enable methanogens to coexist to a competitive environment where hydrogenotrophic microbes (acetogens, methanogenic archaea and sulfate-reducing bacteria) utilize H₂ via different pathways [63–65]. Both in the human gut and in the rumen, microbes grow in association with biofilms that form a polymer-based matrix that provides nutritional and hydraulic advantages for microbial growth and resistance to shear forces [19,66]. Indeed, in our modelling work of human gut fermentation [19], we suggested that, from the different actions the mucus has on colonic fermentation, the mechanism of promoting conditions for microbial aggregation appears as the most relevant factor for attaining the high microbial density and the high level of fibre degradation characteristic of the human gut. Altogether, these factors result in nonlinear behaviours, spatial and temporal variations that promote coexistence and diversity, that, as discussed in dedicated literature on microbial ecology [67–71], render the classical formulation of the competitive exclusion principle [56,72] inapplicable to gut ecosystems.

Finally, mathematical modelling is expected to enhance our understanding of gut ecosystems [66,73]. It is then key that in addition to metabolic aspects, mathematical models of gut fermentation incorporate the multiple aspects that shape microbial dynamics to provide accurate predictions and improve insight on gut metabolism dynamics and its potential modulation. For ruminants, the development of precision livestock technologies provides promising alternatives for integrating real-time data of key animal phenotypes such as feeding behaviour with mathematical models for estimating methane emissions [74] and rumen function indicators at large scale. These tools will be instrumental to support livestock management decisions and guide timely interventions. Similarly, for humans, mathematical models coupled with electronic technologies for online monitoring of gut function [75] might facilitate the diagnosis and the design of personalized therapies for gastrointestinal diseases.

Supporting information

S1 Table [docx]
Methanogens growth media composition.

S2 Table [docx]
Summary of initial OD and pressure measured immediately after primary inoculation.

S3 Table [docx]
qPCR quantification of 16S rRNA genes.

S4 Table [docx]
Calculation of thermodynamic properties of the methanogenesis.