Nitrogen Metabolism Genes from Temperate Marine Sediments

In this study, we analysed metagenomes along with biogeochemical profiles from Skagerrak (SK) and Bothnian Bay (BB) sediments, to trace the prevailing nitrogen pathways. NO3− was present in the top 5 cm below the sediment-water interface at both sites. NH4+ increased with depth below 5 cm where it overlapped with the NO3− zone. Steady-state modelling of NO3− and NH4+ porewater profiles indicates zones of net nitrogen species transformations. Bacterial protease and hydratase genes appeared to make up the bulk of total ammonification genes. Genes involved in ammonia oxidation (amo, hao), denitrification (nir, nor), dissimilatory NO3− reduction to NH4+ (nfr and otr) and in both of the latter two pathways (nar, nap) were also present. Results show ammonia-oxidizing bacteria (AOB) and ammonia-oxidizing archaea (AOA) are similarly abundant in both sediments. Also, denitrification genes appeared more abundant than DNRA genes. 16S rRNA gene analysis showed that the relative abundance of the nitrifying group Nitrosopumilales and other groups involved in nitrification and denitrification (Nitrobacter, Nitrosomonas, Nitrospira, Nitrosococcus and Nitrosomonas) appeared less abundant in SK sediments compared to BB sediments. Beggiatoa and Thiothrix 16S rRNA genes were also present, suggesting chemolithoautotrophic NO3− reduction to NO2− or NH4+ as a possible pathway. Our results show the metabolic potential for ammonification, nitrification, DNRA and denitrification activities in North Sea and Baltic Sea sediments.


