The Impact of Methane on Microbial Communities at Marine Arctic Gas Hydrate Bearing Sediment

Carrier, Vincent; Svenning, Mette M.; Gründger, Friederike; Niemann, Helge; Dessandier, Pierre-Antoine; Panieri, Giuliana; Kalenitchenko, Dimitri

doi:10.3389/fmicb.2020.01932

ORIGINAL RESEARCH article

Front. Microbiol., 24 September 2020

Sec. Aquatic Microbiology

Volume 11 - 2020 | https://doi.org/10.3389/fmicb.2020.01932

The Impact of Methane on Microbial Communities at Marine Arctic Gas Hydrate Bearing Sediment

$\r\nVincent Carrier,*$ Vincent Carrier^1,2*

Mette M. Svenning^1,2*

Friederike Gründger³

Helge Niemann^2,4,5

Pierre-Antoine Dessandier²

Giuliana Panieri²

Dimitri Kalenitchenko²

¹Department of Arctic and Marine Biology, The Arctic University of Norway, Tromsø, Norway
²Centre for Arctic Gas Hydrate, Environment and Climate, The Arctic University of Norway, Tromsø, Norway
³Department of Bioscience, Arctic Research Centre, Aarhus University, Aarhus, Denmark
⁴Department of Marine Microbiology and Biogeochemistry, Royal Netherlands Institute for Sea Research, and Utrecht University, Den Burg, Netherlands
⁵Department of Earth Sciences, Faculty of Geosciences, Utrecht University, Utrecht, Netherlands

Cold seeps are characterized by high biomass, which is supported by the microbial oxidation of the available methane by capable microorganisms. The carbon is subsequently transferred to higher trophic levels. South of Svalbard, five geological mounds shaped by the formation of methane gas hydrates, have been recently located. Methane gas seeping activity has been observed on four of them, and flares were primarily concentrated at their summits. At three of these mounds, and along a distance gradient from their summit to their outskirt, we investigated the eukaryotic and prokaryotic biodiversity linked to 16S and 18S rDNA. Here we show that local methane seepage and other environmental conditions did affect the microbial community structure and composition. We could not demonstrate a community gradient from the summit to the edge of the mounds. Instead, a similar community structure in any methane-rich sediments could be retrieved at any location on these mounds. The oxidation of methane was largely driven by anaerobic methanotrophic Archaea-1 (ANME-1) and the communities also hosted high relative abundances of sulfate reducing bacterial groups although none demonstrated a clear co-occurrence with the predominance of ANME-1. Additional common taxa were observed and their abundances were likely benefiting from the end products of methane oxidation. Among these were sulfide-oxidizing Campilobacterota, organic matter degraders, such as Bathyarchaeota, Woesearchaeota, or thermoplasmatales marine benthic group D, and heterotrophic ciliates and Cercozoa.

Introduction

Cold seep microbial communities thrive where geofluids, characterized by high concentrations of hydrocarbons, in particular methane (CH₄), provide a primary energy source for these organisms (Boetius et al., 2000; Orphan et al., 2002; Niemann et al., 2013). These geofluids and/or free gas migrate upward through faults, cracks, and sediment pores that provide a transport vector from sub-seafloor reservoirs to the seafloor. The origin of methane can be either from the geological cracking of organic matter at high temperature or from biologically mediated decomposition of organic matter (Schoell, 1988; Joye et al., 2010). Under certain thermobaric conditions, CH₄ forms gas hydrates, i.e., an ice-like lattice comprising molecules of CH₄ trapped in crystalline cages of water molecules. The formation or the dissociation of gas hydrates can modify the seafloor morphology, and subsequently can lead to the genesis of pockmarks, craters, and gas domes (Vogt et al., 1994; Hovland and Svensen, 2006; Koch et al., 2015; Portnov et al., 2016; Serov et al., 2017; Waage et al., 2019).

The CH₄ present in the fluid can be oxidized aerobically or anaerobically (Krüger et al., 2005; James et al., 2016). In aerobic environments, the oxidation of CH₄ is driven by methane oxidizing bacteria that utilize oxygen as an electron acceptor. Most of them are associated with the alpha and gamma proteobacteria, but also with Verrucomicrobia or Crenothrix (Hanson and Hanson, 1996; Knief, 2015). Nevertheless, microbial activity at the cold seep seafloor rapidly depletes the available oxygen in marine sediments and limits its penetration depth to a small surface layer, usually of a few millimeters thickness at most (Niemann et al., 2006, 2009; Reeburgh, 2007; Boetius and Wenzhöfer, 2013). In the absence of oxygen, methane is oxidized anaerobically through a process that has been termed the anaerobic oxidation of methane (AOM; Reeburgh, 2007). Anaerobic oxidation of methane is driven by anaerobic methanotrophic Archaea (ANME) and so far, three main ANME clades of phylogenetically distinct groups were detected: ANME-2 and ANME-3 are placed within the methanosarcinales, while ANME-1 forms a distinct group within the Halobacterota (Knittel and Boetius, 2009; Quast et al., 2012; Yilmaz et al., 2014). The phylogenetic dissimilarity of these ANME groups suggests different levels of tolerance to various environmental parameters. Previous study results suggested that ANME-2 might be more sensitive than ANME-1 to high concentrations of sulfide and low concentrations of sulfate (Timmers et al., 2015; Bhattarai et al., 2018). The ANME-2 group would then often be limited to the layers at the sulfate-methane transition zone (SMTZ) and ANME-1 would dominate in more sulfidic sediments, at deeper layers (Knittel et al., 2005; Roalkvam et al., 2012). Nevertheless, ANME-2 groups were also retrieved in sulfide-rich sediments (for example at the Hydrate Ridge; Knittel et al., 2003), insinuating the impact of other factors on the observed stratification of ANME groups. Additional environmental conditions that were suggested to select for differential ANME groups include temperature (Nauhaus et al., 2005; Rossel et al., 2011), salinity (Maignien et al., 2013), or CH₄ flux rates (Girguis et al., 2005; Yanagawa et al., 2011; Marlow et al., 2014).

Most ANME use sulfate, but some were also found to use iron, manganese, and nitrite/nitrate as electron acceptors (Beal et al., 2009; Ettwig et al., 2010, 2016; Hu et al., 2014). Reduction of sulfate at the SMTZ generally requires sulfate reducing bacteria (SRB) and a syntrophic consortium with ANME that are commonly found as AOM drivers (Boetius et al., 2000; Wegener et al., 2015). However, in the last decade, community studies of methanotrophs have shown evidence of free-living ANME cells particularly assigned to the ANME-1 group, but also to the ANME-2 group, that might perform sulfate reduction alone (Orphan et al., 2002; Knittel et al., 2005; Roalkvam et al., 2011; Milucka et al., 2012; Stokke et al., 2012; Gründger et al., 2019).

The AOM coupled with sulfate reduction generates HS^– which can subsequently be oxidized by sulfide-oxidizing bacteria, such as the bacterial mat forming Beggiatoa or Campilobacterota species. Some chemoautotrophs can also be present as intracellular and extracellular symbionts within larger fauna, but also in the eukaryotic euglenozoans and ciliates (Buck et al., 2000; Rinke et al., 2006). Additionally, a higher bacterial and archaeal biomass becomes a trophic basis for grazing megafauna or microbial eukaryotes, including diverse bacterivore ciliates, Cercozoa, and stramenopiles (Werne et al., 2002; Takishita et al., 2007, 2010; Niemann et al., 2013). Potentially parasitic or pathogenic eukaryotes, such as Apicomplexa, Ichthyosporea, and fungi, are also likely to benefit from the denser faunal community (Atkins et al., 2002; Takishita et al., 2006).

