Using Holocene paleo-fire records to estimate carbon stock vulnerabilities in Hudson Bay Lowlands peatlands

The study site is located on the western margin of the HBL, Canada (Fig. 1; the Hudson Plains Ecozone of the Ecological Stratification Working Group 1995). The study site and surrounding area are underlain by Paleozoic carbonates, shales, and siltstones and are close to the boundary with the Archean Superior Province of the Canadian Shield (Ontario Geological Survey 2011). Surficial deposits consist of peat and massive to bedded diamicton (Barnett et al. 2013a; 2013b). Land coverage is almost entirely peatlands and other wetlands, with minor extent of coniferous forest in slight upland regions (Fig. 1; Government of Ontario 2021). Peat depth is generally shallow in the area immediately surrounding the study site, ranging from 100 to 150 cm.

The HBL has a humid microthermal Arctic climate (Martini 2006). Mean annual temperature is −1.3 °C and total annual precipitation is 700 mm at the closest weather station (Lansdowne House; approximately 120 km from the study site; Government of Canada 2021). The highest and lowest average monthly temperatures occur in July and January, respectively (17.2 °C and −22.3 °C). Snowfall >20 cm occurs in 7 months of the year (October to April; 242 cm yr⁻¹ total; Government of Canada 2021).

The study site was accessed by helicopter and is located within the northwestern portion of an area that burned in the summer of 1991 (Fig. 1; identified via satellite imagery [Landsat]; Canadian Forest Service 2020; Hall et al. 2020). The fire was caused by lightning and an estimated 76–100% of the delineated area in Fig. 1 was burned (420 Ha; Canadian Forest Service 2020; Hall et al. 2020). The burned portion of the study site is a shrub-rich bog with hummock and hollow microtopography and is predominantly surrounded by unburned forested bog (Fig. 1). Vegetation consisted of predominantly Picea mariana, Rhododendron groenlandicum, and Sphagnum and Cladina spp. in both the burned and unburned portions of the study site. Chamaedaphne calyculata, Maianthemum trifolium, and Rubus chamaemorus were also present. Eriophorum vaginatum and Larix laricina were also found infrequently surrounding each of the coring localities. Water table depth and pH were measured at two points at the site at the time of sampling (hollow and high hummock) and were 42 and 66 cm and 5.0 and 3.6, respectively. Three cores were collected in August 2018 at the study site, “ROF-F01,” “ROF-F02,” and “ROF-F03” (184 m asl; Fig. 1; Table 1). The cores were taken in different regions of the site (burned-edge, burned-center, unburned) and in different microtopography (low hummock, hollow, high hummock) to capture some of the spatial variability (Table 1). Coring was done using a 1-m Jeglum corer for the first 100 cm; subsequent drives were collected with a 1 m × 15 cm diameter Russian peat corer. Basal sediments underneath the peat at each coring locality consisted of a clay matrix and grains of silt and sand. Cores were stored at 4 °C prior to analysis at the University of Toronto, Canada.

The approach for sampling for chronology and peat properties was the same for each coring locality, with exceptions noted below. Ages were assigned to each core using Bayesian age-depth modeling in the “rbacon” package for R (Blaauw and Christen 2011). Each model was developed using radiocarbon dates on aboveground plant material (Table 2) and the surface age (−68 cal yr BP). A total of seven to nine radiocarbon dates per core were measured to obtain a resolution of roughly one date per 1000 years. Cores were initially dated at 45 cm, 95 cm, and the base of the core, and additional dates between these were measured to obtain the targeted resolution. Conventional ¹⁴C ages were calibrated to “calendar” (cal) years using the IntCal20 calibration curve (Reimer et al. 2020). Radiocarbon ages were measured at the Lalonde AMS Laboratory (Ottawa, Canada). For the coring locality adjacent to the burned area (ROF-F02), where there was no suspected disturbance of the upper peat column from the most recent known fire, Pb-210 activities from bulk peat samples in the upper 50 cm were also used in the Bayesian age-depth model for the location using the “plum” package for R (Table 3 and Fig. S1; Blaauw et al. 2021). Pb-210 samples were measured at the Flett Research Laboratory using alpha spectrometry (Winnipeg, Canada). Select samples across the upper 50 cm of the core were measured to obtain sufficient data to determine exponential decrease in Pb-210 activity to background levels (Table 3 and Fig. S1). Because the coring locality had Pb-210 ages, radiocarbon ages above 45 cm were not measured. We placed a hiatus in a record when the average accumulation rate between radiocarbon ages was much lower than our prior accumulation rate in our Bayesian age-depth model (i.e., << 50 years·cm⁻¹; threshold of 5 × lower = 250 years cm⁻¹).

