Glacial Runoff Promotes Deep Burial of Sulfur Cycling-Associated Microorganisms in Marine Sediments.

Marine fjords with active glacier outlets are hot spots for organic matter burial in the sediments and subsequent microbial mineralization. Here, we investigated controls on microbial community assembly in sub-arctic glacier-influenced (GI) and non-glacier-influenced (NGI) marine sediments in the Godthåbsfjord region, south-western Greenland. We used a correlative approach integrating 16S rRNA gene and dissimilatory sulfite reductase (dsrB) amplicon sequence data over six meters of depth with biogeochemistry, sulfur-cycling activities, and sediment ages. GI sediments were characterized by comparably high sedimentation rates and had "young" sediment ages of <500 years even at 6 m sediment depth. In contrast, NGI stations reached ages of approximately 10,000 years at these depths. Sediment age-depth relationships, sulfate reduction rates (SRR), and C/N ratios were strongly correlated with differences in microbial community composition between GI and NGI sediments, indicating that age and diagenetic state were key drivers of microbial community assembly in subsurface sediments. Similar bacterial and archaeal communities were present in the surface sediments of all stations, whereas only in GI sediments were many surface taxa also abundant through the whole sediment core. The relative abundance of these taxa, including diverse Desulfobacteraceae members, correlated positively with SRRs, indicating their active contributions to sulfur-cycling processes. In contrast, other surface community members, such as Desulfatiglans, Atribacteria, and Chloroflexi, survived the slow sediment burial at NGI stations and dominated in the deepest sediment layers. These taxa are typical for the energy-limited marine deep biosphere and their relative abundances correlated positively with sediment age. In conclusion, our data suggests that high rates of sediment accumulation caused by glacier runoff and associated changes in biogeochemistry, promote persistence of sulfur-cycling activity and burial of a larger fraction of the surface microbial community into the deep subsurface.

Introduction

Arctic fjords with marine-terminating glaciers constitute an important interface for freshwater and sediment influx from land into the sea, thereby influencing the physical and chemical conditions in the coastal marine ecosystems (Svendsen et al., 2002; Etherington et al., 2007). The high influx of sedimentary materials, e.g., minerals, terrigenous organic matter and metals, in glacier-associated fjords has a strong effect on the distributions of the benthic microbial communities (Park et al., 2011; Bourgeois et al., 2016; Buongiorno et al., 2019). Increased water turbidity in close proximity to the glacier can negatively influence surface water primary production (Etherington et al., 2007; Zajączkowski, 2008), which leads to lower organic matter availability in the underlying sediments (Bourgeois et al., 2016). High sediment accumulation rates often seen in such glacier-proximal environments are also limiting benthic life. On the other hand, glacial meltwater also provides an important source of dissolved nutrients, which can stimulate phytoplankton growth beyond the high turbidity zone (Statham et al., 2008; Meire et al., 2015; Sørensen et al., 2015). Consequences of increased primary production together with strong sediment supply are a net CO2 uptake in glaciated fjords as well as rapid burial of fresh detrital phytoplankton biomass to the underlying sediments (Meire et al., 2015; Smith et al., 2015; Bourgeois et al., 2016). Fjord sediments also receive significant amounts of terrigenous organic matter as evidenced by high C/N ratios of the sediment organic matter pool (Goñi et al., 2013; Wehrmann et al., 2014). With these large inputs of both marine and terrigenous organic matter, the ultimate role of glaciated fjord sediments in the global carbon cycle depends on the extent to which the large organic matter inputs are degraded, and thus there is a need to better understand constraints on microbial community structure and degradation potential.

The labile fraction of organic matter in marine sediments is initially degraded by microorganisms with hydrolytic and fermenting capabilities (Müller et al., 2018). These largely unidentified microbial species typically excrete a wide range of enzymes, which enable rapid organic matter turnover even in cold arctic sediments (Arnosti et al., 1998; Teske et al., 2011). Fermentation products released by primary organic matter degrading microorganisms are either further degraded by secondary fermenters or are mineralized completely to CO2 via microbial respiration (Arndt et al., 2013).

Of prime importance in microbial organic matter degradation is the respiratory reduction of sulfate to sulfide, which facilitates up to 69% of the total organic matter mineralization in Arctic fjord sediments (Sørensen et al., 2015). Other key electron acceptors are metals such as iron (oxyhydr)oxides and manganese oxides (Froelich et al., 1979; Thamdrup et al., 1994), which can be introduced in high amounts via glacial runoff to marine sediments and are subjected to redox cycling (Bhatia et al., 2013; Wehrmann et al., 2014; Laufer et al., 2016; Buongiorno et al., 2019). Fe(III) and Mn(IV) also facilitate the oxidation of reduced sulfur compounds (Zopfi et al., 2004; Wehrmann et al., 2014) and thereby fuel a cryptic sulfur cycle in glacier-influenced (GI) sediments (Wehrmann et al., 2017).

In sediments with comparably slow sediment accumulation, which typify non-glaciated Arctic shelf areas (Kuzyk et al., 2013), the rate of organic matter mineralization is highest at the surface and decreases significantly with increasing sediment depth (Lomstein et al., 2012; Kuzyk et al., 2017), reflecting a concomitant decrease in energy available for cell maintenance and growth (Starnawski et al., 2017). This depth gradient results in pronounced compositional changes in the benthic microbial community, which are driven by highly selective survival of microorganisms that are able to subsist in the energy-limited subsurface (Petro et al., 2017; Starnawski et al., 2017; Bird et al., 2019; Marshall et al., 2019). Here, we hypothesized that this strong environmental filtering effect would be attenuated in sediments with high rates of sedimentation (i.e., in coastal sediments affected by glacial runoff) and result in a different pattern of microbial community assembly with depth. Therefore, we investigated microbial diversity and community compositions of GI and non-glacier-influenced (NGI) coastal sediments in the Godthåbsfjord region of Greenland. Compositions of bacterial and archaeal communities and the community of putative sulfite/sulfate-reducing microorganisms (SRM), as analyzed by 16S rRNA and dissimilatory sulfite reductase (dsrB) gene amplicon sequencing, respectively, were compared across sediment samples of different depths and ages. Co-occurrence analyses of operational taxonomic units (OTUs) and correlations with biogeochemical data revealed key environmental factors that were driving the major community differences between GI and NGI sediments. These analyses also identified various uncultivated microorganisms that were associated with sulfur cycling. As hypothesized, our results demonstrate that glacial runoff exerts a strong influence on microbial community assembly processes and community functions in marine sediments.

