Introduction

Arctic soils store large amounts of organic matter (OM) and nutrients in permafrost1,2. When permafrost is degraded, interactions between nutrient availability and OM composition will play an important role in determining the net effects on delivery of dissolved organic carbon (DOC) and inorganic nutrients from upland landscapes to major arctic rivers and ultimately the Arctic Ocean2,3. Predicted and recently observed changes in the concentration and export of OM and nutrients vary geographically across the Arctic and as a function of permafrost structure, underlying parent material, vegetation, and disturbance3,4. Boreal forests have been burning with greater frequency due to longer growing seasons, warmer temperatures, and changing weather patterns5,6,7, thereby adding additional uncertainty to projections of how these ecosystems will respond to climate change. The response of arctic watersheds varies dramatically among regions, with both increases and decreases in DOC concentrations and inorganic nutrients following wildfire, and differing return intervals to baseline conditions8,9,10,11,12,13. Although various studies have documented the effects of fire on stream chemistry, few have evaluated how these changes will impact the processing and export of nutrients from Arctic fluvial systems. Arctic rivers mobilize large quantities of organic matter and nutrients to the Arctic Ocean14 which may be altered due to increasing permafrost thaw and wildfire frequency. Watersheds in the Central Siberian Plateau, a major component (2.5 million km2) of the Yenisei River Basin, provide a platform for understanding the long-term resilience and recovery of stream chemistry following fire in taiga landscapes. Streams of the Central Siberian Plateau drain watersheds with extensive larch forest (Larix gmelinii) growing on continuous permafrost (300–500 m deep) underlain by unglaciated basalt deposits15,16.

Uptake and transformation of nutrients and OM by stream channel processes can result in losses of C and N through mineralization and denitrification, as well as alteration of the timing of delivery of DOC and nitrate (NO3-) to downstream ecosystems17,18,19,20. Uptake rates (usually expressed as uptake velocity, Vf) are commonly used to quantify nutrient demand and the efficiency with which in-stream biota remove and process nutrients at a given concentration21,22. Uptake velocity is typically determined via in-situ nutrient manipulations since uptake is difficult to infer from downstream changes in the ambient pool of nutrients. Generally, as nutrient concentrations increase, uptake efficiency of the microbial community declines, resulting in a larger fraction of nutrients exported downstream17,23. Uptake dynamics have not been evaluated extensively in arctic fluvial systems affected by wildfires10.

Here we examine how years since the last fire affect nutrient concentrations and the quantity and composition of dissolved organic matter (DOM) entering fluvial networks using a space-for-time approach spanning a 100-year well-constrained wildfire chronosequence across 17 watersheds (Fig. 1) (Supplementary Table 1) in the Central Siberian Plateau16. We also determined uptake rates of inorganic N in a subset of these watersheds with individual nutrient pulse additions24 of nitrate (NO3−) and ammonium (NH4+) to determine how fire influences the processing and export of nutrients in arctic fluvial systems. We assess changes to the composition of DOM via optical properties and Fourier transformed ion cyclotron resonance mass spectrometry (FT-ICR MS) to provide an ultra-high resolution assessment of molecular formulae for thousands of individual DOM molecules (>10,000)25 in multiple water samples from each stream. Samples were taken in multiple years during June (freshet) and July (summer low flow).

Figure 1
figure 1

Map of sub-watersheds sampled in the Central Siberian Plateau, blue region shown in insert of Russia. Each watershed is colored with its respective fire history between 3 and >100 years since the last fire, draining into the Kochechum River and Nizhnyaya Tunguska River which are part of the greater Yenisei River Basin. On the fire history key, 10 represents watersheds burned between 3 and 7 years ago, 25 are watersheds that burned between 18 and 25 years ago, 50 represents watersheds between 51 and 57 years ago, 60 represents watersheds burned between 66 and 71years ago, and >100 are watersheds that burned up to 122 years ago. See Supplementary Table 1 for further details on individual watersheds.

Full size image

