In-depth characterization of denitrifier communities across different soil ecosystems in the tundra

In contrast to earlier assumptions, there is now mounting evidence for the role of tundra soils as important sources of the greenhouse gas nitrous oxide (N2O). However, the microorganisms involved in the cycling of N2O in this system remain largely uncharacterized. Since tundra soils are variable sources and sinks of N2O, we aimed at investigating differences in community structure across different soil ecosystems in the tundra. We analysed 1.4 Tb of metagenomic data from soils in northern Finland covering a range of ecosystems from dry upland soils to water-logged fens and obtained 796 manually binned and curated metagenome-assembled genomes (MAGs). We then searched for MAGs harbouring genes involved in denitrification, an important process driving N2O emissions. Communities of potential denitrifiers were dominated by microorganisms with truncated denitrification pathways (i.e., lacking one or more denitrification genes) and differed across soil ecosystems. Upland soils showed a strong N2O sink potential and were dominated by members of the Alphaproteobacteria such as Bradyrhizobium and Reyranella. Fens, which had in general net-zero N2O fluxes, had a high abundance of poorly characterized taxa affiliated with the Chloroflexota lineage Ellin6529 and the Acidobacteriota subdivision Gp23. By coupling an in-depth characterization of microbial communities with in situ measurements of N2O fluxes, our results suggest that the observed spatial patterns of N2O fluxes in the tundra are related to differences in the composition of denitrifier communities.


