The geothermal gradient shapes microbial diversity and processes in natural-gas-bearing sedimentary aquifers

Deep subsurface microorganisms constitute over 80 % of Earth's prokaryotic biomass and play an important role in global biogeochemical cycles. Geochemical processes driven by geothermal heating are key factors influencing their biomass and activities, yet their full breadth remains uncaptured. Here, we investigated the microbial community composition and metabolism in microbial-natural-gas-bearing aquifers at temperatures ranging from 38 to 81 °C, situated above nonmicrobial-gas- and oil-bearing sediments at temperatures exceeding 90 °C. Cultivation-based and molecular gene analyses, including radiotracer measurements, of formation water indicated variations in predominant methanogenic pathways across different temperature regimes of upper aquifers: high potential for hydrogenotrophic–methylotrophic, hydrogenotrophic, and acetoclastic methanogenesis at temperatures of 38, 51–65, and 73–81 °C, respectively. The potential for acetoclastic methanogenesis correlated with elevated acetate concentrations with increasing depth, possibly due to the decomposition of sedimentary organic matter. In addition to acetoclastic methanogenesis, in aquifers at temperatures as high as or higher than 65 °C, acetate is potentially utilized by microorganisms responsible for the dissimilatory reduction of sulfur compounds other than sulfate because of their high relative abundance at greater depths. The stable sulfur isotopic analysis of sulfur compounds in water and oil samples suggested that hydrogen sulfide, generated through the thermal decomposition of sulfur compounds in oil, migrates upward and is subsequently oxidized with iron oxides present in sediments, yielding elemental sulfur and thiosulfate. These compounds are consumed by sulfur-reducing microorganisms, possibly reflecting elevated microbial populations in aquifers at temperatures as high as or higher than 73 °C. These findings reveal previously overlooked geothermal-heat-driven geochemical and microbiological processes involved in carbon and sulfur cycling in the deep sedimentary biosphere.

1 Introduction

The deep subsurface environment harbors a substantial fraction of Earth's prokaryotes (McMahon and Parnell, 2014; Magnabosco et al., 2018), constituting over 80 % of the total prokaryotic biomass (Bar-On et al., 2018). The metabolic activities of these microorganisms play a pivotal role in global biogeochemical cycling, such as carbon, nitrogen, and sulfur (Magnabosco et al., 2018; Aloisi et al., 2006). A prominent example is methane production. Much of the methane hydrate, the largest methane reservoir on Earth, is suggested to be generated by these subsurface microorganisms (Kvenvolden, 1995). Additionally, coal-bed methane and shale gas reservoirs may also contain a significant amount of microbially derived methane (Vinson et al., 2017).

In the absence of light energy, microorganisms inhabiting deep sedimentary environments rely on chemical energy derived from the oxidation of reduced substances in sediments, with organic matter oxidation being particularly critical (Lovley and Chapelle, 1995; Jørgensen and Boetius, 2007; Arndt et al., 2013). The labile components are consumed during burial, limiting the availability of energy sources for microorganisms as sediment age increases (Middelburg, 1989). However, substantial populations of active microbial cells have been observed even in deep buried sediments older than 16×10⁶ years (Schippers et al., 2005). It has been hypothesized that temperature increase during burial stimulates thermal or biological degradation of recalcitrant organic matter, possibly sustaining microbial activities (Parkes et al., 2000). To address the fundamental questions of how microbial cells in the deep biosphere can survive with limited energy sources, this hypothesis has been examined through numerical simulations (Horsfield et al., 2006), laboratory incubation experiments (Parkes et al., 2007), and geochemical analysis (Malinverno and Martinez, 2015) of subseafloor sediments. Thus, while temperature increase during burial is posited to drive subsurface microbial activity, field-observation-based microbiological research remains limited, and the mechanisms are not fully understood.

As sediment compaction progresses with burial, pore size and permeability decrease, reducing the living space, available water, and nutrients, thereby inhibiting microbial growth (Fredrickson et al., 1997). To investigate the impact of temperature increase on the subsurface microbial ecosystem, it is essential to target subsurface environments where these inhibitory factors are minimized. Therefore, we focused on aquifers, which, even at great depths, maintain high porosity and permeability, providing ample living space and water for microorganisms (Fredrickson et al., 1997; McMahon and Chapelle, 1991; Lovley and Chapelle, 1995; Krumholz et al., 1997).