Fires altered concentrations of DOC, dissolved organic nitrogen (DON), and nutrients in streams of the Central Siberian Plateau. Concentrations of DOC (Fig. 2A; p = 5.87 × 10−9, F = 14.33, df=4) and DON (Fig. 2B; p = 2.78 × 10−5, F = 7.61, df=4) both increased with longer time since the last fire, with patterns consistent in June and July. Similar responses between June and July reflect consistency across changing flow regimes shifting from run-off dominated spring freshet in early June, and flows in July that are more representative of low-flow conditions typical of summer. In recently burned sites DOM inputs to the stream were likely reduced due to the combustion of the surface litter and organic layer, which reduces soil standing stocks of C and production of DOM from soil OM7. Post-fire declines in stream water DOC and DON concentrations can also be attributed to deepening of the active layer, which exposes mineral soil2,3 leading to increased adsorption of DOM3,15,26. Changes to concentrations of DOC and DON across the chronosequence indicate an ecosystem recovery time of approximately 50 years post wildfire disturbance, during which time the vegetation (moss, lichens, trees) undergoes regrowth re-establishing OM inputs into the soil. DOC concentrations in central Siberia exhibit a positive relationship with forest productivity and can serve as a proxy for forest C stocks15,16. The stabilization of DOC and DON concentrations after 50 years since the last fire suggests the reestablishment of an equilibrium as described by resilience theory and the ecosystem succession and nutrient retention hypothesis27,28.

Figure 2
figure 2

DOM, nutrient concentrations, and stoichiometric ratios across the burn gradient. Boxplot panels represent (A) DOC, (B) DON, (C) NO3–N concentrations, and (D) molar DOC:DIN ratios, across 17 streams sampled 2011, 2013, and 2016–2017 during June and July across the burn gradient. The paired boxplots correspond to June (dark shade) and July (lighter shade). Boxes represent interquartile range with the median value as the bold line, whiskers represent 1.5 interquartile range, and points are possible outliers. Letters denote significant differences (α = 0.05) where uppercase correspond to June and lowercase to July. Note that July data in 60 years since the last fire was excluded from statistical analysis due to low n. Significant differences were tested using ANOVA for the parametric variables and Kruskal-Wallis test for nonparametric variables. Statistics for DOC June p = 1.8 × 10−8 and July p = 8.3 × 10−4; DON June p = 0.0002 and July p = 0.005; DIN June p = 0.008 and July p = 0.04; DOC:DIN June p = 0.0008 and July p = 0.06. Respective n-values across the burn gradient for June: 16, 16,19, 4, 10; July: 5, 5, 3, 1, 10.

Full size image

Stream dissolved inorganic nitrogen (DIN) responded differently to forest fires than DOM. NO3− concentrations were highest in the most recently burned sites and decreased with increasing recovery time (Fig. 2C, p«0.05); consistent responses to fire were observed in both June and July. Multiple mechanisms have been proposed to explain the elevated NO3− concentrations often observed in recently burned watersheds on permafrost, such as an increase in soil aeration from improved drainage after fire leading to declines in soil denitrification9, increases in mineralization-nitrification of organic nitrogenous compounds29 and decreased nutrient uptake due to the reduction in biomass30,31. In contrast to DOC and DON, concentrations of NO3− demonstrate a shorter recovery time, returning to their pre-disturbance concentrations approximately at 10 years. In larch forests post-fire, growth of the moss layer and OM storage is very slow32 but there is a greater demand for inorganic N from growing shrubs and larch seedlings, which require fire to propagate32,33. This relatively short recovery of NO3− in the Central Siberian Plateau is consistent with earlier studies, which have found recovery intervals of 5–6 years after fire for inorganic N in other boreal streams30. Since the recovery of DOC and DON concentrations lags behind that of NO3−, it appears that the processes which produce, and export DOM serve as the rate-limiting step in determining the resilience and stability of heterotrophic processes in these watersheds. As a result, these Central Siberian Plateau streams at 10 years remain in a state of relatively higher energy limitation compared to their pre-disturbance steady-state (i.e. watersheds that were burned over 100 years ago) as described by molar DOC:DIN ratios (Fig. 2D; p = 9.78 × 10−5). In the study streams, DIN is mostly composed of NO3− where ammonium (NH4+) is often at or below detection limit (~5 µg/L as N). During the post-fire interval when DIN is elevated and DOC is reduced, this balance between energy and nutrient demand (as measured by stoichiometry) is re-established after approximately ten years reflecting a system that returns to a state of inorganic-N limitation and high DOM availability. During this initial recovery period when DON concentrations are relatively low, the predominant in-stream nutrient source may shift from organic to inorganic forms of N. There was no statistically significant difference in ammonium and phosphate concentrations across the burn gradient (Supplementary Table 2).