Materials and Methods

Sediment Sampling

The sediment cores used in this study were collected using a gravity corer in 2013 during a research cruise on board of RV Sanna (Seidenkrantz et al., 2014). Sampling of the cores for this study was described previously (Glombitza et al., 2015). In brief, up to 6 m long gravity cores were recovered from four sites on the open shelf and within the Godthåbsfjord (Nuup Kangerlua) system in South West Greenland (Table 1 and Figure 1) in August 2013. The four stations can be broadly subdivided into two groups, the NGI stations 3 (St3-NGI) and 6 (St6-NGI), and the GI stations 5 (St5-GI) and 8 (St8-GI). St3-NGI (core SA13-ST3-20G) is located outside the fjord on the continental shelf of the Labrador Sea. St6-NGI (core SA13-ST6-40G) is situated in the Kapisigdit Kanderdluat, a side-fjord without glaciers. St5-GI (core SA13-ST5-30G) is located in the main channel of the Godthåbsfjord. Although St5-GI is not directly in front of a glacier, most glacier-derived material is transported toward the Labrador Sea across this site. St8-GI (core SA13-ST8-47G) is in very close proximity to a glacier front at the northernmost outlet of the Greenland ice sheet in the Kangersuneq fjord (Figure 1).

TABLE 1

Description of the sampling stations and cores, including sampling position, water depth and core length.

Station	Core name	Latitude (N)	Longitude (W)	Water depth (m)	Core length (cm)
St3-NGI	SA13-ST3-20G	64°26.7425′	52°47.6486′	498.2	587
St5-GI	SA13-ST5-30G	64°25.3479′	51°30.6209′	622.4	607
St6-NGI	SA13-ST6-40G	64°29.0604′	50°42.3240′	389	562
St8-GI	SA13-ST8-47G	64°40.7078′	50°17.4672′	475.8	569

FIGURE 1

Locations of the four sampling stations in the Godthåbsfjord region. The coordinates of the stations are indicated in Table 1. St3-NGI and St6-NGI, denote Non-glacier-influenced stations 3 and 6. St5-GI and St8-GI, denote Glacier-influenced stations 5 and 8. Label a, main fjord (Nûp Kangerdlua); b, Qôrnup Suvdlua; c, Umánap Suvdlua; d, Kapisigdlit Kanderdluat; e, Kangersuneq. The background map is a NASA Worldview satellite image (Terra MODIS Corrected Reflectance) of August 19, 2013.

The sediment at St3-NGI consists of olive-gray to olive-green mud, weakly laminated at the base and somewhat mottled in the upper part. St6-NGI and St5-GI sediment also consists of olive-gray to olive-green sediments, where radiographic images display a clear but not strong lamination. The core of St5-GI contains some distinct turbedite (from mudflows or slumps) layers; analyses of these layers are not included in the present paper. St8-GI consisted of beige to light gray distinctly laminated sediments, with laminations of some mm to several cm seen in optical light and highly distinctive in radiographic images, and also in grain size data and elemental composition (data not shown).

For molecular analyses, sediment cores were sampled by cutting small holes into the core liner and scraping off the surface with sterile spatulas before collecting sediment samples in duplicates with sterile 5 mL cut-off syringes. The sediment sub-samples were packed in Whirl-Pak bags and immediately frozen at −80°C. Porewater was extracted from the sediment cores as described previously (Glombitza et al., 2015). Data from Rumohr cores shown in Glombitza et al. (2015) were not included in the present paper because corresponding samples for molecular analyses were not available. Rumohr core devices such as those used in Glombitza et al. (2015) are small diameter gravity corers that enable more undisturbed, high quality sampling of upper sediment layers, although less material is recovered due to the smaller diameters.

Analytical Methods

Sulfate, hydrogen sulfide, and dissolved inorganic carbon (DIC) concentration measurements, as well as sulfate reduction rate (SRR) measurements were previously obtained and described from Gravity cores described in Glombitza et al. (2015). In addition, ferrous iron, total organic carbon (TOC), and total nitrogen (TN) were analyzed. For analyses of ferrous iron (Fe2+) concentrations, porewater was amended with 0.5 M HCl (1:1, v/v) and stored at 4°C. One mL of ferrozine solution (1 g L–1 in 50 mM HEPES-buffer at pH 7) was subsequently added to 20 μL of acid-preserved porewater, producing a magenta color reaction (Sørensen, 1982; Phillips and Lovley, 1987). The absorbance was measured at 562 nm (Stookey, 1970) with a spectrophotometer (FLUOstar Omega, BMG Labtech GMBH). To determine TOC and TN, sediment samples were treated with sulfurous acid (5–6% w/w) to remove inorganic carbon. Once dried, 50 mg of acidified sediment were packed into cleaned tin cups and burned in the elemental analyzer (FLASH EA 1112 series, Thermo Scientific). TOC and TN concentrations were calculated from a standard curve with wheat flour, which contains 43.37% of carbon and 2.31% of nitrogen. The C/N ratio was calculated as the molar ratio of TOC to total nitrogen. For determination of methane concentrations, 2 cm3 sediment was transferred to 20 mL GC vials containing 2.5 mL saturated NaCl with excess crystalline salt. The bottles were vigorously shaken to release methane into the headspace of the GC vials and stored upside down at −20°C until further analysis. Methane concentrations in the vial headspace were subsequently measured by gas chromatography (SRI 310C GC, SRI Instruments Europe GmbH) with a flame ionization detector. Ammonium concentrations were determined from 2 mL of porewater. Concentrations were analyzed spectrophotometrically as previously described (Dansk-Standard, 1975; Bower and Holm-Hansen, 1980).

