In situ development of a methanotrophic microbiome in deep-sea sediments

Emission of the greenhouse gas methane from the seabed is globally controlled by marine aerobic and anaerobic methanotrophs gaining energy via methane oxidation. However, the processes involved in the assembly and dynamics of methanotrophic populations in complex natural microbial communities remain unclear. Here we investigated the development of a methanotrophic microbiome following subsurface mud eruptions at Håkon Mosby mud volcano (1250 m water depth). Freshly erupted muds hosted deep-subsurface communities that were dominated by Bathyarchaeota, Atribacteria and Chloroflexi. Methanotrophy was initially limited to a thin surface layer of Methylococcales populations consuming methane aerobically. With increasing distance to the eruptive center, anaerobic methanotrophic archaea, sulfate-reducing Desulfobacterales and thiotrophic Beggiatoaceae developed, and their respective metabolic capabilities dominated the biogeochemical functions of the community. Microbial richness, evenness, and cell numbers of the entire microbial community increased up to tenfold within a few years downstream of the mud flow from the eruptive center. The increasing diversity was accompanied by an up to fourfold increase in sequence abundance of relevant metabolic genes of the anaerobic methanotrophic and thiotrophic guilds. The communities fundamentally changed in their structure and functions as reflected in the metagenome turnover with distance from the eruptive center, and this was reflected in the biogeochemical zonation across the mud volcano caldera. The observed functional succession provides a framework for the response time and recovery of complex methanotrophic communities after disturbances of the deep-sea bed.


Introduction
The ocean seabed is the largest methane reservoir on Earth, comprising this climate-relevant gas in the form of semistable methane hydrates, as gas bubbles or dissolved in porewater. Globally, most of the methane rising from deeper subsurface layers is oxidized by methanotrophic microbial communities before it can reach the hydrosphere [1]. The methanotrophic communities in the seabed are diverse, but dominated by relatively few globally distributed types [2]. The thin oxic surface layer of methanerich sediments is often inhabited by aerobic methanotrophic bacteria of the Methylococcales [2][3][4][5]. Anoxic subsurface layers, where methane and sulfate overlap, are inhabited by consortia of anaerobic methanotrophic archaea (ANME) and their partner bacteria of the sulfate-reducing Desulfobacterales [6][7][8][9]. These methanotrophic communities, also referred to as the microbial methane filter, remove >90% of the methane in undisturbed continental margin sediments [1]. Methanotrophs also play an important role in methane removal at shallow [10,11] and deep-sea gas-emitting seep habitats [12,13]. Hence, only a small fraction of the seabed methane escapes from these sediments to the hydrosphere and atmosphere. However, the microbial methane filter at geologically highly dynamic seeps such as mud volcanoes has a lower efficiency, removing only 10−30% of the rising methane [14,15]. Understanding the causes for these different efficiencies, as well as the time scales needed for the establishment of an efficient methane filter, is crucial in order to assess the consequences of natural and man-made seafloor disturbances, such as rapidly dissociating hydrates [16,17], mud slides, eruptive mud volcanoes [14] or large oil spills [18][19][20].
Here we study the development of a deep-sea microbiome disturbed by seafloor mixing due to gas eruptions and mud slides at the actively gas-emitting Håkon Mosby Mud Volcano (HMMV) on the Norwegian continental slope. Marine mud volcanoes are seabed structures formed by upward migration of subsurface gasses together with fluids and sediments, from hundreds of meters to several kilometers depth by buoyancy and gravitational forces [21]. They are an important source of the greenhouse gas methane, globally emitting an estimated 27 Tg per year [22]. It has been speculated that the reduced efficiency of the microbial methane filter at mud volcanoes could be due to the low availability of electron acceptors, since the sediments are purged with anoxic subsurface fluids rising with the gas [14,23]. Other factors may be fluctuating temperatures, or frequent disturbances by mud mixing, which affect the growth of methanotrophs [24]. To investigate this further we compared the biogeochemistry and microbial community composition between recently disturbed, partially recovered, and undisturbed seafloor, using time-series observations and sampling of the Håkon Mosby in the framework of the deep-sea observatory "LOOME-Long term observations of mud volcano eruptions" (2003−2010). The hypotheses tested were (1) that the subsurface microbial signature of freshly erupted muds disappears with exposure to deep oxygenated seawater, (2) that freshly erupted muds lack complex methanotrophic communities and hence may have a low capacity to remove methane, and (3) that it needs years to develop complex cold-seep communities due to the slow generation times and cold temperatures.

