Shortened Fire Intervals Stimulate Carbon Losses from Heterotrophic Respiration and Reduce Understorey Plant Productivity in Boreal Forests

Fire frequency is increasing with climate warming in the boreal regions of interior Alaska, with short fire return intervals (< 50 years) becoming more common. Recent studies suggest these “reburns” will reduce the insulating surface organic layer (SOL) and seedbanks, inhibiting black spruce regeneration and increasing deciduous cover. These changes are projected to amplify soil warming, increasing mineral soil organic carbon (SOC) decomposition rates, and impair re-establishment of understorey vegetation and the SOL. We examined how reburns changed soil temperature, heterotrophic soil respiration (RH), and understorey gross primary production (GPP), and related these to shifts in vegetation composition and SOL depths. Two distinct burn complexes previously covered by spruce were measured; both included areas burned 1x, 2x, and 3x over 60 years and mature (≈ 90 year old) spruce forests underlain by permafrost. A 2.7 °C increase in annual near-surface soil temperatures from 1x to 3x burns was correlated with a decrease in SOL depths and a 1.9 Mg C ha−1 increase in annual RH efflux. However, near-surface soil warming accounted for ≤ 23% of higher RH efflux; increases in deciduous overstorey vegetation and root biomass with reburning better correlated with RH than soil temperature. Reburning also warmed deeper soils and reduced the biomass and GPP of understory plants, lessening their potential to offset elevated RH and contribute to SOL development. This suggests that reburning led to losses of mineral SOC previously stored in permafrost due to warming soils and changes in vegetation composition, illustrating how burn frequency creates pathways for accelerated regional C loss.


INTRODUCTION
Boreal forests circle the northern hemisphere, spanning 11% of the Earth's terrestrial surface and storing between a third and two-thirds of all carbon (C) found in global forests, largely as soil organic carbon (SOC) (Bonan and Shugart 1989;Pan and others 2011;Bradshaw and Warkentin 2015). In North America, boreal black spruce (Picea mariana (Mill.) BSP) forests play a critical role in SOC storage through the development of a thick, lowdensity surface soil organic layer (SOL). This SOL stores as much C as the overstorey at forest maturity (Van Cleve and others 1983; Kane and Vogel 2009;Brown and Johnstone 2011;Mann and others 2012;Alexander and Mack 2016;Walker and others 2019) and insulates the underlying mineral soil in summer (Zhu and others 2019), reducing the active layer thickness and promoting the development and retention of discontinuous permafrost throughout the region (Bonan and Shugart 1989;Yoshikawa and others 2003). Perennially frozen soils inhibit decomposition (O'Donnell and others 2009;Deluca and Boisvenue 2012), facilitating the accrual of SOC over millennia (Grosse and others 2011), and comprise an estimated 88% of the total boreal SOC pool (1672 Pg C; Tarnocai and others 2009), or nearly twice the mass of the atmospheric C pool (Schuur and others 2008). Process models suggest that permafrost loss may increase SOC vulnerability to decomposition losses and shift ecosystems into net C sources, causing positive feedback to climate warming (Koven and others 2011;McGuire and others 2018).
Black spruce forests have dominated the Alaskan interior for the past 5000 years, expanding during a timeframe where cold winters and stand-replacing fires at % 145 year intervals facilitated a self-replacing successional pattern, along with a mosaic of transitional vegetation communities (Van Cleve and others 1983;Johnstone and Chapin 2006;Lloyd and others 2006;Higuera and others 2009;Johnstone and others 2010b;Brown and John-stone 2011). However, increased fire frequency over the past century has substantially shortened the mean fire return interval to % 105 years (Alaska Dept. of Natural Resources 2005-2019; Todd and Jewkes 2006;Kasischke and others 2010), and some models predict intervals of less than 50 years in some areas by the year 2100 (Balshi and others 2009a). This has the potential to change boreal forest ecosystem structure and soils by altering a long-standing pattern of interactions amongst fire, vegetation succession, and climate (Weber and Flannigan 1997;Balshi and others 2009a;Kasischke and others 2010;McGuire and others 2016;Whitman and others 2019). Multiple short-interval burns (that is, ''reburns'') leading to less than 50 year return intervals can impede black spruce regeneration and resemble more intense longer-interval burns by reducing the aerial seedbank, consuming more of the SOL, and increasing mineral soil exposure (Brown and Johnstone 2012;Whitman and others 2019), conditions favouring the recruitment of faster growing angiosperms (Kasischke and Johnstone 2005;Johnstone and Chapin 2006;Greene and others 2007;Johnstone and others 2010a;Hayes and Buma 2021). Under deciduous cover, the micro-climate and groundlevel vegetation communities critical to building the SOL in mature spruce forests may fail to reestablish (Natalia and others 2008;Alexander and Mack 2016), preventing the re-development of this insulating layer and resulting in longer-term warming of the underlying soil. Experimental and modelling studies indicate this SOL loss will increase the abundance and potential persistence of deciduous communities in upland areas (Johnstone and Chapin 2006;Beck and others 2011; Mekonnen and others 2019; Whitman and others 2019; Johnstone and others 2020; Mack and others 2021), leading to long-term soil warming and permafrost decline (Hinzman and others 2003) and consequently to changes to SOC dynamics in the region (Balshi and others 2009b;Turetsky and others 2011).
Models predict an acceleration of climate warming and reduced fire return intervals for much of the North American boreal forest (Balshi and others 2009a;others 2010a, 2011;Kasischke and others 2010). As fires become larger and more frequent, the frequency of shortinterval reburns will increase, as will the area of forest affected by them (Whitman and others 2019; Buma and others 2020). A concurrent shift to greater deciduous cover across the interior of Alaska would inevitably warm the region's soils (Johnstone and others 2010a; Brown and John-Shortened Fire Intervals Stimulate Carbon Losses from Heterotrophic Respiration stone 2011). This warming has the potential to stimulate soil respiration (RS) rates, that is, the combined respiration of heterotrophic soil organisms (RH) and plant roots, which increase with soil temperature as a first-order kinetic response, often exponentially (Lloyd and Taylor 1994;Fang and Moncrieff 2001). Following fire, root respiration initially decreases with vegetation loss but increases over time with regeneration, while RH can increase with soil warming (Zhuang and others 2002;Bergner and others 2004). As a consequence of more short-interval reburns over coming decades, RH may continue to increase with SOL loss and climate change-driven soil warming. Although most RH is generated from SOC-rich substrates near the soil surface, permafrost thaw increases the overall volume of the active layer (Viereck 1983;O'Donnell and others 2009), providing decomposers with access to a larger, previously unavailable SOC pool (Davidson and Janssens 2006;Schuur and others 2009;Brown and Johnstone 2011;Gibson and others 2019). Furthermore, a shift from spruce to deciduous cover may stimulate decomposition of mineral SOC decomposition by changes to root and litter quality and turnover rates (O'Neill and others 2002; Laganiè re and others 2010; Adamczyk and others 2019).
Increased nutrient availability and ground surface insolation following fire stimulates germination and growth of vegetation, though reestablishment of pre-fire composition may be impaired by more intense or successive short-interval fires (Ott and others 2006;Kurkowski and others 2008;Hayes and Buma 2021). An increased loss of the SOL associated with higher fire intensity has been linked to a decline in seedbank and viable rhizomes/roots from low shrubs and herbaceous vegetation and a shift in moss communities (Schimmel and Granstrom 1996), likely altering patterns of early successional establishment and composition. Despite comprising a small fraction of the total standing biomass in upland spruce forests, this ground-level vegetation has a high turnover rate and is responsible for up to half of all aboveground net primary productivity in mature boreal forests, providing substantial litter inputs to the forest floor (Wardle and others 2003;Nilsson and Wardle 2005;Turetsky and others 2010;Ikawa and others 2015). Consequently, the primary productivity of the understorey may comprise a substantial early successional transitional C sink that builds the SOL and persists to maturity in black spruce forests.
Most studies of early successional primary production in boreal forests focus on ecosystem-scale measurements from eddy-covariance flux towers (Kolari and others 2004;Goulden and others 2011;Ueyama and others 2013) or models (Amiro and others 2000;Zhuang and others 2002;Chertov and others 2009) that are unable to distinguish contributions from differing plant functional groups (for example, vascular/non-vascular) and height class (for example, ground-level vs. taller understorey shrubs vs. emerging overstorey). Direct measurements of net ecosystem exchange from the ground surface (NEE GS ) using transparent flux chambers can be used in conjunction with measurements of RS from opaque flux chambers to estimate gross primary production (GPP) from ground-level vegetation (GPP GV ), allowing for a microscale understanding of the contribution of the ground-level vegetation to C dynamics. Few studies have published GPP GV measurements in the boreal forest using flux chambers (Goulden and Crill 1997;Heijmans and others 2004;Kolari and others 2006) and none to our knowledge have been specifically designed to assess reburns. Understanding this component is critical because successional pathways are likely altered by reburning, influencing which vegetation communities are contributing to photosynthetic uptake and their subsequent inputs to soil C pools via litter and root turnover.
Our overall study objective was to estimate the response of RS, RH, and GPP fluxes to fire frequency in a chronosequence of areas burned once, twice, or thrice over the past 60 years, and to link these to changes in the soil and ground surface environment. We measured these fluxes at two study sites % 200 km apart in the Alaskan interior. Both sites were covered by mature spruce forest prior to the first burn, had a similar pattern of fire event frequency and timing, and contained areas of % 90-year-old mature spruce forest to provide reference measurements. We hypothesized that repeated burns at shorter intervals would: i) thin the SOL, leading to higher soil temperatures that increase RH rates, ii) reduce understorey biomass and GPP GV , primarily from vascular plants, and iii) increase net C loss at the ground surface (that is, NEE GS ).

