Thawing Yedoma permafrost is a neglected nitrous oxide source

In contrast to the well-recognized permafrost carbon (C) feedback to climate change, the fate of permafrost nitrogen (N) after thaw is poorly understood. According to mounting evidence, part of the N liberated from permafrost may be released to the atmosphere as the strong greenhouse gas (GHG) nitrous oxide (N2O). Here, we report post-thaw N2O release from late Pleistocene permafrost deposits called Yedoma, which store a substantial part of permafrost C and N and are highly vulnerable to thaw. While freshly thawed, unvegetated Yedoma in disturbed areas emit little N2O, emissions increase within few years after stabilization, drying and revegetation with grasses to high rates (548 (133–6286) μg N m−2 day−1; median with (range)), exceeding by 1–2 orders of magnitude the typical rates from permafrost-affected soils. Using targeted metagenomics of key N cycling genes, we link the increase in in situ N2O emissions with structural changes of the microbial community responsible for N cycling. Our results highlight the importance of extra N availability from thawing Yedoma permafrost, causing a positive climate feedback from the Arctic in the form of N2O emissions.


Supplementary Tables
Supplementary Table 1  b < a Results of quantitative PCR at the Kurungakh exposure. a, Gene copy number of archaeal amoA. b, Gene copy number of bacterial amoA. c, Gene copy number of total amoA (archaeal + bacterial). d, Gene copy number of 16S rRNA. e, Gene ratio of amoA/16S rRNA. The studied surfaces are arranged according to the distance from the Yedoma cliff border, with intact Holocene cover on the left and earliest thawed revegetated Yedoma on the right. Small grey symbols indicate values for individual samples, large red symbols indicate means, and error bars indicate standard error of mean (n = 3 biologically independent samples). Lower case letters indicate significant differences between studied soils (Kruskal-Wallis test with pairwise comparisons with Dunn's test; p < 0.05). DW = dry weight; VH -Vegetated Holocene cover, BYF -Bare freshly thawed Yedoma, BYE -Bare earlier thawed Yedoma, VYM -Yedoma revegetated with mosses, and VYG -Yedoma revegetated with grasses.

Supplementary Figures
The relative abundance of genes involved in key processes of nitrogen cycle: nitrogen fixation (nifH), nitrification (bacterial and archaeal amoA), nitrate reduction (narG, napA), denitrification (nir (nirK+nirS), norB, nosZ), dissimilatory nitrate reduction to ammonium (DNRA; nrfA) & anammox (hdhA). a, Relative abundances separately for each gene. b, Relative abundances grouped according to metabolic pathway. The studied surfaces are arranged according to the distance from the Yedoma cliff border, with intact Holocene cover on the left and earliest thawed revegetated Yedoma on the right. n = 3 biologically independent samples. Relative abundance of N cycling genes at the Kurungnakh exposure from all functional gene sequences captured with the targeted metagenomics tool. a, Relative abundance of nifH gene. b, Relative abundance of narG gene. c, Relative abundance of napA gene. d, Relative abundance of norB gene. e, Relative abundance of nrfA gene. f, Relative abundance of hdhA gene. g, Ratio of norB/nosZ genes. h, Relative abundance of nir/nrfA genes. The studied surfaces are arranged according to the distance from the Yedoma cliff border, with intact Holocene cover on the top of the Yedoma exposure on the left and earliest thawed revegetated Yedoma on the right side. Small grey symbols indicate values for individual samples, large red symbols indicate means, and error bars indicate standard error of mean (n = 3 biologically independent samples). Lower case letters indicate significant differences between studied soils (Kruskal-Wallis test with pairwise comparisons with Dunn's test; unadjusted p < 0.05). The relative abundance of amoA (including both archaeal and bacterial), nir (including nirK and nirS) and nosZ genes, as well as ratio of (nirK+nirS)/nosZ genes are shown in Fig. 4 of the main manuscript. VH -Vegetated Holocene cover, BYF -Bare freshly thawed Yedoma, BYE -Bare earlier thawed Yedoma, VYM -Yedoma revegetated with mosses, and VYG -Yedoma revegetated with grasses.
Effect of drying (25% reduction in water content), carbon (C) addition (as glucose) and C + nitrogen (N) addition (glucose + nitrate (NO3 ─ )) on nitrous oxide (N2O) production under different headspace conditions. Acetylene was used to inhibit the last step of denitrification (N2O reduction to N2), allowing estimation of the total denitrification rate. Box plots show lower and upper quartiles, median (thick black line), smallest and largest values without outliers (thin black line) and outliers (circles); n = 5 biologically independent samples. Stars indicate significant differences from control of the same headspace treatment (two-sided Wilcoxon signed-rank test; * unadjusted p < 0.05, ** unadjusted p < 0.01). Note the logarithmic scale on y-axes.

Supplementary Notes
Supplementary Note 1. Extended methods.

Study sites
The Kurungnakh exposure (N 72°20', E 126°17') is located on the Eastern shore of disturbed Yedoma revegetated by lower and taller canopies, respectively, of grasses Arctagrostis arundinacea and Puccinella neglecta (Supplementary Fig. S2). For the data analysis, the two bare earlier thawed surfaces 3) and 4) as well as the two revegetated surfaces with grasses 6) and 7) were merged in order to match with the Kurungnakh data.