In the Arctic, gas hydrate bearing domes were observed 50 km south of Svalbard in Storfjordrenna, at ∼390 m below sea level (Serov et al., 2017). They are referred to as pingos, after similar terrestrial features observed in glacial valleys (Mackay, 1998), although they differ by their formation (i.e., gas hydrates instead of regular water ice; Serov et al., 2017). At the water depth of the gas hydrate pingos (GHP; ∼390 m, ∼0.5–2.5°C bottom water T°C), the gas hydrates remain within the gas hydrate stability zone (GHSZ), but are close to its upper limit and are sensitive to even small changes of temperature and pressure (Hong et al., 2018). Hydroacoustic observations have revealed acoustic flares originating from methane gas bubbles in the water column. These were primarily located at the summit on four of the five Storfjordrenna pingos. The dating of methane derived authigenic carbonates suggested that CH₄ seepage has been active for several thousand years (Serov et al., 2017). Visual observations have revealed a higher biomass in sediments of the pingos compared to the surrounding seafloor (Åström et al., 2018). This can be explained by the presence of a carbonate crust induced by AOM, which offers a hard substrate for the attachment of benthic organisms, such as sponges and anemones (Niemann et al., 2005; Cordes et al., 2010; Vaughn Barrie et al., 2011).

Past investigations at the Storfjordrenna pingos have primarily addressed the geochemical conditions (Serov et al., 2017; Hong et al., 2018) or the biodiversity of larger fauna (Sen et al., 2018; Åström et al., 2018), but the microbial community structure remains mostly unknown [with the exception of a biofilm retrieved within deeper sediments at the pingos (Gründger et al., 2019)]. At a circular seep further south, the Haakon Mosby Mud Volcano (HMMV), the composition of the bacterial and archaeal communities varied between concentric zones around the apex of the edifice, i.e., along a methane flux/concentration gradient (Niemann et al., 2006; Lösekann et al., 2007). In Storfjordrenna, the gas flares at the summit of the structures could suggest a similar concentric arrangement of microbial habitats. However, these pingos contrast with HMMV by presenting a multitude of small geological fractures and gas hydrates chaotically distributed around the structures through which methane migrates to the seafloor surface (Hong et al., 2018; Waage et al., 2019).

Our study aimed at determining spatial variations in the microbial community structure along a gradient from the apex to the edge of three pingos. We addressed key environmental factors that are influencing the prokaryotic and eukaryotic community structures and their spatial distribution. Finally, we identified key taxa characteristics for these Arctic CH₄-rich environments, demonstrating the uniqueness of this ecosystem.

Materials and Methods

Study Site

The sampling site was located in the Arctic Ocean at the mouth of Storfjordrenna, 50 km south of Svalbard at approximatively 390 m water depth (Serov et al., 2017). A group of 5 GHPs (∼10 m high, 500 m in width) distributed on the seabed over a 2.5 km² area were recently found¹. Hydroacoustic surveys and real time visually guided observations with a TowCam-Multicore (see footnote 1) System (TC-MC) and Remotely Operated Vehicle (ROV) dives² have revealed acoustic flares of gas bubbles consisting predominantly of CH₄ and being emitted from 4 of the 5 GHPs. One structure (GHP 5) did not show any visible flares or gas hydrate in sediments (Serov et al., 2017). During the sampling campaigns for this study, seep activity at the different sampling sites was assessed by tracking flares through hydro acoustic surveys with a multibeam echosounder (Kongsberg Simrad EM 302) or by visual observations using a TC-MC system configuration (Panieri et al., 2017).

Sampling Procedure

Field campaigns were conducted with RV Helmer Hanssen and sediment cores at the GHP 3 and at GHPs 1 and 5 were taken in June 2016(see footnote 2) and June 2017³, respectively. Cores were taken along a spatial gradient from the apex of the geological feature to its edge. Core IDs MC_900 (apex), MC_902, MC_918, and MC_919 (edge) were taken at GHP 1. Core IDs MC_1061 (apex), MC_1062, MC_1063, and MC_1065 (edge) were collected at GHP 3 while core IDs MC_920 (apex), MC_922, and MC_923 (edge) were collected at GHP 5 (Figure 1). A reference core (core ID 898) was retrieved at one kilometer away from the closest GHP. The use of the multicore system KC Denmark DK8000 integrated with a TC-MC with a real time transmission of images (Daniel et al., 2003) allowed for the collection of six 60 cm long real time visually guided cores. The combined TC-MC was used to visually survey and sample sediments from the study site and the sediment recovery varied between 15 and 40 cm. Exceptionally, core ID BC_1029 was taken using a blade core mounted on a Sperre Subfighter 30k ROV to target directly sediments in close vicinity to a CH₄ gas flare at GHP 3 in June 2016.

FIGURE 1

Figure 1. Geological dome structures referred to as Gas Hydrate Pingo (GHP) located at the mouth of Storfjordrenna, 50 km south of Svalbard. The upper panel gives an overview of this area. The lower panels show the three selected GHPs for this study. White dots represent locations of the different sediment cores at GHPs 1, 3, and 5 using a multicore system. Core MC_898 was sampled as reference site and core BC_1029. GHP 3 was taken using a blade core mounted on a ROV to sample in the vicinity of a methane gas flare.

Porewater Geochemistry

Porewater geochemistry was measured for all cores, and data for BC_1029 and MC_1063 were collected from Hong et al. (2020). CH₄ concentrations were measured with a head space technique and gas chromatography (Thermoscientific Trace 1310) equipped with a flame ionization detector (Hoehler et al., 2000; Panieri et al., 2017). For this, we extruded 3 mL of bulk sediments per 2 cm intervals in all cores which were immediately transferred to a 20 mL headspace vial with 7 mL of NaOH solution (1 M) and two glass beads, and instantly capped. Samples were analyzed subsequently to an equilibration period of 24 h and are represented as concentration per sediment volume. Sediment porosity was determined from weight and volume measurements as presented in Boyce et al. (1973).

Dissolved iron (Fe²⁺), alkalinity, total sulfide (ΣHS), sulfate (SO₄^2–), and dissolved inorganic carbon (DIC) were measured from a neighboring core of the multicore system, recovered during the same sampling round. Using rhizon samplers (Seeberg-Elverfeldt et al., 2005), up to 18 mL of porewater was collected at each cm in the upper 10 cm and at intervals of 4–10 cm in the lower part of the core. Alkalinity and Fe²⁺ were measured onboard by titration and by spectrophotometry, respectively (Hong et al., 2017). SO₄^2– was measured onshore using ion chromatography (Hong et al., 2017), while ΣHS was measured using a spectrophotometer at a wavelength of 670 nm (Cline, 1969). A detailed protocol on the measurement of ΣHS can be found in the Supplementary Material of Hong et al. (2020). Due to equipment availability on the two field cruises, ΣHS and DIC concentrations were not measured for all sediment cores while alkalinity and Fe²⁺ concentrations were only measured for a selection of sediment layers (Supplementary Tables 1–5).

DNA Extraction, Sequencing, and Sequences Analyses

