# Reversibility controls on extreme methane clumped isotope signatures from anaerobic oxidation of methane

Jiarui Liu, Rachel L. Harris, Jeanine L. Ash, James G. Ferry, Sebastian J. E.

Krause, Jabrane Labidi, Divya Prakash, Barbara Sherwood Lollar, Tina

Treude, Oliver Warr, et al.

<figure><img src="image_1.png" /></figure>

► To cite this version:

Jiarui Liu, Rachel L. Harris, Jeanine L. Ash, James G. Ferry, Sebastian J. E. Krause, et al. Reversibility controls on extreme methane clumped isotope signatures from anaerobic oxidation of methane. Geochimica et Cosmochimica Acta, 2023, 348, pp.165-186. 10.1016/j.gca.2023.02.022. insu-04155736

HAL Id: insu-04155736

https://insu.hal.science/insu-04155736

Submitted on 8 Jul 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# Reversibility controls on extreme methane clumped isotope signatures from anaerobic oxidation of methane

<figure><img src="image_4.png" /></figure>

Jiarui Liu<sup>a,*</sup>, Rachel L. Harris<sup>b</sup>, Jeanine L. Ash<sup>c</sup>, James G. Ferry<sup>d</sup>, Sebastian J.E. Krause<sup>a</sup>, Jabrane Labidi<sup>e</sup>, Divya Prakash<sup>f</sup>, Barbara Sherwood Lollar<sup>e,g</sup>, Tina Treude<sup>a,h</sup>, Oliver Warr<sup>g</sup>, Edward D. Young<sup>a</sup>

## ARTICLE INFO

**Article history:**
Received 17 August 2022
Accepted 27 February 2023
Available online 9 March 2023
Associate editor: Filip Meysman

**Keywords:**
Methane isotopologues
Methyl-coenzyme M reductase
Kinetic fractionation
Equilibrium fractionation
Isotopic bond re-ordering

## ABSTRACT

Microbial anaerobic oxidation of methane (AOM) substantially mitigates atmospheric methane emissions on Earth and is a process to consider for astrobiological targets where methane has been detected. The measurement of doubly substituted, or “clumped”, methane isotopes has proven useful in tracing processes of methane formation and oxidation. Both near-equilibrium and extreme disequilibrium methane clumped isotope signatures can be attributed to AOM, but, to date, understanding the mechanistic and environmental controls on those signatures has been lacking. We report measurements of methane clumped isotope compositions of residual methane in AOM-active microbial incubations using sediment slurries from Svalbard and Santa Barbara Channel methane seeps. Incubation experiments of Svalbard sediment slurries resulted in residual methane with very high $\Delta^{13}CH_3D$ and $\Delta^{12}CH_2D_2$ values up to 19.5‰ and 65.1‰, respectively. We found similarly high $\Delta^{13}CH_3D$ and $\Delta^{12}CH_2D_2$ values in fluid samples from the Chamorro Seamount, a serpentinite mud volcano in the Mariana forearc, suggesting that minimal reversibility of AOM intracellular reactions leads to kinetic fractionation of clumped isotopologues. When conditions were consistent with a low thermodynamic drive for AOM, however, methane isotopologues approached intra-species quasi-equilibrium. This was clearly observed in isotope exchange experiments with methyl-coenzyme M reductase (Mcr) and in microbial incubations of the Santa Barbara Channel sediment slurries. Using an isotopologue fractionation model, we highlight the critical role of reversibility in controlling the trajectory of gases in $\Delta^{13}CH_3D$ vs. $\Delta^{12}CH_2D_2$ space during AOM. The near-equilibrium methane isotopologue signatures are generalized as a result of the Mcr-catalyzed intracellular isotope exchange operating under near-threshold free energy conditions, as shown in the deep-biosphere incubations. Our results show that the reversibility of the Mcr-catalyzed reaction is central to understanding the meaning of methane isotopologue ratios affected by microbial production and oxidation.

© 2023 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

Methane, the simplest of all hydrocarbons, is a gas with economic, societal, and environmental relevance. It is the main component of natural gas, a key component of the global carbon cycle (Reeburgh, 2007), and a primary target for understanding the rela-

tive role of geologic processes and potential extinct/extant life elsewhere in the universe (e.g., Yung et al., 2018). As a potent greenhouse gas, methane contributes to climate change and its concentration in the atmosphere has been increasing since the beginning of industrialization (IPCC, 2021). The main “biological” processes forming methane are microbial degradation (microbial-genic) and pyrolysis (thermogenic) of organic matter (e.g., Reeburgh, 2007; Repeta et al., 2016). Abiotic methane can be formed through various types of water-rock interactions (Etiope and

<footnote><sup>a</sup>Department of Earth, Planetary and Space Sciences, University of California, Los Angeles, CA, USA</footnote>

<footnote><sup>b</sup>Department of Organismic and Evolutionary Biology, Harvard University, Cambridge, MA, USA</footnote>

<footnote><sup>c</sup>Department of Earth, Environmental and Planetary Sciences, Rice University, Houston, TX, USA</footnote>

<footnote><sup>d</sup>Department of Biochemistry and Molecular Biology, Pennsylvania State University, University Park, PA, USA</footnote>

<footnote><sup>e</sup>Université de Paris, Institut de Physique du Globe de Paris, CNRS, Paris, France</footnote>

<footnote><sup>f</sup>School of Chemical and Biomolecular Sciences, Southern Illinois University, Carbondale, IL, USA</footnote>

<footnote><sup>g</sup>Department of Earth Sciences, University of Toronto, Toronto, Ontario, Canada</footnote>

<footnote><sup>h</sup>Department of Atmospheric and Oceanic Sciences, University of California, Los Angeles, CA, USA</footnote>

<footnote>* Corresponding author.
E-mail address: jiaruiliu@ucla.edu (J. Liu).</footnote>
Sherwood Lollar, 2013 and references therein) under a range of both low-temperature and high-temperature conditions. Identifying the sources and processes of methane emission is vital for constraining fluxes from major reservoirs contributing to climate change, and more broadly, for understanding the biogeochemical dynamics of the planet's critical zone. Traditionally, microbial and thermogenic methane gases are interpreted by their carbon (δ¹³C) and hydrogen (δD) isotopic compositions (e.g., Schoell, 1983; Whiticar, 1999). However, applying δ¹³C and δD as a diagnostic framework alone has significant limitations, most notably circumstances in which isotopic signatures of methane sources are not sharply demarcated, exhibit some degree of overlap, and such empirical frameworks are inevitably redefined as additional experimental and field data come to light (Sherwood Lollar et al., 2006; Douglas et al., 2017). In addition, empirical fields in δ¹³C vs. δD space are inherently linked to the isotopic signature of the source materials coupled with the formation processes, which in turn are strongly affected by mixing and secondary alteration processes (Etiope and Sherwood Lollar, 2013 and references therein).

Recently, sophisticated instrumentation has become available to resolve the abundances of doubly substituted, "clumped", isotopologues of methane gas, including ¹³CH₃D and ¹²CH₂D₂ (Ono et al., 2014; Stolper et al., 2014b; Young et al., 2016; Eldridge et al., 2019; Gonzalez et al., 2019; Dong et al., 2021; Zhang et al., 2021). The relative abundances of these two mass-18 methane isotopologues are expressed as Δ¹³CH₃D and Δ¹²CH₂D₂ values, representing deviations from a stochastic standard in which distributions of isotopes across all isotopologues are effectively random. The latter is the expected state at sufficiently high temperatures (>1000 K). The relative abundances of doubly substituted methane isotopologues differ from the stochastic distribution at lower temperatures due to the enhanced thermodynamic stability of bonds involving the heavier isotopes (e.g., Ma et al., 2008; Webb and Miller, 2014; Liu and Liu, 2016). In principle, isotopic bond ordering, or clumping, is a tracer of the processes attending the formation of the molecules, coupled with any consumption and/or re-equilibration processes, and therefore gives the potential for their application to discern processes of methane formation and oxidation both on Earth and beyond (Young, 2019). The utility of Δ¹³CH₃D and Δ¹²CH₂D₂ lies in the fact that the values should be sensitive to process, and interpretations of their significance may not dependent upon long-lost isotope exchange partners (e.g., water). In practice, kinetically controlled methane clumped isotope compositions can be sensitive to hydrogen sources of methane that differ in δD (Taenzer et al., 2020).

In some cases, Δ¹³CH₃D has been used to reconstruct the formation temperatures of thermogenic methane (Stolper et al., 2014a; Stolper et al., 2015; Wang et al., 2015; Douglas et al., 2016; Lalk et al., 2022) and the combination of Δ¹³CH₃D and Δ¹²CH₂D₂ has been used to verify that the molecules record equilibrium at a given temperature (Young et al., 2017; Giunta et al., 2019; Gonzalez et al., 2019; Thiagarajan et al., 2020; Xie et al., 2021). The abundances of mass-18 isotopologues of microbial methane, on the other hand, have been shown to be inconsistent with thermodynamic equilibrium, meaning that their distribution does not reflect the temperature at which the methane was formed (Stolper et al., 2015; Wang et al., 2015; Young et al., 2017; Ash et al., 2019; Giunta et al., 2019; Warr et al., 2021). The Δ¹³CH₃D and Δ¹²CH₂D₂ of microbial methane can be modeled as the result of the combination of equilibrium and kinetic steps during the enzymatically facilitated formation of methane (Cao et al., 2019; Young, 2019; Gropp et al., 2021; Gropp et al., 2022; Ono et al., 2022). In particular, the extremely negative Δ¹²CH₂D₂ values are most likely due to the combinatorial effect of accessing multiple

sources of hydrogen with significantly different hydrogen isotopic signatures during the formation of methane (Röckmann et al., 2016; Yeung, 2016; Young, 2019; Taenzer et al., 2020).

While the ability to define a methane isotopologue signature related to microbial methanogenesis is promising, there are caveats. One is that abiotic methane and low-maturity thermogenic methane has been found with deficits in Δ¹²CH₂D₂ compared with equilibrium values in methane gas, which is further confirmed by pyrolysis experiments on n-octadecane (Douglas et al., 2017; Young et al., 2017; Dong et al., 2021; Warr et al., 2021; Xie et al., 2021). On the other hand, abiotic or biotic consumption processes could potentially alter the relative abundances of the mass-18 isotopologues in ways that might mimic the putative microbial methanogenesis signature. Laboratory experiments and *ab initio* modeling indicate that methane oxidation in the atmosphere by OH or Cl radicals is accompanied by large shifts in methane clumping down to very negative Δ¹³CH₃D and Δ¹²CH₂D₂ values due to Rayleigh distillation (Haghnegahdar et al., 2017; Whitehill et al., 2017; Young, 2019). Similarly, both Δ¹³CH₃D and Δ¹²CH₂D₂ values of residual methane in microbial culture were observed to decrease during aerobic oxidation of methane (Wang et al., 2016; Krause et al., 2022).

Photocatalytic oxidation and aerobic microbial oxidation of methane occur in the atmosphere and oxygenated environments. Therefore, these processes are in principle separated spatially from microbial methanogenesis, which operates almost exclusively under anoxic, and thus reducing, conditions. Anaerobic oxidation of methane (AOM), however, is closely related to methanogenesis both spatially and metabolically. For example, concurrent methanogenesis and AOM have been found in near-surface sediments and at the sulfate-methane transition (Xiao et al., 2017; Beulig et al., 2019; Krause and Treude, 2021). Anaerobic methanotrophic (ANME) archaea are close relatives of methanogenic archaea and have the full enzymatic machinery of methanogenesis working in reverse (Timmers et al., 2017). Methanogens belonging to *Methanosarcina* can reverse key reactions of methanogenesis coupled to ferric iron reduction (Yan et al., 2018; Yu et al., 2022). The prospect of a methanotroph (ANME-1) from natural sediments conducting methanogenesis further underscores the potential versatility of methanoarchaeal enzymes (Beulig et al., 2019; Kevorkian et al., 2021). Factors forcing these enzymes to operate in either direction need further investigation. However, even when operating in full AOM mode with the highest net turnover of methane, ANME archaea demonstrate some fraction (~4–15%) of net reversibility, defined as the ratio of the backward to forward fluxes (Treude et al., 2007; Holler et al., 2011; Wegener et al., 2021). Back flux of material during the enzymatic reactions of AOM has implications for isotope fractionation, leading to carbon isotope (¹³C/¹²C) equilibration between methane and intracellular metabolites, in particular when AOM is operating close to thermodynamic limitations where the methyl-coenzyme M reductase (Mcr)-catalyzed reaction is near equilibrium (Hoehler et al., 2000; Holler et al., 2011; Yoshinaga et al., 2014; Chuang et al., 2019; Wegener et al., 2021).

The precise mechanisms that are responsible for Δ¹³CH₃D vs. Δ¹²CH₂D₂ effects of AOM are largely unknown, and further laboratory incubation experiments under different environmental conditions are required. The overarching goal is to develop these intra-methane isotopologue signals as tracers for the origin and processing of methane in general, and AOM in particular. One of the biggest challenges for understanding the isotopologue composition of natural samples is limited knowledge of how isotopic bond re-ordering caused by AOM can overwrite formation-process signatures. Given the high level of reversibility of
Mcr-catalyzed reactions during AOM, Ash et al. (2019), Giunta et al. (2019), and Warr et al. (2021) suggested that AOM could be responsible for thermodynamic isotopic bond-order equilibrium in methane found in sulfate-depleted Baltic Sea sediments, sedimentary basins, and the deep subsurface, respectively. Recently, Ono et al. (2021) found that the methane clumped isotopologue ¹³CH₃D exhibits kinetic fractionation in an ANME-1 sediment-free enrichment culture with a high sulfate concentration (28 mM), where Δ¹³CH₃D values increased by up to 5.8%. Importantly, Wegener et al. (2021) demonstrated a sulfate dependence of the net bulk isotope fractionation by virtue of the thermodynamic drive of the intracellular reactions of AOM. Based on the earlier studies, therefore, it appears that the reversibility of intracellular reactions of AOM is responsible for the observed discrepancy between the kinetic isotopologue fractionation in a sulfate-replete culture (Ono et al., 2021) and the intra-species isotope equilibrium in sulfate-depleted settings (Ash et al., 2019; Giunta et al., 2019; Warr et al., 2021).

To document the range of effects that AOM may have on modifying methane isotopologue signatures, we collected a suite of marine sediments and fracture fluids from diverse marine and terrestrial settings that showed different degrees of ANME abundance and AOM activity (Fig. 1). We performed microbial incubations using these samples to evaluate the relationship between the reversibility of intracellular reactions of AOM and the methane isotopologue fractionation/re-ordering, and to elucidate the unknown fractionation factors for ¹²CH₂D₂/¹²CH₄. We also collected natural fluids from the South Chamorro Seamount where AOM is known to substantially consume abiotic source methane (Wheat et al., 2020), in order to examine whether the isotopologue fractionation differs between laboratory incubations and endogenous AOM in natural environments. In combination, our work aims to address how environmental conditions can potentially affect isotopologue fractionation, as well as rates and magnitudes of isotopic bond re-ordering during AOM. This information will improve our ability to distinguish primary signatures of methane formation from biologically-driven overprinting to maximize the utility of ¹³CH₃-D/¹²CH₄ and ¹²CH₂D₂/¹²CH₄ as tracers of methane formation and subsequent processing.

## 2. Materials and methods

### 2.1. Incubation of methane seep sediment slurry from Svalbard

Marine sediments were collected in summer 2011 with a video-guided multicorer from active methane seeps covered by a dense sulfur-bacteria mat offshore northwestern Svalbard during the R/V Poseidon cruise 419 (Site MUC 12; 79°00.417′ N, 06°54.131′ E, 1235 m water depth). Details of the sampling site, including ex-situ activity of AOM (Fig. 2A) and sulfate reduction and geochemical parameters in replicate sediment cores, can be found in Melaniuk et al. (2022). Sediment used for the present study was collected from the top 10 cm of a replicate multicorer core (10 cm inner diameter) and filled headspace-free into 100- and 250-ml borosilicate bottles with butyl stoppers. The sediment was stored anoxically at 4 °C in the dark for 9 years without methane or sulfate addition. The AOM community in the sediment was reactivated in the laboratory at UCLA in April 2020 by mixing with anoxic artificial seawater medium at a ratio of 1:1 (v:v) (Widdel and Bak, 1992; Laso-Pérez et al., 2018) and a methane gas headspace as the sole amended carbon source in a 2-L glass bottle. The cultivation procedures of Laso-Pérez et al. (2018) were followed.

We performed a set of sulfide analyses to confirm the activity of AOM before starting the actual experiments. Sulfide concentration was measured after Cord-Ruwisch (1985) using a Shimadzu UV-Spectrophotometer (UV-1800). The relative standard deviation (RSD) of sulfide analysis was better than 5%. Sulfide concentration in the supernatant started to increase two months after the addition of substrates (sulfate and methane). Further, sulfide concentration increased ca. 0.07 mM per day between the 134th and 205th days since reactivation. The supernatant of the sediment slurry was replaced with fresh anoxic medium with seawater sulfate concentration (28 mM) when sulfide concentration approached 11 mM. After that, sulfide concentration continued to increase ca. 0.07 mM per day between the 206th and 283rd days since reactivation. Sulfide production from sulfate reduction with non-methane substrates was excluded by the control experiment (see below). The linear increase in sulfide concentration over time

<figure><img src="image_5.png" /><figcaption>**Fig. 1.** Global map showing sampling sites for this study. Marine sediments or fracture fluids were collected for slurry or fluid incubations at (1) Svalbard methane seeps, (2) Santa Barbara Channel methane seeps, (3) Nankai Trough (International Ocean Discovery Program Hole C0023A), and (4) Beatrix Gold Mine, South Africa. Methane in natural fluids was collected from a sub-seafloor borehole observatory (Ocean Drilling Program Hole 1200C) on (5) South Chamorro Seamount, a serpentinite mud volcano in the Mariana forearc.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 2.** Geochemistry of the Svalbard methane seep sediment and slurry. (A) Ex-situ AOM rates determined by ¹⁴CH₄ radiotracer at Site MUC12 offshore Svalbard (Melaniuk et al., 2022). R1–3 represents replicate one through three. (B) Evolution of dissolved sulfide concentration and percentage of methane oxidized in the Svalbard methane seep sediment slurry with a starting sulfate concentration of 20 mM. Linear regressions are shown in panel B.</figcaption></figure>