The effects of fire on landscape recovery are also reflected in DOM composition. As surface flow paths change with landscape recovery, exported OM shifts in quantity, composition, and bioavailability from microbial-derived low molecular weight DOM from deep soils to less processed and more aromatic DOM derived from surface flow paths34. This is reflected in the spectral slope ratios (SR) shifting from low molecular weight to higher molecular weight values with increasing time since fire (Fig. 3A, p = 0.002). Aliphatic compounds (including N-containing aliphatics) increased with increasing recovery time (Fig. 3B; p = 0.04, F = 2.96, df=3) while polyphenolic and condensed aromatic compounds decreased after 60 years since the last fire although the change was not statistically significant (Fig. 3C p = 0.07, F = 2.49, df=3; Fig. 3D p = 0.18, F = 1.72, df=3). With the prolonged absence of fire in boreal forests, decomposition rates decline thus promoting greater storage of OM in soils35, a change in watershed C dynamics that is reflected in the increase of aliphatic compounds. Condensed aromatic compounds, including poly-condensed aromatic compounds derived from charred material that are known as black carbon25,36 were relatively low across our study watersheds, similar to others that have also found no clear relationship between black carbon in streams and fire history37. DOM composition as assessed by ultrahigh-resolution mass spectrometry reflects the recovery of the adjacent landscape, patterns that are not evident using DOM optical properties only (Supplementary Table 3) or the four component PARAFAC model that describes mostly humic-like DOM (Supplementary Table 4).

Figure 3
figure 3

DOM composition of streams across the burn gradient. Boxplots represent (A) slope ratios (SR) from June 2013 and 2016, relative abundance of (B) aliphatic compounds including N-containing aliphatics, (C) polyphenolic compounds, and (D) condensed aromatic compounds. Streams burned 60 years ago for FT-ICR MS were excluded from statistical analyses due to low n. Boxes represent interquartile range with the median value as the bold line, whiskers represent 1.5 interquartile range, and points are possible outliers. Uppercase letters denote significant differences (α = 0.05). Only data from June 2016 for FT-ICR MS compound groups and June 2013 and 2016 for SR. Significant differences were tested using Kruskal-Wallis test for nonparametric variables. P-value for SR p = 0.008, aliphatics 0.047, polyphenolics p = 0.07, and condensed aromatics p = 0.18. Respective n-values across the burn gradient SR: 11, 9, 12, 4, 7; FT ICR: 11, 8, 9, 2, 6.

Full size image

Biogeochemical regimes in stream ecosystems reflect multiple physical, chemical, and biological processes occurring at the watershed scale as well as underlying geology and climate. Wildfire disturbance interacts with each of these drivers to determine stream chemistry and microbial biogeochemical processes. To test hypotheses regarding the potential role of other landscape scale drivers we used a backward stepwise regression to relate stream chemistry and DOM composition to watershed characteristics (area, slope, elevation, aspect, years since the last fire, and season as described by either freshet period or low discharge). Across the different models, years since the last fire and season consistently were the most significant predictors of DOC, DON, and DIN concentrations, DOC:DIN ratios, and DOM composition (Supplementary Table 5). Although season was a strong predictor of solute concentrations, the patterns across the chronosequence remain similar despite the differences in concentrations. Others have also found that time since wildfire disturbance also emerges as a strong predictor variable of stream chemistry even when additional landscape variables are included38.

