The study was performed in Aarhus Bay where organic-rich Holocene marine mud has accumulated up to >10 m thick layers on top of the glacial deposits over the past 8000 years (Flury et al., 2016). The sedimentation rate at station M5 (56°06.20°N and 10°27.47°E; water depth 28 m), based on ¹⁴C-dating of mollusk shells, has been stable for the last 4300 years (Langerhuus et al., 2012). The water column is generally oxic, oxygen penetration depth in the sediment is low (1–5 mm; Rasmussen and Jørgensen, 1992), and bioturbation is visible from sediment color and faunal tracks in the top 3 cm. Sediments for incubation were sampled in December 2015 and May 2016 with a Haps corer (Kanneworff and Nicolaisen, 1972), then transferred into a gravity core liner (diameter 11.8 cm) on deck, brought back to the laboratory, kept at in situ temperature, and processed within 48 h. The temperature of the bottom water, measured by thermometer, was around 10°C at both sampling times and the salinity was 30‰, determined by a handheld refractometer. Rumohr Lot cores (Meischner and Rumohr, 1974) were taken in parallel for geochemical profiles and molecular analysis of the methanogen community, and the sampling intervals were 2 cm in the top 20 cm. Sediment samples for molecular analysis were taken with sterile 5 ml cut-off plastic syringes and preserved at -80°C. Pore-water samples for sulfate and dissolved inorganic carbon (DIC) analysis were extracted from the intact sediment core via Rhizon suction samplers (0.1 μm porous polymer, Rhizosphere Research, Wageningen, Netherlands), pretreated as described previously (Glombitza et al., 2014). Pore-water samples for DIC were kept at 4°C before analysis. For sulfate, the extracted pore-water was flushed with CO₂ to remove hydrogen sulfide, and frozen at -20°C until analysis. For methane sampling, 3 ml of water close to the sediment surface (< 1 cm) was taken with a 20 ml syringe mounted with a long PVC tube, and 2 ml of sediment was collected with cut-off plastic syringes. All samples were transferred immediately to 10 ml glass serum vials (preloaded with 4 ml distilled water and 1.5 g NaCl), which subsequently were closed with a butyl rubber septum and a crimp cap. The serum vials containing sediment were shaken, and preserved at -20°C before analysis.

Diffusive fluxes of CH₄ were estimated from concentration gradients using Fick’s first law:

where J is the diffusive flux of CH₄, C is the concentration, z is depth (cm), and φ is porosity (Boudreau, 1997). D is the diffusion coefficient in free solution for a salinity of 30 and temperature of 10°C: D_CH4 = 1.02 cm² d^-1 (Dale et al., 2008). The factor 1 - ln(φ²) corrects for the lower diffusion coefficient in sediment due to tortuosity (Boudreau, 1997).

Plastic bags made of gas-tight multi-laminar plastic film (Supplementary Figures S1a, b) (Amcor, Denmark) were used for sediment incubation and ¹³CH₄ spiking experiments. This bag incubation method was developed by Hansen et al. (2000), but some modifications were made: (1) the bags were sealed with a black slide binder (LEITZ, Germany) on one side, (2) the silicon gaskets were replaced by butyl rubber gaskets to improve gas tightness (Supplementary Figure S1c), (3) the plastic bags were mounted with a glass outlet custom made to fit the cut-off 2.5 ml syringes (Supplementary Figure S1d), (4) to physically protect the incubation bags and minimize oxygen penetration, they were placed in an outer gas-tight bag filled with N₂. Prior to the isotope spiking experiments, gas-tightness tests were conducted by filling bags with artificial sea water brine of 300‰ salinity (for inhibiting microbial activities). CH₄-saturated artificial sea water was added to two parallel bags, leaving no headspace, and incubated at 15°C. Triplicate samples were taken by syringes with a needle through the stopper and transferred into 3 ml glass vials. Methane in the water samples was measured by the headspace method using gas chromatography (SRI Instruments, United States). Gas-tightness tests were conducted by monitoring CH₄ concentration in the bags over time. Within 140 days, methane concentrations in the bag were very stable with little fluctuations (ANOVA, P > 0.05; Supplementary Figure S2), suggesting that this type of bags can be used for incubation experiments over many days or even months. By this incubation method, sediment is kept in closed bags free of air and headspace and without dilution or stirring. This maintains an environment close to the in situ conditions (Hansen et al., 2000), which is critical for the study of natural process rates.