indicates that the ANME community was successfully reactivated in the sediment slurry. Thus, the homogenized slurry was equally distributed to two 600-ml culture bottles in replicate and filled up completely with an anoxic artificial seawater medium (Laso-Pérez et al., 2018). Sixty ml of the medium volume was replaced with ca. 200 kPa of methane (Airgas), and both incubations were performed under excess sulfate (ca. 10 mmol) over methane (5.4 mmol) (cf. Ono et al., 2021). The sediment slurry was then stored horizontally at 4 °C in the dark. Methane isotopologue compositions were analyzed on days 1, 36, 44, 58, and 71 of the incubation. After the first experiment, we reset the slurry for longer incubation through medium replacement and bubbling with N₂/CO₂ gas. Methane isotopologue compositions were analyzed on days 1, 60, 90, and 168 of the incubation.

To further confirm and quantify the activity of AOM and sulfate reduction, we measured turnover rates with ¹⁴C- and ³⁵S-radiotracer techniques, respectively, in the high-sulfate experiment. AOM rate was determined by injecting 20 µL of ¹⁴C-methane tracer (dissolved in MilliQ water, activity 0.7 kBq, specific activity 185 MBq mmol⁻¹) into two 13 ml headspace-free sub-samples of the sediment slurry. Sulfate reduction rate was determined by injecting 10 µL of carrier-free ³⁵S-sulfate tracer (dissolved in MilliQ water, activity 234 kBq, specific activity 37 TBq mmol⁻¹) into a 10 ml sub-sample of the sediment slurry. The vials were incubated for two days at 4 °C in the dark. After incubation, sulfate reduction activity was terminated by transferring the sample from the glass vial to a centrifuge tube filled with 20 ml 20% zinc acetate followed by freezing at -30 °C. Sulfate reduction rate was analyzed and calculated according to the cold chromium distillation method (Kallmeyer et al., 2004). To terminate AOM activity, samples were transferred to 50 ml glass vials filled with 20 ml 5% sodium hydroxide. The vials were sealed with rubber stoppers immediately and shaken thoroughly. AOM rate was determined according to Treude et al. (2005) and Joye et al. (2004).

After the high-sulfate experiment described above, the supernatant in the slurry was replaced with low-sulfate artificial seawater medium (1 mM) four times. The sulfate concentration in the slurry decreased over time due to medium replacement and microbial consumption. The produced sulfide was maintained at low concentrations (<1 mM) through medium replacement and bubbling with N₂/CO₂ gas. The culture bottle was then filled up completely with the low-sulfate medium, reaching a sulfate concentration of 0.4 mM in the slurry. Sixty ml of the medium volume

was replaced with ca. 200 kPa of methane (Airgas), allowing the incubation to be performed under excess methane (5.4 mmol) over sulfate (<0.2 mmol). The low-sulfate experiment was performed identically to the high-sulfate experiment with sulfate concentration as the only variable at the beginning of each incubation. Towards the end of the low-sulfate experiment, sulfate was depleted (0.01 mM) in the slurry. Thus, an additional concentrated anoxic sulfate solution (2.5 ml) was added to the culture bottle by a syringe without opening the bottle, reaching a final sulfate concentration of 0.4 mM again in the slurry. Sulfate concentration was analyzed by ion chromatography (Metrohm 761). Alkalinity was determined by acid titration with a Metrohm 876 Dosimat Plus (Dale et al., 2015). The concentrations were calibrated against the IAPSO standard seawater (n = 10) with RSD of sulfate and alkalinity analyses better than 2%. In addition, pH was determined by a pH meter (VWR sympHony B10P).

Methane concentration and isotopologue abundance were determined during the incubations on a vacuum line interfaced with a gas chromatograph, and the Panorama mass spectrometer, respectively, as described below in section 2.6. Given the known volume and porosity of the slurry, the amounts of methane that dissolved in the liquid phase were calculated using Henry's law and the Bunsen solubility coefficient (Yamamoto et al., 1976). For porosity analysis, the homogenized slurry was sampled using a syringe. Porosity was measured from the volume of slurry and water weight loss after drying at 60 °C until a constant weight was reached. The total amounts of methane in the culture bottle were therefore the sum of methane in the headspace and the liquid phase.

The liquid phase was also sampled to determine the δ¹³C of dissolved inorganic carbon (DIC) and the δD of water in the slurry (Atekwana and Krishnamurthy, 1998; Kopec et al., 2019). In brief, analysis of δ¹³C_DIC was performed in 12 ml Labco Exetainer vials on a Thermo Scientific GasBench II coupled to a Thermo Finnigan Delta Plus XL isotope-ratio mass spectrometer (IRMS) at the UC Davis Stable Isotope Facility. Isotopic values are reported in δ-notation relative to the Vienna Pee Dee Belemnite (VPDB) standard. Analytical precision for δ¹³C is better than 0.1‰ (1σ) based on replicate analyses of laboratory standards. Analysis of hydrogen isotopic composition of water was conducted at the Stable Isotope Laboratory at Dartmouth College. The δD values of water were measured using an H-Device, in which water was reduced by hot chromium (850 °C), and the resulting hydrogen gas was measured by a Thermo Delta Plus XL IRMS. Isotopic ratios are reported in δ-
notation relative to the Standard Mean Ocean Water (VSMOW) standard. Analytical precision for δD is better than 0.5‰ (1σ) based on replicate analyses of laboratory standards.

A control incubation was set up to assess potential microbial methanogenesis and organoclastic sulfate reduction in the same sediment slurry. The control experiment was prepared by replacing the supernatant with a new methane-free anoxic artificial seawater medium and a 60-ml headspace was left to monitor methane concentration. The slurry was then bubbled thoroughly with N₂/CO₂ (80:20) three times (3 h in total) to get rid of the residual dissolved methane and sulfide from the previous experiments. The slurry was stirred during and between flushing to drive the dissolved methane into the headspace. The control experiment slurry was stored in the dark at 4 °C for three months. The slurry supernatant was collected for sulfide and sulfate concentrations at the start and end of the control incubation. Methane concentration in the gas phase was measured by gas chromatography throughout the three months.

## 2.2. Incubation of methane seep sediment slurries from the Santa Barbara Channel

Marine sediments were collected with four push cores from the Coal Oil Point seep field (Rostocker Seep and Isla Vista Super Seep) in the Santa Barbara Channel in 2017 (Jordan et al., 2020). The top 12.5 cm sediment of the cores was transferred to 300 ml culture bottles at UCLA following the same procedure as detailed in the Svalbard slurry section above. Low sulfate (<1 mM) and high sulfate (>28 mM) concentrations were achieved in the sediment slurry by either dilution steps with anoxic artificial seawater (without sulfate) or additions of concentrated sodium sulfate solution (100 mM). The sediment slurry was incubated in the dark with a 4.4 ml headspace of methane at ca. 100 kPa and 20 °C on a shaker. The methane gas bubble from the sediment slurry was extracted at the end of each experiment and geochemistry analyses were performed as described in the Svalbard slurry section. A similar control experiment followed the Svalbard slurry section with a 4.4 ml Argon headspace.

## 2.3. Incubations of sediment slurry and fracture fluid from the deep biosphere (Nankai Trough and Beatrix Gold Mine)

Hole C0023A (32°22.0018' N, 134°57.9844' E), located in the Nankai Trough on the subduction boundary between the Philippine Sea and Eurasian plates, was drilled during International Ocean Discovery Program (IODP) Expedition 370 in 2016 (Heuer et al., 2017). Marine sediments from 257 m below seafloor (mbsf) were sampled anoxically with rigorous contamination control (see Heuer et al., 2017) and utilized for sediment incubation. In brief, whole-round cores (WRCs) were prepared from recovered sections under super-clean and anoxic conditions onboard the Chikyu Drilling Vessel. X-ray computed tomography was performed on all WRCs for lithological identification and to assess core quality on board. Samples without drilling disturbances were designated for incubation and were subjected to secondary scraping to minimize the potential for contamination.

At South Africa's Beatrix Gold Mine (28 14.1' S, 2647.7' E), fracture fluid was collected in 2016 from the BE326-BH2 borehole, from 1390 m below land surface (mbls) in shaft 3, level 26. A detailed description of the study site and sample collection has been published in Lau et al. (2016). Briefly, fracture fluid was filtered using a 0.2 µm hollow fiber MediaKap®-10 filter. The microorganisms on the filter were anoxically back-flushed into sterile, capped, N₂-sparged 160-ml borosilicate serum vials to a final cell concentration of ~10⁷ cells ml⁻¹. Samples were stored at 4 °C upon return to the surface and used for incubation.

Incubations were set up inside a clean anoxic glove bag at Princeton University. Sterile aluminum foil was placed on the working surface of the glove bag and sterile nitrile gloves were used over the glove bag's butyl rubber gloves to minimize potential contamination of low-biomass samples. The surface of C0023A samples was scraped using a sterile scalpel to remove sediment exposed to gas-tight packaging. A total of 10 g (wet weight) of interior sediment was weighed out from sample source core C0023A-5F02 (top depth 257.3 mbsf) and transferred into combusted 160-ml borosilicate serum vials containing 100 ml artificial sulfate-free seawater medium modified from the recipe by Widdel and Bak (1992). For the BE326-BH2 incubation, 10 ml of fracture fluid was added to 90 ml of artificial sulfate-free seawater medium. Serum vials were sealed with 0.1 N NaOH-boiled butyl rubber stoppers and aluminum crimps, and the headspace was subsequently replaced with 100% methane. Each incubation was supplemented with 10 mM of one of the following electron acceptors: NO₃⁻, NO₂⁻, SO₄²⁻, or Fe³⁺ (in the form of hydrous ferric oxide). A control without added electron acceptors was included to assess endogenous microbial activity from remnant electron acceptors available in the inoculum, and an autoclaved sample was included as a killed control. Serum vials containing sediment slurry or fracture fluid were incubated upside-down at an approximate in-situ temperature of 40 °C. At UCLA, the collected residual methane was purified for isotopologue analysis following the method described in section 2.6.

Electron acceptor depletion was monitored using Dionex IC25 ion chromatography coupled to an MSQ-quadruple mass spectrometer (Thermo Scientific). The sediment slurry or fracture fluid was also subsampled to monitor the δ¹³C of DIC. Briefly, 500 µL aliquots were anoxically transferred to combusted, amber borosilicate serum vials treated with saturated HgCl₂ and aluminum-crimped sealed with butyl rubber stoppers. Following sample transfer, vials were over-pressurized with ultra-high purity N₂ gas, supplemented with 0.5 N H₃PO₄, and heated overnight in a water bath at 70 °C to extract all DIC out of the solution. Isotopic composition of the headspace was analyzed using a Picarro cavity ring-down spectrometer equipped with a G2101-I Isotopic CO₂ analyzer in CO₂ focus mode.

## 2.4. Methyl-coenzyme M reductase (Mcr) experiments

The experiment was performed as described previously with the Mcr enzyme purified from Methanothermobacter marburgensis (Mahlert et al., 2002; Scheller et al., 2010) at Pennsylvania State University. Mcr that was used for this assay had absorbance maxima at 387 nm, corresponding to the enzyme's active (Ni⁺) form. Coenzyme B (HS-CoB) was prepared from the symmetric disulfide CoB-S-S-CoB by reaction with NaBH₄. Methyl-coenzyme M (CH₃-S-CoM) was synthesized from coenzyme M (sodium salt) by methylation with methyl iodide. The complete reaction mixture (4 ml) contained 5 mM CH₃-S-CoM, 2.5 mM HS-CoB and 32.3 mg of Mcr in 50 mM phosphate buffer (pH 7.6). This allowed isotope exchange of the methyl moiety in methyl-coenzyme M and methane to occur. The reaction mixtures were incubated at 20 °C for 48 h or 60 °C for 2–6 h in a stoppered 10-ml serum vial with a methane headspace of ca. 100 kPa. The difference in incubation times is due to more rapid reaction at 60 °C. Control reactions contained the complete reaction mixture minus enzyme. At UCLA, the collected gas was purified for isotopologue analysis following the method detailed in section 2.6.

## 2.5. Methane-bearing natural fluids at the South Chamorro Seamount

South Chamorro Seamount is an active serpentinite mud volcano at ca. 3150 m water depth in the Mariana forearc. In 2001,
Hole 1200C was drilled to 266 mbsf and cased during Ocean Drilling Program (ODP) Leg 195 in the summit knoll of South Chamorro Seamount (Fryer and Salisbury, 2006). An oceanic borehole observatory, commonly called a CORK (Circulation Obviation Retrofit Kit), was deployed for subsequent fluid sampling (Wheat et al., 2008). The discharged fluids were largely altered relative to seawater composition, displaying a high pH up to 12.3 with abundant dissolved methane up to 33 mM (Wheat et al., 2008; Wheat et al., 2020). In January 2009, pristine crustal fluids were collected directly as they discharged at Hole 1200C with a remotely operated vehicle (ROV) HyperDolphin (HPD Dives 941–947) during Cruise NT09-01 (Wheat et al., 2020). Isobaric gas-tight samplers (Seewald et al., 2002) were used for fluid sampling. Immediately upon recovery of the ROV, fluid samples were subsampled for liquid and gas analyses. Gas extraction from fluid samples was conducted using previously reported techniques (Seewald et al., 2002). At UCLA, the collected gas was purified for isotopologue analysis following the method detailed in the next section.

## 2.6. Doubly substituted isotopologue measurements and isotope notation

Methane isotopologue abundances of methane gas samples were measured using the Panorama (Nu Instruments) high-mass-resolution gas-source isotope ratio mass spectrometer housed at UCLA. Details surrounding the purification and measurement of methane gas were previously published (Young et al., 2016; Young et al., 2017) and are briefly summarized here. Methane sample gases were purified on a vacuum line interfaced with a gas chromatograph (GC). Samples were delivered to the vacuum line through a septum by a gas-tight syringe and trapped on silica gel at liquid nitrogen temperature. The Helium carrier gas was then used to flush the sample to the GC. Separation was accomplished with a 3-meter 1/8-inch OD stainless steel column packed with 5 Å molecular sieve, followed in series by a 2-meter 1/8-inch OD stainless steel column packed with HayeSep D porous polymer. Peaks were identified using an in-line, passive thermal conductivity detector (TCD). Once methane collection was complete, the sample was transferred to an evacuated sample tube filled with silica gel at liquid nitrogen temperature. Methane in this tube was introduced to the inlet of the mass spectrometer where it was warmed to 40 °C and expanded into the bellow of the instrument.

The Panorama mass spectrometer was set to a mass resolving power of ~40,000 or greater, allowing the measurement of ion currents for resolved ¹²CH₄, ¹³CH₄, ¹²CH₃D⁺, ¹³CH₃D⁺, and ¹²CH₂D₂⁺. Isotopologues of masses 16 and 17 were measured using Faraday collectors with amplifier resistors of 10¹¹ Ω. Both doubly substituted mass-18 isotopologues, ¹³CH₃D⁺ and ¹²CH₂D₂⁺, were measured with an electron multiplier as the axial collector. The measured ratios of these ion currents yield values for bulk ¹³C/¹²C and D/H as well as for both Δ¹³CH₃D and Δ¹²CH₂D₂. The isotopic compositions of carbon and hydrogen are reported as deviations from the carbon and hydrogen reference materials VPDB and VSMOW. Standard delta notation is used to express the fractional differences in per mil units:

$$ \delta^{13}\text{C} = [(^{13}\text{C}/^{12}\text{C})_{\text{sample}}/(^{13}\text{C}/^{12}\text{C})_{\text{VPDB}} - 1] \times 1000 \quad (1) $$

$$ \delta\text{D} = [(\text{D/H})_{\text{sample}}/(\text{D/H})_{\text{VSMOW}} - 1] \times 1000 \quad (2) $$

The relative abundances of the two mass-18 isotopologues of methane are reported relative to the stochastic reference frame expressed in per mil using the capital delta notation:

$$ \Delta^{13}\text{CH}_3\text{D} = [(^{13}\text{CH}_3\text{D}/^{12}\text{CH}_4)_{\text{sample}}/ \times (^{13}\text{CH}_3\text{D}/^{12}\text{CH}_4)_{\text{stochastic}} - 1] \times 1000 \quad (3) $$

$$ \Delta^{12}\text{CH}_2\text{D}_2 = [(^{12}\text{CH}_2\text{D}_2/^{12}\text{CH}_4)_{\text{sample}}/ \times (^{12}\text{CH}_2\text{D}_2/^{12}\text{CH}_4)_{\text{stochastic}} - 1] \times 1000 \quad (4) $$

Analytical uncertainties are estimated from the long-term reproducibility of the thermogenic Utica Shale gas through the same purification process (n = 22). External precision for δ¹³C, δD, Δ¹³CH₃D, and Δ¹²CH₂D₂ is found to be approximately 0.1‰, 0.3‰, 0.3‰, and 0.7‰, respectively (1σ). The relationship between temperature and both Δ¹³CH₃D and Δ¹²CH₂D₂ has been predicted through ab initio calculations based on the harmonic approximation and can be expressed by the following equations (Young et al., 2017):

$$ \Delta^{13}\text{CH}_3\text{D}(\text{T}) \approx 1000\ln(1 + 0.0355502/\text{T} - 433.038/\text{T}^2 + 1270210.0/\text{T}^3 - 5.94804 \times 10^8/\text{T}^4 + 1.196630 \times 10^{11}/\text{T}^5 - 9.07230 \times 10^{12}/\text{T}^6) \quad (5) $$

$$ \Delta^{12}\text{CH}_2\text{D}_2(\text{T}) \approx 1000\ln(1 + 0.183798/\text{T} - 785.483/\text{T}^2 + 1056280.0/\text{T}^3 + 9.37307 \times 10^7/\text{T}^4 - 8.919480 \times 10^{10}/\text{T}^5 + 9.901730 \times 10^{12}/\text{T}^6) \quad (6) $$

where T is in Kelvin. Eqs. (5) and (6) show that Δ¹³CH₃D and Δ¹²CH₂D₂ values are both positive when methane is formed at thermodynamic equilibrium, and approach 0‰ at high temperatures (>1000 K).

## 2.7. Calculation of isotope fractionation factors using Rayleigh equation

Incubations were performed in glass bottles with frequent shaking and the headspace methane was considered a well-mixed source in a closed system. Although methane in the headspace was extracted once or twice between the start and end of each experiment, the amount of methane extracted was minor (<3%) compared to the amount of methane in the whole bottle, resulting in a negligible deviation from a strictly closed system. Thus, the Rayleigh fractionation equation (Mariotti et al., 1981) was used to extract fractionation factors for the ¹³CH₄, ¹²CH₃D, ¹³CH₃D, and ¹²CH₂D₂ isotopologues relative to ¹²CH₄ resulting from anaerobic fractionation:

$$ \frac{R_t}{R_0} = \left( \frac{[^{12}\text{CH}_4]_t}{[^{12}\text{CH}_4]_0} \right)^{\alpha^{-1}} \quad (7) $$