Results and discussion
We investigated the seabed microbial community in mud flows of the HMMV (72°N, 14°44′E, 1250 m water depth) during research campaigns in August 2009 and September 2010. In this period, the long-term geophysical recordings of the LOOME observatory ( Fig. 1, S1; [25]) measured three major and 12 minor eruptions that occurred every 3−4 weeks. From this eruption pattern, detected by our deployed instruments and visual observations, we were able to infer an average mud flow velocity of 0.4 m day −1 across the HMMV center [25]. The mud velocity was used to convert sample distance to the eruptive central conduit into time. This space-for-time substitution approach used in ecological analyses of disturbances [26] allowed us to infer the spatiotemporal development of the methanotrophic microbiome. After visual seafloor inspection by ROV Quest and AUV Sentry, and biogeochemical sampling, we categorized the center muds into four zones (Figs. 1 and 2, S1; [27]). Zone 1 covered an area of about 50 × 90 m at the HMMV center. It consisted of fresh subsurface muds that were deposited during several gas eruptions recorded in 2009−2010 (Fig. 1c, Fig. S2; [25]). Zone 2 consisted of older subsurface mud flows southeast of the active center. The muds had a smooth, slightly rippled surface and were exposed to seafloor conditions for 1−2 years according to the morphology of the seabed and the measured flow velocity of the muds. Zone 3 were muds > 200 m away from the eruptive center, with thin thiotrophic mats, exposed to cold bottom waters already for 2−5 years. As a fourth zone we sampled the hummocky rim around the active center of HMMV. These sediments are stabilized by layers of hydrate at a few meters sediment depth [28]. They are not physically mixed by the center eruptions, as evident by comparative high-resolution mapping of the structure in 1996 [29,30], and by the dense coverage of long-lived siboglinid tubeworms (Fig. S2F), which were absent in zones 1-3.
Biogeochemistry and microbial methanotrophic rates change with time and distance from the eruptive center The fresh warm subsurface muds exposed during eruptions were saturated with methane as indicated by spontaneous in situ degassing, and had high concentrations of dissolved inorganic carbon (DIC), alkalinity and ammonium, indicative of a deep-subsurface origin of the pore fluid (Table 2, Fig. 2b, c). The fluids originate from a depth of up to 3 km below the seafloor, where the central conduit of HMMV is rooted [28,31]. The methane carbon isotopic signature of the dissolved gas in the warm sediments is similar to that of the surrounding gas hydrates with around δ 13 C = −60‰ (PDB), indicating a mixed thermogenic/biogenic origin in the deep subsurface [32]. Due to holes and cracks from degassing of the center muds (Fig. S2A), sulfate-containing bottom water percolated 5−10 cm into the seafloor. In the exposed surface muds 9 months after the main eruption, we measured methane oxidation (MOx-total aerobic and anaerobic methane oxidation) rates of up to 80 nmol ml −1 day −1 . These are low rates for seep ecosystems [6], especially when considering the high availability of methane and oxygen at the surface, or methane and sulfate in HMMV subsurface sediments. Sulfate reduction (SR) was not detectable, and the concave shape of the sulfate gradient ( Fig. 2b) is explained by the downward diffusion of sulfate against an upward flux of subsurface fluids [33]. In zone 2, which still showed ripples from the mud eruptions (Fig.  S2D), MOx rates were higher than in zone 1, peaking at 100 nmol ml −1 day −1 at the surface (Fig. 2c); SR rates peaked in the top few centimeters, but still no free sulfide was detected in the porewaters. Zone 3 is characterized by mats of sulfide-oxidizing bacteria (Figs. S1B, C, S2D, E). Here, the sulfate profile showed substantial consumption in the upper centimeters, and sulfide concentration reached almost 6 mmol l −1 (Fig. 2d). Integrated SR rates of zone 3 sediments measured during six expeditions between 2001 and 2010 ( Table 2) consistently showed the high rates that are typical for a well-established community of anaerobic methanotrophs (18 ± 4 mmol m −2 d −1 ; mean ± S.E.; n = 18; Table 2; [34]). These rates are around 60-fold higher than average sulfate reduction rates of nonseep impacted shelf sediments [35]. Sulfide production peaked in surface sediments of zone 3 (1.19 ± 0.15 mmol l −l ; n = 110), whereas sulfide concentrations were low in zones 1 and 2 ( Table 2) confirming that AOM was not established in fresh sediments. Porewater analyses as well as methane oxidation and sulfate reduction rate measurements performed on these six expeditions support the lateral zonation across the caldera. This zonation is in accordance with the geophysical model, in which centrally rising mud and gas are laterally transported and eventually sink in upon degassing [25]. Our biogeochemical results strongly suggest that over the period of at least one decade AOM continuously established in muds of the same distance to the eruptive center-i.e. muds of a similar age-despite a continuous mud flow.
Subsurface communities get transported to surface sediments due to mud volcanism We analyzed archived subsurface samples from the prestudy phase in 2003 and additionally assessed the microbial community composition of ten surface (S) and five Fig. 2 Biogeochemistry of surface sediments of HMMV and of a reference site outside of the mud volcano. In the reference sediment (a) outside of the HMMV caldera around 500 m upslope of the mud volcano, dissolved inorganic carbon (DIC), sulfate and alkalinity showed typical background concentrations, ammonium (NH 4 ) and sulfide (H 2 S) were not detected, and there was no measurable methane oxidation (MOx) or sulfate reduction (SR). Sediments at the center of HMMV (b) show an upward transport of sulfate-depleted subsurface fluids enriched in DIC and NH 4 and show low MOx rates. SR is first detectable in sediments of zone 2 (c) and H 2 S production is first detectable in aged sediments of zone 3 (d). Note: MUC827 is a parallel core of Sample 9 (MUC823) originating from the same dive and same area. All profiles were measured on the same expedition in 2010. For details as well as pore-water concentrations and rates of other samples in zones 1−4, see Table 2 and ref. [61] subsurface (D) sediment samples across the same zones sampled after the eruptions. The microbial communities in the warm (10−20°C) subsurface sediments (3.8−2.5 m below seafloor) of the HMMV caldera were characterized by a relatively low archaeal and bacterial alpha diversity, and showed a low community turnover, i.e. replacement of microbial taxa across zones (Figs. 3, 4, Fig. S3). These simple and homogeneous communities support the geophysical model of a uniform, warm subsurface mud layer filling the central chimney of HMMV [25,28,36]. Interestingly, the subsurface communities were most similar to the surface communities of zone 1 (Fig. 3b, c; Table S2), which is in accordance with the observed upward transport and deposition of the subsurface sediments by gas eruptions (Fig. S2A). Subsurface and also surface sediments of Zone 1 contained typical heterotrophic deep biosphere clades, such as Bathyarchaeota (Miscellaneous Crenarchaeotic Group), Chloroflexi and Atribacteria (candidate division JS1) [37,38] as well as Peptococcaceae and methanogenic Methanosaeta (Fig. 4, Fig. S4). The latter four clades were suggested to form a syntrophic network, degrading proteins and fatty acids under methanogenic conditions [39]. Bathyarchaeota, which comprise organisms that also degrade detrital proteins [40], greatly dominated the archaeal community in the freshly deposited muds (Fig. 4, S4). We detected many genes involved in fermentation and methanogenesis (Fig. S4, S5A) in the subsurface metagenomes, while genes for sulfate reduction, sulfur oxidation or methane oxidation were absent or very rare. The subsurface communities in the central mud conduct of HMMV may be fueled by organic compounds transported with rising porewater fluids from the deep subsurface. The subsurface and surface community of zone 1 contained very few ANME-3 sequences (<1% relative sequence abundance), and few genes affiliated with Methylococcales, sulfatereducing bacteria (SRB) and sulfur-oxidizing bacteria (SOB) (Fig. S4), in accordance with the biogeochemical profiles of the freshly exposed muds (Fig. 2, Table 2).
By comparative sequence analysis, we tested if the subsurface clades deposited at the surface, such as Atribacteria [41], would persist in the surface muds with increasing exposure and distance from the central mud conduit. The subsurface-derived clades decreased in relative sequence abundance with increasing distance to the center. In zone 3, these clades were barely detectable in the sediments. The subsurface microbial signature was replaced by that of typical seep communities within around 2 years. Only 3−7% of the operational taxonomic units (OTUs; at 98% 16S rRNA gene V4-V6 sequence identity) were shared between subsurface and surface muds (Fig. 3, Table S2) and even less (1% shared OTUs) with the reference site that was characterized by typical oligotrophic deep-sea sediment organisms, such as Xanthomonadales and Thaumarchaeota [42]. This diversity pattern was also confirmed based on the turnover of gene families in the community metagenomes (Fig. 5a). Only 18−35% of the gene families found in the subsurface metagenome of zone 1 (sample 11) were found in the surface metagenomes, while 66 and 77% gene families were shared between sample 11 and the other two subsurface metagenomes (Table S6). Similarly, the surface metagenomes shared between 52 and 79% gene families among each other (60 ± 9.8 %, mean ± SD, n = 6), but only 18−48 % with the subsurface metagenomes (33 ± 9.4%, mean ± SD, n = 12). This suggests that subsurface-derived genes, and hence subsurface community functions, were rapidly depleted in the surface muds. Subsurface microbes have longer generation times compared to surface bacteria [43][44][45]. In addition, they were likely repressed by the exposure to the seafloor conditions; i.e. cold temperatures of −1°C and high oxygen concentrations (>280 μM; [46]), and were eventually overgrown by others.