The response in stream chemistry due to fires influences the processing and export of DIN. Across all watersheds, DIN uptake velocities (Vf) correlated negatively with DIN concentrations (Fig. 4A, r2 = 0.26, p = 0.005). This pattern has been observed in numerous studies in a variety of biomes, suggesting that the efficiency with which the microbial community removes DIN decreases with increasing DIN concentrations17,23,39. The changes in concentrations of DOC due to fires also influence DIN processing where streams with greater DOC:DIN molar ratios supported the greatest DIN removal (Fig. 4B r2 = 0.39, p = 0.0004). Relationships between uptake and molar ratios19,20 demonstrate the important role of DOM in DIN processing and indicate that the balance and interactions between DOM and DIN can influence the resilience of heterotrophic processes. Moreover, DOM compound groups such as polyphenolics, can enhance DIN removal (Fig. 4D, r2 = 0.36, p = 0.004); increasing relative abundance of polyphenolics leads to higher DIN uptake velocities. The positive relationships between DIN uptake and polyphenolics (i.e. vascular plant derived DOM40) suggests that uptake will be greatest in watersheds that have longer recovery times since the last fires, where inputs of DOM sources from terrestrial plants can be elevated. In addition, polyphenolics which have typically been thought to be aromatic and less bioavailable may be important energy sources in potentially energy limited streams (low primary productivity) despite high ambient DOM41. Alternatively, greater availability of DIN could have enhanced DOM decomposition by relieving any DIN limitations that inhibit the breakdown of aromatic molecules like polyphenolics. Although nutrient pulse additions do not allow us to determine what specific processes are responsible for DIN removal from the water column, the data demonstrate a close relationship between DOM and DIN uptake, supporting the idea that these patterns are due to heterotrophic or assimilative processes that require an energy source such as DOM. With shorter fire return intervals, elevated DIN concentrations, and declines in both DOM quantity and DOC:DIN ratios we can expect greater export of DIN as the strength of watersheds as sinks is greatly reduced leading to larger amounts of nutrients reaching larger rivers and ultimately the Arctic Ocean. DOC:DIN ratios provide a valuable integrator of energy versus nutrient limitation in stream ecosystems, in which energy limitation (low DOC:DIN ratios) or nutrient limitation (high DOC:DIN ratios) can be important in determining the delivery of nitrogen and OM to downstream ecosystems19,20. Downstream exports of inorganic N will likely persist for several years to a decade post fire, until DOM inputs from the watershed are re-established.

Figure 4
figure 4

Drivers of DIN uptake velocity. Uptake velocity (Vf) of NH4+ (circles) and NO3− (triangles) grouped as DIN Vf from streams across the burn gradient (10 red, 25 orange, 50 blue, 60 yellow, and>100 green) related to (A) ambient DIN concentration, (B) molar DOC:DIN ratios, and (C) relative abundance of polyphenolics. Vf values here are from 2016 through 2018. The n-values for NH4+ Vf is 18 and NO3¯ Vf is 10.

Full size image

Across arctic watersheds, the response of stream chemistry to biomass burning is not consistent. Watershed recovery following wildfire across arctic watersheds appears to be a function of many variables, including composition of vegetation and rates of regrowth, nature and extent of OM consumed by fire, and the effects of both of these on active layer depth3,7,9,12,33. Although our study focuses on streams within the Central Siberian Plateau, the results presented here can be applied to a wide array of systems in northern high-latitude and boreal regions. This will better suit predictive models of how Arctic watersheds respond to fire disturbance across the pan-Arctic. The majority of previous studies include only the assessment of stream chemistry before and after a fire8,9,13, making predictions of long-term recovery and resilience difficult to evaluate and estimate. Our study addresses this challenge by taking advantage of a long chronosequence in a space-for-time substitution and combining observational and experimental data to better understand the response of watersheds to major landscape disturbances. Differences in the response to disturbance are likely driven by the role that permafrost structure (i.e. mineral soil layer vs. deep soil OM) and its continuity, hydrology, and vegetation communities (i.e. tundra vs. taiga) play in regulating elemental fluxes and the balance of energy and nutrient demand in arctic fluvial networks. A better understanding of the effects of large-scale perturbations like fires is important to develop a greater mechanistic understanding of how site conditions drive variability in responses, exports, and resiliency to increasing alteration of arctic landscapes.

Methods

Grab samples across chronosequence

We collected stream water samples in June and July from 2016 to 2018, but data from 201111 and 201310 are included in the main body of text. Samples were filtered through pre-combusted Whatman GF/F glass fiber filters and collected in acid washed high density poly-ethylene bottles (DOC, TDN, and nutrients), pre-combusted amber glass vials (optical DOM metrics), and polycarbonate bottles (FT-ICR MS) and rinsed 3 times in the field. Samples were kept frozen except for DOM optical samples which were kept at 4 °C; all samples were transported and analyzed at the Water Quality Analysis Laboratory at the University of New Hampshire, except samples collected for FT-ICR MS which were analyzed at the National High Magnetic Field Laboratory at Florida State University.