Sediment age profiles for cores from St3-NGI, St6-NGI and St5-GI were based on radiocarbon dating and age modeling. Materials and depths were chosen for dating based on availability (Supplementary Table S1). For the 14C age determination, calcareous mollusk shells, benthic foraminifera, and seaweed samples were collected from cores and the 14C concentrations were determined by Accelerator Mass Spectrometry at the AMS 14C Dating Center of Aarhus University. The 14C ages were calibrated using the Marine13 radiocarbon calibration curve (Reimer et al., 2013) with a reservoir correction of ΔR = 140 ± 35 years. Age-depth models for cores from St3-NGI and 6 were calculated using the software Oxcal v4.2 (Ramsey, 2008). The presence of rapidly deposited turbidite sequences in the core from St5-GI made a detailed age model difficult.

In the core retrieved from St8-GI no material usable for 14C dating was found, and thus radiocarbon dating was not possible. Instead, the sediment age was estimated from 210Pb and 226Ra measurements of freeze-dried, homogenized sediment samples. The water content of each sample was determined before and after freeze-drying, and results were reported on a dry-weight basis and salt-corrected for bottom-water salinity at the time of sampling. Core sections were analyzed by gamma spectroscopy using a CANBERRA® Broad Energy Germanium with a P-type detector (model BE3830). Detector efficiency and self-absorption were corrected by counting reference material from the International Atomic Energy Agency (IAEA) within the same geometry. The reproducibility errors, determined by counting the same sample four times, were 5.7 and 3.8%, for 210Pb and 226Ra, respectively. 137Cs was not detected in the 25 analyzed sediment samples. However, due to the very long core, we were not able to analyze all intervals from the entire core and therefore a brief event such as the Chernobyl accident may easily have been missed. 210Pb was also determined from its 210Po daughter isotope using alpha spectroscopy (n = 23 sediment sections). The samples were prepared as previously described (Flynn, 1968). A 209Po tracer, calibrated against a 210Po NIST standard (Isotope Product Laboratories) was employed for quantitation. The reproducibility error was less than 1%.

Sedimentation rates that were used for the age model presented in this study were estimated from the least-squares fit to the natural log of the excess 210Pb (210Pbex) in the core and the output of a one-dimensional two-layer advection diffusion model that accounts for both bioturbation and compaction with depth (Lavelle et al., 1985; Kuzyk et al., 2015; Supplementary Figure S1). Although surface layer mixing that can occur in organic-rich shelf sediments may interfere with tracer-based age model reconstructions (Kuzyk et al., 2013), this was not evident in St-8-GI for the following reasons: the sediments were low in organics and 210Pbex decreased exponentially with depth; and the sediments were highly laminated, which was interpreted as annual deposition layers, with shifts in grain size caused by varying melting rates in summer vs. winter.

To estimate the loss of surface sediment by gravity coring and to correct the age-depth model, we compared the porewater profiles of ammonium, DIC, and the carbon isotope ratio 13C/12C of DIC (data not shown) between Rumohr cores collected during the same cruise and the gravity cores presented in this study. The upper 18 cm of the gravity core retrieved at St3-NGI and the upper 10 cm of the core retrieved at St6-NGI were missing. The depths of potential surface sediment loss at the two GI stations were corrected as the mean values of the two first cores (14 cm) since the Rumohr core casts at the GI stations were not successful. Finally, sediment age was recalculated as “actual age,” i.e., the surface of the seafloor was considered as 0 year old at the time of sampling, based on the age models and the above mentioned depth correction.

DNA Extraction and Preparation of 16S rRNA Gene and dsrB Amplicon Libraries

Approximately 0.5–1 g of sediment was used for DNA extraction according to a previously established protocol (Kjeldsen et al., 2007). Per station, 8–10 samples from different depths, corresponding to approximately 2 samples per meter core were selected for further analyses, with highest resolution at the top of the cores (Supplementary Table S2). Barcoded 16S rRNA gene amplicons were produced with a two-step PCR barcoding approach (Herbold et al., 2015), using the general bacterial and archaeal primers U519F (5′-CAGCMGCCGCGGTAATWC-3′) and 802R (5′-TACNVGGGTATCTAATCC-3′) for initial amplification (Klindworth et al., 2013). These primers were modified with a 16 bp head sequence as described previously (Herbold et al., 2015). The first round of amplification was performed in triplicates with 12.5 μL per reaction volume. The reaction mix contained 1 × Taq buffer (Thermo Scientific), 0.2 mM dNTP mix (Thermo Scientific), 2 mM MgCl2, 0.25 U Taq polymerase (Thermo Scientific), 0.2 μM of each primer, and approximately 1−10 ng DNA. The PCR started with a denaturation at 95°C for 3 min, followed by 30 cycles of 95°C for 30 s, 48°C for 30 s and 72°C for 30 s, and a final elongation at 72°C for 2 min. The subsequent barcoding PCR round (50 μL total volume) was performed with 1 × Taq buffer (Thermo Scientific), 0.2 mM dNTP mix (Thermo Scientific), 2 mM MgCl2, 1 U Taq polymerase (Thermo Scientific), 0.2 μM of each primer, and 2 μL of the pooled triplicate PCR products from the first PCR reaction. The thermal cycling program consisted of an initial denaturation at 95°C for 3 min, 12 cycles of 95°C for 30 s, 52°C for 30 s, and 72°C for 30 s, followed by a final elongation at 72°C for 2 min. The dsrB amplicons were produced with the equimolar primer mixes of DSR1762F and DSR2107R variants (Supplementary Table S3) according to an established protocol (Pelikan et al., 2016). Barcoded amplicons were mixed and further prepared for multiplexed, paired-end MiSeq sequencing (Herbold et al., 2015). Sequence data-sets are available in the NCBI Sequence Read Archive under study accession number PRJNA546002.