Sediment cores were extruded and 2 cm thick layers were transferred in Whirl-Pak^® sterile sampling bags (Nasco, United States) and stored at -80°C. Following the measurements of the different environmental parameters in the laboratory, 55 of these samples were selected for amplicon libraries sequencing. These samples were selected at regular depths (surface, ∼5, ∼10, and ∼15 cm) and at clear geochemical interfaces as detected by porewater geochemical gradients (e.g., SMTZ). In a cold room (4°C), sediments were manually ground in liquid nitrogen using a sterilized mortar. The DNA was extracted using the DNeasy PowerSoil Kit (Qiagen, Germany). The manufacturer protocol was followed, except that the samples were placed in G2 DNA/RNA Enhancer beads tubes (Ampliqon, Denmark) for physical lysis (Jacobsen et al., 2018) instead of the kit lysis tubes. Once the DNA samples were quality checked using electrophoresis gels, the DNA concentrations were measured using a NanoDrop^TM 2000 spectrophotometer (Thermo Fisher Scientific, United States) and normalized before being sent to the IMGM Laboratories GmbH for library preparation and amplicon sequencing. For each sample, eukarya were amplified using 18S rDNA degenerate primers to target the V4 region, and bacteria and archaea were amplified using 16S rDNA degenerate primers to target the V3-V4 region (Supplementary Table 6). Library generation was conducted in accordance with the company’s protocols before being sequenced using a Miseq System (Illumina inc., United States). Paired-end nucleotide reads were deposited at Sequence Read Archive Genebank⁴ as BioProject accession number PRJNA593930.

Paired-end reads were meticulously processed and the workflow was derived from the USEARCH suggested protocol⁵. Pairs were merged before being length trimmed and quality filtered with USEARCH v10.0.240. Thereafter, operational taxonomic units (OTUs) were constructed using the UPARSE-OTU greedy algorithm at 97% pairwise sequence identity. Singleton OTUs were removed and taxonomy was assigned using the method Wang implemented in Mothur to the SILVA database release 138 (Edgar, 2010; Quast et al., 2012; Yilmaz et al., 2014). Sequences that were not classified to their domain were discarded prior to further statistical analyses. Finally, sequences from multicellular organisms are likely detected within the 18S rDNA libraries and therefore OTUs that were assigned to Metazoan groups and unclassified eukaryotes were discarded to focus only on the microbial community.

Statistical Analyses

Archaeal, bacterial, and eukaryotic libraries were rarefied at 8900, 4700, and 1300 sequences, respectively, corresponding to the lowest number of sequences in one sample. Preliminary analyses of the libraries demonstrated a large fraction of OTUs that contained just a few sequences in a sample, especially for the bacterial communities (Figure 6). In this study, we aimed to determine the distribution patterns of key microbes. The inclusion of a large fraction of rarer taxa in the diversity analyses, despite sharp gradients in the dominating OTUs, prevented the visualization of these gradients of community changes. Therefore, only OTUs containing at least 1% of the overall sequences of one sample were kept for further statistical analyses.

For the bacterial and archaeal communities, beta diversity, measuring changes in the composition of communities between different samples, was calculated on the relative abundance of the selected abundant OTUs using the Bray–Curtis dissimilarity index implemented in the Vegan v2.5-5 package on R (Oksanen et al., 2019). Clusters of sediment samples sharing similar OTUs abundance and composition for both domains of life were formed at a dissimilarity index of ca. 0.5–0.6. For each cluster, the relative abundance of each OTU was averaged and used to build a doughnut diagram with the R package ggplot2 v3.2.1. Thereafter, distance-based redundancy analyses (dbRDA) were performed to reveal whether the environmental parameters measured had an impact on the observed community dissimilarity between the different sediment cores. A dissimilarity matrix was built using the Bray–Curtis dissimilarity index. As the environmental parameters differed between the GHPs, and the fact that missing values can affect the outcome of the analyses, the dbRDA were performed and presented for each GHP separately. Environmental parameters were logarithmically transformed and standardized through Z scoring (Legendre and Legendre, 1998). The significance of the resulting axis from the dbRDA was evaluated through permutation tests (n = 999). Both functions for dbRDA and permutations tests are implemented in the Vegan v2.5-5 package on R (Oksanen et al., 2019).

For the eukaryotic libraries, biodiversity analyses were likely affected by the removal of sequences assigned to Metazoa, as in some samples they could represent on average 40% of the sequences. Furthermore, a large fraction of the community structure at the GHPs site was dominated by reads assigned to photosynthetic eukaryotes that might have originated from the sedimentation of phytoplankton cells, undermining any subsequent attempts at describing the structure of the eukaryotic communities thriving at the GHPs and evaluating the impact of environmental factors on the biodiversity (Rey and Rune Skjoldal, 1987). Therefore, a different approach was used for the eukaryotic libraries and we emphasized instead on the contrast of the abundant OTUs composition between the reference site and CH₄-rich sediments. To do so, once sequences assigned to Metazoa or unclassified eukaryotes were removed and eukaryotic libraries were rarefied, OTUs that were abundant at the reference site were subtracted and presented separately. We hypothesized that the remaining abundant OTUs would be indicators of taxonomic groups influenced by local conditions at the GHPs. Analyses on the relative abundances of these taxonomic groups were calculated using the Bray-Curtis dissimilarity index (Oksanen et al., 2019) and clusters of sediment samples were formed at a dissimilarity index of 0.5–0.6.

Benthic Foraminiferal Analyses

We observed that the relative abundances of certain prokaryotic taxonomic groups, including the genus Sulfurimonas, increased in CH₄-rich sediments. To ensure that the changes in relative abundances of these taxonomic groups were caused by the presence of CH₄, we compared results from DNA sequences with an independent proxy for surface CH₄-rich sediments. Agglutinated foraminifera are often observed in Arctic seas (Wollenburg and Mackensen, 1998; Jernas et al., 2018) and are particularly sensitive to cold seeps where they are very rare or even absent (Panieri and Sen Gupta, 2008; Martin et al., 2010; Dessandier et al., 2019). Accordingly, changes in their abundances can be used to assess the impact of CH₄ seepage disturbance on the local biological communities. Foraminiferal samples (0–1 cm sediment depth) from GHP 1 were stored for 14 days at 4°C in a 2 g L^–1 Rose Bengal solution in ethanol 96%, in order to identify the living (Schönfeld et al., 2012), or recently alive individuals (Rose Bengal stained foraminifera; Corliss, 1991). All samples were wet sieved using 63 and 125 μm mesh sieves and dried at 40°C (48 h). We considered “living” individuals as the ones characterized by a pink stain of all chambers in their test, with the exception of the last one. In case of doubt, the test was broken to investigate the staining of the endoplasm (Schönfeld et al., 2012). All benthic foraminiferal specimens from >125 μm size fraction were handpicked, identified, and counted. The density was calculated by dividing the number of agglutinated foraminiferal individuals (Supplementary Table 10) in each core by the surface of the core (5.02 × 10³ m²). The relationship between the density of agglutinated foraminiferal cells and the logarithm of the number of resampled Sulfurimonas sequences was tested using a linear model.

Results

Environmental Characterization and Geochemistry

At the reference site, CH₄ was nearly absent, gas flares were not detected on the echosounder, and CH₄ sediment concentrations did not exceed 4 μM (Figure 2). ΣHS remained undetectable throughout the reference core, while measured concentrations of SO₄^2– slightly decreased from 28 mM at the sediment surface to 26 mM at 11 cm below seafloor (bsf), correlating to seawater concentration of the Barents Sea. The seafloor was muddy and authigenic carbonates were not observed (Supplementary Table 4).