Surface communities developing with distance from the eruptive center
Using the mud flow velocities measured by LOOME observatory and space-by-time substitution, we investigated how the microbiome composition and metabolic capabilities of the surface communities would shift across the different zones with time, and if these changes were in accordance with the changes in biogeochemical rates ( Fig. 2; Table 2). The analysis of 16S rRNA gene amplicons confirmed that community structure was significantly different between the zones, as tested by ANOSIM (R Arch = 0.7, p Arch < 0.01; R Bac = 0.5, p Bac < 0.01, Fig. S3). Also, the communities were more similar between adjacent zones than between zones that were further apart, e.g. zone 1 shared more OTUs with zone 2 than with zone 3 (Table S2). Total cell counts increased with distance to the active center (Figs. 3, 6a). Relative cell abundances of key microbial clades involved in methanotrophy and sulfate reduction also changed between the zones (Fig. 6b), increasing twofold from zone 1 to 2. Further, we observed a strong increase in microbial richness and evenness in zone 3 compared to surface sediments of zone 1 and 2 (Fig. 3). Microbial diversity continued to increase in the zone 4 sample and peaked in the reference site. Around 85% of archaeal and bacterial OTUs of the freshly erupted muds of zone 1 were replaced within 1−2 years of exposure to surface conditions (Fig. 3, Table S2). Yet, the communities of zones 1 and 2 still shared key organisms (e.g. Methylococcales and Desulfobacterales), especially those samples taken at a similar distance to the central conduit (e.g. Samples 3 and 5, Fig. 3). Zone 3 sediments reached total cell numbers of 3 × 10 9 cells ml −1 , similar to the community density in the stable hummocky rim (zone 4). Zone 4 harbored the most complex microbial communities at HMMV, with the highest diversity and evenness (Fig. 3). The overall community development observed at the level of 16S ribosomal RNA gene diversity and turnover was also confirmed by the diversity and turnover of metabolic gene families ( Fig. 5; Figs. S4, S5). Rarefaction curves showed that the expected number of gene families was highest in the consolidated muds of zone 3 and the reference site (Fig. 5d) and generally lower in freshly exposed sediments and subsurface muds. Similarly, the total number of observed gene families was highest in the metagenome of surface sediment from the reference site and low in the sediments of the active center. At the same time the number of unique gene families that exclusively occurred in one metagenome was high in the active center and lowest in the reference site (Fig. 5a). This indicates that with distance from the active center the surface sediments are becoming increasingly diverse, but also increasingly similar on the level of community functions.