To estimate long-term and apparent carbon accumulation rates (aCARs), carbon density, and total carbon mass, and to investigate peat properties (carbon and nitrogen content), dry bulk density measurements were performed at a 1 cm resolution to match the highest resolution of the charcoal record. A total of 2–28 cm³ was taken from each 1 cm interval to obtain >0.5 g of the dried sample for carbon and nitrogen analysis. Each sample was dried to a constant mass overnight at 60 °C and then ground to a fine powder using a MM200 Retch ball mill. Carbon and nitrogen analysis was performed at the Radiochronology Lab at the Centre for Northern Studies, Université Laval (Québec City, Canada), where 70–100 mg per sample was analyzed using a Leco CHN628 analyzer. Standards analyzed alongside samples were within 0.6% relative standard deviation for both carbon and nitrogen (EDTA Lot1061 and 502-309 soil standard). Carbon accumulation rates were calculated using two methods, aCARs and long-term rates of carbon accumulation (LORCAs). LORCAs represent the net carbon storage at the site through its averaged accumulation history, while aCARs are time-stepped estimates approximating net carbon uptake at a given time (Young et al. 2021). aCARs were calculated by multiplying the dry bulk density (in g cm⁻³) by the accumulation rate (cm yr⁻¹) and the proportion of carbon in each sample (Chambers et al. 2010). Cumulative carbon mass was determined by summing the carbon densities (in kg C m⁻³) multiplied by the thickness (in m) for each sample. The cumulative carbon mass (in kg C m⁻²) was then used to calculate LORCA for each coring locality by dividing it by the basal age.

To assess fire frequency through the Holocene, subsamples of 2 cm³ of peat were taken contiguously along each core at a 1–2 cm interval. Each sample was placed in a 5% KOH treatment overnight (>16 h) and then sieved to >150 μm (after Mooney and Tinner 2010; van Bellen et al. 2012). The entire sample was counted using a stereomicroscope with a gridded Petri dish at 32 × magnification. Charcoal fragments were identified by referencing New et al. (2016) (Fig. S2). Peat type was assigned to each sample based on the main vegetation groups observed. Samples were assigned to one of three groups: Sphagnum-dominated peat, mixed ligneous and Sphagnum peat, and mixed ligneous, Sphagnum, and herbaceous peat. To be assigned a group, each of the various plant macrofossil components was at least 25% of the plant macrofossils present in the sample.