Nutrient pulse additions

Short term nutrient pulse additions were performed at several sites across the fire chronosequence in June 2016, June and July 2017, and July 2018. Amendments consisted of NO3− as NaNO3, and NH4+ as (NH4)2SO4, with and without PO43+ as KH2PO4 but we report NH4+ uptake from both as PO43+ had no influence on NH4+ uptake metrics. Amendments were added along with NaCl as a conservative tracer mixed with stream water where we increased background nutrient and Cl concentrations 2–3X above background. Experimental reach lengths ranged between 15 m to 200 m (see Supplementary Table 6 for more details on individual additions) and were selected to exclude large pools, tributaries, or other overland flows. The amendments were released at the top of the experimental reach and samples were collected at a fixed point at the bottom of the experimental reach throughout the break through curve (BTC), where between 20 and 40 samples were collected per experiment, depending on flow conditions and changes in specific conductance monitored with a HANNA HI 991300 (HANNA Instruments, Woonsocket, RI, US). Background samples were collected as described above in duplicates before every experiment at the top and bottom of the experimental reach.

Calculations of uptake metrics follow the break through curve integration method24 where uptake length (Sw) was determined as the negative inverse of the ratio between the natural log of the fraction of background corrected nutrient and Cl retrieved and the reach length distance. Uptake velocity (Vf) was determined as:

$${V}_{f}=\frac{\frac{Q}{w}}{{S}_{w}}$$
(1)

where Q is discharge, w stream width, and Sw the uptake length. Uptake velocity is the best metric to compare across sites and dates as it is normalized for stream depth and velocity22,23,24. Uptake velocity is reported in terms of DIN (combining NO3− and NH4+ Vf) due to low statistical power if reported individually. Discharge was determined using the dilution gaging method42 using NaCl dissolved in stream water and added to the stream as an instantaneous slug addition and conductivity was logged every 5 seconds using a HOBO conductivity data logger (Onset, Bourne, MA).

Water chemistry

All samples were analyzed for DOC and total dissolved nitrogen (TDN) using high temperature catalytic oxidation with a Shimadzu TOC-V CPH/TNM. Nitrate (NO3−) was analyzed with ion chromatography using an Anion/Cations Dionex ICS-1000 with AS40 autosampler. Ammonium (NH4+) was measured using a SmartChem 200 discrete automated colorimetric analyzer using the alkaline phenate standard method. Phosphate (PO43−) was analyzed as soluble reactive phosphorus by discrete colorimetric analysis on a Seal AQ2 (Seal Analytical) by the ascorbic molybdenum blue method (EPA 365.1). DON was determined arithmetically by subtracting DIN (NO3− + NH4+) concentrations from TDN. Samples below detection limit were substituted with half the detection limit; method detection limit for NO3− (4 µg/L as N), NH4+ (5 µg/L as N), PO43− (1 µg/L as P), DOC (0.05 mg/L), TDN (0.035 mg/L). Stream chemistry data was also obtained from the same streams as this study during June 201111 and 201310. Those samples were also collected and analyzed as described here.

Optical properties

All grab samples across chronosequence from July 2013, June 2016, and July 2018 and background nutrient pulse samples from June 2016 and July 2018 were analyzed for DOM optical properties. UV absorbance was measured with a Shimadzu photo diode array detector with HPLC at 1 nm intervals between 200 to 700 nm. Specific UV absorbance at 254 nm (SUVA254) was determined by dividing the UV absorbance at 254 nm by DOC concentration, where SUVA is an index of DOM aromaticity43. Fluorescence index (FI) was determined as the ratio of 470 and 520 nm fluorescence intensities at 370 nm excitation and FI used to identify sources of DOM (terrestrial vs microbial)44. Humification index (HIX) was determined as the ratio between area under the emission spectra 435–480 nm and the sum of the peak area 330–345 nm and 435–490 nm and describes the degree of humification of DOM45. Spectral slope (S) was determined using a non-linear fit of an exponential function to the log transformed absorption spectrum in the ranges of 275–295 nm and 350–400 nm46. Slope ratio (SR) was determined as the ratio of the slope (275–295 nm) and slope (350–400 nm) which have been negatively correlated to DOM molecular weight and aromaticity46.