Development of aerobic methanotrophic populations
In zones 1 and 2, OTUs affiliated with aerobic methanotrophic Methylococcales were abundant in the top cm of surface sediments exposed to the oxygen-rich cold bottom waters, as shown by relative sequence and cell abundances (Figs. 4, 6) and a high number of 16S rRNA and pmoA gene sequences in the metagenome (Fig. S4, S5B). The Methylococcales reached 5 × 10 8 cells ml −1 sediment (3.1 ± 1.9 × Fig. 3 Alpha and beta diversity, and total cell abundance across HMMV sediments. a Archaeal and bacterial diversity was determined using OTUs (Operational taxonomic units at 98% 16S rRNA V4-V6 gene sequence identity, corresponding to the recommended taxonomic threshold for microbial species [74]). OTU data were used to assess observed richness, Inverse Simpson diversity, and Chao1 estimated richness. Note: The axes for archaeal and bacterial values differ by one order of magnitude. Total cell numbers were determined by Acridine Orange Direct Counts and were integrated over the top 10 cm sediment depth. Dashes denote missing data points. b Shift of microbial community structure in HMMV sediments as visualized by nonmetric multidimensional scaling (NMDS) using relative sequence abundance of archaeal and bacterial OTUs. Color indicates the sample origin (Subsurface = gray polygon; Zone 1 = red; Zone 2 = purple; Zone 3 = light blue; Zone 4 = dark blue; Reference site = brown). The percentages of microbial OTUs that are shared between zones (numbers next to arrows) are based on presence−absence data, i.e. showing that only 1% of OTUs that are present in the subsurface are also found in the surface sediment of the reference site. The microbial communities of the subsurface and of zones 1−3 were all significantly different from each other (ANOSIM based on presence/absence data: R = 0.7, p = 0.001). Zone 4 and the reference site could not be included in the ANOSIM as there was only one sample retrieved. c Dendrogram showing hierarchical clustering of the samples based on archaeal and bacterial OTUs. Fig. 4 Succession of microbial clades in sediments of HMMV based on relative abundance of OTUs. Together the bars of one sample add up to 100% archaeal and 100% bacterial relative sequence abundance focusing on functionally relevant clades involved in the methane and sulfur cycle. Each bar in these panels shows the top 2−3 most sequence abundant OTUs and all remaining OTUs (Other) that belong to this subset. A complete table with archaeal and bacterial OTUs, as well as a summary table of relative abundances of key populations (e.g. ANME-3) is available at PANGAEA, see ref. [61] 10 8 cells ml −1 sediment, mean ± SD, n = 6; Table 2) already after a few months of exposure of the subsurface muds (Fig.  S6). This peak in cell abundance of aerobic methylotrophs in fresh muds was observed in previous expeditions [14,15,47], suggesting that these organisms rapidly colonize fresh HMMV muds in general. In the thin sediment surface layer where oxygen is available, Methylococcales can respond faster to the high supply of methane, as their higher energy yield supports faster growth rates compared to those of anaerobic methanotrophs [48][49][50][51]. Also, aerobic methylotrophs could relatively rapidly colonize freshly exposed gassy muds, as they can disperse with bottom waters [2,52]. Assuming that representatives of this group would settle on the freshly deposited muds by sediment resuspension across the mud volcano, and applying a mud transport rate of 0.4 m per day between sediments of zones 1 and 2 as observed in this time period [25], the aerobic methylotrophs could have achieved a net growth rate of 0.01 day −1 corresponding to a doubling time of 60 days. This rate is similar to that calculated for Methylococcales in arctic, boreal swamps (0.02 day −1 ; [53]). In comparison, cold-water anaerobic methanotrophs commonly show growth rates of around (0.003 day −1 ) [48]. This finding supports previous hypotheses, that aerobic methanotrophs dominate surface sediments of emerging methane leaks and young seep systems [24,54,55]. However, spatially, they can only occupy a small niche due to limited oxygen penetration into the seafloor [12,15,54].