FIGURE 2

Figure 2. Geochemical profiles of the different sediment cores collected at the reference site (MC_898) and GHPs 1, 3, and 5, with the exception of BC_1029 that is presented in Figure 3. Profiles for CH₄, SO₄^2–, and ΣHS at these sediment depths are given. Black bars correspond to the sediment layers from which the DNA was extracted and sequenced.

At GHP 1, gas flares and high CH₄ sediment concentrations were suggestive of high CH₄ seepage activity (Figures 1, 2). Dense patches of chemosynthetic organisms, such as siboglinids, as well as carbonate crusts colonized by anemones and sponges, were scattered across GHP 1 (Supplementary Table 4). Concentrations of CH₄ were low in the sediment surface layer, ranging from 0.61 to 6.73 μM, and increased with depth in cores taken at the GHP 1 apex, reaching a maximum of 169 μM at 37 cmbsf in core MC_900 and 1109 μM at 19 cmbsf in core MC_902 (Figure 2). MC_902 was also characterized by a stronger depletion of SO₄^2– with depth than at the reference site as concentration dropped below 5 mM at 15 cmbsf. With the decrease in SO₄^2–, ΣHS concentrations increased, peaking at 4558 and 2078 μM in MC_900 and MC_902, respectively. MC_918 was collected close to the rim of the GHP, where concentrations of CH₄ and ΣHS increased with depth, but at lower concentrations than at cores taken near the apex of GHP 1. The SO₄^2– concentrations values at MC_918 ranged from 27.8 mM at the surface to 25.9 mM at 19 cmbsf. MC_919 was taken outside the GHP, but close to its edge. Here, environmental parameters became more similar to the reference site. Low concentrations of CH₄ (yet still slightly higher than at the reference site) were detectable and SO₄^2– concentrations were only slightly lower than at the reference site and remained above 26 mM within this core.

At GHP 3, BC_1029 had the highest CH₄ concentrations of all sites, reaching up to 12.8 mM at 12 cmbsf (Figure 3). This core was taken in the vicinity of a CH₄ gas flare (Figure 1). The four other cores from GHP 3 had lower CH₄ concentrations than BC_1029 (<15 μM). Still, the cores MC_1061 and MC_1062, located close to the GHP 3 apex, had higher CH₄ concentrations than cores MC_1063 and MC_1065, collected near the edge and outside GHP 3, respectively. SO₄^2– maximum concentrations in the surface sediment layers were in the range of 27-28 mM for all cores, but the SO₄^2– level decreased to 12.21 mM at 12 cmbsf and at 23.83 mM at 14 cmbsf for BC_1029 and MC_1061, respectively. Within other cores taken at GHP 3, the decreasing concentrations of SO₄^2– showed a similar pattern to the reference site. Fe²⁺ concentrations were only measured in two cores (BC_1029 and MC_1063) and showed a sharp decrease at the sediment surface in core BC_1029, but remained high in core MC_1063, where it was depleted only at 20 cmbsf (Supplementary Table 2).

FIGURE 3

Figure 3. Geochemical profiles of the sediment core BC_1029 collected at GHP 3 near a CH₄ gas flare. Profiles for CH₄ and SO₄^2– at these sediment depths are shown. Black bars correspond to the sediment layers from which the DNA was extracted and sequenced.

At GHP 5, similar CH₄ concentrations in the upper sediment layer were measured in core MC_920 and at the apex of GHP 3 (Figure 2). However, gas flares were not visible on the echosounder at the apex of GHP 5. In addition, a CH₄ concentration of ∼18 μM was measured at 19 cmbsf in core MC_922, occurring concomitantly with an increasing concentration of ΣHS. The seafloor was covered with hard surfaces, mostly ice raft debris, and colonized by anemones and sponges (Supplementary Table 5). Complementary information on visual observations at the sampling sites and on concentrations of Fe²⁺, alkalinity, and DIC are available as Supplementary Information (Supplementary Tables 4, 5).

Taxonomy and Abundant OTUs

Once pair-ends reads were quality filtered, 8 129, 36 301, and 8184 OTUs were successfully assigned to the archaeal, bacterial, and eukaryotic domains, respectively (Supplementary Table 7). After rarefaction, within the archaeal OTUs, 87 were found to be abundant in at least one of the sediment layers collected from the reference site or the GHPs. Among the bacteria and eukaryotes, 107 and 140 abundant OTUs were identified, respectively.

Composition Similarities of the Microbial Communities

Based on the beta-diversity dissimilarity analyses in the GHP sediments (Figure 4), six different community clusters designated A1, A2, A3, A4, A5, and A6 were identified for the archaeal domain. Cluster A1 included nearly all surface sediment samples and was dominated by the Crenarchaeota Nitrososphaera and the Nanoarchaeota Woeserchaeales, with 61.4 and 10.0% of the total archaeal community, respectively. Sediment layers associated to the cluster A2 were from different depths, although most were collected between 4 and 10 cmbsf. The cluster A2 was characterized by a stronger dominance of Woeserchaeales (21.86%), Bathyarchaeia (12.14%), Nitrososphaera (11.73%), the marine benthic group D (MBG-D) within the Thermoplasmatota (4.49%), and Asgardarchaeota (3.35%). In addition, 2.1% of the sequences were associated to an unclassified archaeal OTU. The community of the cluster A3 was driven by the MBG-D (20.1%), the Bathyarchaiea (14.5%), and the Woesearchaeales (12.7%). The Asgardaeota groups of Heimdallarchaeia (7.3%) and Lokiarchaeia (1.1%), and the Halobacterota ANME-1 group (1.3%) were also predominant. The cluster A4 had a similar community composition to the cluster A3 and was dominated by the MBG-D (18.5%), the Wosearchaeales (10.8%), and ANME-1 (8.8%). Bathyarchaeia (7.7%) and Asgardarchaeota (5.8%) were also major components of the A4. The clusters A5 and A6 differed from the other groups particularly by a higher relative abundance of sequences associated to ANME groups. The cluster A5, representing sediment layers at the gas flare (core BC_1029) was mainly composed of ANME-1 (17.9 %), ANME-2a/-2b (6.6 %), and ANME-2c (6.1%). Other abundant taxonomic groups included Wosearchaeales (13.2 %) and the MBG-D (4.9%), in addition to the Asgardarchaeota Heimdallarchaeia (6.1%) and Lokiarchaeia (5.3%). The ANME communities of A6 was in contrast to A5 by a stronger dominance of ANME-1 (60.7%; Figure 7), in comparison to the ANME-2a-2b (2.9%) and ANME-2c (4.1%; Figure 5), were also abundant in the cluster A6 representatives from the MBG-D (8.3%).

FIGURE 4

Figure 4. The sediment archaeal communities clustered in six different groups, calculated based on the Bray–Curtis dissimilatory index from the composition of abundant archaeal OTUs (A1, A2, A3, A4, A5, and A6). The averaged taxonomic composition of each cluster is illustrated in a doughnut chart with colors indicating taxa listed in the box. “Non-abundant” includes sequences assigned to archaeal OTUs that were not retrieved in abundance in this study. Finally, the heatmap gives standardized values (Svalue) of depth and of the logarithmic concentrations of CH₄ and HS. Non-available data are represented by striped squares.

FIGURE 5