Study Area
Two upland boreal forest study sites % 200 km apart were selected within the interior of Alaska from areas previously identified via satellite and aerial photography (Figure 1). The ''Dalton'' site (DA) is located south of the Dalton highway, 33 km NW of Livengood, Alaska (65.69°, -149.16°). The ''Steese'' site (ST) is located north of the Steese highway, 10 km west of Central, Alaska (65.59°, -144.90°). Soils were primarily silt loams, accrued from either eolian (Dalton) or fluvial (Steese) deposition processes, minimally weathered due to cold soil temperatures and the relatively young age (10,000-15,000 years) of the soils (Harden and others 1992). Based on historical aerial imagery from the Alaska Large Database and modern remotely sensed fire perimeters assessed by Hayes and Buma (2021) both sites were predominantly covered by black spruce in the early 1950s, and underlain by discontinuous permafrost (V. Romanovsky, 2018;personal communication). Beginning in 1953, a series of severe forest fires (complete aboveground mortality) at both sites produced a mosaic of mature black spruce forest and regener-ating forest areas burned up to three times over the previous 60 years (Table 1), with all burned most recently between 2003-2005. Four burn histories were represented in our study as of the year 2019: once-burned areas (1x: burned 14-16 years earlier), twice-burned areas (2x: burned 14-16 and 52-66 years earlier), thrice-burned areas (3x: burned 14-16, 28-45, and 52-62 years earlier), and mature black spruce forest (Mature: burned % 90-100 years ago) dated using tree ring counts (Hayes and Buma in prep; Table 1).

Plot Selection and Development
Representative flux plot locations were established within 1x, 2x, and 3x burned areas at each site, accounting for variation in slope, (0-5%), hillslope position (plateau, backslope, toeslope), and aspect (NW to NE) within each burn frequency. Most plots were established on well-drained backslopes at the hillier Dalton site, while plots at the flatter Steese site were mostly on plateaus with moderate drainage. Due to the minimal amount of accessible mature forest at each study site, mature forest reference plots (''Mature'') at the Dalton site were located on a foot slope with a SE facing aspect. Permafrost was found at a depth of 15-30 cm in the mineral soil at most mature forest areas. No permafrost was found within 45 cm of the soil surface at any burn frequency, despite the north facing aspects of the burned plots.
Each burn contained three to four plots (20 m diameter) at both sites, for a total of 26 plots (Table S1). These were subdivided into eight subplots (area = 1 2 m) radially distributed at a distance of 10 m from the centre to be used as soil surface C flux measurement locations ( Figure S1). To isolate the response of heterotrophic respiration, every other subplot was ''trenched'' to sever roots by incising a 2-4 mm gap to a depth of 30 + cm encircling each subplot (% 80 9 80 cm) using a flat shovel and serrated knife. Vascular plants were also removed from trenched plots, either by pulling to remove roots where possible or cutting to the stem base to minimise soil disturbance. Moss and lichens were left in place as these were not separable without significant disruption to the soil surface and to measure the GPP of non-vascular vegetation alone. Trenches were re-incised at least one week prior to the first set of measurements and maintained biweekly to prevent root in-growth. Similarly, vegetation regrowth was removed prior to measurement. The distance between ''trenched'' and ''intact'' plots (% 8 m) minimised the potential for proximal root damage of plants within intact plots.

Mineral Soil Root Sampling
Mineral soil samples were collected for root analysis using a 7.25 cm diameter tube corer (AMS Inc, American Falls, ID, USA) from random locations within burn frequencies (n = 3-4 per burn) in August 2018 that did not co-occur with flux measurement plots. At each location, three sub-samples % 10 m apart were extracted to a depth of 45 cm after removing the SOL, and later combined during processing to obtain representative samples. Samples were kept refrigerated until sorted to extract live roots for mass analysis. Live roots were identified and sorted using methods described by Ruess and others (1996).

SOL and Vegetation
The depth of the SOL was assessed using a spade to dig a face along the vertical soil profile and measuring the distance between the interface of mineral and organic horizons and the top of the organic horizon, excluding live moss. Values were estimated as the mean of three measurements collected from random locations at each plot. Vegetation was also measured to genus or species where possible in three categories: ground-level,  1967, 2003 35 1953, 2005 41 3x 1967, 1991, 2005 24 1957, 1974, 2004 27 ** Oldest dated spruce trees in site from Hayes and Buma (in prep); exact fire date unknown.
tall understorey vegetation, and overstorey vegetation. Ground-level vegetation was defined as all vegetation less than 30 cm in height and separated into vascular and non-vascular components. Vascular ground-level biomass was estimated using a method combining ''hits'' within a point-frame, percent ground-cover estimates, and biomass collection as defined by Schuur and others (2007), whereas the non-vascular component was estimated as percent ground-cover only. Tall understorey vegetation root-collar area was estimated by measuring the root-collar diameter of all woody stems greater than 30 cm in height and less than 1 cm diameter at breast-height (1.37 m-DBH) within 1 m of subplot centres. Overstorey vegetation basal area (BA) was estimated by measuring the DBH of all woody vegetation with greater than 1 cm DBH within 2 m from the centre of untrenched subplots (n = 12-16 per burn frequency).