Development of anaerobic methanotrophic communities and sulfur-oxidizing bacteria
Cell counts and sequences showed that anaerobic methaneoxidizing archaea (ANME) and their sulfate-reducing partner bacteria were rare in the freshly exposed center sediments of zone 1 (Figs. 4, 6; Fig. S7). Free-living sulfate-  Tables S5  and S6, respectively. Turnover of gene families between metagenomes is visualized by dendrogram (b) and nonmetric multidimensional scaling ordination (c). d Rarefaction indicates that most gene families, and possibly metabolisms, that are present at HMMV were detected reducing bacteria (SRB) increased tenfold in relative abundance from zone 1 to zone 2, the latter harboring 2.4 ± 1.5 × 10 8 SRB cells ml −1 sediment (Table 2). Zone 3 had similar SRB counts as zone 2, but in addition harbored large numbers of ANME/SRB consortia (Fig. 6, Table 2). The sulfide that was produced by these AOM consortia (Fig. 2d, Table 2) supported the growth and establishment of sulfuroxidizing bacteria (SOB), forming white mats covering zone 3 (Fig. S1, S2D, E). In comparison, the thiotrophs were rare in sequence abundance in sediments of zone 1. Their relative sequence abundance increased between zones 2 and 3, where an increasing amount of oxidative dsrAB genes were detected in the metagenomes (Figs. S4, S8). The archaeal community of zone 3 was dominated by a single ANME-3 OTU. This OTU had a relative sequence abundance of below 0.1% in freshly exposed subsurface muds of zone 1, 2−30% in zone 2 and comprised 63−88% of all archaeal sequences in zone 3, indicating a significant increase of the population with time and distance from the mud and gas conduit, but also suggesting a slow doubling time of 100−200 days. This estimated doubling time of ANME-3 in situ corresponds to a growth rate of 0.003 −0.006 day −1 which is in good agreement to the rates estimated by in vitro experiments with psychrophilic ANME-2 (0.003 day −1 ; [48]). All populations of ANME, SRB and SOB (Figs. 4, 6) and their metagenomic signatures, especially the relative abundances of mcrA and dsrAB genes (Figs. S4, S5A, S8) showed the same pattern. They were rare in subsurface and freshly exposed muds, and became abundant in surface sediments with increasing distance from the center and thus exposure time of the muds. This demonstrates that slow-growing, and initially rare hydrocarbon-consuming microorganisms are able to out-compete others at cold seeps when the conditions are favorable and the time-scale permits [56,57]. In the undisturbed zone 4, ANME-2a and ANME-3 archaea had similar relative sequence abundances (Fig. 4), indicating that these consolidated sediments harbor niches for more were assessed with DAPI (white bars) and compared to bacterial cell abundance (probe: EUB338-I-III, gray bars). Replicate samples were available for zones 1−3 (zone 1: n = 3, zone 2/3: n = 2), shown as multiple bar pairs. b Relative cell abundances based on single-cell counts using CARD-FISH with specific probes for Bacteria (probe EUB338 I-III), Methylococcales (probe MTMC-701 and competitor probes), Desulfosarcina/Desulfococcus (probe DSS658), Archaea (probe Arch915), and ANME-3 (probe ANME3-1249 and helper probes). Bacteria that did not overlap with Methylococcales or DSS are denoted "other Bacteria". Archaea that did not overlap with ANME-3 are referred to as "other Archaea". "Cells without probe signal" were only stained by the nucleic acid stain DAPI and not by general archaeal or bacterial probes. Relative abundances were averaged over the top ten centimeters; all three layers are shown in Figure S7. Note: The relative cell abundances for ANME-3, Desulfosarcina/Desulfococcus and Methylococcales are underestimated as these clades formed cell aggregates that were not included in the single-cell counts. Probe details are given in Table S3; detailed values and cell counts are archived [61] ANME clades, as compared to the disturbed sediments of the caldera, which were highly dominated by ANME-3. In zone 4, very little methane reaches the surface sediment due to active methanotrophic communities at the roots of the tubeworms in~60 cm depth [47], resulting in very low methane oxidation rates at the surface ( Table 2). The low temperature and stability of the sediments as well as the presence of hydrates below 0.5−1 m sediment depth may explain why the investigated surface sediment are of low activity and share 12% OTUs with surface sediment of the nonmethane reference site (Fig. 3, Fig. S3). Together, the observation of community succession on mud flows at HMMV match previous experiments with wood and whale falls showing that deep-sea methanogenic, methanotrophic and thiotrophic clades need years to develop functional communities on allochthonous surfaces and energy supplies [58][59][60].

