Distinct methane-driven biogeochemical regimes in Arctic seafloor gas hydrate mounds

doi:10.21203/rs.3.rs-93738/v1

Archaea mediating anaerobic methane oxidation are key in preventing methane produced in marine sediments from reaching the hydrosphere; however, a complete understanding of how microbial communities in natural settings respond to changes in methane flux remains largely uncharacterized. We investigate microbial communities in gas hydrate-bearing seafloor mounds at Storfjordrenna, offshore Svalbard in the high Arctic, where distinct methane flux regimes ranging from steady-state dynamics, recent increase in subsurface diffusive flux, and gas seepage were identified. Populations of anaerobic methanotrophs and sulfate-reducing bacteria were highest at the seep site, while a recent increase in methane influx was associated with decreased community diversity. Despite high methane fluxes and methanotroph doubling times estimated at 5-9 months, microbial community responses were largely synchronous with the advancement of methane into shallower sediment horizons. Together, these provide a framework for interpreting subseafloor microbial responses to methane escape in a changing Arctic.

Geochemistry

Arctic seafloor

methane escape

gas hydrate mounds

Microbially-generated methane in marine sediments has been estimated at 1013 – 1014 g per year1. Microbial anaerobic methane oxidation (AOM) is responsible for consuming the majority of this methane–up to 90%1–before it can escape to the hydrosphere. This globally-widespread2 microbial methane filter consists of very slow-growing3,4, currently uncultured clades of anaerobic methanotrophic archaea (ANME) and often-symbiotic sulfate-reducing bacteria (SRB) that thrive at sulfate-methane transitions (SMTs), sediment depths where methane is oxidized with sulfate (SR-AOM)5. In contrast to the large areas where SMTs occur within the sediment, at discrete locations of active methane gas release, such as pockmarks and mud volcanoes, over 90% of the methane can escape the microbial filter, even though AOM rates at these sites are considerable (several nmol cm-3 d-1)6.

Methane release from the Arctic seafloor has received significant attention over the past two decades7. Seafloor permafrost methane venting to the atmosphere has been documented along a wide portion of the East Siberian Margin8, the South Kara Sea shelf9, and the upper slope of the Beaufort Sea10. Extensive geophysical surveys have characterized thousands of fault-associated seeps below warming waters along the West Spitsbergen (Svalbard) margin11,12; numerical modeling and U/Th dates from authigenic carbonates revealed that seepage has persisted here for hundreds to thousands of years13,14.

The Storfjordrenna trough mouth fan, ~50 km south of Svalbard, hosts gas hydrate-bearing mounds (GHMs) on the seafloor that are morphologically similar to those described in the Beaufort15 and Kara16 Seas (Fig. 1A). These GHMs lie below water depths of 370-390 meters, which approach the upper limit of gas hydrate stability in this area17. Gas leakage was observed at four of five GHMs, which are thought to have formed from hydrate accumulation and methane gas overpressure following glacial retreat17. Microbial community responses to subsurface methane release, whether driven by tectonic18, climate19 and/or oceanographic20 forcings, are important to constrain because they support macrofaunal communities21 of ecological and economic importance22. However, how these microbial communities respond to changes in methane release over time in Arctic cold seeps remains largely uncharacterized. Such knowledge is of immediate importance as environmental changes, from either natural or anthropogenic causes, could potentially result in a sudden increase of methane flux.

In the Arctic Ocean, abrupt release of methane from gas hydrate dissolution in the central Barents Sea has been hypothesized23, while methane release from the Deepwater Horizon oil spill into deep Gulf of Mexico waters was correlated with the growth of aerobic methane-oxidizing Gammaproteobacteria and oxygen drawdown24. Sediment microbial community responses to fluctuating methane regimes have been characterized at mud volcanoes25,26, and methane has recently been found to shape community structure at Storfjordrenna GHMs27. However, a dynamic understanding of how microbial activity may mitigate methane release in this habitat is currently unknown.

Changes in concentration gradients of porewater sulfate in marine sediments have been used to constrain the timing of subaqueous landslides28 and to indicate irrigation (through bioturbation or ascending gas bubbles29) or migration of upwards-diffusing methane14. Under steady-state conditions with a constant methane flux, sulfate decreases linearly with depth until the SMT is reached30 (Fig. 1B). At locations experiencing increases in methane flux, sulfate deviates from the linear profile, showing an abrupt decrease to <1 mM within a small depth range (Fig. 1C). The sulfate reduction (SR) zone thins in response to methane advancement into shallower layers, which stimulates AOM within them. Reactive transport modeling of the deviation from linearity in porewater sulfate profiles (Fig. 1D) can be used to estimate how long ago methane began to diffuse into shallower sediment zones, provided that other phenomena (advection, seawater irrigation, bioturbation, or mass transport deposits) are minimized or constrained14,28. Other geochemical signatures, combined with observations of free gas and gas hydrates, support a model where episodic methane emission occurs in pulses, with distinctive pre- and post-active stages31.

In this work, we constrain temporal responses of microbial communities as methane migrates upwardly towards shallow sediment horizons. Using samples and data from Storfjordrenna GHMs, where varying regimes of methane transport are evident, we employ interdisciplinary geochemical, numerical, and microbiological approaches. We report shifts in AOM rates, ANME/SRB abundances, and microbial communities concomitant with recent changes in methane flux, showing a tightly-coupled microbial response to intensifying subseafloor methane fluxes.

Identification of distinct regimes of methane transport.

