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

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 ux remains largely uncharacterized. We investigate microbial communities in gas hydrate-bearing seaoor mounds at Storfjordrenna, offshore Svalbard in the high Arctic, where distinct methane ux regimes ranging from steady-state dynamics, recent increase in subsurface diffusive ux, and gas seepage were identied. Populations of anaerobic methanotrophs and sulfate-reducing bacteria were highest at the seep site, while a recent increase in methane inux was associated with decreased community diversity. Despite high methane uxes 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 subseaoor microbial responses to methane escape in a changing Arctic.


Introduction
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 lter 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 lter, even though AOM rates at these sites are considerable (several nmol cm-3 d-1)6. Methane release from the Arctic sea oor has received signi cant attention over the past two decades7. Sea oor 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 sea oor 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 ve 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 ux. 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 uctuating 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 ux, sulfate decreases linearly with depth until the SMT is reached30 (Fig. 1B). At locations experiencing increases in methane ux, sulfate deviates from the linear pro le, 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 pro les ( 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 ux, showing a tightly-coupled microbial response to intensifying subsea oor methane uxes.

Results
Identi cation 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 ares over multi-year surveys and thus represent active seepage 17 . Here, ascending free gas is the dominant carrier of methane, with limited movement in porewater 31 . In contrast, abundant trace elements in porewaters from GHM5 point towards upward uid advection, which may have followed prior seepage activity 31 . Geochemical data, numerical modeling results, microbial communities, and functional gene abundance data are presented for these regimes in gures 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 ux), with gravity core (GC) 1048 representing a low methane ux case from a non-GHM area (Fig. 2). Increasing methane uxes in areas at GHM3 and GHM4 have allowed methane to diffuse into shallower sediment horizons, and sulfate pro les 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 re ect high methane ux, 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, re ecting the precipitation of iron sul de minerals resulting from high rates of sul de production 32 . Authigenic carbonate nodules were retrieved in several cores, and chunks of gas hydrates several cm in diameter were observed between 40-50 cm below sea oor 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 pro les 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 fteen 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 33 . 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 34,35 , share similar distribution patterns. Droplet digital PCR counts of the methane-xing 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 de ned 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 steadystate with respect to sulfate-methane dynamics, and hereafter refer to them as steady-state cores. Sul de pro les mirror the shape of the alkalinity curves, peaking at SMT depths. In addition, macroscopic SMTassociated mucoid bio lms consisting predominantly of ANME-1 36 , 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 uxes 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 pro les otherwise display considerable variability and dsrAB counts are typically low, between 10 4 -10 5 copies per gram bulk sediment (Fig. 2D).
GC1045 was sampled from the southern margin of GHM3, and GC1081 from the center of GHM4 (Fig. 1). Sulfate pro les from these cores show concave-up curvature, indicating that the methane-sulfate dynamics are not at steady-state, but likely re ect a recent increase in methane ux 14 (Fig. 3A). Porewater sulfate pro les show a rapid decrease in concentration down core and SMTs are well established. As estimated by our reduced modeling approach, total methane uxes 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 ux adequately t these data 14 . 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 rst 50 cm of GC1045 allows us to discount the possibility of oxic bottom water intrusion. These cores experiencing increases in methane ux 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 7 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 re ect a larger or more diverse sulfate-reducing community than in GC1081.
PC1029 was recovered from an established patch of frenulate siboglinid tubeworms (Oligobrachia sp. CPL clade, Fig. S2) whose chemosynthetic lifestypes are supported by sul de 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 sea oor at PC1029 (Fig. 4A) may be attributed to seawater in ltration (siboglinid bioirrigation or bubble-driven convection) and/or sul de oxidation from bacterial symbionts 39 . 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 ux 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 pro le with the greatest concentration gradient (10-15 cmbsf) yields a peak AOM rate on the order of 10 3 nmol cm -3 d -1 . This rate estimate would be increased signi cantly by accounting for siboglinid-driven pumping of bottom seawater sulfate, or sul de 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 ux 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 27 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 8 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 25,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 ux, 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 re ect 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 signi cant 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 in ux (Fig. 3D). Across samples for which we were con dently 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 in ux (Fig. 3C). Besides these major groups, other poorly-understood taxa include the Aminicenantes, Anaerolineae, and Phycisphaerae, all thought to be fermentative saccharolytic heterotrophs 41,42,43 .
Dehalococcoidia, also present, contain members capable of reductive dehalogenation 44 . Among all ASVs in our dataset, we identi ed ninety-nine whose relative abundances were signi cantly 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 ux, 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 ux, 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 2 = 0.06, p < 0.001), and when samples from above or below the SMT were considered separately (R 2 = 0.19, p = 0.003; R 2 = 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 pro les: the linear sulfate reduction (SR) zone, the nonlinear SR zone impacted by recent methane in ux, and below-SMT. Community structure also varied signi cantly 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 ux were more similar to below-SMT communities than those in linear SR zones, suggestive of recent adaptations to methane in ux (Fig. S4). Aside from methane regimes and redox zones, communities also varied signi cantly 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 27 that remains unconstrained.

Discussion
The presence of distinct regimes of methane ux 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 slowgrowing anaerobic fermenters of organic matter. Growth of bio lms 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 pro les indicate the methane front began around 290 years ago. As methane ux 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 in ltration driven by bubble-driven convection, bioirrigation by frenulate siboglinids, or reoxidation of sul de 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 uxes 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 45 . 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 uxes from steady-state cores fall within the middle 46 . At the seep site, PC1029, our methane ux estimate on the order of 10 2 mol m -1 yr -1 is several times the maximum of other modeled AOM rates at seep sites 46 , but less than the highest empirically measured AOM rate 6 . We acknowledge that the rate estimated from PC1029 is associated with large uncertainties, as we were not able to satisfactorily t 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 in ltration from the bottom water or sul de 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 uxes may indicate stable conditions over timescales of years 47 .
We now consider ANME doubling times at sites experiencing an increase in methane ux. 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 48 . These values complement the only other published estimate of in situ ANME doubling time, at approximately 100-200 days 47 . We can then estimate per-cell AOM rates across SMTs at an average of 0.65 pmol CH 4 d -1 , comparable to the 0.5-1.8 reported in a bioreactor where AOM was stimulated 49 . The observation that AOM rates and mcrA abundances in areas experiencing increasing methane ux peak at nearly the same depths suggests that the notoriously slow ANME doubling times may not present a signi cant lag (only 0-1 years) in the response of the benthic biogeochemical methane lter, in contrast with ndings from a prior modeling study from coastal ocean sediments with methane produced in-situ through acetotrophic methanogenesis 50 .
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 ndings 5 . 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 51 , and fewer multi-heme cytochromes thought to be involved in direct intercellular electron transfer 52 . Additionally, the lack of mer in ANME-1 53 , and the suggestion of met 54 or fae, hps, and adh 55 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 ux, ANME-1a and ANME-1b were present at nearequal 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 ux (Figs. 2, 3, 4). This may point towards a boom-and-bust cycle where methane in ux 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 re ect 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 pro le. Alternatively, this may indicate inactive relic communities, though ANME-1 may still be capable of AOM 56 or even methanogenesis 57 when starved of sulfate.
Co-occurrences between ANME-1b and SEEP-SRB2 have been reported 2 , 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 una liated with ANME in our dataset-are presumed to oxidize a wide variety of reduced hydrocarbons 35 . The presence of several potentially fermentative and saccharolytic clades like the Atribacteria, Aminicenantes, Anaerolineae, and Phycisphaerae may re ect alternate organic matter-dependent metabolic strategies that are interrupted by ANME and SRB when methane enters sulfate-rich porewaters. Macroscopic ANME-dominated bio lms found at two SMTs in GC1048 and GC1070 36 contained mcrA in concentrations of up to 7.6 x 10 10 g -1 . We hypothesize these bio lms re ect sediment regimes that have experienced steady methane and sulfate supply over many years. ANME bio lms have been described at SMTs in other subsea oor locations, often in fracturedominated cores 58 .
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 2 . Broadly, diversity decreases with depth, but only signi cantly across depths corresponding to peak AOM rates in high methane ux 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 ux (i.e. inhabiting non-steady-state SR zones of cores where methane ux 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 in uenced 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 59 .
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 60 and in a submarine mud volcano offshore Japan 61 . Four JS1 ASVs were identi ed across different regimes of methane ux 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-ux 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 ux, 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 32 .
The presence of differentially abundant ASVs at distinct regimes may re ect the sampling of comparatively shallow sediments at PC1029, and an in uence from macrofauna. Sulfur-oxidizing gammaproteobacterial symbionts of the siboglinid frenulate Oligobrachia have previously been reported in cold seeps 62 . Notably, there is an absence of Oligobrachia and a decreased prevalence of sea oor bacterial mats at GHM5, where steady-state cores were collected 21 . Despite the short (several km) distances between individual GHMs, many interdependent factors, such as physical disturbances, differences in uid ow regimes, and colonization of foundation species provide heterogeneity across seep ecosystems 63 .

Summary and prospective
In summary, our integrated approach allows us to detail regimes of methane transport where 1) steadystate 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 ux 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 ux. 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 64 . A modeling scenario that considers the dynamics of ANME growth may be of use in constraining estimates of marine subsurface methane ux into the hydrosphere. Global microbial methane lter e ciencies of 50-60% 5 have been used in modeling studies 19 , but seep sites display wide heterogeneity 6 . Our nding that mcrA gene copy numbers correlate positively with modeled AOM rates provides some justi cation for coupling these populations and their associated activities (Fig. S3), mirroring the coupling of methane uxes and transcripts of methane cyclers in peat soils 65 . Though microbial community data can provide explanatory power for predicting ecosystem processes, community changes do not always coincide with processes being measured 66 . 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 67 , or even apply these data into biogeochemical models to infer the presence of cryptic cycles 68 . 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 ares 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 69 . 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 Sub ghter 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-

Modeling
For the sites with non-steady-state porewater pro les (GC1045 and GC1081), we applied the same reduced modeling scheme as described previously 14 . Only sulfate and methane were considered in this reduced model, with AOM as the only reaction consuming both constituents. We assigned the diffusion coe cients of sulfate and methane to be 1.58x10 -2 and 3.01x10 -2 m 2 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 tting porewater pro les around the depth of the SMT for all cores 14 . Different lower boundary conditions of methane were assigned, which were constrained by the curvatures of the sulfate pro les. 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-ux lower boundary condition was used for sulfate.
For the sites with steady-state porewater pro les, we applied the comprehensive modeling scheme as described previously 14 . We coupled the CrunchFlow routine for this simulation with pro les 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 pro le, and therefore the same organic matter degradation rate, as in core GC920 from a previous study 14 . Seawater composition was used for upper and initial conditions for all solutes. Except for methane, no-ux 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) 14 . Depth-integrated uxes 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 pro le is also slightly bent as in the pro les of GC1045 and GC1081. To satisfactorily t this observation, we similarly applied a slightly intensifying methane ux when we modeled GC1048. However, the increase in methane ux we assigned for GC1048 is much smaller than in the two sites with obvious non-steady-state sulfate pro les. To eliminate concerns that only three porewater samples were collected above the SMTs which may lead to misinterpretation of the status of pore uid system (i.e., steady-state vs. increasing methane ux), we also t the porewater pro les from the four steady-state sites with model results for scenarios that assume an increase in methane ux (Fig. S5). Despite the scarcity of this sulfate data collected, we conclude it is di cult 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, Ampli cation, Sequencing, and Analysis DNA was extracted from sediments in a clean laminar ow 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. Brie y, V4 regions of bacterial and archaeal 16s rRNA genes were ampli ed in triplicate 25 ul reactions using universal 515-forward and 806-reverse primers 71 modi ed with dual-indexed Illumina sequencing adapters 72 . The thermal cycling protocol of Caporaso et al 2011 71 was followed without modi cations. After con rming ampli cation with agarose gel electrophoresis, triplicate PCR products were pooled and puri ed with a Qiagen QIAquick PCR puri cation kit. Amplicon concentrations were quanti ed with a Qubit uorometer 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 ampli ed and included in the sequencing run.
Working in R version 3.6.1, 16S rRNA amplicon data was processed with DADA2 73 (version 1.12.1) following an established pipeline 74 . Reads were denoised, chimeras removed, and taxonomies classi ed using version 132 of the SILVA nonredundant 16S reference database 75 . Sequences were aligned with DECIPHER 76 (version 2.12), and a phylogenetic tree was constructed using phangorn 77 (version 2.5.5).
Phyloseq 78 (version 1.28.0) was used to combine read count data with sample and taxonomy information. Sequences identi ed as Eukaryotes, Chloroplasts, or Mitochondria were removed, and the "combined" method of decontam 79 (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 80 (version 2.5-6), weighted Unifrac 81 distances calculated from a Hellinger-transformed ASV count table, and PERMANOVA tests were run to assess differences in community structure among groups. DESeq2 82 (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 83 and Luton 84 , 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 ampli cation 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 uorescence distributions using Quantasoft software (Bio-Rad). Threshold uorescence 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/sval ux.