In this study, we targeted a gas field in central Japan, where aquifers exhibit a wide temperature range, spanning approximately 35 to 80 °C, due to a steep geothermal gradient (5 °C per 100 m) (Kato, 2018). We collected formation water (FW) from each aquifer and employed a comprehensive approach, including geochemical analysis, radiotracer measurements, molecular gene analysis, and cultivation experiments to evaluate microbial diversity, community structure, and potential metabolic processes. Furthermore, in this field, high-temperature oil deposits are situated deeper than the series of aquifers. We therefore also examined the impact of oil components on microorganisms in the upper aquifers. The aim of this study is to elucidate the effects of temperature rise and associated geochemical processes, such as the decomposition of sedimentary organic matter and petroleum formation, on microbial diversity, community structure, and metabolic processes, including the conversion of carbon and sulfur compounds.

2 Materials and methods

2.1 Site description and sample collection

The study site is located along the Sea of Japan coast in the city of Tainai, Niigata Prefecture, Japan (Fig. 1a). The reservoirs are distributed in the Neogene and Quaternary formations of the Niigata sedimentary basin. The reservoirs of dissolved natural gases are distributed within the Pliocene to Pleistocene Nishiyama and Pleistocene Haizume Formations and mainly consist of sandstone or conglomerate (Kobayashi, 2000; Fukano et al., 2023). These sediments are interpreted as turbidite sediments deposited under submarine fan to shelf conditions (Takano et al., 2001). Logging analysis of the gas production wells revealed that each reservoir layer was separated by silt–mudstone, and their horizontal continuity was confirmed. Oil and associated gas are deposited along anticlinal structures in the lower part of the Nishiyama Formation and the Pliocene Shiiya Formation (Fukano et al., 2023). The reservoirs of oil deposits consist of sandstone in turbidite sediments deposited under submarine fan and shelf conditions (Takano et al., 2001). The chemical and isotopic compositions of natural gases (Kaneko and Igari, 2020) indicate that methane dissolved in upper aquifers is primarily of microbial origin, whereas that from lower oil deposits primarily originates from oil-associated thermogenic processes (Fig. S1 in the Supplement) based on the classification by Milkov and Etiope (Milkov and Etiope, 2018). In this gas field, gases are dissolved in FW and produced for commercial purposes by pumping gas-associated FW from upper aquifers. Crude oil and gases are also collected from lower oil deposits.

Figure 1 Location of the wells from which water samples were collected for microbiological analyses (a), the depths of the wells shown on the geological cross-section (b), and the measured water temperatures related to those depths (c). The ratio of the vertical to horizontal scale of the cross-section is 2.1:1 (b). The range of the well strainers is indicated by a red bar (b, c), and the average depth (m b.g.s.) of each well strainer is also indicated (b). Because three wells are present at the same location, only one well is shown (b).

The FW samples subjected to microbiological analyses were collected from the separate reservoir layers in the Haizume and Nishiyama Formations through six commercial production wells with depths ranging from approximately 500 to 1430 m b.g.s. (meters below the ground surface; Fig. 1a and b). Crude oil and associated FW samples were collected from the Shiiya Formation through tanks separating oil, FW, and gas produced from multiple wells and used for microbial cell counting and geochemical analysis. The geological information and analytical purpose of each sample are summarized in Table S1 in the Supplement.

FW samples for methanogen activity measurements and cultivation were collected in sterilized glass bottles with butyl rubber stoppers and screw caps. The bottles were purged with N₂ gas before and during sample collection and then filled with FW to prevent the penetration of air. For the molecular analysis, 2 L FW samples were collected and filtered through a 0.2 µm pore size Millipore Express Plus membrane filter (Millipore, Billerica, MA, USA) to harvest microbial cells, which were subsequently stored at −20 °C. The samples used for total cell counts were fixed with formalin at a final concentration of 2 % ($v/v$) immediately after sampling and stored at 4 °C. Zinc acetate solution (2.0 M) was added to the water sample to fix the hydrogen sulfide dissolved in the FW immediately after collection, and the solution was used for the measurement of the hydrogen sulfide concentration.

2.2 Geochemical analysis