where R refers to the ratio of isotopologues of interest in the gas phase (e.g., ¹³CH₃D/¹²CH₄, ¹²CH₂D₂/¹²CH₄), the subscript 0 signifies the initial property, and the ratio of methane abundance at time t relative to the initial methane abundance is commonly referred to as F (fraction remaining). For kinetic fractionations, the fractionation factor α is simply the ratio of the rate constant associated with the reaction of an isotopically substituted isotopologue relative to the rate constant of the unsubstituted isotopologue. These can include ratios of symmetry numbers. For the sake of simplification, the isotopologue fractionation factor for ¹³CH₃D and ¹²CH₂D₂ relative to ¹²CH₄ are termed ¹³Dα and D²α, respectively. By convention, α < 1 signifies that the reaction favors the isotopically light species, leaving the residue enriched in the heavy isotope or isotopologue. On a plot of -lnF vs. ln(R/R₀), the slope (often termed ε) provides the fractionation factor (¹³ε = 1-¹³α or Dε = 1-Dα). Fractionations and their uncertainties (95% confidence interval) are calculated by the weighted least square method (York et al., 2004). Errors for methane concentrations are estimated to be 0.4% of the measured value. Note that the F ratio was only measured in the incubation experiment of the Svalbard methane seep sediments.
At zeroth order, the relationship among isotopologue fractionation factors would follow the rule of the geometric mean (Bigeleisen, 1955), in which a stochastic distribution of isotopes occurs among the isotopologues. This would be the case where there is no energetic preference for forming doubly substituted isotopologues as opposed to randomly distributing isotopes among bonds (e.g., ¹³Dα ≈ ¹³α · Dα). Under these circumstances, the fractionation factors for the doubly substituted species are simply the products of the fractionation factors for the individual isotopes (i.e., the square of the geometric mean of the fractionation factors for doubly substituted species). The deviation from the rule of the geometric mean can be characterized by the kinetic clumped isotopologue fractionation factor, γ, yielding (cf. Wang et al., 2016):

$$^{13D}\alpha = ^{13D}\gamma \cdot ^{13}\alpha \cdot ^{D}\alpha \quad (8)$$

$$^{D2}\alpha = ^{D2}\gamma \cdot ^{D}\alpha \cdot ^{D}\alpha \quad (9)$$

## 2.8. Open system model -- steady-state between transport and oxidation of moving methane

In natural environments, Rayleigh-type methane oxidation is rarely observed because replenishment of methane, facilitated by advection or diffusion, usually applies. We will consider two end-member environmental scenarios: a steady-state open system model presented in this section and a time-dependent closed system model with a methane source and sink in the next section.

In the open system model, methane is transported into and out of the system via advection with concurrent AOM. A steady state between oxidation and supply is assumed to have been attained. A mass balance equation at steady state can be expressed as (Hayes, 2001; Wang et al., 2016):

$$R_0 = (1 - \varphi)R_{SS} + R_{SS}\varphi\alpha \quad (10)$$

where $R_0$ is the initial isotopologue ratio, $R_{SS}$ is the steady-state value, $\varphi$ is the fraction of methane removed via oxidation relative to that which flows through, and $\alpha$ is the kinetic fractionation factor associated with the oxidation reaction. Such a steady state will occur where the supply and sink occur at constant rates. Here, advection is assumed to have no isotope effect (Alperin et al., 1988), whereas AOM has associated fractionation factors ¹³α, Dα, ¹³Dγ, and D²γ. The fractional contribution of AOM to the total sink, $\varphi$, is physically related to the Damköhler number, i.e., Da, the ratio of the residence time in the flow system to the characteristic reaction time (Da = $\tau_{flow} / \tau_{reaction}$). For a first-order reaction, the dimensionless Damköhler number can be simplified to $\tau_{flow} \times k$ where $\tau_{flow}$ is the residence time and $k$ is the reaction rate constant with units of inverse time. The $\varphi$ values are normalized Da values, as in $\varphi = Da/(1+Da)$. Therefore, if $k$ is large in comparison to $\tau$, Da is large and $\varphi$ approaches 1. At this extreme, a balance between production and oxidation is achieved in a closed system as described in the next section. Conversely, if $k$ is small in comparison to $\tau$, $\varphi$ approaches 0 and the isotopologue ratio does not change. Values in between these extrema cause the steady-state values to be intermediate between the initial isotopologue ratios and the values at $\varphi = 1$.

## 2.9. Closed system model -- formation of methane balanced against oxidation

In the closed system model, the time-dependent evolution of the moles of an isotopologue of methane can be described in terms of a constant rate of production and a rate of oxidation that behaves as a first-order reaction with an invariable rate constant (Haghnegahdar et al., 2017). In this simplest model, where the balance

is only between production and consumption, the system can be described as:

$$\frac{dn_i}{dt} = E_i - k_i n_i \quad (11)$$

where $n_i$ is the moles of the isotopic species of interest, $E_i$ is the source term for $i$ (e.g., rate of abiotic or microbial methanogenesis), and $k_i$ is the rate constant for the sink $i$ (i.e., the rate constant for oxidation). The solution to Eq. (11) yields:

$$n_i(t) = n_i^0 e^{-k_i t} + \frac{E_i}{k_i} (1 - e^{-k_i t}) \quad (12)$$

where $n_i^0$ is the initial moles of $i$. The moles of species $i$ at steady state are obtained by evaluation Eq. (12) where $t \to \infty$, and thus $e^{-k_i t} \to 0$, yielding:

$$n_i(t \sim \infty) = \frac{E_i}{k_i} \quad (13)$$

Considering $i$ to be the major isotopologue, the steady-state amount of methane is controlled by the ratio of the production rate and the rate constant for oxidation. The steady-state ratio of two isotopologues can therefore be expressed as:

$$\frac{n_j(t \sim \infty)}{n_i(t \sim \infty)} = \frac{E_j}{E_i} \frac{k_i}{k_j} \quad (14)$$

where $n$ is the moles of the two isotopologues $i$ and $j$. The steady-state isotopologue ratio depends on the isotopologue ratio of the methane produced ($E_j/E_i$) and the reciprocal of the fractionation due to oxidation ($k_i/k_j$), independent of the absolute rates.

## 2.10. Comparison between measured and theoretically estimated isotope fractionation factors

To understand the meaning of the measured isotope fractionation factors from the Svalbard experiment, we used two methods to estimate the fractionation factors with mechanistic implications. Firstly, if the kinetics is dominated by rupture of single carbon-hydrogen bonds comprising the reaction coordinate, the fractionation factors can be estimated by the square root of the inverse ratios of relevant reduced masses, $\sqrt{\mu_i/\mu_j}$. Here the $\mu_i$ values are the reduced masses for ¹³C-H and ¹²C-D as required, $\mu$ is the reduced mass for ¹²C-H, and the derived fractionation factors are those implied by these reduced masses, including the use of the rule of the geometric mean for the doubly-substituted species.

Secondly, we consider the ratio of the logarithm of rotational/vibrational partition functions for a crude estimate for the transition state of methane being dismantled by the Mcr-mediated reaction. The $Q^*$ values are estimated by calculating the partition function for methane at 4 °C with the three asymmetric stretch modes (commonly referred to as $v_3$) removed. The differences among the partition functions for the different isotopologues are evidenced clearly by writing their ratios where the prime subscript refers to the isotopically-substituted species (Bigeleisen and Mayer, 1947; Urey, 1947):

$$\frac{Q'}{Q} = \frac{\sigma}{\sigma'} \prod_i \frac{u'_i}{u_i} \frac{e^{-u'_i/2}}{1 - e^{-u'_i}} \frac{1 - e^{-u_i}}{e^{-u'_i/2}} \quad (15)$$

We used the vibrational frequencies calculated by Lee et al. (1995). For our rough estimate of a suitable transition state bound to Mcr (an imprecise exercise in its own right), we removed the $v_3$ modes from Eq. (15). This approach provides a test of the possibility that it is the highest frequency modes, $v_3$, that are most readily broken and dominate the kinetics (and thus are not counted in the vibrational partition functions). Ratios of vibrational frequencies, or the square root of the inverse ratio of reduced masses, are ratios
of energies associated with a vibration. Values lnQ* are also proportional to energies, but represent the weighted mean across all vibrations not involved in bond rupture (Bigeleisen, 1952). The ratio is therefore an estimate of the relative energy associated with motions in the molecules.

### 3. Results

Methane isotopologue data from laboratory experiments and natural fluids are presented in Fig. 3. The $\Delta^{13}CH_3D$ and $\Delta^{12}CH_2D_2$ values from the Svalbard sediment slurry incubation and Chamorro Seamount natural fluids are higher than thermodynamic equilibrium values at their ambient temperature (>0 °C; Fig. 3A–C). Other incubations, including the Mcr-catalyzed isotope exchange and incubations of the Santa Barbara Channel sediment, Nankai Trough sediment, and Beatrix Gold Mine fracture fluid, show moderate increases in $\Delta^{13}CH_3D$ but minimal increases in $\Delta^{12}CH_2D_2$ (Fig. 3D–F). Their apparent temperatures based on clumped isotopes are similar to the experimental temperatures, ranging from 20 °C to 60 °C. In the following sections, we first present the laboratory incubation data and then show natural fluid data for comparison.

#### 3.1. Incubation of the Svalbard methane seep sediment slurries

Off-shore methane seep sediments from the Svalbard archipelago host high AOM activity. Ex-situ AOM rates determined by $^{14}$C-radiotracer techniques peak at 182 nmol cm$^{-3}$ d$^{-1}$ at 2–5 cm below seafloor (Fig. 2A) (Melaniuk et al., 2022). In the reactivated

sediment slurry with saturated methane and high sulfate concentrations (8–24 mM), dissolved sulfide concentrations increased linearly as methane concentrations decreased linearly (Fig. 2B). In addition, sulfate concentrations decreased, while alkalinity increased during the experiment (see the Supplementary Data). In the control experiment containing a headspace gas of N$_2$ and CO$_2$, no methane was detected and no increase in sulfide concentration was observed, indicating negligible methanogenesis and organoclastic sulfate reduction, respectively (Fig. S1A). Collectively, AOM is active in the sediment slurry and is the only mechanism for methane removal in the experiment. Accordingly, the in-vitro AOM rate determined by methane concentrations in the slurry incubation was 90 nmol of methane per cm$^3$ sediment slurry per day (Fig. 2B), translating into about 180 nmol per cm$^3$ of undiluted sediment per day. Furthermore, the in-vitro AOM and sulfate reduction rates determined by $^{14}$C- and $^{35}$S-radiotracers were 86 and 84 nmol cm$^{-3}$ d$^{-1}$, respectively, in the slurry incubation, which are equivalent to roughly 172 and 168 nmol cm$^{-3}$ d$^{-1}$ for undiluted sediment. The in-vitro AOM rates are consistent with the ex-situ AOM rates of the same sediment measured by $^{14}$C-radiotracer (Fig. 2A) and other in-vitro AOM enrichment cultures (e.g., Ono et al., 2021; Wegener et al., 2021).

In two experiments with high sulfate, the bulk $\delta^{13}C$ and $\delta D$ values of residual methane increased by 9.6‰ and 195.4‰, respectively, with declining methane concentrations (Fig. 4A). Based on a Rayleigh distillation model, the respective carbon and hydrogen isotope fractionation factors, expressed in per mil, were found to be 4.9 ± 0.1‰ and 149.0 ± 1.5‰ in the first experiment, and 7.5 ± 0.1‰ and 154.0 ± 1.3‰ in the second experiment (Fig. 5A–B). With progressive methane consumption, $\Delta^{13}CH_3D$ and $\Delta^{12}CH_2D_2$ values

<figure><img src="image_7.png" /><figcaption>**Fig. 3.** Clumped isotope data of residual methane in AOM incubation experiments and natural fluids. (A–B) Incubation of the Svalbard methane seep sediment slurry with high sulfate concentration (8–24 mM; panel A) and low sulfate concentration (0–0.4 mM; panel B). (C) Methane in natural fluids collected from the South Chamorro Seamount. (D) Exchange reaction between methane and methyl-coenzyme M. (E) Incubation of the Santa Barbara Channel methane seep sediment slurry. (F) Incubations of the Nankai Trough sediment slurry and Beatrix Gold Mine fracture fluid from the deep biosphere. The purple stars represent the initial tank gas. The solid black line depicts theoretical thermodynamic equilibrium abundances of methane isotopologues, along with corresponding temperatures. Error bars are one standard error. Note different scales between panels A–C and D–F.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 4.** Bulk isotope data of residual methane in AOM incubation experiments and natural fluids. (A–B) Incubation of the Svalbard methane seep sediment slurry with high sulfate concentration (8–24 mM; panel A) and low sulfate concentration (0–0.4 mM; panel B). The arrow and dashed line in panel B indicate the evolution of residual methane towards isotope equilibrium among CH₄, HCO₃⁻ and H₂O at 4 °C (Zhang et al., 1995; Turner et al., 2021). (C) Methane in natural fluids collected from the South Chamorro Seamount. (D) Exchange reaction between methane and methyl-coenzyme M. (E) Incubation of the Santa Barbara Channel methane seep sediment slurry. (F) Incubations of the Nankai Trough sediment slurry and Beatrix Gold Mine fracture fluid from the deep biosphere. The purple stars represent the initial tank gas. Uncertainties of δ¹³C and δD values are encompassed by individual data points.</figcaption></figure>

increased by 16.4‰ and 57.6‰, reaching values of 19.5‰ and 65.1‰, respectively (Fig. 3A). These values are substantially higher than the values expected for isotopologue equilibrium at the incubation temperature of 4 °C. The kinetic clumped isotopologue fractionations (γ), defined as the ratio of the mass-18 isotopologue fractionation factors to the product of the bulk carbon and hydrogen fractionation factors, describe deviations from the rule of the geometric mean (Bigeleisen, 1955; Wang et al., 2016). The respective γ values for ¹³CH₃D and ¹²CH₂D₂ defined by these data were 0.985 ± 0.002 and 0.912 ± 0.006 in the first experiment, and 0.985 ± 0.002 and 0.909 ± 0.005 in the second experiment (Fig. 5C–D). With declining methane concentrations, the δ¹³C values of DIC decreased by 6‰ in the first experiment and 10‰ in the second experiment, whereas the δD values of water were more-or-less invariant (see the Supplementary Data).

The same sediment slurry was then incubated with low sulfate concentration (<0.35 mM) under the same conditions (e.g., saturated methane), resulting in a lower thermodynamic driving force for AOM. Methane isotopologue compositions were analyzed on days 1, 11, 58, and 79 of the incubation. In contrast to the high sulfate experiment, the bulk δ¹³C and δD values of residual methane decreased by 10.1‰ and 15.2‰, respectively (Fig. 4B). The Δ¹³CH₃D and Δ¹²CH₂D₂ values, however, increased linearly by 4.6‰ and 14.8‰, respectively (Fig. 3B). There is no substantial difference in Δ¹³CH₃D and Δ¹²CH₂D₂ values of residual methane taken on days 58 and 79, likely signaling a reduction in the AOM rate due to sulfate limitation (0.01 mM). To examine whether the residual methane will approach isotopologue equilibrium or overshoot it, additional sulfate was added on day 83, to increase the concentration again to 0.39 mM in the slurry. Isotopologue ratios

of the residual methane were analyzed on days 92 and 178 of the incubation. Both Δ¹³CH₃D and Δ¹²CH₂D₂ values increased up to 11.4‰ and 31.5‰, respectively (Fig. 3B). Net isotope fractionation factors are not reported here for low sulfate incubation because of the challenge of accurately quantifying methane consumption, limiting our ability to estimate the fraction of methane remaining in the Rayleigh process, and thus causing large errors, but both ¹³α and ¹⁵α are higher than 1.0, in contrast with the high sulfate incubation in which both ¹³α and ¹⁵α are lower than 1.0 (Figs. 4A–B, 5). The δ¹³C values of DIC decreased by 2‰, less pronounced than the high-sulfate experiments, whereas the δD values of water remained constant (see the Supplementary Data).

## 3.2. Incubation of the Santa Barbara Channel methane seep sediment slurries

Analogous to the incubation of Svalbard sediment slurry, we conducted incubation with Santa Barbara Channel sediment slurries at 20 °C. However, the volume and pressure of headspace methane in the Santa Barbara slurry (4.4 ml under ca. 100 kPa) were much lower than those of the Svalbard slurry (60 ml under ca. 200 kPa), presumably leading to a lower thermodynamic driving force for AOM. An earlier study at the same sites of the Santa Barbara Channel found that the ex-situ AOM rates determined by ¹⁴C-radiotracer techniques were relatively high, on the order of tens to hundreds of nmol cm⁻³ d⁻¹ (Treude and Ziebis, 2010). All incubations with methane headspace showed decreases in sulfate concentration and increases in dissolved sulfide concentration and alkalinity over time (see the Supplementary Data). Along with the ex-situ AOM rates reported in the earlier study and the control
<figure><img src="image_9.png" /><figcaption>**Fig. 5.** Evolution of methane isotopologue ratios during methane consumption in the incubation of the Svalbard methane seep sediment slurry with high sulfate concentration (8–24 mM). Data from experiments 1 and 2 are shown in open and filled circles, respectively. F is the fraction of methane remaining and R is the ratio of isotopologues in the gas phase. α and γ represent the kinetic isotope fractionation factor and clumped isotopologue fractionation factor, respectively. The linear regressions and 2σ errors are calculated by the weighted least square method of York et al. (2004).</figcaption></figure>

experiment (Fig. S1B), these trends indicate the activity of AOM. Some δD values of residual methane decreased by up to 15.0‰, whereas the δ¹³C values remained roughly constant over time under both high and low sulfate conditions (Fig. 4E). In the high sulfate slurries (>28 mM), both Δ¹³CH₃D and Δ¹²CH₂D₂ of residual methane increased, ranging from 3.5‰ to 5.4‰ for Δ¹³CH₃D and from 8.5‰ to 12.1‰ for Δ¹²CH₂D₂ (Fig. 3E). In the low sulfate slurries (<1 mM), the Δ¹²CH₂D₂ values remained nearly constant over time (6.7 ± 0.9‰), while the Δ¹³CH₃D increased by 2.7‰, reaching a value of 5.5 ± 0.2‰ relative to stochastic (Fig. 3E). For comparison, when methane molecules are at thermodynamic equilibrium at 20 °C, the respective Δ¹³CH₃D and Δ¹²CH₂D₂ values would be 5.9‰ and 20.3‰ (Young et al., 2017).

### 3.3. Incubation of sediment slurries and fracture fluids from the deep biosphere