The following manipulations were all done under strict anoxic conditions in a glove box (Plas-Lab, Inc., United States), which was continuously flushed with nitrogen gas. Oxygen concentration in the glove box was <0.5 μmol l^-1, as monitored with an oxygen sensor (Revsbech, 1989). Sediments from Haps cores were sliced into six 2 cm segments, five in the top 10 cm, and one around 15 cm. The sliced segments were passed through a 2 mm sieve to remove macrofauna and hard objects such as carbonate shells which could puncture the bag, and then loaded into the gas-tight bag, carefully avoiding headspace and removing gas bubbles. After 2–3 days of stabilization, ¹³C-labeled methane was introduced into each bag by injecting 20–100 μL of ¹³CH₄-containing artificial sea water (excluding sodium sulfate), which resulted in a ¹³C isotope enrichment of ∼2% with minimal change (∼2%) of the total CH₄ concentrations in pore-water. The sediment was homogenized by manual kneading of the bag for 20 min (Hansen et al., 2000). Six bags with sediment were incubated at 10°C in a thermostated incubator (MMM Incucell, Germany). At each sampling time, triplicate samples were taken by squeezing the bag to push the sediment through the glass outlet, where a PVC tube connected the outlet to a cut-off 2.5 ml syringe to ensure tightness (Supplementary Figure S1d). Sediment collected into the syringe was immediately transferred into 10 ml glass serum vials, which were closed with a butyl rubber septum and crimp cap, and preserved at -20°C, as described above.

Methane production and oxidation were calculated from the CH₄ concentration and its isotope composition using an isotope dilution model, which was originally developed to determine rates of NH₄⁺ turnover in anoxic marine sediments (Blackburn, 1979). The model here presumes that the concentration of CH₄ changes during the experiment through time due to constant methane consumption (i.e., oxidation) (r) and production (p), thus:

The model further presumes that the mole fraction of ¹³CH₄ [R_t = ¹³CH₄/(¹³CH₄ + ¹²CH₄)] changes trough time as a result of production of CH₄ from indigenous sources and consumption of CH₄ due to anaerobic oxidation in the sediment. In case the rate of methane production, p, is different from rate of methane consumption, r, the relationship between R_t and time can be described with the following expression:

where R_b is the mole fraction of ¹³CH₄ in CH₄ being produced from indigenous sources in the sediment and R₀ is the mole fraction of ¹³CH₄ at the start of the experiment (see SI for the deduction of Eq. 2, and the equation when p = r). Equation 1 shows that the quantity “p - r” equals the slope of a regression line obtained from plotting C_t against, while Eq. 2 shows that the quantity “-p/(p - r)” equals the slope of a regression line obtained from plotting ln(R_t-R_b) against ln(C_t/C₀). Therefore p and r can each be calculated from the combination of these slopes.

The mineralization rates were measured by two methods, (1) during both cruises, sulfate reduction rates (SRRs) were determined the using the whole-core injection method (Jørgensen, 1978). Two microliters of carrier-free ³⁵S-sulfate solution containing approximately 100 kBq (ca. 2.5 μCi) of radioactivity, were injected in 2 cm intervals at multiple depths in 26-mm-i.d. core. The formation of radiolabeled products was analyzed after incubation for 12–19 h at the in situ temperature. Incubations were terminated by sectioning the core into 2 cm segments that were fixed in weighed vials containing 10 ml of cold 5% Zn-acetate. The sediment was immediately mixed and stored at -20°C up to 4 weeks until further analysis. After thawing, the total reduced inorganic sulfur (TRIS = H₂S, S₀, FeS, and FeS₂) was separated from the ³⁵SO₄^2- by a single step cold chromium distillation (Kallmeyer et al., 2004). The radioactivity of ³⁵SO₄^2- and ³⁵S-TRIS was counted in 15 ml scintillation cocktail (Gold Star, Meridian Biotechnologies, United Kingdom) on a TriCarb 2900TR liquid scintillation analyzer (Packard Instrument Company, Germany). The ex situ SRRs were then calculated accordingly with a correction factor of 1.06 (Jørgensen, 1978; Flury et al., 2016); (2) In May 2016, accumulation of ammonium (NH₄⁺) and DIC were also monitored in the bags to estimate mineralization rates. Sediment subsamples of 2.5 ml were taken from each bag and centrifuged, NH₄⁺, DIC, and calcium were analyzed in the supernatant. The Ca²⁺ concentrations in sediment during the incubation were quite stable with little change over time (ANOVA, P > 0.05). The production rates of NH₄⁺ and DIC were calculated by plotting concentrations against time and using linear fitting, and production rates of NH₄⁺ were multiplied by a factor of 1.3 to correct for adsorption, according to Canfield et al. (1993a). The rates of organic carbon mineralization rates from SRR were calculated with a 1:2 stoichiometry between SRR and organic C oxidation (Flury et al., 2016).