To investigate microbial communities across regimes of methane dynamics in Storfjordrenna GHMs, we collected sediment cores from GHMs where seepage was either observed or absent (Fig. 1, Table S1). GHM3 and GHM4 show persistent hydroacoustic gas flares over multi-year surveys and thus represent active seepage¹⁷. Here, ascending free gas is the dominant carrier of methane, with limited movement in porewater³¹. In contrast, abundant trace elements in porewaters from GHM5 point towards upward fluid advection, which may have followed prior seepage activity³¹. Geochemical data, numerical modeling results, microbial communities, and functional gene abundance data are presented for these regimes in figures 2 through 4, grouped according to different methane transport scenarios. First, we present four cores from areas exhibiting steady-state methane-sulfate dynamics (i.e. no temporal changes in methane flux), with gravity core (GC) 1048 representing a low methane flux case from a non-GHM area (Fig. 2). Increasing methane fluxes in areas at GHM3 and GHM4 have allowed methane to diffuse into shallower sediment horizons, and sulfate profiles at these sites document the departure from a steady-state condition (Fig. 3). Finally, at an active seep site, downcore changes from a push core precisely taken by a remote-operated vehicle (ROV) from GHM3 capture biogeochemical signatures that reflect high methane flux, gas bubble emission, and advective delivery of seawater into the upper sediment column (Fig. 4). Corresponding data are provided in supplementary materials.

Field descriptions and general patterns.

Black-colored glaciomarine sediments were recovered in all cores, reflecting the precipitation of iron sulfide minerals resulting from high rates of sulfide production³². Authigenic carbonate nodules were retrieved in several cores, and chunks of gas hydrates several cm in diameter were observed between 40-50 cm below seafloor in a replicate of push core (PC) 1029. Cores PC1029 and GC1081 were taken from areas of gas seepage indicated by the white polygons in Fig. 1. Core recovery lengths ranging from 102 to 335 cm captured SMTs in all cores except for PC1029 (Table S1). Across all cores, alkalinity profiles that contrast with the decline in sulfate provide further support of AOM as the dominant sink for sulfate (Figs. 2A, 3A, 4A). In situ methane values are probably higher than those reported, as gas samples were taken from cores at atmospheric pressure. No bubbles or moussy sediment texture was observed in the recovered cores, discounting the possibility of degassing upon core retrieval.

Bacterial and archaeal 16S rRNA gene sequencing recovered 3.12 million sequences and 16,470 amplicon sequence variants (ASVs) after contaminants were removed. Bubble plots (Figs. 2C, 3C, 4B, left panels) show the fifteen most abundant taxonomic classes in the dataset, each of which individually constitute 1% or more of the total sequences, and combined account for 83.6% of reads in the dataset. The three most common ASVs, which alone comprise 22.2% of all sequences, belong to the class JS1 (phylum Atribacteria) which are thought to ferment organic matter³³. Two other dominant classes, Deltaproteobacteria and Methanomicrobia, are subdivided into genera of sulfate-reducing bacteria (SRB) and anaerobic methanotrophs (ANME), respectively (Figs. 2C, 3C, 4B, right panels). Respectively, ANME and SRB make up 10% and 12% of total sequences in this dataset. ANME are most dominant at or near SMTs, but correlate most strongly with peak AOM rates derived from numerical modeling. Two clades of sulfate-reducing bacteria that commonly associate with ANME at seeps, SEEP-SRB1 and SEEP-SRB2³⁴^,³⁵, share similar distribution patterns. Droplet digital PCR counts of the methane-fixing methyl-coenzyme reductase gene mcrA and dissimilatory sulfate reduction gene dsrAB span several orders of magnitude across cores and depths (Figs. 2D, 3D, 4C). Methods are detailed at the end of the manuscript.

Regimes of methane dynamics defined by geochemical and microbial data.

Steady-state biogeochemical equilibrium.

Three gravity cores from GHM5 (and one nearby, in the case of GC1048) showed linear decreases in sulfate, with methane present only below the SMT (Fig. 2A). We classify these sites as being at steady-state with respect to sulfate-methane dynamics, and hereafter refer to them as steady-state cores. Sulfide profiles mirror the shape of the alkalinity curves, peaking at SMT depths. In addition, macroscopic SMT-associated mucoid biofilms consisting predominantly of ANME-1³⁶, were observed in split cores at 63 and 68 cm in GC1070 and at 305 cm in GC1048 (Fig. 2A). Modeled AOM rates peak at or several cm above the SMT, and show maxima of < 20 nmol cm^-3 d^-1 (Fig. 2B); depth-integrated fluxes are below 1 mol m^-2 yr^-1 (Fig. S1). Though both ANME-1 subclades are nearly equally abundant in all samples from this study, ANME-1a are more dominant in steady-state cores (Fig. 2C). In GC1068, highest mcrA counts are seen at the peak AOM depth, though gene abundance profiles otherwise display considerable variability and dsrAB counts are typically low, between 10⁴–10⁵ copies per gram bulk sediment (Fig. 2D).

Increasing methane flux.

GC1045 was sampled from the southern margin of GHM3, and GC1081 from the center of GHM4 (Fig. 1). Sulfate profiles from these cores show concave-up curvature, indicating that the methane-sulfate dynamics are not at steady-state, but likely reflect a recent increase in methane flux¹⁴ (Fig. 3A). Porewater sulfate profiles show a rapid decrease in concentration down core and SMTs are well established. As estimated by our reduced modeling approach, total methane fluxes throughout these two cores have increased over the past two decades (Table S2). Modeling scenarios were constructed on a prior dataset of several porewater species from Storfjordrenna in an attempt to account for other processes, but only a scenario applying contrasts in methane flux adequately fit these data¹⁴. No fractures, mass transport deposits, porosity changes, or evidence of bioturbation were found in the gravity cores analyzed, and a buildup of ammonium to 60 µM (data not shown) in the first 50 cm of GC1045 allows us to discount the possibility of oxic bottom water intrusion. These cores experiencing increases in methane flux were initially assumed to be at steady-state before methane migrated upwards through the sediment column. Fluxes are integrated from all modeled AOM rates assuming AOM as the only sink for sulfate. Ammonium concentrations from GC1045 do not exceed 200 µM, and assuming a C:N ratio of 6.1, a maximum of only 0.61 mM of sulfate would be consumed by organoclastic sulfate reduction. Following these constraints, our model estimates peak AOM rates of ~200 nmol cm^-3 d^-1, which are an order of magnitude higher than those derived for steady-state cores (Fig. 3B). At the time of sampling, peak AOM rates were positioned at or up to 20 cm above the SMT depth (Fig. 3B), though running the model backwards or forward in time reveals an upward migration at a linear rate of 10 cm per year (Fig. 3B, Table S2). A more thorough description of this modeling approach is provided in the Supplemental Information.