Marine sediments and fracture fluids were sampled from the deep subsurface of the Nankai Trough and the Beatrix Gold Mine, respectively. Metagenomics, metatranscriptomics, and fluorescent in situ Hybridization (FISH) of the Beatrix Gold Mine fluids were previously published and the results demonstrated that ANMEs are present and active in the incubations (Lau et al., 2016; Harris et al., 2018; Harris et al., 2021). In parallel to the natural abundance incubations described in Section 2.3, long-term (350 day) high-pressure (40 MPa) ¹³C-CH₄ tracer incubations were performed on Nankai Trough sediment slurries (see Supplementary Methods and Fig. S2). Based on the production of ¹³C-DIC (Fig. S3) and δ¹³-C_DIC measurements (Fig. S4), trace AOM activity was statistically distinguishable above background, albeit at very low rates (on the order of pmol cm⁻³ day⁻¹). These findings agree with ¹⁴C-CH₄ radiotracer experiments reported by Beulig et al. (2022), who found potential AOM activity near detection limits in the Nan-

kai Trough sediments due to a high background of abiotic ¹⁴C-CH₄ conversion to ¹⁴C-CO₂ in the medium controls. These results collectively highlight the difficulty in quantifying AOM activity in these oligotrophic, low-biomass, deep biosphere sediments (see Heuer et al., 2020). Thus, only potential AOM activity has been observed and no explicit discussion can be made in this respect.

The bulk isotope data of incubation experiments are consistent with the sluggish nature of microbial activity in the deep sub-seafloor. Here, nearly all headspace methane showed no substantial changes in δ¹³C and δD values (<2‰; Fig. 4F). However, substantial changes were found for clumped isotope ratios. In the Nankai Trough slurries with nitrate or nitrite as the added electron acceptors, only Δ¹³CH₃D increased, reaching a value of 5.7 ± 0.1‰ as incubation time progressed (Figs. 3F, S5C). Similarly, when iron (III) in the form of hydrous ferric oxide was added as the electron acceptor in the Beatrix Gold Mine incubation, Δ¹³CH₃D increased up to 5.7 ± 0.1‰, while there was less increase in Δ¹³CH₃D when sulfate or no electron acceptor was added (Figs. 3F, S5D). Little variation in Δ¹²CH₂D₂ values was found in these experiments (Fig. 3F).

### 3.4. Incubation with methyl-coenzyme M reductase (Mcr) enzyme

We performed experiments with Mcr purified from Methanothermobacter marburgensis that catalyzes the exchange of the methyl moiety of methyl-coenzyme M with methane. The δ¹³C and δD of residual methane decreased by up to 2.3‰ and 6.4‰ in all cases, respectively (Fig. 4D). Both Δ¹³CH₃D and Δ¹²CH₂D₂ increased to different degrees relative to the initial methane isotopologue ratios (Fig. 3D). In detail, at 60 °C, the incubations lasted for 2–6 h in which a small increase of 0.5‰ in Δ¹³CH₃D and 2.0‰ in Δ¹²CH₂D₂ were observed. At 20 °C, the incubation lasted for 48 h and the Δ¹³CH₃D of residual methane increased substantially by
4.1‰, reaching a value of 5.6 ± 0.1‰ relative to stochastic, while Δ¹²CH₂D₂ showed a modest increase of only 2.6‰, reaching a value of 8.5 ± 0.5‰.

### 3.5. Natural methane-rich fluids from the South Chamorro Seamount

Three fluid samples were collected from South Chamorro Seamount during HyperDolphin Dives 941, 945, and 947 in 2009. The general geochemical data were previously published (Wheat et al., 2020). In brief, the stagnant fluids collected during Dive 941 are rich in sulfate (12 mM) and methane (33 mM), whereas discharging fluids from Dives 945 and 947 have low sulfate concentrations (0.8 mM) and relatively low methane concentrations (23 mM). Further, sulfate and methane concentrations in the fluids were negatively correlated with aqueous sulfide concentration and alkalinity, indicating that the fluids were altered by AOM to different degrees (Wheat et al., 2020). The bulk δ¹³C and δD values range from -37.0‰ to -26.1‰ and -110.3‰ to -100.6‰, respectively (Fig. 4C). The respective Δ¹³CH₃D and Δ¹²CH₂D₂ values are 7.5‰ and 25.2‰ in the high-sulfate, high-methane fluid, but increase to 12.6‰ and 42.6‰ in low-sulfate, low-methane fluids (Fig. 3C). These values are analogous to the results from the Svalbard slurry incubations (Fig. 3A–B).

## 4. Discussion

### 4.1. Kinetic fractionations of methane clumped isotopologue during AOM

AOM was highly active in the sediment slurry from the Svalbard methane seep (Fig. 2). The high-sulfate incubation experiment showed progressive enrichment of methane isotopologues containing ¹³C and D isotopes relative to the initial gas (Figs. 4A, 5), indicating a kinetic isotope fractionation. Although the kinetic fractionation for D/H ratios is typical, the carbon isotope fractionation is slightly lower than previous studies, which are in the range of 11‰ to 39‰, suggesting a higher degree of reversibility of the carbon pathway in our AOM incubations (Holler et al., 2009; Wegener et al., 2021). The γ value for ¹³CH₃D of 0.985 is lower than those obtained for other methane oxidation reactions such as aerobic oxidation and oxidation by OH and Cl radicals, while the γ values for ¹²CH₂D₂ of 0.909–0.912 are similar to those obtained for other methane oxidation reactions (Wang et al., 2016; Haghnegahdar et al., 2017; Whitehill et al., 2017; Ono et al., 2021; Krause et al., 2022). Lower absolute magnitudes of γ values imply larger magnitudes of clumped isotope effects and larger deviations in the clumped isotope fractionation factor relative to the product of the bulk isotope fractionation factors. Ono et al. (2021) proposed that the different γ values could be linked to the transition state structure of the reaction. For example, the C–H bond length of methyl-coenzyme M reductase (2.6 Å) and soluble methane monooxygenase (1.3 Å), the key enzymes involved in AOM and aerobic methanotrophy, respectively, are different, necessitating changes in the ¹³C–D stretching vibrational mode and thus perhaps the difference in kinetic clumped isotopologue effects (Ono et al., 2021).

The extremely high Δ¹³CH₃D and Δ¹²CH₂D₂ values up to 12.6‰ and 42.6‰, respectively, from the Chamorro Seamount, a serpentinite mud volcano in the Mariana forearc (Fig. 6C–F), as well as the high values seen in the Svalbard incubations (Fig. 6A), appear to be signatures of AOM. Based on samples from experiments and natural environments for which measurements currently exist, AOM is the only known process that drives methane isotopologue abundances substantially above equilibrium predictions for both

Δ¹³CH₃D and Δ¹²CH₂D₂. In the laboratory, we treat AOM as a Rayleigh process in which methane is consumed only by AOM in a closed system. This explains the trends in Δ¹³CH₃D vs. Δ¹²CH₂D₂ space for the Svalbard incubations (Fig. 6A–B). For the Chamorro Seamount, an open system that involves flow, is likely to be more suitable as a model for isotopologue fractionation. There, Wheat et al. (2020) demonstrated that methane production was abiotically fueled by hydrogen production from serpentinization and carbonate dissolution in the deep-sourced fluids, and the fluid compositions were largely altered by AOM as the fluid migrates upwards and encounters deposited pelagic sediment underlying South Chamorro Seamount. Although the pH of the collected fluids is 12, the pH of deep subsurface niches where ANME archaea live and AOM occurs can be lower, as observed in neighboring serpentinite mud volcanoes (Wheat et al., 2020). Earlier studies have shown that microorganisms can tolerate high-pH fluids and perform AOM at the Chamorro Seamount (Takai et al., 2005; Curtis et al., 2013; Kawagucci et al., 2018). These AOM-affected fluids make it to the seafloor where they discharge. Therefore, the environment where fluids are traversing through the seamount represents an open system, in which methane is transported into the subduction channel via advection and removed by both advection and AOM, likely resulting in a steady state, or at least a transient steady state. A simple mass balance model can be used to investigate the effects of simultaneous advection and oxidation (see Materials and Methods). Here, we assume that the serpentinite-sourced abiotic methane is in isotopologue equilibrium at its formation temperature (>250 °C; Wheat et al., 2020). Such near-zero clumped isotope compositions are consistent with the results from abiotic silane decomposition experiments at 300–600 °C reported in Young et al. (2017). We consider that the impact of diffusion in fluid transport on the Δ¹³CH₃D and Δ¹²CH₂D₂ values is negligible within the subduction channel given the significantly smaller-scale molecular diffusion operates relative to advective processes (Wheat et al., 2020). Thus, AOM is assumed to be the only process for isotope fractionation.

The combined flow and oxidation steady-state models shown in Fig. 6C–D largely reproduce the observed large Δ¹³CH₃D and ¹²CH₂D₂ values at the Chamorro Seamount. Note that to achieve the best fit, a slightly lower γ value for ¹³CH₃D of ~0.980 or a higher γ value for ¹²CH₂D₂ of ~0.935 than obtained from the experiments is needed. While the open system model is more suitable to describe these fluids, we also considered the scenario of a closed system in which the formation of methane is balanced against oxidation. In this model, isotopologue abundances vary with time, eventually reaching a steady state that is sensitive to the precise γ values (Fig. 6E–F). We set the ratio of the steady-state moles to initial moles of methane, (E/k)/n₀, to 1 since the modeled trajectories remain largely the same as (E/k)/n₀ varies. Here, ¹³CH₃D γ values of 0.980–0.985 and ¹²CH₂D₂ γ values of 0.912–0.935 are needed to achieve the best fit in the closed system model. These estimated γ values for the Chamorro Seamount fluids from either the open- or closed-system model are largely consistent with those obtained from the Svalbard slurry incubations.

One caveat of these two models is the, to date, largely unknown bulk isotope fractionation factors. In Fig. 6, we set ¹³α and ¹⁵α to be 0.995 and 0.851, respectively, based on experimental observations made in the Svalbard experiment. We then test the sensitivity of the models by varying the bulk isotope fractionation factors to allow depletion in heavy isotopes of methane (i.e., ¹³α and ¹⁵α > 1). The model outputs remain largely the same, though with different curvatures (Fig. S6). This demonstrates that the trajectories in Δ¹³CH₃D vs. Δ¹²CH₂D₂ space are much more sensitive to the relationship between the bulk α values (i.e., the γ factors) than
<figure><img src="image_10.png" /><figcaption>**Fig. 6.** Kinetically-driven methane isotopologue data of residual methane in AOM incubation experiments and methane in natural fluids. (A–B) Residual methane during incubation of the Svalbard methane seep sediment slurry with high sulfate concentration (8–24 mM; panel A) and low sulfate concentration (0–0.4 mM; panel B). (C–F) Methane in natural fluids from the South Chamorro Seamount. Grey lines depict simulation outputs of the Rayleigh fractionation model (A), open-system flow and steady-state model (C–D), and closed-system time-dependent shift model (E–F). The fractionation factors used in each model are shown in each diagram. The purple stars represent the initial tank gas in panel A, and the source methane in panels C–F, respectively. This source methane is adopted from Young et al. (2017), representing abiotic methane gas formed through Si₅C₁₂H₃₆ decomposition experiments, which resembles methane production through serpentinization. In panel A, points along the Rayleigh fractionation line are marked at intervals of 0.1 in *f*, the fraction of initial methane remaining, while in panel C, points along the open-system model are marked at intervals of 0.1 in *φ*, the fraction of methane removed via oxidation. The modeled trajectories in panels E–F remain largely the same as (E/k)/n₀ increase from 10⁻¹⁰ to 10¹⁰ (data not shown), and we use 1 as an arbitrary representation. The solid black line depicts theoretical thermodynamic equilibrium abundances of methane isotopologues, along with corresponding temperatures. Uncertainties of Δ¹³CH₃D and Δ¹²CH₂D₂ values are encompassed by individual data points.</figcaption></figure>

to the absolute ¹³α and Dα values in these two models (Figs. 6, S6). Given the uncertainties in model inputs and limited data points, we stress that the estimated γ values are associated with large uncertainties but, overall, are significantly lower than unity.

Taken together, both open and closed system models display similar positive trajectories in Δ¹³CH₃D vs. Δ¹²CH₂D₂ space, similar to the Rayleigh distillation model (Fig. 6), suggesting that extremely positive Δ¹³CH₃D and Δ¹²CH₂D₂ values are potentially a
clumped isotope signature of kinetically-dominated AOM. These positive values are the result of $\gamma < 1$ for both mass-18 isotopologues. Indeed, if $\gamma = 1$ for both rare isotopologues, and all else being equal, $\Delta^{13}CH_3D$ and $\Delta^{12}CH_2D_2$ would both decrease rather than increase with oxidation. In detail, the Rayleigh fractionation equation (Eq. (7)) can be simplified to (Wang et al., 2016):

$$ \Delta^{13}CH_3D = \Delta^{13}CH_3D_{init} + (^{13}D\gamma \cdot ^{13}\alpha \cdot ^{D}\alpha - ^{13}\alpha - ^{D}\alpha + 1) \cdot \ln f \quad (16) $$

$$ \Delta^{12}CH_2D_2 = \Delta^{12}CH_2D_{2init} + (^{D2}\gamma \cdot ^{D}\alpha^2 - 2 \cdot ^{D}\alpha + 1) \cdot \ln f \quad (17) $$

where $\Delta^{13}CH_3D_{init}$ and $\Delta^{12}CH_2D_{2init}$ are the initial isotopologue compositions. Therefore, the threshold $^{13}D\gamma$ and $^{D2}\gamma$ for constant $\Delta^{13}CH_3D$ and $\Delta^{12}CH_2D_2$ values are $(^{13}\alpha + ^{D}\alpha - 1)/(^{13}\alpha \cdot ^{D}\alpha)$ and $(2 \cdot ^{D}\alpha - 1)^{D}\alpha^2$, respectively. Using the measured $^{13}\alpha$ and $^{D}\alpha$ from the high-sulfate Svalbard incubation, the calculated $^{13}D\gamma$ and $^{D2}\gamma$ thresholds are found to be 0.999 and 0.969, respectively. Any measured $\gamma$ factors lower than the threshold conditions mean that both $\Delta^{13}CH_3D$ and $\Delta^{12}CH_2D_2$ values should evolve to more positive values in a closed system, and vice versa. Although the $\gamma$ factors appear to be linked to the transition state structure of the enzyme-mediated reaction (Ono et al., 2021), the exact structures of transition states and their role in kinetic clumped isotopologue effect remain unknown. Further theoretical and experimental investigations are needed to elucidate the fundamental mechanism.

We note that with progressive AOM from the stagnant fluids to the discharging fluids at the Chamorro Seamount as indicated by alkalinity, methane, sulfate, and sulfide concentrations (Wheat et al., 2020), the $\delta^{13}C$ of residual methane decreases by 10.9‰ and $\delta D$ increases by 9.7‰ in three fluid samples (Fig. 4C), resembling the ANME incubation under low-sulfate conditions (<1 mM) in Wegener et al. (2021). These authors determined radiotracer-based forward and backward AOM reaction rates at different sulfate concentrations and found that the net reversibility of the AOM reaction increases with decreasing sulfate concentration. They concluded that the net bulk isotope fractionations of AOM reflect a combination of expression of kinetic isotope effects (KIEs) and equilibrium isotope effects (EIEs) (Yoshinaga et al., 2014; Ono et al., 2021; Wegener et al., 2021). Indeed, when the Svalbard sediment slurry was incubated with even lower sulfate concentrations (<0.35 mM), both $\delta^{13}C$ and $\delta D$ of residual methane decreased over time (Fig. 4B), suggestive of a greater expression of EIEs in both carbon and hydrogen pathways of AOM due to a lower net thermodynamic drive and thus relatively higher net reversibility of AOM. At 4 °C, the carbon and hydrogen equilibrium isotope fractionation factors between $HCO_3^- (aq)-CH_4(g)$ and $H_2O(l)-CH_4(g)$ are 1.091 and 1.216, respectively (Zhang et al., 1995; Turner et al., 2021). When inter-species isotope equilibrium is achieved, the $\delta^{13}C$ and $\delta D$ of methane should be around -95‰ and -241‰ based on the measured $\delta^{13}C_{DIC}$ and $\delta D_{H_2O}$ values of -13.1‰ and -76.7‰ in the low-sulfate experiment.

Despite the fact that the residual methane in the low-sulfate incubation moves towards isotope equilibrium between both $CH_4-CO_2$ and $CH_4-H_2O$ molecules (Fig. 4B), $\Delta^{13}CH_3D$ and $\Delta^{12}CH_2D_2$ increase linearly and fall on nearly the same Rayleigh fractionation line in $\Delta^{13}CH_3D$ vs. $\Delta^{12}CH_2D_2$ space as the high-sulfate incubation (Fig. 6B). In particular, the $\Delta^{13}CH_3D$ reaches 7.4‰ when sulfate is depleted (0.01 mM), which is higher than the value expected for intra-species isotope equilibrium at the experimental temperature ($\Delta^{13}CH_{3Eq.} = 6.5‰$ at 4 °C). This high $\Delta^{13}CH_3D$ value suggests a kinetic clumped isotopologue fractionation. We surmise that rather than representing unidirectional consumption of methane, the isotopic effects we observed appear to reflect a convolution of transition-state kinetics and equilibrium isotopic effects. Furthermore, both $\Delta^{13}CH_3D$ and $\Delta^{12}CH_2D_2$ increase linearly again after new sulfate is added to the slurry

(0.39 mM), indicating that the final $\Delta^{13}CH_3D$ and $\Delta^{12}CH_2D_2$ values in a methane-rich closed system depend on the electron acceptor concentration and availability. The last measured $\Delta^{13}CH_3D$ and $\Delta^{12}CH_2D_2$ values correspond to apparent temperatures lower than -34 °C (Fig. 6B), far below the incubation temperature of 4 °C. Collectively, our results from the Svalbard sediment slurry and the Chamorro Seamount show clear evidence for kinetic clumped isotopologue fractionation, with the bulk isotope ratios either increasing or decreasing with progressive methane consumption depending upon the availability of the electron acceptor.

## 4.2. Equilibrium between methane isotopologues during AOM

It has been proposed that AOM drives methane isotopologue abundances to thermodynamic isotopic bond-order equilibrium in a range of marine and continental settings (Ash et al., 2019; Giunta et al., 2019; Tyne et al., 2021; Warr et al., 2021; Ono et al., 2022). The initial step of AOM is the Mcr-catalyzed reversal of the final reaction in methanogenic pathways as shown for the enzyme from M. marburgensis (Eq. (18)):