Conclusion
Here we studied the development of a deep-sea methanotrophic microbiome on mud flows of a methane-emitting mud volcano in situ. We were able to sample the communities of freshly exposed subsurface muds, and compare them to their source community, as well as to increasingly developed seep and nonseep communities outside the caldera. Changes in biogeochemical rates with increasing distance to the eruptive center of the mud volcano were mirrored by changes in the corresponding metabolic genes and cell counts of the respective clades. At the level of the whole microbial metagenome, both 16S rRNA sequence turnover as well as the diversity of metabolic gene families showed a pattern of increasing complexity with increasing development of the methanotrophic assemblages, supporting a rich and diverse bacterial and archaeal community. We were able to confirm our initial hypothesis based on biogeochemical measurements that subsurface communities of bacteria and archaea are replaced by pioneering aerobic methanotrophs and later complex anaerobic methanotrophic and thiotrophic communities. Even when electron donors and acceptors are not limiting, the succession of benthic deep-sea bacterial and archaeal populations may need years, before the typically high diversity and evenness of deep-sea sediment communities is reached. Our findings indicate that loss of seafloor integrity-in this case by gas eruptions and mud mixing-and thereby the local decline of active and complex methanotrophic communities can explain the low efficiency of methane consumption that is globally observed at active mud volcanoes. Over several years, a seep microbiome can develop from initially rare populations to a complex community, in this case study evidenced by increasing cell and sequence numbers, and increasing diversity, of aerobic and anaerobic methanotrophs and thiotrophs. The observed functional succession provides insights into the response time and recovery of complex microbial communities to natural and anthropogenic disturbances in the deep sea.

