Patterns of Ecosystem Structure and Wildfire Carbon Combustion Across Six Ecoregions of the North American Boreal Forest

Increases in fire frequency, extent, and severity are expected to strongly impact the structure and function of boreal forest ecosystems. An important function of the boreal forest is its ability to sequester and store carbon (C). Increasing disturbance from wildfires, emitting large amounts of C to the atmosphere, may create a positive feedback to climate warming. Variation in ecosystem structure and function throughout the boreal forest is important for predicting the effects of climate warming and changing fire regimes on C dynamics. In this study, we compiled data on soil characteristics, stand structure, pre-fire C pools, C loss from fire, and the potential drivers of these C metrics from 527 sites distributed across six ecoregions of North America’s western boreal forests. We assessed structural and functional differences between these fire-prone ecoregions using data from 417 recently burned sites (2004–2015) and estimated ecoregion-specific relationships between soil characteristics and depth from 167 of these sites plus an additional 110 sites (27 burned, 83 unburned). We found that northern boreal ecoregions were generally older, stored and emitted proportionally more belowground than aboveground C, and exhibited lower rates of C accumulation over time than southern ecoregions. We present ecoregion-specific estimates of depth-wise soil characteristics that are important for predicting C combustion from fire. As climate continues to warm and disturbance from wildfires increases, the C dynamics of these fire-prone ecoregions are likely to change with significant implications for the global C cycle and its feedbacks to climate change.


