Influence of minor hydrocarbon seepage on sulfur cycling in marine subsurface sediments

All hydrocarbon (HC) reservoirs leak to some extent. When small quantities of HCs escape offshore reservoirs and migrate through overlying organic-poor marine sediments towards the surface, these HCs are often completely metabolized by microorganisms before reaching the sediment–water interface. However, these low and often unnoticed upward HC fluxes still influence the geochemistry of surrounding sediments and potentially stimulate metabolic activity by microbial populations in shallow subseafloor environments.

In this study, we examined how localized HC seepage affects microbial sulfate reduction in organic-poor sediments from the SW Barents Sea, focusing on three sampling areas overlying known HC deposits and two pristine seabed reference areas. The analysis of 50 gravity cores revealed significant variability in the predicted sulfate depletion depth, ranging from 3 to 12 m below the seafloor. We observed nearly linear pore water sulfate and alkalinity profiles, alongside low rates of sulfate reduction (pmol cm³ d⁻¹). Metagenomic and metatranscriptomic data indicated metabolic potential and activity towards sulfate reduction and anaerobic oxidation of methane (AOM). Expression of functional marker genes (aprAB, dsrAB, mcrA) revealed metabolic activities by sulfate-reducing Desulfobacterota and methanotrophic ANME-1 archaea sustained by HC traces in the sediment. Furthermore, in concomitance to AOM processes, we found that archaea amongst the classes Lokiarchaeia and Thorarchaeia were expressing genes involved in sulfur reduction (hydB, hydG).

Overall, our study demonstrates that the gradient in pore water geochemistry, the rates of sulfate reduction processes, and the genetic features of microbial populations actively involved in sulfate-driven AOM processes are all affected by low and often unnoticed upward HC fluxes. This slight HC seepage resulted in sulfate depletion at shallower depths compared to what would be expected at a comparable non-affected site and produced concomitant biogeochemical signatures in the shallow subsurface that enable the inference of deeply buried reservoirs.

1 Introduction

All hydrocarbon (HC) reservoirs leak to some degree (Hunt, 1995; Yergin, 2009; Heggland, 1998), with over 80 % of all seeps occurring directly above the reservoirs (Ciotoli et al. 2020). While large leakages result in conspicuous manifestations at the seafloor (e.g. Cramm et al., 2021), minor seeps, characterized by low, primarily diffusive HC fluxes, usually remain invisible at the sediment–water interface (SWI), as there are no visible signs of seepage, e.g. carbonate crusts or chemosynthetic communities such as sulfur bacteria mats, clams, and tubeworms. In cases where HC fluxes are so low that hydrocarbons are fully mineralized and do not reach the SWI, sediment geochemistry and microbiology at minor seepage sites will be altered (Rasheed et al., 2013; Abrams, 2020; Riedinger et al., 2014; Joye et al., 2004). Even if HC fluxes were higher and if hydrocarbons do indeed reach the SWI, they are often overlooked (Abrams, 2020; Joye, 2020).

HC seepage triggers both geochemical and microbiological processes, changing physicochemical properties, such as sediment porosity, mineral and pore water composition, and microbial community composition and activity (Abrams, 2020; Hvoslef et al., 1996; Joye, 2020). Likewise, the presence of free gas in methane-containing pore fluids alters sediment porosity and compressibility (Jang et al., 2018), which can be detected through seismic profiling due to the resulting impedance contrast (Rovere et al., 2020). These changes may also affect the sediment's geochemical properties, such as clay texture, organic matter (OM) reactivity, and metal content (Chen et al., 2023). HC influx and subsequent HC consumption will also alter pore water geochemistry and pH, promoting ion exchange reactions, dissolution and precipitation (Smrzka et al., 2019; Turchyn et al., 2021; Torres et al., 1996), and structural conversion of minerals (Jiang, 2012). As a consequence, specific trace metals and metalloids (e.g. Cr, Cu, As, Se, Sb) that can be metabolized by specific microorganisms (Raab and Feldmann, 2003) tend to become enriched in the sediment pore water at seepage sites (Rasheed et al., 2013; Guseva et al., 2021). Finally, most offshore sediments are composed of pelagic deposits, i.e. clays, silts, or biogenic oozes, containing only low amounts of bioavailable OM that become increasingly recalcitrant with depth of burial (Middelburg, 2018). Thus, the influx of HCs in such sedimentary systems represents an additional source of electron donors, providing chemical energy for microorganisms and potentially fuelling various metabolic processes. The metabolic potential to break down and assimilate HCs has been identified in a large number of microorganisms, e.g. reviewed in Joye (2020), whose activities may induce concomitant effects in the pore water, organic, and mineral fraction of the sediment (Kim et al., 2004).

The quantitatively most important anaerobic OM degradation process in the upper layers of marine sediments, where oxygen is depleted but sulfate is still present in relatively high concentrations, is sulfate reduction (Jørgensen et al., 2024; Bowles et al., 2014; Kasten and Jørgensen, 2000; Jørgensen, 1982). Due to its abundance in seawater, sulfate diffuses into the sediment pore water, where it is used as an electron acceptor by sulfate-reducing microorganisms composed of sulfate-reducing bacteria and ANME (Widdel et al., 2010). Sulfate reduction (SR) can be divided into two main pathways, namely organoclastic and methanotrophic SR. The respective overall reactions are as follows:

$$\begin{array}{}\text{(R1)}& {\displaystyle }\begin{array}{rl}& \text{organoclastic SR:}\\ & \mathrm{2}{\mathrm{CH}}_{\mathrm{2}}\mathrm{O}+{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}\to \mathrm{2}{\mathrm{HCO}}_{\mathrm{3}}^{-}+{\mathrm{H}}_{\mathrm{2}}\mathrm{S},\end{array}\text{(R2)}& {\displaystyle }\begin{array}{rl}& \text{methanotrophic SR:}\\ & {\mathrm{CH}}_{\mathrm{4}}+{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}\to {\mathrm{HS}}^{-}+{\mathrm{HCO}}_{\mathrm{3}}^{-}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}.\end{array}\end{array}$$

Due to the abundance of easily biodegradable OM and sulfate close to the sediment–water interface, organoclastic SR normally prevails in anaerobic OM degradation within the uppermost part of the sediment column, leading to a decrease in sulfate concentration. Methanotrophic SR becomes predominant only when biogenic or thermogenic methane diffusing from deeper sediments reaches the zone where sulfate is still available but approaches depletion (Martens and Berner, 1974). The narrow zone where both methane and sulfate are available is referred to as the sulfate–methane transition zone (SMTZ). Within the SMTZ, mainly methane, potentially along with other OM compounds (Beulig et al., 2019; Jørgensen et al., 2019b), is oxidized by microorganisms, while sulfate is reduced (Eq. 2).