$$ CH_4 + \text{CoM-S-S-CoB} \rightleftharpoons CH_3-S-CoM + \text{HS-CoB} \quad (\Delta G^\circ = + 30 \text{ kJ mol}^{-1}) \quad (18) $$

where CoM-S-S-CoB is the heterodisulfide of coenzymes M and B, $CH_3-S-CoM$ is methyl-coenzyme M, and HS-CoB is the reduced form of coenzyme B (Scheller et al., 2010). The reversibility of this reaction is evidently central to the mechanism for methane isotopologue equilibration by AOM. The in-vitro experiments that comprise exchange of methane and the methyl moiety in methyl-coenzyme M ($CH_3-S-CoM$) characterize the role of the Mcr enzyme in reaction Eq. (18). At 60 °C, the optimal temperature for the thermophile-derived enzyme, both $\Delta^{13}CH_3D$ and $\Delta^{12}CH_2D_2$ values of residual methane increase slightly (Fig. 3D). At 20 °C, however, $\Delta^{13}CH_3D$ approaches the equilibrium value, while there is little change in $\Delta^{12}CH_2D_2$. This implies that carbon plays a larger role than hydrogen in the isotopologue selectivity for reaction with the Mcr enzyme at a temperature below optimal.

Incubations of slurry from the Santa Barbara Channel, with a lower amount of methane, further inform the assessment of the role of thermodynamic driving forces in determining the isotopologue signatures of AOM. There is a trend towards $\Delta^{13}CH_3D$ and/or $\Delta^{12}CH_2D_2$ values expected for isotopologue equilibrium at experimental temperature (Fig. 3E). Two of the high-sulfate incubations of the Santa Barbara sediments resemble the high-temperature Mcr experiments under the optimal condition, whereas the low-sulfate incubations resemble the low-temperature Mcr experiments (Fig. 3D-E). Importantly, the amounts of methane injected into the Santa Barbara slurry (0.5 mmol) are much lower than those into the Svalbard slurry (5.4 mmol), leading to a much lower thermodynamic drive for AOM for the former and an overall lower rate of reaction. Thus, the reversibility and isotope exchange during AOM in the Santa Barbara slurry may be much higher than the Svalbard slurry under either high- or low-sulfate conditions. This vigorous exchange ultimately allows methane isotopologues to react with Mcr reversibly, achieving partial intra-species isotope exchange in the Santa Barbara slurry incubations. The high reversibility may be exacerbated by the fact that the ANME community in coastal sandy Santa Barbara sediments with high depositional rates is less active with less biomass compared to that from the deep-sea fine-grained Svalbard methane seep sediments. The lower sedimentation rates at the Svalbard seep focus AOM in the same sediments for an extended period, allowing the ANME community to grow and enrich.

All experiments of the deep biosphere sediment and fluid incubations suggest that there is a trend toward equilibrium values for
$\Delta^{13}CH_3D$ but little change in $\Delta^{12}CH_2D_2$ within experimental time-scales, regardless of which electron acceptor is present (Fig. 3F). In detail, a positive correlation between incubation time (up to 500 days) and increase in $\Delta^{13}CH_3D$ values is found in incubation experiments of the Nankai Trough sediments (Fig. S5C), consistent with the well-documented sluggish nature of microbial activity in the deep subseafloor (Hoehler and Jørgensen, 2013; Heuer et al., 2020).

One may postulate that the $\Delta^{13}CH_3D$ values in our incubations that appear to approach equilibrium values would in fact continue to increase beyond the equilibrium value given sufficient time. However, this is probably not the case. The respective maximum observed for $\Delta^{13}CH_3D$ values is 5.56‰ (2 days), 5.50‰ (3 days), 5.70‰ (500 days), and 5.66‰ (350 days) in the Mcr experiment, and the incubations of the Santa Barbara Channel sediment, the Nankai Trough sediment, and the Beatrix Gold Mine fracture fluid (Fig. 3D–F). The apparent temperatures based on $\Delta^{13}CH_3D$ are 32 °C, 33 °C, 27 °C, and 28 °C, respectively (30 ± 3 °C, n = 4). These are all close to the ambient temperature under which these different incubations took place. Although it cannot be categorically ruled out that these $\Delta^{13}CH_3D$ values reflected ambient temperatures by happenstance during an ongoing kinetically driven trend, it would be fortuitous that this would occur four times in independent experimentation. Therefore, we conclude that the $\Delta^{13}CH_3D$ values increase over time towards equilibrium and indeed remain at thermodynamic equilibrium despite a large range of incubation times of 2 to 500 days in the above four experiments.

Given the observed intra-species isotope quasi-equilibrium, one may similarly expect inter-species isotope exchange to occur as seen in the low-sulfate Svalbard experiment (Fig. 4B). The residual methane in most of the Santa Barbara and deep biosphere incubations, however, showed no substantial changes in $\delta^{13}C$ and $\delta\text{D}$ values (Fig. 4E–F). Changes in bulk and clumped isotope compositions are decoupled in these incubations. While more experiments are needed to explicitly confirm the conditions under which this decoupling occurs, one tentative hypothesis is that, under certain conditions, the Mcr enzyme only mediates the exchange of isotopes among methane molecules without expressed bulk isotope fractionation, allowing $^{13}CH_3D$ and/or $^{12}CH_2D_2$ to be re-ordered to the incubation temperature via bond breaking and reformation. This hypothesis is analogous to the heating experiments in which methane intra-species equilibration is promoted by the presence of a metal catalyst at temperatures above 200 °C (Ono et al., 2014; Stolper et al., 2014b; Young et al., 2017). In our case, the Mcr enzyme may act as the catalyst and promote isotopic bond re-ordering with limited isotope exchange between $CH_4-CO_2$ and $CH_4-H_2O$ molecules, as methane molecules are disassembled and reassembled. This hypothesis requires further experimental investigation.

### 4.3. A general model for methane isotopologue fractionation due to AOM

#### 4.3.1. Model construction

To understand the mechanism of methane isotopologue fractionation and/or bond re-ordering during AOM, the methane isotopologue fractionation factors obtained in AOM laboratory incubation experiments can be used to construct a geochemical model with several simplifying assumptions for fractionation due to AOM. The result is a parsimonious model that explains most methane isotopologue data affected by anaerobic oxidation by ANMEs collected thus far. A crucial aspect of the model is the important role that reversibility has in controlling the trajectory of gases in $\Delta^{13}CH_3D$ vs. $\Delta^{12}CH_2D_2$ space during AOM.

The model addresses the most important single step, the potentially reversible cleaving of $CH_4$ molecules catalyzed by the Mcr enzyme to form co-enzyme $CH_3$ and H in the ANME cells (e.g., Scheller et al., 2010). The relevant reaction can therefore be written succinctly as

$$CH_4 \overset{k_f}{\underset{k_r}{\rightleftharpoons}} CH_3 + H \quad (19)$$

where the forward and reverse rate constants, $k_f$ and $k_r$, are controlled by the kinetics of binding and bond rupture/formation in the first step of AOM. While there are a total of 280 isotopic isomer equivalents of the reaction in Eq. (19), with one exception (see below), we will focus on a simplified reaction network of isotopologues (as opposed to all isotopic isomers). This is necessitated by the fact that the measured fractionation factors convolve a number of isotopic effects, including the symmetry effects of all isomers. The sixteen reactions in terms of 12 isotopologues are

$$1. CH_4 \overset{k_1}{\underset{k_2}{\rightleftharpoons}} CH_3 + H$$

$$2. CH_3D \overset{k_3}{\underset{k_4}{\rightleftharpoons}} CH_3 + D$$

$$3. CH_3D \overset{k_3}{\underset{k_4}{\rightleftharpoons}} CH_2D + H$$

$$4. CH_2D_2 \overset{k_5}{\underset{k_6}{\rightleftharpoons}} CH_2D + D$$

$$5. CH_2D_2 \overset{k_5}{\underset{k_6}{\rightleftharpoons}} CHD_2 + H$$

$$6. ^{13}CH_4 \overset{k_6}{\underset{k_7}{\rightleftharpoons}} ^{13}CH_3 + H$$

$$7. ^{13}CH_3D \overset{k_7}{\underset{k_8}{\rightleftharpoons}} ^{13}CH_3 + D$$

$$8. ^{13}CH_3D \overset{k_8}{\underset{k_9}{\rightleftharpoons}} ^{13}CH_2D + H \quad (20)$$

$$9. CH_3 + H \overset{k_9}{\underset{k_{10}}{\rightleftharpoons}} CH_4$$

$$10. CH_3 + D \overset{k_{10}}{\underset{k_{11}}{\rightleftharpoons}} CH_3D$$

$$11. CH_2D + H \overset{k_{12}}{\underset{k_{13}}{\rightleftharpoons}} CH_3D$$

$$12. CH_2D + D \overset{k_{13}}{\underset{k_{14}}{\rightleftharpoons}} CH_2D_2$$

$$13. CHD_2 + H \overset{k_{15}}{\underset{k_{16}}{\rightleftharpoons}} CH_2D_2$$

$$14. ^{13}CH_3 + H \overset{k_{14}}{\underset{k_{15}}{\rightleftharpoons}} ^{13}CH_4$$

$$15. ^{13}CH_3 + D \overset{k_{15}}{\underset{k_{16}}{\rightleftharpoons}} ^{13}CH_3D$$

$$16. ^{13}CH_2D + H \overset{k_{16}}{\underset{k_{17}}{\rightleftharpoons}} ^{13}CH_3D$$

where $k_f$ are the relevant rate constants (first order for reactions (1) through (8), 2nd order for (9)–(16)), and again we emphasize that the different isomers are embedded in the fractionation factors extracted from experiments.

Isotopologue fractionation factors for the first step in deconstructing methane isotopologues are in effect ratios of rate constants 1–8. We measured these fractionation factors for the Svalbard incubations in which the loss of methane gas due to AOM is consistent with an irreversible Rayleigh process. The derived fractionation factors, $\alpha$, are listed in Table 1 together with relevant comparisons to be discussed. Here, the 3rd column is the square root of the inverse ratios of relevant reduced masses ($^{12}C$-H,

<caption>**Table 1**</caption>

Measured and estimated fractionation factors relative to $^{12}CH_4$. The experimental data are from the Svalbard incubations. $\mu$ is reduced mass. Q* is rotational/vibrational partition functions for a crude estimate for the transition state of methane being dismantled by the Mcr-mediated reaction. The prime superscript refers to the isotopically-substituted species.

<table><thead><tr><th>Species/<sup>12</sup>CH<sub>4</sub></th><th>Experiment</th><th>√μ/μ'</th><th>lnQ*/lnQ*</th></tr></thead><tbody><tr><td><sup>13</sup>CH<sub>4</sub></td><td>0.995</td><td>0.997</td><td>0.997</td></tr><tr><td><sup>12</sup>CH<sub>3</sub>D</td><td>0.851</td><td>0.734</td><td>0.832</td></tr><tr><td><sup>13</sup>CH<sub>3</sub>D</td><td>0.834</td><td>0.732</td><td>0.827</td></tr><tr><td><sup>12</sup>CH<sub>2</sub>D<sub>2</sub></td><td>0.660</td><td>0.539</td><td>0.694</td></tr></tbody></table>
¹³C-H, ¹²C-D, etc.). These are the fractionation factors expected where the kinetics is dominated by rupture of single carbon-hydrogen bonds comprising the reaction coordinate. Note that the experimentally-derived fractionation factors are considerably closer to unity than those implied by the reduced masses.

The 4th column in the table of fractionation factors is the ratio of the logarithm of rotational/vibrational partition functions for a crude estimate for the transition state of methane being dismantled by the Mcr-mediated reaction. We note that the measured fractionation factors are crudely similar, but not identical, to the lnQ*/lnQ* ratios; the ratios of logs of partition functions provide better rationalizations for the measured fractionation factors than do reduced masses along individual bonds. We note further that this is evidently a characteristic of AOM as opposed to abiotic oxidation in the atmosphere and perhaps bacterial aerobic oxidation, where fractionation factors are significantly closer to the square root of the inverse ratios of relevant reduced masses (Haghnegahdar et al., 2017; Whitehill et al., 2017; Krause et al., 2022). For example, the values for α¹³C/¹²C, αD/H, α¹³CH₃D/CH₄, and α¹²CH₂D₂/CH₄ of methane oxidation by OH radical in the atmosphere are 0.994, 0.758, 0.752, and 0.521, respectively (Haghnegahdar et al., 2017), which by inspection are seen to be closer to expectations from reduced masses along a C-H reaction coordinate. We suggest that larger fractionation nearer to unity, and resembling plausible ratios of energies of methane transition states, reflect significant binding of methane to the catalyzing enzymes (e.g., Mcr), in the case of AOM, leading to a more dispersed dependence of bond rupture on overall rotational/vibrational energy for the different isotopologues.

The model formulated to accommodate the measured fractionation factors comprises the set of rate constants, with branching ratios, and fractionation factors shown in Table 2. Several assumptions are needed based on theoretical consideration because some inputs are not readily available from the experiments. Taking reactions (2) and (3) as an example, the breakage of a ¹²C-D bond in ¹²CH₃D (primary kinetic isotope effect for D/H, α₂) is associated with a larger kinetic isotope effect than breaking a ¹²C-H bond in ¹²CH₃D (secondary kinetic isotope effect, α₃). Thus, we set α₂ to be 0.5, which is close to the primary KIE (0.41 ± 0.04) observed in the methane activation experiment using the Mcr enzyme (Scheller et al., 2013). The primary KIEs for doubly substituted species are derived following the rule of the geometric mean with measured γ values (Eqs. (8)–(9)). The branching ratios are set to

<caption>**Table 2**</caption>

Rate constants (k) and isotope fractionation factors (x) for the simplified reaction network of isotopologues. Dαeq-p and Dαeq-s are the equilibrium fractionation factors for primary and secondary fractionations, respectively. D²γeq(-p/s) and ¹³Dγeq(-p/s) are the fractionation factors representing equilibrium abundances of ¹²CH₂D₂ and ¹³CH₃D in thermodynamic equilibrium. The reader is referred to the text for interpretation of the table.

<table><thead><tr><th>Reaction</th><th>k</th><th>α</th></tr></thead><tbody><tr><td>1. CH<sub>4</sub> → CH<sub>3</sub> + H</td><td>k<sub>f</sub></td><td>α<sub>1</sub> = 1.000</td></tr><tr><td>2. CH<sub>3</sub>D → CH<sub>3</sub> + D</td><td>k<sub>2</sub> = 1/4 k<sub>f</sub> α<sub>2</sub></td><td>α<sub>2</sub> = 0.500</td></tr><tr><td>3. CH<sub>3</sub>D → CH<sub>2</sub>D + H</td><td>k<sub>3</sub> = 1/4 k<sub>f</sub> α<sub>3</sub></td><td>α<sub>3</sub> = 0.968</td></tr><tr><td>4. CH<sub>2</sub>D<sub>2</sub> → CH<sub>2</sub>D + D</td><td>k<sub>4</sub> = 1/2 k<sub>f</sub> α<sub>4</sub></td><td>α<sub>4</sub> = 0.441</td></tr><tr><td>5. CH<sub>2</sub>D<sub>2</sub> → CHD<sub>2</sub> + H</td><td>k<sub>5</sub> = 1/2 k<sub>f</sub> α<sub>5</sub></td><td>α<sub>5</sub> = 0.879</td></tr><tr><td>6. ¹³CH<sub>4</sub> → ¹³CH<sub>3</sub> + H</td><td>k<sub>6</sub> = k<sub>f</sub> α<sub>6</sub></td><td>α<sub>6</sub> = 0.995</td></tr><tr><td>7. ¹³CH<sub>3</sub>D → ¹³CH<sub>3</sub> + D</td><td>k<sub>7</sub> = 1/4 k<sub>f</sub> α<sub>7</sub></td><td>α<sub>7</sub> = 0.490</td></tr><tr><td>8. ¹³CH<sub>3</sub>D → ¹³CH<sub>2</sub>D + H</td><td>k<sub>8</sub> = 1/4 k<sub>f</sub> α<sub>8</sub></td><td>α<sub>8</sub> = 0.949</td></tr><tr><td>9. CH<sub>3</sub> + H → CH<sub>4</sub></td><td>k<sub>9</sub> = 1/2 k<sub>f</sub> α<sub>9</sub></td><td>α<sub>9</sub> = α<sub>1</sub></td></tr><tr><td>10. CH<sub>3</sub> + D → CH<sub>3</sub>D</td><td>k<sub>10</sub> = 1/2 k<sub>f</sub> α<sub>10</sub></td><td>α<sub>10</sub> = α<sub>2</sub> • Dαeq-p</td></tr><tr><td>11. CH<sub>2</sub>D + H → CH<sub>3</sub>D</td><td>k<sub>11</sub> = 1/2 k<sub>f</sub> α<sub>11</sub></td><td>α<sub>11</sub> = α<sub>3</sub> • Dαeq-s</td></tr><tr><td>12. CH<sub>2</sub>D + D → CH<sub>2</sub>D<sub>2</sub></td><td>k<sub>12</sub> = 1/2 k<sub>f</sub> α<sub>12</sub></td><td>α<sub>12</sub> = α<sub>4</sub> • D²γeq-p • Dαeq-p • Dαeq-s</td></tr><tr><td>13. CHD<sub>2</sub> + H → CH<sub>2</sub>D<sub>2</sub></td><td>k<sub>13</sub> = k<sub>f</sub> α<sub>13</sub></td><td>α<sub>13</sub> = α<sub>5</sub> • D²γeq-s • Dαeq-s • Dαeq-s</td></tr><tr><td>14. ¹³CH<sub>3</sub> + H → ¹³CH<sub>4</sub></td><td>k<sub>14</sub> = 1/2 k<sub>f</sub> α<sub>14</sub></td><td>α<sub>14</sub> = α<sub>6</sub> • ¹³αeq</td></tr><tr><td>15. ¹³CH<sub>3</sub> + D → ¹³CH<sub>3</sub>D</td><td>k<sub>15</sub> = 1/2 k<sub>f</sub> α<sub>15</sub></td><td>α<sub>15</sub> = α<sub>7</sub> • ¹³Dγeq-p • ¹³αeq • Dαeq-p</td></tr><tr><td>16. ¹³CH<sub>2</sub>D + H → ¹³CH<sub>3</sub>D</td><td>k<sub>16</sub> = k<sub>f</sub> α<sub>16</sub></td><td>α<sub>16</sub> = α<sub>8</sub> • ¹³Dγeq-s • ¹³αeq • Dαeq-s</td></tr></tbody></table>

