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

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 4 ), 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 4 forms gas hydrates, i.e., an icelike lattice comprising molecules of CH 4 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 4 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 4 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(Niemann et al., , 2009Reeburgh, 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 Rossel et al., 2011), salinity (Maignien et al., 2013), or CH 4 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., 2010Ettwig et al., , 2016Hu 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., 2007Takishita et al., , 2010Niemann 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 . 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 . 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 4 seepage has been active for several thousand years . 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 Hong et al., 2018) or the biodiversity of larger fauna Å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 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 4 -rich environments, demonstrating the uniqueness of this ecosystem.

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 . A group of 5 GHPs (∼10 m high, 500 m in width) distributed on the seabed over a 2.5 km 2 area were recently found 1 . Hydroacoustic surveys and real time visually guided observations with a TowCam-Multicore (see footnote 1) System (TC-MC) and Remotely Operated Vehicle (ROV) dives 2 have revealed acoustic flares of gas bubbles consisting predominantly of CH 4 and being emitted from 4 of the 5 GHPs. One structure (GHP 5) did not show any visible flares or gas hydrate in sediments . 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 .

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 3 , 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 4 gas flare at GHP 3 in June 2016.

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 4 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 2+ ), alkalinity, total sulfide ( HS), sulfate (SO 4 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 2+ were measured onboard by titration and by spectrophotometry, respectively (Hong et al., 2017). SO 4 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 2+ 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 R 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) 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. 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 4 as BioProject accession number PRJNA593930.
Paired-end reads were meticulously processed and the workflow was derived from the USEARCH suggested protocol 5 . 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 4 -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 4 -rich sediments. To ensure that the changes in relative abundances of these taxonomic groups were caused by the presence of CH 4 , we compared results from DNA sequences with an independent proxy for surface CH 4rich 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 4 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 3 m 2 ). 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.

Environmental Characterization and Geochemistry
At the reference site, CH 4 was nearly absent, gas flares were not detected on the echosounder, and CH 4 sediment concentrations did not exceed 4 µM (Figure 2). HS remained undetectable throughout the reference core, while measured concentrations of SO 4 2− slightly decreased from 28 mM at the sediment surface to 26 mM at 11 cm below seafloor (bsf), correlating 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 4 , SO 4 2-, and HS at these sediment depths are given. Black bars correspond to the sediment layers from which the DNA was extracted and sequenced.
Frontiers in Microbiology | www.frontiersin.org 6 September 2020 | Volume 11 | Article 1932 to seawater concentration of the Barents Sea. The seafloor was muddy and authigenic carbonates were not observed (Supplementary Table 4). At GHP 1, gas flares and high CH 4 sediment concentrations were suggestive of high CH 4 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 4 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 4 2− with depth than at the reference site as concentration dropped below 5 mM at 15 cmbsf. With the decrease in SO 4 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 4 and HS increased with depth, but at lower concentrations than at cores taken near the apex of GHP 1. The SO 4 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 4 (yet still slightly higher than at the reference site) were detectable and SO 4 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 4 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 4 gas flare (Figure 1). The four other cores from GHP 3 had lower CH 4 concentrations than BC_1029 (<15 µM). Still, the cores MC_1061 and MC_1062, located close to the GHP 3 apex, had higher CH 4 concentrations than cores MC_1063 and MC_1065, collected near the edge and outside GHP 3, respectively. SO 4 2− maximum concentrations in the surface sediment layers were in the range of 27-28 mM for all cores, but the SO 4 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 4 2− showed a similar pattern to the reference site. Fe 2+ 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).
At GHP 5, similar CH 4 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 4 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 2+ , alkalinity, and DIC are available as Supplementary Information (Supplementary Tables 4, 5).
FIGURE 3 | Geochemical profiles of the sediment core BC_1029 collected at GHP 3 near a CH 4 gas flare. Profiles for CH 4 and SO 4 2at these sediment depths are shown. Black bars correspond to the sediment layers from which the DNA was extracted and sequenced.

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.
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 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 4 and HS. Non-available data are represented by striped squares.
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.

Distribution and Co-occurrence of the Domain Clusters
The community clusters showed particular patterns of cooccurrence 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 4 concentrations drove the dissimilarity between the other community types. At all GHPs, A2 and B2 correlated negatively with CH 4 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 2+ and SO 4 2− .

Sulfurimonas and Agglutinated Foraminifera Relationship
In general, we found lower numbers of agglutinated foraminifera at habitats characterized by higher densities of the sulfideoxidizing 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).

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 . 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 4 -rich sediments hold distinct community clusters (A3-A6, B3-B5, E2-E4) while communities retrieved in CH 4 -poor sediments were more similar to the reference site (Figure 8). According to our hypothesis, CH 4 -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 4 -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 4 -rich sediments (Figure 6; Niemann et al., 2013;Bomberg et al., 2015), agglutinated foraminifera are known to be sensitive to CH 4 -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 4 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.
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.

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 4rich 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 4 -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 4 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 4 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 4 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 4 2− and high HS − environments (Timmers et al., 2015). These different tolerances to the presence of SO 4 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 4 seepage (Vigneron et al., 2013(Vigneron et al., , 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 4 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 4 seeping activity from the sediments and subsequently CH 4 oxidation rates in the water column. This seasonality in CH 4 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 4 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 4 -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 4 2− with depth in CH 4 -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 4 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 4 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 4 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 4 oxidizers is contrasting to most seep sites where they were found when both CH 4 and O 2 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 4 -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 4 -rich sediments than at other sediment layers. Similarly to the distribution of ANME groups, these bacterial groups showed different relative abundances between CH 4 -rich sediments collected at the gas flare to the other samples. While most communities in CH 4 -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 4 , 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 wellknown 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 ammoniaoxidizing 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 4 -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 4 -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 4 -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 4 concentrations. We thereby demonstrated that in the absence of CH 4 , eukaryotic communities across the GHPs have similar composition than to the reference site. In contrast, clusters E2-E4 were retrieved in or near CH 4 -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 4 -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 4 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 4 -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 4 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 4 . 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 4 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 4 -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 4 -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 4 .

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).