In contrast to the steady-state cores, the dominance of Deltaproteobacteria, JS1, and Methanomicrobia are more apparent in GC1045 and GC1081, and ANME-1b are the most abundant ANME genus (Fig. 3C). Counts of mcrA are higher than in steady-state samples, reaching maxima around 10⁷ copies per gram at SMTs in both cores (Fig. 3D). Though dsrAB counts are also comparable at these SMTs, higher dsrAB abundances at shallower depths in GC1045 likely reflect a larger or more diverse sulfate-reducing community than in GC1081.

Active methane seepage.

PC1029 was recovered from an established patch of frenulate siboglinid tubeworms (Oligobrachia sp. CPL clade, Fig. S2) whose chemosynthetic lifestypes are supported by sulfide generated from SR-AOM at sites with high methane discharge^21,37,38. Observations of vigorous gas bubbling and recovery of gas hydrate support the inference that the site was experiencing high methane seepage at the time of sampling. Sulfate concentrations at near-seawater values up to 10 cm below seafloor at PC1029 (Fig. 4A) may be attributed to seawater infiltration (siboglinid bioirrigation or bubble-driven convection) and/or sulfide oxidation from bacterial symbionts³⁹. Further downcore, the incomplete drawdown of sulfate and high methane concentrations suggest that sulfate-coupled AOM is an ongoing process, pointing towards a high methane flux at the center of GHM3. As processes other than sulfate diffusion from seawater are not accounted for in our model parameterization, we are unable to precisely calculate AOM rates from PC1029. Our rough estimation of the AOM rate based on the part of the sulfate profile with the greatest concentration gradient (10-15 cmbsf) yields a peak AOM rate on the order of 10³ nmol cm^-3 d^-1. This rate estimate would be increased significantly by accounting for siboglinid-driven pumping of bottom seawater sulfate, or sulfide reoxidation mediated by their endosymbionts. Nevertheless, this estimated AOM rate is an order of magnitude higher than the rates calculated for cores experiencing increases in methane flux shown in Figure 3B.

Microbial communities from PC1029 show higher percent abundances of several classes, notably Bacteroidia and Gammaproteobacteria, than in cores representing other regimes of methane dynamics (Fig. 4B). ANME-2 are present at 1-3 cm in PC1029, but ANME-1b are dominant at depths with lower sulfate concentrations (Fig. 4B). ANME-1a are nearly absent, agreeing with recent observations from this seep location²⁷ and contrasting the two other regimes. Near-equal abundances of SEEP-SRB1 and SEEP-SRB2 at PC1029 are reminiscent of GC1081, the other core from a seep site (Figs. 4B, 3C). Highest mcrA concentrations, exceeding 10⁸ copies per gram bulk sediment, were recovered in PC1029, even in depths with high sulfate, low methane, and low alkalinity (Fig. 4C); these values are comparable to ANME cell counts reported from other seep sites²⁵^,⁴⁰. Counts of dsrAB were over an order of magnitude lower than mcrA throughout the core, but still higher than in nearly all other samples from different regimes.

Response times of ANME and SRB to methane pulses as inferred from porewater modeling.

In areas experiencing increases in methane flux, modeled AOM peaks migrate upward with time (Fig. 3B). Different depths, and thus microbial communities inhabiting them, can be assigned by the time they experienced (or are expected to experience) this upward-moving peak in AOM. At GC1045 and GC1081, AOM fronts migrate upward at roughly 10 cm per year (Fig. 3B, Table S2); thus GC1045 communities from 66, 76, 86, and 110 cm depths respectively correspond to AOM peaks at the time of collection, and one, two, and over four years before. Highest relative abundances of SEEP-SRB1 and total ANME, as well as highest concentrations of mcrA, are seen in the community sampled at 76 cm (Fig. 3C, D), suggesting these taxa dominate microbial communities after about a year following methane migration into this sediment horizon. In contrast, ANME-1b are highest in the community from 66 cm, which may reflect a quicker growth or a preference for lower methane concentrations compared to ANME-1a. In the timespan from one to four years after the AOM pulse has passed through, mcrA abundances decreased by nearly three orders of magnitude, but dsrAB by less than one. After four years, the AOM pulse moves onward and microbial communities are starved of sulfate, ANME and SRB populations respectively decrease from 46 to 1.1% and 22 to 1.8% of the total community.

GC1081 communities from 56.5, 69, and 86 cm correspond to maximum AOM rates from the time of sampling, one and a half, and three years ago, respectively, while the community at 49 cm is associated with high (but not yet peaking) AOM rates (Fig. 3B). In contrast to communities from GC1045, ANME percent abundances do not decrease as quickly, and SEEP-SRB1 increases with depth (Fig. 3C). A similar trend of ANME-1b growth preceding ANME-1a is noticed, but surprisingly ANME-1b are present in high relative abundance at 24 cm, where AOM rates are not expected to be significant until two years after sampling. Concentrations of mcrA and dsrAB both roughly correspond to the present-day AOM pulse, showing no lag time with respect to methane influx (Fig. 3D). Across samples for which we were confidently able to calculate AOM rates, mcrA gene abundances correlate positively with rates when plotted on a log-log scale (Fig. S3).

Microbial community diversity and analysis.