be 1/4 for reactions (2) and (7), 3/4 for reactions (3) and (8), and 1/2 for reactions (4) and (5), as only 1/4 of the time the D is removed from ¹²CH₃D or ¹³CH₃D in its destruction while 3/4 of the time an H is removed. Accordingly, the secondary KIEs are calculated from the measured fractionation factor, Dα, and the branching ratios such that Dα = 1/4 α₂ + 3/4 α₃. We further consider bulk and clumped isotope equilibrium in the reversible reactions. At equilibrium, the ratio of the reverse and forward KIEs is equal to the EIE of the reaction. Thus, the KIEs for the reverse reactions are derived using the forward KIEs above and the EIEs calculated by Gropp et al. (2021). When the reaction is irreversible, only the forward KIEs are expressed, whereas the EIEs are expected to be more dominant when the reaction is highly reversible. In the latter case, we have allowed for exchange equilibrium among methane molecules to result from methane formation. This is done by including D²γeq and ¹³Dγeq, the fractionation factors representing equilibrium abundances of ¹²CH₂D₂ and ¹³CH₃D in thermodynamic equilibrium. In what follows we use an equilibration temperature for AOM of 20 °C. We verified the branching ratios in Table 2 using the full 280 reactions involving 50 isotopomer species and simple rate constants composed of square roots of inverse ratios of reduced masses (see below).

Reversibility is determined by the ratio Φ = k<sub>f</sub>/k<sub>f</sub>. Where Φ = 0, oxidation is “irreversible”, corresponding to Rayleigh fractionation with purely kinetic fractionation factors. We assume that the high-sulfate Svalbard incubation represents this irreversible end-member. Where Φ > 0, the influences of reversibility become immediately evident for even small values (and being mindful that k<sub>f</sub> refers to a 2nd order reaction while k<sub>f</sub> is for a first-order reaction, with the commensurate change in units for the rate constants). As Φ → 1, the oxidation at the Mcr step approaches full reversibility and equilibrium relative abundances of methane isotopologues result from the reaction with only negligible decreases in methane abundance. In all cases, the systems reach steady state, given sufficient time, with the steady-state abundance of methane, and its isotopologue composition, determined by the value for Φ.

### 4.3.2. Example applications

**Irreversible Rayleigh fractionation:** This example application of the model is trivial in so far as it is the model assumed to derive the fractionation factors. Nonetheless, it is useful to illustrate that the Svalbard sediment incubation data are properly accounted for using the model. Fig. 7A shows the path in Δ¹³CH₃D vs. Δ¹²CH₂D₂ space predicted by the model compared with the calibration data for the case of Rayleigh fractionation where Φ = 0. The model also fits the bulk isotope data well for similar values of F, the fraction of methane remaining (Fig. 8A).

**Reversibility:** Some data for methane known to have been processed by AOM exhibit significant variations in Δ¹³CH₃D at roughly constant Δ¹²CH₂D₂ and minimal variations in bulk isotope ratios. Our model explains these data as being the result of more reversibility in the Mcr-mediated transformation between CH₄ and CH₃ + H. We illustrate this here by comparing our results with the sediment slurry incubation data from the Santa Barbara Channel.

Most of the data exhibit a “horizontal” trend in the mass-18 isotopologue diagram that can be reproduced with our model using k<sub>f</sub>/k<sub>f</sub> ~ 5 × 10⁻⁶, resulting in a predicted steady-state concentration of methane that is about 99% of the initial value, and with negligible changes in δD and δ¹³C, in broad agreement with the data (Figs. 7B, 8B). As the k<sub>f</sub>/k<sub>f</sub> increases to 10⁻¹, the model leads to increases in both Δ¹³CH₃D and Δ¹²CH₂D₂ (Fig. 7B). This at least qualitative agreement between the AOM data from the Santa Barbara Channel and the model for significant reversibility in the first AOM step mediated by Mcr raises the specter of a potential monitor of Mcr-catalyzed intracellular isotope exchange operating at
<figure><img src="image_11.png" /><figcaption>**Fig. 7.** Methane clumped isotope data and model outputs of the general model for methane isotopologue fractionation due to AOM. (A) Comparison of our simplified model with incubations of Svalbard methane seep sediment slurry showing the model fits the data. (B) Comparison between Santa Barbara Channel methane seep incubations and our simplified model with inclusion of low turnover from CH₄ to CH₃ of about 1%, defined as $[CH₃]/[CH₄]$ at steady state. (C) Comparison between deep biosphere incubations and our simplified model. The *in-situ* methane gas from the Beatrix Gold Mine is shown for comparison (Young et al., 2017). The model in panels A–C uses 16 reactions involving 12 isotopologues and measured fractionation factors of the Svalbard incubation (Table 2). (D) Comparison of enzymatically-mediated exchange of methane molecules with our model prediction using the full 280 reactions involving 50 isotopomer species and simple rate constants composed of square roots of inverse ratios of reduced masses. Reversibility ($\Phi$) is determined by the ratio of the reverse rate constant to the forward rate constant ($k_f/k_j$). Where $\Phi$ is closer to one, oxidation is more reversible, and vice versa. The solid black line depicts theoretical thermodynamic equilibrium abundances of methane isotopologues, along with corresponding temperatures. Error bars are one standard error.</figcaption></figure>

near-threshold free energy geochemical settings, such as the deep biosphere (Fig. 7B–C).

Indeed, near-equilibrium methane isotopologue ratio of *in-situ* methane was found in the Beatrix Gold Mine with minor negative excursion in $\Delta^{12}CH_2D_2$ relative to equilibrium (Fig. 7C; Young et al., 2017). Here, methane was produced abiotically with negative $\Delta^{12}CH_2D_2$ values that are progressively erased by the incursion of biological activity that drives the isotopic bond ordering in methane towards equilibrium (Young et al., 2017). The board consistency in isotopologue ratios between the *in-situ* gas and laboratory incubations of the Beatrix samples further confirms the modeled trajectory (Fig. 7C).

We checked that the modeled trajectories are robust, and can be derived using kinetic rate constants defined by reduced masses alone combined with the full 280 reactions for all isomers (see below). We used the measured fractionation factors from the Svalbard incubation as model input for the Santa Barbara Channel and deep biosphere incubations as well. However, this may not be the case. To test whether the model is sensitive to the input fractionation factors, all KIEs and EIEs are set to 1.0 with the $^{D^2}\gamma_{eq}$ and $^{13D}\gamma_{eq}$ unchanged. This sensitivity test indicates that the modeled trajec-

tory remains the same and is not sensitive to the input fractionation factors when the reversibility is relatively high (Fig. S7).

*Exchange mediated by Mcr:* In the last example we make use of a related but different model to illustrate the path in $\Delta^{13}CH_3D$ vs. $\Delta^{12}CH_2D_2$ space due to enzymatically-mediated exchange of isotopes among methane molecules. The Mcr experiments provide the data against which to test this model. We do not use precisely the model shown in Table 2 in this case. The reason is that pure exchange is likely to be dominated by symmetry effects, and the fractionation factors obtained from the Svalbard experiments are an unresolvable convolution of these effects with other factors. To overcome this issue, this model is analogous to the simple model above, but in this case all isomers are included, producing 50 species related by 280 reactions, to explain the kinetics. Rather than the measured fractionation factors, we here use generic kinetic fractionation factors corresponding to square roots of the inverse ratios of relevant reduced masses ($^{12}C$-H, $^{13}C$-H, $^{12}C$-D, etc.).

We impose a value for $\Phi$ of >0.1, making the result dependent only on the relative number of isomers for each species and the two equilibrium fractionation factors $^{D^2}\gamma_{eq}$ and $^{13D}\gamma_{eq}$. Fig. 7D shows our prediction dominated by symmetry with the 60 °C
<figure><img src="image_12.png" /><figcaption>**Fig. 8.** Methane bulk isotope data and model outputs of the general model for methane isotopologue fractionation due to AOM. (A) Comparison of our simplified model with incubations of Svalbard methane seep sediment slurry showing the model fits the data. (B) Comparison between Santa Barbara Channel methane seep incubations and our simplified model. (C) Comparison between deep biosphere incubations and our simplified model. The model in panels A–C uses 16 reactions involving 12 isotopologues and measured fractionation factors of the Svalbard incubation (Table 2). (D) Comparison of enzymatically-mediated exchange of methane molecules with our model prediction using the full 280 reactions involving 50 isotopomer species and simple rate constants composed of square roots of inverse ratios of reduced masses. Reversibility ($\Phi$) is determined by the ratio of the reverse rate constant to the forward rate constant ($k_f/k_r$).</figcaption></figure>

exchange experiments in which the Mcr enzyme operates at its optimal condition. Here again, the agreement is reasonably good, suggesting that the inferences about the trajectory in $\Delta^{13}CH_3D$ vs. $\Delta^{12}CH_2D_2$ space due to nearly complete reversibility are correct. However, at 20 °C, below the optimal condition for the Mcr enzyme, a much lower $\Phi$ is needed to reproduce the data.

**Limitation:** Due to the complexity of isotopomer species involved in the reactions, the models presented here only focus on the first and, arguably, the most important step of AOM but do not include other intracellular pathways. As demonstrated by Wegener et al. (2021), a multistep isotope model for AOM including KIEs and EIEs is useful to explain the opposing bulk isotope effects (cf. Northrop, 1981). Such multistep isotopologue models have also been used to simulate clumped isotope effects of microbial methanogenesis (Cao et al., 2019; Young, 2019; Gropp et al., 2022; Ono et al., 2022). Indeed, this one-step model fails to yield inter-species isotope equilibrium for bulk isotopes with high reversibility, though equilibrium fractionation factors for both bulk and clumped isotopes are parameterized in the model. Therefore, we stress that the inputs and outputs of our model should be considered as illustrative approximations. Furthermore, the rate of AOM itself was not directly measured in our experiments other than the Svalbard slurry incubations. The lack of rate information limits the ability to quantify the degree of reversibility on the isotopologue ratios in our study. Nonetheless, the models presented

here allow multiple scenarios to be tested and suggest that high reversibility and low reversibility have characteristic trajectories in $\Delta^{13}CH_3D$ vs. $\Delta^{12}CH_2D_2$ space that are exhibited by the data.

### 4.4. Kinetically-driven AOM versus equilibrium-driven AOM

Taking all of our results together, we find that the $\Delta^{13}CH_3D$ and $\Delta^{12}CH_2D_2$ signatures of AOM can be understood in terms of reversibility at the Mcr step. The degree of reversibility is presumably controlled by the tenets of disequilibrium thermodynamics in which rates of reaction are proportional to reaction affinity (derivations in chemical potentials from equilibrium). Kinetic clumped isotope fractionation is found in the Svalbard slurry in which large amounts of methane are provided to the active ANME community. This kinetic signature is also seen in natural fluids from the South Chamorro Seamount, indicating that the kinetic isotopologue fractionation does not differ substantially between laboratory incubations and natural environments (Fig. 6). When less methane is provided, $\Delta^{13}CH_3D$ tends to increase towards thermodynamic near-equilibrium values with minimal increase in $\Delta^{12}CH_2D_2$ (Fig. 3E). For example, there is a striking difference in the isotopologue evolution of residual methane between the Svalbard slurry, in which kinetics drives $\Delta^{13}CH_3D$ and $\Delta^{12}CH_2D_2$ to extreme values, and the Santa Barbara slurry incubations in which changes occur mainly in $\Delta^{13}CH_3D$ and appear to approach equilibrium, even
<figure><img src="image_13.png" /><figcaption>**Fig. 9.** A schematic representation of methane clumped isotope signatures of diverse methane sources and sinks, as well as kinetically- and equilibrium-driven AOM. Circle symbols depict residual methane during incubations of the Svalbard methane seep sediment slurry, while methane data of natural fluids from the South Chamorro Seamount are shown as square symbols. Rhombus symbols illustrate residual methane from the Mcr exchange experiment and incubations of the Santa Barbara Channel sediment slurry. Deep biosphere incubations are shown as triangle symbols. The solid black line depicts theoretical thermodynamic equilibrium abundances of methane isotopologues, along with corresponding temperatures. The purple star represents the tank gas used in the slurry incubation experiments, and the zone of equilibrium-driven AOM is based on an initial gas of thermogenic origin. For comparison, two modeled trajectories for methane oxidation with OH radical (γ for ¹³CH₃D = 0.999, γ for ¹²CH₂D₂ = 0.907; Haghnegahdar et al., 2017), the major methane sink in air, are shown: the solid line is the Rayleigh distillation model (cf. Fig. 6A), the dashed line is the steady-state closed-system model ((E/k)/n₀ = 1; cf. Fig. 6E). The trajectories for aerobic oxidation of methane are similar to those of CH₄ + OH reaction (Krause et al., 2022). The zonation of methane sources is adopted from Young et al. (2017), Young (2019), Dong et al. (2021), etc. All the plotted data are from this study. The reader is referred to the text for interpretation of the figure.</figcaption></figure>

though the sediments are both from methane seeps (Fig. 3). The simple kinetic model shown here reproduces the salient features of this behavior as being the result of reversibility (Fig. 7). Another reason for this discrepancy is the biomass of ANME archaea under different conditions, such that in slurry that experienced high methane concentrations over a longer period of time, the biomass may be higher and vice-versa. We therefore infer that AOM with high sulfate concentrations may show either kinetic or equilibrium behavior depending on the local ANME community, the availability of methane, and thus the overall thermodynamic drive.

We conclude that both Δ¹³CH₃D and Δ¹²CH₂D₂ values reflect kinetic isotope fractionation and increase beyond intra-methane thermodynamic equilibrium when the net AOM reaction has low reversibility due to high rates, whereas largely reversible exchange between methane isotopologues facilitated by the Mcr enzyme under a lower thermodynamic drive leads to methane isotopic bond re-ordering, driving the abundances of methane isotopologues ultimately to thermodynamic quasi-equilibrium values. This observation is consistent with reaction affinity as a primary determinant for the behavior of the system. These departures from equilibrium, or reaction affinities, are controlled by environmental factors, including but not limited to, temperature, concentrations of methane and electron acceptors, and perhaps the type of electron acceptors. Although it seems that ¹³CH₃D and ¹²CH₂D₂ equilibration are ultimately controlled by the reversibility of AOM, the driver for different trajectories in Δ¹³CH₃D vs. Δ¹²CH₂D₂ space is still elusive. Future efforts should explicitly focus on the environmental controls of the enzymatic activity of intracellular pathways and the reversibility of AOM, and their intrinsic link to methane isotopologue ratios.

### 4.5. Implications

A schematic representation of methane clumped isotope signatures of diverse methane sources, as well as kinetically- and equilibrium-driven AOM, is shown in Fig. 9. Extremely low Δ¹²CH₂D₂ values down to ca. -40% have been observed for microbial methane under experimental and natural conditions, as well as thermogenic and abiotic methane generated by laboratory experiments (e.g., Young et al., 2017; Young, 2019; Taenzer et al., 2020; Dong et al., 2021). Besides equilibrium-driven AOM, isotopologue data near thermodynamic equilibrium have been attributed to high-maturity thermogenic or high-temperature abiotic methane, as well as low-rate methanogenesis (Stolper et al., 2015; Wang et al., 2015; Young et al., 2017; Ash et al., 2019; Giunta et al., 2019; Labidi et al., 2020; Warr et al., 2021; Xie et al., 2021; Gropp et al., 2022; Ono et al., 2022). Recently, Ono et al. (2022) demonstrated that methanogenesis does not occur or only proceeds at extremely slow rates at low pH₂ in energy-limited deep sedimentary environments using an isotopologue flow network model, and therefore inferred that near-equilibrium methane isotopologue signals in deep marine sediments are produced by the catalytic reversibility of the Mcr enzyme, likely from ANME archaea performing either AOM or net methanogenesis. This inference is consistent with our incubation experiments in which Mcr-catalyzed isotope exchange occurs.

To date, extremely positive Δ¹³CH₃D and Δ¹²CH₂D₂ values have been observed, as demonstrated here, only under conditions where AOM occurs with a low degree of reversibility. This observation raises an intriguing question whether methane clumped isotope signature of AOM could have the potential to be diagnostic of
AOM on Earth and perhaps throughout the solar system where methane has been detected (Webster et al., 2015; Waite et al., 2017; Yung et al., 2018; Thompson et al., 2022). For example, the in-situ discovery of isotopically depleted methane gas on Mars by Curiosity (Webster et al., 2015; House et al., 2022) indicates that methane is produced either biogenically or abiogenically (Atreya et al., 2007; Yung et al., 2018), and its emission results in the presence of methane in the atmosphere of Mars (Formisano et al., 2004). On the other hand, although the atmospheric methane sink on Mars has been quantified (Atreya et al., 2007; Yung et al., 2018), abiotic loss mechanisms would suggest a mean atmospheric residence time of ~300 years – substantially longer than observed methane lifetimes on the order of months to years (Lefèvre, 2019). The conditions for a subsurface microbial sink on Mars are still elusive. Nevertheless, the reaction transport model by Marlow et al. (2014) demonstrated that AOM could be a feasible metabolism on ancient Mars, for example, in sulfate-rich groundwater with methane from serpentinization by-products, and acid-sulfate fluids with methane produced from basalt alteration. The surface of Mars also includes deposits of iron oxides and manganese oxides, Fe- and Mn-dependent AOM therefore has also been proposed as a possible metabolism for recent and ancient Mars (House et al., 2011). Applying these inferences to our understanding of methane cycling on Mars is speculative. Nonetheless, a robust in-situ approach is needed to trace processes of methane cycling on future missions such as Mars Life Explorer (National Academies of Sciences, Engineering, and Medicine, 2022).

Mars is one of the most promising candidates for the application of methane isotope clumping. The reason is that doubly substituted isotopologues of methane remove the difficulties associated with using bulk ¹³C/¹²C and D/H elsewhere in the solar system, where the geochemical context necessary for interpreting these ratios is hampered by incompletely characterized chemical cycles (Lefèvre, 2019; Young, 2019; House et al., 2022). In Δ¹³CH₃D vs. Δ¹²CH₂D₂ space, microbial, thermogenic, and abiotic methane gases are either below or near the thermodynamic equilibrium predictions (Fig. 9). Once formed, methane could be oxidized photochemically or by microbial oxidation. These sinks have different Δ¹³CH₃D vs. Δ¹²CH₂D₂ signatures, depending upon whether a steady state is achieved or not. The details depend on the precise γ values, but in general, by analogy with Earth's atmosphere, oxidation in the atmosphere (e.g., by OH radicals) should ultimately lead to either high Δ¹²CH₂D₂ of order tens to hundreds of per mil at less variable Δ¹³CH₃D values in residual methane at steady state, or extremely low values for both where a steady state between methane formation and oxidation has not been achieved (Fig. 9) (Haghnegahdar et al., 2017). The trajectories for microbial aerobic oxidation of methane are like those of the CH₄ + OH reaction (Krause et al., 2022). AOM, on the other hand, can lead to extremely high Δ¹²CH₂D₂ and Δ¹³CH₃D values up to 65.1‰ and 19.5‰, respectively, as shown here. Indeed, these signatures are essentially orthogonal. Although mixing scenarios could mimic the extreme signatures from AOM, those mixing scenarios require extremely large differences in bulk carbon and hydrogen isotopic values of the two endmembers of methane (e.g., some 80‰ and 300‰ differences in endmember δ¹³C and δD values, respectively), possibly allowing kinetically-driven AOM to be distinguished from mixing with concurrent measurement of bulk isotopes.