The inorganic chemical and total organic carbon contents in the FW samples were measured as described previously (Katayama et al., 2015), and the analytical methods used for each compound are summarized in Table S2. The concentrations of formate, acetate, propionate, and butyrate were measured using a high-performance liquid chromatograph system (Prominence HPLC, Shimadzu Co. Ltd., Kyoto, Japan) equipped with an electrical conductivity detector. The concentration of the fixed hydrogen sulfide sample was measured colorimetrically with a SmartSpecTM Plus spectrophotometer (Bio-Rad, CA, USA) using ferrous hydrogen sulfide reagents (HACH; Loveland, CO, USA). Hydrocarbon components in FW samples from upper aquifers were extracted and concentrated with dichloromethane and qualitatively analyzed using a gas chromatograph–mass spectrometer (5973 GC–MS; Agilent Technologies, CA). The crude oil sample from the oil deposits was diluted with hexane and qualitatively analyzed via a GC–MS. The concentration of phenols in the FW sample was measured via the 4-aminoantipyrine spectrophotometric method. The concentrations of methanol, toluene, and xylene in the FW sample were measured using a GC–MS (Agilent 5973; Agilent Technologies) instrument equipped with a headspace sampler (Agilent 7697A, Agilent Technologies). The number of replication and the standard deviation value for each measurement are provided in Table S3.

The stable oxygen and hydrogen isotopic compositions of the FW samples were measured via wavelength-scanned cavity ring-down spectroscopy (WS-CRDS) via a laser spectroscopy analyzer (Picarro L2120-i). Each isotopic ratio was measured six times per sample, with standard deviations of δD and δ¹⁸O less than 2.0 ‰ and 0.5 ‰, respectively.

Stable sulfur isotope analysis of sulfate and hydrogen sulfide dissolved in FW was performed as described previously (Yamanaka et al., 2013; Katayama et al., 2019). For stable isotope analysis of elemental sulfur (dissolved in FW), hydrogen sulfide dissolved in FW samples was first removed by precipitation using zinc acetate solution, followed by filtration using filter paper to trap the generated zinc sulfide. Liquid–liquid extraction of elemental sulfur was performed from the filtrate with carbon disulfide followed by evaporation of carbon disulfide in a sand bath (∼60 °C), and crystallized elemental sulfur was measured by an elemental analyzer–isotope ratio mass spectrometer (IsoPrime EA; GV Instruments, Manchester, UK). For stable sulfur isotope analysis of sulfur compounds in oil, the compounds were oxidized and decomposed in a high-pressure oxygen combustion unit (Parr Bomb), precipitated as barium sulfate, wrapped in a tin capsule containing V₂O₅, and measured using the elemental analyzer–isotope ratio mass spectrometer. Each sulfur isotope ratio represents the average of duplicate measurements, with deviations less than 0.4 ‰.

2.3 Methanogenic activity measurement using a radiotracer

The methane production rates in the FW samples were measured using radiotracer experiments as described previously (Katayama et al., 2015). Briefly, a 20 mL water sample in a 50 mL serum vial sealed with a butyl rubber stopper and aluminum crimp was injected with one of the radiotracers: ¹⁴C-bicarbonate (10 mL, 199 kBq), [2-¹⁴C]-acetate (10 mL, 81 kBq), or ¹⁴C-methanol. The vials were incubated under an atmosphere of ${\mathrm{N}}_{\mathrm{2}}/{\mathrm{CO}}_{\mathrm{2}}$ (80:20, $v/v$) for 14, 28, or 42 d. The activity measurements were conducted in triplicate for each of the three cultivation periods and for each of the three radiotracers. The incubation temperatures were set close to the measured water temperature of each sample (i.e., 559 m b.g.s., 37 °C; 754 m b.g.s., 50 °C; 845 m b.g.s., 55 °C; 1049 m b.g.s., 65 °C; 1115 m b.g.s., 75 °C; 1373 m b.g.s., 80 °C). The produced ¹⁴CH₄ was oxidized to ¹⁴CO₂ and measured using a Tri-Carb 3100TR liquid scintillation counter (Perkin Elmer). The ¹⁴C activity at time zero was used as a control. The methane production rate was calculated using the equation $a_{\text{p}}/\left(a_{\text{r}}t\right)C$, where a_p and a_r are the activities of the product and added reactant, respectively; t is the incubation period; and C is the in situ concentration of the reactant.

2.4 Direct cell counts

A total of 20 mL of formation water sample was filtered through a 0.2 µm pore size Isopore membrane filter (Millipore), stained with SYBR Green solution (10 µg mL⁻¹) for 10 min, and observed under an epifluorescence microscope (BX51; Olympus, Tokyo, Japan). The detection limit for cell abundance was 4.3×10² cells mL⁻¹.

2.5 Quantitative PCR