Charcoal counts were analyzed to identify fire peaks using “CharAnalysis” (Higuera 2009). Each core was rescaled to the median sample resolution prior to analysis to reduce bias from a changing vertical accretion rate in the site’s accumulation history (Table 1). If a core had a suspected hiatus, it was split into two components (i.e., above and below the boundary; Section A and B; Table 1). The rescaled charcoal accumulation rate (CHAR) was split into three components: background (C_back), noise (C_noise), and fire (C_fire). The background component is attributed to changes in long-term regional charcoal production and transport, post-fire deposition, and the potential movement of charcoal within the acrotelm (Higuera et al. 2009; van Bellen et al. 2012). Both the noise and fire components are part of the peaks (C_peak) identified above C_back, where C_noise is attributed to random variability due to external inputs, and the C_fire represents charcoal produced due to local burning at or near the site (Gavin et al. 2006; Higuera et al. 2010). We defined C_peak as a residual, assuming that there is an additive relationship. We identified C_back using LOWESS smoothing (robust to outliers) and defined the smoothing window to have >30 samples per window (Higuera et al. 2010). We used a local threshold to separate C_fire and C_noise to account for variation in the record with decomposition when the core intervals analyzed were greater than the smoothing window (Table 1; Higuera et al. 2010). C_fire was considered to represent a local peat fire when it was above the 99% threshold of the C_noise distribution identified using Gaussian mixture models (Gavin et al. 2006). We utilized a signal-to-noise index (SNI) to evaluate each of the peaks in our record (after Kelly et al. 2011). A SNI greater than 3 for a given peak was the threshold to consider a peak well separated from C_noise. To separate peaks that directly impacted the peat column and those that impacted vegetation inputs, we used C:N ratio, as N is enriched relative to C during peat smoldering (Zaccone et al. 2014).

Relationships between FRI and aCAR through the Holocene were investigated using the combined datasets of the three coring localities. FRI and aCAR values were assigned to one of the following intervals: 0–500 cal yr BP, 500–4200 cal yr BP, and 4200–7200 cal yr BP. We selected 500 cal yr BP as a boundary because aCARs shift to much higher values associated with acrotelm peat that cannot be directly compared to deeper values (Young et al. 2019). Further, there could be potentially higher amounts of charcoal movement through fluctuating water tables in acrotelm peat. Therefore, these values were omitted from the statistical analysis of the difference between aCAR values through time and their relation to changing FRIs. The boundary between the other two intervals is the stage boundary between the Middle and Late Holocene (4200 cal yr BP; Walker et al. 2019). For FRIs that overlapped the boundaries, they were assigned to the interval where the majority of the FRI occurred. A Mann–Whitney U test was performed to test differences between aCARs between the Middle (4200–7400 cal yr BP; N = 74, 3 coring locations) and Late Holocene (500–4200 cal yr BP N = 123, 3 coring locations). This test was selected because the aCAR data are non-normal (Shapiro–Wilk, p < 0.05). Statistical tests were not performed on the FRI data, as sample size for each interval was low, even in the combined record. Median fire frequency values were used in the comparison between time periods to minimize the influence of outliers. We estimated carbon release for the Middle and Late Holocene time bins (in kg C m⁻²) by multiplying the median fire frequency (1/FRI; in fires·yr⁻¹) and time bin interval (yr) by the annual average emission per unit of burned area within the Hudson Plains (1430 g C m⁻²·fire⁻¹; Balshi et al. 2009). This annual average emission value was selected as it is a modeled estimate for the entire region and also lies conservatively within estimates of carbon loss from fire in boreal peatlands, which is expected in more humid regional climates such as the HBL where fire is more likely restricted to aboveground biomass (0–6000 g C m⁻²·fire⁻¹; Zoltai et al. 1998; Turetsky and Wieder 2001; Lukenbach et al. 2015b; Wilkinson et al. 2018).

Holocene-scale fire frequencies were used to place our study region in the wider context of the relationship between fire frequency and LORCAs in peatlands across the boreal region. A Holocene-scale fire frequency (in fires 1000·yr⁻¹) was calculated for each coring location by dividing the total number of fire peaks by the basal age. A Holocene FRI was then calculated as the inverse of fire frequency. We estimated carbon release at each coring location for the entire record (in kg C m⁻²) by multiplying the median fire frequency (1/FRI; in fires·yr⁻¹) and the core basal age (yr) by the annual average emission per unit of burned area within the Hudson Plains as above (1430 g C m⁻²; Balshi et al. 2009). To compare our study site to other boreal regions, we compiled fire frequency and LORCA data from other boreal peatlands (Kuhry 1994; Turunen et al. 2001; van Bellen et al. 2012; Magnan et al. 2020) and investigated their relationship using linear regression analysis. A linear relationship was explored as this was previously found to be significant (Kuhry 1994; Moore and Robinson 2000). Fire frequency was square-root transformed prior to linear regression analysis to normalize the dataset (Shapiro–Wilk, p > 0.05).