Background
Nitrous oxide (N 2 O) is a greenhouse gas (GHG) that has approximately 300 times the global warming potential of carbon dioxide on a 100-year scale [1]. Atmospheric N 2 O concentrations have increased by nearly 20% since preindustrial times, with soils-both natural and anthropogenic-accounting for up to 70% of the global emissions [2]. Despite being nitrogen (N) limited and enduring low temperatures throughout most of the year, tundra soils are increasingly recognized as important sources of N 2 O [3][4][5][6][7]. The relative contribution of tundra soils to global GHG emissions is predicted to increase in the future [8,9], as the warming rate at high latitude environments is more than twice as high than in other regions [10].
Microbial denitrification is an important source of N 2 O [11]. Denitrification is a series of enzymatic steps in which nitrate (NO 3 -wide range of archaea, bacteria, and some fungi, most of which are facultative anaerobes that switch to N oxides as electron acceptor when oxygen becomes limiting [12]. Denitrification is a modular community process performed in synergy by different microbial taxa that execute only a subset of the complete denitrification pathway [12,13]. With the growing number of microbial genomes sequenced in recent years, it has become evident that only a fraction of the microorganisms involved in the denitrification pathway encode the enzymatic machinery needed for complete denitrification [14,15]. Compared to high N 2 O-emitting systems such as agricultural and tropical soils, our knowledge of denitrifier communities in tundra soils is limited. As denitrification leads to the loss of N to the atmosphere, it enhances the N-limited status of tundra systems thus impacting both microbial and plant communities [16,17]. Investigations of denitrifier diversity in the tundra have been largely limited to gene-centric surveys using microarrays, amplicon sequencing, qPCR, and read-based metagenomics, which provide limited information on the taxonomic identity and genomic composition of community members. These studies have shown that denitrifier communities in the tundra are dominated by members of the phyla Proteobacteria, Actinobacteria, and Bacteroidetes, and that the potential for complete denitrification is usually present at the community level [18][19][20][21][22]. However, it is not known whether the complete denitrification potential occurs within discrete microbial populations or is widespread throughout populations of truncated denitrifiers lacking one or more denitrification genes. In addition, tundra soils encompass many different ecosystems, some of which are notorious N 2 O sources (e.g. bare peat surfaces [3]). N 2 O consumption is usually favoured in wetlands, where low NO 3 availability due to anoxia promotes the reduction of N 2 O to N 2 [23]. In upland soils, N 2 O fluxes vary in both time and space. Strong N 2 O sinks have been observed specially in sparsely vegetated upland soils [7], but the microbial processes underlying N 2 O consumption in these systems are largely unknown [24]. Altogether, these large differences in N 2 O fluxes across tundra ecosystems indicate differences in the structure of microbial communities, but a comprehensive understanding of the microorganisms driving N 2 O fluxes in tundra soils is lacking.
Modelling N 2 O emissions based on microbial community structure is challenging. N 2 O fluxes are characterized by a high temporal and spatial heterogeneity driven by several environmental constraints related to soil pH, N, moisture, and oxygen content [11]. In addition, our knowledge of the regulation of the denitrification process is largely based on the activity of model organisms such as the complete denitrifier Paracoccus denitrificans [25].
It has been suggested that incomplete denitrifiers that contain Nir and Nor but lack Nos contribute substantially to soil N 2 O emissions [26], while non-denitrifying N 2 O reducers, i.e., microorganisms that contain Nos but lack Nir, can represent an important N 2 O sink [27][28][29]. Furthermore, the partitioning of metabolic pathways across different populations with truncated pathwaysalso known as metabolic handoffs [30]-has been linked to higher efficiencies in substrate consumption compared to complete pathways [15,31]. However, it remains largely unclear how populations of truncated denitrifiers with different sets of denitrification genes interact with each other and the environment impacting in situ N 2 O emissions.
The paucity of in-depth knowledge on denitrifying communities in the tundra impairs our ability to model current and future N 2 O fluxes from this biome. A better understanding of the ecological, metabolic, and functional traits of denitrifiers is thus critical for improving current models and mitigating N 2 O emissions [32]. This invariably relies on the characterization of the so-called uncultured majority, i.e., microorganisms that have not been cultured to date but which comprise a high proportion of the microbial diversity in complex ecosystems [33,34]. Genome-resolved metagenomics is a powerful tool to access the genomes of uncultured microorganisms and has provided important insights into carbon cycling processes in tundra soils [35][36][37]. However, this approach has not yet been applied to investigate the mechanisms driving N 2 O fluxes in the tundra. Here, we used genomeresolved metagenomics to investigate the diversity and metabolic capabilities of denitrifiers across different tundra soil ecosystems characterised by a high variability in net N 2 O fluxes in an area of mountain tundra in Kilpisjärvi, northern Finland.

Study area and sampling
The Saana Nature Reserve (69.04°N, 20.79°E) is located in Kilpisjärvi, northern Finland (Fig. 1a). The area is part of the mountain tundra biome and is characterized by a mean annual temperature of − 1.9 °C and annual precipitation of 487 mm [38]. Sampling was performed across 43 sites during the peak of the growing season in the northern hemisphere. Our study sites are distributed across Mount Saana and Mount Korkea-Jehkas and the valley in between (Fig. 1b), and include barren soils (n = 2), heathlands (dominated by evergreen and deciduous shrubs) (n = 18), meadows (dominated by graminoids and forbs) (n = 7), and fens (n = 16) (Fig. 1c). Elevation across the sampling sites varies from 586.6 to 904.5 m.a.s.l. (Additional file 1: Table S1). Fen sites were sampled in July 2018 and all other sites in July 2017. Samples were obtained with a soil corer sterilized with 70% ethanol and, when possible, cores were split into organic and mineral samples using a sterilized spatula. In total, 69 samples (41 organic and 28 mineral) were obtained from the 43 sites (Fig. 1c, Additional file 1: Table S1). Samples were transferred to a whirl-pack bag and immediately frozen in dry ice. Samples were transported frozen to the laboratory at the University of Helsinki and kept at − 80 °C until analyses.

Soil physicochemical characterization and in situ measurement of GHG fluxes
Soil pH, moisture, and soil organic matter (SOM) content were measured from the 69 samples according to Finnish (SFS) and international (ISO) standards (SFS 300, ISO 10390, and SFS 3008). Carbon (C) and N content were measured using a Vario Micro Cube machine (Elementar, Langenselbold, Germany). In situ ecosystem-level N 2 O and methane (CH 4 ) fluxes were measured from the 43 sites using a static, non-steady state, non-flow-through system composed of a darkened acrylic chamber (20 cm diameter, 25 cm height) [4,39]. Measurements were conducted between 2nd July and 2nd August 2018, between 10 am and 5 pm. Simultaneous measurement of GHG fluxes and sampling for metagenomic sequencing was not possible due to limited resources and logistic constraints. At each site, five 25 mL gas samples were taken during a 50-min chamber closure and transferred to evacuated Exetainer vials (Labco, Lampeter, UK). Gas samples were analysed using an Agilent 7890B gas chromatograph (Agilent Technologies, Santa Clara, CA, USA) equipped with an autosampler (Gilson, Middleton, WI, USA) and a flame ionization detector for CH 4 and an electron capture detector for N 2 O. Gas concentrations were calculated from the gas chromatograph peak areas based on standard curves with a CH 4 concentration of 0-100 ppm and a N 2 O concentration of 0-5000 ppb.
Differences in physicochemical composition and rates of GHG fluxes across soil ecosystems and depths were assessed using one-way analysis of variance (ANOVA) followed by Tukey's HSD test with the lm and TukeyHSD functions in R v3.6.3 [40]. The relationship between soil ecosystem, depth, and physicochemical properties was also verified using a multivariate approach consisting of principal component analysis (PCA) and permutational ANOVA (PERMANOVA) with the package vegan v2.5-6 in R v3.6.3 (functions rda and adonis, respectively) [40,41]. C, N, and C:N ratio were not included in the multivariate dataset due to a high amount of missing data, and moisture and SOM were log-transformed prior to analysis. Due to the limited number of samples from barren sites, these were not included in the ANOVA and PER-MANOVA procedures.

Metagenome sequencing and processing of raw data
Total DNA and RNA were co-extracted as previously described [42]. Briefly, extraction was performed on 0.5 g of soil using a hexadecyltrimethyl ammonium bromide (CTAB), phenol-chloroform, and bead-beating protocol. DNA was purified using the AllPrep DNA Mini Kit (QIAGEN, Hilden, Germany) and quantified using the Qubit dsDNA BR Assay Kit (ThermoFisher Scientific, Waltham, MA, USA). Library preparation for Illumina metagenome sequencing was performed using the Nextera XT DNA Library Preparation Kit (Illumina, San Diego, CA, USA). Metagenomes were obtained for the 69 samples across two paired-end NextSeq (132-170 bp) and one NovaSeq (2 × 151 bp) runs. Two samples were additionally sequenced with Nanopore MinION. For this, libraries were prepared using the SQK-LSK109 Ligation Sequencing Kit with the long fragment buffer (Oxford Nanopore Technologies, Oxford, UK) and the NEBNext Companion Module for Oxford Nanopore Technologies Ligation Sequencing Kit (New England Biolabs). Each sample was sequenced for 48 h on one R9.4 flow cell.
We obtained more than 9 billion Illumina (1.4 Tb) and 7 million Nanopore (21.5 Gb) reads from the 69 soil metagenomes (mean: 19.9 Gb, minimum: 0.7 Gb, maximum: 82.9 Gb) (Additional file 1: Table S1). The quality of the raw Illumina data was verified with fastQC v0.11.9 [43] and multiQC v1.8 [44]. Cutadapt v1.16 [45] was then used to trim sequencing adapters and low-quality base calls (q < 20) and to filter out short reads (< 50 bp). Nanopore data were basecalled with GPU guppy v4.0.11 using the high-accuracy model and applying a minimum quality score of 7. The quality of the basecalled Nanopore data was assessed with pycoQC v2.5.0.21 [46] and adapters were trimmed with Porechop v0.2.4 [47].  Table S1. e In situ ecosystem-level nitrous oxide (N 2 O) and methane (CH 4 ) fluxes measured from the 43 sites using a static, non-steady state, non-flow-through system. Negative values represent net uptake and positive net emissions. For clarity, one outlier measurement from a meadow site (660 µg N 2 O m -2 day -1 ) was removed. In panels d and e, ecosystems followed by different letters are significantly different (one-way ANOVA, p < 0.05). Samples from barren soils were not included in the ANOVA procedure due to the limited number of samples (ND: not determined)

Taxonomic profiling
Taxonomic profiles of the microbial communities were obtained using a read-based approach, i.e., based on unassembled Illumina data. Due to differences in sequencing depth across the samples, the dataset was resampled to 2,000,000 reads per sample with seqtk v1.3 [48]. Reads matching the SSU rRNA gene were identified with METAXA v2.2 [49] and classified against the SILVA database release 138.1 [50] in mothur v1.44.3 [51] using the Wang's Naïve Bayesian Classifier [52] and a 80% confidence cut-off. Differences in community structure across soil ecosystems and depths were assessed using non-metric multidimensional scaling (NMDS) and PERMANOVA with the package vegan v2.5-6 in R v3.6.3 (functions metaMDS and adonis, respectively) [40,41]. The relationship between community structure, soil physicochemical properties, and elevation was also assessed using PERMANOVA and distance-based redundancy analysis (db-RDA) with forward selection with the package vegan v2.5-6 in R v3.6.3 (functions adonis and capscale/ordistep, respectively) [40,41]. The physicochemical dataset included only pH, moisture, and SOM due to a high amount of missing data for the other variables, and barren sites were not included in the PERMANOVA procedure due to the limited number of samples. Moisture and SOM were log-transformed prior to analysis. Relationships between the abundance of individual genera and N 2 O flux rates were assessed using linear regression in R v3.6.3 [40].

Metagenome assembling and binning
Metagenome assembling of the Illumina data was performed as two co-assemblies. One co-assembly comprised the upland soils (barren, heathland, and meadow; n = 47) and the other the fen samples (n = 22). For each co-assembly, reads from the respective samples were pooled and assembled with MEGAHIT v1.1.1.2 [53]. Assembling of the Nanopore data was done for each sample individually with metaFlye v2.7.1 [54], and contigs were corrected based on Illumina data from the respective sample with bowtie v2.3.5 [55], SAMtools v1.9 [56], and pilon v1.23 [57]. Quality assessment of the (co-) assemblies was obtained with metaQUAST v5.0.2 [58]. Binning of metagenome-assembled genomes (MAGs) was done separately for each Illumina and Nanopore (co-)assembly with anvi'o v6.2 [59] after discarding contigs shorter than 2500 bp. The two Illumina co-assemblies and the two individual Nanopore assemblies yielded more than 4 million contigs longer than 2500 bp, with a total assembly size of 21.1 Gb. Gene calls were predicted with prodigal v2.6.3 [60]. Single-copy genes were identified with HMMER v.3.2.1 [61] and classified with DIAMOND v0.9.14 [62] against the Genome Taxonomy Database (GTDB) release 04-RS89 [63,64]. Illumina reads were mapped to the contigs with bowtie v2.3.5 [55] and SAM files were sorted and indexed using SAMtools v1.9 [56]. The co-assemblies covered a significant fraction of the original metagenomic data, with an average read recruitment rate of 54.6% across samples (minimum: 22.9%, maximum: 75.8%). Due to their large sizes, Illumina co-assemblies were split into 100 smaller clusters based on differential coverage and tetranucleotide frequency with CONCOCT v1.0.0 [65]. Contigs were then manually sorted into bins based on the same composition and coverage metrics using the anvi-interactive interface in anvi'o v6.2 [59]. Nanopore contigs were binned directly without pre-clustering. Bins that were ≥ 50% complete according to the presence of single-copy genes were further refined using the anvi-refine interface in anvi'o v6.2 [59]. In addition to taxonomic signal (based on singlecopy genes classified against GTDB), either differential coverage or tetranucleotide frequency was used to identify and remove outlying contigs. The former was used for bins with a large variation in contig coverage across samples, and the latter for those with marked differences in GC content across contigs. Medium-and high-quality bins (≥ 50% complete and < 10% redundant according to the MIMAG standard [66]) were renamed as MAGs and kept for downstream analyses.

Gene-centric analyses
Functional profiles of the microbial communities were obtained using a gene-centric approach based on assembled data. For each (co-)assembly, gene calls were translated to amino acid sequences and searched against the KOfam hidden Markov model (HMM) database with KofamScan v1.3.0 [67]. Only matches with scores above the pre-computed family-specific thresholds were kept. Genes putatively identified as denitrification genes (nirK, nirS, norB, and nosZ) were submitted to further analyses to identify false positives consisting of distant homologues that are not involved in denitrification. Amino acid sequences were aligned with MAFFT v7.429 [68] and alignments were visualized with Unipro UGENE v38.1 [69]. Sequences were then inspected for the presence of conserved residues at positions associated with the binding of co-factors and active sites: nirK, Cu-binding and active sites [70]; nirS, c-heme and d 1 -heme binding sites [71]; norB, binding of the catalytic centres cyt b, b 3 , and Fe b [72]; nosZ: binding of the Cu Z and Cu A centres [72]. Sequences which did not contain the correct amino acid at these positions were removed. Finally, resulting amino acid sequences were aligned with MAFFT v7.429 [68] along with reference sequences from the genome of cultured denitrifiers [14] and a maximum-likelihood tree was computed with FastTree v2.1.11 [73] using the LG + GAMMA model. Annotation of denitrification genes was also performed for previously published genomes retrieved from GenBank. These included a set of 1529 MAGs obtained from soils in Stordalen Mire, northern Sweden [37], and all (n = 69) genomes of Acidobacteriota strains and candidate taxa (accessed on 9 October 2020).
The abundance of functional genes was computed based on read coverage with CoverM v0.6.1 [74]. For this, Illumina reads were mapped to the contigs with minimap v2.17 [75] and coverage was normalized to reads per kilobase million (RPKM). Differences in functional community structure were assessed using NMDS, PERMANOVA, and db-RDA as described above for the taxonomic profiles. Differences in the abundance of individual genes across soil ecosystems were assessed using ANOVA followed by Tukey's HSD test with the lm and TukeyHSD functions in R v3.6.3 [40]. Due to the limited number of samples from barren sites, these were not included in the ANOVA and PERMANOVA procedures. Relationships between the abundance of denitrification genes and N 2 O flux rates were assessed using linear regression in R v3.6.3 [40].

Phylogenomic analyses of MAGs and metabolic reconstruction
Phylogenetic placement of MAGs was done based on 122 archaeal and 120 bacterial single-copy genes with GTDB-Tk v1.3.0 [76] and the GTDB release 05-RS95 [63,64]. Acidobacteriota MAGs containing denitrification genes were submitted to further phylogenomic analyses alongside all genomes of Acidobacteriota strains and candidate taxa available on GenBank (n = 69; accessed on 9 October 2020). For this, the amino acid sequence of 23 ribosomal proteins was retrieved for each genome with anvi'o v6.2 [59] and aligned with MUSCLE v3.8.1551 [77]. A maximum likelihood tree was then computed based on the concatenated alignments with FastTree v2.1.11 using the LG + GAMMA model [73]. Escherichia coli ATCC 11,775 was used to root the tree.
For metabolic reconstruction, MAGs were annotated against the KOfam HMM database [67] with HMMER v.3.2.1 [61] using the pre-computed score thresholds of each HMM profile. The anvi-estimate-metabolism program in anvi'o v6.2 [59] was then used to predict the metabolic capabilities of the MAGs. A metabolic pathway was considered present in MAGs containing at least 75% of the genes involved in the pathway. Carbohydrateactive enzymes (CAZymes) were annotated with dbCAN v.2.0 based on the dbCAN v7 HMM database [78]. Only hits with an e-value < 1 × 10 −14 and coverage > 0.35 were considered.

MAG dereplication and read recruitment analysis
Prior to read recruitment analyses, Illumina and Nanopore MAGs were dereplicated based on a 99% average nucleotide identity (ANI) threshold with fastANI v1.3 [79] to remove redundancy (i.e., MAGs that were recovered multiple times across the different assemblies). Read recruitment analyses were then performed with CoverM v0.6.1 [74]. For this, Illumina reads were mapped to the set of non-redundant MAGs with minimap v2.17 [75] and relative abundances were calculated as a proportion of the reads mapping to each MAG.

Environmental characterization and in situ GHG fluxes
Our sampling design in Kilpisjärvi included two soil depths across four ecosystems that are characteristic of the tundra biome (barren soils, heathlands, meadows, and fens) (Fig. 1a-c). In previous studies, we have established in the area a systematic fine-scale sampling of microclimate, soil conditions, and vegetation in topographically distinct environments [42,80,81]. Local variation in topography and soil properties creates a mosaic of habitats characterized by contrasting ecological conditions. This makes the study setting ideal to investigate species-environment relationships and ecosystem functioning in the tundra [42,82,83].
Physicochemical composition varied across samples (Additional file 1: Table S1). Soil ecosystem and depth explained a significant fraction of the variation in pH, gravimetric soil moisture, and SOM across samples (PERMANOVA, R 2 = 0.91, p = 0.01) (Additional file 2: Fig. S1). Samples from the organic layer had higher moisture and SOM than samples from the mineral layer (one-way ANOVA, R 2 = 0.44-0.71, p < 0.001), while pH did not vary significantly between soil layers (oneway ANOVA, p > 0.05). Soil physicochemical properties did not differ across soil ecosystems in the mineral layer (one-way ANOVA, p > 0.05). In the organic layer, however, fens were characterized by higher pH, moisture, and N content (one-way ANOVA, R 2 = 0.10-0.66, p < 0.001) and, together with the meadows, lower C:N ratio (one-way ANOVA, R 2 = 0.55, p < 0.001) (Fig. 1d).
In situ measurements showed a high sink-source variability in net N 2 O fluxes across the ecosystems (Fig. 1e). Although the average N 2 O flux across all sites was small (net consumption of 6 µg N 2 O m -2 day -1 ), high N 2 O emission at rates of up to 660 µg N 2 O m -2 day -1 was observed at the meadow sites. Likewise, strong N 2 O consumption (up to − 435 µg N 2 O m -2 day -1 ) was observed particularly at the heathland and meadow sites. Net CH 4 emissions were observed exclusively at the fen sites.

Differences in microbial community structure across soils ecosystems
Read-based analyses of unassembled SSU rRNA gene sequences showed that microbial community composition differed across the ecosystems, with fen soils harbouring contrasting microbial communities compared to the other ecosystems (PERMANOVA, R 2 = 0.35, p < 0.001) (Additional file 2: Fig. S2a). No differences in community structure were observed between soil depths or the interaction between soil ecosystem and depth (PERMANOVA, p > 0.05). A significant relationship was observed between community structure and soil physicochemical properties (pH, gravimetric soil moisture, and SOM; PERMANOVA, R 2 = 0.28, p < 0.001), but not elevation (PERMANOVA, p > 0.05). Due to the significant overlap between soil ecosystem and physicochemical composition (Additional file 2: Fig. S1), we used db-RDA with forward selection to investigate in more detail the links between community structure and the environment. The best model explaining community structure comprised soil ecosystem and pH (db-RDA, R 2 = 0.40, p < 0.001). Addition of elevation did not improve the model (db-RDA, p > 0.05).
Among previously described (i.e., not unclassified) taxa, microbial communities in barren, heathland, and meadow soils were dominated by aerobic and facultative anaerobic heterotrophs such as Acidipila/Silvibacterium, Bryobacter, Granulicella, Acidothermus, Conexibacter, Mycobacterium, Mucilaginibacter, Bradyrhizobium, and Roseiarcus (Additional file 2: Fig. S2b). On the other hand, fen soils were dominated by methanogenic archaea from the genera Methanobacterium and Methanosaeta and anaerobic bacteria such as Thermoanaerobaculum, Desulfobacca, and Smithella, but also the putative aerobic heterotroph Candidatus Koribacter. We did not observe a significant relationship between the abundance of individual microbial genera and N 2 O flux rates (linear regression, p > 0.05).
Communities from different ecosystems also differed in their functional potential (Additional file 2: Fig. S2c). Denitrification genes (nirK, nirS, norB, and nosZ) were in general more abundant in the meadows and fens (oneway ANOVA, R 2 = 0.48-0.76, p < 0.001) (Additional file 2: Fig. S2d). Fen soils, which had the highest gravimetric soil moisture content, also had a higher abundance of genes involved in sulfate reduction (dsrA and dsrB) and methanogenesis (mcrA and mcrB) (one-way ANOVA, R 2 = 0.59-0.90, p < 0.001), indicating the prevalence of anoxic and reductive soil conditions in these wet sites. We did not observe a significant relationship between N 2 O flux rates and neither the abundance of individual denitrification genes nor the ratio between nosZ and nirK + nirS abundances (linear regression, p > 0.05).
However, the ratio between nosZ and nirK + nirS abundances was higher in the meadows (one-way ANOVA, R 2 = 0.29, p < 0.001) (Additional file 2: Fig. S2d), which indicates a higher potential for N 2 O consumption in this ecosystem.
To investigate their distribution across the different soil ecosystems, MAGs were dereplicated based on a 99% ANI threshold, yielding a set of 761 non-redundant MAGs (Fig. 2). On average, 15.8% of the reads from each sample were recruited by the set of non-redundant MAGs (minimum: 7.6%, maximum: 30.5%). In agreement with the read-based assessment, we observed differences in MAG composition across the soil ecosystems, with only 50 MAGs shared between the heathland, meadow, and fen soils (Additional file 2: Fig. S4a). Fen soils harboured the highest number of MAGs, with an average of 155 MAGs per sample (Additional file 2: Fig. S4b). Although barren and fen soils had similar taxonomic richness according to the read-based estimates, only a small number of MAGs was detected in the barren soils (average of four MAGs per sample). This is likely a result of limited sampling and sequencing of this ecosystem, which consisted of four samples and a total of 7.9 Gb of metagenomic data (Additional file 1: Table S1). The number of MAGs in heathland and meadow soils was similar (average of 47 and 63 MAGs per sample, respectively) (Additional file 2: Fig. S4b). In general, barren, heathland, and meadow soils were dominated by the same set of MAGs (Additional file 2: Fig. S4c). These included members of the Acidobacteriota (Sulfotelmatobacter and unclassified genera in the class Acidobacteriae), Actinobacteriota (Mycobacterium and unclassified genera in the family Streptosporangiaceae), and Proteobacteria (Alphaproteobacteria: Reyranella, Bradyrhizobium, and unclassified Xanthobacteraceae; Gammaproteobacteria: unclassified Steroidobacteraceae). On the other hand, fen soils were dominated by MAGs that were not assigned to formally described genera, including lineages of Acidobacteriota (family Koribacteraceae), Actinobacteriota (family Solirubrobacteraceae), Chloroflexota (class Ellin6529), Desulfobacterota (order Desulfobaccales), and Halobacterota.

Microorganisms from tundra soils have truncated denitrification pathways
To gain insights into the microorganisms involved with the cycling of N 2 O in tundra soils, we traced the curated denitrification genes to the set of recovered MAGs. Denitrification genes were found in 110 of the 796 MAGs (13.8%) (Additional file 1: Table S2). These were affiliated with the archaeal phylum Thermoproteota and many bacterial phyla such as Proteobacteria (classes Gamma-and Alphaproteobacteria), Acidobacteriota, Bacteroidota, Actinobacteriota, Chloroflexota, and Verrucomicrobiota (Fig. 3a). However, only 17 MAGs were assigned to a validly described genera (Additional file 1: Table S2). These included members of the Acidobacteriota (Solibacter, Sulfotelmatobacter, Terracidiphilus, and Gaiella), Myxococcota (Anaeromyxobacter), Planctomycetota (Singulisphaera), Proteobacteria (Alphaproteobacteria: Bauldia, Bradyrhizobium, Methylocella, and Reyranella; Gammaproteobacteria: Gallionella and Rhizobacter), and Verrucomicrobiota (Lacunisphaera and Opitutus). On average, 1.8% of the reads in each sample were recruited by all denitrifiers combined (minimum: 0.4%, maximum: 6.1%). In general, Genes involved in denitrification were found exclusively in MAGs with truncated denitrification pathways, i.e., MAGs missing one or more genes involved in the complete denitrification process (Fig. 3a). Of the 110 MAGs harbouring denitrification genes, the vast majority (n = 104) encoded only one of the Nir, Nor, and Nos enzymes and no MAG encoded all the three enzymes required for the reduction of NO 2 to N 2 . Unsurprisingly, co-occurrence of genes encoding the three enzymes was also not observed in any of the other genomic bins of lower quality that were discarded from the final MAG dataset (i.e., bins that were < 50% complete and/or > 10% redundant). To verify if microorganisms with truncated denitrification pathways are common in other tundra systems, we expanded our analysis to 1529 MAGs recovered from permafrost peatland, bog, and fen soils in Stordalen Mire, northern Sweden [37]. Among these, 225 MAGs (14.7%) contained denitrification genes (Additional file 2: Fig. S5). MAGs encompassed a similar taxonomic profile as observed in the Kilpisjärvi dataset, and MAGs with truncated denitrification pathways were also the norm in Stordalen Mire soils. Only one MAG, assigned to the Gammaproteobacteria genus Janthinobacterium, encoded all the Nir, Nor, and Nos enzymes required for the reduction of NO 2 to N 2 . Metabolic potential for denitrification in tundra soils. a Distribution of denitrification genes across 110 metagenome-assembled genomes (MAGs) recovered from tundra soils in Kilpisjärvi, northern Finland. Genes encoding the nitrite (nirK/nirS), nitric oxide (norB), and nitrous oxide (nosZ) reductases were annotated using a three-step approach including (1) identification using hidden Markov models from the KOfam database, (2) manual inspection for the presence of conserved residues at positions associated with the binding of co-factors and active sites, and (3) phylogenetic analyses along with sequences from archaeal and bacterial genomes (Additional file 2: Fig. S6). More information about the MAGs can be found in Additional file 1:

Microorganisms affiliated with the Chloroflexota lineage Ellin6529 are the main denitrifiers stricto sensu in fen soils
The reduction of NO 2 to NO, performed by microorganisms harbouring the nirK or nirS genes, is the hallmark step of denitrification and is often referred to as denitrification stricto sensu as it involves the conversion of a soluble substrate to a gaseous product thus leading to the removal of N from the system [12]. Of the 110 Kilpisjärvi MAGs harbouring genes involved in denitrification, 46 contained nirK/nirS genes and are thus potential denitrifiers stricto sensu (Fig. 3a). These belonged mainly to the bacterial phyla Chloroflexota, Actinobacteriota, and Proteobacteria (classes Alpha-and Gammaproteobacteria). Most MAGs (n = 43) contained the nirK gene, which encodes the copper-containing form of Nir (Additional file 2: Fig. S6a). The nirS gene encoding the cytochrome cd 1 -containing form of Nir was present in four Gammaproteobacteria MAGs (Additional file 2: Fig. S6b), including one MAG that contained both genes.
The composition of potential denitrifier stricto sensu communities differed across the ecosystems (Fig. 3b). MAGs belonging to the Alphaproteobacteria class of the Proteobacteria were the most abundant in the barren, heathland, and meadow soils, particularly the MAG KUL-0154 assigned to the genera Bradyrhizobium (Fig. 4). Two other Alphaproteobacteria MAGs that do not correspond to formally described genera in the families Acetobacteraceae and Beijerinckiaceae (KUL-0057 and KUL-0056, respectively) were also found at high abundances. In addition, one Actinobacteriota MAG assigned to an uncharacterized genus in the family Gaiellaceae (KWL-0073), was abundant in the meadow soils. On the other hand, fen communities were dominated by MAGs belonging to the phylum Chloroflexota (Fig. 3b), which included seven MAGs assigned to the class-level lineage Ellin6529 (Fig. 4).
None of the Ellin6529 MAGs that were dominant in the fen communities contained the key genes involved in autotrophic carbon fixation, dissimilatory sulfate reduction, dissimilatory nitrate reduction to ammonia, and nitrogen fixation (Additional file 1: Table S3). Analysis of genes encoding terminal oxidases involved in the aerobic respiratory electron chain revealed that all seven Ellin6529 MAGs harboured the coxABC genes encoding the aa3-type cytochrome c oxidase. Four MAGs also contained the cydAB genes encoding the cytochrome bd ubiquinol oxidase, a terminal oxidase with high affinity for oxygen that also plays a role in preventing the inactivation of oxygen-sensitive enzymes and protecting against oxidative and nitrosative stress, toxic compounds such as cyanide, and other stress conditions such as high temperature and high pH [84,85]. The dominant MAGs in the barren, heathland, and meadow soils encoded a different set of aerobic terminal oxidases. In addition to the cydAB genes, the MAGs KUL-0057 and KUL-0154 also contained the cyoABCD genes encoding the cytochrome o ubiquinol oxidase, which is the main terminal oxidase under highly aerobic conditions [86], and KUL-0057 also contained genes encoding the cbb3-type cytochrome c oxidase, a terminal oxidase with high affinity for oxygen [87]. Genes involved in the Calvin cycle (e.g., rbcL, rbcS, and prkB) were found in the Bradyrhizobium MAG (KUL-0154), and none of the key genes for autotrophic carbon fixation pathways were present in the other Alphaproteobacteria MAGs that were dominant in the barren, heathland, and meadow soils.

Acidobacteriota with the potential to reduce NO and N 2 O are abundant in the fens
The stepwise reduction of NO to N 2 O and N 2 carried out by microorganisms containing the norB and nosZ genes, respectively, represents the final step of denitrification and the main biotic control on N 2 O emissions. Soil denitrification rates depend on multiple environmental conditions such as adequate moisture and inorganic N availability, but whether it results in the emission of N 2 O or N 2 is ultimately linked to a balance between the activity of NO and N 2 O reducers [11,15]. norB and nosZ genes were identified in 62 and 9 Kilpisjärvi MAGs, respectively, belonging mostly to the phyla Actinobacteriota, Bacteroidota, Acidobacteriota, and Proteobacteria (class Alphaproteobacteria) (Fig. 3a). Apart from one Gemmatimonadota and one Acidobacteriota MAG, norB-and nosZ-containing MAGs were almost exclusively non-denitrifiers stricto sensu, i.e., they did not harbour the nirK/nirS genes involved in the reduction of NO 2 to NO. Most MAGs (n = 48) harboured a norB gene encoding the monomeric, quinol-dependent form of Nor (qNor), while the remaining MAGs (n = 8) encoded the cytochrome c-dependent Nor (cNor) (Additional file 2: Fig. S6c). In regards to the nosZ gene, most MAGs (n = 6) contained sequences affiliated with the clade II (also known as atypical) NosZ [14,15,27] (Additional file 2: Fig. S6d). Only four MAGs contained both the norB and nosZ genes and thus have the potential to reduce NO completely to N 2 (Fig. 3a).
As observed for the denitrifier stricto sensu communities, the communities of potential NO and N 2 O reducers also differed between the ecosystems (Fig. 3b). MAGs assigned to the Alphaproteobacteria class of the Proteobacteria were the most abundant in the barren, heathland, and meadow soils. In particular, the MAG KWL-0112 assigned to the genera Reyranella was the dominant norB-containing MAG, while KUL-0116 (belonging to an uncharacterized genus in the family Acetobacteraceae) was the dominant MAG harbouring the Pessi et al. Environmental Microbiome (2022) 17:30 nosZ gene (Fig. 4). On the other hand, fen communities were dominated by Acidobacteriota MAGs (Fig. 3b), particularly the norB-and nosZ-containing MAG KWL-0326 affiliated with the class Thermoanaerobaculia (Fig. 4). This MAG contained the same set of genes encoding aerobic terminal oxidases as found in the nirK-containing Fig. 4 Relative abundance of metagenome-assembled genomes (MAGs) harbouring denitrification genes across different soil ecosystems in the tundra. MAGs were recovered from soils in Kilpisjärvi, northern Finland, and annotated for genes encoding the nitrite (nirK/nirS), nitric oxide (norB), and nitrous oxide (nosZ) reductases using a three-step approach (see methods). Relative abundances were computed as the proportion of reads mapping to each MAG. MAG taxonomy is based on the Genome Taxonomy Database (GTDB) release 05-RS95 Ellin6529 MAGs that were dominant in the fen sites, namely coxABC and cydAB (Additional file 1: Table S3). No genes involved in carbon fixation, dissimilatory sulfate reduction, dissimilatory nitrate reduction to ammonia, and nitrogen fixation were found in any of the dominant norB-and nosZ-containing MAGs.
To elucidate the phylogenetic placement of the Acidobacteriota MAGs and to verify if the potential for NO and N 2 O reduction is present in other members of this phylum, we analysed all available genomes of Acidobacteriota strains and candidate taxa available on Gen-Bank. This revealed that genes encoding the Nir and Nos enzymes are widespread across the phylum Acidobacteriota (Fig. 5). Genes encoding the Nor enzyme were present in all but one of the six Acidobacteriota subdivisions with genomes from cultured representatives. This included the strains Acidobacterium ailaaui PMMR2 (subdivision Gp1), Acidipila sp. 4G-K13 (Gp1), Silvibacterium bohemicum DSM 103,733 and S. bohemicum S15 (Gp1), Acidobacteriaceae bacterium URHE0068 (Gp1), Edaphobacter aggregans DSM 19,364 (Gp1), Luteitalea pratensis DSM 100,886 (Gp6), Geothrix fermentans DSM 14,018 (Gp8), and Thermoanaerobaculum aquaticum MP-01 (Gp23), as well as the candidate taxa Candidatus Koribacter versatilis Ellin345 (Gp1), Candidatus Sulfotelmatomonas gaucii SbA5 (Gp1), and Candidatus Solibacter usitatus Ellin6076 (Gp3). On the other hand, genes encoding the Nos enzyme were found only in members of the subdivisions Gp6 and Gp23.

Discussion
The 796 MAGs obtained in the present study by a manual binning and curation effort represent one of the largest genomic catalogues of microorganisms from tundra soils to date. Earlier gene-centric investigations have revealed the potential for complete denitrification in tundra soils [22,88], however, these approaches fail to reveal the wider genomic context of the genes involved in this pathway. By applying the genome-resolved metagenomics approach, we traced denitrification genes to specific microbial populations, thereby allowing a detailed investigation of the genomic makeup of potential denitrifiers in tundra soils. This approach also enabled us to access the genomes of uncultured, poorly characterized taxa, which comprise the majority of the microorganisms in soils and other complex ecosystems [33,34].
Our genome-resolved survey revealed that denitrification across different tundra soil ecosystems is dominated by microorganisms with truncated denitrification pathways (i.e., harbouring only a subset of the genes required for complete denitrification), most of which represent poorly characterized taxa without cultured representatives. The congruence of these findings in both our original dataset of northern Finland soils and a re-analysis of a comprehensive metagenomic dataset from soils in northern Sweden [37] suggests that truncated denitrification pathways are not a methodological artifact arising from the metabolic reconstruction of fragmented genomes. Indeed, recent genome-resolved investigations have shown that cross-feeding between microorganisms with truncated metabolic pathways, also known as metabolic handoffs, are the norm across a wide range of ecosystems such as grassland soil, aquifer sediment, groundwater, and the ocean, and not only in relation to denitrification but other redox transformations as well [30,89,90]. Although it has been established that denitrification is a community effort performed by different microbial populations [12][13][14][15], these genome-resolved metagenomic studies are beginning to reveal a more indepth, ecosystem-centric representation of the denitrification pathway. In addition to their predominance in genomic databases [14], it appears that truncated denitrifiers are also dominant within defined ecosystems across various terrestrial and aquatic biomes, including the tundra. It has been suggested that the partitioning of metabolic pathways across different populations via metabolic handoffs is advantageous as it eliminates competition between enzymes accelerating substrate consumption [15,31] and provides flexibility and resilience to the communities in face of environmental disturbances [30]. We further hypothesize that the predominance of denitrification pathways characterized mostly by metabolic handoffs in tundra soils could be related to N limitation. If metabolic handoffs enable a more effective substrate consumption as previously suggested [15,31], truncated denitrification pathways would be favoured in tundra soils which are mostly N limited but undergo rapid surges in N availability, e.g., during the spring melting season [91].
Tundra ecosystems are typically heterogeneous. Previous studies in the Kilpisjärvi region have shown that soil properties such as pH and moisture do not have any strong relationship with the macrotopography of the area (50-500 m scale). Instead, environmental variation is controlled by the fine-scale mesotopographic variation of the relief (2-20 m scale), resulting in a mosaic of different soil ecosystems with contrasting vegetation [42,[80][81][82][83]. Our results agreed with this observation and showed that denitrifier communities in the tundra differ between drier upland ecosystems (barren, heathland, and meadow soils) and waterlogged fens. This is likely related to differences in soil moisture affecting oxygen availability in these ecosystems. The dominant denitrifier populations in the oxic dry upland soils, related to the genera Bradyrhizobium, Reyranella, and other uncharacterized genera in the class Alphaproteobacteria, encoded aerobic terminal oxidases that are active under highly aerobic conditions as well as oxidases with high oxygen affinity [86,87]. The former likely provides an adaptive advantage in these soils by allowing rapid aerobic growth under standard conditions of high oxygen availability, and the latter would sustain growth in microoxic niches within the soil matrix and during periods of reduced oxygen availability (e.g., during the spring melting season).
On the other hand, fen soils are continuously inundated because they are located at lower topographic positions where the water table is permanently at or near the soil surface. The result is a mostly anoxic environment due Metabolic potential for denitrification among members of the phylum Acidobacteriota. Phylogenomic analysis of 85 Acidobacteriota metagenome-assembled genomes (MAGs) containing denitrification genes recovered from tundra soils in Kilpisjärvi (northern Finland) and Stordalen Mire (northern Sweden), and 69 genomes of Acidobacteriota strains and candidate taxa. Maximum likelihood tree based on concatenated alignments of 23 ribosomal proteins and rooted with Escherichia coli ATCC 11775 (not shown). Genes encoding the nitrite (nirK), nitric oxide (norB), and nitrous oxide (nosZ) reductases were annotated using a three-step approach (see Methods) to the slow rate at which oxygen diffuses into the waterlogged soil, favouring reduced rather than oxidized soil chemistry. In line with this, we found a predominance of anaerobic processes in the fens, including a higher abundance of genes involved in denitrification, sulfate reduction, and methanogenesis, the latter supported by in situ measurements showing net CH 4 emission at the fen sites. Communities of potential denitrifiers in the fen soils were dominated by somewhat enigmatic taxa, namely potential NO 2 reducers affiliated with the class Ellin6529 of the Chloroflexota and NO/N 2 O reducers assigned to the subdivision Gp23 of the Acidobacteriota. Both groups are major members of microbial communities in soils worldwide [92], and RNA-based investigations have shown that they are active in tundra soils during both summer and winter seasons [42,93]. Thermoanaerobaculum aquaticum MP-01, the only cultivated member of the Acidobacteriota subdivision Gp23, is a strictly anaerobic bacterium that has been shown to use Fe and Mn, but not NO 3 nor NO 2 -, as electron acceptors in anaerobic respiration [94]. However, studies investigating the use of nitrogen oxides in anaerobic respiration usually provide soluble NO 3 -or NO 2 as electron acceptors, not the gases NO and N 2 O, which bias against truncated denitrifiers that do not contain the narG and nirK/nirS genes [95]. Ellin6529-formerly G04-were first detected by culture-independent methods in alpine tundra wet meadow soil in the Colorado Rocky Mountains, USA [96], and later isolated in a study targeting slow-growing and mini-colony forming bacteria from Australian agricultural soil [97]. However, their ecological, physiological, and metabolic preferences remain largely unknown. Their genomic composition and high abundance in the water-logged, anoxic fen soils suggest that the Ellin6529 and Gp23 populations found in this study are likely able to grow anaerobically with the use of NO and N 2 O as electron acceptors. However, it is known that in addition to their role in anaerobic respiration, NO and N 2 O reduction can be used as a detoxification mechanism or as electron sink for metabolism. For example, the aerobe Gemmatimonas aurantica T-27 is not able to grow on N 2 O alone, but can use N 2 O as electron acceptor transiently when oxygen is depleted [98].
In addition to microbial community structure, differences in N 2 O fluxes observed between upland and fen soils also appear to be linked to soil moisture. Some of the drier upland sites investigated were hotspots of N 2 O consumption. This is particularly interesting for the acidic heathland soils, as low pH is known to impair the expression of the NosZ enzyme thus promoting N 2 O emission [99,100]. On the other hand, fens had close to net-zero N 2 O fluxes, which is in line with previous observations for water-saturated soils both in the tundra [7] and worldwide [11,13]. This has been linked to lower rates of N mineralization and nitrification in anoxic ecosystems, which limit the availability of NO 3 and NO 2 and promote complete denitrification, resulting in N 2 as end product rather than N 2 O. Indeed, supplementing fen soils in the tundra with NO 3 and NO 2 has shown to promote N 2 O emissions [101]. Moreover, climate change models predict lowering of the water table in high-latitude wetlands, which could lead to increased N 2 O emissions from these ecosystems which contain substantial amounts of both C and N bound to the soil organic matter [102,103].

Conclusions
A better understanding of denitrification is paramount for our ability to model N 2 O emissions and mitigate climate change. High-latitude environments in particular have experienced amplified warming in recent decades, a trend that is likely to continue in the coming centuries. As mechanisms of GHG emissions are very climate sensitive, the contribution of tundra soils to global GHG atmospheric levels is thus predicted to increase in the future leading to a positive feedback loop. Compared with CO 2 and CH 4 , measurements of N 2 O fluxes in tundra soils are sparse and are rarely coupled with a characterization of the microorganisms involved, making the magnitude and drivers of N 2 O fluxes across the polar regions uncertain. While microorganisms with truncated denitrification pathways appear to dominate the denitrifier communities investigated here, the potential for complete denitrification was present at the ecosystem level. In addition to a better monitoring of N 2 O emissions throughout the tundra biome, our results suggest that a better understanding of the contribution of tundra soil to global N 2 O levels relies on the elucidation of the regulatory mechanisms of metabolic handoffs in communities dominated by truncated denitrifiers.
Additional file 2. Fig. S1. Physicochemical composition of tundra soils in Kilpisjärvi, northern Finland. Fig. S2. The microbial diversity of Kilpisjärvi soils as seen using a gene-centric approach. Fig. S3. Genome-resolved metagenomics of tundra soils. Fig. S4. Overview of the microbial diversity in Kilpisjärvi soils based on a genome-resolved approach. Fig. S5. Metabolic potential for denitrification in Stordalen Mire soils. Fig. S6. Phylogeny of a) nirK, b) nirS, c) norB, and d) nosZ sequences from metagenomeassembled genomes (MAGs) recovered from tundra soils in Kilpisjärvi, northern Finland.