Anaerobic methanotrophic communities thrive in deep submarine permafrost.

Thawing submarine permafrost is a source of methane to the subsurface biosphere. Methane oxidation in submarine permafrost sediments has been proposed, but the responsible microorganisms remain uncharacterized. We analyzed archaeal communities and identified distinct anaerobic methanotrophic assemblages of marine and terrestrial origin (ANME-2a/b, ANME-2d) both in frozen and completely thawed submarine permafrost sediments. Besides archaea potentially involved in anaerobic oxidation of methane (AOM) we found a large diversity of archaea mainly belonging to Bathyarchaeota, Thaumarchaeota, and Euryarchaeota. Methane concentrations and δ13C-methane signatures distinguish horizons of potential AOM coupled either to sulfate reduction in a sulfate-methane transition zone (SMTZ) or to the reduction of other electron acceptors, such as iron, manganese or nitrate. Analysis of functional marker genes (mcrA) and fluorescence in situ hybridization (FISH) corroborate potential activity of AOM communities in submarine permafrost sediments at low temperatures. Modeled potential AOM consumes 72-100% of submarine permafrost methane and up to 1.2 Tg of carbon per year for the total expected area of submarine permafrost. This is comparable with AOM habitats such as cold seeps. We thus propose that AOM is active where submarine permafrost thaws, which should be included in global methane budgets.

Introduction

Terrestrial permafrost landscapes, which developed during glacial cold periods, are known to be a large reservoir of organic carbon (~1300 Pg)[1]. Permafrost thaw and the following microbial production of carbon dioxide and methane from liberated organic matter may act as a positive feedback to climate warming[2]. Methane has 34 times higher global warming potential over a 100 year period[3] and is a more critical greenhouse gas than carbon dioxide. The amount and the release rates of methane are not well constrained although they are critical for evaluating future climate change.

Several Arctic sources of methane have been identified, including methane bursts during soil freezing[4], thermokarst lakes[5], lakes and ponds[6], wetlands[7], gas hydrates[8] and submarine permafrost[9]. Submarine permafrost on continental shelves of the Arctic Ocean is a consequence of the inundation of terrestrial permafrost by sea water during the Holocene marine transgression[10,11]. Coastal submarine permafrost froze under subaerial terrestrial conditions in alluvial/fluvial settings and has remained frozen since then[12,13]. Submarine permafrost is much more susceptible to thawing than permafrost on land, because of overlying warm marine water causing diffusion of salt water into the sediments from the top, and because of geothermal heat flux from below[9]. Based on thermal modeling of permafrost development over glacial/interglacial cycles, submarine permafrost is likely to have persisted in deep sediment layers for hundreds of millennia[14]. Over those millennial timescales submarine permafrost reaches its freezing point, between −2 and −1 °C, after which it slowly starts degrading due to the influence of high salt concentration[10,15]. Our knowledge of evolving carbon pools and carbon turnover in submarine permafrost is, however, scarce.

Arctic submarine permafrost regions have the potential to emit large amounts of methane to the atmosphere[16], however flux measurements are controversial. First estimates suggest a release of 8.0 Tg of methane per year from the East Siberian Arctic Shelf (comprising Laptev Sea, East Siberian Sea and the Russian part of the Chukchi Sea), which would equal the release of methane from the entire world ocean[9]. In contrast, recent estimates for the Laptev and East Siberian Seas are almost three times lower (2.9 Tg y−1)[17]. Methane is mainly released from the sediment but the source is unclear[18]. It is expected that permafrost degradation can create pathways for the release of gas captured within submarine permafrost sediment or within gas hydrates underneath[16]. Low δ13C-CH4 values indicate a biogenic or thermogenic origin of methane released from submarine permafrost[13,19]. A previous study demonstrated a drop of methane concentrations at the permafrost thaw front deep inside the sediments, which coincided with increasing δ13C-CH4 values from −70‰ to −35‰ that points to anaerobic oxidation of methane (AOM)[13]. Sulfate penetrating from the seabed into the sediment created a deep sulfate-methane transition zone (SMTZ) similar to other geological formations in the subsurface[20], with sulfate being a potential electron acceptor for AOM[21].

Most studies on AOM in the marine system have focused on ANME archaea that couple the oxidation of methane to the reduction of sulfate via a syntrophic lifestyle with sulfate-reducing bacteria (SRB)[22-24]. Three main marine clades, ANME-1a/b, ANME-2a/b/c, and ANME-3 have been identified so far[25]. These clades are associated with sulfate-reducing bacteria of the genera Desulfosarcina, Desulfococcus, Desulfobulbus and Desulfofervidus, which belong to the class Deltaproteobacteria[25] and of the phylum Thermodesulfobacteria[26]. Besides the marine ANMEs there is evidence of terrestrial AOM driven by the ANME-2d clade[27,28]. ANME-2d sequences were detected in wetland and permafrost habitats[29-31] so they might mitigate the release of methane in these environments. A recent study on wetlands showed that even low amounts of sulfate (µM) are sufficient to couple AOM with sulfate reduction in terrestrial environments[32]. In addition, alternative electron acceptors for AOM such as iron[33,34], manganese[33] as well as nitrate[27,35] and humic substances[36,37] have been suggested. To date, the microbial mitigation of methane emissions from permafrost environments is only evident in aerobic soils and sediments[38,39], while communities that oxidize methane anaerobically in thawing permafrost in the deep subsurface remain unexplored[8]. We hypothesize that active marine and terrestrial ANME clades exist in thawing submarine permafrost similar to microorganisms found in terrestrial enrichments[27,35] or marine SMTZ[40,41], and that these clades could function as an efficient methane filter for a diffusive methane release.

We investigated AOM in two deep submarine permafrost cores (48 m and 71 m) from the Siberian Shelf with different stages of submarine permafrost thaw[42] that were inundated for about 540 and 2,500 years, respectively[11,13]. To identify microbial assemblages involved in AOM, we performed high-throughput sequencing of archaeal and bacterial 16S rRNA genes and of the functional marker gene methyl co-enzyme reductase (mcrA) that is encoded by both methanogens and anaerobic methanotrophs. We quantified mcrA and dissimilatory sulfate reductase (dsrB) gene markers by quantitative PCR (qPCR) and correlated their abundance with the 16S rRNA gene abundance of AOM-related microorganisms. We visualized AOM-related microorganisms by catalyzed reporter deposition fluorescence in situ hybridization (CARD-FISH) and constructed ANME-specific 16S rRNA gene clone libraries to resolve their diversity in submarine permafrost sediments. Pore water geochemistry i. e., concentrations of methane, nitrate, sulfate, iron, and manganese, was analyzed to verify potential horizons of methane consumption. Analysis of stable carbon isotope geochemistry of methane and microbial membrane lipid biomarker were performed to evaluate the activity of AOM communities in situ.