Methane was measured using a Gas Chromatograph (GC) (310C, SRI Instruments, United States) with a flame ionization detector. The carbon isotope composition of methane was determined using a coupled pre-concentration GC/Isotope Ratio Mass Spectrometer system (GC/IRMS) (Thermo Fisher Scientific, Germany) (Rice et al., 2001). The δ¹³C values were reported vs. the Vienna Pee Dee Belemnite standard (VPDB). DIC was measured on the same GC/IRMS system without pre-concentration (Torres et al., 2005). Sulfate was measured by an IC-2500 ion chromatography system (Dionex Corporation, United States), NH₄⁺ was determined by a FLUOstar Omega spectrophotometer (BMG LABTECH, Germany) based on the salicylate-hypochlorite method (Bower and Holm-Hansen, 1980), and calcium was measured by inductively coupled plasma optical emission spectrometry (ICP-OES) (Perkin Elmer Optima 2000DV, United States).

Aliquots of frozen sediment samples (0.2 g) from the parallel Rumohr Lot core were used for total DNA extraction. Sediment samples were thawed at room temperature in a lysis buffer mixture containing 0.65 ml sodium phosphate buffer solution (112.9 mM Na₂HPO₄, 7.1 mM NaH₂PO₄) and 0.2 ml SDS solution (500 mM Tris–HCl, 100 mM NaCl, 10 w.-% sodium dodecyl sulfate, pH 8.0), as well as 0.25 ml zirconia beads (0.1 mm diameter, BioSpec, United States). The thawed mixture was subjected to bead beating at 50 oscillations s^-1 for 1 min using a TissueLyser LT 2500 (Qiagen, Germany), followed by incubation in a thermomixer with 600 rpm at 50°C. After lysis, the mixture was centrifuged for 10 min at 19,000 × g at 4°C. The supernatant was extracted with equal volumes of phenol:chloroform:isoamyl alcohol (25:24:1, vol:vol:vol; Sigma Aldrich, United States), followed by extraction with equal volumes of chloroform: isoamyl alcohol (24:1, vol:vol; Sigma Aldrich, United States). DNA was precipitated with 1 ml polyethylene glycol 8,000 (Sigma Aldrich, United States) at 4°C overnight, centrifuged at 19,000 × g for 30 min. The precipitates were washed with ice cold 70 vol% ethanol solution, dried up in the air and then dissolved in 400 μl TE buffer (10 mM Tris, 1 mM EDTA, pH 8.0) and stored at -20°C before use.