While intriguing, we presently point out that clumped isotopes for methane would be an approach with a high potential for false negatives, as only anaerobic oxidation of methane, and only under conditions of low reversibility, would reasonably produce a sufficiently diagnostic signal. Another caveat here is that significant technical development would be necessary to measure methane that has an atmospheric mixing ratio on the order of a few ppb. Nonetheless, along with previous isotopologue studies on microbial

methanogenesis (Stolper et al., 2015; Wang et al., 2015; Young et al., 2017; Taenzer et al., 2020), our work on the effects of AOM on doubly substituted isotopologues of methane suggests that future missions could benefit from additional research and technique development for in-situ mission measurements of methane clumped isotopes on Mars, Enceladus, and other solar system bodies where the methane cycling is a key (bio)geochemical tracer.

## 5. Conclusions

We performed clumped isotopologue analysis on methane samples from laboratory microbial incubations and natural environments to elucidate the range of clumped isotope effects caused by AOM under different conditions. We found extremely high Δ¹³CH₃D and Δ¹²CH₂D₂ values up to 19.5‰ and 65.1‰, respectively, in microbial incubations of sediment slurries from the Svalbard methane seep and in natural fluids from the South Chamorro Seamount. We found near-equilibrium methane clumped isotope compositions in the Mcr exchange experiment and microbial incubations of sediment slurries and/or fracture fluids from the Santa Barbara Channel and deep biosphere. Thermodynamic drive and the reversibility of AOM intracellular reactions are invoked to interpret the observed differences in Δ¹³CH₃D and Δ¹²CH₂D₂ signatures and the interpretations are consistent with results from a simple, one-step isotopologue fractionation model. Low reversibility of the AOM reaction means a pronounced expression of kinetic fractionation of clumped isotopologues, whereas the Mcr-catalyzed intracellular isotope exchange may be more prominent with higher reversibility. The latter drives the isotopic bond ordering in methane towards equilibrium and may be responsible for the widely observed near-equilibrium clumped isotope signatures in substrate-limited anoxic natural environments where microbial methane production and/or oxidation occur.

## Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

We acknowledge C. Geoffrey Wheat and Jeffrey S. Seewald for providing methane samples from the Chamorro Seamount. Haolan Tang is acknowledged for Panorama Mass Spectrometer training, assistance, and consultation. We are grateful to Evert Duin for the gift of Mcr. We acknowledge Carolyn M. Harris for assistance with water hydrogen isotope analysis and David Yousavich for assistance with sulfate analysis. Formative discussions with Edwin A. Schauble, Jiawen Li, William D. Leavitt, Shuhei Ono, Ao Zhang, and Jingdong Gao helped shape this manuscript. We thank two anonymous reviewers for their helpful and constructive reviews of this paper. The research was funded by the NASA FINESS Fellowship 80NSSC21K1529 (to J. Liu and T. Treude) and the Alfred P. Sloan Foundation under the auspices of the Deep Carbon Observatory (to E.D. Young and R.L. Harris). We acknowledge additional funding for U. Toronto from the NSERC Discovery grant (to B. Sherwood Lolar) and the CIFAR Earth 4D program. R.L. Harris acknowledges financial support from the NASA Award 80NSSC19K1427. J.L. Ash was funded by a C-DEBI Postdoctoral Fellowship during this work.

## Appendix A. Supplementary material

The supplementary material includes supplementary figures S1–S7, supplementary methods for the ¹³CH₄ tracer microcosm
experiment at high hydrostatic pressure, and the research data.
Supplementary data to this article can be found online at
https://doi.org/10.1016/j.gca.2023.02.022.

References

Alperin, M.J., Reeburgh, W.S., Whiticar, M.J., 1988. Carbon and hydrogen isotope
fractionation resulting from anaerobic methane oxidation. Global Biogeochem.
Cycles 2, 279–288.

Ash, J.L., Egger, M., Treude, T., Kohl, I., Cragg, B., Parkes, R.J., Slomp, C.P., Sherwood, L. B., Young, E.D., 2019. Exchange catalysis during anaerobic methanotrophy revealed by ¹²CH₂D₂ and ¹³CH₃D in methane. Geochem. Persp. Lett. 10, 26–30.

Atekwana, E.A., Krishnamurthy, R.V., 1998. Seasonal variations of dissolved inorganic carbon and δ¹³C of surface waters: application of a modified gas evolution technique. J. Hydrol. 205, 265–278.

Atreya, S.K., Mahaffy, P.R., Wong, A.-S., 2007. Methane and related trace species on Mars: Origin, loss, implications for life, and habitability. Planet. Space Sci. 55, 358–369.

Beulig, F., Roy, H., McGlynn, S.E., Jørgensen, B.B., 2019. Cryptic CH₄ cycling in the sulfate-methane transition of marine sediments apparently mediated by ANME-1 archaea. ISME J. 13, 250–262.

Beulig, F., Schubert, F., Adhikari, R.R., Glombitza, C., Heuer, V.B., Hinrichs, K.U., Homola, K.L., Inagaki, F., Jørgensen, B.B., Kallmeyer, J., Krause, S.J.E., Morono, Y., Sauvage, J., Spivack, A.J., Treude, T., 2022. Rapid metabolism fosters microbial survival in the deep, hot subseafloor biosphere. Nat. Commun. 13, 312.

Bigeleisen, J., 1952. The Effects of Isotopic Substitution on the rates of Chemical Reactions. J. Phys. Chem. 56, 823–828.

Bigeleisen, J., 1955. Statistical Mechanics of Isotopic Systems with Small Quantum Corrections. I. General Considerations and the Rule of the Geometric Mean. J. Chem. Phys. 23, 2264–2267.

Bigeleisen, J., Mayer, M.G., 1947. Calculation of Equilibrium Constants for Isotopic Exchange Reactions. J. Chem. Phys. 15, 261–267.

Cao, X., Bao, H., Peng, Y., 2019. A kinetic model for isotopologue signatures of methane generated by biotic and abiotic CO₂ methanation. Geochim. Cosmochim. Acta 249, 59–75.

Chuang, P.-C., Yang, T.F., Wallmann, K., Matsumoto, R., Hu, C.-Y., Chen, H.-W., Lin, S., Sun, C.-H., Li, H.-C., Wang, Y., Dale, A.W., 2019. Carbon isotope exchange during anaerobic oxidation of methane (AOM) in sediments of the northeastern South China Sea. Geochim. Cosmochim. Acta 246, 138–155.

Cord-Ruwisch, R., 1985. A quick method for the determination of dissolved and precipitated sulfides in cultures of sulfate-reducing bacteria. J. Microbiol. Methods 4, 33–36.

Curtis, A.C., Wheat, C.G., Fryer, P., Moyer, C.L., 2013. Mariana Forearc Serpentinite Mud Volcanoes Harbor Novel Communities of Extremophilic Archaea. Geomicrobiol J. 30, 430–441.

Dale, A.W., Sommer, S., Lomnitz, U., Montes, I., Treude, T., Liebetrau, V., Gier, J., Hensen, C., Dengler, M., Stolpovsky, K., Bryant, L.D., Wallmann, K., 2015. Organic carbon production, mineralisation and preservation on the Peruvian margin. Biogeosciences 12, 1537–1559.

Dong, G., Xie, H., Formolo, M., Lawson, M., Sessions, A., Eiler, J., 2021. Clumped isotope effects of thermogenic methane formation: Insights from pyrolysis of hydrocarbons. Geochim. Cosmochim. Acta 303, 159–183.

Douglas, P.M.J., Stolper, D.A., Eiler, J.M., Sessions, A.L., Lawson, M., Shuai, Y., Bishop, A., Podlaha, O.G., Ferreira, A.A., Santos Neto, E.V., Niemann, M., Steen, A.S., Huang, L., Chimiak, L., Valentine, D.L., Fiebig, J., Luhmann, A.J., Seyfried, W.E., Etiope, G., Schoell, M., Inskeep, W.P., Moran, J.J., Kitchen, N., 2017. Methane clumped isotopes: Progress and potential for a new isotopic tracer. Org Geochim. 113, 262–282.

Douglas, P.M.J., Stolper, D.A., Smith, D.A., Walter Anthony, K.M., Paull, C.K., Dallimore, S., Wik, M., Crill, P.M., Winterdahl, M., Eiler, J.M., Sessions, A.L., 2016. Diverse origins of Arctic and Subarctic methane point source emissions identified with multiply-substituted isotopologues. Geochim. Cosmochim. Acta 188, 163–188.

Eldridge, D.L., Korol, R., Lloyd, M.K., Turner, A.C., Webb, M.A., Miller III, T.F., Stolper, D.A., 2019. Comparison of Experimental vs Theoretical Abundances of ¹³CH₃D and ¹²CH₂D₂ for Isotopically Equilibrated Systems from 1 to 500 °C. ACS Earth Space Chem. 3, 2747–2764.

Etiope, G., Sherwood Lollar, B., 2013. Abiotic methane on Earth. Rev. Geophys. 51, 276–299.

Formisano, V., Atreya, S., Encrenaz, T., Ignatiev, N., Giuranna, M., 2004. Detection of Methane in the Atmosphere of Mars. Science 306, 1758–1761.

Fryer, P.B., Salisbury, M.H., 2006. Leg 195 synthesis: Site 1200-Serpentinite seamounts of the Izu-Bonin/Mariana convergent plate margin (ODP Leg 125 and 195 drilling results). Proc. ODP Sci. Res. 195, 1–30.

Giunta, T., Young, E.D., Warr, O., Kohl, I., Ash, J.L., Martini, A., Mundle, S.O.C., Rumble, D., Pérez-Rodríguez, I., Wasley, M., LaRowe, D.E., Gilbert, A., Sherwood, L.B., 2019. Methane sources and sinks in continental sedimentary systems: New insights from paired clumped isotopologues ¹³CH₃D and ¹²CH₂D₂. Geochim. Cosmochim. Acta 245, 327–351.

Gonzalez, Y., Nelson, D.D., Shorter, J.H., McManus, J.B., Dyroff, C., Formolo, M., Wang, D.T., Western, C.M., Ono, S., 2019. Precise Measurements of ¹²CH₂D₂ by Tunable Infrared Laser Direct Absorption Spectroscopy. Anal. Chem. 91, 14967–14974.

Gropp, J., Iron, M.A., Halevy, I., 2021. Theoretical estimates of equilibrium carbon and hydrogen isotope effects in microbial methane production and anaerobic oxidation of methane. Geochim. Cosmochim. Acta 295, 237–264.

Gropp, J., Jin, Q., Halevy, I., 2022. Controls on the isotopic composition of microbial methane. Sci. Adv. 8, eabm5713.

Haghnegahdar, M.A., Schauble, E.A., Young, E.D., 2017. A model for ¹²CH₂D₂ and ¹³CH₃D as complementary tracers for the budget of atmospheric CH₄. Global Biogeochem. Cycles 31, 1387–1407.

Harris, R.L., Lau, M.C.Y., Cadar, A., Bartlett, D.H., Cason, E., Heerden, E., Onstott, T.C., Hotopp, J.C.D., 2018. Draft Genome Sequence of "Candidatus Bathydayaetha" Archaeon BE326-BA-RLH, an Uncultured Denitrifier and Putative Anaerobic Methanotroph from South Africa's Deep Continental Biosphere. Microbiol. Resour. Announce. 7, e01295-01218.

Harris, R.L., Vetter, M.C.Y.L., van Heerden, E., Cason, E., Vermeulen, J.-G., Taneja, A., Kieft, T.L., DeCoste, C.J., Laevsky, G.S., Onstott, T.C., 2021. FISH-TAMB, a Fixation-Free mRNA Fluorescent Labeling Technique to Target Transcriptionally Active Members in Microbial Communities. Microb. Ecol.

Hayes, J.M., 2001. Fractionation of Carbon and Hydrogen Isotopes in Biosynthetic Processes. Rev. Mineral. Geochem. 43, 225–277.

Heuer, V.B., Inagaki, F., Morono, Y., Kubo, Y., Maeda, L., the Expedition 370 Scientists, 2017. Temperature Limit of the Deep Biosphere off Muroto. Proceedings of the International Ocean Discovery Program, Volume 370.

Heuer, V.B., Inagaki, F., Morono, Y., Kubo, Y., Spivack, A.J., Viehweger, B., Treude, T., Beulig, F., Schubotz, F., Tonai, S., Bowden, S.A., Cramm, M., Henkel, S., Hirose, T., Homola, K., Hoshino, T., Ijiri, A., Imachi, H., Kamiya, N., Kaneko, M., Lagostina, L., Manners, H., McClelland, H.-L., Metcalfe, K., Okutsu, N., Pan, D., Raudsepp, M.J., Sauvage, J., Tsang, M.-Y., Wang, D.T., Whitaker, E., Yamamoto, Y., Yang, K., Maeda, L., Adhikari, R.R., Glombitza, C., Hamada, Y., Kallmeyer, J., Wendt, J., Wörmer, L., Yamada, Y., Kinoshita, M., Hinrichs, K.-U., 2020. Temperature limits to deep subseafloor life in the Nankai Trough subduction zone. Science 370, 1230–1234.

Hoehler, T.M., Borowski, W.S., Alperin, M.J., Rodriguez, N.M., Paull, C.K., 2000. Model, stable isotope, and radiotracer characterization of anaerobic methane oxidation in gas hydrate-bearing sediments of the Blake Ridge. Proc. Ocean Drilling Program, Sci. Results 164, 79–85.

Hoehler, T.M., Jørgensen, B.B., 2013. Microbial life under extreme energy limitation. Nat. Rev. Microbiol. 11, 83–94.

Holler, T., Wegener, G., Knittel, K., Boetius, A., Brunner, B., Kuypers, M.M.M., Widdel, F., 2009. Substantial ¹³C/¹²C and D/H fractionation during anaerobic oxidation of methane by marine consortia enriched in vitro. Environ. Microbiol. Rep. 1, 370–376.

Holler, T., Wegener, G., Niemann, H., Deusner, C., Ferdelman, T.G., Boetius, A., Brunner, B., Widdel, F., 2011. Carbon and sulfur back flux during anaerobic microbial oxidation of methane and coupled sulfate reduction. P. Natl. Acad. Sci. USA 108, E1484–E1490.

House, C.H., Beal, E.J., Orphan, V.J., 2011. The Apparent Involvement of ANMEs in Mineral Dependent Methane Oxidation, as an Analog for Possible Martian Methanotrophy. Life 1, 19–33.

House, C.H., Wong, G.M., Webster, C.R., Flesch, G.J., Franz, H.B., Stern, J.C., Pavlov, A., Atreya, S.K., Eigenbrode, J.L., Gilbert, A., Hofmann, A.E., Millan, M., Steele, A., Glavin, D.P., Malespin, C.A., Mahaffy, P.R., 2022. Depleted carbon isotope compositions observed at Gale crater, Mars. P. Natl. Acad. Sci. USA 119. e2115651119.

IPCC, 2021. Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change. IPCC, Geneva, Switzerland.

Jordan, S.F.A., Treude, T., Leifer, I., Janßen, R., Werner, J., Schulz-Vogt, H., Schmale, O., 2020. Bubble-mediated transport of benthic microorganisms into the water column: Identification of methanotrophs and implication of seepage intensity on transport efficiency. Sci. Rep. 10, 4682.

Joye, S.B., Boetius, A., Orcutt, B.N., Montoya, J.P., Schulz, H.N., Erickson, M.J., Lugo, S. K., 2004. The anaerobic oxidation of methane and sulfate reduction in sediments from Gulf of Mexico cold seeps. Chem. Geol. 205, 219–238.

Kallmeyer, J., Ferdelman, T.G., Weber, A., Fossing, H., Jørgensen, B.B., 2004. A cold chromium distillation procedure for radiolabeled sulfide applied to sulfate reduction measurements. Limnol. Oceanogr. Methods 2, 171–180.

Kawagucci, S., Miyazaki, J., Morono, Y., Seewald, J.S., Wheat, C.G., Takai, K., 2018. Cool, alkaline serpentinite formation fluid regime with scarce microbial habitability and possible abiotic synthesis beneath the South Chamorro Seamount. Progress in Earth and Planetary Science 5, 74.

Kovrikian, R.T., Callahan, S., Winstead, R., Lloyd, K.G., 2021. ANME-1 archaea may drive methane accumulation and removal in estuarine sediments. Environ. Microbiol. Rep. 13, 185–194.

Kopec, B.G., Feng, X., Posmentier, E.S., Sonder, L.J., 2019. Seasonal Deuterium Excess Variations of Precipitation at Summit, Greenland, and their Climatological Significance. J. Geophys. Res. Atmos. 124, 72–91.

Krause, S.J.E., Liu, J., Young, E.D., Treude, T., 2022. Δ¹³CH₃D and Δ¹²CH₂D₂ signatures of methane aerobically oxidized by Methyllosinus trichosporium with implications for deciphering the provenance of methane gases. Earth Planet. Sci. Lett. 593, 117681.

Krause, S.J.E., Treude, T., 2021. Deciphering cryptic methane cycling: Coupling of methylotrophic methanogenesis and anaerobic oxidation of methane in hypersaline coastal wetland sediment. Geochim. Cosmochim. Acta 302, 160–174.

Labidi, J., Young, E.D., Giunta, T., Kohl, I.E., Seewald, J., Tang, H., Lilley, M.D., Früh-Green, G.L., 2020. Methane thermometry in deep-sea hydrothermal systems:
Evidence for re-ordering of doubly-substituted isotopologues during fluid cooling. Geochim. Cosmochim. Acta 288, 248–261.

Lalk, E., Pape, T., Gruen, D.S., Kaul, N., Karolewski, J.S., Bohrmann, G., Ono, S., 2022. Clumped methane isotopologue-based temperature estimates for sources of methane in marine gas hydrates and associated vent gases. Geochim. Cosmochim. Acta 327, 276–297.