Sequence Data Processing

16S rRNA gene and dsrB amplicon raw reads were demultiplexed, filtered and clustered as described previously (Herbold et al., 2015; Pelikan et al., 2016) using fastq-join (Aronesty, 2013) to merge reads and UPARSE version 8.1.1861 (Edgar, 2013) to generate OTUs. Phylum/class-level classification of 16S rRNA-OTUs was performed with the Ribosomal Database Project naïve Bayesian classifier in MOTHUR (Wang et al., 2007; Schloss et al., 2009), using the SILVA database v.128 (Quast et al., 2013) as a reference. The dsrB-OTUs were classified by phylogenetic placement of representative sequences into a DsrAB reference tree (Müller et al., 2015) that was updated with novel sequences from diverse candidate phyla (Parks et al., 2017; Anantharaman et al., 2018; Hausmann et al., 2018). This DsrAB reference tree was constructed by de-replicating novel DsrA and DsrB sequences with less than 100% similarity to any DsrAB sequence in the original reference database and aligning them to the reference alignments of DsrA and DsrB (Müller et al., 2015) using MAFFT (Katoh et al., 2002). The combined DsrA and DsrB alignments were then concatenated and sequences with a total length of less than 500 amino acids were removed. The concatenated DsrAB alignment was clustered at 70% sequence identity with usearch (Edgar, 2010), and alignment positions were kept if they were conserved in at least 10% of all sequences in the 70%-clustered alignment (56 sequences). The unclustered DsrAB alignment (2985 sequences) was then filtered for conserved alignment positions using seqmagick[1] and was used to generate a maximum likelihood tree with FastTree (Price et al., 2010). 16S rRNA gene and dsrB OTU tables were processed in R using native functions (R Core Team, 2015) and the R software package phyloseq (McMurdie and Holmes, 2013). OTU counts were rarefied, i.e., sub-sampled at the smallest library size (dsrB: 1521; 16S rRNA gene: 4517) and transformed into relative abundances for all further analyses, except for network analyses, which were performed with the unrarefied OTU count matrices (Friedman and Alm, 2012). The DsrAB-alignment, -tree and -taxonomy files used in this study are provided as Supplementary Files S1–S3.

Statistical Analyses

Shannon alpha diversity was calculated using R (R Core Team, 2015). Beta diversity analyses were performed with the R software package vegan (Oksanen et al., 2017), including calculations of Bray-Curtis distances with the function ‘vegdist()’ and non-metric multidimensional scaling ordination analysis with the function “nMDS().” Environmental variables were tested for effects on the overall community composition by Mantel tests using the native R function “mantel().” Obtained p-values were corrected for multiple testing with the native R function “p-adjust()” using the Benjamini–Hochberg correction method. Correlations of individual OTUs with environmental variables were calculated with the native R function cor() using the Spearman correlation coefficient. p-Values were generated by permuting the values of each environmental variable followed by correlation of individual OTU abundances with the permuted environmental variable. This process was repeated 1000 times and the obtained p-values were corrected as described above.

Correlation network analyses were performed separately for 16S rRNA- and dsrB-OTUs to highlight potential synergistic interactions between microbial community members (Weiss et al., 2016). Species co-occurrence networks were calculated using SparCC (Friedman and Alm, 2012) based on count matrices of all OTUs with >10 reads in at least 7 out of 34 samples. Use of more abundant and prevalent OTUs increases sensitivity of the network analyses (Berry and Widder, 2014). p-Values were generated as described above. Only positive OTU correlations >0.5 were considered. Attributes of individual OTUs, i.e., sampling station and sediment age at which the OTU was found at the highest relative abundance, were assigned to OTUs in R and networks were visualized in Cytoscape (Shannon et al., 2003). Significant OTU clusters, i.e., significantly more interactions between OTUs within the community cluster than with OTUs outside the community cluster, were defined by Mann–Whitney U-tests using the Cytoscape plugin “clusterONE” (Nepusz et al., 2012).

Phylogenetic Analysis

Representative sequences of 16S rRNA-OTUs were aligned with the SINA aligner (Pruesse et al., 2012) using the SILVA database v.128 (Quast et al., 2013) as a reference. Sequences that were closely related to 16S rRNA-OTUs were extracted from the SILVA database and used to construct a reference tree with FastTree (Price et al., 2010). Subsequently, 16S rRNA-OTU sequences were placed into the reference tree using the EPA algorithm (Berger et al., 2011) in RAxML (Stamatakis, 2014). The placement trees of 16S rRNA-OTUs and dsrB-OTUs (utilized for dsrB-OTU classification) were visualized with iTOL (Letunic and Bork, 2007).

Results

Depth Profiles of Sediment Age and Porewater Chemistry Differ Substantially Between Non-glacier-influenced and Glacier-Influenced Sediments

