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 seepage17. Here, ascending free gas is the dominant carrier of methane, with limited movement in porewater31. In contrast, abundant trace elements in porewaters from GHM5 point towards upward fluid advection, which may have followed prior seepage activity31. 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 production32. 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 matter33. 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-SRB234,35, 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-136, 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 104 –105 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 flux14 (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 data14. 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 107 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 discharge21,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 symbionts39. 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 103 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 location27 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 108 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 sites25,40. 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 heterotrophs41,42,43. Dehalococcoidia, also present, contain members capable of reductive dehalogenation44. 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 (R2 = 0.06, p < 0.001), and when samples from above or below the SMT were considered separately (R2 = 0.19, p = 0.003; R2 = 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 heterogeneity27 that remains unconstrained.