Carbon Flux Measurements
Fluxes of CO 2 from the soil surface and plants less than 30 cm in height (FCO 2 ) were measured using clear acrylic flux chambers (30 9 30 9 30 cm) placed directly on the ground surface at the centre of each subplot via both manual and automated dynamic closed-chamber flux systems. A 20 cm wide tarpaulin ''skirt'' was affixed to the base of chambers and weighed down using a length of chain, preventing leakage and air exchange and minimising disturbance of soil and vegetation (Street and others 2007). Chamber headspace air was mixed using a low-speed fan and circulated through an infrared CO 2 analyser via Bev-A-Lineä tubing at 0.9 LPM before returning to the chamber. Flux was estimated from the regression slope defining change in chamber headspace CO 2 concentration over the first 30 s of each measurement set. We avoided measurements when winds were gusting and speeds were more than about 10 km h -1 to avoid turbulence effects on fluxes. Manual FCO 2 measurements were assessed using a portable single-chamber system equipped with a LiCOR LI-830 (Licor Biosciences, Lincoln NE, USA) CO 2 analyser. Two consecutive sets of CO 2 measurements spanning 45 s apiece were logged at a rate of one-per-second using the LI-8x0 software (LiCOR Biosciences, v.1.0.2). During the first set, the chamber was covered with an opaque cover to block sunlight and prevent photosynthesis, assessing only soil and plant respiration (''Dark''). The cover was then removed prior to the second set, allowing sunlight to enter the chamber (''Sun'') to assess instantaneous net ecosystem exchange (NEE). We defined positive NEE as representing a net efflux of C to the atmosphere from the ground surface. Manual measurements were collected weekly in 2018 between June and September (DOY 175-265) from the Dalton site only. In 2019, FCO 2 was measured biweekly from the Dalton and Steese sites between May and September .
Automated measurements were assessed using a non-portable flux chamber system equipped with a LiCOR LI-850 CO 2 analyser plumbed to eight chambers placed on each subplot within a plot and logged using a CR1000x datalogger (Campbell Scientific, Logan UT, USA). Pneumatically actuated lids closed for four minutes during flux measurements collected every 15 min in a clockwise sequence, allowing each chamber to be measured every 2 h. ''Dark'' equivalent measurements were collected during the night hours (00:00-3:00) while ''Sun'' measurements spanned daytime hours. Due to mobility constraints, this system was only used to measure flux from two specific plots at the Dalton site (one 1x and one 2x plot) for 2-3 consecutive days. Measurements were collected on two occasions for each plot between August to September in 2018 and three times between May to August in 2019. For analysis, all automated measurements were averaged to one dark and sun measurement per site visit. Mean daily values were similar to manual flux measurements.

Partitioning of C Fluxes
The difference between dark and sun measurements of FCO 2 for each subplot were used to estimate GPP ( Figure 2). Intact plots accounted for FCO 2 from both total respiration of soil (roots + heterotrophic organisms) and ground-level vegetation (RS GV ), while trenched plots excluded roots and vascular vegetation, accounting for FCO 2 from only heterotrophic organisms and non-vascular plants (mostly mosses) and lichen (RH NV ). Similarly, estimates of GPP from intact subplots assessed all ground-level vegetation (GPP GV ; vascular + non-vascular), while trenched subplots assessed GPP from non-vascular vegetation and lichens only (GPP NV ).
Trenching and removing vegetation from plots may introduce sources of error into measurement: RH may increase from root decay, lack of fresh root exudates and new growth may reduce RH, and reduced evapotranspiration rates and ground-cover from ground vegetation may influence water availability and soil temperatures (Bond-Lamberty and others 2004a). However, comparisons of Shortened Fire Intervals Stimulate Carbon Losses from Heterotrophic Respiration studies that use trenching versus different methods to partition RH from RS show minimal differences in results, indicating that the method is reasonably robust (Bond-Lamberty and others 2004b).

Ancillary Measurements
During flux measurements, in-chamber air temperatures were measured using an Omega type T thermocouple probe within chambers. Soil temperature and moisture values for subplots were measured to a depth of 10 cm using an Omega type T thermocouple and thermometer and a Campbell Scientific Hydrosense II (Campbell Scientific, Logan UT, USA), allowing values to stabilise over 2 min. For manual measurements only, photosynthetically active radiation (PAR) was measured inside the chamber above vegetation using an Apogee SQ-616 PAR meter and logged using Apogee Connect software (Apogee Instruments, Inc., N Logan UT, USA).
Soil and air temperatures were also recorded hourly from permanent weather stations installed in a single representative plot for each burn frequency at each site. Weather stations were comprised of thermocouples measuring soil temperature at 10 cm and air temperature at 30 cm above the ground connected to Campbell Scientific dataloggers (CR10x). In 2018, only the DA-1x and DA-2x plots had weather stations. In 2019, weather stations were in all burn frequencies at both sites but not in the mature forest areas. To supplement missing data, daily mean soil temperatures at 10 and 50 cm depths for each burn frequency and mature forest were also provided by Romanovsky and others (2020; unpublished results) using a U30-NRC HOBO data logger and S-TMB-M temperature sensors (Onset Computer Corp, Cape Cod, Massachusetts) calibrated using an ice bath, improving the accuracy for temperatures near 0°C to ± 0.02°C. More information on the methods can be found in Romanovsky and others (2020).

Data Analysis and Modelling
Differences amongst burn frequency and site for RS GV , RH NV , GPP GV , and GPP NV were tested in a linear mixed model (LMM; Tables 2, S2) framework using R version 4.1.1 (R Core Team 2020) and the LME4 package (v1.1-27.1; Bates and others 2015). Preliminary analyses produced heteroscedastic residuals with a right-tail distribu- tion and extreme values in both tails. To mitigate this, we employed log or inverse hyperbolic sine (IHS) transformations where appropriate to meet assumptions for subsequent mixed model analyses.
The IHS function performs similarly to a log transformation but allows for use of negative values and is less sensitive to values approaching zero (Burbidge and others 1988).
Even after transformation, a handful of extreme values introduced heteroscedasticity into models. Since these were ''real'' values produced by processes outside those measured in the study, values were ''windsorized'' to 99% of the distribution using the R package DescTools (v0.99.42; Andri and others 2021). This improved the homogeneity of variance in model residuals and reduced inflated effect sizes while retaining all measurements.
Adjusted log-likelihood scores (AIC) were used to aid in random effect selection. A random intercept was included for subplots to control for repeated measures. Plot was tested as a random effect but ultimately within-plot variance (that is, between subplots) was not significantly different from between plot variance (that is, spatial variance of FCO 2 at the subplot level equalled that between plots) and inflated the AIC score.