Results

Isotopic evidence of AOM and potential electron acceptors

Methane concentrations were low (average 48 µM) in the ice-bonded permafrost section of the Cape Mamontov Klyk core (hereafter C2) inundated for ~2,500 years[13] and peaked at approximately 52 m below seafloor (bsf) with a concentration of 990 µM (Fig. 1A). In contrast, methane concentrations in the ice-bonded permafrost section of the Buor Khaya core (hereafter BK2) were on average 8 times higher than in C2 (Fig. 2A). δ13C values of methane in the ice-bonded permafrost of BK2 showed values between −71 and −53‰ typical for a biological origin[43] (Fig. 2A). Above the ice-bonded permafrost, δ13C values of methane ranged from −37 to −30‰ within a few meters on top of the thaw front and stayed constant to a depth of about 12 mbsf[13]. Above 12 mbsf δ13C-values of methane showed high dynamics with values ranging between −67 to −32‰ (Fig. 2A). δ13C values of methane in the C2 core were only measurable in parts of the ice-bonded permafrost and showed 13C values between −72‰ to −37‰ (Fig. 1A). Independently of one another, both cores showed opposing dynamics of methane concentration and δ13C values with a clear signal of AOM at the SMTZ of BK2. Estimated release rates of methane from the ice-bonded permafrost of BK2 together with the modeled fraction of methane that was oxidized yielded potential AOM rates of 1.7–2.1 nmol cm−3 d−1 ± 0.9–1.2 nmol cm−3 d−1 [13] in the SMTZ.

Figure 1

Pore water profiles of methane, nitrate, sulfate, manganese and relative archaeal abundance in the Mamontov Klyk core C2. (A) Methane concentrations are shown in black and corresponding δ13C-values in magenta[19]. (B) Nitrate concentrations are shown in blue, sulfate concentrations in red, manganese concentrations in green, and iron concentrations in light blue. (C) Purple dots on the left show depth location in the core. Relative abundance of archaeal 16S rRNA gene sequences are shown as bar plots. Colors of bars refer to the taxa in the legend below figure. Numbers on the right refer to the exact depth in the core. ANME = ANaerobic MEthanotrophic archaea, DHVEG-1 = Deep Hydrothermal Vent Euryarchaeotal Group 1, DSPEG = Deep Submarine Permafrost Euryarchaeotal Group, MBG-B = Marine Benthic Group B, MBG-D = Marine Benthic Group D, MCG = Miscellaneous Crenarchaeotal Group. The shaded area represents the permafrost degradation zone, the red area a sulfate-methane transition zone, and the light blue areas ice-bonded permafrost. All uncolored areas of the plots correspond to unfrozen submarine permafrost and marine sediments.

Figure 2

Pore water profiles of methane, nitrate, sulfate, iron, manganese and relative archaeal abundance in the Buor Khaya core BK2. (A) Methane concentrations are shown in black and corresponding δ13C-values in magenta[13]. (B) Nitrate concentrations are shown in blue, sulfate concentrations in red, manganese concentrations in green, and iron concentrations in light blue. (C) Purple dots on the left show depth location in the core. Relative abundance of archaeal 16S rRNA gene sequences are shown as bar plots. Colors of bars refer to the taxa in the legend below figure. Numbers on the right refer to the exact depth in the core. ANME = ANaerobic MEthanotrophic archaea, DHVEG-1 = Deep Hydrothermal Vent Euryarchaeotal Group 1, DSPEG = Deep Submarine Permafrost Euryarchaeotal Group, MBG-B = Marine Benthic Group B, MBG-D = Marine Benthic Group D, MCG = Miscellaneous Crenarchaeotal Group. The red area represents a sulfate-methane transition zone and the light blue area ice-bonded permafrost. All uncolored areas of the plots correspond to unfrozen submarine permafrost and marine sediments. n.d. – not detected.

Sulfate from sea water penetrated into the ice-bonded permafrost section of core C2 and into the thaw front of core BK2 indicating the presence of a SMTZ (Figs 1B and 2B). Sulfate concentrations of 12 to 24 mM in the unfrozen part of the C2 core decreased with depth and penetrated almost 10 m into the ice-bonded permafrost. Here, sulfate concentrations were still high with up to 12.5 mM. Below 40 mbsf sulfate concentrations were low (<0.4 mM) with the exception of ~52 mbsf (~1.5–2 mM) and the lower permafrost thaw front at ~58 mbsf (2.3 to 5.5 mM). Sulfate in the marine-affected layers of BK2 showed low concentrations (<0.5 mM) in the first 5 mbsf and brackish concentration of 2 to 7 mM between 5 mbsf and the ice-bonded permafrost. Inside the ice-bonded permafrost sulfate concentrations were low (average 0.08 mM) (Fig. 2B). In core C2, other potential terminal electron acceptors (TEA) such as nitrate, manganese, and iron were detected in the permafrost sediment layers (Fig. 1B). The latter two were only detected as total ions with no information on the oxidative state. Nitrate and manganese concentrations displayed high fluctuations in ice-bonded layers showing marine influence[42] in the sediment i.e. between 29 and 43 mbsf. A high concentration of nitrate and manganese was also observed at the lower boundary of the ice-bonded permafrost (58.5 mbsf). In the ice-bonded permafrost part of BK2, manganese, iron and nitrate were close to or even below the detection limit. In the unfrozen part, manganese and iron concentrations were an order of magnitude higher and iron concentrations showed an increase towards the surface (Fig. 2B).

Archaeal community structures of degrading permafrost