The three most abundant classes in our dataset, the Methanomicrobia, Deltaproteobacteria, and JS1, a class of Atribacteria, are especially dominant in communities from cores experiencing recent methane influx (Fig. 3C). Besides these major groups, other poorly-understood taxa include the Aminicenantes, Anaerolineae, and Phycisphaerae, all thought to be fermentative saccharolytic heterotrophs⁴¹^,⁴²^,⁴³. Dehalococcoidia, also present, contain members capable of reductive dehalogenation⁴⁴. Among all ASVs in our dataset, we identified ninety-nine whose relative abundances were significantly different across communities when grouped according to methane dynamics (Fig. 5). These ASVs on average comprise 12.7% of the sequences in communities associated with active seepage, 1.6% of communities experiencing methane flux, and 9.7% of steady-state communities. When compared to the other two regimes, communities from sites displaying steady-state sulfate-methane dynamics contained higher abundances of several ASVs belonging to Aminicenantia, Dehalococcoidia, and Woesearchaeota (Fig. 5). In addition, one ANME-1a ASV was higher in this group, though only three of the 41 ANME ASVs in the entire dataset were differentially abundant across methane regimes. Several ASVs belonging to SEEP-SRB1 and Desulfatiglans also displayed variation among regimes (Fig. 5).

In communities from cores experiencing high methane flux, Shannon-Weiner alpha diversity indices decrease as depths approach peak model-derived AOM rates (Fig. 6A). Linear regressions show no such decrease in diversity across AOM peaks from steady-state areas (Fig. 6B). Interestingly, in samples from core GC1069, highest diversity is seen at depths of peak AOM, while the opposite is apparent in GC1070 (Fig. 6B).

Differences in community structure are evident across regimes of methane dynamics, with communities from PC1029 showing especially clear separation from those in steady-state cores (Fig. 7A). These distinctions were still observed even when seep samples were omitted (R² = 0.06, p < 0.001), and when samples from above or below the SMT were considered separately (R² = 0.19, p = 0.003; R²= 0.22, p = 0.001 respectively). In addition, we classify samples according to three geochemical zones they inhabit based on the shapes of porewater sulfate profiles: the linear sulfate reduction (SR) zone, the nonlinear SR zone impacted by recent methane influx, and below-SMT. Community structure also varied significantly across these redox zones, though PERMANOVA tests revealed only 9.6% of the variance among communities could be explained by redox zone (Fig. 7B) as compared to 19.7% by methane regime. Though containing high relative abundances of ANME, communities in nonlinear SR regions of cores experiencing increasing methane flux were more similar to below-SMT communities than those in linear SR zones, suggestive of recent adaptations to methane influx (Fig. S4). Aside from methane regimes and redox zones, communities also varied significantly by the GHM and core they were sampled from (32.5% and 22.7% of variance), suggesting these Arctic GHM communities contain a high degree of biogeographic heterogeneity²⁷ that remains unconstrained.

The presence of distinct regimes of methane flux at Storfjordrenna GHMs allows us to examine concomitant changes in inferred AOM activity and microbial community composition. We conceptually summarize results from integrated geochemical, numerical, and microbiological analyses that characterize three distinct biogeochemical regimes corresponding to changes in methane supply across all seven cores (Fig. 8).

In panel A, steady-state sulfate-methane dynamics occur when methane consumption is balanced by sulfate diffusion from seawater. ANME and SRB abundances often do not peak around the SMT, and these populations are accompanied by several other microbial groups, many of which may be slow-growing anaerobic fermenters of organic matter. Growth of biofilms at SMTs may also be supported over long timescales, if given steady supplies of methane and sulfate. In panel B, AOM draws down sulfate and a SMT is established, though it has not reached steady-state dynamics. At GC1045, diffusion modeling of sulfate porewater profiles indicate the methane front began around 290 years ago. As methane flux increases, this pulse diffuses upwards, stimulating AOM in increasingly shallow sediment layers. AOM rates are approximately an order of magnitude higher than in panel A, supporting the growth of ANME/SRB and decreasing microbial community diversity. In panel C, gas seepage and the presence of hydrates at PC1029 indicate methane is at or above saturation in porewaters throughout the core. Sulfate is delivered into the sediment column through seawater infiltration driven by bubble-driven convection, bioirrigation by frenulate siboglinids, or reoxidation of sulfide by their endosymbionts. High sulfate concentrations predict high (though unconstrained) AOM rates, supporting large populations of ANME and SRB based on respective counts of mcrA and dsrAB. Though direct cell counts were not obtained, the presence of several other abundant bacterial and archaeal classes suggests these shallow sediments support high cell densities overall.

Model-derived methane fluxes from Storfjordrenna cores GC1045 and 1081 are an order of magnitude higher than those from seeps associated with pockmark features at Vestnessa Ridge, west of Svalbard⁴⁵. When compared to other estimates across continental margins worldwide, the values we report for these two cores are high but well within reported ranges, while fluxes from steady-state cores fall within the middle⁴⁶. At the seep site, PC1029, our methane flux estimate on the order of 10² mol m^-1 yr^-1is several times the maximum of other modeled AOM rates at seep sites⁴⁶, but less than the highest empirically measured AOM rate⁶. We acknowledge that the rate estimated from PC1029 is associated with large uncertainties, as we were not able to satisfactorily fit the modeled curve to empirically-derived sulfate data under the current setup of the model. The model currently does not consider gaseous phase transport or bioturbation, which would enhance gaseous methane transport from deeper sediments, nor does it include sulfate infiltration from the bottom water or sulfide oxidation, which may provide additional substrates for SR-AOM. Though the timing of seepage at the center of GHM3 is unconstrained, the large populations of anaerobic methanotrophs and sulfate reducers supported by high methane fluxes may indicate stable conditions over timescales of years⁴⁷.