In situ N2O fluxes
In situ nitrous oxide (N2O) fluxes were measured by the static chamber technique 8 , twice in July 2016 on Kurungnakh and once in July 2017 in Duvanny Yar. Although N2O emissions are known to commonly occur as high pulses or hot moments often outside the peak summer period 9, 10 , our own studies on bare peat surfaces have shown that temporal pattern of the emissions in these hot spots follows closely the temperature variations, peaking at high summer 11 . This (besides obvious logistical reasons) guided our selection of timing for the field campaign since we wanted to catch the peak emissions and thus maximize the emission range to reveal the drivers of spatial variability. The two measurement rounds in Kurungnakh confirm that the spatial patterns observed in N2O emissions were robust with respect to timing: the correlation between the two measurement occasions was very high (Pearson correlation with log-transformed data, R = 0.92, p < 0.001). Another challenge related to soil N2O measurements is the large spatial variability 10

Soil pore gas sampling and N2O flux estimation with a regression function
Soil pore gas at 10 cm depth was sampled at the revegetated Yedoma surfaces with a stainless-steel probe equipped with a three-way stopcock to allow tight connection with the sampling syringe. Water saturated conditions or high bulk density of the soil prevented soil pore gas sampling at most of the bare Yedoma surfaces as well as bare sand. The gas samples were stored and analyzed as described above. We used the dependence between soil pore N2O concentration at the depth of 10 cm and N2O flux observed at revegetated Yedoma surfaces of the two study sites (Supplementary Figure S5). Both parameters were logtransformed before fitting a linear regression to the data. The regression was able to explain 80% of the variability in the N2O fluxes from revegetated sites.

Soil sampling and analysis
After the chamber measurements, soil samples were taken from the topsoil (0-10 cm) of each chamber plot using a steel corer (i.d. 8.5 cm). Roots and other plant parts and stones were removed, after which the soils were homogenized by sieving (mineral soils; 5 mm mesh size) or by hand-mixing (organic soils).
Bulk density was determined from volumetric soil samples after drying until constant weight at 60 or 105 °C for organic and mineral soils, respectively. The dried soil was used in determination of particle density and content of soil organic matter (SOM), carbon (C) and nitrogen (N). Particle density was determined from dry, finely ground soil by a pycnometric method. For elemental analysis, dried soil was homogenized in a ball mill (Retsch MM301, Haan, Germany). Total content of C, organic C and N, as well as δ 13 C of SOC and δ 15 N in the bulk soil were analyzed with an elemental analyzer (Thermo Finnigan Flash EA 1112 Series, San Jose, CA, USA). For organic C analysis, inorganic C was removed from a subsample with the acid fumigation method 14 . Water-filled pore space (WFPS) was calculated from VWC measured in situ, using bulk density and particle density determined as described above. Soil pH was measured from slurries with a soil:H2O ratio of 1:4 ratio (w/v).