The upper layer of core C2, inundated for 2,500 years, was dominated by sequences of typical marine archaeal taxa. They were represented by aerobic ammonia-oxidizing Thaumarchaeota of the genus Nitrosopumilus (~44%), Lokiarchaeota found at hydrothermal vents (~15%), marine benthic group D (6%) and other marine taxa such as marine group II (~2%, Fig. 1C). At 18 mbsf Bathyarchaeota sequences of the MCG-6/pGrfC26 clade dominated (~53%) the archaeal community. This clade includes members able to degrade polymeric carbohydrates in the sulfate reduction zone[44]. In the same depth a shift towards taxa that are commonly found in terrestrial environments such as the ammonia-oxidizing Nitrososphaerales (~25%) was observed. Here we also detected a so far unknown clade, which we designated as Deep Submarine Permafrost Euryarchaeotal Group I (DSPEG I, ~12%). This new clade is closely related to Methanomassiliicoccales of the superclass Diaforarchaea. DSPEG I increased in relative sequence abundance with depth (Figs 1C and 3). DSPEG I and a sister clade (designated DSPEG II) were most dominant in the sediment horizons between 29 and 44.8 mbsf, where they represented between 62 to 84% of all archaeal sequences (Fig. 1C). Both groups showed less than 82% sequence similarity values to the next cultured representative Methanomassiliicoccus luminyensis[45] placing them in a putative new order. Other abundant sequences belong to Methanoflorentaceae (~0.1–5%), ANME-2d (~0.2–6%), Methanosarcina (~0.2–9%), Bathyarchaeota MCG-6 and MCG-8 (~0.8–11%), as well as Nitrososphaerales (~3.7–22%). At 45 mbsf there was a pronounced shift towards Nitrososphaerales sequences (~65%). In the sediment layers between 49–52 mbsf, where most of the electron acceptors showed low concentrations, Bathyarchaeota of the MCG-6/pGrfC26 clade dominated (~41–78%). In these sediment horizons, many methanogen-related sequences appeared, such as Methanoregula (~0.1–6% between 50.1–52.3 mbsf), Methanosarcina (~2–18%), Methanobacterium (~3–10% 50.1–52.3 mbsf) and Methanomassiliicoccaceae (0.7–14% between 48.8–50.6 mbsf). Between 50.1 and 51.8 mbsf Woesearchaeota of the YLA114 clade were present (~3–13%). In core C2 ANME-2d sequences increased in relative abundance from less than 1% at 48.80 mbsf to 18% at 52.3 mbsf and up to 70% at 52.8 mbsf. The high abundance of anaerobic methane-oxidizers coincides around the highest peak of the methane concentration profile in the ice-bonded permafrost (284 µM) at 52 mbsf (Fig. 1). In the same depth, about 1.1% of sequences are related to Methanosarcinaceae/ANME-3 as described by reference[24] (Table S1). A second, much smaller peak of the methane concentration profile (10 µM) at 55.6 mbsf showed 16% ANME-2d sequences. At the same depth Methanosarcina represented the most abundant group (~47% of all archaeal sequences) and DSPEG I and II related sequences were the second most abundant (~21%, Fig. 1C). Towards the lower boundary of the ice-bonded permafrost at 58.7 mbsf, the DSPEG I and II increased again in sequence abundance (~33%). The sediment horizon directly underneath the permafrost contained some sequences related to marine environments, such as Lokiarchaeota (~9%), and ANME-2a/b (2.5%). The deepest sediment sample showed mostly Bathyarchaeota sequences of subgroups MCG-6 and MCG-8 (~77%) and DSPEG I and II, ANME-2d, Methanoflorentaceae and Nitrososphaerales (~6%, 4%, 3% and 2%, respectively) sequences.

Figure 3

Phylogenetic affiliation of submarine permafrost archaeal sequences based on 16S rRNA gene. Taxonomic cluster in red contain sequences of submarine permafrost. Numbers in brackets show the number of OTUs per cluster. The scale bar represents 10 percent sequence divergence. The tree was rooted with the DPANN superphlyum.

The more recently inundated core BK2 contained a mixture of marine and terrestrial archaeal sequences in the depth between 4.7 and 12 mbsf. Marine clades were represented by ANME-2a/2b (~0.5–35%), Methanomicrobia of clade F99a103 (~0.1–16%) that were first discovered at a hydrothermal chimney, ANME-2c sister clade (~16%, 8.1 mbsf), Marine Benthic Group D (~0.1–10%), Lokiarchaeota (~0.7–9), and Nitrosopumilales (~0.3–4%). Terrestrial clades were represented by Methanosarcinales/Methanosarcina (~0.7–15%), DSPEG I and II (~0.1–9%), and Nitrososphaerales (~0.2–5%). The high numbers of ANME-2a/b sequences (Fig. 2C) at depths of ~5 and ~9 mbsf correlated with highest δ13C-methane values in the uppermost 12 m that were characterized by highly variable methane δ13C-values (−67 to −31‰, Fig. 2C). Still, Bathyarchaeota-related sequences of the clades MCG-6 and MCG-8 (~40–71%) dominated throughout the layers above the ice-bonded permafrost. In contrast, they drastically decreased (~2%) in the upper layer of the ice-bonded permafrost. The sediment horizon between 24.0 and 24.7 mbsf was characterized by opposing gradients of sulfate and methane. In the upper part of this layer that is strongly marine influenced ANME-2a/b made up about 26% of all archaeal sequences while in the lower and still partially ice-bonded part of the permafrost thaw front both ANME-2a/b and ANME-2d clades were detected in high abundances (19% and 28%, respectively, Fig. 2C). In the ice-bonded permafrost of core BK2, archaea were only detected in about 45.46 mbsf. Here, Methanosarcina (~37%) and DSPEG I and II (~46%) related sequences dominated. A detailed description of the bacterial community structure and diversity of both cores was described elsewhere[42].

Diversity of ANME-related groups in degrading permafrost

Both cores displayed layers where ANME sequences dominated the archaeal community. The two major groups, ANME-2a/b and ANME-2d showed a remarkably high diversity with 23 and 17 operational taxonomic units (OTUs), respectively (Figs 3 and 4A). ANME-2d separated into 4 different clusters (Fig. 4B), with the majority of all sequences falling into cluster 1. This cluster also contained Candidatus Methanoperedens nitroreducens BLZ1[46], known to either use nitrate[27] or iron[47] as electron acceptor for methane oxidation. Besides the two dominant ANME-groups, we also detected sequences in a sister clade of ANME-2c (designated ANME-2e) with 5 OTUs (Figs 3 and 4 for details see also Figs S1 and S2). Other sequences of this group were retrieved from cold seep sites and mud volcanoes. Moreover, we detected 1 OTU that fell into the ANME-3 cluster (Figs 3 and 4, for details see also Fig. S2). Clone libraries of the SMTZ, using specific ANME-2a primers revealed mainly ANME-2a/b sequences but also ANME-2d as well as one sequence of the ANME-2e cluster (Fig. 4A, for details see also Fig. S1). Corroborative CARD-FISH analysis of the SMTZ detected ANME-2a consortia with the specific ANME-2a-647 probe (Fig. S3), while detection with a specific probe for Desulfosarcina/Desulfococcus showed no signal in any consortium.