We now consider ANME doubling times at sites experiencing an increase in methane flux. If we interpret the downcore increases in mcrA concentrations approaching the SMT (Fig. 3D) as methane-fueled ANME growth, doubling times of 147 days (GC1081, 23 to 56.5 cm) to 261 days (GC1045, 66 to 76 cm) can be derived by assuming one copy of mcrA per ANME genome⁴⁸. These values complement the only other published estimate of in situ ANME doubling time, at approximately 100-200 days⁴⁷. We can then estimate per-cell AOM rates across SMTs at an average of 0.65 pmol CH₄ d^-1, comparable to the 0.5–1.8 reported in a bioreactor where AOM was stimulated⁴⁹. The observation that AOM rates and mcrA abundances in areas experiencing increasing methane flux peak at nearly the same depths suggests that the notoriously slow ANME doubling times may not present a significant lag (only 0-1 years) in the response of the benthic biogeochemical methane filter, in contrast with findings from a prior modeling study from coastal ocean sediments with methane produced in-situ through acetotrophic methanogenesis⁵⁰.

At Storfjordrenna GHMs, ANME-1 is the most abundant anaerobic methanotroph in nearly all communities (Figs. 2C, 3C, 4B). However, our observations of ANME-2 at sulfate-rich surface sediments in PC1029 (Fig. 4B) agree with previous findings⁵. At all other depths and locations, the reasons for ANME-1 dominance at Storfjordrenna GHMs is unclear. Genomic explanations may include the lack of an energetically expensive nifDHK nitrogenase in ANME-1⁵¹, and fewer multi-heme cytochromes thought to be involved in direct intercellular electron transfer⁵². Additionally, the lack of mer in ANME-1⁵³, and the suggestion of met⁵⁴ or fae, hps, and adh⁵⁵ as functional replacements in the reverse methanogenesis pathway suggest that separate niches could exist for these two clades, but this remains an open question. In cores experiencing increased methane flux, ANME-1a and ANME-1b were present at near-equal abundances, with ANME-1b sequences more abundant at shallower depths. In contrast, the 1a subclade was more dominant in three of four steady-state cores (Figs. 2C, 3C). ANME percent abundances and mcrA concentrations are mostly higher in cores experiencing increasing methane flux (Figs. 2, 3, 4). This may point towards a boom-and-bust cycle where methane influx into shallower sediment layers quickly stimulates a large but ultimately unsustainable methanotrophic population, which may decline as sulfate is drawn down or as other community members establish.

In several instances, high abundances of ANME or concentrations of mcrA are seen at depths above those where methane is expected (in GC1081 and GC1048) or below the SMT (GC1069). These may reflect the limited resolution of our alkalinity and sulfate measurements. At GC1068 and GC1069, any cessation in the methane supply would allow sulfate to diffuse into a deeper depth without affecting the linearity of the sulfate profile. Alternatively, this may indicate inactive relic communities, though ANME-1 may still be capable of AOM⁵⁶ or even methanogenesis⁵⁷ when starved of sulfate.

Co-occurrences between ANME-1b and SEEP-SRB2 have been reported², and their relative abundances appear to mirror each other in GHM4 samples (Fig. 3C). Both clades of SEEP-SRB, as well as Desulfatiglans–the most common SRB genus unaffiliated with ANME in our dataset–are presumed to oxidize a wide variety of reduced hydrocarbons³⁵. The presence of several potentially fermentative and saccharolytic clades like the Atribacteria, Aminicenantes, Anaerolineae, and Phycisphaerae may reflect alternate organic matter-dependent metabolic strategies that are interrupted by ANME and SRB when methane enters sulfate-rich porewaters. Macroscopic ANME-dominated biofilms found at two SMTs in GC1048 and GC1070³⁶ contained mcrA in concentrations of up to 7.6 x 10¹⁰ g^-1. We hypothesize these biofilms reflect sediment regimes that have experienced steady methane and sulfate supply over many years. ANME biofilms have been described at SMTs in other subseafloor locations, often in fracture-dominated cores⁵⁸.

Microbial communities inhabiting Storfjordrenna GHMs show lower richness and evenness than most other reported communities from methane seeps, sulfate-methane transition zones, and marine subsurface environments². Broadly, diversity decreases with depth, but only significantly across depths corresponding to peak AOM rates in high methane flux areas (Fig. 6). This decrease is still evident when other metrics of richness (number of ASVs) and evenness (inverse Simpson) are used (data not shown). In communities recently impacted by methane flux (i.e. inhabiting non-steady-state SR zones of cores where methane flux is increasing), this decrease in diversity and convergence towards a community type found below SMTs may be associated with certain taxa being outcompeted by ANME/SRB on timescales of years as the methane front progresses. Cell generation times can decrease by several orders of magnitude across the SMT, below which community assembly can be influenced more by the slow growth of few taxa, such as Atribacteria, capable of thriving in an energy-limited environment than by evolutionary adaptation during burial⁵⁹.

Below-SMT communities are dominated by Atribacteria of the JS1 class (Figs. 3C, 4B), while similar observations have been reported in methane-rich deep Antarctic marine sediments⁶⁰ and in a submarine mud volcano offshore Japan⁶¹. Four JS1 ASVs were identified across different regimes of methane flux and positions above or below the SMT (Fig. 5), though interestingly, one of them (ASV91) was preferentially abundant in above-SMT steady-state and below-SMT increasing-flux communities, evidence towards its persistence during methane migration into shallower sediment horizons. Despite steady-state communities showing larger numbers of differentially abundant ASVs, two Calditrichia (genus Caldithrix) were higher in communities experiencing increased methane flux, while four Campylobacteria from the genera Sulfurimonas and Sulfurovum and seven Gammaproteobacteria were associated with active methane seepage (Fig. 5). Sulfurovum is capable of oxidizing elemental sulfur or thiosulfate using oxygen or nitrate as electron acceptors³².