Excitation emission matrices (EEM) were conducted on room temperature water at excitations ranging from 240 to 450 nm at 1 nm intervals and emission ranging from 350 to 550 nm at 2.5 nm intervals using a Jobin Yvon Horiba Fluoromax-3 fluorometer. For 2016 samples EEM emission range was extended from 300 to 600 nm, but for purposes of Parallel Factor Analysis (PARAFAC) modeling these were clipped to match samples from previous years.

Parallel Factor Analysis (PARAFAC) was performed on 132 samples of stream water and soil pore water collected in 2011, 2013, and 2016. PARAFAC is a multivariate technique that decomposes fluorescence spectra into individual fluorescence components47. Modeling of PARAFAC components followed established procedures and was normalized to reduce collinearity with reversal of the normalization prior to exporting the model48. Prior to PARAFAC modeling, data in the region of first-order and second-order scatter was removed.

FT-ICR MS analysis

Samples analyzed for FT-ICR MS were collected in June 2016 and July 2018 only. DOM was obtained for FT-ICR MS by a solid phase extraction (SPE) method49. Briefly, each sample was acidified to pH 2 and passed through a Bond Elute PPL (Agilent Technologies) cartridge. Residual salts were rinsed from the cartridge with acidified (pH 2) Milli-Q water. The stationary phase was then dried with a stream of N2 and the DOM was eluted with 100% MeOH (JT Baker) to a final concentration of 50 µg/mL C and stored at 4 °C prior to analysis. Negatively charged ions were produced via direct infusion electrospray ionization at a flow rate of 700 nL/min and analyzed with a custom-built FT-ICR MS equipped with a 21-tesla superconducting magnet50.

Each signal>6σ RMS baseline noise was assigned a molecular formula with EnviroOrg©™ software51 developed at the NHMFL. The mass measurement accuracy did not exceed 200 ppb. Each molecular formula was classified based on stoichiometry; condensed aromatic (modified aromaticity index (AImod) ≥ 0.67), polyphenolic (0.67> AImod > 0.5), unsaturated, low oxygen (AImod < 0.5, H/C < 1.5, O/C < 0.5), unsaturated, high oxygen (AImod < 0.5, H/C < 1.5, O/C ≥ 0.5), aliphatic (H/C ≥ 1.5, N = 0), peptide-like (H/C ≥ 1.5, N > 0)52,53. Although FT-ICR MS allows for the precise assignment of molecular formulae to peaks that may represent multiple isomers, they describe the underlying molecular compounds comprising DOM, therefore the term compound may be used when describing the peaks detected by FT-ICR MS. Signal magnitude was converted to relative abundances by dividing the signal magnitude of an individual peak by the sum of all assigned signals. Subsequently, the percent contribution of each compound category was calculated as the percent that the relative abundance in each compound category contributed to the summed abundance of all assigned formulae.

Statistical analyses

Normality was tested for all parameters using the Shapiro-Wilk normality test (α ≤ 0.05). Differences across parameters that follow the assumptions of normality were determined using a one-way ANOVA and nonparametric Kruskal-Wallis rank sum test for parameters that violate the assumptions of normality. We used simple linear regressions to determine relationships between DIN Vf and DIN concentrations, stoichiometric ratios, and DOM compound groups. DIN concentrations and DOC:DIN ratios were log transformed to normalize data. To determine if watershed characteristics were potential predictors of patterns in stream chemistry across the fire chronosequence, we used stepwise multiple regression analysis with backward selection. Watershed characteristics (area, elevation, slope, aspect), season (as freshet period or low discharge), and years since the last fire were loaded as independent variables and DOC, DON, DIN, DOC:DIN, relative abundance of combined aliphatics and N-containing aliphatics, condensed aromatics, and polyphenolics, and slope ratio as the dependent variables. Final models were used as those with lowest AIC scores. For all statistical analyses α=0.05 and were conducted in RStudio (version 1.2.1335, RStudio, Inc. Team, Boston, MA 2016) using the stat, rcompanion, and MASS packages. The PARAFAC model was validated using split-half analysis using the drEEM toolbox48 in MATLAB (MATLAB R2015b; The Mathworks, Natick, USA, 2008).