Figure 5. Physico-chemical and biological vertical profiles of the cores MC_900 and MC_902 (GHP 1) where a SMTZ was identified. For each core, the left figure presents the CH₄, SO₄^2–, and HS concentration profiles. The centered and right figures present the relative abundances of archaeal methane oxidizers (ANME groups) and sulfate-reducing bacteria groups, respectively.

FIGURE 6

Figure 6. The sediment bacterial communities clustered in five different groups, calculated based on the Bray–Curtis dissimilatory index from the composition of abundant bacterial OTUs (B1, B2, B3, B4, and B5). The averaged taxonomic composition of each cluster is illustrated in a doughnut chart with colors indicating taxa listed in the box. “Other” relates to sequences that are assigned to OTUs abundant throughout the whole communities, but not within the illustrated cluster. As for the group “non-abundant”, it includes sequences assigned to bacterial OTUs that were not retrieved in abundance in this study. Finally, the heatmap gives standardized values (Svalue) of depth and of the logarithmic concentrations of CH₄ and HS. Non-available data are represented by striped squares.

FIGURE 7

Figure 7. Relative abundance of the taxonomic groups that contain sequences associated to abundant eukaryotic OTUs. Overall, 140 abundant OTUs were retrieved, but 39 of them, abundant at the reference site, were mostly associated to taxa that are suggested to be allochthones and have fallen from surface waters. Therefore, for each sediment community library, these 39 OTUs were separated and are shown in the left bar charts (OTUs shared), while the remaining GHP OTUs are shown to the right (OTUs GHP). Bray-Curtis dissimilarity hierarchical clustering of the microbial communities at selected sediment depths was based on the GHP OTUs and separated these in four different clusters (E1, E2, E3, and E4). The eukaryotic communities from MC_902 at 15 cmbsf, from MC_900 at 29 cmbsf, from MC_919 at 35 cmbsf, and from MC_922 at 15 cmbsf strongly diverged and were therefore not included in these clusters. Leaves correspond to the core ID of the sediment layer and its depth [CoreID – depth (cm)]. “Other” corresponds to the relative abundance of sequences that were associated to taxonomic groups that are not illustrated in the figure.

For the bacterial domain, five community clusters, designated B1, B2, B3, B4, and B5, were identified for the GHP sediments (Figure 6). The rare biosphere represented by the non-abundant OTUs composed of a large fraction of all the bacterial communities and particularly for the clusters B1 and B2. Within these two clusters, the rare biosphere composed of an average of 76.7% of the bacterial sequences. Among the abundant OTUs, sequences within the cluster B1 were mostly assigned to the Gammaproteobacteria (5.2%), the Verrucomicrobiota (3.2%), and the Campilobacterota Sulfurimonadaceae (1.6%). The cluster B2 had stronger presence of Desulfobacterota Desulfobacterales (4.5%), including sequences associated to the cold seeps clade SEEP-SRB1 and Desulfobulbales (3.1%), in addition of Bacteroidota (4.7%). Cluster B3 represented sediment communities retrieved at the gas flare (core BC_1029) and was dominated by the Campilobacterota Sulfurovaceae (20.2%) and Sulfurimonadaceae (13.5%). Throughout all cores, the Sulfurovaceae and Sulfurimonadaceae were strictly represented by the genera Sulfurovum and Sulfurimonas, respectively. Additionally, B3 was characterized by the occurrence of Dissulfuribacterales (1.5%), mainly due to an OTU of the SEEP-SRB2 group, Desulfatiglandales (1.6%) and Desulfobacterales (2.9%). Remaining abundant taxa of the cluster B3 were assigned to the Bacteroidota (7.8%) and the Chloroflexi Anaerolinaeae (3.2%). It is also to be noted the presence of the Gammaproteobacteria Methylococcales in cluster B3 (1.2%). Communities within the cluster B4 primarily hosted sequences assigned to the Desulfobacterota (9.0%), largely included within the Desulfobacterales (5.2%), and the Chloroflexi Anaerolinaeae (6.1%) and Dehalococcoidia (3.6%). Additionally, abundant OTUs characterizing the cluster B4 were assigned to the Bacteroidota (4.3%), the Caldatribacteriota Japan Sea 1 (JS1) clade (4.8%), the Campilobacterota Sulfurimonadaceae (2.2%), and Sulfurovaceae (1.6%). The cluster B5 was dominated by the Desulfobacterota (19.6%), including representatives of Desulfobacterales (6.0%), Desulfatiglandales (4.4%), and Dissulfuribacterales (9.1%). One OTU assigned to SEEP-SRB2 and two OTUs assigned to SEEP-SRB1 composed 9.1 and 4.7% of the overall sequences, respectively (Figure 5). In comparison to other bacterial clusters, the cluster B5 was also characterized by a higher relative abundance of the Caldatribacteriota JS1 (10.7%) in addition to the Chloroflexi Dehalococcoidia (11.7%) and Anaerolinaeae (4.4%). Sediment samples clustering within the groups B3–B5 were mostly dominated by abundant OTUs, as the rare biosphere composed ca. 43% of the overall sequences.

With the eukaryotic primers, 39 abundant OTUs were retrieved at the reference core MC_898 and they composed from 20 to 100% of the sequences in all sediment communities. In Figure 7, these 39 OTUs are presented separately from the 101 eukaryotic OTUs retrieved exclusively at the GHPs. Beta diversity in the relative abundances of the taxonomic groups of these 101 OTUs retrieved in sediment communities at the GHPs site resulted in four clusters, designated as E1, E2, E3, and E4 (Figure 7). The proportions of sequences assigned to these OTUs varied between clusters, where an average of 8.5, 19.8, 29.8, and 24.1% of the sequences for the clusters E1, E2, E3, and E4 were assigned to them, respectively. In cluster E1, Cercozoa and ciliates corresponded respectively to 2.7 and 1.5% of the overall sequences. Within the cluster E2, these groups were more abundant, and their relative abundances increased to 8.2% for the Cercozoa and to 5.4% for the ciliates. For clusters B1 and B2, sequences were primarily assigned to an unclassified group of Cercozoa, while the class Spirotrichea primarily dominated the ciliates. Within the cluster E3, the taxonomic diversity was higher than for E1 or E2. Other cercozoan groups, such as Granofilosea, Phytomyxea, and Thecofilosea, in addition to the ciliates classes Armophorea and Conthreep, are frequently seen in higher abundances. In addition to Cercozoa and ciliates, abundant taxa exclusive to these sediment layers included representatives of the Holozoa, uncultivated marine stramenopiles (MAST) groups 6 and 12, in addition to the fungi (Ascomycota, Basidiomycota, and Chrytridiomycota). Sediment samples clustering within the cluster E4 were characterized by a higher proportion of Apicomplexa among the OTUs_.

Similarly to the distribution of the 101 eukaryotic OTUs presented above, the community structure of the 39 OTUs also thriving at the reference site varied between the clusters. The relative abundances of Ochrophyta were lower in clusters E2, E3, and E4, which are more predominantly composed by Cercozoa and ciliates. Finally, alpha diversity metrics that were used to assess biodiversity richness and evenness and the taxonomic composition for all domains within each sediment community are available as Supplementary Information (Supplementary Tables 8, 9 and Figures 1–3).

Distribution and Co-occurrence of the Domain Clusters