INTRODUCTION
The boreal forest is one of the largest biomes on earth, covering almost 1.9 billion hectares and encompassing ∼30% of the global forested area (Brandt et al., 2013). From a global perspective, the most significant and critical function of the boreal forest is its ability to sequester and store carbon (C), as it contains approximately one third of terrestrial C stocks (Pan et al., 2011). This globally important biome is becoming increasingly vulnerable to change as the climate continues to warm.
One of the most rapid pathways through which boreal forests can be altered is with changes to the fire regime. Climate warming and drying has led to an intensification of boreal forest fires (Balshi et al., 2009;Flannigan et al., 2009), with large increases in the average annual area burned over recent decades (Coops et al., 2018). Further increases, up to a factor of five, are expected by the end of the century (Balshi et al., 2009;Boulanger et al., 2014). Increasing fire extent, frequency, and severity could alter boreal forest C storage (Bond-Lamberty et al., 2007;Walker et al., 2019): from the historical net accumulation of C from the atmosphere over multiple fire cycles, to a net loss, which in turn would cause a positive feedback to global climate warming (Oris et al., 2013;Li et al., 2017). Understanding the structure and function of these fire-prone ecoregions is needed to quantify the role of fire in the global C cycle and its feedbacks to climate change.
The accumulation of C in boreal ecosystems increases with stand age or time after fire, but the rate of accumulation depends on vegetation composition and environmental conditions that impact net primary productivity and decomposition (Luyssaert et al., 2008;Jonsson and Wardle, 2010). Black spruce [Picea mariana (Mill.) B.S.P.] is one of the dominant tree species in boreal North America and is generally less productive and slower growing than jack pine (Pinus banksiana Lamb.) or deciduous trees (Subedi and Sharma, 2013;Alexander and Mack, 2016). The majority of organic C sequestered in black spruce boreal forests resides in the soil organic layer (SOL; Ping et al., 2008;Johnson et al., 2011). Combustion of the SOL dominates C loss during fire, typically accounting for 80-90% of total C emitted in boreal forests of Alaska and Canada (Boby et al., 2010;Rogers et al., 2014;Walker et al., 2018b). Despite its importance, the greatest uncertainty in modeling C combustion from wildfire is estimating the SOL component (French et al., 2004), which requires information on SOL bulk density (g cm −3 ) and C fraction (%) to model C stocks (g C m −2 ). Both the rate of C accumulation and SOL characteristics are likely to be ecoregionspecific as a consequence of differences in long-term climate, geological and biogeographical history, soil development, and parent materials (Houle et al., 2017).
Fire severity can be defined as the loss of above-and belowground organic material that occurs as a direct consequence of fire (Keeley, 2009). Given that the majority of combustion takes place in the SOL (Boby et al., 2010;Rogers et al., 2014;Walker et al., 2018b), estimates of SOL burn depth and post-fire residual SOL depth have been used as metrics of fire severity (Greene et al., 2004(Greene et al., , 2007Johnstone and Kasischke, 2005;Johnstone and Chapin, 2006;Walker et al., 2018a). Total C combusted from fire is also a widely used metric of fire severity (Boby et al., 2010;Walker et al., 2018b), but variations in above-and belowground C combustion, C combusted relative to pre-fire C, and the proportion of belowground C combusted relative to total C combusted are also important metrics for understanding fire severity and its ecological impacts. In this study, we compiled data on pre-fire C pools, C loss from fire, and the potential drivers of these C metrics from 417 recently (2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) burned sites within 18 fires. We also collected depth-wise measurements of soil characteristics from 167 of these sites plus an additional 110 sites (27 recently burned and 83 control sites with no historical record of burning in the last 50-75 years). Sites were established by five independent research projects, sampled one year after fire, and distributed across six ecoregions of North America's western boreal forests (Figure 1 and Supplementary Table S1). Using this unique and spatially extensive dataset, we addressed the following questions: (1) How do pre-fire aboveground (stand age, density, basal area, proportion of black spruce, biomass, and C pools) and belowground (SOL depth and C pools) forest attributes differ among fire-prone ecoregions? (2) Do aboveground, belowground, and total C pools increase over time at a similar rate among ecoregions? (3) In the SOL profile, what is the relationship between depth and (i) bulk density, (ii) C fraction, and (iii) C stocks, and does this vary between ecoregions in a way that is dependent on total SOL thickness? (4) Do fire severity and C combustion (i.e., residual SOL depth, burn depth, above-and belowground C combustion, proportional combustion) vary among ecoregions?

Study Areas and Data Acquisition
We collated data from 527 sites established by five independent research projects located in six ecoregions in the western North American boreal forest (Figure 1 and Supplementary Table S1). Sites were established in the northern ecoregions of Alaska Boreal Interior, Boreal Cordillera, Taiga Plains, and Taiga Shield, and in southern ecoregions of Softwood Shield and Boreal Plains. Ecoregions differ in geological history, soils development, parent materials, and mean annual temperatures and precipitation (US EPA, 2015). We collated site-level data on pre-fire stand structure, stand age, topography, and pre-and post-fire aboveand belowground C pools from 417 sites across 18 wildfires. These sites were established by four of the five research projects (Supplementary Table S1), sampled one year after fire, and included in this project if they had all of the measured variables in Table 1. We also collated depth-wise measurements of soil characteristics from 167 of the 417 sites plus an additional 110 sites (27 burned, 83 unburned) across 12 wildfires and eight spatially independent unburned areas. Three out of four of the previously mentioned research projects recorded these depthwise soil measurements and additional sites were included from a fifth research project (Supplementary Table S1). Sites were chosen to be representative of burned or unburned forests within  Table S1 for details on sites used for ecoregion models and sites used for soil property models.
each ecoregion using remote sensing imagery and fire history records and stratifying by a combination of drainage conditions or fire severity. See references in Supplementary Table S1 for details on study-specific methods of site selection. Across all research projects, calculations largely followed the methods described in Walker et al., 2018b, but see references in Supplementary Table S1 for details. Briefly, each site was assigned a moisture class based on topography-controlled drainage and adjusted for soil texture and presence of permafrost, ranging from xeric to subhygric, using the method outlined by Johnstone et al. (2008). Stand age at the time of fire was based on tree ring counts from five to ten dominant trees per site. Stem counts, diameter at breast height (DBH) measurements, or study-specific allometric equations were used to calculate tree density (stems m −2 ), basal area (m 2 ha −1 ), aboveground biomass (g m −2 ), and aboveground C content (g C m −2 ). Tree combustion estimates were used to quantify aboveground C combustion These were estimated for each tree as either total percent burned or combustion of structural classes (i.e., foliage, fine branches, large branches, bark). Residual SOL depth was measured at five to 25 points per site. Soil organic layer (SOL) burn depth was estimated based on the height of black spruce adventitious roots above the residual SOL measured on 10 trees per site and calibrated using depth of black spruce adventitious roots in unburned sties (as per Boby et al., 2010;Walker et al., 2018a). In 59 out of the 417 sites, black spruce trees were absent and burn depth was estimated by subtracting the residual SOL from the SOL depth in unburned sites of the same moistureclass (as per Walker et al., 2018a). Pre-fire SOL depth was calculated as the sum of the residual SOL and SOL burn depth. We compiled site-level estimates of residual SOL C, pre-fire SOL C, and SOL C combusted. These were estimated by sampling five to ten soil profiles per site, assessing depth-wise bulk density (g cm −3 ), C content (%), and C stocks (g C m −2 ), and modeling C stock as a function of SOL depth. Based on the above projectspecific estimates, we calculated total C combusted as the sum of above-and belowground C combustion, proportion of prefire C combusted as total C combusted divided by total prefire C, and proportion of total C combusted attributed to the belowground C pool as belowground C combustion divided by total C combusted.
In addition to the site-level variables, we acquired depthwise measurements of soil characteristics from 167 of these sites plus an additional 110 sites (27 burned, 83 unburned) (Supplementary Table S1). The additional 27 burned sites had no data on aboveground C pools and the additional 83 Pre-fire basal area m 2 ha −1 10.6 ± 1.9 ab 9.5 ± 0.8 a 15.4 ± 0.9 ab 8.4 ± 0.5 a 28.5 ± 5. unburned sites had no data on C combustion and thus could not be included in our analyses of ecoregion comparisons. Five to ten SOL profiles were collected per site, and bulk density (g cm −3 ), C content (%), and C stock (g C m −2 ) were assessed from either pre-determined depth increments (e.g., 5 cm) or depth measurements of horizons (e.g., brown moss, fibric, humic). To ensure that all measurements were acquired from the SOL and not mineral soil, we excluded all depthwise SOL samples (n = 3794) with bulk density >0.75 g cm −3 and C content <20%. This resulted in 2685 measurements from 1075 residual or unburned SOL profiles nested within 283 sites (Supplementary Table S1). We adjusted each residual SOL profile depth measurements by the site-level estimate of burn depth. Given that each study grouped SOL depthwise measurements uniquely, we binned the data into 10 cm increments. For example, if the depth of the soil sample or horizon in the pre-fire or unburned SOL profile was <10 cm, we gave it a depth of 10 cm, if the depth was ≥10 cm and <20 cm it was given a depth of 20 cm and so on to a maximum depth of 60 cm.

Statistical Analyses
All statistical analyses were performed using R statistical software version 3.5.2 (R Development Core Team, 2018).
To assess if the structure of these fire-prone ecoregions differed (Q1), we fit linear (LMM) or generalized linear mixed effects models (GLMM) in the package "lme4" (Bates et al., 2015) or the package "glmmTMB" (Brooks et al., 2017). Each model included the fixed effect of ecoregion (six levels) and random effects of fire (18 levels) to allow for different intercepts for each fire scar in order to account for the spatial non-independence of sites within fire scars. Note that because ecoregions cannot be replicated or spatially interspersed, our study is not randomized with replication at the ecoregion level but is randomized with replication within each fire in each ecoregion. We tested response variables of moisture class, pre-fire stand age, density, basal area, proportion of black spruce, aboveground biomass, aboveground C pool, SOL depth, belowground C pool, total C pool, and the proportion of the total C contained in the belowground pool. For stand age, we used a LMM with a normal distribution. For moisture class, we used a GLMM with a Poisson distribution. For each of the proportional response variables we first added a constant of 0.000001 to any zero values and subtracted this constant from any values of one in order to fit a GLMM with a beta distribution and a log-link function. We used GLMM for the remaining variables with a gamma distribution and a loglink function. For these and all mixed models that follow, the significance of fixed effects were determined using maximum likelihood ratio tests comparing the full model to a reduced model and confirmed using Akaike information criterion (AIC, AIC < 2.0; Zuur et al., 2009). Because we had multiple response variables with the same experimental units we corrected the p-values obtained from likelihood ratio tests comparing the full model to a reduced model using the false discovery rates (FDR; Benjamini and Hochberg, 1995). When a significant effect of ecoregion was detected based on the adjusted p-values (p-value <0.05) we used Tukey-Kramer post hoc analysis for multiple comparisons in the R package "emmeans" (Lenth et al., 2019) to test for pairwise differences in marginal means between ecoregions. For these and all LMM or GLMM that follow we assessed bias in model fit by visually inspecting residual versus fitted values.
To assess if aboveground, belowground, and total C pools increased over time at a similar rate between ecoregions (Q2), we fit LMM with response variables of pre-fire aboveground C, belowground C, and total C, all of which were natural log transformed to ensure normality. For each model, we included fixed effects of stand age at the time of fire (i.e., time after establishment), ecoregion (six levels), and their first order interaction and a random effect of fire (18 levels). The importance of fixed effects and p-value adjustments proceeded as described above. When a significant interaction with ecoregion was found, we used a Tukey-Kramer post hoc analysis to complete pairwise comparisons of marginal mean intercepts and slopes between ecoregions.
To determine how (i) bulk density (g cm −3 ), (ii) C content (%), and (iii) C stocks (g C m −2 ) change with depth in the SOL profile and if these vary with total pre-fire or unburned SOL thickness and ecoregion (Q3), we fit LMM for each soil characteristic. Bulk density and C stocks were natural log transformed for normality. We grouped ecoregions into four categories to ensure sufficient sample size within each unique combination of ecoregion and total pre-fire or unburned SOL thickness category. Taiga Plains (n = 290 soil profiles) and Taiga Shield (n = 192 soil profiles) were left as is, but Alaska Boreal Interior and Boreal Cordillera were grouped as "Alaska" (n = 311 profiles) and the Boreal Plains and Softwood Shield were grouped as "Saskatchewan" (n = 282 soil profiles). We assigned each site to a SOL thickness category (total SOL; 0-20 cm, 20-40 cm, and >40 cm) and modeled soil characteristics using fixed effects of depth within soil profile (depth), ecoregion, total pre-fire or unburned SOL thickness, and the first-order interactions between depth and ecoregion and depth and SOL thickness. For these LMM, we used random effects of soil profile (1075 levels), nested within site (283 levels), nested within fire or unburned area (20 levels). The importance of fixed effects and p-value adjustments proceeded as above. Differences in marginal mean intercepts and slopes of bulk density, C content and C stocks as a function of depth were compared between ecoregions within each SOL thickness class and between SOL thickness classes within each ecoregion using Tukey-Kramer post hoc analysis.
To assess if metrics of fire severity and C combustion varied among ecoregions (Q4), we fit LMM or GLMM, with response variables of residual SOL depth, residual SOL C, burn depth, aboveground, belowground, and total C pools, the proportion of total pre-fire C combusted, and the proportional contribution of the belowground C pool to total C combusted. For each model, we included random effects of fire scar (18 levels) and a fixed effect of ecoregion (six levels). For burn depth, we used a LMM with a normal distribution. For each of the proportional response variables, we fit a GLMM with a beta distribution and a loglink function. We used GLMM for the remaining variables with a gamma distribution and a log-link function. The importance of fixed effects and p-value adjustments proceeded as described  above. When a significant effect of ecoregion was detected, we used Tukey-Kramer post hoc analysis to test for pairwise differences in marginal means between ecoregions.

RESULTS
In assessing the pre-fire structure of these fire-prone ecosystems (Q1), we found that stand age, proportion of black spruce, stand basal area, aboveground biomass, aboveground C pools, SOL depths, belowground C pools, total C pools, and the proportion of total C stored belowground differed between ecoregions (Figure 1, Table 1 and Supplementary Table S2). Only site drainage class and pre-fire stand density did not differ between ecoregions (Table 1). In general, the northern Taiga and Alaskan ecoregions were more similar to one another than to the southern boreal ecoregions (Boreal Plains and Softwood Shield). Compared to the southern ecoregions, sites in the northern boreal ecoregions burned at an older age (Figure 2A), were more dominated by black spruce (Figure 2B), had lower aboveground C pools (Figure 2C), and stored proportionally more of the total ecosystem C belowground than aboveground ( Table 1). In addition to these north versus south differences, we also observed distinct characteristics of the Taiga ecoregions, where pre-fire SOL depth, belowground C pools (Figure 2E), and total prefire C pools were higher compared to the Alaskan or southern boreal ecoregions.
Ecoregion differences were also apparent in the rates at which C accumulated (Q2). We found significant interactions between stand age and ecoregion in predicting each of prefire aboveground, belowground, and total C (Figures 3A-C, Table 2 and Supplementary Table S3). Pre-fire aboveground C increased slightly in the Taiga Plains and Taiga Shield, but we found no significant relationships between aboveground C pools and time of stand establishment in the Alaskan (Alaska Boreal Interior and Boreal Cordillera) or southern boreal (Boreal Plains and Softwood Shield) ecoregions ( Figure 3A and Table 2). Pre-fire belowground C pool increased over time in the Boreal Plains but did not change with time after fire in the other ecoregions (Figures 3B,C and Table 2). Total C pools, the majority of which is from belowground (Table 1 and Figure 3C), significantly increased over time in the Taiga Shield and Boreal Plains ( Figure 3A and Table 2). The rate of increase in each of aboveground, belowground, and total C pools was much higher, albeit not significantly different due to large standard errors, in the southern boreal ecoregions compared to the northern. FIGURE 3 | Ecoregion specific effects of time after the previous fire, based on stand age at the time of the most recent fire, on natural log transformed pre-fire (A) aboveground C pools, (B) belowground C pools, and (C) total C pools (g C m −2 ). Colors represent ecoregions and lines represent model fits of ecoregion specific significant relationships and 95% confidence intervals. Non-significant relationships are not shown. See Supplementary Table S3 for model output and Table 2 for estimated intercepts and slopes.
In examining the relationships between soil characteristics and depth (Q3), we found that bulk density (g cm −3 ), C fraction (%), and C stock (g C m −2 ) varied with depth and that these relationships differed between ecoregions and total thickness of the pre-fire or unburned SOL profile (Figure 4, Table 3 and Supplementary Table S4). We found that bulk density (g cm −3 ) always increased with depth and within each total SOL thickness categories this increase was always greatest for the Taiga Plains. When comparing within each ecoregion, increases with depth were greatest when total SOL thickness was 0-20 cm (Figure 4 and Table 3). Carbon content (%) decreased with depth in all ecoregions and SOL thickness category, but this decrease was greatest for all ecoregions when SOL thickness was <20 cm. Within each total SOL thickness category the decrease in C content with depth was greatest in the Taiga Plains. In all cases, C stocks (g C m −2 ) increased with depth and this rate of increase was greatest for shallow SOL (<20 cm) in the Taiga Plains (Figure 4 and Table 3). Despite having high C stocks, the Taiga Shield had the smallest increase in cumulative C stocks with depth across all SOL thickness categories (Figure 4 and Table 3).
We found ecoregion differences in most of our metrics of fire severity (Q4). Specifically, burn depth, residual SOL depth and C, aboveground C combustion, proportion of C combustion from belowground, and proportion of total prefire C combusted were different between ecoregions (Figure 2, Table 1 and Supplementary Table S2). Burn depth (Figure 2E) was slightly higher in the Alaskan ecoregions and residual belowground C pools ( Figure 2F) were highest in the Taiga ecoregions, but these differences did not result in any significant ecoregion differences in belowground or total C combustion ( Table 1). The most notable difference in fire severity between ecoregions, was that aboveground C combustion was higher ( Figure 2F) and the proportion of total C combustion coming from belowground was lower ( Figure 2G) in the southern compared to the northern ecoregions.

DISCUSSION
Through compiling a spatially extensive dataset spanning six ecoregions in North America's western boreal forests, we were able to assess how the structural and functional attributes that relate to above-and belowground C stocks of these fire-prone ecoregions differ. Our study sites captured a broad gradient in pre-fire conditions of tree productivity, stand age, and ecosystem C storage as well as burn depth and C combustion from fire both within and among ecoregions. Given the breadth of our study sites, and that they were chosen to be representative of burned forests within each ecoregion, our results should allow future research to estimate C accumulation, storage, and combustion throughout the western North American boreal forest.
Our regional analysis highlighted pronounced ecoregion differences associated with latitude in boreal forest structural and functional attributes relating to C sequestration and storage. In general, southern ecoregions (Boreal Plains and Softwood Shield) were more similar to one another than to the northern ecoregions, suggesting that differences in longterm climate between these regions are more important than ecoregion differences in geological and biogeographical history, soil development, and parent materials, for predicting ecosystem C dynamics. The age distributions of our examined field sites TABLE 2 | Marginal mean intercept and slope estimates and standard errors of pre-fire aboveground carbon (C), belowground C, and total C as a function of stand age for each ecoregion. suggest that there is a higher probability of younger stands burning in southern compared to northern boreal ecoregions. Northern ecoregions also had smaller aboveground C pools, larger belowground C pools, and lower rates of C accumulation. Taken together, these differences suggest that sufficient biomass for fire ignition and spread occurs at an earlier stage in stand development in the southern compared to the northern boreal ecoregions. However, we did not account for potential differences in lightning rates between ecoregions, which is correlated to temperature and precipitation and a dominant driver of total burned area in boreal North America (Veraverbeke et al., 2017). Although C accumulates slowly in northern boreal forests, the relative low severity and frequency of past fires has allowed soil C pools to accumulate and act as a reservoir for long-term C storage (Harden et al., 2000;Chapin et al., 2006). This has resulted in deep SOL and large belowground C pools in the Alaskan and Taiga ecoregions. These soil C sinks are particularly high in the Taiga ecoregions where peatlands, bogs, and fens are common (Tarnocai et al., 2009;Hossain et al., 2015). However, the soil C pools in northern ecoregions may be more likely to combust with a changing fire regime  or changing land use (Turetsky et al., 2011a) than they have in the past. The higher frequency of fires in southern boreal regions, in combination with storing proportionally more C above-than belowground and having lower belowground C pools, suggests that the longterm C sink of southern boreal forests is weak compared to the northern ecoregions.
Warmer temperatures increase rates of organic soil decomposition, which release nutrients and increase aboveground C production and total net primary productivity (Chapin et al., 2002;Day et al., 2019). This leads to increased litter fall, root turnover, and an overall higher yield in C input per year in boreal forest soils (Deluca and Boisvenue, 2012). As the southern ecoregions have higher mean annual temperatures and a greater number of growing degree days (>5 • C) than the northern ecoregions, these long-term climate parameters are likely important determinants of the higher aboveground C storage and accumulation rates we observed in the southern ecoregions. Similar shifts in C distribution from below-to aboveground pools as temperatures increase have been observed in black spruce forests of Interior Alaska (Kane and Vogel, 2009).
Whether stands established post-fire or post-timber harvest also influence boreal C dynamics (Seedre et al., 2014;Dieleman et al., 2020). Boreal forest stands that originate following harvest tend to have higher rates of aboveground C recovery than those that establish after fire due to reforestation practices or because harvesting activities usually occur on more productive sites (Ilisson and Chen, 2009;Seedre et al., 2014). In the southern boreal ecoregions, almost half of the sites (19 of 44) established following harvest rather than fire, which might also explain why these ecoregions had higher aboveground C pools compared to the northern ecoregions .
Differences in species composition also influence C storage and accumulation rates; black spruce are generally less productive and slower growing than jack pine or deciduous trees (Subedi and Sharma, 2013;Alexander and Mack, 2016). Black spruce was dominant in the majority of sites in northern ecoregions, but in the southern boreal ecoregions, jack pine dominated in 18 out of 44 sites. As the frequency and severity of boreal wildfires continue increasing in these northern ecoregions, stands previously dominated by black spruce are likely to transition to jack pine or deciduous dominance Whitman et al., 2018;Hart et al., 2019). This change in species composition in combination with direct impacts of climate warming on C storage and accumulation rates, suggests that the future C dynamics of northern boreal forests will resemble those of southern boreal forests.
As expected, soil bulk density and C stocks increased and C content decreased with depth in the SOL profile and the rate at which these soil characteristics changed with depth were dependent on ecoregion and total pre-fire or unburned SOL thickness. The differences in these rates were generally greater between total SOL thickness categories than between ecoregions, suggesting that soil development is more closely related to local drainage conditions than to ecoregion differences. Interestingly, within each total SOL thickness category the Taiga Plains ecoregion had the highest bulk density and C content, potentially due to extensive lacustrine deposits from the last major glaciation that are found throughout this ecoregion (Ecosystem Classification Group et al., 2009). These results highlight the importance of using regionally specific estimates of SOL C storage for predicting C combustion. FIGURE 4 | Soil properties of bulk density (g cm −3 ), carbon content (%) and carbon pool (g C m −2 ) as a function of sample depth, total plot SOL thickness, and ecoregion. Lines represent significant ecoregion specific relationships between soil properties and depth in the soil profile and shading represents 95% confidence intervals. Note that points are jittered horizontally at each 10 cm depth increment (width = 2 cm). See Supplementary Table S4 for model results and Table 3 for Tukey-Kramer post-hoc pairwise comparisons of marginal mean intercepts and slopes compared across total SOL thickness categories within ecoregions and across ecoregions within total SOL thickness categories.
Fire severity, measured as the proportion of total ecosystem C combusted, was close to 50% in Alaska and approximately 30% in the other ecoregions (Table 1). Lower proportional combustion in the southern ecoregions is likely due to more C stored aboveground in tree stems that is unavailable for combustion. In contrast, we attribute the lower proportional 3 | Linear mixed effect model estimated marginal mean intercepts and slopes (± standard error) for the effect of sample depth (depth) on the response variables of bulk density, carbon, and carbon pool for each combination of ecoregion (four levels) and total soil organic layer thickness (three levels). Note that Bulk density and Carbon Stock were natural log transformed for the model and the back transformed values are presented here. Bold values indicate nonsignificant effects (p-value > 0.05) and non-bolded values are significant (p-value < 0.05). Superscript letters represent significant differences between ecoregions within each total soil organic layer thickness category, and superscript numbers represent significant differences between total soil organic layer thickness categories within each ecoregion. The models were fit on 2596 soil increments from 1041 soil profiles in 277 plots in 20 burn scars/unburned areas sampled in four different projects. See Supplementary Table S4 for original model results and Figure 4 for a graphical depiction.
combustion in the Taiga Plains and the Taiga Shield to deep SOL and large belowground C pools (Table 1) in poorly drained areas where deep soils are protected from burning by permafrost or saturated water table conditions (Turetsky et al., 2011b). Burn depth, residual SOL depth and C, and the proportion of total C combustion from belowground were lowest in the southern ecoregions, emphasizing the high C combustion from more rapidly accumulating aboveground biomass and shallow pre-fire SOL depths in these ecoregions. The fact that burn depth differed but belowground C combustion did not highlights the importance of local or regional conditions in modeling C combustion. As climate continues to warm and boreal wildfire regimes continue to intensify, the structure and function of the boreal forest is likely to change. Notably, the boreal net ecosystem C balance could switch from net C sink to a net C source . Understanding C combustion from fires is critical for forecasting such changes and the results we present in this study should allow future research to estimate C accumulation, storage, and combustion throughout the western North American boreal forest. Our findings that C storage, accumulation rates, and sources of C combustion from fire vary latitudinally suggests that the C dynamics of northern boreal forests are likely to shift and resemble the dynamics of southern boreal forests in the future. Specifically, as climate continues to warm northern regions are likely to start accumulating C more rapidly, store more C aboveground, have smaller belowground C stocks, burn more frequently, and emit proportionally more C from the aboveground component than from belowground. Understanding these changes to the structure and function of fire-prone boreal ecoregions are needed to quantify the role of fire in the global C cycle and its feedbacks to climate change.