Flux Response Models
In order to estimate the cumulative annual C flux from each burn frequency over the entire snowfree period (DOY 135-285), we created statistical models that predict the response of respiration and GPP to soil and air temperature and PAR (flux response models). Hourly estimates of flux for each subplot were summed, accounting for diel patterns of flux associated with changes in temperature and PAR. To achieve this, we first predicted hourly Shortened Fire Intervals Stimulate Carbon Losses from Heterotrophic Respiration estimates of soil temperature, air temperature, and PAR for each subplot between DOY 100-300 using a LMM with manually measured values as a dependent variable and continuous weather station temperature values as predictors. Our in-plot weather stations provided soil and air temperatures, while PAR was estimated from solar radiation data collected from nearby regional weather stations. Predictions were bootstrapped (500x) to provide robust estimates and standard errors.
To ensure our models were not overpredicting flux from colder temperatures outside the range and period of measured values (DOY < 140 and > 265), we included estimates of FCO 2 measurements for each subplot at soil/air temperatures and PAR values predicted for DOY 105 and 295 (Ts 10 : -2 to -5°C). Soil respiration rates were estimated by a LMM with a random intercept for subplots using IHS-transformed soil respiration and temperature data pooled between our study and values aggregated by Natalia and others (2008) from studies of similar upland boreal forests during the winter. No relevant literature values were available to provide seasonal end-points for GPP, so we estimated these using a general linear model that pooled the IHS-transformed GPP dataset across all grouping levels into two categories: intact and trenched plots. Untransformed air temperature and log-transformed PAR values were used as dependent variables to predict values for each subplot.
Flux response models for RS GV and RH NV were developed using IHS-transformed dark chamber FCO 2 measurements as the dependent variable and the IHS of soil temperature at 10 cm (TS 10 ), air temperature (TA), and burn frequency (as a factor) as fixed covariate predictors (Tables 3a, S3a). Models for GPP GV and GPP NV were developed using IHS-transformed GPP estimates as the dependent variable and TA, the natural log of PAR, and burn frequency (as a factor) as fixed covariate predictors (Tables 3b, S3b). Model estimates were bootstrapped (500x) to produce robust error estimates and predictions.

Temperature Sensitivity (Q 10 ) Models
The response of RS GV and RH NV to changes in soil temperature over the growing season were also evaluated as a Q 10 function, defined as the rate of change in flux over a 10°C range in soil temperature. The Q 10 function provides a widely recognised metric to compare the apparent temperature sensitivity of soil substrates and vegetation amongst groups. Values were calculated using the b 1 coefficients from LMMs defined by a relationship be-tween log-transformed RS GV or RH NV values and IHS-transformed TS 10 values (Eq. 2; Tables 4, S4) using the emtrends function of the emmeans package in R (v1.6.3; Length 2020). These models focus on apparent soil temperature sensitivity alone, and differ from the flux response models above which included other covariates.
A full list of acronyms used in this paper can be found in the supplementary materials.

Site and Vegetation Characteristics
Ground-cover varied considerably across burn frequencies and exhibited a patchy mosaic of mosses, leaf litter, lichen, and exposed mineral or organic soil ( Figure 3A) and clusters of vascular ground vegetation (Figures 3B, S2; Table S5); < 30 cm in height). Mature plots at both sites were underlain by a continuous cover of mostly pleurocarpous mosses with minimal litter, no bare soil, and 0.5-0.7 Mg C ha -1 of vascular ground vegetation biomass. In turn, burns were mostly covered by a variety of acrocarpous mosses, with only ST-1x plots retaining any pleurocarps (% 5%). At the Dalton site, acrocarp cover increased with burn frequency from 33% in DA-1x plots to 79% cover in DA-3x plots. The opposite pattern was observed for vascular ground vegetation: biomass in DA-1x plots was nearly double that of mature plots but declined with burn frequency to less than 0.1 Mg C ha -1 in DA-3x plots. At the Steese site, acrocarp cover was continuous in ST-1x plots (72%) and slightly lower and more heterogeneous in ST-2x plots (65%). Both ST-1x and DA-1x plots were accompanied by the highest vascular ground vegetation biomass recorded in the study (1.0-1.2 Mg C ha -1 ). The ST-3x plots were the least vegetated plots at ground-level, containing only 18% acrocarp cover and less than 0.1 Mg C ha -1 of vascular ground vegetation, but was covered by a large amount of litter (72%). Exposed soil was only significantly present at the DA-1x (47%) and ST-2x plots (27%).
Mature plots at both sites contained the highest overstorey basal area in the study (% 25 m 2 ha -1 ), covered almost exclusively by black spruce (Table S5). At the Dalton site, 1x plots lacked an overstorey but contained a tall understorey (Table S5; Figure 3B) comprised by roughly equal amounts by black spruce and mixed-deciduous regeneration (white birch, willows, and alder). A mostly deciduous tall understorey in DA-2x plots filled gaps between clumps of an emerging mixeddeciduous overstorey (3.1 m 2 ha -1 ) that was better established in DA-3x plots (5.3 m 2 ha -1 ), predominantly by paper birch. Vegetation communities amongst burn frequencies differed at the Steese site; compared to DA-1x plots, ST-1x plots retained a small surviving mature spruce overstorey component (2.4 m 2 ha -1 ) and contained twice as much tall understorey black spruce regeneration. Steese 3x plots had an aspen-dominated overstorey with the highest basal area of any burn (9.2 m 2 ha -1 ), twice that of DA-3x plots. Total live root biomass in the mineral soil to a 45 cm depth was substantially lower in 1x plots (0.9-1.7 Mg C ha -1 ) than the mature plots (2.3-2.4 Mg C ha -1 ) at both sites ( Figure 3C; Table S6). Values were significantly higher in both the DA-2x (1.9 Mg C ha -1 ) and DA-3x (3.6 Mg C ha -1 ) plots than DA-1x plots, while root densities were similar between ST-1x and 2x plots and highest in ST-3x plots (4.6 Mg C ha -1 ).  The SOL depth halved with successive fires at the Dalton site, declining from 17 cm in mature plots to < 1 cm in DA-3x plots ( Figure 3D; Table S6). A similar pattern at the Steese site showed less change in ST-2x plots (6.4 cm) compared with DA-2x plots (4.6 cm), but 3x plots reached the same minimal depth for both sites. Increasing depth of the SOL was a strong predictor of increasing vascular ground vegetation biomass and declining deciduous overstorey BA in burned plots (r 2 = 0.71/0.53, p £ 0.001; Figures S4a, S4b; Table S7). In turn, increasing deciduous overstorey BA was a strong predictor of higher live root biomass in the mineral soil and leaf litter cover (r 2 = 0.86/0.51, p £ 0.001; Figures S4c, S4d; Table S7).