The community clusters showed particular patterns of co-occurrence between each domain, especially for the prokaryotes (Supplementary Figure 4). For instance, 12 of the 13 sediments communities within cluster A1 were associated with the bacterial cluster B1. The pairs A2/B2, A3-A4/B4, and A6/B5 were also commonly co-occurring. However, concomitance patterns between prokaryotic and eukaryotic clusters were less supported. Still, the eukaryotic cluster E1 usually fell together with the pairs A1/B1 or A2/B2. The clusters E2 and E3, instead, coincided with the pairs A3-A4/B4 and A6/B5, respectively. The paired clusters A1/B1 were retrieved at the surface of nearly all sediment cores while the clusters A2/B2 generally corresponded to the subsurface communities at the reference site and at cores taken toward the edge of a GHP. Pairs of A3/B4 or A4/B4 occurred below the sediment surface at the apex of GHP 1 (cores MC_900 and MC_902) and of GHP 3 (MC_1061 and MC_1062). The pair A6/B5 occurred in subsurface sediments at the apex of GHP 1, but also toward the outskirt of the GHPs at the surface of MC_918 (GHP 1) and in subsurface sediments of GHP 5 (core MC_922). Finally, the microbial communities retrieved at the gas flare (core BC_1029) of the GHP 3 could not be related to other communities at the GHPs site for all domains of life and clustered separately. Communities from all sediment depths at BC_1029 clustered within A5, B3, and E4.

Impact of Environmental Conditions on the Microbial Community Structure

The community clusters for the two prokaryotic domains demonstrated a profile primarily related to sediment depth and methane availability (Figure 8). The impact of measured environmental parameters on the dissimilarity between the different prokaryotic communities, observed through the formation of six archaeal and five bacterial community types, was assessed through dbRDA. Overall, the unconstrained proportions of the two principal axes (RDA 1 and 2) explained 43.71–62.52% of the dissimilarity between the different prokaryotic communities and were all significant (Figure 8). Depth correlated negatively with the prokaryotic community types A1 and B1 while CH₄ concentrations drove the dissimilarity between the other community types. At all GHPs, A2 and B2 correlated negatively with CH₄ concentrations, while A3-4-5-6 and B3-4-5 correlated positively. At GHP 1, these community types were also impacted by higher concentrations of ΣHS, while types A2 and B2 thrived in sediments richer in Fe²⁺ and SO₄^2–.

FIGURE 8

Figure 8. The impact of different environmental parameters on the archaeal and bacterial community structure within the sediments of the different GHPs assessed through dbRDA. A distance matrix was calculated based on the Bray-Curtis dissimilarity index from the composition of abundant archaeal and bacterial OTUs for each GHP. The correlation between the environmental variables and the built distance matrices are presented by biplots. The unconstrained proportion for each axis explaining the variability in a distance matrix is presented in percentage along the axis. Permutation tests were used to assess the solidity of the analyses and axes with a * were found significant.

Sulfurimonas and Agglutinated Foraminifera Relationship

In general, we found lower numbers of agglutinated foraminifera at habitats characterized by higher densities of the sulfide-oxidizing Suflurimonas. The relationship between the logarithm of the number of resampled Sulfurimonas sequences at the sediment surface and the density of agglutinated foraminiferal species showed a significant (F = 43.122, p-value = 0.007183), negative, and linear correlation (Supplementary Figure 5).

Discussion

Community Types Distribution Across the Pingos

Our first objective was to test the hypothesis that variations in the community structure occur along a radial gradient from the apex of the GHPs, which was expected to concentrate most of the gas seeping activity (Serov et al., 2017). Investigating the microbial communities thriving along spatial and depth pingos gradients led to the distinction of different community clusters for each domain of life (Figures 4, 6, 7). CH₄-rich sediments hold distinct community clusters (A3–A6, B3–B5, E2–E4) while communities retrieved in CH₄-poor sediments were more similar to the reference site (Figure 8). According to our hypothesis, CH₄-rich sediments were recovered from coring locations close to the apex of GHP 1 (MC_900 and MC_902) or GHP 3 (BC_1029) where active gas flares were visible. However, we did also find high dissolved methane concentrations sediments hosting the CH₄-rich community clusters we have described at the edge of GHPs 1 (MC_918) and 5 (MC_922; Supplementary Figure 4). This unpredicted spatial distribution of the different microbial community types at the GHPs was further supported through the observed significant negative correlation between the relative abundance of Sulfurimonas and the density of agglutinated foraminifera on the seafloor (Supplementary Figure 5). While Sulfurimonas is a genus that is often retrieved in higher relative abundances in CH₄-rich sediments (Figure 6; Niemann et al., 2013; Bomberg et al., 2015), agglutinated foraminifera are known to be sensitive to CH₄-rich environments (Panieri and Sen Gupta, 2008; Martin et al., 2010; Dessandier et al., 2019).

The use of these two independent methods further confirmed that there was no radial gradient at the GHPs. This contrasted thereby with earlier studies on active mud volcanoes where the community composition and the nature of the dominating methane oxidizers varied along concentric zones around the apex of the structure (Niemann et al., 2006; Lösekann et al., 2007; Lee et al., 2019). Instead, across the GHPs, community types were scarcely distributed and mainly depth and the availability of CH₄ appeared to drive the transition between them (Figure 8). Furthermore, changes in community composition at the GHPs occurred on a smaller scale than at the HMMV, where the identified concentric zones extended over tenth to hundreds meters. In our study, nearby sediments cores MC_918 and MC_919, or BC_1029 and MC_1061, were less than 40 m apart, but the first hosted a community type dominated by ANME-1 while the latter was more similar to the reference site. This suggests that the microbial community spatial succession at these pingos is still not yet fully grasped. Thereby, further investigations on the variability of the microbial community composition should be addressed at a higher site resolution.

Microbial Biodiversity Across the Study Area

Our second objective was to describe the microbial biodiversity at the GHPs and to identify key taxa influenced by this CH₄-rich environment. Overall, the communities presented different assemblages, depending on their vertical positioning in the sediment matrix; i.e., surface, a few cm below the seafloor, in CH₄-rich sediments, or at the gas flare (BC_1029). The variability in the structure of eukaryotic communities and the nature and quantities of Foraminifera at the GHPs were analyzed differently than for prokaryotes. We therefore discuss the composition of the prokaryotic and eukaryotic communities within the different sediment habitats separately.

Prokaryotes

Sediment characterized by a CH₄ depletion and ΣHS increase hosted a microbial community dominated by ANME and SRBs, strongly suggesting ongoing AOM. The archaeal community was primarily dominated across all GHPs by the anaerobic CH₄ oxidizing group ANME-1 (Figures 4, 5). Interestingly, methanotrophic communities primarily driven by ANME-1 have been less frequently observed than by ANME-2, or were found only in deeper sediments (Girguis et al., 2005; Ruff et al., 2015; Gründger et al., 2019). Our understanding of the factors favoring the growth of the different ANME groups is still limited. Their tolerance levels to various environmental factors and the impact of CH₄ flux rates on their growth rate have been two common orientations used by studies to investigate their biogeography. Within the first orientation, it is suggested that the ANME-1 would be more tolerant to broader ranges of environmental conditions, and could predominate over ANME-2 in low SO₄^2– and high HS^– environments (Timmers et al., 2015). These different tolerances to the presence of SO₄^2– and HS^– has been suggested to explain vertical successions in dominance of these groups along different SMTZ (Roalkvam et al., 2011; Biddle et al., 2012; Ruff et al., 2015). However, at the GHPs, although ANME-2 and ANME-3 were also detected, their relative abundances remained low, and there was no clear vertical transition in the nature of the dominant ANME group along the SMTZ in cores MC_900 and MC_902. This could suggest that other factors at the GHPs favor the growth of ANME-1 and/or inhibit the proliferation of ANME-2. Within the second orientation, observations were made at the Hydrate Ridge or the Gulf of Mexico that ANME-2 groups were more commonly retrieved in areas with highly active CH₄ seepage (Vigneron et al., 2013, 2019). In our study, although ANME-1 still predominated the methanotrophic community near the gas flare (BC_1029), the relative abundances of ANME-2 groups were in contrast higher than in other clusters. However, this hypothesis would contradict previous observations where ANME-2 demonstrated higher growth rates than ANME-1 at low CH₄ flux rates (Girguis et al., 2005). Beyond these two hypotheses presented above, the hydrographic conditions above the GHPs could also induce an additional set of environmental constraints, as the bottom-water temperature seasonally varies (Ferré et al., 2020). This creates fluctuations in both CH₄ seeping activity from the sediments and subsequently CH₄ oxidation rates in the water column. This seasonality in CH₄ seepage activity could potentially also impact the selection of the ANME groups. The biogeography of ANME groups remains therefore still unclear. With its five GHPs presenting different CH₄ flux history and its multiple ecological niches, the GHPs, combined with the usage of appropriate tools for sampling sediments at a higher precision, present thereby an ideal site to provide further insights into the distribution of ANME groups.