A goal of the present study was to identify environmental factors (biogeochemical data is partially described in Glombitza et al., 2015) that shape the microbial community compositions and interactions in NGI and GI sediments. St3-NGI and St6-NGI are located on the open shelf and within the Godthåbsfjord, respectively, and were both characterized by a strong gradient of sediment age due to comparably low sedimentation rates, with maximum ages of the gravity cores close to 10,000 years (Figure 2). Furthermore, NGI sediments had high TOC and TN concentrations, as well as low C/N ratios down to 500 cm sediment depth. SRR decreased with depth at both stations and sulfate became depleted in the bottom of the core from St3-NGI, but not in the core from St6-NGI. Important to note, SRR were determined in greater resolution from Rumohr cores in Glombitza et al. (2015), whereby high SRR were especially more evident in the upper 25 cmbsf of stations St3-NGI and St6-NGI. For instance, SRR reached 35 nmol SO42– cm–3 d–1 in the upper sediments of St3-NGI (Glombitza et al., 2015). At St3-NGI, hydrogen sulfide concentrations gradually increased with depth and decreased again below a depth of 400 cmbsf, coinciding with the appearance of methane in the porewater. At St6-NGI, hydrogen sulfide was present in lower amounts and methane did not accumulate at any depth.

FIGURE 2

Physicochemical sediment properties in non-glacier-influenced (NGI) and glacier-influenced (GI) sediment cores. Colors indicate the sampling station. St3-NGI and St6-NGI, Non-glacier-influenced stations 3 and 6. St5-GI and St8-GI, Glacier-influenced stations 5 and 8. Note that the scales are different for each physicochemical parameter. Data on SRR and DIC, as well as, sulfate and sulfide concentrations were taken from Glombitza et al. (2015).

In comparison, St5-GI and St8-GI were both characterized by high sedimentation rates as indicated by the low sediment ages, i.e., around 200 years at the bottom of the core at St8-GI and around 500 years at the bottom of the core at St5-GI (Figure 2). GI stations had low TOC and TN concentrations, however, total carbon burial flux into GI stations was higher. Furthermore, GI stations had high C/N ratios of up to 48. The porewater contained dissolved Fe2+ in the upper 250 cmbsf at both GI stations. SRR remained detectable throughout most of the sediment depths in the GI stations, still reaching up to 1 nmol SO42– cm–3 d–1 at depths around 400 cmbsf at both GI sites. Despite substantial SRR, hydrogen sulfide did not accumulate any depths in GI sediments.

Glacier Runoff and Sediment Age Are Strong Determinants of Microbial Community Composition

In total, 6755 16S rRNA-OTUs and 1094 dsrB-OTUs were obtained by amplicon sequencing. NGI and GI sediments clearly differed in 16S rRNA- and dsrB-OTU compositions and these differences increased with depth at NGI stations (Supplementary Figures S2A,B). Mantel correlations between beta-diversity and environmental factors revealed that the 16S rRNA- and dsrB-OTU compositions were mostly impacted by sediment age and C/N ratio of organic matter (Table 2). Further differences in 16S rRNA gene and dsrB-OTU compositions were explained to a lesser extent by sediment structure (i.e., density and porosity), TOC, TN, SRR, and Fe2+ concentrations (the latter two parameters were only significant for the 16S rRNA gene community) (Table 2). Gradually increasing sediment age with depth at NGI stations (Figure 2) was associated with gradual changes in 16S rRNA gene and dsrB beta-diversity with depth (Supplementary Figures S2A,B). 16S rRNA and dsrB alpha-diversity at the NGI stations gradually decreased with depth (Supplementary Figures S2C,D). In contrast, alpha-diversity in GI stations remained rather high throughout the cores (Supplementary Figures S2C,D). Relative abundance patterns of most phyla/classes and DsrAB families at St5-GI and St8-GI did not follow a gradual change in compositions with increasing sediment depth like at St3-NGI and St6-NGI, but remained rather constant, with some fluctuations among taxa (Figure 3). Alpha-, Delta-, and Gamma-proteobacteria, Campylobacterota and notably Cyanobacteria were overall more abundant at GI sediments as compared to NGI sediments (Figure 3A). Furthermore, the dsrB-containing community in GI sediments had higher relative abundances of Desulfobacteraceae, Desulfobulbaceae, uncultured DsrAB family-level lineages 4, 7, and 9 and unclassified DsrAB sequences from the Firmicutes group and the Nitrospira supercluster (Figure 3B). At NGI stations, several phyla/classes, i.e., Acidobacteria, Bacteroidetes, Deltaproteobacteria, Dependentiae (TM6), Gammaproteobacteria, Omnitrophica (OP3), Planctomycetes, and Woesearchaeota (DHVEG-6) decreased in relative abundance with depth, particularly at St3-NGI (Figure 3A). The phyla Atribacter (JS1), Aerophobetes (BHI80-139), Aminicenantes (OP8), Alphaproteobacteria, Betaproteobacteria, and Chloroflexi increased in relative 16S rRNA gene abundances with depth at both NGI stations (Figure 3A). We note that although some Betaproteobacteria are known to be present in marine environments (Orcutt et al., 2011), some sequences detected here, e.g., Herbaspirillum or Ralstonia spp. commonly represent laboratory-derived sequences (Salter et al., 2014). Their low-level presence should therefore be treated with caution, although they did not affect the overall conclusions of this study. The relative abundances of the following dsrB-containing groups decreased with depth at NGI stations: Desulfobacteraceae, Syntrophobacteraceae, the uncultured family-level lineages 7 and 9, and uncultured bacteria within the Environmental supercluster 1 (Figure 3B). In contrast, representatives of the uncultured family-level lineage 3, as well as uncultured bacteria within the Deltaproteobacteria supercluster and the Firmicutes group, increased in relative abundances with depth at the NGI stations.

TABLE 2

Mantel correlations between 16S rRNA gene and dsrB community compositions and physicochemical parameters.

	16S rRNA gene		dsrB
	Mantel statistic	p-Value	Mantel statistic	p-Value