Soil Moisture, PAR, and Temperature
Seasonal volumetric soil moisture values (Figure 4A; Table 5) were at their peak in mid-May after snowmelt, ranging from 11 to 33%) in all burns and mature plots. Values declined rapidly by mid to late June (3-12%) and tended to recover slightly by late-summer to autumn (5-15%). Soil moisture values in DA-2x, DA-3x, and ST-3x plots were about half that of other burns after initial decline and changed minimally until autumn, while values in the mature plots continued to decline, reaching values similar to 3x plots in autumn. At both sites, 3x plots had slightly higher bulk densities and a very thin SOL, indicating reduced pore space available for moisture storage and a lack of insulating cover to reduce evaporative losses. Although soil moisture did not differ in general between sites, mineral soils at the Steese site had higher bulk densities (+ 11-32%; Table S6) than the Dalton site. Soil at the flatter Steese site also showed moderate mottling within 20 cm of the surface in all burns, particularly in ST-1x and ST-2x plots, but not at the more highly sloped  PAR and soil temperature are normalized between continuous weather station and site level instantaneous measurements (see ''Methods'' section). lower in 3x plots (9-11 mol m -2 day -1 ) where shading from a developing deciduous overstorey reduced insolation.
Mean daily soil temperatures at 10 cm depth (TS 10 ) increased rapidly after snowmelt and plateaued between July and mid-August, declining thereafter ( Figure 4C). Peak temperatures were lowest in the mature plots (9-10°C) and increased consistently with burn frequency to % 13°C in 3x plots. Mean TS 10 values in burns (Table 5), averaged by plot over the growing season, were strongly correlated with declining SOL depths (r 2 = 0.66-0.80, n = 20, p £ 0.001; Figure 5A; Table S8). Even greater differences in soil temperatures emerged with soil depth: the difference in mean annual soil temperature at 50 cm depth (TS 50 ) during the growing season between reburns and 1x plots (Table 5; 3-6°C) was comparatively much higher than the differences for TS 10 (2-3°C). The TS 50 of the mature plots remained below 0°C and soil sample cores retrieved during late-summer frequently contained frozen soil and ice lenses between 30 and 50 cm depths, indicating the presence of permafrost. Dalton 3x plots had surface temperatures (TS 10 ) similar to those in DA-2x plots but deeper soil temperatures (TS 50 ) that were similar to those from DA-1x plots. This was likely due to the more northerly hillslope aspect and hillslope shading of the plot, which received 10% less estimated incoming annual solar radiation relative to the mature plots at the Dalton site (Points Solar Radiation Tool, ArcGIS 10.3).

Differences in Fluxes Amongst Burn Frequencies
Measured rates of FCO 2 ( Figure 6; Table 6) were highest in the mature plots (RS GV : 3.3-8.4 mg C m -2 min -1 ; RH NV : 1.9-3.8 mg C m -2 min -1 ), declined significantly after the first burn (RS GV : 1.1-4.7 mg C m -2 min -1 ; RH NV : 0.6-2.9 mg C m -2 min -1 ), and increased substantially in reburns (RS GV : 0.9-6.7 mg C m -2 min -1 ; RH NV : 0.9-4.7 mg C m -2 min -1 ). Values were lowest in spring and autumn and highest in summer for all burn frequencies. Burns with the lowest respiration rates (DA-1x and ST-2x plots) also expressed the least seasonal variation (RS GV : CV = 23-28%; RH NV : CV = 10-23%) and contained no deciduous overstorey vegetation, while those with the highest rates (DA-2x plots and 3x plots) had the highest seasonal variability (RS GV : CV = 40-47%; RH NV : CV = 30-33%) and deciduous overstorey BA. The mature plots also had low seasonal variability (RS GV : CV = 19-26%; RH NV : CV = 18-19%) despite high FCO 2 rates. This pattern was stronger for Data aggregated by site and burn frequency. The standard error of the mean (SEM) was estimated using the standard deviation of the mean annual values for each subplot.
Shortened Fire Intervals Stimulate Carbon Losses from Heterotrophic Respiration RS GV than RH NV , presumably since RS GV is more sensitive to over/understorey phenology via inclusion of autotrophic respiration. Overall mean differences in respiration between sites were small (Tables 2, S2) and largely attributable to high rates in ST-3x plots.
Rates of RH NV averaged to the plot scale over the growing season were positively correlated with plot scale measurements of overstorey deciduous BA, live root biomass, and % leaf litter cover (r 2 = 0.55/ 0.52/0.38, p £ 0.01; Figure 5A; Table S8), and negatively correlated with SOL depth (r 2 = 0.36, p £ 0.01). Associations were generally similar but weaker for RS GV . Q 10 was also a very strong predictor of both RS GV and RH NV (r 2 = 0.81/0.86, p £ 0.0001; Figure 5B; Table S8).
Relative to FCO 2 , GPP GV ( Figure 6) in burns peaked earlier in the season (July) and tapered more rapidly towards September. Peak values were highest in 1x plots (3.1-3.2 mg C m -2 min -1 ) and declined in 2x plots (2.2-2.5 mg C m -2 min -1 ) and especially 3x plots (1.1-1.2 mg C m -2 min -1 ). In mature plots, values were more consistent throughout the growing season, with lower peaks (1.7-1.9 mg C m -2 min -1 ). In comparison, GPP NV was much lower than GPP GV , varied little over the growing season, and was generally similar amongst burns (0.4-0.8 mg C m -2 min -1 ). Mean GPP GV averaged to the plot scale over the growing season was strongly positively correlated to SOL depth and  Table S8) and negatively correlated to overstorey deciduous vegetation and live root biomass (r 2 = 0.75/0.64, p £ 0.0001). GPP NV was significantly reduced by litter cover (r 2 = 0.38, p £ 0.01).
All linear mixed models showed a significant effect for burn frequency for all C fluxes and Q 10 values (Tables 2, 3, 4). Both RS GV and RH NV showed a significant interaction effect between site and burn site-based differences (Table 2a). This was largely due to differing responses for burns, with 1x and 3x plots having higher respiration values at the Steese site. In the flux response models, site-based differences became insignificant, except where interacting with soil temperature and burn frequency (Table 3).

Flux Response Models and Estimated Annual Fluxes
Soil temperature at 10 cm depth (TS 10 ), air temperature (TA), and PAR were all important predictors of FCO 2 in our flux response models (Table 3). Seasonal patterns of predicted daily FCO 2 values ( Figure S3) closely mirrored observed values between DOY 140-260 (r 2 : RS GV /RH NV = 0.77/ 0.73, GPP GV /GPP NV = 0.53/0.26), while values predicted outside the measured range were similar to those used from other studies to inform the models (Natali and others 2008). Using the IHS function to transform both the dependent variable and primary independent covariate (TS 10 ) prevented the over-estimation of values typical of exponential models when soil temperatures approached zero (Bond-Lamberty and others 2004a), and unlike log transformations returned real numbers when temperatures were negative. We estimated the annual cumulative flux of C from RS GV , RH NV , and GPP by summing hourly flux predicted for each subplot during the snow-free period (DOY 135-285; Table S9), excluding winter estimates. Total annual FCO 2 efflux in mature plots (RS GV : 9.4-9.8 Mg C ha -1 ; RH NV : 4.9-5.2 Mg C ha -1 ) declined substantially after the first burn (RS GV : 3.5-6.2 Mg C ha -1 ; RH NV : 1.9-3.2 Mg C ha -1 ) and was highest in 3x plots (RS GV : 5.9-6.9 Mg C ha -1 ; RH NV : 3.3-5.4 Mg C ha -1 ), particularly for RH NV . Conversely, total annual GPP GV was highest in 1x plots (2.6-2.8 Mg C ha -1 ) and declined substantially in 3x plots (0.8-1.2 Mg C ha -1 ). Despite much lower PAR values and similar vegetation communities, annual GPP GV in the mature plots (1.7-2.4 Mg C ha -1 ) was similar to 2x plots (2.1-2.2 Mg C ha -1 ), indicating that other site factors were more important to differences amongst burn frequencies. Annual GPP NV was lowest in ST-3x plots (0.3 Mg C ha -1 ) but did not differ significantly amongst other burns (0.4-0.7 Mg C ha -1 ). Values were consistently higher in the mature plots (0.7-0.8 Mg C ha -1 ) where moss cover was highest, highlighting the importance of mosses to ground-level understorey GPP under the low-light conditions of black spruce forest canopies.