The presence of differentially abundant ASVs at distinct regimes may reflect the sampling of comparatively shallow sediments at PC1029, and an influence from macrofauna. Sulfur-oxidizing gammaproteobacterial symbionts of the siboglinid frenulate Oligobrachia have previously been reported in cold seeps⁶². Notably, there is an absence of Oligobrachia and a decreased prevalence of seafloor bacterial mats at GHM5, where steady-state cores were collected²¹. Despite the short (several km) distances between individual GHMs, many interdependent factors, such as physical disturbances, differences in fluid flow regimes, and colonization of foundation species provide heterogeneity across seep ecosystems⁶³.

Summary and prospective

In summary, our integrated approach allows us to detail regimes of methane transport where 1) steady-state sulfate-methane dynamics supports moderate rates of AOM at SMTs, low ANME/SRB populations, and a diverse community of organic matter degraders; 2) as methane flux increases, diffusion of methane into shallower sediment horizons stimulates ANME growth therein with little lag time; 3) seepage and sulfate transport into shallow sediments support high populations of ANME and SRB. Cold seeps are dynamic systems that undergo temporal perturbations in methane flux. These results highlight the importance of framing microbial community data and estimates of their metabolic processes within a spatially and temporally constrained geochemical context to more thoroughly understand microbial contributions in structuring habitats and mediating biogeochemical cycles.

The incorporation of genomic data into reactive transport models describing other microbially-mediated processes has demonstrated utility in predicting subsurface microbial responses⁶⁴. A modeling scenario that considers the dynamics of ANME growth may be of use in constraining estimates of marine subsurface methane flux into the hydrosphere. Global microbial methane filter efficiencies of 50-60%⁵ have been used in modeling studies¹⁹, but seep sites display wide heterogeneity⁶. Our finding that mcrA gene copy numbers correlate positively with modeled AOM rates provides some justification for coupling these populations and their associated activities (Fig. S3), mirroring the coupling of methane fluxes and transcripts of methane cyclers in peat soils⁶⁵. Though microbial community data can provide explanatory power for predicting ecosystem processes, community changes do not always coincide with processes being measured⁶⁶. At higher resolutions, -omics strategies capable of characterizing functional genes, transcripts, pathways, and draft genomes can link sequence data with processes and characterize ecosystem changes⁶⁷, or even apply these data into biogeochemical models to infer the presence of cryptic cycles⁶⁸. Further studies could apply the framework discussed here towards interpreting the biogeochemistry of seep ecosystems at other locations, or to other microbially-mediated cycles constrained by distinct mechanisms of solute transport.

Fieldwork and Sample Collection

Samples and data were collected aboard the RV Helmer Hanssen on CAGE cruise 16-5, from June 16^th to July 4^th, 2016. Bathymetric data were acquired with the RV Helmer Hanssen’s shipboard Kongsberg Simrad EM 302 multibeam echo sounder using a max. frequency of 34 KHz and a max. swath as a function of depth of 5.5. Gas flares were detected with single (split)-beam EK60 and multibeam EM302 echosounders using 18 and 38 KHz transducers.

Gravity core (GC) 1045 was recovered from the south slope of GHM3, and GCs 1068–1070 from three locations at GHM5. GC1081 was collected from gas seepage area at GHM4. Once recovered, the plastic liner containing the core was removed from the barrel, sectioned into 1 m segments, labeled, and split in half with a table saw to obtain working and archive halves. Core halves were stored horizontally at 4˚C. Following sectioning, Rhizons were used to sample porewater on archive halves. Sediment headspace gas samples for methane measurements were collected from depressurized cores, and thus should be considered underestimates for in situ concentrations. 5 ml bulk sediment was collected with cutoff plastic syringes from the working half of the core, transferred to 20 ml headspace glass vials with 5 ml 1M NaOH and 2 glass beads, capped with rubber septa and aluminum crimpers, and stored at 2˚C. Total alkalinity (TA) was titrated onboard less than a few hours after the syringes were detached from the Rhizons. Depending on the expected TA, we used 0.1 to 0.5 ml of porewater for titration in an open beaker with constant stirring. pH was manually recorded with every addition of 0.0012M HCl. 7-10 measurements were performed for every sample. TA was calculated from the recorded pH and amount of acid added using the Gran function. Details of the calculation were reported previously⁶⁹. Increases in porewater alkalinity determined by onboard titrations were used to roughly constrain the SMT depths (within 30 cm) for sampling purposes.

Sediment microbiology samples of 2 cm depth were then taken every 5-10 cm near the SMT and every 20-50 cm above and below it. Less than 12 hours after cores were collected, ethanol-sanitized spatulas were used to scrape away the outer several mm of sediment from the working core half, and ~100 g from the interior of each sample was placed into a sterile whirlpak bag (VWR) and immediately frozen at -80˚C.

Replicate PVC push cores for geochemical and microbiological sampling were collected ~30 cm from the seep at GHM3 using a Sperre Subfighter 30k remotely operated vehicle (ROV) equipped with a raptor arm from the Centre for Autonomous Marine Operations and Systems (AMOS). Recovery ranged from 23 to 50 cm. Rhizons were used to extract porewater from one core, and microbiology samples were extruded on deck from the other in 2-cm sections using an ethanol-sanitized spatula. These were placed into sterile bags and frozen immediately at -80˚C. Deep-frozen sediment samples were shipped from UiT-Tromsø to Oregon State University (OSU) in a Cryo Solutions MVE Doble 47 dry shipper and were subsequently stored at -80˚C.

Geochemistry