DATA AVAILABILITY STATEMENT
The data used in this manuscript is archived at the Oak Ridge National Laboratory Distributed Active Archive Center (ORNL DAAC), doi: 10.3334/ORNLDAAC/1744 (Walker et al., 2020).

AUTHOR CONTRIBUTIONS
XW and MM conceived the study. XW compiled datasets, analyzed the data, and led the writing in collaboration with MM. All authors established the sites and contributed to field data, read, edited, and contributed to the writing of this manuscript.

FUNDING
This writing of this manuscript and synthesis of data was supported by funding the NASA Arctic Boreal and Vulnerability Experiment (ABoVE) Legacy Carbon grant NNX15AT71A awarded to MM. The original field studies were supported by funding in the United States from NSF DEB RAPID grant #1542150 to MM, NASA Rapid Response grant NNX15AD58G and NASA ABoVE grant NNX15AT83A to LB-C, NASA ABoVE grant NNX15AU56A to BR, SV, and MT, Joint Fire Science Program grant 05-1-2-06 to JJ, NSF grant 0445458 to MM, NSF support to the Bonanza Creek LTER (DEB-0423442); and in Canada from NSERC Discovery Grant funding to JJ and MT; Government of the Northwest Territories Cumulative Impacts Monitoring Program Funding project #170 to JB; NSERC PDFs to ND and CD; GNWT logistical and financial support through the Laurier-GNWT Partnership Agreement; Polar Knowledge Canada's Northern Science Training Program funding awarded to Canadian field assistants; SV acknowledges Vidi grant support from the Netherlands Organisation for Scientific Research (NWO).