Differences in Environmental Variables and C Flux Within Burn Frequency Between Sites
Overall, patterns of vegetation, SOL thickness, root biomass, soil temperatures, and FCO 2 were similar between sites for mature plots and 3x plots (Figures 3, 4, 5, S2; Tables S3, S5, S6). Inter-site differences for the same burn frequency were present for 1x and 2x plots, potentially due to differences in burn severity in the most recent fires (14-16 years ago). Despite sharing a similar loss of SOL material and tall understorey, ST-1x plots retained some mature spruce and a near-continuous moss cover, while DA-1x plots showed greater declines in respiration, complete overstorey mortality, and widespread patches of exposed soil. Steese 2x plots more closely resembled DA-1x plots: no overstorey and Data aggregated by site and burn frequency. SEM = the standard error of the mean. Post-hoc tests use a compact letter display to denote significant differences (p £ 0.05) amongst burn frequencies within sites using a Holm-Bonferroni correction for multiple comparisons.
significant spruce regeneration, analogous cover of ground vegetation and exposed soil, and similar rates of respiration and GPP GV . In comparison, DA-2x plots had a thinner SOL, a developing deciduous overstorey, and similar respiration rates and root biomass as DA-3x plots. Near-surface soil temperatures (TS 10 ) were nearly identical in DA-2x plots and 3x plots at both sites. However, lower TS 50 values in DA-3x plots were associated with slower growing paper birch, while higher TS 50 values in ST-3x plots were coupled with a faster growing aspen overstorey and higher RH NV rates and Q 10 values. Despite these differences, general patterns of C flux for burn frequencies emerged when predictions from the flux response models were pooled between sites (Figure 7; Table S9). Soil respiration was halved after the first burn and subsequently increased with burn frequency. This pattern was stronger for RH NV , increasing by 29% in 2x plots and 80% in 3x plots, while RS GV remained similar between 1x and 2x plots but increased by 40% by the 3 rd burn. The ratio of RH NV :RS GV for cumulative annual flux was % 52% in both mature plots and 1x plots and increased to % 69% in both 2x and 3x plots. Changes to GPP from the ground surface displayed the opposite pattern; GPP GV values increased by 26% after the first burn and declined to their lowest rates in 3x plots. GPP NV declined by after the first burn and then again slightly by the 3rd burn, and comprised an increasing proportion of GPP GV with burn frequency when averaged across sites (1x = 18%, 2x = 26%, 3x = 48%).
All burns displayed a net positive annual exchange of C from the ground surface to the atmosphere (NEE GS ; Figure 7; Table S9) that was lowest in 1x plots (2.1 Mg C ha -1 ). Values increased with burn frequency to 2.4 Mg C ha -1 in the 2x plots and 5.5 Mg C ha -1 in 3x plots, approaching those from the mature plots (7.7 Mg C ha -1 ).

DISCUSSION
Our study results support the hypothesis that increased fire frequency in boreal forest regions has accelerated mineral soil C loss via heterotrophic respiration. Reburned areas had less SOL material to insulate the mineral soil, and thus the soil tended to be warmer and emit more C as annual soil CO 2 efflux. This tendency differed somewhat between our two sites; at the Dalton site, RS GV and RH NV increased substantially in both the 2x/3x plots relative to 1x plots. At the Steese site, RS GV and RH NV were lower in the 2x than the 1x plots, while the 3x plots had the highest values of any burn at either site. Compared across both sites, RS GV and especially RH NV tended to increase with burn frequency. These trends mostly followed near-surface and deeper soil temperatures that also increased substantially with burn frequency, especially in the 3x plots, and higher soil temperatures were correlated to decreasing SOL depth (Figure 5A).
Ground vegetation was a notable modifier of soil respiration and GPP across the burn sequences. At both sites, GPP GV values were much higher in the 1x and 2x plots than in the 3x plots, nearly offsetting C losses (RS GV ) in DA-1x and ST-2x plots ( Figure 6; Table S9). These differences were primarily due to greater ground vegetation coverage in the 1x and 2x plots, which was correlated to increasing SOL depth ( Figure 5A). Consequently, twice-burned areas have the potential to be similar to areas burned either once or three times, but the occurrence of a third burn within the past 60 years appeared to be a critical threshold to triggering an Figure 7. Model predicted cumulative annual flux of CO 2 by soil and vegetation for each subplot over the snow-free period (DOY 135-285) amongst burn frequencies for RS GV , RH NV , GPP GV , and GPP NV (Flux), and net ecosystem exchange for all ground-level vegetation and soil (NEE GS ). Solid lines represent changes to mean C efflux from RS GV , RH NV , and NEE GS while the dashed lines represents changes to mean C influx for GPP. Line ribbons represents ± 1 SEM. Dots fitted using horizontal jitter for enhanced visualisation. Triangles represent intact subplots while circles represent trenched subplots. elevated soil respiration response and inhibiting establishment of vascular ground vegetation. This reburn-induced threshold shift in ground-level understorey GPP and soil respiration dynamics can be attributed to changes in both the thermal dynamics and biogeochemical properties of the underlying soil and to the overlying vegetation communities, including both ground-level and overstorey vegetation cover and associated root biomass.

Soil Temperature as a Driver of Elevated Soil Respiration Rates in Reburns
Many studies of artificial soil warming in forests credit increases in FCO 2 rates to temperature enhanced decomposition (Eliasson and others 2005; Coucheney and others 2013; Melillo and others 2017; Lim and others 2019). In our flux response models, near-surface soil temperature at 10 cm (TS 10 ) over the growing season was an effective predictor of FCO 2 (Table 2), accounting for 43-47% of variance in FCO 2 . Relative to 1x plots, reburned plots had higher TS 10 (+ 2-3°C) and RH NV values (+ 28-80%), suggesting that elevated soil temperatures enhanced decomposition rates. However, differences in soil temperature only partially explained differences in RH NV amongst burn frequencies. The interaction term between TS 10 and RH NV in the flux response model accounted for only 7% of total variance, and model estimates of RH NV made by equalising soil and air temperatures in all burn frequencies to those from 1x plots only reduced the differences between reburns and 1x plots by 10-23%. Differences in air temperatures amongst burn frequencies were minimal, and soil moisture, another common covariate, was a weak overall predictor and failed to account for differences in RH NV . Furthermore, SOL losses in reburns reduced potential decomposer contributions from a large C pool otherwise available in 1x plots. Consequently, factors other than near-surface soil temperatures were contributing to elevated RH NV in reburns.
An additional mechanism to account for elevated RH NV in reburns may have originated from deeper soils. The difference in mean soil temperatures between warmer reburns and cooler 1x plots (Table S5) was twice as large at 50 cm depths (TS 50 ) than at 10 cm (TS 10 ). Additionally, Q 10 values increased with fire frequency and were strongly correlated with RH NV (Figure 5A), but not nearsurface soil temperature. Over a growing season, the Q 10 of surface fluxes reflects both plant phenology and the changing availability of soil C to microbes at different soil depths (Davidson and others 2006). Graf and others (2008) observed that Q 10 increased substantially with soil depth, potentially as a response to increased temperature sensitivity of deeper SOC to decomposition (Davidson and Janssens 2006). In reburns, greater Q 10 values for RH NV and a higher ratio of RH:RS than mature and 1x plots may have reflected the integration of more temperature sensitive C and warmer deep soil temperatures.
Using an empirical model based on field measurements, Gibson and others (2019) predicted that increased active layer depth and soil warming after fire in a black spruce peatland raised RH rates in deeper soils by four-fold. O'Neill and others (2002) also found warmer deep soil temperatures and greater Q 10 values in aspen forests than black spruce forests seven years after a fire. Our results suggest that burn frequency accentuated this burnsoil warming tendency, both from a reduction in the SOL layer depth during combustion and a shift in vegetation toward angiosperms. For example, higher RH NV rates and Q 10 values in DA-2x plots than DA-3x plots may have been due to warmer mean growing season soil temperatures at 50 cm (8.4°C vs. 4.0°C), even though TS 10 values were similar for both (% 10.5°C). However, this deep soil temperature relationship to FCO 2 did not hold for other burns. Rates of RH NV in ST-2x plots were similar to DA-1x plots and much lower than ST-1x plots ( Figure 6; Tables 3, S9), despite being warmer than both at 50 cm (6.7°C vs. 3.0-3.4°C; Table 5). Furthermore, RH NV rates declined substantially after the first burn at both sites despite warmer TS 50 values (+ 3°C) and deeper active layers than in mature plots. Together, these observations suggest that in addition to deep soil warming, other reburn-specific factors such as vegetation species composition and biomass contributed to an elevated RH NV response.