Laso-Pérez, R., Krukenberg, V., Musat, F., Wegener, G., 2018. Establishing anaerobic hydrocarbon-degrading enrichment cultures of microorganisms under strictly anoxic conditions. Nat. Protoc. 13, 1310–1330.

Lau, M.C.Y., Kieft, T.L., Kuloyo, O., Linage-Alvarez, B., van Heerden, E., Lindsay, M.R., Magnabosco, C., Wang, W., Wiggins, J.B., Guo, L., Perlman, D.H., Kyin, S., Shwe, H. H., Harris, R.L., Oh, Y., Yi, M.J., Purtschert, R., Slater, G.F., Ono, S., Wei, S., Li, H., Sherwood, L.B., Onstott, T.C., 2016. An oligotrophic deep-subsurface community dependent on syntrophy is dominated by sulfur-driven autotrophic denitrifiers. P. Natl. Acad. Sci. USA 113, E7927–E7936.

Lee, T.J., Martin, J.M.L., Taylor, P.R., 1995. An accurate ab initio quartic force field and vibrational frequencies for CH₄ and isotopomers. J. Chem. Phys. 102, 254–261.

Lefèvre, F., 2019. The Enigma of Methane on Mars: In: Cavalazzi, B., Westall, F. (Eds.), Biosignatures for Astrobiology. Springer International Publishing, Cham, pp. 253–266.

Liu, Q., Liu, Y., 2016. Clumped-isotope signatures at equilibrium of CH₄, NH₃, H₂O, H₂S and SO₂. Geochim. Cosmochim. Acta 175, 252–270.

Ma, Q., Wu, S., Tang, Y., 2008. Formation and abundance of doubly-substituted methane isotopologues (¹³CH₃D) in natural gas systems. Geochim. Cosmochim. Acta 72, 5446–5456.

Mahler, F., Grabarbe, W., Kahnt, J., Thauer, R.K., Duin, E.C., 2002. The nickel enzyme methyl-coenzyme M reductase from methanogenic archaea: in vitro interconversions among the EPR detectable MCR-red1 and MCR-red2 states. J. Biol. Inorg. Chem. 7, 101–112.

Mariotti, A., Germon, J.C., Hubert, P., Kaiser, P., Letolle, R., Tardieu, A., Tardieu, P., 1981. Experimental determination of nitrogen kinetic isotope fractionation: Some principles; illustration for the denitrification and nitrification processes. Plant Soil 62, 413–430.

Marlow, J.J., LaRowe, D.E., Ehlmann, B.L., Amend, J.P., Orphan, V.J., 2014. The Potential for Biologically Catalyzed Anaerobic Methane Oxidation on Ancient Mars. Astrobiology 14, 292–307.

Melaniuk, K., Sztybor, K., Treude, T., Sommer, S., Rasmussen, T.L., 2022. Influence of methane seepage on isotopic signatures in living deep-sea benthic foraminifera, 79° N. Sci. Rep. 12, 1169.

National Academies of Sciences, Engineering, and Medicine. 2022. Origins, Worlds, and Life: A Decadal Strategy for Planetary Science and Astrobiology 2023–2032. The National Academies Press, Washington, DC.

Northrop, D.B., 1981. The Expression of Isotope Effects on Enzyme-Catalyzed Reactions. Annu. Rev. Biochem. 50, 103–131.

Ono, S., Rhim, J.H., Gruen, D.S., Taubner, H., Kölling, M., Wegener, G., 2021. Clumped isotopologue fractionation by microbial cultures performing the anaerobic oxidation of methane. Geochim. Cosmochim. Acta 293, 70–85.

Ono, S., Rhim, J.H., Ryberg, E.C., 2022. Rate limits and isotopologue fractionations for microbial methanogenesis examined with combined pathway protein cost and isotopologue flow network models. Geochim. Cosmochim. Acta 325, 296–315.

Ono, S., Wang, D.T., Gruen, D.S., Sherwood, L.B., Zahniser, M.S., McManus, B.J., Nelson, D.D., 2014. Measurement of a Doubly Substituted Methane Isotopologue, ¹³CH₃D, by Tunable Infrared Laser Direct Absorption Spectroscopy. Anal. Chem. 86, 6487–6494.

Reeburg, W.S., 2007. Global Methane Biogeochemistry. In: Holland, H.D., Turekian, K.K. (Eds.), Treatise on Geochemistry. Pergamon, Oxford, pp. 1–32.

Repeta, D.J., Ferrón, S., Sosa, O.A., Johnson, G.C., Repeta, L.D., Acker, M., DeLong, E.F., Karl, D.M., 2016. Marine methane paradox explained by bacterial degradation of dissolved organic matter. Nat. Geosci. 9, 884–887.

Röckmann, T., Popa, M.E., Krol, M.C., Hofmann, M.E.G., 2016. Statistical clumped isotope signatures. Scientific Reports 6, 31947.

Scheller, S., Goenrich, M., Boecher, R., Thauer, R.K., Jaun, B., 2010. The key nickel enzyme of methanogenesis catalyses the anaerobic oxidation of methane. Nature 465, 606–608.

Scheller, S., Goenrich, M., Thauer, R.K., Jaun, B., 2013. Methyl-Coenzyme M Reductase from Methanogenic Archaea: Isotope Effects on the Formation and Anaerobic Oxidation of Methane. J. Am. Chem. Soc. 135, 14975–14984.

Schoell, M., 1983. Genetic Characterization of Natural Gases. AAPG Bull. 67, 2225–2238.

Seewald, J.S., Doherty, K.W., Hammar, T.R., Liberatore, S.P., 2002. A new gas-tight isobaric sampler for hydrothermal fluids. Deep-Sea Res. Pt. I 49, 189–196.

Sherwood Lollar, B., Lacrampe-Couloume, G., Slater, G.F., Ward, J., Moser, D.P., Ghirring, T.M., Lin, L.H., Onstott, T.C., 2006. Unravelling abiogenic and biogenic sources of methane in the Earth's deep subsurface. Chem. Geol. 226, 328–339.

Stolper, D.A., Lawson, M., Davis, C.L., Ferreira, A.A., Neto, E.V.S., Ellis, G.S., Lewan, M. D., Martini, A.M., Tang, Y., Schoell, M., Sessions, A.L., Eiler, J.M., 2014a. Formation temperatures of thermogenic and biogenic methane. Science 344, 1500–1503.

Stolper, D.A., Martini, A.M., Clog, M., Douglas, P.M., Shusta, S.S., Valentine, D.L., Sessions, A.L., Eiler, J.M., 2015. Distinguishing and understanding thermogenic and biogenic sources of methane using multiply substituted isotopologues. Geochim. Cosmochim. Acta 161, 219–247.

Stolper, D.A., Sessions, A.L., Ferreira, A.A., Santos Neto, E.V., Schimmelmann, A., Shusta, S.S., Valentine, D.L., Eiler, J.M., 2014b. Combined ¹³C-D and D-D clumping in methane: Methods and preliminary results. Geochim. Cosmochim. Acta 126, 169–191.

Taenzer, L., Labidi, J., Masterson, A.L., Feng, X., Rumble, D., Young, E.D., Leavitt, W.D., 2020. Low Δ¹²CH₂D₂ values in microbalgentic methane result from combinatorial isotope effects. Geochim. Cosmochim. Acta 285, 225–236.

Takai, K., Moyer, C.L., Miyazaki, M., Nogi, Y., Hirayama, H., Nealson, K.H., Horikoshi, K., 2005. Marinobacter alkaliphilus sp. nov., a novel alkaliphilic bacterium isolated from subseafloor alkaline serpentine mud from Ocean Drilling Program Site 1200 at South Chamorro Seamount, Mariana Forearc. extremophiles 9, 17–27.

Thiagarajan, N., Kitchen, N., Xie, H., Ponton, C., Lawson, M., Formolo, M., Eiler, J., 2020. Identifying thermogenic and microbial methane in deep water Gulf of Mexico Reservoirs. Geochim. Cosmochim. Acta 275, 188–208.

Thompson, M.A., Krissansen-Totton, J., Wogan, N., Telus, M., Fortney, J.J., 2022. The case and context for atmospheric methane as an exoplanet biosignature. P. Natl. Acad. Sci. USA 119.

Timmers, P.H.A., Welte, C.U., Koehorst, J.J., Plugge, C.M., Jetten, M.S.M., Stams, A.J.M., 2017. Reverse Methanogenesis and Respiration in Methanotrophic Archaea. Archaea 2017, 22.

Treude, T., Krüger, M., Boetius, A., Jørgensen, B.B., 2005. Environmental control on anaerobic oxidation of methane in the gassy sediments of Eckernförde Bay (German Baltic). Limnol. Oceanogr. 50, 1771–1786.

Treude, T., Orphan, V., Knittel, K., Gieseke, A., House, C.H., Boetius, A., 2007. Consumption of Methane and CO₂ by Methanotrophic Microbial Mats from Gas Seeps of the Anoxic Black Sea. Appl. Environ. Microbiol. 73, 2271–2283.

Treude, T., Ziebis, W., 2010. Methane oxidation in permeable sediments at hydrocarbon seeps in the Santa Barbara Channel, California. Biogeosciences 7, 3095–3108.

Turner, A.C., Korol, R., Eldridge, D.L., Bill, M., Conrad, M.E., Miller, T.F., Stolper, D.A., 2021. Experimental and theoretical determinations of hydrogen isotopic equilibrium in the system CH₄–H₂–H₂O from 3 to 200 °C. Geochim. Cosmochim. Acta 314, 223–269.

Tyne, R.L., Barry, P.H., Lawson, M., Byrne, D.J., Warr, O., Xie, H., Hillegonds, D.J., 2021. Rapid microbial methanogenesis during CO₂ storage in hydrocarbon reservoirs. Nature 600, 670–674.

Urey, H.C., 1947. The thermodynamic properties of isotopic substances. J. Chem. Soc. (Resumed), 562–581.

Waite, J.H., Glein, C.R., Perryman, R.S., Teolis, B.D., Magee, B.A., Miller, G., Grimes, J., Perry, M.E., Miller, K.E., Bouquet, A., Lunine, J.I., Brockwell, T., Bolton, S.J., 2017. Cassini finds molecular hydrogen in the Enceladus plume: Evidence for hydrothermal processes. Science 356, 155–159.

Wang, D.T., Gruen, D.S., Lollar, B.S., Hinrichs, K.-U., Stewart, L.C., Holden, J.F., Hristov, A.N., Pohlman, J.W., Morrill, P.L., Könneke, M., Novelche, K.B., Reeves, E.P., Sutcliffe, C.N., Ritter, D.J., Seewald, J.S., McIntosh, J.C., Hemond, H.F., Kubo, M.D., Cardace, D., Hoehler, T.M., Ono, S., 2015. Nonequilibrium clumped isotope signals in microbial methane. Science 348, 428–431.

Wang, D.T., Welander, P.V., Ono, S., 2016. Fractionation of the methane isotopologues ¹³CH₄, ¹²CH₃D, and ¹³CH₃D during aerobic oxidation of methane by Methylococcus capsulatus (Bath). Geochim. Cosmochim. Acta 192, 186–202.

Warr, O., Young, E.D., Giunta, T., Kohl, I.E., Ash, J.L., Sherwood, L.B., 2021. High-resolution, long-term isotopic and isotopologue variation identifies the sources and sinks of methane in a deep subsurface carbon cycle. Geochim. Cosmochim. Acta 294, 315–334.

Webb, M.A., Miller, T.F., 2014. Position-Specific and Clumped Stable Isotope Studies: Comparison of the Urey and Path-Integral Approaches for Carbon Dioxide, Nitrous Oxide, Methane, and Propane. Chem. A Eur. J. 118, 467–474.

Webster, C.R., Mahaffy, P.R., Atreya, S.K., Flesch, G.J., Mischna, M.A., Meslin, P.-Y., Farley, K.A., Conrad, P.G., Christensen, L.E., Pavlov, A.A., Martín-Torres, J., Zorzano, M.-P., McConnochie, T.H., Owen, T., Eigenbrode, J.L., Glavin, D.P., Steele, A., Malespin, C.A., Archer, P.D., Sutter, B., Coll, P., Freissinet, C., McKay, C.P., Moores, J.E., Schwenzer, S.P., Bridges, J.C., Navarro-Gonzalez, R., Gellert, R., Lemmon, M.T., 2015. Mars methane detection and variability at Gale crater. Science 347, 415–417.

Wegener, G., Gropp, J., Taubner, H., Halevy, I., Elvert, M., 2021. Sulfate-dependent reversibility of intracellular reactions explains the opposing isotope effects in the anaerobic oxidation of methane. Sci. Adv. 7, eabe4939.

Wheat, C.G., Fryer, P., Fisher, A.T., Hulme, S., Jannasch, H., Mottl, M.J., Becker, K., 2008. Borehole observations of fluid flow from South Chamorro Seamount, an active serpentinite mud volcano in the Mariana forearc. Earth Planet. Sci. Lett. 267, 401–409.

Wheat, C.G., Seewald, J.S., Takai, K., 2020. Fluid transport and reaction processes within a serpentinite mud volcano: South Chamorro Seamount. Geochim. Cosmochim. Acta 269, 413–428.

Whitehill, A.R., Joelsson, L.M.T., Schmidt, J.A., Wang, D.T., Johnson, M.S., Ono, S., 2017. Clumped isotope effects during OH and Cl oxidation of methane. Geochim. Cosmochim. Acta 196, 307–325.

Whiticar, M.J., 1999. Carbon and hydrogen isotope systematics of bacterial formation and oxidation of methane. Chem. Geol. 161, 291–314.

Widdel, F., Bak, F., 1992. Gram-Negative Mesophilic Sulfate-Reducing Bacteria: In: Balows, A., Trüper, H.G., Dworkin, M., Harder, W., Schleifer, K.-H. (Eds.), The Prokaryotes: A Handbook on the Biology of Bacteria: Ecolphysiology, Isolation, Identification, Applications. Springer New York, New York, NY, pp. 3352–3378.

Xiao, K.-Q., Beulig, F., Kjeldsen, K.U., Jørgensen, B.B., Risgaard-Petersen, N., 2017. Concurrent Methane Production and Oxidation in Surface Sediment from Aarhus Bay, Denmark. Front. Microbiol. 8.
Xie, H., Dong, G., Formolo, M., Lawson, M., Liu, J., Cong, F., Mangenot, X., Shuai, Y., Ponton, C., Eiler, J., 2021. The evolution of intra- and inter-molecular isotope equilibria in natural gases with thermal maturation. *Geochim. Cosmochim. Acta* 307, 22–41.

Yamamoto, S., Alcauskas, J.B., Crozier, T.E., 1976. Solubility of methane in distilled water and seawater. *J. Chem. Eng. Data* 21, 78–80.

Yan, Z., Joshi, P., Gorski, C.A., Ferry, J.G., 2018. A biochemical framework for anaerobic oxidation of methane driven by Fe (III)-dependent respiration. *Nat. commun.* 9.

Yeung, L.Y., 2016. Combinatorial effects on clumped isotopes and their significance in biogeochemistry. *Geochim. Cosmochim. Acta* 172, 22–38.

York, D., Evensen, N.M., Martínez, M.L., Delgado, J.D.B., 2004. Unified equations for the slope, intercept, and standard errors of the best straight line. *American Journal of Physics* 72, 367–375.

Yoshinaga, M.Y., Holler, T., Goldhammer, T., Wegener, G., Pohlman, J.W., Brunner, B., Kuypers, M.M.M., Hinrichs, K.-U., Elvert, M., 2014. Carbon isotope equilibration during sulphate-limited anaerobic oxidation of methane. *Nat. Geosci.* 7, 190–194.

Young, E.D., 2019. A Two-Dimensional Perspective on CH₄ Isotope Clumping: Distinguishing Process from Source. In: Orcutt, B.N., Daniel, I., Dasgupta, R. (Eds.), Deep Carbon: Past to Present. Cambridge University Press, Cambridge, pp. 388–414.

Young, E.D., Kohl, I.E., Lollar, B.S., Etiope, G., Rumble, D., Li, S., Haghnegahdar, M.A., Schauble, E.A., McCain, K.A., Foustoukos, D.I., Sutcliffe, C., Warr, O., Ballentine, C. J., Onstott, T.C., Hosgormez, H., Neubeck, A., Marques, J.M., Pérez-Rodríguez, I., Rowe, A.R., LaRowe, D.E., Magnabosco, C., Yeung, L.Y., Ash, J.L., Bryndzia, L.T.,

2017. The relative abundances of resolved ¹²CH₂D₂ and ¹³CH₃D and mechanisms controlling isotopic bond ordering in abiotic and biotic methane gases. *Geochim. Cosmochim. Acta* 203, 235–264.

Young, E.D., Rumble, D., Freedman, P., Mills, M., 2016. A large-radius high-mass-resolution multiple-collector isotope ratio mass spectrometer for analysis of rare isotopologues of O₂, N₂, CH₄ and other gases. *Int. J. Mass Spectrom.* 401, 1–10.

Yu, L., He, D., Yang, L., Rensing, C., Zeng, R.J., Zhou, S., 2022. Anaerobic methane oxidation coupled to ferrihydrite reduction by Methanosarcina barkeri. *Sci. Total Environ.* 844, 157235.

Yung, Y.L., Chen, P., Nealson, K., Atreya, S., Beckett, P., Blank, J.G., Ehlmann, B., Eiler, J., Etiope, G., Ferry, J.G., Forget, F., Gao, P., Hu, R., Kleinböhl, A., Klusman, R., Lefèvre, F., Miller, C., Mischna, M., Mumma, M., Newman, S., Oehler, D., Okumura, M., Oremland, R., Orphan, V., Popa, R., Russell, M., Shen, L., Sherwood, L.B., Staehle, R., Stamenković, V., Stolper, D., Templeton, A., Vandaele, A.C., Viscardy, S., Webster, C.R., Wennberg, P.O., Wong, M.L., Worden, J., 2018. Methane on Mars and Habitability: Challenges and Responses. *Astrobiology* 18, 1221–1242.

Zhang, J., Quay, P.D., Wilbur, D.O., 1995. Carbon isotope fractionation during gas-water exchange and dissolution of CO₂. *Geochim. Cosmochim. Acta* 59, 107–114.

Zhang, N., Snyder, G.T., Lin, M., Nakagawa, M., Gilbert, A., Yoshida, N., Matsumoto, R., Sekine, Y., 2021. Doubly substituted isotopologues of methane hydrate (¹³CH₃D and ¹²CH₂D₂): Implications for methane clumped isotope effects, source apportionments and global hydrate reservoirs. *Geochim. Cosmochim. Acta* 315, 127–151.