Abundances of archaeal 16S rRNA genes were determined on LightCycler 480 (Roche, US) using the primer sets Arch 806F (ATT AGA TAC CCS BGT AGT CC) (Takai and Horikoshi, 2000) and Arch 958R (YCC GGC GTT GAM TCC AAT) (DeLong, 1992). A fragment covering the hypervariable regions (V3–V6) of the Archaeal 16S rRNA gene were PCR amplified using the primers Arch344Fmod (5′-GGGYGCAGCAGKCGMGAA-3′) and Arch915R (5′-GTGCTCCCCCGCCAATTCCT-3′) (Goodfellow and Stackebrandt, 1991). The forward primer was modified from the initial primer Arch349F (5′-GYGCASCAGKCGMGAAW-3′; Takai and Horikoshi, 2000) by extending the 5′-end and reducing degeneracy without losing coverage as determined by in silico tests using TestPrime 1.0 (Klindworth et al., 2013). PCR was carried out in 25 μl reaction volume with GeneAmp 9700 Thermal Cycler (Applied Biosystems, United States) with the following cycling parameters, initial pre-denaturation at 98°C for 2 min, followed by 25 cycles of denaturation at 95°C for 30 s, annealing at 60°C for 45 s, elongation at 72°C for 1 min, plus a final extension at 72°C for 10 min. The PCR mixture contained 0.125 μl TaKaRa ExTaq Hot Start polymerase (TaKaRa, Japan), 2.5 μl of 10× Ex Taq buffer, 2 μl dNTP solution (2.5 mM), 1 μl of each primer (10 μM), 2 μl BSA (10 μg μl^-1), 1 μl of template DNA, and 15.375 μl de-ionized ultrapure water. In a second PCR, all PCR products were supplied with the forward and reverse Illumina^® adapter overhang sequences, which are provided in Illumina^® published “16S metagenomic sequencing Library preparation protocol”¹. The resultant PCR products were purified and supplied with indices using the Illumina^® Nextera XT Index Kit. Negative controls without template were included to test for reagent contamination in each set of PCR reactions. PCR products were evaluated on 2% agarose gels, and purified using the Agencourt AMPure XP Kit (Beckman Coulter, Inc., United States) before used as template for the next PCR. The final purified PCR products were quantified by a Qubit 2.0 fluorometer with the dsDNA HS Assay Kit (Life Technologies, United States), and then pooled in an equimolar ratio. The pooled library was sequenced on an Illumina^® MiSeq system using a 600 cycle MiSeq v3 Reagent Kit (Illumina, United States), which produces two 300-bp long paired-end reads. An average of 39,026 (ranging from 16,477 to 66,068) reads were generated for all samples, and more than 99.7% of them classified as archaea, thus indicating good specificity of the primers.

The paired-end MiSeq reads were processed using MOTHUR 1.36.1 (Schloss et al., 2009) by following the MiSeq SOP pipeline (Kozich et al., 2013). Merged reads with length shorter than 540 bp or longer than 599 bp, or homopolymeric stretches longer than 7 bp were all removed from the dataset in the initial quality filtering. Reads were further denoised by pre-clustering and the chimeric sequences were checked and filtered out by UCHIME (Edgar et al., 2011). The filtered sequences were aligned and classified based on the SILVA SSU REF NR v.123 database (Quast et al., 2013; Yilmaz et al., 2014). Sequences not aligning within the expected region flanked by the PCR primers were removed from further analysis. Methanogen taxonomic lineages were classified into different physiological categories depending on their presumed substrate preference (Canfield et al., 2005; Jabłoński et al., 2015): hydrogenotrophic (H₂) using H₂/CO₂, acetoclastic (Acetate) using acetate, methylotrophic (Methyl) using non-competitive methylated substrates (MAs, methanol, etc.), mixotrophic (Mix) using more than one type of substrates, and unclassified (N). The sequences generated in this study were deposited in the NCBI Sequence Read Archive (SRA) under project no. PRJNA359637.

Profiles of geochemical parameters showed similar pattern in cores from the two cruises in December 2015 and May 2016 (Figure 1). Methane pore-water concentrations increased with depth, from 0.6 to 51 μM in the top 20 cm, where the sulfate concentration was >15 mM. There was a slight gradient of CH₄ from the sediment to the overlying seawater (Figure 1, enlarged section). The gradient was used to calculate diffusive methane fluxes across the sediment-water interface using Fick’s first law. The estimated methane flux was 0.22 nmol cm^-2 d^-1 during December and 2.13 nmol cm^-2 d^-1 during May (Table 1). This is a minimal estimate, as the effect of bioturbation was not taken into account. The sulfate–methane transition zone was shallow, at around 50 cm, but there was a long tailing of methane up through the sulfate zone. The δ¹³C of methane decreased with depth from -30‰ to -80‰ in both profiles, with a small shift to more negative values in the top 20 cm for the May samples.