Sulfate content in porewater was analyzed by a Dionex ICS1100 ion chromatography (IC) at the Geological Survey of Norway (NGU). An IonPac AS23 column was equipped on the IC with the eluent (4.5 mM NaCO₃ and 0.8 mM NaHCO₃) flow set to be 1 mL/min. Due to a dilution issue when analyzing sulfate concentrations with IC, measured values were corrected by assuming a constant chloride concentration of 556 mM across the samples. From our previous knowledge of chloride concentration in the region, the concentration can be at most 4% apart from the concentration we assigned for correction (this translates to a few mM uncertainty in the sulfate concentration, though error bars are barely visible beyond widths of points). However, this correction does not affect our interpretation of methane dynamics based on sulfate profiles, as we observed concomitant increases in alkalinity from TA measurements.

Total sulfide (ΣHS) concentrations were measured spectrophotometrically following the Cline method⁷⁰. Samples were preserved onboard with 23.8 mM Zn(OAc)₂ solution onboard < 30 minutes after the syringes were disconnected from the Rhizons. The samples were then kept frozen until shore-based analysis. Details of the analyses were also given in Latour et al⁶⁹. Depending on the factor of dilution, the detection limit is around tens of μM. To determine the concentration of dissolved methane in the porewater of the sediment matrix, a conventional headspace method was applied¹⁷.

Gas measurements were performed using a Thermo Scientific Trace 1310 gas chromatograph equipped with a flame ionization detector (GC-FID) and a Thermo Scientific TG-BOND alumina (Na₂SO₄) column (30 m × 0.53 mm × 10 μm).

Modeling

For the sites with non-steady-state porewater profiles (GC1045 and GC1081), we applied the same reduced modeling scheme as described previously¹⁴. Only sulfate and methane were considered in this reduced model, with AOM as the only reaction consuming both constituents. We assigned the diffusion coefficients of sulfate and methane to be 1.58x10^-2 and 3.01x10^-2 m² yr^-1, respectively, when corrected for tortuosity and temperature. We simulated a 60-meter sediment column, which is the bottom of the gas hydrate stability zone in the area. AOM rates were controlled by the theoretical maximum rate, methane and sulfate concentrations, and the lower boundary condition of methane. The theoretical maximum AOM rate was determined to be 2 mol m^-3 yr^-1, based on fitting porewater profiles around the depth of the SMT for all cores¹⁴. Different lower boundary conditions of methane were assigned, which were constrained by the curvatures of the sulfate profiles. A higher concentration of methane for the lower boundary condition results in a more abrupt change in the sulfate concentration gradient, and vice versa. We assigned seawater sulfate and methane concentrations for the initial and upper boundary conditions; a no-flux lower boundary condition was used for sulfate.

For the sites with steady-state porewater profiles, we applied the comprehensive modeling scheme as described previously¹⁴. We coupled the CrunchFlow routine for this simulation with profiles of sulfate, ΣHS, TA, calcium, magnesium, and ammonium to constrain the reaction network. Because ammonium concentration measurements were not available for these sites, we assumed a similar profile, and therefore the same organic matter degradation rate, as in core GC920 from a previous study¹⁴. Seawater composition was used for upper and initial conditions for all solutes. Except for methane, no-flux lower boundary conditions were assigned for all solutes. Methane was produced at the deepest cell in the model from an imaginary mineral. Such a setup is to simulate a methane source that was not produced in situ in the sediments, and to overcome the limitation of the software package used that required all solutes to have the same boundary condition. A detailed description of the reaction network can be found in Hong et al. (2017)¹⁴. Depth-integrated fluxes were calculated by summing the product of modeled AOM rates and the cell thickness, considering all modeled AOM rates from above the methanogenic zone. GC1048 is a special case compared with the other steady-state sites: Its sulfate profile is also slightly bent as in the profiles of GC1045 and GC1081. To satisfactorily fit this observation, we similarly applied a slightly intensifying methane flux when we modeled GC1048. However, the increase in methane flux we assigned for GC1048 is much smaller than in the two sites with obvious non-steady-state sulfate profiles. To eliminate concerns that only three porewater samples were collected above the SMTs which may lead to misinterpretation of the status of pore fluid system (i.e., steady-state vs. increasing methane flux), we also fit the porewater profiles from the four steady-state sites with model results for scenarios that assume an increase in methane flux (Fig. S5). Despite the scarcity of this sulfate data collected, we conclude it is difficult to misinterpret these as steady-state sites. For both modeling approaches, more detail including assumptions, equations, and setup is provided in supplementary information.

DNA Extraction, Amplification, Sequencing, and Analysis

DNA was extracted from sediments in a clean laminar flow hood using a Qiagen DNeasy PowerSoil kit following the manufacturer’s protocol. The Earth Microbiome Project 16S Illumina Protocol was used to prepare amplicons for sequencing. Briefly, V4 regions of bacterial and archaeal 16s rRNA genes were amplified in triplicate 25 ul reactions using universal 515-forward and 806-reverse primers⁷¹ modified with dual-indexed Illumina sequencing adapters⁷². The thermal cycling protocol of Caporaso et al 2011⁷¹ was followed without modifications. After confirming amplification with agarose gel electrophoresis, triplicate PCR products were pooled and purified with a Qiagen QIAquick PCR purification kit. Amplicon concentrations were quantified with a Qubit fluorometer using the Qubit dsDNA high sensitivity assay kit and pooled in equimolar amounts. Illumina Miseq V2 paired-end 250 bp sequencing was performed by technicians at Oregon State University’s Center for Genome Research and Biocomputing (CGRB). Two sediment-free DNA extraction blanks were amplified and included in the sequencing run.

Working in R version 3.6.1, 16S rRNA amplicon data was processed with DADA2⁷³ (version 1.12.1) following an established pipeline⁷⁴. Reads were denoised, chimeras removed, and taxonomies classified using version 132 of the SILVA nonredundant 16S reference database⁷⁵. Sequences were aligned with DECIPHER⁷⁶ (version 2.12), and a phylogenetic tree was constructed using phangorn⁷⁷ (version 2.5.5). Phyloseq⁷⁸ (version 1.28.0) was used to combine read count data with sample and taxonomy information. Sequences identified as Eukaryotes, Chloroplasts, or Mitochondria were removed, and the “combined” method of decontam⁷⁹ (version 1.4.0) was then used to identify and remove 81 contaminant ASVs. In addition, after noticing the presence of Micrococcus in one blank sample, all four ASVs from this genus were manually removed. In total, the removed ASVs comprised 1.05% of the reads in the dataset.