Thus, the respective methane and sulfate fluxes exert control on the depth of the SMTZ (Borowski et al., 1996, 1999; Jørgensen and Kasten, 2006), which, under steady-state conditions, is usually a reflection of methanogenesis rates (Henrichs and Reeburgh, 1987) or fluxes of thermogenic methane (Hu et al., 2017). However, the SMTZ is highly dynamic and responds relatively rapidly (i.e. within decades) to shifts in pore water geochemistry (Sultan et al., 2016; Hong et al., 2016). Thus, microbially active subseafloor sediments exhibit significant variability in the depth of the SMTZ across different oceans, from 1 to 10 m on the continental shelf (Egger et al., 2018; Treude et al., 2020; Treude et al., 2014; Riedinger et al., 2017; Riedinger et al., 2014; März et al., 2008; Niewöhner et al., 1998). On the continental slope, the SMTZ can be down to 100 m below seafloor (m b.s.f.) in extreme cases (Nunoura et al., 2009) but is usually not deeper than 25 m b.s.f.

The subseafloor biosphere hosts phylogenetically diverse and metabolically active microbial communities (Parkes et al., 2014; D'Hondt et al., 2004). Their metabolic activity depends mainly on biogeochemical cycling of nitrogen, phosphate, and sulfur. At HC seepage sites, microbial communities are well defined in terms of taxonomy and highly specialized in terms of metabolic functions, and they are referred to as “seep biomes” (Ruff et al., 2015). It has been known for many decades that methanogenic hydrocarbon degradation is carried out jointly by bacteria and archaea in underground hydrocarbon reservoirs. New studies suggest that Candidatus methanoliparum may be able to independently degrade long-chain alkanes and produce methane (Wegener et al., 2022; Zhou et al., 2021; Zhang et al., 2024). Seep biomes show spatial heterogeneity but are distinct from non-seep environments (Pop Ristova et al., 2015). The reported key microbial guilds include methane and oil oxidizers, sulfate reducers, and nitrogen fixers. Taxa found near HC seeps involve members of the bacterial phyla Desulfobacterota and Chloroflexota that can, respectively, use simple and halogenated HCs as electron acceptors (Kleindienst et al., 2014; Zhang et al., 2023), whereas archaeal groups, like the ANME or Asgardarchaeota, consume methane or make use of HCs or sulfur compounds (Firrincieli et al., 2021; MacLeod et al., 2019).

Our study aims to determine in what ways minor HC seepage modifies the diversity and activity of anaerobic methanotrophs and sulfate-reducing bacteria in surface sediments of the SW Barents Sea. We hypothesize that minor HC seepage influences pore water geochemical properties, even in near-surface sediment, and that resident microbial communities respond quickly in terms of composition and activity to the geochemical alteration of their habitats.

For this, we collected a total of 50 gravity cores, of which 40 cores are from areas affected by light HC seepage from underlying HC reservoirs (i.e. HC-affected sites), but without any visible surface manifestations like carbonate crusts or chemosynthetic communities such as sulfur bacteria mats, clams, and tubeworms, and 10 cores are from seabed areas without seepage (i.e. reference sites). We combine detailed analyses of pore water sulfate and alkalinity with potential and modelled sulfate reduction rates and with dissolved methane (CH₄) and carbon dioxide (CO₂) concentrations, and we supplement them with metagenomic (i.e. total DNA) and metatranscriptomic (i.e. total RNA) data on taxonomic and functional marker genes in order to trace metabolic activity by microbial sulfate-reducing consortia in their geochemical context.

2 Material and methods

2.1 Geological context of the Barents Sea, sampling sites

The Barents Sea is an epicontinental shelf sea surrounded by the marginal Norwegian Sea to the southwest, the Arctic Ocean to the north, and the Russian archipelago Novaya Zemlya to the east (Fig. 1a). Since its formation during the Caledonian orogeny (i.e. 490–390 Ma), the Barents Sea Basin experienced several phases of uplift and subsidence, with subsequent tilting and erosion. These geologic events resulted in the characteristic structure of the Barents Sea in the form of a succession of structural highs and basins (Gabrielsen et al., 1990; Larssen et al., 2002).

Figure 1 Sampling locations. (a) Map of the Norwegian coast. The square marks the sampling area. (b) A total of 50 cores were retrieved from three HC-affected sites (Zone 1, Zone 10, Zone 12) and two reference sites (Ref 5, Ref 6) (bathymetric map: EMODnet project (European Commision), 2022). (c) Bathymetric map of the three HC-affected sites. (d) Each site was sampled in a grid of cores, Zone 1 is enlarged here as an example. The missing samples in the middle are due to the increased occurrence of dropstones in this area, preventing the deployment of a gravity corer.

HC exploration in the Barents Sea has a history stretching back several decades (Doré, 1995; Johansen et al., 1993). The sediments have been characterized in numerous studies (e.g. Sættem et al., 1991; Elverhøi and Solheim, 1983) and are mainly composed of organic-poor silty clays (Knies and Martinez, 2009; Nickel et al., 2013) originating from eroded sedimentary and igneous rocks.

During an expedition to the southern Loppa High region of the Barents Sea (Fig. 1a) from 29 October to 8 November 2021, we collected 50 gravity cores (up to 3 m in length) at about 350 m water depth in three areas with known underlying HC reservoirs (i.e. Zone 1, Zone 10, Zone 12) and in two reference areas (i.e. Ref 5, Ref 6) without underlying HC reservoirs (Fig. 1b), labelled hereafter as HC-affected sites and reference sites, respectively. The cores were subsampled for sediment and pore water on board the research vessel.

2.2 Pore water and sediment sampling

Pore water samples were extracted from all cores at a resolution of 10 cm, using rhizons (Rhizosphere Research Products) (Seeberg-Elverfeldt et al., 2005), collecting a minimum volume of ca. 5 mL over 24 h. The pore water was filtered through 0.2 µm pore size cellulose ester (MCE) membrane syringe filters (Merck MF-Millipore) and aliquoted as follows: 1.5 mL transferred into a plastic screw cap vial without further treatment for anion and cation measurements, 2 mL into a glass vial amended with 50 µL saturated HgCl₂ solution to inhibit microbial activity and sealed without headspace for alkalinity measurements (Edenborn et al., 1985), and 1.5 mL into a plastic screw cap vial mixed with 200 µL ZnCl₂ (20 % weight vol⁻¹) to precipitate dissolved hydrogen sulfide as ZnS for subsequent sulfide measurements. All samples were stored at +4 °C until analysis in the home lab.

Immediately after core recovery, the lowermost 20 cm of each sediment core was cut off and subsampled for DNA and RNA extraction. The sediment was pushed out of the liner, and the outer 1–2 cm of sediment rims was removed with a sterile spatula. Aliquots of the remaining sediment were transferred into gas-tight aluminium foil bags, flushed with nitrogen gas, heat-sealed, and stored at −80 °C for later extraction of various biomolecules.

2.3 Pore water sulfate, sulfide, and alkalinity

Anion concentrations were determined using an ion chromatography (IC) system equipped with a SykroGel A × 300 AB-A01 column (all Sykam). The eluent contained 7.3 mg L⁻¹ NaSCN and 636 mg L⁻¹ NaCO₃. The pump rate was set to 1 mL min⁻¹, and the injected sample volume was 50 µL. All samples were measured in triplicate, and the results were averaged. Because of their high salt concentrations, all samples had to be diluted 1:50 with Milli-Q water before injection. The detection limit of the IC system is 0.5 mg L⁻¹. The average standard deviation of replicate measurements was always better than 3 %.

Sulfide concentrations were measured according to the protocol of Cline (1969). Aliquots with absorption values greater than 1 were diluted and remeasured. Sulfide concentrations were calculated by comparing adsorption values measured for pore water samples to those of an Na₂S standard. The detection limit of the method is 0.1 µM. All samples were measured in triplicate, and the results were averaged, with standard deviation ≤ 3.5 %.

Pore water alkalinity was determined using the Visocolor HE alkalinity AL 7 kit (Macherey-Nagel GmbH). As we only had a total of ca. 2 mL of pore water per sample, we scaled down the volume of sample to 500 µL in order to perform all measurements in triplicate. For this, 10 µL of indicator was added to 500 µL of pore water and titrated with the Visocolor solution until we observed a change in colour from blue to orange. Alkalinity of the pore water samples was calculated based on the titration solution volume. Replicates differed by less than 3 % on average.

2.4 Sulfate reduction rates and modelled sulfate reduction rates

For sulfate reduction measurements, each gravity core was subsampled in triplicate every 40–50 cm and at the bottom end. Structurally intact sediment plugs (volume ca. 3.5 cm; Environmental Metagenomics, Research Center One Health Ruhr of the University Alliance Ruhr, Faculty of Chemistry, University of Duisburg-Essen, Essen, Germany) were recovered by inserting 5 mL glass barrels fitted with a syringe plunger into the sediment. Following the whole core injection method (Jørgensen, 1978), 15 µL of radioactive ³⁵SO₄–sulfate tracer (200 kBq) was injected into the glass barrels containing the samples retrieved on the same day (Fossing et al., 2000). The glass barrels were closed with butyl rubber stoppers, and the samples were incubated at in situ temperature (ca. 4 °C) for 24 h. Each sample was then transferred into a previously weighed 50 mL centrifuge tube containing 10 mL of 20 % zinc acetate solution to stop all microbial activity and to precipitate all volatile hydrogen sulfide as zinc sulfide. Each vial was then thoroughly shaken to break down all sediment aggregates and stored frozen until analysis. Blank samples were prepared by adding the radiolabelled sulfate tracer just before transferring the sample into the zinc acetate solution. The tracer blank was prepared in the same way without sediment.

Sulfate reduction rates (SRRs) were quantified in the home lab, using single-step cold chromium distillation (Kallmeyer et al., 2004). The sample vials were thawed and centrifuged, and the supernatant was carefully removed. A small aliquot of supernatant was kept for quantification of total radioactivity, and the rest was discarded. The sediment pellet was then quantitatively transferred into a distillation flask. Due to the clayey consistency of the sediments, we had to re-freeze the vials and transfer the frozen sediment pellet from the centrifuge tube to the distillation flask. SRRs were calculated according to the following formula:

$$\begin{array}{}\text{(1)}& \mathrm{\text{sulfate reduction rate:}}\mathrm{SRR}={\displaystyle \frac{{\mathrm{SO}}_{\mathrm{4}}\cdot \mathit{\phi }\cdot a_{\mathrm{tris}}\cdot \mathrm{1.06}\times {\mathrm{10}}^{\mathrm{6}}}{a_{\mathrm{tot}}\cdot t}},\end{array}$$

where SRR is the sulfate reduction rate (pmol cm⁻³ d⁻¹), SO₄ is the sulfate concentration in the pore water (mmol L⁻¹), φ is the porosity of the sediment set to 0.7, a_tris is the radioactivity of total reduced inorganic sulfur (cpm), a_tot is the total radioactivity used in the spiking (cpm), t is the incubation time (d), 1.06 is the correction factor for isotopic fractionation, and 10⁶ is the unit conversion factor (from µmol L⁻¹ to pmol L⁻¹).

Assuming steady-state conditions, the net rate of sulfate production, or consumption, can be modelled based on the measured sulfate concentration profiles by the diffusion–reaction modelling software PROFILE (Berg et al., 1998). This software considers three different kinds of vertical transport (i.e. diffusion, bioturbation, and irrigation) together with the flux across the SWI. The final output corresponds to the minimum number of intervals, as defined by their respective sulfate production or consumption rates, that is necessary to reproduce the measured pore water profiles. For the modelling calculation, we used the coefficient of sulfate diffusion in water at 4 °C $D_{\mathrm{Sulfate}}=\mathrm{0.56}\times {\mathrm{10}}^{-\mathrm{5}}$ cm s⁻¹ (Iversen and Jørgensen, 1993), a porosity (φ) of 0.7, and a seawater sulfate concentration of 28 mmol cm⁻³. Because the sediment appeared to be free of macrobenthos, the biodiffusivity and irrigation coefficients were set to D_B=0 cm s⁻¹ and α=0, respectively.

2.5 Methane and carbon dioxide concentrations

For dissolved CH₄ and CO₂ concentrations, we subsampled the lowermost 20 cm of sediment from each gravity core directly after recovery. Using 5 mL cut-off syringes, we extracted ca. 4 cm of sediment per core. The first and last 0.5 cm of sediment in the syringe was discarded, while the central 3 cm was transferred into a 10 mL glass crimp vial containing a saturated NaCl solution. The vials were immediately sealed with thick butyl rubber stoppers, crimped, and stored upside down without headspace until analysis in the home lab.

Prior to measurement, we introduced 3 mL of ultrapure helium gas as headspace while withdrawing the same amount of NaCl solution from the vial. To equilibrate dissolved gases with the headspace, the content of the vials was mixed at 220 rpm on an orbital shaker for 18 h and further vortexed to break up the remaining small clayey aggregates. We extracted ca. 350 µL of headspace from the vial using a gas-tight syringe. Prior to measurement, ca. 100 µL of the sampled gas was flushed through the injection needle, and 250 µL of gas sample was then injected into a 7890A Gas Chromatograph system equipped with a flame ionization detector (FID), a thermal conductivity detector (TCD), and an HP PLOT Q column (all Agilent). Oven temperature was set to 50 °C, flow rate was set to 17.2 mL min⁻¹, and pressure was set to 13 psi. The detectors both worked at 200 °C with flow rates of 40 mL min⁻¹ (FID) and 15 mL min⁻¹ (TCD). The system was calibrated using 250 µL of analytic pure standards, injecting CO₂ concentrations of 310 and 5270 ppm and CH₄ concentrations of 10 and 5170 ppm. The initial CH₄ and CO₂ concentrations were converted from parts per million (ppm) to molar concentrations by applying the ideal gas law.

2.6 DNA/RNA extractions, sequencing libraries, and functional marker genes

DNA was extracted for a total of 12 samples from eight different cores corresponding to 4, 2, and 2 samples from HC-affected sites (i.e. Zone 1, Zone 10, and Zone 12), respectively, along with 3 and 1 samples from Ref 5 and Ref 6, respectively. Samples were subjected to DNA extraction using the DNEasy PowerMax Soil Kit (QIAGEN, Germany) following the manufacturer's instructions with few adaptations. For 2.5 g of starting sediment material, the final elution was performed in two steps using 2.5 mL of solution C6 twice and was successively incubated and centrifuged. Eluted DNA was then concentrated by adding 510 µL of a 3 M sodium acetate (pH 5.2)/glycogen (0.4 µg µL⁻¹) solution and incubated for 1/,h at −20 °C and then centrifuged to pellet DNA. After ethanol cleaning and air-drying, DNA was resuspended in 100 µL of 10 mM Tris buffer. Finally, DNA was cleaned and concentrated using the Zymo DNA Clean & Concentrator (Zymo, UK) following the manufacturer's instructions and was eluted in a final elution volume of 35 µL.

DNA sequencing was achieved by using 32.5 µL of purified and concentrated DNA for metagenomic library preparation with the QIAseq FX DNA Library UDI A/B kit (QIAGEN, Germany), irrespective of the individual sample concentrations. The genomic DNA was enzymatically fragmented for 10 min, and DNA libraries were prepared with seven PCR cycles for the PCR-less procedure and without a dedicated PCR step for the PCR-free procedure. The average library sizes were ±400 base pairs. Prepared libraries were quantified using a Qubit fluorometer (Thermo Fisher, USA) and quality-checked on a Bioanalyzer (Agilent, USA). Sequencing was performed on an Illumina NextSeq 2000 platform using 2×151 bp read length at the Luxembourg Centre for Systems Biomedicine (LCSB) Sequencing Platform, aiming at an average of 10 G base pairs per sequencing library.

RNA was extracted from 2×5 g of sediment per sample, using the RNeasy PowerSoil Total RNA Kit (QIAGEN) according to the manufacturer's instructions, yielding 30 µL of RNA-containing extract, and quantified on a Qubit 4 fluorometer, using the RNA high-sensitivity kit (Thermo Fisher Scientific). RNA extracts were directly subjected to library preparation in technical duplicates, using the Revelo™ RNA-Seq High Sensitivity kit (Tecan Life Sciences) with 8 µL of RNA extracts as a template. This trio kit includes first- and second-strand cDNA synthesis, adaptor ligation, and library amplification through 10 PCR cycles. Prepared libraries were quantified using a Qubit 4 fluorometer (Thermo Fisher Scientific) and quality-checked on a Bioanalyzer (Agilent) to determine the molarity and overall size distribution of RNA molecule fragments. The average library fragment size was ±300 bps, and the yields were between 4 and 12 ng µL⁻¹. Metatranscriptomic libraries displaying high and low molarities were pooled into two different batches, both sequenced by the commercial company EMBL (Heidelberg, Germany) on a NextSeq 2000 with a P3 kit and an Illumina micro kit (2×150 bps), respectively.

De novo assembly of metagenomic reads was performed (Table S1 in the Supplement) using a customized Snakemake workflow (Mölder et al., 2021), including quality control, binning into contigs, and gene prediction and annotation, as described in Bornemann et al. (2023). Paired-end reads were interleaved, Illumina adapters and controls were removed, and reads were quality-based trimmed and de-interleaved into forward and reverse read files. Bacterial and archaeal scaffolds were concatenated and mapped back to raw reads to retrieve coverage, % GC, and contig lengths. Genes were then predicted by extracting open reading frames (ORFs) from metagenome-assembled scaffolds using Prodigal (Hyatt et al., 2010). Metatranscriptomic reads were quality-filtered and trimmed using Trimmomatic (Bolger et al., 2014). Functional profiles were obtained by mapping the filtered and trimmed reads to their corresponding metagenomic databases using BWA-MEM (Li and Durbin, 2009). ORFs were retrieved using Bakta (Schwengers et al., 2021) and combined with the mapping using the featureCounts tool (Liao et al., 2014).

Taxonomic identifications of all ORFs were further integrated with functional annotations by performing BLASTp searches against an aggregated database of predicted proteins (Orsi et al., 2018), using the DIAMOND protein aligner v. 0.9.24 (Buchfink et al., 2015). Cut-off values for assigning the best match to specific taxa present in the aggregated database were performed at a minimum bit score of 50, minimum amino acid similarity of 60, and an alignment length of 50 residues. We used this approach to draw conclusions about metabolic traits derived specifically from phylum to class taxonomic levels only (Vuillemin et al., 2020a, b).

Here, we focus on ORFs encoding genes involved in sulfur cycling and methane-related processes only. To trace metabolic potential in anaerobic respiration of sulfur compounds (i.e. sulfate, sulfite, thiosulfate, sulfur, polysulfide), we looked for the presence of ORF-encoding genes with similarities to adenylylsulfate reductase (apr), dissimilatory sulfite reductase (dsr), anaerobic sulfite reductase (asr), thiosulfate/polysulfide reductase (phs/psr), and sulfhydrogenase (hyd) (Vuillemin et al., 2022). For methane production and anaerobic consumption, we targeted ORFs encoding proteins of the methyl coenzyme M reductase (mcr). To identify and characterize metabolically active taxa, we targeted the small subunit ribosomal ribonucleic acid (SSU rRNA) as a taxonomic marker and the filamenting temperature-sensitive mutant A-Z (FtsAZ) for cell division.

To confirm taxonomic assignments and metabolic activity of taxa involved in sulfate reduction and anaerobic oxidation of methane (AOM) processes, we performed a phylogenetic analysis of the RNA polymerase sigma factor (rpoD) and ribosomal protein S3 (rps3) gene proteins, aprAB, dsrAB, mcrA, and hydB, that we could detect in the transcriptomes (Vuillemin et al., 2020a). All ORFs annotated to these genes were aligned against their top two BLASTp hits in the NCBI database, using MUSCLE (Edgar, 2004). Conserved regions of the alignments were selected using Gblocks 0.91b (Castresana, 2000) with the following settings: allowing smaller final blocks, gap positions within the final blocks, and less strict flanking positions (Dereeper et al., 2008). Phylogenetic analysis of the resulting amino acid alignments of the predicted proteins were conducted in SeaView version 5.0.5 (Gouy et al., 2010), using PhyLM maximum likelihood (Guindon et al., 2010), with BLOSUM62 as the evolutionary model and 100 bootstrap replicates.

3 Results

3.1 Pore water geochemical profiles

All sulfate concentrations in the uppermost profiles roughly match the sulfate concentration of the Barents Sea bottom waters (ca. 28 mM). In all cores, sulfate concentrations decrease almost linearly with depth, with linear regressions showing correlation coefficients of R²>0.85 (Figs. 2a–b and S1 in the Supplement). The slopes of the regression lines can be divided into two groups corresponding to one group encompassing reference sites Ref 5 and Ref 6 and another group including the three HC-affected sites Zone 1, Zone 10, and Zone 12. Based on a t-test using the software PAST v. 3.20 (Hammer et al., 2001), the slopes of these two groups are statistically significantly different (t=5.574, $p=\mathrm{1.25}\times {\mathrm{10}}^{-\mathrm{6}}$) from each other (Fig. 2). The sulfate profiles from HC-affected sites show a steeper decrease with depth than those from reference sites. On average, sulfate concentrations decrease by 10 mmol m⁻¹, which would correspond to complete sulfate depletion at roughly 3 m b.s.f. at HC-affected sites. At reference sites, the average slope is 3 mmol m⁻¹ with few variations, which implies that sulfate would be depleted at ca. 12 m b.s.f.

Figure 2 Pore water linear regression for sulfate concentrations, alkalinity, and sulfide concentrations, along with some representative profiles. (a, b) Regression lines for sulfate profiles (derived from measured pore water concentrations) are linear (R²>0.85) and exhibit decreasing trends with depth. The regression lines shown in panel (a) represent HC-affected sites (Zone 1, Zone 10, Zone 12), which display a steeper slope compared to those in the panel (b) reference sites (Ref 5, Ref 6). (c, d) Regression lines for pore water alkalinity (calculated from measured pore water concentrations) are linear (R²≥0.85) and increase with depth. The regression lines for HC-affected sites in panel (c) show a more gradual increase compared to the reference sites in panel (d). (e) Sulfide concentrations measured across all cores from HC-affected sites (represented by dots), with their corresponding regression lines shown where sulfide concentrations were quantifiable. (f) Subset of representative profiles selected for each sampling site. The uppermost and lowermost rows depict cores with the highest and lowest R² values, respectively, from HC-affected and reference sites.

All alkalinity profiles display a clear linear increase with depth (R²≥0.85). From 5 mmol L⁻¹ at the SWI, pore water alkalinity increases differently with depth, corresponding to 8, 5, and 6 mmol L⁻¹ m⁻¹ at HC-affected sites Zone 1, Zone 10, and Zone 12 and ≤ 1 mmol L⁻¹ m⁻¹ at references sites Ref 5 and Ref 6, respectively (Figs. 2c–d and S2).

Based on the slopes of pore water profiles, the calculated alkalinity fluxes are larger at HC-affected sites than at reference sites. The same is true for sulfate fluxes, but differences across HC-affected sites and reference sites are less pronounced (Fig. S3). In comparison, variations observed within each sampling site are small.

Sulfide could hardly be detected in pore water samples collected from reference sites, most values being below 1 µM. In contrast, sulfide was detected in most of the cores from HC-affected sites. Especially at Zone 1 and Zone 10, sulfide concentrations increase linearly with depth (Fig. 2e). The majority (ca. 70 %) of the cores from Zone 1 show an increase of 2.5 mM m⁻¹ in sulfide concentrations. One core (Z01C03) stands out, reaching sulfide concentrations of ca. 11 mM at 1 m b.s.f. Several cores (ca. 57 %) from Zone 10 show a binary distribution of sulfide concentrations, corresponding either to linear profiles reaching ca. 1 mM sulfide at 1 m b.s.f. or to a steeper increase reaching ca. 3.7 mM at the same depth. At Zone 12, sulfide is detectable in only 2 of the 17 cores, with a much smaller increase in concentrations compared to Zone 1 and Zone 10.

3.2 Measured and modelled sulfate reduction rates

SRRs were detectable in only 29 out of the 50 cores (Fig. 3), with rates between 4 and 320 pmol cm⁻³ d⁻¹. In 403 of the 455 samples across the 50 cores analysed (i.e. 89 %), SRRs were below the detection limit (1 pmol cm⁻¹ d⁻¹). The occurrence of very low but detectable SRRs increased with sediment depth, which mostly stands true for HC-affected sites (Fig. 3).

Figure 3 Sulfate reduction rates (SRRs) and percentages of samples with sulfate reduction rates above the detection limit. Low SRRs became increasingly detectable with increasing sediment depth. The number of samples with measurable SRRs (% grey bars) is higher in sediments from HC-affected sites (dots) compared to reference sites (triangles); NB due to a lack of samples, the two uppermost intervals were merged.

Diffusion–reaction modelling with PROFILE (Berg et al., 1998) suggests steady SRRs in 47 of the 50 cores (Fig. S1), implying that SRRs remain constant over the entire length of the cores. The modelled SRR profiles show SR activity (i.e. > 0 pmol cm⁻³ d⁻¹) in 33 cores, of which 6 display rates ≤ 0.01 pmol cm⁻³ d⁻¹ and 27 have a supposed SRR ≥ 1 pmol cm⁻³ d⁻¹. For three specific cores (i.e. Z01C02, Z01C09, Z12C17), PROFILE proposes two intervals with different SRRs, suggesting that SR processes may only be active in the upper 50 cm b.s.f. and inactive in the sediment below.

The modelled and measured SRRs are indicative of net and gross turnover rates, respectively. The gross rate reflects microbial turnover, whereas the net turnover rate includes the recycling of reduced sulfur compounds (Berg et al., 1998). The modelled (i.e. net) rates for most samples are higher than the measured (i.e. gross) rates, which is rather unusual and can be explained by very low values whose absolute changes have a drastic effect on the calculated values. The modelled and measured SRRs are nevertheless within the same order of magnitude and are thus comparable.

3.3 Methane and carbon dioxide concentrations

All samples contain measurable methane concentrations, ranging from 0.08 to 19.79 µmol L⁻¹. All samples from reference sites contain less than 5 µmol L⁻¹ of methane, with detection of any higher concentrations being restricted to some samples from HC-affected sites and only below 100 cm b.s.f. (Fig. 4). CO₂ was detectable in all samples, with concentrations ranging from 29 to 512 µmol L⁻¹, but without any specific trend with depth (Fig. 4). The reason is probably that any CO₂ from whichever source enters the carbonate equilibrium and the gas is removed from the pore water, both in HC-affected sites and reference sites. The low methane and carbon dioxide concentrations prevented any measurement of their stable isotopic signatures.

Figure 4 Concentrations of dissolved methane and carbon dioxide. (a) At HC-affected sites, methane concentrations are highest in samples originating from 100 cm b.l.f. or deeper, whereas methane concentrations at reference sites remain below 5 µmol L⁻¹. (b) CO₂ concentrations across samples do not display any specific trend with sediment depth.

3.4 Statistical correlations

To identify geochemical differences between HC-affected and reference sites, we performed a Pearson correlation analysis. Correlation coefficients calculated for each separate gravity core (Table S2) indicate that alkalinity is strongly negatively correlated ($p\ge -\mathrm{0.85}$) with sulfate concentrations in sediments from HC-affected sites. The same applies to quantifiable sulfate and sulfide concentrations ($p\ge -\mathrm{0.7}$). For reference sites, alkalinity and sulfate concentrations are also strongly correlated (p≈0.8), whereas sulfide concentrations are below the detection limit and have hence been omitted for the correlation analysis.

We established ratios of sulfate to alkalinity (HCO${}_{\mathrm{3}}^{-}$) fluxes (Fig. 5). As about 90 % of seawater alkalinity can be contributed to bicarbonate (89.8 % HCO${}_{\mathrm{3}}^{-}$, 2.9 % CO${}_{\mathrm{3}}^{\mathrm{2}-}$) (Kerr et al., 2021), we set alkalinity synonymous with bicarbonate concentration. We compared the measured fluxes to those expected from organoclastic (1:2 based on Eq. 1) and methanotrophic (1:1 based on Eq. 2) SRR activities. At reference sites, the average ratio of sulfate to alkalinity fluxes equals 2:1, whereas it is 1.25:1 (5:4) at HC-affected sites, implying either that HCO${}_{\mathrm{3}}^{-}$ is being removed from the pore water or that pore water sulfate is eventually being replenished via sulfide oxidation.

Figure 5 Relationship between sulfate and alkalinity fluxes at HC-affected sites and reference sites. The majority of samples from HC-affected sites (dots) have lower alkalinity than sulfate fluxes. This is even more pronounced at reference sites (triangles). Black lines signify theoretical ratios for methanotrophic (1:1) and organoclastic (1:2) SR.

3.5 Functional marker genes specific to the SMTZ

From the 12 sequenced metagenomes, we extracted a total of 307 234 ORFs with multiple annotations that corresponded to 821 850 predicted genes. From the 5 sequenced metatranscriptomes, we extracted a total of 31 093 ORFs corresponding to 47 316 predicted genes (Table S2).

For each of the targeted marker genes, we plotted the relative percentage of total predicted genes and their corresponding taxonomic assignments at the phylum and class level (Fig. 6). Taxonomic assignments of SSU rRNA and FstAZ proteins show that metabolically active and growing microbial populations at HC-affected sites are affiliated with the phyla Chloroflexota, Desulfobacterota, and Euryarchaeota and more specifically with anaerobic methanotrophs (i.e. ANMEs). The phyla Desulfobacterota and Chloroflexota are actively engaged in processes of sulfate reduction, as shown by the expression of ORFs encoding the aprAB and dsrAB genes (Fig. 7a–b). Taxonomic assignments of expressed ORFs encoding mcrA genes show an increased abundance of methanotrophic ANMEs (Fig. 7c) over methanogenic Halobacterota. ORFs encoding hydB-G genes are expressed by taxa related to the classes Lokiarchaeia and Thorarchaeia (Fig. 7d) among the phylum Asgardarchaeota, indicating metabolic capability to use reduced sulfur compounds (e.g. polysulfide, elemental sulfur). However, their overall expression levels appear limited, and these Asgardarchaeota do not seem to actively divide (Fig. 6). At reference sites, Chloroflexota appear to be the only phylum capable of growing and dividing in these organic-poor sediments (Fig. 6). However, concomitant expression of functional marker genes involved in sulfate reduction remains below the limit of detection. We acknowledge that the non-detection of certain genes may also result from limited sequencing depths.

Figure 6 Open reading frames encoding functional marker genes involved in sulfate reduction and production and anaerobic consumption of methane and their taxonomic assignments at the phylum to class level. The relative abundances of certain phyla increase at HC-affected sites compared to reference sites. Concomitant expression of apr and dsr genes by Desulfobacterota, phs genes by Chloroflexota, hyd genes by Asgardarchaeota, and mcr genes by ANMEs and Halobacterota are consistent with methanotrophic sulfate-reducing activities. apr: adenylylsulfate reductase; dsr: dissimilatory sulfite reductase; asr: anaerobic sulfite reductase; phs/psr:  thiosulfate and polysulfide reductase; hyd: sulfhydrogenase; mcr: methyl coenzyme M reductase; MG: metagenomics; MT: metatranscriptomics; BDL: below detection limit.

Phylogenetic trees plotted for transcript ORFs assigned to the rps3 and rpoD taxonomic marker genes (Figs. S4–S5) provided concomitant evidence for metabolic activity by the phyla Desulfobacterota, Chloroflexota (class Dehalococcoidia), Acidobacteriota (class Aminicenentia), Asgardarchaeota (class Thorarchaeia), and Halobacterota (class ANME-1).

4 Discussion

4.1 Pore water geochemical profiles correlate to meta-omics features

Pore water profiles reflect the spatial distribution of net turnover processes. The shape of the profile allows us to distinguish intervals in which transport is purely diffusive (i.e. linear profile) from those in which there is either net production (i.e. convex profile) or consumption (i.e. concave profile) of chemical species (Schulz, 2006; Hesse and Schacht, 2011).

Based on the largely linear profiles for pore water alkalinity, sulfate, and sulfide (Figs. 2–3, S1–S2), we infer the absence or near-absence of net sulfate consumption in the depth interval covered by the cores, a sink for sulfate, and a source of sulfide and alkalinity at greater depth. Microbially mediated AOM with sulfate as an electron acceptor yields hydrogen sulfide, bicarbonate, and water (Eq. 2), providing an explanation for both sink and source (Martens and Berner, 1974; Iversen and Jorgensen, 1985), assuming sufficient supply of methane at greater depth. The occurrence of AOM at deeper depth is further supported by the relationship between sulfate and alkalinity fluxes (Fig. 5) and negligibly small methane concentrations (< 10 µM L⁻¹) at the bottom end of all cores from HC-affected sites and reference sites (Fig. 4), indicating almost complete consumption of methane originating from biogenic or thermogenic sources below our sampling interval.

Linear concentration profiles are only obtained in regions with no net consumption or production of the respective compound (Schulz, 2006). Despite the apparent linearity of sulfate profiles (Fig. 2), we observed very low SRRs (≤ 300 pmol cm⁻³ d⁻¹), confirming previous reports of sporadic SRRs in Barents Sea sediments with similar rates above the SMTZ (Nickel et al., 2012). In our case, sulfate reduction is mainly detectable at HC-affected sites (Fig. 3), with a greater prevalence in deeper sediments (Fig. 5). With seepage supplying electron donors, like methane and/or other HCs, from deeper reservoirs to near-surface sediment, sulfate reduction is fuelled by those HCs rather than from the buried and more recalcitrant OM. Inconspicuous supplies of methane and/or other HCs to the overlying sediment can stimulate SRRs by providing alternative electron donors. Such sporadic energy supply promotes geochemical conditions that stimulate benthic HC-degrading microbial communities. Although these communities can be physiologically and phylogenetically diverse, they are mainly composed of anaerobic guilds of methanogens, methanotrophs, sulfate reducers, and fermenters (Teske, 2019).

The analysis of functional marker genes supports the notion of sporadic HC supply, providing complementary insights into the metabolic potential and activity of specific microbial populations in the sediment. The occurrence and abundance of aprAB, dsrAB, and mcrA genes supports interpretations of geochemical data in terms of sulfate-driven AOM processes and metabolic activity (Fig. 6). Taxonomic profiling of populations expressing ORFs encoding genes involved in sulfate reduction (i.e. aprAB, dsrAB) reveals Desulfobacterota (family Desulfobacteraceae) and Chloroflexota (class Dehalococcoidia) as the prevalent bacterial phyla constituting the consortium of sulfate-reducing bacteria at HC-affected sites (Figs. 7, S4–S5). The aprAB and dsrAB phylogenetic trees clearly show that none of these genes correspond to their reverse function, and the ORFs annotated as aprAB and dsrAB are predicted to function in dissimilatory sulfate reduction. Moreover, with the sediment being fully anoxic, any detection of sulfur oxidation via reversible aprAB or reverse dsrAB was unexpected, as such activities are mostly reported for aerobic Gammaproteobacteria (Vuillemin et al., 2022).

Figure 7 Phylogenetic trees of aprAB, dsrAB, mcrA, and hydB-G protein-encoding genes. PhyML maximum likelihood trees of transcript ORFs encoding (a) the adenylylsulfate reductase (714 amino acids) subunit alpha (aprA) and beta (aprB); (b) the dissimilatory sulfite reductase (493 amino acids) subunit alpha (dsrA) and beta (dsrB); (c) the methyl coenzyme M reductase (698 amino acids) subunit alpha (mcrA); and (d) the sulfhydrogenase, or sulfur reductase (286 amino acids), subunit beta (hydB) and gamma (hydG). Phylogenetic trees are based on 100 bootstrap replicates with BLOSUM62 as the evolutive model. Boldface type signifies sequence accession numbers to the NCBI database

Furthermore, ORFs encoding genes involved in anaerobic production and consumption of methane (i.e. mcrA) clearly show that the clade ANME-I exhibits concomitant, although low, metabolic activity towards anaerobic oxidation of methane (Fig. 7c); namely they express ORFs encoding the reversible mcrA gene. The higher expression levels of methanotrophic clade ANME-I over methanogenic Halobacterota (Fig. 6) also suggest that metabolic activity by ANMEs could be sustained by an alternative source of methane (Dong et al., 2020), although cryptic production of biogenic methane by the ANME-1 clade has also been observed in experiments replicating SMTZ conditions (Beulig et al., 2019). Because methane concentrations measured at both HC-affected sites and reference sites were similarly low (Fig. 4), it appears unlikely that AOM processes are sustained by in situ production of biogenic methane only (Dong et al., 2020). Finally, the classes Lokiarchaeia and Thorarchaeia in the phylum Asgardarchaeota expressed ORFs encoding beta and gamma subunits of hydB-G genes (Fig. 7d), and thereby appear metabolically active, using various sulfur and sulfide species, potentially produced as by-products of sulfate reduction.

In contrast, Chloroflexota (class Dehalococcoidia) was the only phylum at reference sites showing sufficient metabolic activity adequate to achieve cell division (i.e. FstAZ). In the absence of a univocal signal towards active dissimilatory sulfate reduction by Dehalococcoidia at HC-affected sites (Fig. 7b), we infer that these Chloroflexota could persist and gain energy through fermentation processes in organic-lean marine sediments instead (Vuillemin et al., 2020a). Due to this metabolic versatility, Chloroflexota were found to be active at both HC-affected and reference sites, which ranks them as poor microbial indicators for minor HC seepage (Iasakov et al., 2022).

Together, these microbial patterns extracted from meta-omics corroborate the geochemical data, arguing in favour of discrete AOM processes concomitant to sulfate reduction at HC-affected sites and for limited in situ production of biogenic methane at both HC-affected sites and reference sites. Alongside an active AOM consortium composed of Desulfobacterota and ANME-1, concomitant production of reduced sulfur compounds could apparently stimulate metabolic activity by Lokiarchaeia and Thorarchaeia (Figs. 6 and 7).

4.2 Geochemical indicators for minor HC seepage

The combined use of meta-omics microbial patterns with the determination of SR and related geochemical features and fluxes in upper sediments enables us to distinguish between HC-affected sites and reference sites to a certain extent. Sediments of the SW Barents Sea are organic-poor (TOC < 0.5 %) and thus characterized by a lack of electron donors (Knies and Martinez, 2009) and low microbial activity (Nickel et al., 2012, 2013) with consequent sulfate penetration to greater sediment depths. In organic-rich sediments, any signal inherent to minor HC seepage fuelling SR would be difficult to detect, as it would be covered by organoclastic SR processes. The organic-lean nature of Barents Sea sediments should enable the detection of discrete HC seeps because the slightest supply of electron donors leads to a direct increase in the metabolic activity of HC-degrading microbial populations. This slight increase in metabolic activity can then be discriminated from the background signal observed at reference sites. Diffusion of HCs was thus expected to promote sulfate consumption, resulting in its depletion in the pore water at shallower sediment depths. This, combined with variable fluxes of thermogenic and biogenic methane, can account for variations in the depth of the SMTZ.

At different locations of the SW Barents Sea, the depth of the SMTZ, or rather the sulfate depletion depth, was postulated, for instance, at 37 m b.s.f. in the Hammerfest/Loppa High and the Tromsø Basin/Ingøydjupet area (Nickel et al., 2012) and between 3.5 and 29.2 m b.s.f. in the Ingøydjupet area (Argentino et al., 2021a). These values are consistent with the global average depth of the SMTZ (12.8 ± 12.1 m b.s.f.) on the continental slope (Egger et al., 2018). For comparison, at the Vestnesa Ridge, which is under the influence of gas hydrates, the depth of the SMTZ fluctuates massively between seepage (ca. 1 m b.s.f.) and non-seepage (ca. 3 to 75 m b.s.f.) sites (Hong et al., 2016). Similarly, Barents Sea sediments with and without methane seeps show an SMTZ located between a few centimetres and about 2 m b.s.f. (Argentino et al., 2021b).

Linear extrapolation of sulfate profiles (Fig. 2) indicates complete sulfate depletion, indicative of an SMTZ, at 12 and 3 m b.s.f. for reference and HC-affected sites, respectively. We suggest that the difference in sulfate penetration depth is caused by upward migration of HCs. The sulfate and alkalinity fluxes inherent to organoclastic (Eq. 1) and methanotrophic (Eq. 2) SR have respective theoretical ratios of 1:2 and 1:1, while mixed SR processes have intermediate values. Current ratios of sulfate to alkalinity fluxes (Fig. 5) at reference sites imply either that HCO${}_{\mathrm{3}}^{-}$ is being removed from the pore water or that pore water sulfate is eventually being replenished via sulfide oxidation. From other fieldwork in the area, we know that the sediment is fully anoxic below the uppermost few centimetres. Phylogenetic analysis of aprAB and dsrAB genes clearly shows that none of the identified amino acid sequences correspond to their reverse function. Furthermore, we did not identify any dsrC sequences, which potentially work in cellular sulfur oxidation (Ruiz-Blas et al., 2024). While dissolved inorganic carbon (DIC) can be assimilated via acetogenesis, methanogenesis, or methanotrophy (Kellermann et al., 2012), authigenic carbonate minerals are known to precipitate from the pore water as by-products of AOM processes (Turchyn et al., 2021), sometimes even forming carbonate hardground at the seabed. However, from the low level of detection of AOM-related transcripts, we inferred only a minor imprint on the geochemistry and therefore assume that bicarbonate ions originating from AOM processes at greater depths are diffusing upward and eventually precipitate as calcium carbonate. Up to 80 % of the sulfide produced through SR can also be converted back into sulfate through reoxidation of reduced sulfur compounds (Jørgensen et al., 2004) and disproportionation (Pellerin et al., 2015; Jørgensen et al., 2019a). Finally, the dissolution of sulfate minerals, such as barite (BaSO₄), may also act as an additional source of sulfate (Griffith and Paytan, 2012). Altogether, this suggests either the reoxidation of sulfide or the removal of bicarbonate from the pore water, e.g. by mineral precipitation, or both. The former is further supported by concomitant sulfate reduction, detected by radiotracer incubations (Fig. 3) and microbial gene expression patterns (Fig. 6), the latter by the ratios of sulfate to alkalinity flux (Fig. 5). Thus, we could identify several biogeochemical features that argue for ongoing sulfate reduction (most likely AOM-related) at HC-affected sites below our sampled depth range and which allow their discrimination from reference sites. Given the sulfate concentration profiles and the small amounts of methane detected in the samples, it becomes clear that AOM is not a quantitatively important process in the cored depth interval, producing an erratic signal in HC-affected sediments.

5 Conclusions

Our study, combining modelled and measured SRR and pore water geochemistry with functional marker genes, revealed remarkable differences between samples originating from a pristine seafloor and those from areas affected by small and inconspicuous HC seepage. The linear shape of pore water profiles was indicative of purely diffusive transport under steady-state conditions and allowed estimation of the depth of sulfate depletion, which appeared markedly different between HC-affected and reference sites. At the HC-affected sites, we found higher pore water fluxes of sulfate, sulfide, and alkalinity and could also demonstrate active microbial SR by turnover measurements and modelling. Furthermore, meta-omics provided additional evidence for enhanced metabolic activity related to microbial sulfate reduction coupled with minor AOM processes at HC-affected sites, which could be clearly discriminated from the background signal inherent to fermentative activities and minor in situ production of methane by the subseafloor biosphere. These observations demonstrate that minor HC seepage, which is not even visible at the seafloor, significantly influences sedimentary biogeochemical cycles driven by microbial populations resident in organic-lean sediments.