Vegetation Drivers of Soil C Flux
Increased burn frequency appeared to alter the successional trajectory away from black spruce and towards an aspen or birch forest, with potential effects on SOC dynamics caused by both temperature and vegetation type. Laganiè re and others (2012) also found warmer soil temperatures in aspen and mixed-wood stands insufficient to account for higher RH rates relative to a black spruce forest, and attributed differences to faster turnover of SOC stimulated by higher-quality aspen litter inputs. Developing upland deciduous boreal forests can increase the soil supply of labile C and plant-available nutrients relative to a developing spruce forest (Flanagan and Van Cleve 1983;Alban and Pastor 1993;Ruess and others 2003;Johnstone and others 2010b). Deciduous forests can also produce twice as much higher-quality leaf-litter and contain three times greater live fine root densities that extend deeper into the mineral soil (Ruess and others 1996;Steele and others 1997;Laganiè re and others 2010). Plant roots, particularly from deciduous trees, have been associated with an increase in fungal decomposer activity and SOC decomposition rates in boreal forests (Bradley and Fyles 1995;Adamczyk and others 2019). Increased labile C and inorganic N inputs from litter leachates and root exudates can stimulate decomposition of mineral SOC via microbial priming (Bradley and Fyles 1995;Lé garé and others 2005;Fontaine and others 2007; Laganiè re and others 2010), particularly in the active layers of high-latitude ecosystems (Wild and others 2016;Keuper and others 2020) and in mineral soils (Karhu and others 2016).
In our study, reburn plots with the highest RH NV rates, Q 10 values, and soil temperatures also contained the greatest biomass of roots in mineral soil and deciduous overstorey vegetation. Positive relationships between these variables ( Figure 5) suggest that a developing deciduous overstorey and a deeper and warmer active layer in reburns were encouraging root growth and increasing RH, likely by expanding C substrate availability and stimulating the decomposition of deeper SOC. This was particularly evident in ST-3x plots, which had an aspen-dominated overstorey, the warmest TS 50 values, and the highest root biomass and RH rates in the study. Models of rhizosphere priming in thawed permafrost soils indicate that increased plant root production might amplify soil respiration rates by 12% (Keuper and others 2020). This suggests that root priming may have contributed to higher RH NV rates in reburns, and that changes to forest cover and C substrate availability from reburning may influence the expression of this process across the permafrost region.