Gross and net N transformation rates
Nitrogen transformation rates were determined in the field laboratory from freshly sampled soils. For determination of gross N mineralization and nitrification rates, we used the pool dilution method, which is based on labelling the product pool (NH4 + for mineralization, NO3for nitrification) with the heavy N isotope, 15 N, and following the dilution in the isotopic label with time 17,18 . However, due to low mineral N content and high N immobilization of the studied soils, we were not able to detect the dilution of the 15 N stable isotope as a result of gross nitrification in any of the studied soils, and also determination of gross mineralization failed in some soil types or replicates. However, we could use the same experiment to calculate the net N mineralization and nitrification rates as described below. The net rates represent the net balance of consumption and production of inorganic N species and, thus, describe well the N saturation of the system. The initial N addition amounted to 2.1─2.6 mg N (kg DW) ─1 depending on the soil moisture content. Since the N addition was targeted to the product, not the precursor pool of each net N transformation process, we have good reasons to believe it did not significantly affect the net transformation rates. However, it may have Gross N mineralization rates per mass of dry soil were calculated based on a pool dilution model [17][18][19] . Net N mineralization rates were calculated by dividing the difference of the total mineral N content (NH4 + and NO3 -) between the first and second sampling points with the incubation time. Net ammonification and nitrification rates were calculated similarly from the change in NH4 + and NO3contents, respectively.

Experiment 1: Nitrous oxide production under different headspace conditions
Rates of N2O production and total denitrification were determined from Kurungnakh soil samples, which had been kept frozen during storage and shipment to the University of Eastern Finland. After thawing the soils, the soils were pre-incubated at 4 o C before homogenization by hand mixing and removal of visible roots. For the incubation experiments, 10 g fresh weight of organic soil (Holocene cover) and 25 g fresh weight (FW) of mineral soils (bare and revegetated Yedoma and sand) were weighed into 120 ml glass flasks. The flasks were closed with parafilm with small holes to allow gas exchange but prevent loss of moisture from soil and kept at 10 o C for three days to allow the soils to acclimatize to the incubation temperature.
The soils were incubated at field moisture content at 10 o C for six days in five replicates under three different headspace treatments: 1) oxic, 2) anoxic and 3) anoxic with acetylene. in the gas samples were determined with GC as described above. Due to intense accumulation of N2O for some soil types and treatments, a linearity test was performed to confirm the linearity of the GC-ECD up to 2000 ppm N2O (data not shown). We calculated the rate of N2O production per mass of dry soil for each sampling interval.
For oxic treatment, flux was calculated from the slope of a linear regression fitted to the first four sampling points. Within this time range, the N2O production rate was constant, after which it evened out in many of the samples. For the anoxic treatments, we report the maximum N2O production between two sampling points, because we often observed N2O consumption (treatment 2 without acetylene) or steady state (treatment 3 with acetylene) after initial N2O production, indicating reduction of N2O to N2.

Experiment 2: Response of nitrous oxide production to different moisture conditions and carbon and nitrogen sources
The aim of the second incubation experiment was to investigate the factors limiting N2O production in freshly thawed Yedoma. We dried the freshly thawed Yedoma to reduce the water content by 25%, weighed 17 g FW to incubation flasks, and incubated under three different headspace treatments as described above. For each headspace treatment, we applied three different amendments within a volume of 250 μl per flask: control (milli-Q H2O addition), addition of C (glucose; 67 μg C (g DW) -1 , equal to 0.3% of SOC), or addition of C (as above) and NO3 ─ (4.7 μg N (g DW) -1 , equal to 0.3% of TN). The GC analysis and calculation followed the procedure described above for Experiment 1.