Sediment age [actual years]	0.380	0.003	0.538	0.002
Sulfate [mmol L–1]	0.179	0.016	0.390	0.002
DIC [mmol L–1]	0.108	0.047	0.207	0.010
H2S [μmol L–1]	0.227	0.003	0.418	0.002
Fe(II) [μmol L–1]	0.262	0.003	–	–
Density [g cc–1]	0.234	0.003	0.349	0.002
Porosity	0.213	0.003	0.364	0.002
TOC [μmol g dw–1]	0.228	0.005	0.329	0.002
C/N ratio	0.471	0.003	0.478	0.002
SRR [nmol Sulfate cm–3 d–1]	0.260	0.003	–	–
Methane [μmol L–1]	–	–	0.310	0.002

FIGURE 3

Microbial community composition in non-glacier-influenced (NGI) and glacier-influenced (GI) sediment cores. Changes in 16S rRNA phylum/class (A) and DsrAB-family (B) relative abundances with sediment depth are shown. Only phyla/classes and DsrAB-families with a relative abundance greater than 1% are shown. St3-NGI and St6-NGI, Non-glacier-influenced stations 3 and 6. St5-GI and St8-GI, Glacier-influenced stations 5 and 8. DS, Deltaproteobacteria supercluster; ES1, Environmental supercluster 1; FG, Firmicutes group; NS, Nitrospira supercluster.

OTUs That Were Positively Correlated to SRR Have High Inter-Species Connectivity in Young Sediments

Correlation network analyses of 16S rRNA- and dsrB-OTUs revealed two nearly separated 16S rRNA-OTU clusters and two completely separated dsrB-OTU clusters, which were structured along a sediment age gradient (Figures 4A,B). One cluster included OTUs that were most abundant in old NGI sediments and the other one included OTUs that were most abundant in young GI and NGI sediments. These separate network clusters in “young” and “old” sediments largely overlapped with regions of particularly high inter-species OTU correlations. These “young” and “old” sediment OTU clusters were separated at a sediment age of around 300–400 years (Figures 4C,D), which corresponds to a sediment depth of 30–40 and 20−30 cmbsf at NGI stations 3 and 6, respectively.

FIGURE 4

Co-occurrence of abundant OTUs across non-glacier-influenced (NGI) and glacier-influenced (GI) sediment cores. Inter-species correlations are indicated for 16S rRNA gene (A) and dsrB (B) OTUs. Only edges with p ≤ 0.01 and R2 ≥ 0.5 are shown. OTUs are colored and shaped according to the approximated sediment age and sampling station at which they were found at the highest relative abundance, respectively. The orange and green border color of OTUs indicates significant correlations to sulfate reduction rates and sediment age, respectively. OTUs that are connected by black and purple edges formed significant community clusters in old and young sediments, respectively. St3-NGI and St6-NGI, Non-glacier-influenced stations 3 and 6. St5-GI and St8-GI, Glacier-influenced stations 5 and 8. (C,D) Age of the sediment layer at which individual OTUs from (A,B) were found at the highest relative abundance across all samples and stations. Each dot represents an OTU. Black and purple background colors indicate the affiliation to significant community clusters determined for the inter-species correlation networks, respectively.

Individual OTUs in the 16S rRNA- and dsrB-OTU networks were additionally subjected to correlation analyses with environmental parameters (Supplementary Figures S3, S4). The majority of OTUs that constituted the “young” and “old” sediment clusters were positively correlated with SRR and sediment age, respectively (Figures 4A,B). Most OTUs that correlated positively with SRR also correlated negatively with sediment age and vice versa (Supplementary Figures S3, S4). OTUs that positively correlated to SRR showed distinct distributions in relative abundances, i.e., highest abundances in the surface sediments of NGI stations and mostly ubiquitous distributions throughout the whole core at GI stations (Figure 5). Many of these 16S rRNA-OTUs (n = 10) and most of the dsrB-OTUs belonged to the family Desulfobacteraceae (Supplementary Figures S5, S6). Other SRR-correlated 16S rRNA-OTUs belonged to the families Desulfobulbaceae, Desulfarculaceae, and Syntrophobacteraceae and to the phyla/classes, Acidobacteria, Actinobacteria, Alphaproteobacteria, Bacteroidetes, Ignavibacteria, Gammaproteobacteria, Planctomycetes, and Woesarchaeota (Supplementary Figure S5). Besides the prevalence of Desulfobacteraceae, dsrB-OTUs positively correlated to SRR were affiliated with the uncultured family-level lineages 7 and 9 (Supplementary Figure S6). The dsrB-OTUs 40 and 41 were also affiliated with the Environmental Supercluster 1. OTU 40 belongs to a sequence cluster of uncultured bacteria, and is related to the metagenome-derived genome RBG_13_60_13 (accession number GCA_001796685.1) of a Chloroflexi bacterium (Supplementary Figure S6).

FIGURE 5

Relative abundances of 16S rRNA- and dsrB-OTUs with significant correlations to sulfate reduction rates in non-glacier-influenced (NGI) and glacier-influenced (GI) sediment cores. The column annotation indicates the approximate sampling depth in centimeters below seafloor (cmbsf). The color range from blue to red indicates the relative abundance of OTUs. Phyla/classes and DsrAB-families that were represented by more than one OTU are indicated in bold. St3-NGI and St6-NGI, Non-glacier-influenced stations 3 and 6. St5-GI and St8-GI, Glacier-influenced stations 5 and 8.

Operational taxonomic units positively correlated with sediment age were affiliated with diverse taxa (Figure 6). The phyla/classes that were represented by at least two age-correlated 16S rRNA-OTUs were Aerophobetes (BHI80-13), Alphaproteobacteria, Aminicenantes (OP8), Atribacter (JS1), Chloroflexi, Deltaproteobacteria, Euryarchaeota (Marine Benthic Group D), Omnitrophica (OP3), and Planctomycetes (Supplementary Figure S5). Three dsrB-OTUs that positively correlated with sediment age were affiliated to the family Desulfobacteraceae (Supplementary Figure S6). The dsrB-OTU 16 belonged to a group of uncultured bacteria in the Deltaproteobacteria supercluster, OTU 58 could only be assigned to the Firmicutes group, and OTU 1 was affiliated with the family Syntrophobacteraceae.