Re-burning Alters GPP and NEE from Ground-Level Vegetation
Burn frequency altered vegetation communities and net soil C flux, expressed as a large decrease in both the biomass and GPP of ground vegetation in reburns ( Figures 2C, 6), primarily as a function of changes to vascular vegetation ( Figure 5C). Mature plots contained a near-continuous pleurocarp moss cover, with GPP NV rates 33% higher than in burned plots and similar to those observed in a mature black spruce forest by Gaumont-Guay and others (2009). Values differed minimally amongst acrocarp-dominant burns except for ST-3x plots, which had substantially more ground covered by leaf litter than other burns ( Figure 3A). This is consistent with observations by Jean and others (2017), who found acrocarp colonisers to be dominant in both spruce and deciduous early successional boreal forests and higher leaf litter cover in aspen forests such as those in ST-3x plots.
Relative to mature forest, 1x plots experienced considerable increases to both annual GPP GV flux (+ 35%) and the biomass of ground-level vascular vegetation (+ 71%; Tables S2, S7), potentially linked to greater light availability (Mean PAR: % 285 vs. % 78 lmol m -2 s -1 ; Figure 4B; Table 5). However, GPP GV and biomass declined somewhat in 2x plots, despite greater light levels (% 360 lmol m -2 s -1 ); even 3x plots had higher light levels (% 130 lmol m -2 s -1 ) than mature plots, yet less than half the primary productivity rates. The response of GPP GV to PAR in our study plateaued around 250-300 lmol m -2 s -1 (Figure S5), similar to levels noted by Kolari and others (2006) for boreal dwarf shrubs, reflecting an adaptation to partial shade. The plateau for bryophytes (GPP GV ) was even lower (% 100 lmol m -2 s -1 ), indicating that photosynthesis was effectively saturated at PAR levels found in 1x plots. Declining in GPP GV rates in burns were more strongly correlated with an increase in deciduous overstorey basal area and root biomass (r 2 = 0.75 and 0.63, p < 0.001; Figure 5; Table S8). This suggests that while increased insolation in 1x plots contributed to greater primary productivity and accrual of biomass at the ground-level, competition for resources with trees and taller shrubs was likely a larger contributing factor amongst burns.
Other studies estimating GPP from ground-level vegetation using flux chambers in the boreal forest have generally focused on mid-successional and mature boreal forests (Goulden and Crill 1997;Heijmans and others 2004;Kolari and others 2006;Gaumont-Guay and others 2009). More apt comparisons with early successional sites can be found in studies using eddy-covariance flux towers. Goulden and others (2011) estimated changes to total ecosystem GPP (GPP ECO ) in areas burned at various intervals in a post-fire chronosequence study of Manitoba boreal forest. The 6-year postfire site (GPP ECO = 3.8 Mg C ha -1 year -1 ) contained no overstorey and a ground vegetation cover similar to our 1x plots (GPP GV = 2.6 Mg C ha -1 year -1 ). Similarly, flux-tower estimates of GPP ECO by Kolari and others (2004) in a 4-year-old Shortened Fire Intervals Stimulate Carbon Losses from Heterotrophic Respiration Scots pine clear-cut without an overstorey (2.9 Mg C ha -1 year -1 ) were nearly identical to our 1x burn. Although these systems likely differ in some ways from our once-burned areas, a comparison among results suggests that GPP GV alone comprises the majority of GPP ECO during early boreal forest succession. In contrast, a decline in GPP GV in 2x plots (2.0 Mg C ha -1 year -1 ) and especially 3x plots (0.9 Mg C ha -1 year -1 ) highlights that reburning may substantially reduce uptake by ground-level vegetation communities.
These reburn-specific differences in ground vegetation succession dynamics may be instrumental to SOL re-development. Both GPP GV and the biomass of vascular ground vegetation were more strongly positively correlated to SOL depth than soil respiration and highest in 1x plots (Figures 5A,  S4; Tables S7, S8). Assuming % 55% of annual GPP by vascular ground vegetation alone (GPP Vasc = GPP GV -GPP NV ) is returned to the atmosphere during autotrophic respiration (Campioli and others 2015), approximately 1.0 Mg C ha -1 is allocated to new biomass annually in 1x plots (that is, net primary production), similar to the existing aboveground biomass present in 1x plots (Figure 3C; Tables S5, S9). Combined with high annual turnover rates observed for boreal dwarf shrubs (Nilsson and Wardle 2005), this suggests that turnover of roots and litterfall from vascular ground vegetation accrued during early post-fire succession may play a substantial role in feeding SOL regeneration. Conversely, the contribution of non-vascular vegetation to overall GPP from ground vegetation was much lower than from vascular vegetation: the ratio of GPP Vasc :GPP NV was % 2:1 in mature plots, % 5:1 in 1x plots, % 3:1 in 2x plots, and % 1:1 in 3x plots.
Pleurocarpous mosses are often identified as a primary component contributing to SOL development in spruce forests (Van Cleve and others 1983; Natalia and others 2008; Alexander and Mack 2016), but our results suggest that enhanced growth of vascular ground vegetation in onceburned areas (that is, forest burned at historically ''typical'' intervals) may be equally important to developing the SOL. A decline in this capacity by reburning reduces the potential for the ground vegetation component to contribute to NEE and SOL development. While all burns had a positive net efflux of C from the ground surface (NEE GS : that is, soil + ground-vegetation; Figure 7), rates were lowest from 1x plots (2.0 Mg C ha -1 year -1 ) and similar to flux-tower derived NEE rates estimated by Kolari and others (2004) in a 4-year-old Scots pine clear-cut with no overstorey (2.2 Mg C ha -1 year -1 ). Reburning increased NEE GS , with values peaking in 3x plots (4.7 Mg C ha -1 year -1 ) due to high RH NV rates and minimal vascular ground vegetation biomass. By comparison, our estimate of NEE GS rates in mature plots (5.9 Mg C ha -1 year -1 ) were also similar to that estimated by Kolari and others (2006) using flux chambers in 40-year-old Scots pine forest (5.3 Mg C ha -1 year -1 ), though both our GPP and RS values were substantially higher.
The post-fire recruitment of deciduous trees and shrubs generally increases with SOL consumption, which is positively correlated to fire intensity (Boby and others 2010). This shifts competitive advantages away from spruce to angiosperms, with early colonizers often predicting dominant future stand dynamics (Johnstone and others 2010b; Alexander and others 2012). We found that an increase in deciduous overstorey biomass in burns was negatively correlated to SOL depth and vascular ground vegetation biomass (Table S8). Deciduous vegetation both more effectively shades out competition and may introduce allelopathic compounds that inhibit understorey growth (Natalia and others 2008). Although SOL retention following fire provides a competitive advantage to spruce for regeneration, it is unclear whether it actively promotes re-establishment of vascular ground vegetation. Conversely, it is unclear to what extent reduced establishment of ground vegetation in 3x plots was due to SOL loss or from competition with a developing deciduous tree and shrub component. Numerous factors influence post-fire succession of overstorey and understorey vegetation, including loss of seed source and rhizome/root mass typically present in the SOL that may otherwise support resprouting after less intense fires (Schimmel and Granstrom 1996;Johnstone and others 2010b). Our study was not designed to assess this recruitment, or differences between fire frequency and intensity, and its role in forest NEE.
The nature of observational studies, such as ours, make them inherently vulnerable to confounding effects. Regional and local differences in pre-fire vegetation communities, drainage, and weather can influence both fire intensity and the likelihood of reburning at shorter intervals (Johnstone and others 2010b; Whitman and others 2018). For example, warmer soils in some reburned plots could have preceded recent fires, and promoted reburning relative to cooler areas. It is also likely that some of the differences in post-fire vegetation and SOL depth between sites in once and twiceburned areas were due to a lower intensity fire in the most recent burn at the Steese site, potentially due to wetter conditions. Effective assessment of fire severity using remote-sensed imagery is challenging for black spruce stands in hilly, high-latitude areas such as the interior of Alaska (Verbyla and others 2008), and the lack of precision and infrequency of older imagery makes it difficult to parse out pre-fire vegetation at smaller-scales necessary to our analysis, especially for multiple fires dating back over six decades.
As such, we cannot disentangle the potential influence of pre-fire conditions and fire intensity from burn frequency. Regardless, our study provides robust results from two distant field sites with similar fire histories and varying landscape conditions. These indicate that: (1) ground-level understorey vascular vegetation is highly productive in once and even twice-burned areas and likely important to near-surface C inputs; and (2) a third burn in a short period of time (for example, < 60 years) can substantially reduce the contribution of vascular ground vegetation to ecosystem GPP more than a decade after fire. Combined with elevated RH rates in reburns, this accentuates the potential for net ecosystem C losses. However, a study of post-fire chronosequence succession by Mack and others (2021) across the interior of Alaska found that a shift from mature spruce to deciduous forest can store up to fivetimes more C in aboveground biomass at maturity, potentially negating these early successional C losses. Their study considered the role of SOL loss but did not evaluate mineral soil C dynamics. Consequently, a better understanding is needed on how fire frequency and SOL consumption might interact to affect ground vegetation recruitment, SOL development, and NEE.

CONCLUSIONS
Short-interval reburning of high-latitude boreal forests increased rates of soil respiration, primarily as heterotrophic respiration, while decreasing ecosystem C uptake as GPP from ground-level understorey vegetation. These changes accompanied a substantial loss of the SOL and a shift from an anticipated successional pathway towards a spruce-dominated forest with a widespread pleurocarpous moss/perennial understorey towards a deciduous overstorey (birch or aspen) with minimal vascular understorey and patchy acrocarpous mosses. The most frequently burned areas also had the highest near-surface soil temperatures, but as a covariate, temperature differences did not adequately explain their greater C flux relative to the less frequently burned areas. Instead, patterns of C flux were better explained by differences in SOL depth, deciduous vegetation cover, and live root biomass.
Ultimately, our results indicate that boreal forest SOC decomposition rates are likely to rise both as an overall effect of rising temperatures associated with climate change and especially as a function of the reburned landscape. The historical paradigm of spruce predominance in the Alaskan interior boreal forest is rapidly shifting with a climate changedriven increase in fire frequency (Beck and others 2011;Mann and others 2012), while developing deciduous forests in reburns are re-partitioning C pool allocation between soil and vegetation components (Mack and others 2021). As such, future research investigating the response of ecosystem C dynamics to short-interval fires in high-latitude environments should focus on feedbacks between understorey vegetation and SOL development, vertical stratification of soil temperatures, and biogeochemical processes influencing labile C dynamics.
Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.