We utilized the “Global Charcoal Database” to create compilations of eastern and western boreal fire histories to compare the site investigated in this study. Sites within the boreal region and with radiocarbon ages (Fig. 1; boreal region defined after Levavasseur et al. 2012) were extracted and analyzed using the “paleofire” package in R (Blarquez et al. 2014). The boundary of 82 °W was used to separate western and eastern boreal sites, based on paleoclimate data from sites in western Ontario (Lake 239; Fig. 1; Moos and Cumming 2012). Each record was transformed prior to compilation by scaling the records between 0 and 1, Box-Cox normalization, and converting to a z-score with the base period set to the last 8000 years (Power et al. 2010). All sites were lake sediment records, which provide a more regional signal than peat charcoal records, and all had influx data. The composite curve for each region was created using the “pfComposite” function in the “paleofire” package (Blarquez et al. 2014). The composite was the mean z-score in 500 year bins, with bootstrap resampling (1000 iterations). We selected 500 year bins to match the lowest sample resolution in the included studies (# of samples divided by the length of record; see Tables S1 and S2).

Peat accumulation initiated at all three coring localities by 5600 cal yr BP or earlier (Table 1; Fig. 2). All three coring locations are interpreted to have rapidly transitioned to either treed poor fen or bog conditions after initiation, with the presence of significant Sphagnum moss and ligneous (woody) plant remains throughout the cores, including the consistent occurrence of conifer needles (Fig. 3, Vitt 2006; Loisel et al. 2014). Rapid succession (i.e., from marsh and/or rich fen conditions or from forest to wetland through paludification) can occur in relatively drier localities controlled by local topography and drainage (Klinger and Short 1996; Glaser et al. 2004). Further, paludification can be accelerated by fire, and since there are high charcoal concentrations at the base of two cores, fire at the site in this study may have also promoted terrestrial to wetland succession variably across the study region and explain the large spread in basal ages (Simard et al. 2007; Magnan et al. 2018).

At all three coring localities, peat properties (BD and C:N ratio) generally vary with either peat type and (or) charcoal concentration, while aCARs are predominantly related to shifts in vertical accretion rate in the age-depth model (Figs. 2 and 3). Further, shifts in broad peat types do not correspond to changes in charcoal concentrations, suggesting that the Holocene fire records for the study region are reflective of changes in regional climatic conditions and not vegetation change (Fig. 3). At ROF-F01 (edge of burned area), shifts in C:N ratio broadly correspond to shifts in peat type. A hiatus was placed in the record at 51 cm, between the two radiocarbon ages where the deposition time was greater than 400 yr·cm⁻¹ (Fig. 2). The highest nitrogen content is above the hiatus in the record, reaching 1.4–1.6 % between 180 and 400 cal yr BP (Fig. 3). High nitrogen content and the corresponding low C:N ratio in relation to the hiatus indicate peat smoldering (Zaccone et al. 2014) and support the placement of the hiatus. At ROF-F02 (adjacent to burned area), C:N ratios are more variable, along with peat type (Fig. 3). The highest nitrogen values at the coring location correspond to the highest charcoal concentration in the record, also suggesting that N was enriched with peat smoldering (Fig. 3; Zaccone et al. 2014). At ROF-F03 (center of burned area), the lower portion of the record (>4000 cal yr BP) has higher peat bulk density that corresponds to lower C and N content and relatively high charcoal concentration (Fig. 3). Unlike the other two coring locations, herbaceous material as a major component of the peat type is restricted to lower portion of the core (>5800 cal yr BP). Further, Sphagnum peat predominates in the uppermost portion (above 65 cm or 1500 cal yr BP). High nitrogen content occurs between 2300 and 3300 cal yr BP, reaching values of up to 1.9% (Fig. 3). A fire peak lies within this interval, but charcoal concentration is low (Fig. 3). All cores have high aCARs associated with high vertical accretion rates in the uppermost acrotelm peat (Figs. 2 and 3).

Differences between LORCA, carbon density, and total carbon mass at the coring locations correspond to differences in basal age and depth, microtopography, and vegetation cover (Table 1). Mean carbon density ranged from 47 to 57 kg C m⁻³ for the three coring localities, with higher carbon densities in the shrub-rich portions of the burned area with older basal ages and that were collected from hummocks (i.e., ROF-F01 and ROF-F03). Total carbon mass of the peat profiles was the highest at the oldest coring locality (ROF-F03) and lowest at the site with the shallowest peat depth (ROF-F01). LORCAs for the sites were between 8 and 11 g C m⁻²·yr⁻¹, with the highest occurring in the treed core locality in a hollow (ROF-F02; Table 1).

The rescaled CHARs (median resolution) range from 0 to 8 pieces cm⁻²·yr⁻¹ for all coring locations (Fig. 4). A total of five to nine fire peaks per core were identified. All peaks across the three cores are well separated from C_noise, since the median SNI was >3 globally, locally for more than 85% of the samples in each core, and for each significant peak (Table 1; Kelly et al. 2011). One fire peak was also placed at the top of the hiatus of ROF-F01 independent of the peak determination in CharAnalysis, as we interpret this hiatus to be caused by peat combustion. The upper boundary of the hiatus also corresponds to fire peaks in the two other records when accounting for age error (i.e., within the 95% CI; Fig. 4). At each coring locality, 3–4 peaks were associated with a shift in C:N ratio that was lower than the core mean, suggesting peat smoldering occurred (Table 1 and Fig. 3). Further, increased CHARs are recorded in the uppermost portion of the record at every coring locality, but fire peaks are only significant for the cores within the burned area (i.e., ROF-F01 and ROF-F03; Figs. 1 and 4). Only the age of the uppermost significant peak at ROF-F03 matches the year burned (1991). However, the cores within the burned area do not have Pb-210 ages to constrain very recent vertical accretion rates, as we did not perform Pb-210 dating where disturbance from the most recent fire occurred (Fig. 2). The age range surrounding the increased CHARs that were not significant at the unburned coring locality (ROF-F02) also contains the year burned (1991) and is constrained by Pb-210 ages (Fig. 4). Holocene-scale FRIs (0–8000 cal yr BP) and their corresponding fire frequencies were similar for the three cores, ranging from 907 to 990 years and 1.0 to 1.1 fires 1000·yr⁻¹, respectively (Table 1). As a result of fire at the site, we estimate that 8–12 kg C m⁻² was lost to the atmosphere over the site’s accumulation history (Table 1).

The fire records at the coring localities generally record similar timing in fire peaks, with at least 70% of the fire peaks overlapping in age with at least one other peak (Fig. 4). However, some peaks are independent despite the localities being < 100 m apart from each other, even when age error is considered (Figs. 1 and 4). Local-scale variability in charcoal records has been demonstrated on replicate and site-level records in other peatlands, as well as at various charcoal size fractions, suggesting that peatland charcoal records are highly localized to the core location (e.g., van Bellen et al. 2012; Cui et al. 2020; Beaulne et al. 2021). However, the FRIs through the Holocene for each of the core localities fit within the upper end of modern estimates of FRIs for the HBL (200–1000 years; Stocks et al. 2002; Boulanger et al. 2014; Coops et al. 2018). Therefore, the paleo-fire return intervals recorded at the site match regional expectations.

Overall, more fires occurred in the Middle Holocene in comparison to the Late Holocene. When individual FRIs from across the three coring locations were pooled to investigate Middle (4200–7200 cal yr BP) and Late Holocene (500–4200 cal yr BP) differences, the median FRI for the Middle Holocene (912 years) was lower than the Late Holocene (1650 years), corresponding to a higher median fire frequency (1.1 vs. 0.6 fires 1000·yr⁻¹; Figs. 5 and 6) and subsequent carbon loss (6 kg C m⁻² vs. 3 kg C m⁻²) across the interval. Fire frequency patterns at the site match the trend in the western boreal compilation, where there were higher CHARs within lake records during the Middle Holocene, signaling more fire on the landscape during that time (Fig. 6).

The relationship between changing fire frequencies and peatland carbon accumulation differed when investigating at a local vs. regional scale. aCARs were not significantly different between the Middle and Late Holocene at the study site in the combined core dataset (W = 4352, p = 0.609; Fig. 5), despite differing median fire frequency between the time bins. However, when LORCA and fire frequency from this study are plotted with other boreal sites that represent a larger dataset across a larger fire frequency and moisture gradient, there is a significant relationship, with decreasing Holocene-scale carbon accumulation rates with increasing fire frequency (Fig. 7).

Median fire frequency recorded from the study region is related to differing hydroclimatic conditions between the Middle and Late Holocene. Higher median fire frequency in the Middle Holocene corresponds to warmer conditions with minimal precipitation change in the HBL (<0.5 °C; Fig. 6; pollen-inferred estimates; Hargan et al. 2020), which has been linked to relatively drier conditions in testate amoeba records because of higher evaporative demand (Fig. 6; Davies et al. 2021). Further, the pattern at the study site is also comparable to Lake 239 record in northwestern Ontario, where despite higher precipitation within the Middle Holocene, higher evaporative demand and likely more instances of prolonged summer drought led to higher fire frequency surrounding the site (12.5 fires 1000·yr⁻¹ in the Middle Holocene and <5 fires 1000 yr⁻¹ in the Late Holocene; Moos and Cumming 2011; 2012; Fig. 6). Because a difference in median fire frequency is recorded between these time periods at the study region, fire frequency is apparently impacted by relatively small changes in temperature and evapotranspiration budgets and therefore supports that some HBL peatlands may be increasingly vulnerable to fire because of anthropogenic climate change.

Different fire frequency patterns between the site and the eastern boreal compilation support that regional factors are causing alternate fire regimes. In the Quebec region, the Middle Holocene was characterized by higher relative humidity associated with less frequent intrusion of cool and dry Pacific air masses (Edwards et al. 1996). Less frequent intrusions meant there were fewer periods of summer drought that could promote fire (Carcaillet and Richard 2000). Since the western HBL has a higher median fire frequency in the Middle Holocene, cool and dry Pacific air masses likely still influenced drought in the region during this time and support that continental air mass movement changes and drought frequency are critical to constrain to predict peatland fire frequency in the future (e.g., Marcisz et al. 2017). However, fire frequency is still comparatively much lower than modern and Holocene fire values for peatlands in western Canadian boreal region, attesting to the site’s location at the transition from a continental to more humid climate (i.e., western boreal fire frequency in peatlands is 2–8 fires 1000·yr⁻¹; Kuhry 1994; Zoltai et al. 1998; Wieder et al. 2009).

FRIs shorten over the last 500 years and are smaller than both the Middle and Late Holocene intervals (Figs. 6 and 7), which we attribute to recent drying as well as decomposition and acrotelm processes. Some northern peatlands have been drying over the last two centuries with climatic and land use changes and are expected to dry out further with anthropogenic climate change and therefore promote increased fire frequency (Turetsky et al. 2015; Swindles et al. 2019). Although these factors may in part explain the trends in fire frequency at the site, the median fire frequency recorded from the pooled dataset in the last 500 years (3.3 fires 1000·yr⁻¹) is over three times higher than the warmer Middle Holocene value, suggesting that factors other than climate are needed to explain large differences. Decomposition can cause the broadening of charcoal peaks and may disguise smaller peaks as a result (Camill et al. 2009). Further, fire severity may impact whether charcoal layers are formed or not in the peat column (Turetsky et al. 2004). Therefore, our Middle Holocene estimate could be a conservative estimate of fire frequency in this interval. Charcoal fragments could also potentially move in the acrotelm, which may cause artificial peaks and increase the apparent number of fires (Turetsky et al. 2004). However, we aimed to minimize this through determining background charcoal in CharAnalysis (van Bellen et al. 2012). Overall, due to these processes, the acrotelm and catotelm fire frequencies are challenging to compare.

aCARs at the site were not significantly different between the Middle and Late Holocene and do not correspond to changes in median fire frequency, suggesting that carbon losses may have been offset by other local factors (Fig. 5). The Middle Holocene has also been linked to higher biomass production with longer growing degree days above zero in northern peatlands (Yu et al. 2010). Higher vegetation inputs into the site between fires may have therefore compensated for some of the loss. Sphagnum macrofossils are also consistently found throughout each core, suggesting rapid recovery from fire (Fig. 3, Magnan et al. 2012). Further, charcoal peaks are not always associated with evidence of peat smoldering (Fig. 3 and Table 1). Therefore, Sphagnum moss was able to protect the deeper soils by preventing water loss in most cases (Benscoter and Wieder 2003; Benscoter et al. 2011). Fire may have also offset carbon losses at the site by stabilizing a fraction of the organic matter through charcoal production (Heinemeyer et al. 2018; Flanagan et al. 2020; Pellegrini et al. 2021). Chemical alteration of peat at the site during fire is supported by the low C:N ratios and high nitrogen content that are linked to some of the charcoal peaks in the cores (Zaccone et al. 2014). Therefore, the records in this study demonstrate the complex interactions between fire and carbon storage that are dependent on local factors, including whether fire directly or indirectly impacts the peat column.

Fire frequency and carbon accumulation relationships vary at local scales. This variation implies that while fire can play a role in removing carbon inputs, regional compilations are clearly needed to predict the relative role of fire on the global peatland carbon balance. Within local studies in boreal peatlands, both positive and negative relationships between fire frequency and LORCAs have been demonstrated (Kuhry 1994; Camill et al. 2009; van Bellen et al. 2012; Magnan et al. 2020) and, depending on the region, local factors including hydrological conditions (van Bellen et al. 2012), vegetation composition (Camill et al. 2009), and decomposition prevention (Heinemeyer et al. 2018) compensated for any impacts on carbon accumulation rates. When sites are pooled from boreal regions that represent a larger gradient in fire frequency and moisture balance conditions seen across the Holocene than at an individual site, we demonstrate that there is a significant negative relationship between fire frequency and LORCA, with considerable variation (Fig. 7). The variation supports that some peatlands will be more resilient to carbon loss from fire than others (Ingram et al. 2019). Resilience has been shown to come from vegetation composition (i.e., Sphagnum moss abundance and recovery; Magnan et al. 2012; Kettridge et al. 2015) and hydrogeological setting (Hokanson et al. 2016; Ingram et al. 2019). At our study site, LORCA values are lower than predicted based on the sites fire frequency in the linear regression model (Fig. 7). Further, the estimated total carbon loss from each coring locality does not fully explain lower LORCAs at the site from median HBL values (18.3 g C m⁻²·yr⁻¹; Packalen and Finkelstein 2014; Table 1). Therefore, the dry conditions at the site needed for higher tree cover also likely contributed further to lowering carbon stocks (Magnan et al. 2020; Beaulne et al. 2021). Therefore, as expected, sites with dry conditions and high tree and shrub cover are more vulnerable to fire and will likely have an even more diminished long-term carbon storage if more frequent fire occurs.

Although only subtle changes in median fire frequency were recorded at the site with modest mean summer temperature changes within the Middle Holocene (Fig. 6; Hargan et al. 2020), these changes also have implications for carbon release across the Holocene. When comparing the median Middle and Late Holocene FRI of the study site, emissions were potentially doubled per unit area between the Middle and Late Holocene, despite minimal annual temperature increases (Hargan et al. 2020). Therefore, characterizing the FRI in HBL peatlands has the potential to constrain the atmospheric inputs from fire on a regional and global scale through the Holocene. However, continued work capturing spatial variability in peatland fire frequency changes across the region is needed to support scaling fluxes from fire (Loisel et al. 2017).