DNA was extracted from the membrane filter filtered through formation water samples using a PowerWater kit (MoBio Laboratories, CA, USA) according to the manufacturer's protocol. Quantitative PCR targeting bacterial and archaeal 16S rRNA genes and the mcrA gene (encoding a methyl-coenzyme M reductase alpha subunit) in formation water samples was performed as described previously (Katayama et al., 2016).

2.6 Sequencing of the 16S rRNA and mcrA genes

The 16S rRNA genes (including variable regions 4 and 5) were amplified and sequenced as described previously (Katayama et al., 2016). In brief, the primers Univ515F (5'-GTGYCAGCMGCCGCGGTA) and Univ926R (5'-CCGYCAATTCMTTTRAGTT) were used for PCR amplification. Amplicons were sequenced on an Illumina MiSeq platform (Illumina, Inc., San Diego, CA, USA) using a 250 bp paired-end protocol at the J-Bio21 Center, Nippon Steel and Sumikin Eco-Tech Corporation (Tsukuba, Japan).

The mcrA genes were amplified, cloned, and sequenced from the FW and methanogen culture samples (as described below) as described previously (Katayama et al., 2016) with the following modifications. The primers MLf and MLr (Luton et al., 2002) were used for PCR amplification, and clones were subjected to Sanger sequencing using Eurofin Sanger sequencing services (Tokyo, Japan).

2.7 Sequence analysis

The 16S rRNA gene amplicon reads were analyzed using the Fastx toolkit (version 0.0.14) and QIIME 2 (version 2022.2). The quality-filtered sequences with an average length of 375 bp were clustered as amplicon sequence variants (ASVs) based on the DADA2 algorithm and classified using a Bayesian classifier based on the Silva taxonomy SSU Ref 138.1 dataset (Pruesse et al., 2007) with a confidence threshold of 80 %. In this study, the nomenclature of prokaryotic lineages was based on this SILVA taxonomy.