Introduction
Nitrogen is one of the essential elements of life and its cycle is driven mostly by microbial activities. Anthropogenic and physical processes also contribute to different sources of nitrogen (Delwiche 1970;Söderlund and Svensson 1976). Elucidating the pathways of nitrification and denitrification has been a topic of interest since the late nineteenth century when Winogradsky isolated Nitrosococcus Winogradsky (Winogradsky 1890). In marine sediments, nitrogen is cycled as part of the redox zonation in the suboxic zone (Froelich et al. 1979;Jørgensen 1989;Middelburg et al. 1993;Herbert 1999). Ammonium (NH 4 + ) is mainly liberated into the porewater by ammonification involving multiple steps of microbial breakdown of proteins, peptides and amino acids by proteolytic enzymes and deaminases (Herbert 1999 and references therein; Fig. 1). Also, the breakdown of nucleic acids can lead to release of urea (Therkildsen et al. 1996(Therkildsen et al. , 1997. Whether certain steps in the ammonification process in marine sediments dominate more than others remain poorly understood. Electronic supplementary material The online version of this article (doi:10.1007/s10126-017-9741-0) contains supplementary material, which is available to authorized users.
With respect to nitrification, whereby ammonia (NH 4 + / NH 3 ) is oxidized to NO 3 − , it is still unclear whether ammonia-oxidizing archaea (AOA) (Könneke et al. 2005;Wuchter et al. 2006;Francis et al. 2007) or ammoniaoxidizing bacteria (AOB) (Freitag et al. 2006) play a more important role in the first step of this pathway in marine sediments (Schleper and Nicol 2010). AOA (e.g. Nitrosopumilus) and AOB (e.g. Nitrosospira, Nitrosococcus, Nitrosomonas) mediate oxidation of NH 3 to nitrite (NO 2 − ) using the coppercontaining enzyme ammonia monooxygenase (AMO). While AOB use the protein hydroxylamine oxidoreductase (HAO) as an intermediate in the ammonia oxidation process, it is not yet clear what that intermediate protein is in the AOA process (Norton et al. 2001;Walker et al. 2010; Fig. 1). In the Baltic Sea, AOA have been observed to be abundant in the water column (Pitcher et al. 2011;Bale et al. 2013) and in sediments (Jørgensen 1989). It also remains debatable which groups of nitrate-oxidizing bacteria (NOB) including Nitrobacter, Nitrospira, and Nitrospina (Teske et al. 1994;Spieck and Bock 2005;Lücker et al. 2010) are likely to contribute most to oxidation of NO 2 − to NO 3 − in marine sediments.
Besides ammonification, dissimilatory nitrate (NO 3 − ) reduction to NH 4 + (DNRA) can also serve as a source of NH 4 + for primary producers and benthic organisms, yet it remains an understudied pathway (Giblin et al. 2013). In this step, respiratory and fermentative microorganisms reduce NO 3 − to NH 4 + (Fig. 1). In the first step of denitrification, NO 3 − is reduced to NO 2 − (Nar proteins). Once NO 2 − is produced, it can proceed to the next step of the denitrification pathway, or NO 2 − can be converted to NH 4 + (via NrfA, NirBD proteins) via the DNRA pathway (Fig. 1). A number of bacteria including Beggiatoa are capable of DNRA (Baggs and Phillipot 2011;Giblin et al. 2013).
In comparison to DNRA, the denitrification pathway is a more intensively studied pathway in marine sediments. Denitrification includes the conversion of NO 3 − to NO 2 − , NO 2 − to nitric oxide (NO) (Nir proteins), NO to nitrous oxide (N 2 O) (Nor proteins) and N 2 O to N 2 (Nos proteins) (Teske et al. 1994;Spieck and Bock 2005;Lücker et al. 2010;Pauleta et al. 2013 and references therein). Bacteria that are capable of completing different parts of the reaction are phylogenetically diverse but most are members of the Proteobacteria and are facultative aerobes (Shapleigh 2013 (Trubitsyn et al. 2013). In spite of our better knowledge about the denitrification pathway, the environmental factors that determine the balance between DNRA and denitrification are far from being understood. In some estuary sediments, DNRA has been found to be the dominant process influencing the fate of NO 3 − as opposed to denitrification (Soonmo and Gardner 2002;Giblin et al. 2010). Continuing to determine the abundances of microbial groups throughout the Baltic Sea and other areas could help elucidate their level of importance with respect to nitrogen cycling in these areas.
A number of studies have identified AOA and AOB in the water column and sediments using primers for 16S rRNA and the amoA gene via cloning and sequencing or by using lipid biomarker techniques (Francis et al. 2005(Francis et al. , 2007Dang et al. 2008;Pitcher et al. 2011, Bale et al. 2013). Similar techniques have been applied to study denitrifying communities in marine sediments Kerkhof 1998, 1999;Braker et al. 2000;Michotey et al. 2000). Furthermore, metagenomics has been extensively used as an initial step in inferring the functional and metabolic potential of microbial communities in environmental studies (Xie et al. 2011;Dini Anderote et al. 2012;Kimes et al. 2013;Scott et al. 2014). Thureborn et al. (2013) studied metagenomes associated with nitrification and  Cabello et al. (2009) denitrification in three anoxic sediment samples from Landsort Deep (Baltic Sea). However, a combined metagenomic and biogeochemical approach to study ammonification, nitrification, DNRA and denitrification pathways has not been reported for suboxic surface sediments of the North Sea and Baltic Sea.
To better understand which microorganisms and metabolic pathways could be contributing most to nitrification, ammonification, DNRA and denitrification, we studied N-cycling gene distributions in three sediment samples using an Illumina technique. Sediments came from the Skagerrak (SK), located at the North Sea-Baltic Sea, and from the Bothnian Bay (BB) located in the Baltic Sea. Computational analysis of the sequencing data using the Metagenomic RAST Analysis Server (MG-RAST) pipeline allowed for the analysis of an unprecedented number of sequences in these sediments providing a more complete representation of microbial communities and their genes. Porewater concentrations of SK and BB ferruginous sediments showed on-going NO 3 − reduction and NH 4 + production providing the opportunity to relate these pathways to biogeochemical zones in the surface sediments.

Study Site
Sediment cores of up to 38 cm in length (10 cm in diameter) were taken during the RV Meteor cruise No. M86-1 in November 2011 using a multi-corer device (Oktopus GmbH, Kiel, Germany). The two sampling locations considered in this study were BB site At4 (65°26.71′ N/23°17.92′ E) and SK Site Geo 2a (58°29.513′ N/9°35.855′ E) (Fig. 2). Samples were taken from a water depth of 75 m at site At4 and from a water depth of 554 m at site Geo 2a. After collection, cores for microbiological analyses from each site were frozen at −20°C on board the ship and later transferred to −80°C in the shorebased laboratory. For a more detailed description of the study sites refer to .

Porewater Analysis
Parallel cores were recovered for geochemical analyses and porewater was extracted from 23 depths from SK cores and 16 depths from BB cores directly after recovery by using Rhizon samplers (Rhizosphere Research Products B.V, Wageningen, The Netherlands) (Seeberg-Elverfeldt et al. 2005). Porewater samples for NO 3 − and NH 4 + were kept frozen until later analysis at the IOW (Leibniz Institute for Baltic Sea Research Warnemünde) and measured after the methods of Grasshoff et al. (1999) by using a continuous flow nutrient analyser (QuAAtro, Seal Analytical GmbH, Norderstedt, Germany).
Additional geochemical measurements of these cores were made and are discussed in .

Pore Water Modelling
Interpretation of interstitial water profiles of dissolved NO 3 − and NH 4 + at sites Geo 2a and At4 were carried out using the PROFILE (Berg et al. 1998) and the REC (Lettmann et al. 2012) models, considering steady-state conditions. Porewater profiles for the dissolved species used in the modelling can be found in . Local or non-local irrigation was neglected in the interpretation of porewater profiles. The diffusion coefficients in free solution at in situ salinity and temperature were calculated according to Boudreau (1997) and Schulz and Zabel (2006). The molecular diffusivity in the sediment was corrected for tortuosity according to Iversen and Jørgensen (1993), considering sediment porosities after Flemming and Delafontaine (2000). Both models calculate net and not gross process rates.

DNA Extraction
Frozen cores were sliced into 1 cm or 2 cm diameter discs with a band saw (K330S, Paul Kolbe GmbH, Elchingen, Germany) with a WIKUS blade (WIKUS DIAGRIT S Nr. 572 D254 VA, WIKUS-Sägenfabrik, Spagenburg, Germany). The blade was cleaned and sterilized with ethanol (70%) after cutting each slice. The sediment that was in contact with the blade was removed, and only the interior parts of the frozen disc were sectioned into aliquots for DNA extraction. DNA was extracted from two BB samples (BB 3-4 cm and BB 6-7 cm) and one SK sample (SK 6-8 cm). Extractions were made from~0.5-0.6 g of sediment based on the method of Lueders et al. (2004). Following extraction, samples were pooled and the nucleic acid pellet was dissolved in 200 μl of RNase/DNase free water. Humic substances and phenols from the crude extracts were removed with the Zymo PCR Inhibitor Removal kit (Zymo Research, Freiburg, Germany) following the manufacturer instructions. The total DNA was purified by digesting the RNA with 10 μl of 10 mg/ml of Roche RNase A (Sigma-Aldrich, Munich, Germany) for 5 min at room temperature and the Zymo Genomic DNA Clean and Concentrator kit (Zymo Research, Freiburg, Germany). The following modification was made to the manufacturer method: following the first elution through the column, the eluate was not discarded. Instead the eluate was added to a new column. Both columns were processed according to manufacturer instructions and total DNA eluted with water from both columns. Extracts were checked for nucleic acid concentration by NanoDrop Spectrophotometer ND-1000 (PeqLab-VWR International GmbH, Erlangen, Germany) and using the Invitrogen Quanti-iT PicoGreen dsDNA Assay Kit (Thermo Fisher Scientific, Darmstadt, Germany) following the manufacturer instructions.

Metagenome Sequencing
DNA shotgun libraries were generated using the Nextera DNA Library preparation kit following the manufacturer instructions (Illumina, San Diego, USA). The metagenomes were sequenced in a 112-bp paired end single indexed run with the Genome Analyzer IIx (Illumina, San Diego, USA).

Processing of Short Reads
Quality filtering of the metagenomic reads was performed with trimmomatic (0.25) (Bolger et al. 2014). This included the removal of Illumina adaptors and sequences that were shorter than 50 bp. Additionally, sliding window clipping was performed to remove read spaces that had an average quality below 15 (Phred score 33) within a 4-bp window. Samples were analysed using the Metagenomic Analysis Server (MG-RAST) (Meyer et al. 2008) using default settings. Sequences were normalized via transformation, standardization and multiple sample scaling as described in the MG-RAST manual (Wilke et al. 2015). Briefly, for normalization, raw abundances were transformed by a (log2 + 1) transformation, standardized and scaled from 0 to 1 (Meyer et al. 2008).
To compare 16S rRNA sequence abundances, the MG-RAST M5RNA database, which integrates SILVA, Greengenes and RDP databases, was used. An e value cutoff of 1e-05, minimum identity 60% cut-off and a minimum alignment cut-off of 15 were the default parameters used in the analysis. The best hit classification option was selected. Before analysis, BB forward and reverse sequences were clustered as one group and SK forward and reverse sequences as another group using the PCoA option. Following classification, the average of forward and reverse abundances was compared between samples. Normalized abundances of AOA were combined and compared to normalized abundances of AOB.
Protein-coding genes were annotated against the SEED Level (function) Subsystems of MG-RAST using the default parameters described above. The hierarchical classification option was selected. Ward Clustering and Bray Curtis distance options were selected for the heat map analysis. Normalization and standardization of forward and reverse sequences was accomplished using the R statistical computing system (version 3.2.1; R Core Team 2013), R package BmatR( version 0.9) (Braithwaite and Keegan 2013) and accessory apps built on the BmatR^package for MG-RAST (Keegan 2015). Following classification, the average of forward and reverse abundances were compared between samples.
A one-way ANOVA (alpha 0.05) and unpaired t tests (alpha 0.05) were used to determine if the abundances of proteincoding genes between and within samples were significant (StatPlus, Microsoft Excel 2011). When comparing between sites, BB samples were grouped together for the analysis. For with respect to the locations of Sites S1 through S9 of the transect studied by Canfield et al. (1993) the ANOVA analysis, the following gene abundances were compared: SK DNRA, SK denitrification, SK both (referring to genes belonging to both DNRA and denitrification pathways), BB DNRA, BB denitrification and BB both. To determine differences in gene abundances between sites in terms of ammonification, DNRA, denitrification and genes belonging to both DNRA and denitrification pathways, the following steps were carried out: (i) an F test was used to determine equal or unequal variance, (ii) differences in the number of genes being compared between sites was corrected by randomly sampling the BB gene abundance data to match the number of genes in the SK and (iii) a two-tailed distribution was used to determine significance. When comparing between BB34 and BB67 samples, BB AOA vs. BB AOB groups and SK AOA vs. SK AOB groups, step (ii) was excluded because the number of genes were similar between these samples or groups.

Characteristics of Skagerrak and Bothnian Bay Metagenomes
Following Illumina sequencing of two BB samples and one SK sample, approximately 4.9 Gb of sequencing data were generated consisting of 19-21 million (BB samples) and 4 million (SK sample) reads with an average read size of 112 bp. The metagenome data was deposited in the MG-RAST database (Table S1). To our knowledge, this is the largest metagenomic dataset to date of Baltic Sea sediments. An average of 700-800 k genes were annotated as non-rRNA proteins for BB samples and 155 k genes for the SK sample. For comparative purposes, reads were normalized to rule out the effect of inter-sample variations on the read abundances as described in the BMethods^section above.

Overview of Nitrogen Metabolism Genes Detected in Metagenomes
Protein-coding genes involved in the nitrogen cycle detected in our samples were not the most abundant of all genes in samples (Fig. 3). Furthermore, genes for ammonification, DNRA, nitrification and denitrification pathways were detected in metagenomes and are summarized in Fig. 4. One-way ANOVA analysis of DNRA, denitrification gene abundances and genes involved in both pathways revealed a significant difference between samples (P < 0.05).

Biogeochemistry in Surface Sediments
The vertical profiles of dissolved NO 3 − and NH 4 + observed at SK (Fig. 5a) and BB (Fig. 5b) indicate different zones of net production or consumption. They follow the expected porewater redox zonation for sediments (e.g. Froelich et al. 1979). Based on steady-state modelling (Berg et al. 1998;Lettmann et al. 2012), the gradients of measured dissolved NO 3 − and NH 4 + were interpreted quantitatively to estimate the zones of net (de-) nitrification and ammonium production. NO 3 − and NH 4 + models for SK (Fig. 5a) and for BB (Fig. 5b) yield similar rates and vertical zones of net transformations, with the major difference being the use of discrete zones in the PROFILE model (Berg et al. 1998), whereas continuous smooth rate changes are used in the REC model (Lettmann et al. 2012).
At the SK site, NO 3 − occurs below 0.5 cm and is present down to about 5 cm (Fig. 5a). Calculated rates of NO 3 − consumption from concentration gradients show a maximum consumption at 5 cm ( Fig. 5a) due to NO 3 − reduction. NH 4 + occurs below 2.5 cm within the zone of NO 3 − reduction and increases constantly below this depth (Fig. 5a). The rate calculations indicate net NH 4 + consumption occurs below the surface and a second consumption peak below~12 cm, and net production below~17 cm (Fig. 5a).
At the BB site, NO 3 − increases immediately below the surface and is present to 5 cm (Fig. 5b). Modelling results indicate NO 3 − production in the top~4 cm (Fig. 5b).
Nitrification at the sediment-water interface is the likely source of NO 3 − . Below 4 cm, NO 3 − is consumed by NO 3 − reduction. NH 4 + appears at 0.5 cm and increases strongly below 5 cm but only a minor increase occurs below this depth (Fig. 5b). Modelling results indicate NH 4 + production mainly between 5 and 10 cm, while its production is masked by NH 4 + consumption by nitrification in the top 4 cm (Fig. 5b).

Ammonification
The relative abundance of protein-coding genes involved in ammonification appeared abundant in SK and BB samples (Table 1). However, an unpaired t test analysis did not reveal any significant difference in ammonification gene abundances between samples (P > 0.05). Ammonification genes were composed of protease, hydratase, peptidase, urease and deaminase genes (Fig. 6).
The genes amo and the hao genes were present in all samples (≤0.02 normalized abundances). AMO and HAO are involved in the first step of nitrification in AOB; however, no hao homologue has been found in archaeal genomes (Stahl and de la Torre 2012). We also include additional supporting pyrosequencing results that show evidence for the high abundance of Thaumarchaeota in BB sediments relative to other archaea (Fig. S1).
Fermentative bacteria such as Clostridium and Bacillus spp. can use NO 2 − as an electron acceptor during fermentative growth, employing the use of cytoplasmic (NADH) nitrite reductase NirBD proteins (Cabello et al. 2009), thereby contributing to production of NH 4 + (Baggs and Phillipot 2011). 16S rRNA abundances showed they were more abundant in BB samples (Fig. 7b). Protein-coding genes for NirB (nitrite reductase [NAD(P)H] large subunit) and NirD were also present in all samples (Fig. 8a). Comparison of DNRA proteincoding gene abundances (≤0.014 normalized abundance) (Fig. 8a) to denitrification protein-coding gene abundances (≥0.01 normalized abundance) (Fig. 8b) within samples showed genes for DNRA were significantly less abundant (unpaired t test; P < 0.05) (Fig. 10a, b). Fig. 3 Heat map showing normalized protein-coding gene abundances predicted to be involved in various metabolic pathways detected in SK and BB metagenomes. The most abundant genes (red) and least abundant genes (grey). BB34 refers to sample 3-4 cmbsf (centimetre below seafloor), BB67 refers to sample 6-7 cmbsf and SK68 refers to sample 6-8 cmbsf. Values for samples were scaled from 0 (minimum value) to 1 (maximum value) using a uniform scaling method implemented in the MG-RAST pipeline. Forward and reverse abundances were averaged and a single value is reported per sample

Denitrification Pathway
16S rRNA genes of Thiothrix were present in all samples (Fig. 7b). At both sites, protein-coding genes for denitrification were found. In BB samples, the abundances of protein-coding genes involved in denitrification were similar between the two depths (unpaired t test; P > 0.05) (Fig. 8b). Similar abundances of protein-coding genes were also observed between SK and BB samples (unpaired t test; P > 0.05). Genes including nirK (copper-containing nitrite reductase), nirS (cytochrome cd1 nitrite reductase) (Helen et al. 2016) and nirV and nirN genes were detected in all samples (Fig. 8b). While most α, β and γ-Proteobacteria use NirK to reduce NO 2 − to NO, other bacteria use NirS as an alternative enzyme and NirN is a homologue of NirS (van Spanning 2011). Additionally, NirV is sometimes associated with NirK, and may incorporate copper into the redox centre of NirK (van Spanning 2011). Nor genes (norQD) were present (Fig. 8b) along with nitric oxide reductase (NOR) subunit C and B genes (EC 1.7.99.7). The catalytic NOR enzyme encoded by norCB can be co-expressed with accessory genes norQ, norD (although their function remains elusive). Additionally, nos genes (nosRDFYLX) and a nitrous oxide reductase gene (EC 1.7.99.6) were present in all samples (Fig. 8b). NosFY and D are ABC transporters that are usually linked to expression of NosZ, the main multi-copper enzyme involved in N 2 O reduction, and NosR is required for transcription of NosZ (Spiro 2012). Additionally, NosL, a membrane anchored copper protein, and NosX, a periplasmic flavoprotein, may serve as accessory proteins to NosZ (Spiro 2012). We found that with the exception of nirK, none of the other genes that were abundant in our DNRA or denitrification heat maps are known to occur in multiple copies. When the abundance of Fig. 4 Ammonification, DNRA, nitrification and denitrification pathways based on proteincoding genes detected in SK and BB samples. Enzymes (blue) and taxa (orange) that could potentially be involved in each pathway are written next to each arrow. Question marks indicate uncertainty in the organism's involvement in the pathway despite detection. Figure modified after Cabello et al. (2009) this gene was left out of the analysis, the results of the heat map remained the same.

Genes Involved in Denitrification and DNRA Pathways
Genes in the Nar and Nap families that are involved in the first step of the denitrification and DNRA pathways were detected in all samples but were not significantly different between samples (unpaired t test; P > 0.05) (Fig. 8c). The narGH genes, which encode the alpha and beta subunit of the Nar enzyme, were the most abundant genes detected out of all denitrification and DNRA protein-coding genes (Fig. 8c). The narI gene (gamma subunit), the chaperone encoding narJ (delta subunit) gene, nitrate response regulator protein (narQ/P) and nitrate/nitrite sensor protein (narX/L) were also detected but at lower abundances ( Fig. 8c) (Blasco et al. 1998;Cabello et al. 2009). Among the Nap family of genes, a napC (a membrane tetraheme c-type cytochrome), napA (NO 3 − reductase) and accessory proteins napG and H were detected. The NapC protein is an essential component of certain δ and γ-Proteobacteria (including Desulfovibrio spp. and MR-1) and donates electrons to NapA (Potter and Cole 1999;Brondijk et al. 2002;Simon 2002;Chen and Wang 2015). However, an alternative pathway in these interactions can also include NapG and H (Baggs and Phillipot 2011).

Geochemistry
Modelling results suggest ammonification and nitrification could be active processes in the BB (3-4 cm) sample. While ammonification occurs most likely throughout the SK core due to organic matter degradation, it is masked by even higher rates of nitrification within the top 5 cm (Fig. 5b). A double peak in NH 4 + consumption may result from non-steady-state conditions, most likely due to bioturbation. In addition, modelling results indicate denitrification and DNRA could take place in the SK (8-10 cm) and BB (6-7 cm) samples (Fig. 5a, c). Sediments with a high organic carbon input and nitrogen limitation are predicted to favour DNRA over denitrification based on previous studies (Giblin et al. 2013;Algar and Vallino 2014). Both the SK and BB receive high inputs of organic matter (sedimentation rate~4 mm year −1 SK;~1.1-1.6 mm year −1 BB; van Weering et al. 1987;Mattila et al. 2006). However, the quality of organic matter could also be an important factor. In the SK and BB, total organic carbon (TOC) is available throughout both sediments ) and could stimulate heterotrophic activity. However, a more refractory organic matter in the BB   could also slow down heterotrophic activity. Total nitrogen (TN) is also readily available in both sediments  resulting in a C/N ratio of~10 and 30 for SK and BB, respectively. Due to the high input of organic matter and nitrogen in both sites, denitrification could possibly dominate over DNRA.

Metagenome
When interpreting results derived from DNA in sediments, it is important to keep in mind certain physical processes that could influence the distribution of genes. One process is sedimentation, which occurs as new sediment layers are Urea_decomposition----Urease alpha subunit (EC 3.5.1.5) 17 Urease_subunits----Urease alpha subunit (EC 3.5.1.5) 18 Inositol_catabolism----Epi-inositol hydrolase (EC 3.7.1.-) 19 Predicted_carbohydrate_hydrolases----COG2152 predicted glycoside hydrolase 20 Acetyl-CoA_fermentation_to_Butyrate----3-hydroxybutyryl-CoA dehydratase (EC 4.2.1.55) 21 Acetyl-CoA_fermentation_to_Butyrate----Enoyl-CoA hydratase (EC 4.2.1.17) 22 Novel_non-oxidative_pathway_of_Uracil_catabolism----Urease alpha subunit (EC 3.5.1.5) 23 Proteasome_bacterial----ATP-dependent Clp protease ATP-binding subunit ClpX 24 Proteasome_bacterial----ATP-dependent Clp protease proteolytic subunit (EC 3.4.21.92) 25 Proteasome_bacterial----ATP-dependent hsl protease ATP-binding subunit HslU 26 Proteasome_bacterial No. refers to protein-coding gene number as shown in Fig. 6. Protein-coding genes with the highest abundances are underlined and in bold continuously deposited over time, thereby shifting the biogeochemical zones. Another process that could shift biogeochemical zones is storm events. During storms, sediments from the surface could be resuspended in the water column, bringing deeper zones to the surface. Thus, extra-cellular DNA or cellular DNA could end up preserved in a different biogeochemical layer. Another important sedimentary process to be aware of is bioturbation. Sedimentary layers could become homogenized by bioturbation, thereby influencing the distribution (and hence abundance) of genes in sediments (Laverock et al. 2014). Although the metagenome data do not directly indicate activity with respect to a particular organism or pathway, the presence of protein-coding genes in the sediment represents a metabolic potential. In this case, the presence of these genes represents a metabolic potential with respect to nitrogen metabolism. Since many dissimilatory bacteria are involved in amino acid and protein degradation thereby contributing to ammonification, perhaps this explains the higher abundance of ammonification genes over other types of nitrogen metabolism genes in both samples. In sediments of Landsort Deep (Baltic Sea), where DNRA and denitrification has been shown to be important, metagenome results also show that ammonification-coding genes are more abundant compared to genes involved in other nitrogen pathways (Thureborn et al. 2013). Moreover, bacterial proteases and hydratases could potentially contribute most to the release of NH 4 + from proteins and amino acids during ammonification in samples from both sites (Table 1).
In other marine sediment studies, the abundance of AOA and AOB varies and the factors that determine whether one is more abundant than the other are not clear and may depend on many variables including salinity (Caffrey et al. 2007;Mosier and Francis 2008;Santoro et al. 2008), NH 4 + availability (Smith et al. 2014) and spatial and temporal variations (Beman et al. 2012;Smith et al. 2015). From the 16S rRNA metagenomic results, it appears that mostly Thaumarchaeota and AOB could contribute to nitrification at both sites (Fig. 7a). Thaumarchaeota have been found to be abundant (Thureborn et al. 2013) and to play an important role in ammonia oxidation at other locations in the Baltic (Labrenz et al. 2010;Feike et al. 2012). Furthermore, DNA pyrosequencing results of these samples indicate Thaumarchaeota/Nitrosopumilus are abundant relative to other archaea and are the major archaeal type in BB sediments (Fig. S1). While a diversity of archaea is present in SK sediments, the BB appears to be dominated by Thaumarchaeota . Based on metagenomic comparisons between AOA and AOB, however, it remains unclear whether AOA or AOB are more important with respect to NH 3 -oxidation in BB and SK sediments. Fig. 6 Heat map showing normalized protein-coding gene abundances predicted to be involved in ammonification metabolism detected in SK and BB metagenomes. The most abundant genes (red) and least abundant genes (grey). BB34 refers to sample 3-4 cmbsf (centimetre below seafloor), BB67 refers to sample 6-7 cmbsf and SK68 refers to sample 6-8 cmbsf). Numbers on the y-axis correspond to the genes listed in Table 1. Values for samples were scaled from 0 (minimum value) to 1 (maximum value) using a uniform scaling method implemented in the MG-RAST pipeline. Forward and reverse abundances were averaged and a single value is reported per sample In terms of denitrification, Thiothrix could have access to both thiosulphates and NO 3 − in SK and BB sediments, allowing it to carry out thiosulphate oxidation coupled to NO 3 − reduction (Meyer et al. 2007). Porewater results of BB samples showed a decrease in sulphate concentration, which overlapped with the zone of NO 3 − reduction but no sulphide was present . Thiothrix (including the isolate known to carry out this chemolithoautotrophic reaction) can live in freshwater and marine habitats (Trubitsyn et al. 2013), possibly explaining the presence of Thiothrix 16S rRNA genes in SK and BB samples. The presence of protein-coding genes related to sulphur oxidation and denitrification in both samples supports the idea that one way in which denitrification could occur in BB (6-7 cm) and SK (8-10 cm) samples is via a chemolithoautotrophic pathway.
Other organisms such as Nitrosomonas and Nitrobacter, which can carry out only certain steps of the denitrification process, could potentially be involved in reduction of NO 3 − to NO 2 − in BB (6-7 cm) and SK (8-10 cm) using nar genes. Nitrosomonas could also contribute to NO 3 − reduction to N 2 using nap, nir and nor genes in these same samples (Fig. 4). Some marine and freshwater Beggiatoa have the ability to denitrify (Sweerts et al. 1990;Muβmann et al. 2007), and DNRA is well documented in Beggiatoa  and references therein). Beggiatoa have been shown to be capable of surviving independently from external sources of sulphur and NO 3 − for up to 2 weeks in laboratory experiments . Protein-coding genes detected in all samples such as otr ( Fig. 8b; MacGregor et al. 2013) and sulphide oxidation (soxADZBX) genes (Fig. 9)  For similar reasons, it is unclear if Clostridium and Bacillus could be involved in DNRA in these sediments. Although we could not infer whether the above taxa could potentially be involved in fermentative DNRA, the presence of small and large NADH nitrite reductase subunit genes (Fig. 8b) Fig. 8 Heat maps showing normalized putative protein-coding gene abundances predicted to be involved in a DNRA, b denitrification or c both, detected in SK and BB metagenomes. BB34 refers to sample 3-4 cmbsf (centimetre below sea-floor), BB67 refers to sample 6-7 cmbsf and SK68 refers to sample 6-8 cmbsf. Values for samples were scaled from 0 (minimum value) to 1 (maximum value) using a uniform scaling method implemented in the MG-RAST pipeline. Forward and reverse abundances were averaged and a single value is reported per sample the idea that DNRA could also occur via heterotrophic fermentation in these sediments.
When genes for denitrification, DNRA or genes involved in both processes were compared, genes involved in denitrification were significantly greater in abundance than those involved in DNRA (Fig. 10a, b). A study by Trimmer et al. (2013) found rates of denitrification to be high at Skagerrak sites (S4, S6, S8 and S9) neighbouring our site (Fig. 2) and the potential for DNRA to be negligible to moderate. Based on the presence of specific taxa, protein-coding genes and geochemistry results, denitrification could potentially be a more important pathway in the suboxic zone at both sites.