Figure 4

Phylogenetic affiliation of ANME sequences based on 16S rRNA gene. (A) The first number in the parenthesis represents the number of OTUs from the Illumina sequencing, while the second number represents the number of sequences of the clone library with ANME-specific primers. Methanopyrus kandleri was used as outgroup. The scale bar represents 10 percent sequence divergence. (B) The ring diagram represents the whole relative abundance of all ANME-2d sequences in both submarine permafrost cores.

Canonical correspondence analysis (CCA) illustrated that ANME-2a/b and ANME-2e clustered with other marine-related archaea and are mainly influenced by salinity, iron and sulfate (Fig. 5). Spearman rank correlations revealed that the abundance of ANME-2a/b significantly correlated with iron concentrations (R = 0.65, p < 0.03). ANME-2d sequences clustered together with sequences of Methanosarcina and Bathyarchaeota of the MCG-6/pGrfC26 clade and fall in the terrestrial methane-influenced quadrant of the plot. However Spearman correlation revealed no significant correlation.

Figure 5

Canonical correspondence analysis (CCA) of environmental factors and archaeal taxa that contributed with more than 1% in any of the depth. Inlet: Environmental factors are plotted as triplot with scaling 3. Samples of core C2 are projected as black dots and samples of BK2 as red dots. The archaeal taxa are shown as blue dots. Percentages at axes represent the ‘Eigenvalues’ that explain the variability for the first two axes.

Functional groups responsible for methane cycling

Quantitative detection of the methyl co-enzyme reductase subunit A (mcrA) as functional marker for methane production (methanogens) and AOM (ANME) showed a decrease with depth in the unfrozen permafrost of both cores (0 to 29 mbsf for C2, and 0 to 24 mbsf for BK2, respectively, Fig. S4A and B). In the ice-bonded permafrost of core C2 we detected an increase in mcrA copy numbers with the highest copy numbers (0.3–1.2 × 105 ± 0.4 - × 103 g−1 wet weight) between ~45 and 52 mbsf. At 53 mbsf total mcrA copy numbers decreased, while 16S rRNA gene sequences at the same depth were almost exclusively affiliated with ANME-2d (Figs 1C and S4A). 454 amplicon analyses of the mcrA sequences showed a relative increase of ANME-2d-related mcrA genes towards the layers with highest ANME-2d-related 16S rRNA gene sequences between 52.3 and 55.6 mbsf (Figs 2C and S5). Quantification of ANME-2d with recently designed mcrA primers in the ice-bonded permafrost[48] showed up to 1.2 ± 0.5 × 105 copy numbers g−1 wet weight at approximately 57 mbsf and an average copy number of 5.1 ± 1.5 × 104 g−1 wet weight over the whole core (Fig. S6).

Analysis of the dissimilatory sulfate reductase subunit B (dsrB) for SRB mainly followed the depth distribution of Desulfosporosinus 16S rRNA sequences in the C2 core, with a few exceptions (29.0, 30.0, 34.4 and 38.5 mbsf) where an increase in the dsrB could not be related to any known sulfate-reducers (Fig. S7A).

McrA copy numbers in BK2 slightly increased at the SMTZ and may be underestimated due to a primer bias against ANME mcrA genes[48], while 16S rRNA gene sequences showed 25–47% ANME-related sequences at the SMTZ (Fig. S4B). The dsrB copy numbers slightly increased at the SMTZ, although SRB represented < 1% of the total bacterial 16S rRNA gene sequences (Fig. S7B). Phylogenetic reconstruction of SRB at the SMTZ indicated that all sequences of the upper part belonged to typical marine clades SEEP-SRB1 and another seep-associated Desulfobacterium anilini-group (red OTUs, Fig. S8). In the lower part of the SMTZ that is dominated by ANME-2d sequences, seep-associated SRB were below detection limit, while only Desulfosporosinus-related sequences (blue OTU, Fig. S8) could be observed.

Lipid biomarker measurements revealed very low concentrations that prevented a precise isotopic analysis to conclude activity of AOM microbial communities. However, calculated branched and isoprenoid tetraether (BIT) indices showed values of terrestrial origin ranging from 0.99 to 1[49]. Moreover calculated methane indices (MI) for isoprenoid glycerol dialkyl glycerol tetraethers (GDGTs) showed values between 0.85 to 0.99 in the unfrozen sediment layer and 1 at the SMTZ[50]. The ratio of archaeal to bacterial ether lipids showed dominant archaeal values in the SMTZ (Fig. S9). For further discussion on lipid biomarker see supplementary material.

Methane oxidation rates in submarine permafrost

The fraction of produced methane that, based on carbon isotopic signature change, was oxidized varied between 72–86% for the C2 core and 79–100% for the BK2 core, respectively. Methane release rates depend on the permafrost degradation rate and the methane concentration in the pore space of the core section. Assuming that our cores are representative for the whole submarine permafrost area with an estimated distribution of 3 million square kilometer[2], modeling for the C2 core resulted in 0.0001–0.0984 Tg C y−1 (average 0.004–0.005 Tg C y−1, highest and lowest oxidation, respectively) and 0.0242–1.1889 Tg C y−1 (average 0.1712–0.2167 Tg C y−1, highest and lowest oxidation, respectively) for the BK2 core.

Discussion

Atmospheric methane concentration has been increasing again for the last decade, but the sources and mechanisms for this increase are not fully understood[51]. One potentially large but highly controversial source of methane is submarine permafrost, which faces drastic changes due to global temperature increase and associated Arctic sea ice reduction[52]. Our study provides multiple lines of evidence that ANME communities are present in submarine permafrost layers where methane is being consumed. We thus propose the microbial mitigation of methane release from thawing deep submarine permafrost on the Siberian Arctic Shelf. We detected both marine and terrestrial ANME clades likely involved in the AOM process at various depths, not only at the permafrost thaw front, but also in still ice-bonded permafrost that undergoes degradation. Our study indicates that AOM occur at temperatures below 0 °C[42] (Table S2).

