Methane from microbial hydrogenolysis of sediment organic matter before the Great Oxidation Event

Model A: abundances and δ¹³C of short alkanes

Considering the cleavage at position m (between the no. m and no. m + 1 carbon atoms) of an n-alkyl chain with n carbon atoms (1 ≤ m < n), the hydrogenolysis reaction equation is

Similarly, the reaction equation for cleavage on an n-alkane molecule is

Branched alkyl chains are omitted for simplicity. The cleavage reactions are consecutive; the products R-C_mH_2m+1, H-C_mH_2m+1, and H-C_n – mH_2(n–m)+1 (written below as RC_m, HC_m, and HC_n−m) are still subject to hydrogenolysis until all C–C bonds break down with CH₄ as the ultimate product. The cleavage rate (r) of a kerogen side chain (RC_m) or an alkane molecule (HC_m) is proportional to the C–C chain length and concentration c:

Here, k is the reaction constant to break any C–C bond. For simplicity, the difference in k between a middle and an end bond in the C–C chains is not considered. The net reaction rate of a kerogen side chain with a length of m carbon atoms is

Here, t is time and N is the maximum chain length of the reaction system. The first term of the right-hand side accounts for the cleavage of the kerogen side chain, and the second term accounts for the generation of residual shorter side chains from the cleavage of the longer side chain. The net reaction rate of a normal alkane with m carbon atoms is

The first term of the right-hand side accounts for the cleavage of this alkane, the second term accounts for the generation of the alkane from the kerogen side chain, and the third term accounts for the generation from the cracking of alkanes longer than this alkane. The factor in the last term is necessary because HC_m is generated from cuttings at the m and the n − m positions of HC_n.

Equations (6) and (7) show that shorter chains become more enriched as the cracking goes on. On the one hand, longer chains are more prone to cracking because they have more C–C bonds; on the other hand, shorter chains are the products of longer chains.

Both primary and secondary ¹³C KIEs on the rate constant k are considered. Suppose that there is a ¹³C substitution at position j (1 ≤  j ≤ n); if j = m − 1 or j = m, then there is a primary ¹³C KIE, or, if j = m − 2 or j = m + 1, there is a secondary ¹³C KIE on the cleavage at the m position.

A normal distribution is used for the molar distribution of the lengths of initial kerogen side chains, with a minimum chain length of n_min, a maximum chain length of n_max, a mean length of (n_min + n_max)/2, and a standard deviation of σ. When the mean length is large enough (>15), the isotopic compositions of gas products are insensitive to the initial kerogen side chain length distribution. For initial values, a δ¹³C value of −35‰ is applied. The initial chain length is in a normal distribution with a peak of C₁₇ and a standard deviation of σ = 2 carbon atoms. The initial alkane concentrations are assumed to be 0.

For simplicity, we assume that since there is no isotopic fractionation within or between the alkyl chains at the beginning of hydrogenolysis, the probability of ¹³C substitution at any position of any side chain is identical and determined by the initial carbon isotopic composition δ¹³C. Multiple ¹³C substitutions on a C–C chain are omitted because consideration of multiple substitutions would drastically increase the modelling complexity. This approximation is valid when the C–C chain is not too long. For example, the ratio between the probabilities of double and single ¹³C substitution in a C₂₀ chain is ({left[{left(begin{array}{c}20 2end{array}right)}{{left(frac{{,}^{13}{{{{{rm{C}}}}}}}{{,}^{12}{{{{{rm{C}}}}}}}right)}}^{2}right]}/{left[{left(begin{array}{c}20 1end{array}right)}{left(frac{{,}^{13}{{{{{rm{C}}}}}}}{{,}^{12}{{{{{rm{C}}}}}}}right)}right]}) ≈ 10% for ¹³C/¹²C ~ 0.01. Such a chain is long enough that the δ¹³C of gas products is insensitive to C–C chain length. Numerical simulation was conducted with Mathworks MATLAB 2020a.

Model B: bulk and clumped isotopic fractionations of CH₄

Conversion of methylene in a long C–C chain to methane is generalised into two steps:

The first step (step a) is the conversion of the methylene group R-CH₂-R′ to a methyl group (RCH₃) by accepting a capping hydrogen atom from the hydrogen donor (activated H₂); the second step (step b) is the conversion of the methyl group to methane by accepting another capping hydrogen atom. This scheme is highly generalised, and each step may involve multiple elementary biochemical reaction steps, such as the binding of H₂ and long alkyl chains to the enzyme, activation of H–H and C–C bonds, and release of the short alkane products from the enzyme. It is beyond the scope of this work to discuss the detailed biochemical reaction steps. But the cleavage and formation of chemical bonds in these steps should be constrained by the observed isotopic patterns.

Due to the computational complexity, we did not use the random scission model (Model A) in the simulation involving clumped isotopic fractionation, as explained in the following. A conventional kinetic model of the decomposition of organic matter without considering the constraints of C–C chain lengths is a zero-dimensional problem. Modelling the random cutting of long C–C chains without considering isotopes is a one-dimensional problem, and modelling bulk carbon isotopic fractionation during random cutting (Model A) is a two-dimensional problem. If ¹³C–¹³C coupling is included in random cutting, the modelling is a three-dimensional problem; a complex Monte Carlo method has been applied to deal with this problem¹⁹. If the ¹³C–D or D–D coupling is included in Model B, as we wish, it is a problem above the fourth dimension. The complexity of programming and the difficulty of computation make the model unattainable; even if it is achievable, solving this problem is far beyond the scope of this work.