Sanger sequences of the mcrA genes derived from the FW and culture samples were aligned using MAFFT version 7 (Katoh and Standley, 2013). Nucleotide sequences with >85 % sequence identity were treated as operational taxonomic units (OTUs) according to the cutoff values estimated by Yang et al. (2014). For each OTU, the most abundant sequence was selected as the representative sequence, which was translated to amino acids in silico and then searched against BLAST (http://blast.ncbi.nlm.nih.gov/Blast.cgi, last access: 18 October 2023).

The 16S rRNA gene amplicon data were submitted to the DDBJ Sequence Read Archive database under accession number PRJDB16863. The GenBank/EMBL/DDBJ accession numbers for the mcrA gene sequences are LC783983–LC783991.

2.8 Cultivation of methanogens

The basal media used for the methanogenic cultures consisted of 10 mM NH₄Cl, 1 mM KH₂PO₄, 15 mM MgCl₂⋅6H₂O, 350 mM NaCl, 1 mM CaCl₂⋅2H₂O, 30 mM NaHCO₃, 1 mL L⁻¹ selenium and tungsten solution, 1 mL L⁻¹ trace element solution, 1 mL L⁻¹ vitamin solution, 1 mL L⁻¹ resazurin solution (1 mL L⁻¹), and 0.1 mM titanium(III) citrate (as a reducing agent) (Katayama et al., 2020). Using an anaerobic chamber, 40 mL of the FW sample was mixed with 20 mL of this saline mineral medium and then dispensed into 110 mL serum vials. The vials were sealed with butyl rubber septa and aluminum crimps under an atmosphere of ${\mathrm{N}}_{\mathrm{2}}/{\mathrm{CO}}_{\mathrm{2}}$ (80:20, $v/v$; 0.1 MPa). The culture was supplemented with either ${\mathrm{H}}_{\mathrm{2}}/{\mathrm{CO}}_{\mathrm{2}}$ (80 : 20, $v/v$; 0.1 MPa) plus formate (20 mM) plus yeast extract (0.02 %) (for culturing hydrogenotrophic methanogens), acetate (20 mM) (for culturing acetoclastic methanogens), or methanol (10 mM) plus trimethylamine (10 mM) (for culturing methylotrophic methanogens). The incubation temperatures were set at a value close to the measured water temperature of each sample (as described for the radiotracer measurement methods). Methane production was measured using a GC equipped with a thermal conductivity detector (TCD). After methane production was terminated, 4 mL of the enrichment culture was harvested by centrifugation, and mcrA gene sequencing analysis was conducted as described above.

3 Results

3.1 Geochemistry of the formation water

The geochemical properties of the FW samples are summarized in Tables S1 and S3. The water temperatures ranged from 38 to 81 °C in the FW samples from upper aquifers (Fig. 1c) and from 67 to 96 °C in the samples from lower oil deposits (Table S1). The FW from the upper aquifers was characterized by higher concentrations of bicarbonate and ammonium ions, comparable concentrations of chloride and sodium ions, and almost no sulfate in comparison with present-day seawater (Table S3). The depletion of electron acceptors for microbes, such as nitrate, nitrite, and sulfate, resulted in a low redox potential of FW lower than −180 mV.

The crude oil sample from the lower oil deposits contained n-alkanes with carbon numbers of 9–35, whereas such alkanes were not detected in the FW samples from the upper aquifers (Fig. S2). On the other hand, volatile hydrocarbons such as toluene and xylene were detected in the FW samples from the upper aquifers and were most likely derived from oil. In the 1049 and 1115 m b.g.s. samples, elemental sulfur was also detected. The concentrations of phenols dissolved in the FW samples increased from 1115 m b.g.s. to greater depths (above >73 °C), whereas toluene or xylene did not (Table S3). Similar increasing trends at greater depths were also observed for formate, acetate, propionate, and total organic carbon. A steep increase in acetate concentration at >70 °C has also previously been observed in deep subseafloor sediments (Heuer et al., 2020; Parkes et al., 2007), possibly due to the decomposition of sedimentary organic matter along with elevated geothermal heating (Parkes et al., 2007). Considering that phenols are prominent pyrolysis products of terrestrial kerogens according to previous studies (Larter and Senftle, 1985; van de Meent et al., 1980), it is conceivable that the presence of these compounds is also due to the decomposition of sedimentary organic matter associated with geothermal heating.

3.2 Stable isotopic analysis of sulfur compounds

Figure 2 illustrates the depth distribution of stable sulfur isotope ratios in different sulfur species, namely, sulfide, elemental sulfur, and sulfate in FW and sulfur compounds in oil. Elemental sulfur exhibited the lowest δ³⁴S (−3.3 ‰ to −2.0 ‰), followed by hydrogen sulfide (+1.5 ‰ to +9.4 ‰), sulfur compounds in oil (+8.6 ‰ to +9.6 ‰), and sulfate (+12.1 ‰ to +19.4 ‰). Differences in δ³⁴S values were not clearly observed between samples from upper aquifers and lower oil deposits for each sulfur species.

Figure 2 Depth distribution of stable sulfur isotopic ratios in sulfur species in the formation water (FW, filled symbol) and oil (open symbol) samples. The arrows indicate the possible changes in the ratios due to isotopic fractionation associated with the following reactions: thermal decomposition of oil sulfur to hydrogen sulfide (A), oxidation of upward-migrated hydrogen sulfide (broken arrow) with iron oxides to elemental sulfur (B), and biological dissimilatory reduction of elemental sulfur to hydrogen sulfide (C) (for details, see the Discussion section).

3.3 Hydrogen and oxygen isotope ratios of the formation water samples

To assess the potential mixing of formation water (FW) from upper aquifers and lower oil deposits, isotopic values (δ¹⁸O and δD) of water samples collected from specific layers were determined. In the δ¹⁸O–H₂O vs. δD–H₂O diagram (Fig. S3), FW samples collected from aquifers above 1115 m b.g.s. plot along the meteoric water line in this region (Kato and Kajiwara, 1986; Waseda and Nakai, 1983). In addition, the chlorine ion concentration of the FW sample at the uppermost depth of the 554 m b.g.s. sample was lower than that of the present seawater (Table S3), suggesting the mixing of surface water in this sample. In samples below 1373 m b.g.s., there was a trend of increasing δ¹⁸O–H₂O concentration with increasing depth, possibly due to the isotopic exchange between formation water and sedimentary minerals at elevated temperatures (Clayton et al., 1966). The difference in water isotopic ratios, coupled with the presence of aquiclude silt–mudstone layers between (Takano et al., 2001), indicates no mixing of formation water between the upper aquifers and lower oil deposits.

3.4 Potential for methanogenic activity

The ¹⁴C tracer experiments detected the potential activity of methanogenesis from CO₂, acetate, and methanol (hydrogenotrophic, acetoclastic, and methylotrophic pathways, respectively) in all six FW samples from the dissolved-natural-gas-bearing aquifers (Fig. 3a). The depth-related patterns of the activity rates clearly differed among these three methanogenic pathways. Methylotrophic methanogenesis peaked at the shallowest sample (559 m b.g.s.) and decreased by at least 3 orders of magnitude in deeper samples. The hydrogenotrophic methanogenic rates were greater in the three shallowest samples (<845 m b.g.s.) than in the deeper samples (>1049 m b.g.s.), and the maximum rates were observed in 845 m b.g.s. The activity rates of acetoclastic methanogenesis were at least 2 orders of magnitude greater in deeper samples (1115 and 1373 m b.g.s.).

Figure 3 Depth (temperature)-related changes in potential methanogenic activity (a), microbial population (b), and prokaryotic (c) and methanogen (d) community compositions based on the 16S rRNA (c) and mcrA (d) gene sequences in formation water samples from the upper aquifers. The growth temperatures (e) are expressed by the detection of mcrA gene sequences from both original formation water samples and methanogen cultures at different temperatures (for details, see Table S4). Related species of mcrA gene sequences and amino acids sequence identities (in parentheses) are shown (e).

3.5 Enumeration of total microbial cells and of the 16S rRNA and mcrA genes

The number of microbial cells in the FW samples ranged from 1.1×10³ to 4.8×10⁴ cells mL⁻¹ (Fig. 3b). The highest and lowest numbers were measured in the shallowest 559 and 845 m b.g.s. samples, respectively. No microbial cells were observed in the oil-associated FW sample (1733 m b.g.s.), in which water temperature was measured to be 96 °C (Table S3). Therefore, molecular gene sequencing analysis was not conducted on oil-associated FW samples from lower oil deposits.

The populations of bacteria, archaea, and methanogens were measured via quantitative real-time PCR (Fig. 3b). The bacterial 16S rRNA gene copy numbers ranged from 10²–10⁴ mL⁻¹, whereas those of archaea ranged from 10²–10³ mL⁻¹. Bacterial and archaeal gene copy numbers were highest in the deepest sample (1373 m b.g.s.). On the other hand, in this sample, the mcrA gene copy number was below the detection limits. The mcrA gene copy numbers tended to decrease with increasing depth.

3.6 Microbial community compositions

16S rRNA gene amplicon analysis was performed to examine the compositions of the prokaryotic and methanogenic communities. After quality filtering, the reads yielded 88 255–108 428 reads per sample. The major taxonomic groups (>10 % of the total reads in at least one sample) were the phyla or classes Bacillota (formerly Firmicutes), Gammaproteobacteria, Deferribacterota, Desulfobacterota, and Methanobacteriota (formerly Euryarchaeota) (Fig. 3c). The members of Bacillota were abundant in two deeper samples (1115 and 1373 m b.g.s.), accounting for >85 % of the total sequences. In the two shallower samples (559 and 754 m b.g.s.), Gammaproteobacteria and Desulfobacterota sequences were abundant, whereas Gammaproteobacteria and Deferribacterota sequences were abundant in the two intermediate-depth samples (845 and 1049 m b.g.s.). In the samples deeper than 1000 m b.g.s., only one ASV accounted for ≥50 % of the total sequences in the sample. These genes were phylogenetically related to Deferribacter desulfuricans (100 % sequence similarity, in the 1049 m b.g.s. sample), Thermacetogenium phaeum (93 % sequence similarity, in the 1115 m b.g.s. sample), and Thermanaeromonas toyohensis (97 % sequence similarity, in the 1373 m b.g.s. sample). These species can commonly utilize sulfur compounds, such as thiosulfate and elemental sulfur, as electron acceptors and acetate as an electron donor (Takai et al., 2003; Hattori et al., 2000; Mori et al., 2002) Additionally, these genera are known to utilize a variety of electron acceptors: Deferribacter can use nitrate, iron(III), or manganese (Slobodkina et al., 2009); Thermacetogenium can use sulfate (Hattori et al., 2000); and Thermanaeromonas can use nitrate, nitrite, sulfate, or fumarate (Gam et al., 2016).

The archaeal sequences accounted for 0.04 %–18.99 % of the total sequence. The majority (>3 % of the total archaea) were assigned to Halobacterota, Methanobacteriota, and Thermoplasmatota, whereas the minor sequences were Ca. Asgardarchaeota, Ca. Hadarchaeota, and Nanobdellota (formerly Ca. Nanoarchaeota). The sequences of hydrogenotrophic methanogens, such as Methanobacterium and Methanocalculus, were abundant in shallower-depth samples (<1050 m b.g.s.), whereas Archaeoglobaceae were abundant in deeper samples (1115 and 1373 m b.g.s.). Members of the Archaeoglobaceae family are known to reduce sulfur compounds (Liu et al., 2012). On the other hand, recent research has reported Mcr-dependent methanogenesis in Archaeoglobaceae from a hot spring (Buessecker et al., 2023). Based on this literature, acetate conversion to methane is possible through methanogenic pathways deduced from the reconstructed Archaeoglobaceae genomes (Buessecker et al., 2023). Other potential candidates for methanogens or hydrocarbon-oxidizing archaea known to possess the Mcr or Acr (alkyl-S-CoM reductase) gene in their genomes, such as Ca. Methanomethyliales, Ca. Nezhaarchaeales, Ca. Korarchaeia, Ca. Methanodesulfokores, Ca. Bathyarchaeia, and Ca. Helarchaeales (Mei et al., 2023), were not detected in any of the samples.

3.7 Methanogen diversity in FW and culture samples based on the mcrA genes

A total of 29–46 clones per sample of the mcrA gene were grouped into eight OTUs (Fig. 3d). Among these, five OTUs were also detected in the methanogenic cultures of the FW samples (Table S4). Similar to the results of the 16S rRNA gene sequencing analysis, methanogen diversity changed with depth. Methanobacterium and Methanocalculus were dominant in the samples from depths above 1050 m b.g.s., whereas only thermophilic Methanothermobacter was recovered in the deeper samples at 1115 and 1373 m b.g.s., resulting in an overall predominance of hydrogenotrophic methanogens throughout the depths. In 559 and 754 m b.g.s., the methylotrophic methanogen Methanolobus was also detected as a major proportion. Acetoclastic methanogens, such as Methanosaeta and Methanosarcina, were not detected in any of the samples. Additionally, unlike the results from the 16S rRNA gene sequencing analysis, Archaeoglobaceae were not detected in the samples from 1115 and 1373 m b.g.s. The absence of mcrA gene sequences for Archaeoglobaceae, which might participate in acetoclastic methanogenesis in these samples (as described above), highlights the bias and limitation of PCR-based mcrA gene sequencing analysis. This issue has been documented by Vítězová et al. (2021), who indicated that mcrA gene analysis fails to detect specific methanogens, with their identification depending on the primers used. It could also explain the observed discrepancy between the highest ¹⁴C acetoclastic activity and the lowest mcrA gene copy numbers at these high-temperature depths.

These depth-related distributions of methanogens were roughly correlated with their growth temperature, as observed in the enrichment cultures (Fig. 3d and e; Table S4). Methane production from H₂ and formate was observed in the enrichment cultures at depths from 37 to 80 °C, except at 70 °C, whereas no methane production was observed in acetate-amended cultures from any FW samples (data not shown). In cultures supplemented with methylated compounds, methane was detected only in 559 and 754 m b.g.s. samples, i.e., at 37 and 50 °C, respectively.

4 Discussion

This study examined how microbial diversity and processes are shaped by geothermal heat and associated geochemical processes in deep subsurface environments, encompassing a range of temperatures from 38 to 81 °C. The geochemical and hydrological analyses suggested the upward migration of volatile hydrocarbons derived from lower oil deposits to upper aquifers (Fig. S2 and Table S3), whereas no evidence indicated the mixing of water or less volatile n-alkanes between the upper aquifers and lower oil deposits (Figs. S2 and S3). Radiotracer measurements provided clear insight into the variations in potential methanogenesis across different temperature regimes (Fig. 3a). Molecular and cultivation experiments indicate that these variations align with the growth temperature ranges of major indigenous methanogens, except for high-temperature depths (≥70 °C) (Fig. 3d). Our previous study also suggested that the upper temperature limit for growth in indigenous methanogens may reflect the potential for methanogenesis in deep subseafloor sediments (Katayama et al., 2022). The correlation between the high potential for acetoclastic methanogenesis and increased acetate concentration suggests that acetate generation, possibly driven by geothermal heating (Parkes et al., 2007), stimulates acetoclastic methanogenesis in high-temperature depths (≥70 °C). Although syntrophic acetate oxidation (SAO) coupled with hydrogenotrophic methanogenesis is generally observed in high-temperature environments (Conrad, 2023), the potential for this SAO methanogenesis is unlikely at our study site. This inference is drawn from the observation that the levels of ¹⁴C activity from acetate were as much as 2 orders of magnitude greater than those from carbonate in high-temperature depths (≥70 °C). If SAO was the primary activity, the ¹⁴C activity from acetoclastic methanogenesis would be much lower than that from hydrogenotrophic methanogenesis. This is because ¹⁴C-CO₂ produced from ¹⁴C-acetate via SAO is diluted into the bicarbonate pool in the incubated FW sample, leading to the monitoring of ¹⁴C-CH₄ production from CO₂ with an extremely low level of ¹⁴C. We speculate that the predominant Archaeoglobaceae may play a role in converting ¹⁴C-acetate to ¹⁴C-methane at greater depths (1115 and 1373 m b.g.s.).

In addition to driving the acetoclastic methanogenesis observed in the radiotracer experiments, elevated concentrations of acetate at greater depths (1115 and 1373 m b.g.s.) may primarily serve for the dissimilatory reduction of sulfur compounds other than sulfate, such as elemental sulfur and thiosulfate. This is inferred from the high relative abundance of sulfur-reducing, rather than sulfate-reducing, microorganisms at those depths. This raises the question of how sulfur compounds are generated to sustain these sulfur-reducing microbes, possibly reflecting the relatively high microbial population at greater depths. Based on the sulfur isotopic analysis and previous literature, we suggest that elemental sulfur and/or thiosulfate are abiotically generated from sulfur compounds within crude oils according to the following reactions. Hydrogen sulfide can be generated through the thermal decomposition of sulfur compounds in crude oil (Zhu et al., 2017). The occurrence of this reaction is supported by the lower δ³⁴S values of hydrogen sulfide in the FW sample from the lower oil deposits compared with those of sulfur compounds in oil (arrow A in Fig. 2). This generated hydrogen sulfide, together with other oil-derived volatile hydrocarbons, such as toluene and xylene, may migrate upward to upper aquifers and undergo chemical oxidation facilitated by the presence of Fe(III)-bearing minerals in the sediments. This process leads to the formation of elemental sulfur or thiosulfate. The occurrence of this reaction is supported by the lower δ³⁴S values observed for elemental sulfur than for hydrogen sulfide in the FW sample from lower oil deposits (arrow B in Fig. 2). Previous studies have indicated the occurrence of chemical sulfide oxidation with metal oxides, such as FeOOH and magnetite (${\mathrm{Fe}}^{\mathrm{2}+}{\mathrm{Fe}}_{\mathrm{2}}^{\mathrm{3}+}{\mathrm{O}}_{\mathrm{4}}$), in marine sediments (Holmkvist et al., 2011; Bottrell et al., 2008). Given that FeOOH is more readily reduced than magnetite during sediment burial, resulting in the preservation of magnetite (Canfield et al., 1992), magnetite likely serves as the primary agent responsible for sulfide oxidation at this study site. Despite the predominance of sulfur-reducing microorganisms in the upper aquifers as opposed to the absence of microbial cells in the lower oil deposits, the differences in δ³⁴S values of hydrogen sulfide and elemental sulfur are not substantial between the upper and lower locations. This suggests that biological dissimilatory sulfur reduction (arrow C in Fig. 2) is much less important than chemical sulfide oxidation (arrow B in Fig. 2) with respect to its influence on δ³⁴S, which should be linked to differences in the sulfur flux of the biological and chemical processes in the upper aquifers.

A previous study indicated that the degree of generation of hydrogen sulfide from crude oil depends on various factors, including the concentration of sulfur compounds in the oil, oil biodegradation coupled with microbial sulfate reduction, and oil densification during burial (Zhu et al., 2017). At the study site, the concentration of sulfur compounds in oil source rocks, as represented by the sulfur-to-carbon atomic ratio ($\mathrm{S}/\mathrm{C}$), was relatively low, ranging from 0.011–0.049, with an average of 0.025 (Suzuki et al., 1995). For comparison, high-sulfur oils exhibit an $\mathrm{S}/\mathrm{C}$ ratio greater than 0.06 (Orr, 1986). This implies a low level of upward migration of hydrogen sulfide from oil deposits to upper aquifers. Nevertheless, our findings indicate that the shifts in microbial population, diversity, and function may be due to the chemical reactions induced by geothermal heating, which underscores the potential for more dynamic chemical and microbial interactions that drive sulfur and carbon cycling in other deeply buried sedimentary environments.