Pore water methane concentrations in the two cores were in the typical range of deep sediments and soils[20,53]. The permafrost thaw front at the BK2 site, inundated about 540 years ago reflects a deep SMTZ. We identified a steep decline of methane concentration above the SMTZ. This is in accordance with core (IID-13)[54] in close proximity that also showed a large decrease (106 to 103 nM) in methane concentrations at the thaw boundary while δ13C-values of methane from unfrozen layers above (~15 to 4 m) of this core are not reported. Altogether the study of Sapart et al., (2017)[54] finds no evidence for AOM, which might be due to regional influences such as river water, land surface run off and warming by river discharge. Our δ13C-values of methane in the unfrozen part clearly show a large shift that can only be explained by microbial oxidation. The calculated potential rates based on methane release and the fraction of methane oxidized (1.7–2.1 nmol cm−3 d−1 ± 0.9–1.2 nmol cm−3 d−1 13)[13] are typical for margin SMTZ and exceed those of subsurface SMTZ[25].

At the permafrost thaw front of BK2, we detected marine and terrestrial ANME clades that potentially mitigate the methane release into overlying sediment layers. The upper unfrozen part of the SMTZ represents a typical marine sulfate-dependent AOM community of ANME-2a/b (Fig. 2C) affiliated with SEEP-SRB1/seep-associated Desulfobacterium anilini-group (Fig. S7A). These sulfate-dependent AOM communities are of marine origin and are often found in mud vulcanoes, and in methane and hydrocarbon seeps[55]. Visualization of ANME-2a in situ using CARD-FISH further supports their occurrence at the permafrost thaw front of BK2.The molecular detection of associated SRB is, however, not conclusive. This might be due to SEEP-SRB partners that are not targeted by the DSS685 probe[56] or by yet unidentified bacterial partners[57]. The lower part of the SMTZ was characterized by a clear transition to a terrestrial AOM community closely related to Candidatus Methanoperedens nitroreducens (ANME-2d) that might be capable of using alternative TEA such as nitrate and iron[27,47], although no information on the oxidative state of iron is available. A recent genomic study proposed that these organisms have the genetic repertoire for an independent AOM process without a bacterial partner[46]. Besides the analysis of AOM assemblages via a DNA approach we also found high ether lipid MI values in the SMTZ of BK2 that further supports a potential involvement of AOM communities[50]. The high abundances of marine-derived ANME-2a/b sequences in several depths above the SMTZ of BK2 with BIT indices towards an exclusively terrestrial origin show that terrestrial and marine sediments were only slightly mixed in the upper meters of BK2 (Fig. 2). The terrestrial sediments were thus influenced by sea water penetrating down to the permafrost table and thereby transport of marine organisms into deeper layers occurred (Fig. 2). This is consistent with profiles of other environmental parameters (pH, temperature, isotopes and ion concentrations) and the stratification and composition of the bacterial communities[42].