Reaction equation Eq. 8 is expanded to the scheme in Fig. 3a to quantify the five most abundant isotopologues in methane (three or more substitutions such as ¹³CH₂D₂ or ¹²CHD₃ are ignored due to their low abundances). For the subscripts in Fig. 3a (m, i, and j in r_amij or r_bmij), the first digit (m = 0 or 1) is the number of ¹³C atoms involved in the reaction, the second digit (i = 0, 1, or 2) is the number of deuterium atoms connected in the methylene or methyl group, and the third digit (j = 0 or 1) is the number of deuterium atoms in the hydrogen donor.

Clumped isotopic compositions of methylene and methane are defined as the following:

Note that the isotopic compositions here are expressed in decimals; they should be multiplied by 1000 to give per mil values.

The deuterium isotope ratio between the hydrogen donor (denoted with subscript B) and the methylene group (subscript A) is expressed as:

For each reaction step in Fig. 3a, the corresponding rate constants are denoted as k_amij for step a or k_bmij for step b. Kinetic fractionation factors α_kamij = k_amij/k_a000 and α_kbmij = k_bmij/k_b000 define KIEs. Note that a DKIE is often expressed as k_H/k_D, which is the reciprocal of the α_k nomenclature here. A DKIE may be primary or secondary; a primary DKIE results in α_ka001 ≠ 1 and α_kb001 ≠ 1, and a secondary one results in α_ka010 ≠ 1 and α_kb010 ≠ 1. Kinetic clumped isotope fractionation factors γ_amij = α_kamij/(α_ka100^mα_ka010ⁱα_ka001^j) and γ_bmij = α_kbmij/ (α_kb100^mα_kb010ⁱα_kb001^j) define the excessive KIE due to isotope clumping in steps a and b, respectively³⁰.

Conversion of the reactant R-CH₂-R′ is defined as 1 − f, where f is the residual fraction of R-CH₂-R′:

Considering the isotope abundance of D << H, the analytical solution of isotopic compositions at the beginning of the reaction (f = 1) is derived as:

If the abundance of the hydrogen donor is excessive and thus approximately constant, then the analytical solution at the end of the reaction (f = 0) is

If the KIE is absent (all α_k and γ factors equal to unity) and the clumped isotopic compositions in the precursor are 0, then clumped isotopic compositions become the following:

With ({alpha }_{{{{{{rm{A}}}}}}}^{{{{{{rm{B}}}}}}}) = 0.354 from δD_A = −126‰ and δD_B = −691‰ given in the text, Δ¹²CH₂D₂ = −76‰ is obtained. Both Δ¹³CH₃D and Δ¹²CH₂D₂ are more depleted of clumped isotopes than the reported values (4–6‰ for Δ¹³CH₃D and −10 to 5‰ for Δ¹²CH₂D₂)^10,11, indicating that ¹³C–D clumping in the methylene precursor should be considered to explain Δ¹³CH₃D, and DKIE should be considered to explain Δ¹²CH₂D₂.

Numerical simulations are carried out to find parameters satisfying the following:

1. (1)

   The δ¹³C_CH4, δD_CH4, Δ¹³CH₃D, and Δ¹²CH₂D₂ values at higher conversions of organic precursors are within the reported ranges.

2. (2)

   The δD_A and δD_B values are close to the derived values from Fig. 1c (−691 and −126‰, respectively) to show that δD_CH4 is mainly determined by δD of the precursors rather than by DKIE during the hydrogenolysis.

A value of ΔR¹³CHDR′ =6‰ in the organic precursor is applied so that the final Δ¹³CH₃D is in the range of reported values. This ¹³C–D clumping in the precursor is acceptable, considering that Δ¹³CH₃D = 5.6‰ has been reported for biogenic gas^10,11. A Δ¹²CH₂D₂ value close to the observed value but much higher than the stochastic one (Eq. 14) requires γ_a011 > 1 or γ_b011 > 1, as shown by the Δ¹²CH₂D₂ expression in Eq. (13). With this prerequisite, either an inverse primary DKIE (1° DKIE, α_ka001 > 1, α_kb001 > 1) or an inverse secondary DKIE (2° DKIE, α_ka010 > 1, α_kb010 > 1) is necessary, and through numerical simulation, we found that only the inverse 1° DKIE satisfies the above-mentioned δD_A, δD_B, and Δ¹²CH₂D₂ values.

Two scenarios (one is the pure stochastic condition, the other is with an inverse 1° DKIE) are modelled (Fig. 3). The parameters are listed in Table 1. For comparison, analytical solutions at the beginning and end of reactions from Eqs. (12) and (13) are presented. The numerical and analytical solutions are nearly identical at the beginning of conversion. There are small differences between the numerical and analytical solutions at the endpoint because the abundance of the hydrogen donor is not extremely excessive. A weak ¹³C fractionation between the organic precursor and the methane product is obtained with the KIE parameters (Fig. 3b). With such a weak ¹³C KIE, Δ¹³CH₃D is nearly constant for reaction extent (Fig. 3c). Note that we applied an inverse ¹³C KIE, as required by the δ¹³C distribution of the alkane gases (Method 1, Model A). The δD and Δ¹²CH₂D₂ values are independent of ¹³C KIE. Both the bulk and clumped isotopic compositions of methane within the range of reported values are obtained at the organic precursor conversion of 0.65–0.70 as constrained by Fig. 2.

Source: Ecology - nature.com