Blanks and other samples with less than 8,931 reads were removed, and alpha diversity metrics (ASV richness, Chao1, Shannon, and Simpson indices) were then determined. Using vegan⁸⁰ (version 2.5-6), weighted Unifrac⁸¹ distances calculated from a Hellinger-transformed ASV count table, and PERMANOVA tests were run to assess differences in community structure among groups. DESeq2⁸² (version 1.24.0) was used to identify differentially abundant ASVs among three discrete regimes of methane dynamics. Each above-SMT methane regime was compared against the other two combined. Below-SMT samples only included two regimes, because all active seepage samples from PC1029 had sulfate concentrations above 1 mM. In this core, where AOM rates could only be roughly estimated, we used a peak AOM rate depth of 13 cm, which corresponded to the steepest decline in porewater sulfate.

Droplet Digital PCR

Droplet digital PCR (ddPCR) was used to quantify abundances of functional genes dsrAB and mcrA using primer pairs described by Kondo⁸³ and Luton⁸⁴, respectively. Reactions of 22 ul volume were prepared in a clean PCR hood in 96-well plates using 1x Bio-Rad QX200 ddPCR EvaGreen Supermix, 200 nM primers, and and 0.88 ul of tenfold-diluted genomic DNA. Droplets were generated on a QX200 AutoDG Droplet Generator using automated droplet generation oil for EvaGreen Supermix (Bio-Rad). Thermal cycling was performed immediately afterwards on a Veriti 96-well thermal cycler. Protocols began with a single initialization step at 95˚C for 5 minutes and then proceeded to 40 cycles of denaturation at 95˚C for 30 seconds, annealing for 1 minute (at a temperature of 53 for mcrA and 58 for dsrAB), and for mcrA only, an extension at 72˚C for 75 seconds. Signal stabilization steps (4˚C for 5 minutes, then 90˚C for 5 minutes) were then performed before maintaining a 4˚C hold. To ensure uniform heating of all droplets, the ramp rate for all amplification cycles was set to 2˚C/minute. Reactions were kept at 4˚C overnight and read with the Bio-Rad QX200 Droplet Reader the following morning. Droplet generation and reading were performed by the lead author at OSU’s CGRB core facility. Normalization was performed by inspecting fluorescence distributions using Quantasoft software (Bio-Rad). Threshold fluorescence values were manually imposed by visually inspecting distributions of DNA extraction blank and no-template-added control samples. Amplicon copy numbers per well were then converted to copies per gram wet sediment.

Data Availability

Raw sequence data and associated metadata have been deposited to the NCBI Sequence Read Archive (BioProject #PRJNA533183). Other data used in this study (geochemical data, modeled rates, and phyloseq objects) have been made publicly available at https://github.com/sklasek/svalflux.

Acknowledgements

We thank the officers and crew of R/V Helmer Hanssen on the CAGE cruise 16-5, cruise leader Michael Carroll, chief engineers Bjørn Runar Olsen, Pedro De La Torre, Frode Volden, and researcher Stein Nornæs from the Centre for Autonomous Marine Operations and Systems (AMOS) for helping with sample acquisition and ROV operation. We also thank Stefan Bünz for providing bathymetric data acquired by the Norwegian Centre of Excellence, Centre for Arctic gas hydrates, Environment, and Climate (CAGE). Mark Dasenko, Anne-Marie Girard Pohjanpelto, and Jessica Nixon at the Oregon State University Center for Genome Resources and Biocomputing (CGRB) provided support with DNA sequencing and droplet digital PCR. Arunima Sen gave helpful interpretations of tubeworm observations and behavior. This work is supported by the Research Council of Norway (RCN) through its Centres of Excellence funding scheme project no. 223259, and also by the U.S. Department of Energy (DE-FE0013531). S.R. was supported by the US National Science Foundation Research Experience for Undergraduate (REU) program. M.E.T. acknowledges a fellowship from the Hanse-Wissenschaftskolleg (HWK), Germany. W.H. acknowledges the support from the RCN-funded project NORCRUST (grant no. 255150).

Author Contributions

SAK and WH designed the study, participated in fieldwork, generated data, analyzed results, and wrote the manuscript. FG participated in fieldwork and contributed methane measurements. AP participated in fieldwork and provided bathymetric mapping. SR and KH assisted with lab work. MET and FSC analyzed and interpreted results and edited the manuscript.

Additional Information

Competing financial interests: The authors declare no competing financial interests.

There is NO Competing Interest.

SupplementalsKlaseksvalflux.pdf
Supplemental information for Klasek et al.
SupplementalsKlaseksvalflux.pdf
Supplemental information for Klasek et al.
TableS2.xlsx
Supplementary Table S2
TableS2.xlsx
Supplementary Table S2

Distinct methane-driven biogeochemical regimes in Arctic seafloor gas hydrate mounds

Status:

Version 1

Abstract

Figures

Introduction

Results

Identification of distinct regimes of methane transport.

Field descriptions and general patterns.

Regimes of methane dynamics defined by geochemical and microbial data.

Response times of ANME and SRB to methane pulses as inferred from porewater modeling.

Microbial community diversity and analysis.

Discussion

Methods

Fieldwork and Sample Collection

Geochemistry

Modeling

DNA Extraction, Amplification, Sequencing, and Analysis

Droplet Digital PCR

Declarations

Data Availability

Acknowledgements

Author Contributions

Additional Information

References

Additional Declarations

Supplementary Files

Status:

Version 1