Conclusions
In this study, we are able to provide a model for ammonification, nitrification, NO 3 − reduction and denitrification processes in the SK and BB sediments based on the presence of corresponding genes (Fig. 4). Proteases and hydratases appeared to make up the bulk of ammonification genes at both sites. Genes associated with aerobic ammonia oxidation (amo and hao) were present and suggest AOA/AOB contribute to aerobic ammonia oxidation at the sediment-water interface at SK and BB. However, it remains unclear which one may have a more important role in both sediments. In addition, the Protein-Coding Genes Normalized abundance Fig. 9 Heat map showing normalized protein-coding gene abundances predicted to be involved in sulphur oxidation metabolism detected in SK and BB metagenomes. BB34 refers to sample 3-4 cmbsf (centimetre below sea-floor), BB67 refers to sample 6-7 cmbsf and SK68 refers to sample 6-8 cmbsf. Values for samples were scaled from 0 (minimum value) to 1 (maximum value) using a uniform scaling method implemented in the MG-RAST pipeline. Forward and reverse abundances were averaged and a single value is reported per sample presence of nrfA, nirBD and otr and NADH nitrite-reductase genes implies that DNRA could contribute to NH 4 + production where NO 3 − is available either via a respiratory or fermentative pathway. Genes for sulphide oxidation (soxADZBX) could allow for Thiothrix and Beggiatoa to carry out chemolithoautotrophic NO 3 − reduction coupled to thiosulphate or sulphide oxidation either via denitrification or DNRA. Nitrosomonas and Nitrobacter could contribute to NO 3 − reduction using nar, nap, nir and nor genes near the surface where NO 3 − is consumed. Biogeochemical and metagenomic results suggest denitrification could play the more important role in both sediments. Overall, these results show that protein-coding genes for these cycles are potentially operative in suboxic marine sediments at these sites. Our study offers the first in-depth metagenomic characterization of nitrogen cycling and associated genes of suboxic SK and BB sediments.