Sampling sites
Surface sediment samples (0−10 cm) at HMMV are exposed to the cold Arctic bottom water, and generally have an ambient temperature of −1°C. They were recovered either by TV-guided Multicorer or by push cores using the remotely operated vehicle Quest (Marum, University Bremen) ( Table 1). Subsurface sediments of all zones (>2 m below sea floor) were obtained by gravity corer. Their temperature at 4 m below the seafloor ranged from around 15°C in the center to around 3°C at the hummocky rim [36]. After recovery, sediments were immediately subsampled in a refrigerated container (0°C) and further processed for biogeochemical analyses or preserved at −20°C for later DNA analyses. Further details to the geographic locations, dates of sampling, and all contextual data are provided in the supporting information Table 1 and in the public archive for Earth and environmental data PAN-GAEA; see ref. [61].

Biogeochemistry
Porewater and turnover rates were measured in surface sediment cores obtained in 2010 using methods described previously [62]. We show four profiles in detail (Fig. 2 Table 2 and can be accessed from the data publisher PANGAEA; see ref. [61]. Biogeochemical parameters of sediments from 2003 and 2009 have been reported previously [14,15]. The pore water was extracted with Rhizons in 1 cm resolution and immediately fixed in 5% zinc acetate (ZnAc) solution for sulfate, and sulfide analyses. The total sulfide concentrations (H 2 S + HS − + S 2− ) were determined using the diamine complexation method [63]. DIC and alkalinity were measured using the flow injection method (detector VWR scientific model 1054) [64]. Nutrients were determined with a Skalar Continuous-Flow Analyzer [65]. Sulfate reduction (SR) and anaerobic oxidation of methane (AOM) were measured ex situ by the whole core injection method [66] as previously described [67,68]. Refer to SI for details on biogeochemical analyses. All cores except the reference site were degassing methane after retrieval; hence we could not measure true in situ methane concentration in Table 1 Sample overview Sample ID Sample  n.d.  Cell counts are given as cells ×10 8 ml −1 sediment pore waters. It was previously estimated that the in situ concentrations of methane in the gassy HMMV center could reach 100 mM [23].
16S rRNA gene V4-V6 amplicon pyrosequencing DNA extraction was done in duplicates using 1 g sediment each and a commercially available extraction kit (Ultra-Clean Soil DNA Isolation Kit, MoBio, Carlsbad, CA We assessed the quality, size and concentration of PCR products on a Perkin Elmer Caliper GX. Reads were demultiplexed and barcodes removed for submission. Sequence reads were submitted to a rigorous quality control procedure using mothur v30 [69] and a routine [2,70] that included denoising of the flow grams [71], single-linkage preclustering [72] and the removal of chimeras [73]. Sequences were clustered at 98% ribosomal RNA gene V4-V6 sequence identity-corresponding to the recommended taxonomic threshold for microbial species [74]-and were taxonomically assigned using the SILVA taxonomy (SSURef v119, 07-2014 [75]). Further information about sequencing datasets and contextual data are available at PANGAEA [61].

Analyses of V4-V6 amplicon data
Relative abundance of archaeal and bacterial OTUs (operational taxonomic units clustered at 98% sequence identity) is based on the original mothur output (Table S1).
To calculate Inverse Simpson diversity indices and Chao1 Richness [76] the OTU abundance tables were rarefied to account for unequal sampling effort using 300 (Archaea) and 1000 (Bacteria) randomly chosen sequences without replacement using mothur. Bray−Curtis dissimilarities [77] between all samples were calculated and used for twodimensional nonmetric multidimensional scaling (NMDS) ordinations with 20 random starts [78]. All analyses were carried out with the R statistical environment and the packages vegan [79], labdsv [80], as well as with custom R scripts (for details see SI).

Shotgun metagenomics
DNA extraction using 3 g sediment (pooled from 0 to 10 cm sediment depth) was performed manually as previously described [81]. DNA was sheared using a Covaris and libraries were constructed with the Nugen Ovation Ultralow Library protocol and were amplified for 10−11 cycles. The amplified product was visualized on an Agilent DNA1000 chip or Caliper HiSens Bioanalyzer assay.
Libraries were pooled at equimolar concentrations based on these results and size selected using a Sage PippinPrep 2% cassette. The final library pool had an average insert size of 170 bp with~25−30 bp partial overlap between pairs of reads. It was quantified using a Kapa Biosystems qPCR library quantification kit, then sequenced on the Illumina HiSeq1000 in a 2 × 101 paired-end sequencing run using dedicated read indexing. The samples were demultiplexed with CASAVA 1.8.2. Details on the library output are given in Table S4. Further information about sequencing datasets and contextual data are available at PANGAEA [61,82].
Ribosomal and metabolic gene reconstruction from metagenomic data 16S rRNA and metabolic gene abundances as well as gene reconstructions were generated using a novel, modified version of the phyloFlash pipeline (https://github.com/ HRGV/phyloFlash) called funcFlash for metabolic genes. In brief, the generated reads were mapped with BBMap at minimal global nucleotide identity of 70% against curated nucleotide databases: The SILVA SSURef v119 database [75], a published dsrAB gene database (http://www. microbial-ecology.net/db_download/dsr_v3.zip, [83]) and two newly generated pmoA and mcrA gene databases that are publicly available at PANGAEA; see ref. [61]. The mapped read pairs were counted when at least one read had a positive mapping. Full-length (>70% of the target length) genes were assembled with SPAdes [84] or reconstructed with EMIRGE [85]. The mapping of reconstructed metabolic genes to curated databases using funcFlash can be used to distinguish, whether a gene of interest is affiliated with organisms performing the reductive or oxidative pathway. In our case this new analysis allowed us to distinguish between dsrAB genes from sulfate reducers and from sulfur oxidizers, as well as between mcrA genes from methanogens and from methanotrophs.

Multivariate analyses of metabolic gene families from metagenomic data
We used metagenomic data from seven sites to investigate richness, abundance, and turnover of gene families across the different zones of HMMV. Each metagenome was filtered to remove low-quality and/or short reads. Raw reads were merged into paired reads with BBMerge [86]. To enable the comparison of gene diversity across sites, each metagenome was subsampled to 10 6 paired reads using the BBMap "reformat" tool. Subsampled reads were analyzed with humann2 [87]. Here, the reads were subjected to a translated nucleotide search against species-level clusters of nonredundant gene families of the UniRef50 database [88] using DIAMOND [89]. Gene abundances were normalized according to gene length, and then conjoined to obtain a gene × metagenome table, analogous to an OTU × sample table, with gene families as rows and metagenomes as columns. This matrix allowed us to investigate diversity using metabolic gene families (functions) rather than the commonly used ribosomal genes (taxonomy). The gene × metagenome table was subjected to multivariate analyses (data reduction, hypothesis testing, visualization) based on the R packages vegan, labdsv, UpSetR [90] and customized R scripts. Prior to the analyses we removed gene families with less than ten read hits cumulated across all metagenomes, to focus on abundant gene families and to minimize the influence of spurious hits. NMDS ordinations, calculated percentages of shared gene families and the UpsetR diagram are based on a presence/absence matrix. Diversity indices and rarefaction curves were calculated with a subsampling-based iterative approach using abundance information. Refer to SI for details.
Nucleotide sequence accession and contextual data availability Cell counts and catalyzed reporter deposition fluorescence in situ hybridization (CARD-FISH) Total numbers of single cells were determined using acridine orange direct counts according to the protocol published elsewhere [91]. CARD-FISH was performed as previously described [54] with the following modifications. 4−6 µl of 25-fold diluted sediment were used for filtration. Archaeal cell walls were permeabilized with 0.1 M HCl for 2 min to detect ANME-3 cells, or Proteinase K solution (15 µg ml −1 (Merck, Darmstadt, Germany) in 0.05 M EDTA (pH 8), 0.1 M Tris-HCl (pH 8), 0.5 M NaCl) for 2−4 min at room temperature for all other archaea. Bacterial cell walls were permeabilized with lysozyme solution (1000 kU/ml) for 60 min at 37°C. Cells were stained with DAPI (1 µg/ ml), embedded in mounting medium and counted in 40−60 independent microscopic fields using an Axiophot II epifluorescence microscope (Carl Zeiss, Jena, Germany). Cell numbers of dense aggregates were estimated semiquantitatively as previously described [47]. A complete list of oligonucleotide probes, helpers, and competitors used in this study is provided (Table S3). A summary of the results from all 260 CARD-FISH experiments is publicly available at PANGAEA; see ref. [61]. source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons. org/licenses/by/4.0/.