FIGURE 6

Relative abundances of 16S rRNA- and dsrB-OTUs with significant correlations to sediment age in non-glacier-influenced (NGI) and glacier-influenced (GI) sediment cores. The column annotation indicates the approximate sampling depth in centimeters below seafloor (cmbsf). The color range from blue to red indicates the relative abundance of OTUs. Phyla/classes and DsrAB-families that were represented by more than one OTU are indicated in bold. St3-NGI and St6-NGI, Non-glacier-influenced stations 3 and 6. St5-GI and St8-GI, Glacier-influenced stations 5 and 8.

Discussion

Glacier-Runoff Affects Age-Depth Relationships and Microbial Community Assembly in Marine Sediments

Rates and amounts of glacial inputs into sedimentary environments of fjords have a major impact on their biogeochemistry (Glombitza et al., 2015). Here, we have compared microbial community structure between NGI and GI sediments, and highlighted the environmental factors that underlie the observed differences in community assembly with sediment depth. Sulfate is the key terminal electron acceptor in marine sediments of the Godthåbsfjord (Sørensen et al., 2015), but the depth distributions of SRR and microbial community structures differed between NGI and GI sediments.

In the NGI sediments, steep gradients of SRR (Figure 2) presented here and previously (Glombitza et al., 2015), indicated that most of the labile organic matter deposited from marine primary production was mineralized near the seafloor surface, which is typically observed for marine shelf sediments (Flury et al., 2016). Nevertheless, methane and ammonium formation in deep sections of St3-NGI also indicated ongoing mineralization of organic matter with depth, with additional methane also likely diffusing upward from deeper layers and being consumed by anaerobic methane oxidation in the deeper sections of the core (Knab et al., 2008). Microbial communities in the NGI sediments became less diverse with depth and increasingly distinct from the surface communities (Supplementary Figure S2). Furthermore, total microbial cell counts reached nearly 109 cells cm–3 in the upper 10 cm of NGI sediments and decreased to under 108 cells cm–3 in the deepest sections (Glombitza et al., 2015). Such shifts in community composition and cell abundances with depth can be attributed to the progressing geochemical stratification of the sediment and decreasing available energy with increasing sediment age (Petro et al., 2017), and was also evidenced by lower TOC in deeper sections of St3-NGI.

In contrast, the two GI sediment cores were characterized by higher sedimentation rates, low porosity (St8-GI), young ages, low TOC concentrations, low TN, and high C/N ratios (Figure 2); the latter being attributed to influx of terrestrial organic matter (Meyers, 1994; Goñi et al., 2013; Wehrmann et al., 2014). Terrestrial organic matter of particulate phase transported by the glaciers is mostly old, diagenetically altered, and likely poorly available for microbial degradation (Wehrmann et al., 2014). Therefore, the strong correlation of C/N ratio differences with the microbial community structures might mainly reflect the high rate of sedimentation in GI sediments. The on average higher SRR throughout the GI sediment cores were possibly sustained by low amounts of reactive organic matter that was deposited from algal blooms in nutrient-rich waters of glaciated fjords (Bourgeois et al., 2016). Ongoing microbial activity with depth in GI sediments was reflected by total cell counts, which remained above 108 cells cm–3 throughout all depths of GI sediments (Glombitza et al., 2015). Glacial runoff also contains considerable amounts of iron and manganese (Bhatia et al., 2013; Wehrmann et al., 2014). Accumulation of dissolved Fe2+ suggested that Fe(III) reduction substantially contributed to organic matter mineralization in upper GI station sediments (Wehrmann et al., 2014). Lack of sulfide accumulation with depth indicated immediate re-oxidation and/or scavenging of sulfide produced from the sulfate reduction activity (Figure 2; Wehrmann et al., 2017). In agreement with previous studies (Park et al., 2011; Bourgeois et al., 2016; Buongiorno et al., 2019), differences in organic matter availability and electron acceptor concentrations are suggested to have had a major influence on the composition of the seafloor microbial community in glaciated fjords.

Identities and Potential Functional Interactions of Sulfur Cycling-Associated Taxa in “Young” NGI and GI Sediments

16S rRNA gene and dsrB correlation network analyses both revealed two main OTU interaction clusters (Figure 4). In one cluster most OTUs were positively correlated to SRR but not sediment age, while in the other cluster most OTUs were positively correlated to sediment age but not SRR. The relative abundances of 16S rRNA- and dsrB-OTUs that positively correlated with SRR were highest in “young” GI and NGI sediments with active sulfur cycling (up to about 400 years of age). The majority of these OTUs was affiliated with the family Desulfobacteraceae, as well as other bona fide deltaproteobacterial SRM taxa from marine sediments (Wasmund et al., 2017). In addition, several SRR-correlated dsrB-OTUs were affiliated with uncultured DsrAB lineages. Some of these lineages contain metagenome-derived sequences of uncultivated bacteria from the phyla Acidobacteria (DsrAB family-level lineage 9), Planctomycetes or Chloroflexi (Wasmund et al., 2016; Anantharaman et al., 2018). Interestingly, several 16S rRNA-OTUs that positively correlated with SRR were also affiliated with these phyla, supporting their putative involvement in sulfite/sulfate reduction or sulfur disproportionation in the Godthåbsfjord sediments.