Furthermore, to mediate AOM, ANME groups require an electron acceptor, such as sulfate, and have therefore been frequently observed in consortia with microorganisms capable of reducing these compounds. The ANME-1 group have regularly been assigned to the uncultured groups of SEEP-SRB1 and SEEP-SRB2, where both are detected in CH₄-rich sediments at the GHPs. In our study, the relative abundance of Desulfobacterota was higher in microbial communities dominated by ANME groups (Figure 6). Furthermore, the decreasing concentration of SO₄^2– with depth in CH₄-rich sediments, combined with an increasing availability of ΣHS, strongly suggested the use of sulfate as the electron acceptor for AOM. However, across all the GHPs, there was no positive correlation between the relative abundance of ANME-1 and a particular SRB group, either SEEP-SRB 1 or 2, further supporting the hypothesis that ANME-1 could metabolize CH₄ alone (Figures 4–6). Indeed, it was observed that ANME-1 could perform both AOM and sulfate reduction within the same cell (Milucka et al., 2012) and the detection of F420-dependent sulfite reductase in ANME-1 communities may be part of this novel pathway (Vigneron et al., 2019). Nevertheless, a previous study could not find a correlation of ANME-1 and the abundance of dissimilatory sulfite reductase, an essential enzyme for active SRB (Vigneron et al., 2019), demonstrating that ANME-1 may not be able to perform SR. Finally, a different explanation of the absent correlation between ANME-1 and SRB groups could be due to the usage of intercellular wires forming cell-to-cell connections for electron transfers, a hypothesis supported by the detection of genes expressing for extracellular cytochrome production, between distanced ANME-1 and SRB cells (Wegener et al., 2015). Our results, based on the sulfate and sulfide profiles, advocate an anaerobic oxidation of CH₄ supported by the reduction of sulfate, but the role of Desulfobacterota and its relation with the ANME groups remain unclear.

While AOM is mediated by ANME in anaerobic environment, methanotrophy in an aerobic environment is primarily performed by distinct bacterial groups (Hanson and Hanson, 1996; Knief, 2015). In our study, higher concentrations of CH₄ than at the reference site were detected at the surface of some sediment cores collected at the GHPs. However, despite the availability of oxygen suggested by the presence of aerobic taxonomic groups, aerobic bacterial methanotrophs were barely detected. We retrieved abundant Verrucomicrobiales OTUs at the surface of most sediment cores, but their assigned family Rubritaleaceae is not known to include aerobic methanotrophs. Aerobic methanotrophs (Methyloccocales) could only be detected at the surface of BC_1029, collected near the gas flare, but this OTU was composed of only 1.2% of all bacterial sequences. Surprisingly, the apparent rarity of aerobic CH₄ oxidizers is contrasting to most seep sites where they were found when both CH₄ and O₂ are present (Lösekann et al., 2007; Roalkvam et al., 2011; Ruff et al., 2015). Nevertheless, we cannot disregard that the near absence of aerobic methanotrophs in our amplicon libraries could be caused by the choice of primers used (McDonald et al., 2008). Different approaches, including the use of primers targeting functional genes such as pmoA, would be required to improve the study of the biodiversity of aerobic methane oxidizers. Finally, CH₄-rich sediments also harbored higher relative abundances of other groups, but which are likely not directly involved in the AOM. Chloroflexi, the Caldatribacteriota JS1, and Campilobacterota groups were also in higher abundance in CH₄-rich sediments than at other sediment layers. Similarly to the distribution of ANME groups, these bacterial groups showed different relative abundances between CH₄-rich sediments collected at the gas flare to the other samples. While most communities in CH₄-rich sediments demonstrated high proportions of Chloroflexi and JS1, the bacterial communities at the gas flare was primarily dominated by sulfide oxidizing bacteria (Figure 6). More precisely, two Campilobacterota genera mediating the oxidation of sulfur, sulfide or thiosulfate, Sulfurimonas and Sulfuvorum, were found in abundance. These genera are commonly found in abundance near hydrothermal plumes and in diffusive flow sediments, as well as at cold seeps (Yamamoto and Takai, 2011; Adams et al., 2013), while sulfide oxidization in marine sediments tends to be driven primarily by Alphaproteobacteria or Gamma proteobacteria (Lenk et al., 2011). In our study, similar observations suggest that these bacteria play an important role in sulfur cycling and largely dominated the bacterial communities at the gas flare, in comparison to the other sites.

In the absence of CH₄, the sediment microbial composition at the GHPs was highly similar to the reference site and was primarily driven by depth (Figure 8). Depth is likely influencing the shape of microbial communities at the GHPs through the presence or absence of oxygen, a parameter well-known to shape the structure of microbial communities in sediments (Fenchel and Finlay, 2008). Surface sediments were primarily dominated by the aerobic ammonia-oxidizing archaea (AOA) Nitrosopumilaceae that plays, along with ammonia-oxidizing bacteria, an important role in the transformation of nitrogen compounds in marine systems, including cold seeps or at hydrothermal vents (Könneke et al., 2005; Dang et al., 2009; Miyazaki et al., 2009; Stahl and de la Torre, 2012). In deeper sediments, the archaeal community (A2) was dominated by Woesearchaeales and Bathyarchaeia (Figure 4). The most abundant OTU of the 38 associated to the Woeserchaeia across all clusters was found predominantly at nearly all sediment layers below the seafloor, including in the CH₄-rich sediments. As oxygen availability is suggested to be the main factor determining the nature of the thriving Wosearchaeales (Liu et al., 2018), its detection in deeper sediments likely suggest an anoxic ecotype that may be involved in a fermentation-based lifestyle (Castelle et al., 2015). Bathyarcheia, previously known as the Miscellaneous Crenarchaeotal Group (MCG), and the thermoplasmatales MBG-D are globally abundant in marine sediments. The detection of protein-degrading enzymes suggest a role in organic matter anaerobic degradation (Webster et al., 2010; Kubo et al., 2012; Lloyd et al., 2013). The relative abundance of OTUs assigned to the Desulfobacterota within the bacterial community increased with depth, but remained lower than in CH₄-rich sediments (cluster B5). The presence of these Desulfobacterota groups are common in marine sediments as they play a major role in mineralizing organic matter through sulfate reduction (Jørgensen, 1982; Abu Laban et al., 2015; Robador et al., 2016).