Nucleic acid extraction and purification
For molecular studies, we used three biological replicates from five surface types on Kurungnakh, representing different stages of thermokarst and post-thaw succession: vegetated Holocene cover, freshly or earlier thawed Yedoma, revegetated Yedoma with mosses or with grasses. We extracted DNA from these samples in three technical replicates as described previously 21 . In brief, DNA samples were extracted with phenol/chloroform/IAA extraction followed by NaCl/PEG3000 precipitation and purification with PVP/Sephadex/Sepharose columns. The lab replicates were then pooled by surface prior to purification and stored in deep freeze at -80 °C until the analysis.

Quantification of amoA and 16S rRNA genes
Quantitative PCR (qPCR) of archaeal and bacterial amoA and 16S rRNA genes was performed using the reaction and cycling conditions as described previously and summarized in Supplementary Table S6 [21][22][23][24] . All reactions were performed in duplicates. Quantification of archaeal and bacterial amoA genes was based on 10-fold dilutions (10 1 ─10 8 ) of amoA gene fragments from N. viennensis EN76 and N. multiformis, respectively, amplified from cloned gene fragment with vector specific primers. For archaeal amoA genes, the qPCR efficiency was 75.5%, with a detection limit between 6.55*10 2 and 6.55*10 7 genes per reaction (25 μl). For bacterial amoA genes, the efficiency was 91.1%, with a detection limit between 1.88*10 3 and 1.88*10 8 genes per reaction. For 16S rRNA the efficiency was 83.1% with a detection limit between 2.4*10 4 and 2.5*10 9 genes per reaction. The specificity of qPCR amplification products was verified by melting-curve analysis and gel electrophoresis.

Captured metagenomics of N cycling microbes
To detect the changes in N-cycling-relevant microbial community structure with permafrost thaw and post-thaw succession, we studied the relative abundances of the key N cycling genes using a captured metagenomic tool. The method has been validated and tested for overall performance and specificity of the probes used for sequence capture 25 , and it has been successfully applied for studying N cycling genes in bioreactors treating aquaculture effluents 26 . This method is designed for targeting and sequencing the organisms carrying the key N cycling genes involved in the following processes: N2 fixation (nifH), nitrification (amoA), NO3 ─ reduction (narG, napA), denitrification (nirK, nirS, norB, nosZ), dissimilatory nitrate reduction to ammonium (DNRA) (nrfA) and anammox (hdhA) using gene specific probes following the NimbleGen SeqCap EZ protocol by Roche NimbleGen, Inc.
To develop sequence capturing tool, all known sequences for the targeted nitrogen cycling genes were collected from publicly available (NCBI, WGS, Fungene) and private databases with gene specific probes into a local server (CSC, Espoo, Finland). These functional genes were then searched for selected areas of genes with nhmmer 27 and tblastx 27,28 with manually curated and aligned databases of known and isolated clades of each gene. These target databases (TDBs) were then submitted into SeqCap pipeline with default parameters for designing up to six unique 50 mer probes for each cluster (with 90% cut-off) of each gene.
Information of these probes were then used to design NimbleGen SeqCap probe set (approximately 2M probes) for Captured metagenomics. Performance and specificity of these probes for sequence capture was validated and tested 25

Statistical analysis
All statistics were conducted with R, version 3.6.1 30

Rate and volume of thermal erosion and related N mobilization at the Kurungnakh exposure
We calculated the rate and extent of thermal erosion for the south-eastern part of the riverbank on Kurungnakh-Sise Island (length of the section 1.7 km) in QGIS v 3.14 (with SAGA and GRASS) via the following steps: We further approximated how much N was annually liberated at the Kurungnakh exposure as a result of permafrost thaw (mN; kg m ─1 year ─1 ). First, we calculated the annually eroded volume per unit area of the cliff based on the total eroded volume (Vdegraded = 223 502 m 3 ), the width of the cliff (width = 1720 m) and the mean rate of cliff retreat (rate = 3.7 m year -1 ).
We used the total ground-ice content (IC = 82 vol %) 37 , and assumed a negligible N content in the ground-ice. We multiplied the remaining sediment volume (18 vol %) with the bulk density (BD; 1.22 g cm ─3 ) and total N content (%N; 0.16 %) detected in freshly thawed Yedoma in this study (Supplementary Table S2 Table S2).