Despite high rates of sulfate reduction, sulfide did not accumulate in the GI sediments (Figure 2) likely due to its reaction with metals resulting in its oxidation and/or precipitation (Wehrmann et al., 2014). Several 16S rRNA-OTUs that correlated positively with SRR were affiliated with taxa containing sulfur-oxidizing microorganisms such as the candidate genus PHOS-HE36 (phylum Ignavibacteriae) (Koenig et al., 2005), the Woeseiaceae/JTB255 sediment group (Gammaproteobacteria) (Dyksma et al., 2016; Mußmann et al., 2017), and the Rhodobacteriaceae (Alphaproteobacteria) (Lenk et al., 2012; Thrash et al., 2017; Supplementary Figure S5). In addition, OTU 30 affiliated with the sulfur-oxidizing genus Sulfurovum, was solely responsible for the high relative abundance of Campylobacterota at St5-GI. We hypothesize that high SRR and chemical oxidation of sulfide by metals may support significant populations of sulfur-oxidizing or sulfur-disproportionating taxa in deep GI sediments, although future work would be required to substantiate this hypothesis, e.g., detection of mRNA transcripts for sulfur-dissimilating enzymes.

We also identified numerous SRR-correlated OTUs that were affiliated with taxa that are not known to have sulfur-based energy metabolisms. While these OTUs could indeed represent unknown sulfur-cycling microorganisms, they may also be degraders of organic matter that fuel sulfate reduction with fermentation products. For instance, some of these 16S rRNA-OTUs were affiliated with BD2-2 (phylum Bacteroidetes), Phycisphaera (phylum Planctomycetes) or OM1 (phylum Actinobacteria). Representatives of these phyla hydrolyze and ferment organic polymers in marine sediments and consequently might have trophic associations with SRM (Schauer et al., 2011; Webster et al., 2011; Baker et al., 2015; Trembath-Reichert et al., 2016). These associations may therefore explain their cooccurrences with SRM detected here.

Trophic and metabolic interactions with primary-degraders could also help explain why the overall dsrB-harboring communities shifted in composition with depth and at different sites. These shifts could be interpreted as surprising because volatile fatty acids (VFAs) such as acetate and formate are thought to be the primary energy and carbon sources of many typical sulfate-reducers (Beulig et al., 2018), yet their concentrations remained relatively stable in concentrations among sites and depths (Glombitza et al., 2015). Additional contributing reasons for shifts among dsrB-harboring taxa could be due to differences in VFA uptake and activation energetics by different sulfate-reducers (Glombitza et al., 2015); differences in nutrient availability (März et al., 2018); switching to alternative sulfur dissimilatory pathways such as disproportionation (Riedinger et al., 2017) or to different metabolic modes such as fermentation (Marshall et al., 2019); or the possibility to use alternative organic and/or energy sources other than VFAs.

Assembly of the Deep Subsurface Microbial Biosphere in NGI Sediments

Operational taxonomic unit correlation analysis showed that OTU clusters and thus the microbial communities of young and old sediment zones in NGI sediments became “disconnected” at a sediment age of about 300–400 years (Figure 4), corresponding to a sediment depth of approximately 30 cmbsf at NGI stations. This is in line with observations that a considerable shift in microbial community structures occur in marine sediments below the zone of bioturbation, which was suggested to be the main site of assembly of the subsurface community (Jochum et al., 2017; Starnawski et al., 2017). Most OTUs that positively correlated with sediment ages in NGI sediments were affiliated with lineages known to harbor members that selectively persist from the surface into deep subsurface sediments, e.g., Chloroflexi, Aerophobetes, Atribacter (JS1), Aminicenantes (OP8), Alphaproteobacteria and Deltaproteobacteria (Supplementary Figure S5; Orcutt et al., 2011; Wang et al., 2016). Members of these taxa, such as the genus Desulfatiglans (Jochum et al., 2018), the deltaproteobacterial candidate lineage SEEP-SRB1 (Schreiber et al., 2010) or the euryarchaeal Marine Benthic Group D (Lloyd et al., 2013; Kaster et al., 2014; Wasmund et al., 2014; Nobu et al., 2016; Robbins et al., 2016; Wang et al., 2016) are postulated to have traits such as fermentation, sulfate reduction or acetogenesis to support the maintenance of basic cellular functions even under extreme energy-limited conditions in most subsurface sediment environments (Petro et al., 2017).

Conclusion

Coastal marine ecosystems in arctic and sub-arctic oceans are poised to be increasingly impacted by melting of glaciers caused by climate change. In this comparative study, we found that discharge from marine-terminating glaciers had a strong control over the depth-dependent microbial community assembly in sediments of a sub-arctic fjord. Increasing differences in the benthic community composition between GI and NGI sites with depth were largely explained by sediment age. High sedimentation rates at GI stations enabled a complex community of sulfur-cycling-associated microorganisms, including both putative SRM and sulfide oxidizers, to continuously thrive at high relative abundances from the surface deep into the subsurface. Similar communities of sulfur-cycling-associated microorganisms were also present in surface sediments at NGI stations. However, with increasing depth the surface communities were largely replaced by microorganisms that positively correlated with sediment age. Lower sedimentation rates at the NGI sites thus resulted in slow burial and highly selective survival of microorganisms adapted to the energy-limited subsurface (Petro et al., 2017). In summary, our results suggest that increased glacier runoff and the associated high sedimentation rates allow processes that are typically predominant in surface sediments such as sulfide oxidation and associated community members to be rapidly buried and maintained at high abundances in deep subsurface sediments.

Data Availability Statement

The datasets generated for this study can be found in the NCBI Sequence Read Archive under study accession number PRJNA546002.

Author Contributions

ClP, HR, KK, and AL designed the research. ClP generated and analyzed the sequencing data. M-SS, HR, and ChP collected the sediment cores. MJ, HR, and KK obtained the samples. MJ performed most of the biogeochemical analyses. KK performed the DNA extractions. M-SS, ChP, and ZK calculated the sediment ages. ClP, KW, and AL wrote the manuscript. All authors revised the manuscript.

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.