Eukaryotes

Microbial eukaryotic communities at cold seeps have received less attention than the prokaryotes, despite their active role as part of bacterial mat type habitats for instance, or the capacity of some to harbor sulfur oxidizing bacteria (Buck and Barry, 1998; Buck et al., 2000). In this study, we investigated protists and fungi based on the V4 region of the 18S rDNA and have identified key taxonomic groups thriving at the GHP sites. Across all sediments layers, large fractions of sequences that clustered into OTUs were assigned to Metazoa and to sedimenting allochthonous cells. The removal of these sequences likely affected the following analyses of the GHPs eukaryotic communities. Therefore, we assessed separately the 39 OTUs proliferating at the reference site from the 101 OTUs found in abundance only at the GHPs to highlight eukaryotic taxa thriving in CH₄-rich sediments. Communities clustering in E1 demonstrated high similarity to the reference site and could be retrieved at different distances from the apex of all GHPs, but were limited to sediments characterized by low CH₄ concentrations. We thereby demonstrated that in the absence of CH₄, eukaryotic communities across the GHPs have similar composition than to the reference site. In contrast, clusters E2–E4 were retrieved in or near CH₄-rich sediments and demonstrated higher relative abundances of OTUs that are absent or barely found at the reference site.

Within the cluster E2, these OTUs were primarily assigned to ciliates and Cercozoa. We also noted that within the 39 subtracted OTUs, the fraction of alveolates and Cercozoa increased and even surpassed their relative abundances in E3 and E4. Higher densities of prokaryotes involved directly or indirectly in AOM can be a food source for these potential heterotrophic eukaryotes, but their growth in communities clustering in E3 and E4 may be limited by the toxicity of sulfidic conditions (Massana et al., 1994; Coyne et al., 2013). Communities clustering in E3 occurred primarily in CH₄-rich sediments with the prokaryotic communities of clusters A4-A6/B5, composed of taxa involved in AOM and sulfate reduction (Supplementary Figure 4). Nearly all communities within the cluster E3 hosted the highest relative abundances of sequences associated to the 101 OTUs that are exclusively found in abundance at the GHPs (Figure 7). The contrast in these relative abundances, in comparison to the cluster E1, demonstrates the impact of CH₄ on the eukaryotic diversity. The assignment of these OTUs was strongly heterogeneous as several taxonomic groups, such as the Protoalveolata Syndiniales, were present only in few communities (Figure 7). The eukaryotic communities within this cluster were also characterized by the emergence of fungal taxonomic groups. Communities within the cluster E3 were especially affected by the proportion of sequences associated to Metazoa, as on average 40% of the sequences were assigned to this taxon and had to be removed. Thereby, while the heterogeneity in the structure of the communities clustering in E3 could be caused by local conditions, we cannot rule out that it may be due to limitations in the coverage of the eukaryotic biodiversity. Communities at the gas flare (BC_1029), similarly as for the prokaryotes, hosted a distinctive eukaryotic biodiversity clustering exclusively in E4. Among the sequences assigned to the OTUs exclusively abundant at the GHPs within E4, most were primarily assigned to Apicomplexa (up to 15%). Apicomplexa are parasitic alveolates, but the nature of potential hosts at the gas flare remains unknown. Overall, our results demonstrated that changes in the eukaryotic biodiversity occur in CH₄-rich sediments. Using different approaches, such as targeting specific genes or using blocking primers, may provide a more accurate profile of eukaryotic biodiversity at the GHPs. These investigations would further improve our understanding on the role of these protists and fungi at the GHPs site on the microbial community, the biogeochemical cycles, and on food web structures.

Overall, our approach suggests that CH₄ and oxygen are two key factors influencing the microbial community structure. Nevertheless, communities within a cluster had up to approximately 60% similarity and the dendrograms (Figures 4, 6, 7) present additional sub-clusters at higher thresholds. It advocates therefore for additional factors influencing the distribution patterns of the microbial taxonomic groups at the GHPs site. Thereby, our study revealed that the GHP ecosystem has to be considered in further investigations as a myriad of ecological niches. In this perspective, the distance between the cores (approx. 20 m) at a GHP is likely too long to investigate gradual changes in microbial communities in relation to fluxes of CH₄. Designing an approach at a small scale may better fill these gaps of knowledge.

Summary and Concluding Remarks

This study shows that both prokaryotic and eukaryotic communities at the GHPs formed a unique structure influenced by the complex distribution of CH₄ seepage. The distribution of the community types presented similar chaotic patterns and methane oxidizing communities could be retrieved at different locations over a GHP. In CH₄-rich sediments, AOM seemed to be primarily driven by a single OTU associated to ANME-1 and had no correlation with a group of SRB. This further supports the hypotheses that ANME-1 can mediate AOM alone or use different sources of electron receptors. Our approach also illustrated that at the GHPs site, metabolites of AOM, such as sulfide and organic compounds, likely explain the predominance of additional taxa, including the Campilobacterota, the thermoplasmatales MBG-D, and the Bathyarchaeia. Eukaryotic communities in the CH₄-rich sediments had a dominance of heterotrophic ciliates and Cercozoa, likely benefiting from the higher abundances of prokaryotes as a food source. The retrieval of these taxa, distributed specifically among the GHPs, suggests a complex functional microbial system supported by, or contributing to, the local oxidation of CH₄.

Data Availability Statement

The datasets generated for this study can be found in the Sequence Read Archive Genebank as BioProject accession number PRJNA593930.

Author Contributions

VC, MS, FG, and HN initially designed the project. VC, FG, MS, HN, and P-AD contributed to the sampling. VC and P-AD performed the laboratory manipulations, sequences analyses, and statistics with advice from DK, MS, FG, and HN. VC wrote the manuscript with input from DK, MS, HN, FG, P-AD, and GP. All authors contributed to the article and approved the submitted version.

Funding

The study was funded by the Research Council of Norway through the Centre of Excellence for Arctic Gas Hydrate, Environment, and Climate (grant number: 223259).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors are grateful for the help from the crew of the R/V Helmer Hanssen and would like to acknowledge the contribution of the chief scientists during the CAGE 16-5 and 17-2 campaigns. The authors would also like to acknowledge the technical assistance of Wei-Li Hong, Haoyi Yao, Simone Sauer, Matteus Lindgren, and Pavel Serov in the analyses of the physico-chemical parameters. Finally, the authors would like to recognize members of the Centre for Arctic Gas Hydrate and Environment (The Arctic University of Norway, Tromsø) for their expertise. Antje Boetius and Frank Wenzhöfer (AWI, Bremeherhaven and MPI, Bremen, Germany) are thanked for their support with blade corers. The computations were performed on resources provided by UNINETT Sigma2 – the National Infrastructure for High Performance Computing and Data Storage in Norway, account nos. NN9639K and NS9593K. Finally, the publication charges for this article have been funded by a grant from the publication fund of UiT The Arctic University of Norway.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2020.01932/full#supplementary-material

Footnotes

References

Abu Laban, N., Tan, B., Dao, A., and Foght, J. (2015). Draft Genome Sequence of Uncultivated Toluene-Degrading Desulfobulbaceae Bacterium Tol-SR, Obtained by Stable Isotope Probing Using [13C6] Toluene. Genom. Announc. 3, e1423–e1414. doi: 10.1128/genomeA.01423-1414