In core C2, the layer between 29 and 43 mbsf is characterized by alternating frozen and partly thawed sediments (Fig. 1) and showed high fluctuations in the concentrations of nitrate and manganese. These fluctuations indicate degradation of thawing permafrost (shaded area, Fig. 1) and an increase in microbial activity resulting in the consumption of labile carbon pools[58]. Active processes in ice-bonded permafrost close to thawing (mean temperature: −1.2 ± 0.2 °C, Table S2) can occur in liquid water films surrounding mineral particles, which form a network in which microbial activity is expected[59]. Microbial activity may thereby be supported through sulfate as additional electron acceptor[60] penetrating into the ice-bonded sediment. Indeed, sulfate concentrations decreased and dsrB copy numbers increased with depth pointing towards active sulfate reduction and effective anaerobic organic matter decomposition[61]. In C2, the SMTZ occurred below the actual permafrost thaw front but inside the ice-bonded permafrost and was characterized by high manganese and nitrate concentration. Nitrate and manganese concentrations at the SMTZ and in the lower thaw boundary were in the higher µM to mM ranges that exceeded those typically reported for subsurface environments[20] by an order of magnitude (http://publications.iodp.org/). So far it is not clear why nitrate as favorable electron acceptor is not quickly consumed under anaerobic conditions. Nevertheless, the high concentrations of nitrate and manganese could promote AOM with alternative TEA[33,35] and shows ongoing degradation of ice-bonded permafrost. Also, the detection of ANME-2a/b in this layer shows their migration into permafrost by downward marine water intrusion[42]. Unlike in core BK2, the ice-bonded layer of C2 was almost free of methane. At the same time ANME-2d occurred almost entirely throughout this layer. The relatively low copy numbers of ANME-2d detected with specific mcrA primers and compared with total cell counts[42] are still in the range of North Sea and River sediment[48]. The highest abundance of ANME-2d coincided with the highest methane concentration, which we consistently detected by several molecular approaches (16S rRNA, mcrA and ANME-2d-specific mcrA). Taken all this together we suggest that ANME-2d members are responsible for AOM in the ice-bonded permafrost before it completely thaws. An alternative explanation for the low methane concentrations and the low abundance of ANMEs in most of the ice-bonded layers are a relic of AOM communities that were active under different environmental conditions in the past. Whether methane was trapped in the permafrost during freezing or it was produced by microbial degradation of organic matter under recent in situ conditions cannot be resolved since radiocarbon analysis would give similar results in both cases.

Even though sulfate reduction might be relevant for organic matter mineralization[61], links to sulfate-dependent AOM were not observed in the core C2. While C2 exhibited high copy numbers of dsrB and of SRB-related sequences, sulfate reducers were almost exclusively linked to Desulfosporosinus that have not been observed in AOM consortia so far. This genus belongs to the phylum Firmicutes and has been found in natural terrestrial environments such as peatlands, aquifer and permafrost[62,63]. Other TEA than sulfate, such as nitrate, iron and manganese, could also be related to AOM, and showed relatively low concentrations at the highest occurrence of ANME-2d sequences at 52 mbsf. This serves further as an indication of methane consumption during the process of AOM as known from physiological studies of ANME-2d enrichments[27,47] but direct evidence is missing. Finally, besides the detected electron acceptors, other TEA such as humic acids could serve in the AOM process. Humic acids were shown to be involved in AOM[36,37] in peatlands where they were detected in high concentrations. Humic acids are produced during organic matter degradation and soil formation[64] and could thus play a role in thawing permafrost, too.

Two clades, which we named DSPEG I and DSPEG II, mainly occurred in submarine permafrost layers that showed relatively high concentrations of iron and manganese in the pore water (Fig. 1B,C). This could point towards an involvement in iron (III) and manganese (IV) reduction within the anaerobic oxidation of organic matter, in addition to sulfate reduction[62]. Since the oxidative state of both metals has not been determined it is not clear whether these metals are used as TEA. Future analysis should focus on oxidative states of metals in these environments to further clarify a potential role in organic matter degradation and AOM. Nevertheless, Spearman rank analysis showed significant correlations between DSPEG I and DSPEG II and manganese concentrations (R = 0.58, p < 0.006 and R = 0.53, p < 0.02, respectively) and negative correlations with methane concentrations (R = −0.57, p < 0.02 and R = −0.64, p < 0.004, respectively) as illustrated in the CCA (Fig. 5). We propose that these two groups reflect indicator taxa (Table S3) for degrading permafrost. This is also supported by high occurrence of DSPEG I and DSPEG II (~19 to 28%, Fig. 1C) at the lower boundary of the ice-bonded permafrost, representing bottom-up permafrost thaw. Here, high concentrations of sulfate, nitrate, and manganese show an upwardly directed thaw process (Fig. 1B). Both groups were also dominant (47%) in the deepest sample of BK2, while degradation and an upward thaw cannot be concluded from our data. Environmental sequences affiliated with the DSPEG groups were mainly retrieved from cold environments[65,66] and pristine aquifers[67], which further support an active role in low temperature habitats.

Taken together our molecular and biogeochemical data from two submarine permafrost cores indicate several microbial assemblages that have the potential to prevent the release of trapped or recently produced methane into the overlying unfrozen sediment following submarine permafrost thaw. Therefore, we challenge the assumption that high methane emissions reported for the Siberian Arctic Shelves originate from degrading submarine permafrost itself[9] and suggest different mechanisms to be responsible, such as diffusion or ebullition through discontinuities in permafrost or the release from gas hydrates[8,68] at a limited spatial scale. Microbial assemblages in deep permafrost environments are usually associated with slow growth rates[69] and low abundances[70], and their activity is difficult to measure. New approaches such as BONCAT-FISH[57] have the potential for more direct detection of active microorganism and the analysis of their genomic potential.

The calculated fraction of methane that was oxidized in the SMTZ of BK2 showed high efficiencies, pointing towards an effective biological methane filter. While methane oxidation within the intact ice-bonded permafrost section of BK2 is unlikely, C2 showed several layers with heavier stable isotopes (≥−45‰) and high fractions (72 to 86%) of methane that were oxidized. Actual methane oxidation rates may be even higher, since methane production of freshly available organic material is not taken into account. Still, on a global scale the estimated submarine permafrost area (~3 million km²)[2], could consume 0.0001 to 1.1889 Tg y−1 of methane from newly degraded permafrost assuming AOM activities similar to those in our cores. This only account for the methane released, since it is difficult to determine the total sediment volume in which AOM activity takes place. Submarine permafrost is thus comparable to AOM in other environments with high methane fluxes such as seep sites (<10 Tg C y−1)[71], while AOM in wetlands (200 Tg C y−1)[32] and marine SMTZ (<50 Tg C y−1)[71] clearly show higher consumption rates. The latter two span areas that are 6 to 40 times larger than those of submarine permafrost.

Our study provides first molecular evidence of microbial communities in thawing submarine permafrost that are likely involved in AOM processes. In addition, many archaeal taxa such as the newly designated DSPEG groups, a large diversity of Bathyarchaeota, and Thaumarchaeota closely related to nitrogen cycling organisms are detected. Their function is unknown and need further investment to understand their contribution in organic matter degradation of permafrost thaw processes.

Materials and Methods

Site description and sampling

Sediment cores were drilled along two transects from terrestrial permafrost to offshore, submarine permafrost in the Siberian Laptev Sea: Mamontov Klyk (2005) and Buor Khaya (2012) (Table S2). The outermost submarine permafrost cores influenced longest by inundation were chosen for chemical and microbial investigation. The Mamontov Klyk study site was located in the western Laptev Sea (Fig. S10). The core drilled 11.5 km offshore is characterized by three different lithostratigraphical units and contains two ice-bonded permafrost layers between 29.5–30.4 mbsf and 34.3–58.7 mbsf, respectively[12]. The core had a total length of 71 m and contained sandy loam sediment with on average 0.38% organic carbon and an average C/N ratio of 14.

The study site of Buor Khaya is located in the central Laptev Sea on the western part of the Buor Khaya peninsula (Fig. S10). The submarine core contained again three lithological units with an ice-bonded permafrost layer between 24.7–47.6 mbsf[13]. The core had a total length of 47.7 m. The retrieved sediment material consisted of fine sand and contained about 1% organic carbon with a C/N ratio of 14.

The drilling was performed with a mobile drilling rig (URB-2A-2/ URB-4T) and is described elsewhere[12,13]. Cores were kept frozen at −20 °C and transported to the laboratories in Germany under frozen conditions. The frozen cores were split along the vertical axis under aseptic conditions. One half was used as an archive, whereas the other half was split into quarters for microbiological and for geochemical, sedimentological and micropalaeontological analyses.

Geochemical analysis

The cores were sectioned at different intervals. The core of the Mamontov Klyk site (C2) was separated into 118 horizons that were used for pore water analyses in 2007, 2 years after the drilling took place. The core was constantly stored at −20 °C between drilling and subsampling. The Buor Khaya site (BK2) was divided into 80 horizons. After thawing of subsamples in 2013, one year after the drilling, pore water was collected using Rhizones© with an effective pore diameter of 0.1 µm. The concentrations of sulfate and nitrate were determined via a KOH eluent and a latex particle separation column on a Dionex DX-320 ion chromatograph. Total manganese and iron ions were measured by a Perkin-Elmer ‘Inductively Coupled Plasma Optical Emission Spectrometry’ (ICP-OES) Optima 3000 XL.

For methane measurements, 3 g of frozen material were retrieved with ice screws and immediately immersed in 20 ml serum vials containing a saturated NaCl solution (315 g l−1). Serum vials were sealed with butyl-rubber stoppers and a crimp seal. Headspace gas was measured in triplicates with different setups. In brief, gas was analyzed with an Agilent 6890[19] or 7890A[13] gas chromatograph equipped with a flame ionization detector and with a carbon plot capillary column or HP-Plot Q (Porapak-Q) column. The temperature of the oven, injector and detector were 40, 120, and 160 °C, respectively. In both cases, helium was used as carrier gas. The amount of gas in the vials was calculated from headspace concentrations, gas pressure and solubility, and the volume of liquid in the bottles. Methane concentrations are reported relative to sediment pore water volume, regardless of whether present as ice or water, based on calculated total sediment water content. The δ13C-CH4 was determined with an isotope ratio mass spectrometer (Finnigan Delta Plus) equipped with a PreCon and a GC/C III interface (Thermo Fisher Scientific). The precision of replicated measurements was better than 0.5‰ VPDB. Methane δ13C-signatures were linked to the VPDB scale using internal (−43.8‰ VPDB) and external (RM8561; −73.27‰ VPDB) standards measured at least every 10 analyses.

DNA extraction, 16S rRNA Illumina sequencing and analysis

In 2014, genomic DNA of 4.7–13 g homogenized sediment from different depth (C2: 19 samples and BK2: 10 samples, Table S2) was extracted after the protocol of Zhou et al., 1996[72]. DNA concentrations were quantified with Nanophotometer® P360 (Implen GmbH) and Qubit® 2.0 Flurometer (Thermo Fisher Scientific).

The 16S rRNA gene for bacteria was amplified with the primer combination S-D-Bact-0341-b-S-17 and S-D-Bact-0785-a-A-21. The 16S rRNA gene for archaea was amplified in a nested approach with the primer combination S-D-Arch-0020-a-S-19 and S-D-Arch-0958-a-A-19 in the first PCR for 40 cycles and S-D-Arch-0349-a-S-17 and S-D-Arch-0786-a-A-20 in the second PCR for 35 cycles, respectively. The primers were labelled with different combinations of barcodes that are listed together with primer sequences in Tables S4 and S5. The PCR mix contained 1 × PCR buffer (Tris•Cl, KCl, (NH4)2SO4, 15 mM MgCl2; pH 8.7) (QIAGEN), 0.5 µM of each primer (Biomers), 0.2 mM of each deoxynucleoside (Thermo Fisher Scientific), and 0.025 U µl−1 hot start polymerase (QIAGEN). The thermocycler conditions were 95 °C for 5 min (denaturation), followed by 40 cycles of 95 °C for 1 min (denaturation), 56 °C for 45 sec (annealing) and 72 °C for 1 minute and 30 sec (elongation), concluded with a final elongation step at 72 °C for 10 min. PCR products were purified with a Hi Yield® Gel/PCR DNA fragment extraction kit (Süd-Laborbedarf). We performed technical duplicates of one DNA sample and PCR products of 3 individual runs per sample were combined. PCR products of different samples were pooled in equimolar concentrations and reduced to a final volume of 10 µl with a concentration of 200 ng µl−1 in a vacuum centrifuge Concentrator plus (Eppendorf).

Illumina sequencing has been performed by GATC Biotech AG using 300 bp paired-end mode. For better performance due to different sequencing lengths, we used 15% PhiX control v3 library.

The quality and taxonomic classification of the Illumina sequences were analyzed with a customized QIIME pipeline[73]. For details see Supplementary material and methods.

Construction of ANME-specific clone libraries

DNA of permafrost thaw horizons with methane peaks were investigated to amplify ANME communities with specific primers. Therefore we used an ANME-specific probe as reverse primer and a universal archaeal primer as forward primer (Table S4). PCR mixes (25 µl) contained 1 × PCR buffer, 0.2 mM dNTP’s, 2 mM MgCl2, 0.08 mg ml−1 bovine serum albumin, and 0.02 units HotStart Taq Plus Polymerase (QIAGEN) and were performed under the following conditions: initial denaturation 94 °C for 10 min, 30 cycles of denaturation 94 °C for 30 sec, annealing 59 °C for 1 min, extension 72 °C for 3 min, and a final elongation at 72 °C for 10 min. PCR amplicons were purified with the HiYield® Gel/PCR DNA Extraction Kit (Süd-Laborbedarf). Purified PCR products were cloned with the TOPO® TA Cloning® Kit (Invitrogen). Positive clones were directly sequenced by Sanger sequencing (GATC Biotech).

Phylogenetic reconstruction

For phylogenetic analyses of 16S rRNA gene of archaeal and bacterial sequences, the ARB software package was used[74]. After manual refinement of representative sequences of OTUs achieved by the Illumina sequencing and partial sequences of clone libraries against the alignment of the SILVA 16S rRNA gene SSU reference database release 115[75], phylogenetic trees were calculated. Phylogenetic reconstructions were based on the maximum-likelihood algorithm (PHYLIP-ML, 100 bootstraps) implemented in ARB with reference sequences (>1400 bp for bacteria and >850 bp for archaea) and the implement archaeal or bacterial position variability filters. Our own partial sequences were added to the tree using the maximum parsimony algorithm without allowing changes in tree topology.

454 pyrosequencing of functional genes and analysis

The mcrA fragment was amplified using the primer set mlas and mcrA-rev (Table S4). In a touchdown PCR denaturation, annealing and elongation time was set at 1 min. The PCR conditions were as follows: initial denaturation at 95 °C for 3 min, 15 cycles with a stepwise temperature decrement from 65 °C to 50 °C, followed by 15 cycles with an annealing temperature of 55 °C and a final elongation at 72 °C for 10 min. Tagging of amplicons with unique multiple identifier (MID) tags (Table S6) for 454 sequencing was conducted in a second PCR using amplicons from the touchdown PCR as template and 15 cycles with a constant annealing temperature of 55 °C. PCR reactions were performed in several separate reactions and pooled until reaching at least 150 ng final product. We used the same chemicals as for 16S rRNA gene amplifications. PCR products were purified with a Hi Yield® Gel/PCR DNA fragment extraction kit (Süd-Laborbedarf). Amplicons were quantified with the Qubit 2.0 system using the ds DNA HS assay kit (Invitrogen), mixed in equimolar amounts and sequenced from both directions by Eurofins Genomics using Roche/454 GS FLX++ technology.

454 mcrA sequences were analyzed with the mothur software[76] by a customized standard operating procedure. For details see Supplementary material and methods.

Quantification of functional genes

Quantitative PCR was performed using the CFX Connect™ Real-Time PCR Detection System (Bio-Rad Laboratories). Each reaction contained iTaq™ Universal SYBR® Green Supermix (12.5 µl per reaction of 2× concentrate, Bio-Rad Laboratories), PCR primers (0.5 µl containing 20 µM each), sterile water (6.5 µl) and DNA template (5 µl) added to a final volume of 25 µl. Primers targeting the functional genes mcrA and dsrB were used (Table S4). The PCR reactions comprised an initial denaturation (5 min at 95 °C), followed by 40 cycles of 5 sec at 95 °C, 30 s at the specific annealing temperature (see Table S4), 10 sec at 72 °C and a plate read step at 80 °C for 3 s. Melt curve analysis from 65 to 95 °C with a 0.5 °C temperature increment per cycle (5 sec) was conducted at the end of each run to check for non-specific amplification of DNA. The qPCR assays were calibrated using known amounts of PCR amplified and cloned gene fragments from corresponding taxa (standards of pure cultures) in the range of 106–101 copies µl−1. Prior to qPCR analysis, DNA templates were diluted 5- to 100-fold and triplicates were analyzed for each sample. The PCR efficiency based on the standard curve was calculated using the BioRad CFX Manager software and varied between 88 and 100%, depending on the standard. All cycle data were collected using the single threshold cq determination mode.

Quantification of ANME-2d

The DNA was diluted 1:100 to prevent inhibition of amplification by environmental compounds, e.g. humic acids, and to keep the cq value of the samples within detectable limits. The abundance of ANME-2d archaea was quantified by qPCR with specific mcrA primers (Table S4). The PCR mix consisted of 10 µl KAPA HiFi SYBR green mix (KAPA Biosystems), primers (0.04 µl containing 100 mM each), MgCl2 solution (0.5 µl of 50 mM), bovine serum albumin (BSA) (0.3 µl), sterile water (9.02 µl) and 1 µl DNA template. qPCR amplifications were performed by 10 min 98 °C initialization; 40 cycles of 5 sec denaturation at 95 °C, 30 sec of annealing at 60 °C, 1 min of elongation at 72 °C, 2 sec at 82 °C for fluorescence detection; and finally a melting curve from 55 °C to 95 °C with a 0.5 °C increment every 5 sec. For quantification, a tenfold dilution series of mcrA156F/mcrA345R product cloned into a pGEM-T easy plasmid of a known copy number was used as a standard[48]. The PCR product specificity was checked by melt analysis and compared to the standard with a melting temperature of 82 °C.

Microbial lipid biomarker extraction

Different fractions of lipids were extracted from sediment of the SMTZ. For a detailed method description see supplementary material.

Catalyzed reporter deposition fluorescence hybridization (CARD-FISH)

For CARD-FISH, sediments were fixed for 1 h at room temperature with 4% formaldehyde (Fluka). For amplification, a fluorescine-labelled tyramide was used. The protocol for CARD-FISH followed the description published earlier[77]. For details see supplementary material and methods.

Statistical analysis of archaeal community

Statistical analysis was performed with the PAST 3.14 software[78]. Canonical correlation analysis (CCA) was performed with 7 environmental parameters (chloride, pH, conductivity, methane, sulfate, manganese, iron, and nitrate) as explanatory variables and relative abundance of archaeal genera to families. Explanatory variables were standardized by log10 transformation prior to computation, with the exception of pH and conductivity. Chloride, pH and conductivity were used to account for the marine influence of seawater penetration into the permafrost. The significance of canonical axes was tested via permutation computed for N = 999. Spearman correlation analysis was performed with the implemented tools in PAST 3.14.

Modeling of methane oxidation rates in submarine permafrost

We used the highest change in stable isotope signature of dissolved methane in both cores to calculate the fraction of produce CH4 that got oxidized (f). For that we used the formula[79]:δo and δp are the δ13C values (in ‰) of CH4 in the oxidized submarine permafrost layers (>−37‰) and of CH4 produced or trapped in subjacent submarine permafrost layers (<−52‰), respectively. αox and αtrans are the fractionation factors for anaerobic oxidation of CH4 and CH4 transport, respectively.

We assumed that CH4 is transported by diffusion and used a fractionation factor for soil of αtrans = 1.001[80]. Fractionation factors (αox) for sulfate-dependent (S-)AOM in marine enrichments range between 1.009–1.039[81]. Therefore, we used the highest and lowest value for calculations. Since the submarine permafrost is a former freshwater system we also used fractionation factors for S-AOM (αox = 1.030)[82], iron-dependent (Fe-) AOM (αox = 1.031)[83], nitrate-dependent (N-) AOM (αox = 1.032)[84] and extracellular electron transfer (EEL-) AOM (αox = 1.0174)[85], latter likely reflecting humic acids as potential electron acceptor. Calculated fractions of oxidized methane for the low marine S-AOM and the EEL-AOM revealed values over 1 and seems to not represent suitable AOM fractionation factors and were not considered further.

To compare the AOM activity of submarine permafrost to that of other habitats, we took the lowest, highest and average methane concentration of ice-bonded permafrost for both cores and calculated the methane release rate due to permafrost degradation rate for the cores. The permafrost degradation rate is calculated by the mean coastal erosion rate and the permafrost table depth[12,13] and gave values of 0.6 cm per year for C2, and 5.3 cm per year for BK2[86], respectively. Assuming the AOM activity of our cores to be representative for submarine permafrost, we further extrapolated methane concentration to the whole area of submarine permafrost (3 million square kilometers)[2]. As last step we used the highest and lowest fraction of methane that got oxidized to calculate specific rates for both cores.

Data deposition

Sequences of the submarine permafrost metagenome have been deposited at the NCBI Sequence Read Archive (SRA) with the Project number BioProject ID# PRJNA352907, under the accession numbers SRR5183846-SRR5183871 for archaeal 16S rRNA gene sequences, SRR5184420-SRR5184446 for bacterial 16S rRNA gene sequences, and SRR5019811-SRR5019818 for mcrA sequences. ANME-specific, partial 16S rRNA gene sequences of clone libraries are available in Genbank, EMBL and DDBJ under the accession numbers KY613956-KY613991.