The initial concentration of CH₄ in each bag did not increase with depth as it did in CH₄ profiles, which was due to loss of CH₄ during slicing and sieving. During incubations, the concentrations of CH₄ decreased with time for all depths except for the 0–2 cm (Figure 2). In this top section, CH₄ increased from 1.7 to 1.9 μM over the first 4 days, and then decreased with a similar trend as at other depths. This is consistent with the idea that that some active substrates were available in the top 0–2 cm for methanogenesis, which became depleted rapidly in a few days. Therefore, at 0–2 cm the modeling and calculation of methane turnover rates was separated in two phases, 0–4 days and 6–50 days. The δ¹³C values decreased with time in all incubations (Figure 2), as evidence of methane production from a non-spiked carbon source (Seifert et al., 2006). Similar patterns of CH₄ concentrations and δ¹³C values were displayed during bag incubations in May (Supplementary Figure S3).

As assumed by the isotope dilution model, the correlation coefficients were all above 0.96 by linearly fitting CH₄ concentrations with time (Figure 3). The validity of the model is more significantly corroborated by the linearity demonstrated by the plots of ln(R - R_b) against ln(C_t/C₀), and the correlation coefficients were all above 0.95. Similar trends and good fitting were also observed in experiments repeated in May 2016 (Supplementary Figure S4), and in the second phase of 0–2 cm (Supplementary Figure S5). Based on the modeling data, production and oxidation rates of CH₄ were calculated (Figure 4). Both production and oxidation rates peaked in the top layer with initial rates of over 200 pmol cm^-3 d^-1, and then decreased steeply to below 100 pmol cm^-3 d^-1 in all other layers, with no clear depth trend. Also during the second phase in 0–2 cm, rates were <100 pmol cm^-3 d^-1. The depth-integrated rates of CH₄ production and oxidation were slightly higher in December 2015 compared to May 2016 (Table 1), which were all around 1 nmol cm^-2 d^-1 in 0–16 cm.

The highest SRRs of up to 100 nmol cm^-3 d^-1 were detected at 0–2 cm in December 2015, and then decreased gradually to below 20 nmol cm^-3 d^-1 at 17.5 cm (Figure 5, left panel). In May 2016, a similar distribution of SRR with depth was observed, while the SRR were slightly higher at each depth compared to December 2015 (Figure 5, left panel). The rates of organic carbon mineralization rates calculated from SRR gave values close to production rates of DIC (Figure 5, right panel, and Supplementary Figure S6), except at 0–2 cm, where the DIC production rate (290 nmol cm^-3 d^-1) was higher than the mineralization rate calculated only from SRR. The depth-integrated rates of SRR and carbon mineralization were two to three orders of magnitude higher than both flux and turnover rates of CH₄ in 0–16 cm (Table 1).

The in situ archaeal communities and their depth distribution were similar at the phylum level during the two sampling occasions (Supplementary Figure S7), and dominated by Woesearchaeota, Euryarchaeota, Thaumarchaeota, and Bathyarchaeota (formerly known as the Miscellaneous Crenarchaeotal Group) (Meng et al., 2014). The relative abundances of Bathyarchaeota were constant with depth at around 10% of all Archaea, but no methanogens were classified into this phylum due to lack of taxonomy information at lower levels. Methanogenic archaea belonging to Euryarchaeota were present throughout the surface sediment, where their relative abundances (reads classified as methanogens/all reads) were generally less than 1.4% of all Archaea (Figure 6). Thus, the absolute abundances of methanogens, estimated by multiplying their relative abundances with total archaeal 16S rRNA gene copies, were two orders of magnitude lower than the total archaeal communities. Both the relative abundances and absolute abundances of methanogens decreased with depth, except at around 13 cm (Figure 6). Most methanogen-related sequences in the surface sediments belonged to unclassified Methanomicrobiales, Methanococcoides, and Methermicoccus (Supplementary Figure S8). The distribution of methanogens classified according to their inferred metabolic type is summarized in Figure 7, and hydrogenotrophic and methylotrophic methanogens dominated in all samples, while putative acetoclastic methanogens appeared only in a few samples with low abundances. A minor proportion of methanogens could not be assigned to a specific metabolic type.