Coherent signature of warming-induced extreme sub-continental boreal wildfire activity 4800 and 1100 years BP

Climate changes are expected to progressively increase extreme wildfire frequency in forests. Finding past analogs for periods of extreme biomass burning would provide valuable insights regarding what the effects of warming might be for tree species distribution, ecosystem integrity, atmospheric greenhouse gas balance, and human safety. Here, we used a network of 42 lake-sediment charcoal records across a ∼2000 km transect in eastern boreal North America to infer widespread periods of wildfire activity in association with past climate conditions. The reconstructed fluctuations in biomass burning are broadly consistent with variations in ethane concentration in Greenland polar ice cores. Biomass burning fluctuations also significantly co-varied with Greenland temperatures estimated from ice cores, at least for the past 6000 years. Our retrospective analysis of past fire activity allowed us to identify two fire periods centered around 4800 and 1100 BP, coinciding with large-scale warming in northern latitudes and having respectively affected an estimated ∼71% and ∼57% of the study area. These two periods co-occurred with widespread decreases in mean fire-return intervals. The two periods are likely the best analogs for what could be anticipated in terms of impacts of fire on ecosystem services provided by these forests in coming decades.


Introduction
The Northern Hemisphere has recently experienced extreme wildfire seasons, which have burned extensive forest areas across the boreal biome. Notably, recordbreaking warm and dry conditions during the summer of 2018 spurred historic wildfire outbreaks in Sweden, ranging across more than 24 310 ha of lands from north of the Arctic Circle to the southernmost county of Scania (San-Miguel-Ayanz et al 2019). The province of British Columbia, Canada, faced the worst wildfire year on record during the summer of 2017, with 1.2 million ha having burned (Ansmann et al 2018, Hu et al 2019, Natural Resources Canada 2018. During the summer of 2010, intense drought and high temperatures caused several hundred wildfires in Russia, resulting in about 5 million hectares of burnt areas (Viatte et al 2013). Such extreme fire activity can have important economic and social consequences, especially when it occurs in populated areas. Communities must face evacuations, mental and physical health problems, damage to infrastructure, and disruptions to business and industry (Bowman et al 2011, Viatte et al 2013, Landis et al 2018. As an example, the 590 000 ha Fort McMurray Horse River wildfire in Alberta (Canada) during 2016 incurred the worst insured losses in Canadian history, destroying 2400 buildings (Statistics Canada 2017, Sankey 2018). Again in Canada, massive wildfires during the summer of 2013 threatened a Cree community in northern Quebec, forcing its evacuation (Erni et al 2017, Sankey 2018.
Human-caused climate changes, associated with warming and drying at high-latitudes, are expected to progressively increase the frequency of such extreme wildfires in the future, with important consequences for species distribution, ecosystem integrity, atmospheric greenhouse gas balance, and human safety (Flannigan et al 2009, Gauthier et al 2015. Fire suppression agencies will be increasingly overwhelmed (Bowman et al 2017, Wotton et al 2017. Recently, an analysis of fire weather climatology from 1979 to 2015, which used the North American Regional Reanalysis dataset and the Canadian Fire Weather Index System, has suggested an ongoing significant, increasing trend in fire-season length in eastern Canada (Jain et al 2017). To date, detection of temporal changes in the frequency (or in the return interval) of major wildfire activity remains a major challenge, owing to the short temporal coverage of circumboreal fire data (usually limited to post-1950), the large interannual and spatial variability in area that is burned, and evolving wildfire monitoring and reporting methods, which tend to mask long-term and subtle changes that can be a consequence of changing climate (Krezek-Hanes et al 2011, Waito et al 2015. The need for a comprehensive assessment of historical wildfire, climate, and vegetation at all scales, including enhanced datasets from historical fire maps and synthesis of climate, vegetation, and fire interactions from past climates, including past warm periods, was identified as national priorities for building the capacity of wildland fire science in the 2020s (Sankey 2018).
Lake-sediment charcoal records are used to reconstruct long-term fire trends and analyze the sensitivity of fire activity to various ecosystem and climatic processes, and can ultimately improve our ability to project future fire impacts (Heyerdahl et al 2008, Girardin et al 2013b, Waito et al 2018. Charcoal particles that are produced during fire are transported by wind and runoff across landscapes and deposited in lake sediments. They have been used to infer trends in typical (average) biomass burning and changes in fire-return intervals (FRIs) from local to hemispheric scale (Marlon et al 2012, Pellatt et al 2015, Waito et al 2015. Here we extend the utility of these records by demonstrating their applicability in distinguishing synchronous fires or fire periods at the sub-continental scale in an attempt to link widespread pulses of biomass burning and rapid warming during preindustrial times. We illustrate our case with the use of a network of 42 lake-sediment charcoal records that cover the last 8400 years in central and eastern boreal North America. The retrospective analysis of past fire activity allowed us to identify signals from two highfire activity periods across the study area, dating from about 4800 and 1100 calibrated years before present (BP), and having occurred during climatic warming at high latitudes.

Data and methods
Forest fire data Charcoal samples from 42 lake sediment cores were obtained from 14 studies sharing a common protocol of charcoal extraction (see the references listed in table S1 for details, and figure S1 is available online at stacks. iop.org/ERL/14/124042/mmedia for study locations). The network extends throughout the eastern North American boreal forest across a west-to-east transect of more than 2000 km, which covers three Canadian provinces (figure 1(a); table S1). Forests in this area developed during the Holocene, which began with the end of the most recent glaciation about 11 650 years BP (Dyke 2004, Walker et al 2018. During this period, ice sheet margins retreated northeastward. As the glacier retreated, and following the disappearance of the proglacial Objibway Lake ca. 8200 BP, forest density increased with fire-adapted conifer tree species Picea becoming, and remaining, the dominant feature in boreal landscapes; together with Abies, broadleaved trees such as Acer were dominant in the mixedwood forests south of the study area (figure S2; Blarquez et al 2015). The average length of the sediment cores was 320 cm (table S1); in most cases, sampling had reached the contact between organic and mineral sediments, allowing fire histories to be traced back to the very beginning of the post-glacial afforestation period. All cores were sliced at contiguous 0.5 or 1.0 cm intervals to discern fire events at high temporal resolution. The analysis of a network of fire activity records of this magnitude allowed us to achieve unequaled spatial and temporal resolution in this part of North America (i.e. Girardin et al 2013a, Blarquez and Aleman 2016). Sediment accumulation chronologies were generated based on AMS radiocarbon dating of terrestrial plant macroremains and bulk gyttja samples, and 210 Pb when available (table S1). 210 Pb measurements were performed from the uppermost 10-20 cm of the cores. 210 Pb values were inferred by measuring the activity of the daughter product, 210 Po, by alpha spectrometry assuming an equal concentration between the two isotopes. The 210 Pb concentrations were interpreted using the constant rate of supply model of 210 Pb accumulation (Appleby and Oldfield 1978). 14 C dates were calibrated using Bayesian models (Parnell et al 2008;Bchron R package) or 'classical' models (Blaauw 2010; CLAM 2.2 program) based on the IntCal13.14 C data set, as presented in original studies (Hua et al 2013, Reimer et al 2013. It has been shown that the two models (Bayesian and 'classical') derive relatively equivalent chronologies (Blaauw 2010) for the studied time period (0-8450 BP; Wang et al 2019) and lakes (figure S3). Charcoal particles were measured and the resulting data were converted into charcoal accumulation rates (CHAR raw ; mm 2 cm −2 yr −1 ) that were based on numerical age-depth models to reconstruct past fire frequencies for each lake (table S1). The dates of local fire events were extracted from the CHAR raw series using CharAnalysis v1.1 software (Higuera et al 2010). CHAR raw series were interpolated (CHAR interpolated ) to equal time steps based on the median time-resolution computed from the sediments for each lake (table S1). CHAR interpolated series were decomposed into a lowfrequency component (CHAR background ) and a highfrequency component (CHAR peak ). CHAR background results from long-distance burning or redeposition processes of charcoal particles that are unrelated to local fire occurrences. CHAR background was estimated by applying the LOWESS-smoother robust technique, and subtracted to CHAR interpolated to obtain the CHAR peak component. CHAR peak was decomposed into two subpopulations: CHAR noise representing variability in sediment mixing and sampling, together with analytical and naturally occurring noise; and CHAR fire representing significant charcoal peaks , Higuera et al 2010. For each peak, we used a Gaussian mixture model to identify the CHAR noise distribution according to a locally-defined threshold. Signal-to-noise indices (Kelly et al 2011) and goodness-of-fit (Brossier et al 2014) were used to evaluate the effectiveness of the discrimination between CHAR fire and CHAR noise and to assess peak detection accuracy by comparing the empirical and fitted noise distributions, respectively (figure S4). Each CHAR peak that exceeded the threshold at the 99th percentile of the noise distribution was assumed to originate from a local fire event (Gavin et al 2006, Higuera et al 2011). Given the short fire rotations in some parts of the study area (e.g. northwestern Quebec, figure 1), we recognize that some peaks may include more than one fire .
Forest fire data from the Canadian National Fire Database (Stocks et al 2002) were used for this study to described current fire activity. The database contains information on the location, date of detection, size (ha), and cause (lightning or human) of all fires that were recorded in Canada, for the period 1959-2016. Starting from the 0.5°×0.5°grid covering the study area (figure 1), we extracted all fires and computed the annual area burned for each cell and year.
FRI and biomass burning FRI is herein defined as the time that has elapsed between two successive fires that are deduced from CHAR fire , and which are assumed to have occurred within a given landscape area. For each lake, fire dates were converted to FRIs by subtracting each date from the previous one.
CHAR raw series were submitted to a min-max normalization to remove biases that could result from differences in taphonomic processes between lakes. Normalized charcoal accumulation (CHAR norm ) was used as a proxy for biomass burning . All FRI and CHAR norm series were extended at both edges using a mirror reflection method in order to cover the beginning and the end of each series.

Weibull smoothing
To avoid overestimating the importance of single fire events, the FRI distributions were smoothed independently for each lake using the two-parameter Weibull probability density distribution (Johnson and Gutsell 1994, Cyr et al 2009) and a 5-observation moving window. The b and c parameters, which respectively corresponded to the shape and scale parameters of the Weibull distribution, were estimated using the fitdist function of the 'fitdistrplus' R package (Delignette-Muller and Dutang 2015). For each lake, we used a gamma (Γ) function (Johnson and Gutsell 1994) to calculate the mean of the Weibull distribution, where within each 5-observation window, the mean FRI was calculated as: Bootstrap replications (1000 iterations) were used to construct equi-tailed (1-2σ) confidence intervals for the estimated mean FRI within each window. For each window, the mean of the 1000 resulting values was defined as the Weibull-modeled mean FRI. Weibull-modeled mean FRIs were then reported along the time series for each lake (Cyr et al 2009).

Kriging
Ordinary kriging (Cressie 1990) was used to spatially extrapolate changes in FRI and CHAR norm during the Holocene to the entire study area at a 0.5°×0.5°r esolution. Kriging is a linear approach to interpolation/extrapolation of responses, where values across the field surface are determined in accordance with a given covariance structure, which is determined from those observations that are present (Cowtan and Way 2014). The reconstructed values vary smoothly and match the observed values at the coordinates of the observations; they approach the global mean as the distance from the nearest observation increases.
For each lake, values of the Weibull-modeled mean FRI and of the CHAR norm were extracted at a 50 year time resolution from the corresponding distributions. A buffer with a 250 km radius was established around the lakes to delimit the perimeter of the study area, upon which a 0.5°×0.5°grid was superimposed. This process was repeated at each time step as the pool of lakes changed with time. Ordinary kriging was realized with nugget=0 and range=150 km, which were kept constant between time steps. The sill parameter was left unspecified in order for the algorithm to optimize its value at each time step. An exponential model was used to fit the semivariograms. Kriging was performed using the autoKrige function of the 'automap' R package (Hiemstra 2009). The resulting FRI krig (which approximates the average FRI, in years, for a specific grid cell during a specific time period) and CHAR krig grids were produced at a 50 year time resolution, and mapped at a 500 year time resolution.
For each 50 year time point, the mean value of all grid cells, along with the corresponding 95% confidence interval, were extracted and reported along the time series to produce the temporal variability of FRI krig and CHAR krig . To do so, we examined the local spatial structure of each 50 year time point by using Moran's I spatial correlograms (moran.test function; R Core Team 2014), and calculated 95% confidence interval after correcting for the effective degrees of freedom (n′) based on lag 1 autocorrelation estimates of Moran's I (see Dale and Fortin 2002). The advantage of averaging the values of the grid cells, rather than the values of the lakes, is that regions where more lakes were sampled were not over-represented.

Extreme biomass burning detection
Extreme biomass burning was defined as periods during which CHAR krig exceeded a detection threshold that was computed from the long-term median (MED) and the median of absolute distances to the median (MAD): with i=1, K , n, for a grid-cell's CHAR krig record x of size n time points (Hampel 1985) (figures S5 and S6). The parameter z is a threshold that was selected. The number of extreme events under analysis was arbitrarily set to approximate the highest 15th percentile (a function of the distance to the median z, which varies by increments of 0.5). For this study, we set z=3. We preferred to use this extreme event detection method on CHAR krig records rather than using local fire events that were detected using CharAnalysis software. Major fires that were identified by CHAR raw peak records can lag actual fire dates by several decades; this is a common occurrence that is related to the timing of charcoal transport and deposition (Higuera et al 2005). The fire events are therefore diluted over time, making the interpolation technique unreliable for such identification of individual events as it averages over noise present in their timing. Alternatively, exceedance of the detection CHAR krig threshold indicates high inputs of CHAR raw into a lake, but does not inform on the occurrence of a single fire.
The sample proportion p, representing the fraction of grid cells k (an integer0) of a given population n′ (an integer>0) that was identified positively as recording an extreme event, was computed for each sampled 50 year period from 8450 to 0 years BP. Binomial proportion confidence intervals (95% CIs) were computed using a Bayesian calculation with an uninformative prior distribution (Brown et al 2001). This uncertainty estimate is important since n′ is unevenly distributed in time, and because small n′ or values of k that are close to 0 or 1 can lead to large errors in the estimation of p. The Bayesian approach that was used here (Jeffreys prior with Beta [1/2, 1/2]) has been shown to be particularly robust under small samples sizes, while also being well-suited for uncertainty estimation under large sample sizes.
The extreme event detection method was repeated on modern area burned data that were extracted from the Canadian National Fire Database, and for each cell over the 1959-2016 period. A year in a given grid cell was deemed as being extreme when the sum of its area burned exceeded the detection threshold set by equation (1) for the particular annual area burned dataset for that particular cell.
Relationship with temperature and Greenland ethane levels Seasonally unbiased and physically constrained Greenland Summit temperatures that covered the Holocene period, and which were inferred from argon and nitrogen isotopes in air trapped within a Greenland ice core (GISP2; Kobashi et al 2017), were used to infer associations between high-latitude temperature variability and biomass burning over the period 8400-0 BP. The Greenland temperature reconstruction has been shown to be representative of hemispheric signals (Kobashi et al 2017). It should be noted that sampling density was not constant through the Holocene, such that uncertainty ranges of the reconstructed temperatures varied with time (Kobashi et al 2017). The relationship between paleoatmospheric ethane in Greenland polar ice cores (Nicewonger et al 2018) and biomass burning was also examined.
Relationships between variables under study were examined using non-parametric stationary bootstrapped correlations (Mudelsee 2003) and wavelet coherence (WTC) analyses (Grinsted et al 2004). All analyses accounted for the presence of serial autocorrelation in time series when testing for significance of relationships. The non-parametric stationary bootstrapped correlation technique resamples blocks of data pairs to account for the presence of serial (auto-) correlation in the time series. The confidence intervals (CI) allow testing whether the correlation between two serially dependent time series is significant. When the confidence interval contains zero, the hypothesis of 'no correlation' cannot be rejected at the 95% level. Where indicated, data were ranked prior to analysis to satisfy normality requirements for the correlation analysis. In contrast, WTC examines the significance and magnitude of coherence (i.e. correlation) between two time series and, furthermore, reveals information about the phase relationship (delays) in a lag-period space (at a time and frequency localization). This analysis of phase relationships is relevant in a context in which the direct correlations between variables may be difficult to assess due to dating uncertainties and offset paleorecords. For each time-series, the required frequency spectrum was constructed using continuous wavelet transform analysis. We used the Paul wavelet (order 4) as it is not very localized in frequency space, while allowing signals that are relatively aperiodic to be detected by the analysis (Moore et al 2007). The statistical significance level (90%) of the WTC against red noise backgrounds was estimated using Monte Carlo generated noise with 1000 surrogate data set pairs (having the same first-order autocorrelation coefficients as the input datasets; see Grinsted et al 2004). WTC was executed using the R package 'biwavelet' version 0.20.11 (Gouhier et al 2018). For the purposes of the correlation and WTC analyses, the paleoclimate data were resampled to a lower resolution of 50 yeartime steps that was compatible with the CHAR records, using spline estimation; for the correlation with ethane levels in Greenland, biomass burning data were resampled to the time resolution of the ethane data (Software Autosignal).
In our next step in the analysis, spatial correlation maps (Von Storch and Zwiers 1999) were created by analysing the relationship between biomass burning, as computed from modern fire statistics, and gridded temperature data from the Climate Research Unit (CRU TS 4.2 data, resolution of 0.5°latitude × 0.5°l ongitude). Temperature data were averaged at the grid level over the April-September season. These analyses were conducted using the KNMI Climate Explorer online research tool; all data were ranked prior to analysis (van Oldenborgh et al 2008). Additionally, spatial correlation maps were created after analysing the relationship between the Greenland temperature reconstruction and CHAR krig at each grid cell.

Results and discussion
Multi-millennial spatiotemporal variability in fire activity Modern fire rotations (i.e. the time required to burn the equivalent of a specified area; period 1959-1999) range from 25 years to over 500 years across the study area ( figure 1(a); Boulanger et al 2014). Of note is that the area has seen an increase in the proportion (p) of its forests that have been affected by extreme biomass burning (i.e. area burned exceeded an empirical extreme detection threshold, see Data and Methods) during the AD 1980s, peaking in the early AD 1990s and then declining ( figure 1(b); proportions were computed using the Canadian National Fire Database, hereafter pĈ NFD ; see Data and Methods). This phenomenon was observed in analyses of both human-ignited and lightning-ignited fires (figure S7). The years 1988 and 1995, for example, had particularly high pĈ NFD (respectively 52% and 57% of grid cells), distributed in all three provinces (figure S8).
We were interested in determining whether there were past analogs of such widespread biomass burning that could be reconstructed from lake-sediment charcoal records. We conducted spatially-interpolated reconstructions for the whole study area of the normalized charcoal (CHAR krig ) accumulation rate, which is representative of the amount of biomass burning, and of the mean FRI (FRI krig ), which is defined as the average time that has elapsed between successive fire events for a specified 0.5°× 0.5°grid cell during a specified 50 year time period.
Holocene distributions of FRI krig (figure 2) and CHAR krig (figure 3) exhibited high variability through time. From ca. 8200 to 7800 BP, the study area experienced low FRI krig and high biomass burning (figure 4). Following a period of fluctuating FRI krig , low biomass burning lasted until about 6000 BP. FRI krig was longer in the western sector and shorter in the eastern sector (figure 2). Between 6000 and 4500 BP, FRI krig decreased and biomass burning increased (figure 2), as has been previously observed in North America (Hu et al 2006, Figure 2. Spatial interpolation of the fire-return interval (FRI krig ) for 42 lakes from 8000 to 0 years before present (BP). Fire-return interval is defined as the elapsed time between two fires that are assumed to have occurred within a landscape area. FRI was extracted for each lake at 500 year intervals, and spatially interpolated using ordinary kriging. The estimated distribution of the Laurentian Ice sheet (Dyke et al 2003) is also illustrated using a conversion from radiocarbon estimated dates to calibrated calendar ages (Reimer et al 2013). Yellow to red colors indicate areas that experienced frequent fires. Ali et al 2009, Girardin et al 2013a. Between 4500 and 2500 BP, fire activity was at its highest (figures 2-4). The amount of biomass burning fluctuated but remained relatively high, on average (figure 4(c)). FRI krig then gradually increased until 1000 BP, especially in the northcentral part of the study area (figure 2), as has been shown by previous studies (Girardin et al 2013b. From 1000 to 250 BP, FRI krig also increased in the south-central part of the study area (figure 2), while remaining relatively high and constant elsewhere. During this period, biomass burning decreased. After 250 BP, we observed an increase in biomass burning and a decrease in FRI krig ( figure 4). A positive trend in biomass burning during 1750-1950 AD, roughly coinciding with the transition from the Pre-Industrial period to the Anthropocene (~1950 AD), was found across 49% of the area, but mostly limited to the western half (figure S9). Notable fluctuations in biomass burning since 1000 BP are broadly consistent with independent estimates of biomass burning emissions that were recently obtained from paleoatmospheric ethane in Greenland polar ice cores (figures 4(c) and S10; data from Nicewonger et al 2018). Biomass burning is one of two important sources of ethane in the preindustrial atmosphere (Nicewonger et al 2018). This demonstrates the relevance of the information that can be deduced from lake-sediment data for large-scale environmental issues: the correlation between the two independent records implies that fire activity in the study area is an important source for fire emissions transport into Greenland (Kehrwald et al 2012).
Embedded within this long-term variability are two periods of very-high biomass burning, each lasting about two centuries (figure 4(d)): from 4950 to 4750 BP and from 1250 to 1050 BP. In both cases, the pp roportions of cells with CHAR krig values exceeding local extreme detection thresholds (hereafter pĈ HAR ; see Data and Methods) were significantly elevated above overall background variability ( figure 4(a)). These CHAR krig periods are not likely to be related to single large or severe fire events (i.e. biomass consumed per unit area), given that they occur simultaneously with decreases of FRI krig (figure 4(a)). This observation implies that multiple fires marked both periods over the affected regions. During the period of 4950 to 4750 BP, pĈ HAR covered 71% of the study area, mostly in the eastern sector (see inset maps of figure 4(d)). CHAR krig values also exceeded the Figure 3. Spatial interpolation of normalized charcoal (CHAR krig ) accumulations for 42 lakes from 8000 to 0 years before present (BP). CHAR krig was extracted for each lake at 500 year intervals, and spatially interpolated using ordinary kriging. The estimated distribution of the Laurentian Ice sheet (Dyke et al 2003) is also illustrated using a conversion from radiocarbon estimated dates to calibrated calendar ages (Reimer et al 2013). Light blue to red colors indicate areas that underwent high biomass burning. detection threshold in 15 individual lake-sediment charcoal records (i.e. 45% of lakes) that were located in eastern and western parts of the study area (see inset maps of figure 4(d)). From 1250 to 1050 BP, pĈ HAR covered 57% of the study area ( figure 4(d)). CHAR krig values exceeding the threshold were also found in 15 of (a) FRI was extracted for each lake at 50 year intervals, spatially interpolated using ordinary kriging, and here averaged across the whole study area. The shading is indicative of the 95% confidence intervals (95% CI) for the averaged quantities, after adjusting the sample size to an effective degree of freedom. The number of lake-sediment charcoal records at a given time is shown in (b). (c) Same as (a) but for the normalized charcoal (CHAR norm ) accumulation. The squares illustrate the distribution of Greenland ice core ethane levels from (Nicewonger et al 2018) and are used as an independent fire proxy for benchmarking (also see figure S10 for a zoom of the last 1000 years). The two records are significantly correlated with r=0.62 and bootstrap 95% CI [0.16, 0.86]. (d) Percent of the study area for which grid-cell CHAR krig exceeded locally-defined extreme detection thresholds at 50 year intervals (pĈ HAR ), with 95% confidence intervals. The gray vertical bars delineate two periods during which a high percent of the study area underwent extreme biomass burning, i.e. 4950-4750 and 1250-1050 years BP. The inset maps show the distribution of those grid cells (outline in black) in which CHAR krig (color shading) exceeded extreme detection thresholds. Also, circles were added to illustrate the distribution of sampled lakes in which CHAR interpolated values, of the specified period, exceeded the extreme detection threshold specific to each lake. the sampled lakes (38% of lakes; figure 4(d)). For this particular period, biomass burning occurred across many regions, stretching from the prairie-boreal ecotone to the easternmost boreal forests.
Aside from these two periods, the period 3100-2850 BP also stands out as having had a high fire activity (figure 4). In particular, we find FRIs among the lowest recorded. It should be noted that although an increase in CHAR krig occurred at the beginning of the Industrial Era, the proportion of cells exceeding the detection limit during the modern period was within the range of background variability ( figure 4(d)).
Signatures of warming-induced extreme fire activity Biomass burning is the result of additive effects of temperature, seasonal drought severity, duration of the snow-free period (equivalent to fire season length), and fuel availability (Molinari et al 2018, Chaste et al 2019). To our knowledge, evidence for a past-temperature influence on spatially widespread periods of high fire activity is lacking for our study area. To examine whether this effect was plausible, we analyzed the correlation between the metric pĈ NFD during 1959-2016 (data from figure 1(b)) and April-September mean land temperatures (see Data and Methods). We found that pĈ NFD co-varied in a positive manner with large-scale mean April-September warming at high latitudes, including over Greenland ( figure 5). This positive covariance between fire-affected areas and Greenland temperature sets the basis for exploring relationships between past widespread boreal fire activity and multi-decadal temperature reconstructions covering the Holocene from Greenland ice cores (Kobashi et al 2017;figure 4(e)). Consistent with the analysis of modern fire and temperature data, wavelet coherence (WTC) analysis between the Greenland  The 90% significance level against red noise is shown as a thick contour. The gray line delineates the 'cone of influence' where zero padding has reduced the variance. The relative phase relationship is shown by arrows, with in-phase pointing right, antiphase pointing left, temperature leading pĈ HAR by 90°(1/4k) pointing down, and pĈ HAR leading temperature by 90°(1/4k) pointing up. In-phase and antiphase relationships may be interpreted respectively as positive and negative correlations in time and period (Grinsted et al 2004). temperature reconstruction and the pĈ HAR record indicated areas of significant and almost in-phase correlation at periodicities of>800 years per cycle since 6000 BP (i.e. pĈ HAR generally lagging behind temperature; figure 6). High power was also noticeable during periods approximating those of high biomass burning, i.e. ∼4250 to 5250 and ∼1750 to 1250 BP at periodicities of about 400 years per cycle (albeit given the implied coarseness of the temporal resolution of the fire record, information that can be obtained is likely to be limited at such frequency). Non-parametric stationary bootstrapped correlations that were computed for the last 6000 years BP corroborated the WTC correlation for that period with r=0.20 and 95% CI  Porter et al 2019) remains to be elucidated. It could be hypothesized that relatively low biomass burning from 7500 to 6500 BP was the result of cold air masses associated with the remaining ice sheet in northern Quebec (Dyke 2004, Marsicek et al 2018, Liu et al 2014 (figures 1 and S13), humid conditions and high lake levels across large areas of eastern North America (Hély et al 2010, their figure 3(e)), and higher biomass of less fire-prone deciduous tree species, notably Betula and Populus (Girardin et al 2013a, Blarquez and Aleman 2016, figure S2). All these factors could have contributed to reduce spreading of crown fires.
It should be noted that the relationship between CHAR krig and temperature is not homogeneous across the study area ( figure 7). It is a dominant feature of the easternmost sectors, and during the last 2000 years or so, it encompassed a larger part of the study area ( figure 7(c), also see figures S11 and S12). Note also, that the pattern of temperature and CHAR krig correlation over the last 2000 years (figure 7(c)) is similar to the pattern of linear trends of CHAR krig from 200 to 0 BP (figure S9). It is likewise similar to the pattern of correlation between pĈ NFD and gridded mean April-September temperatures for 1959-2016 (figure 5). From these observations, it appears to be a region in the center of the study area, southeast of James Bay, which emerges from the ensemble where CHAR krig temporal trajectories do not co-occur with high latitude temperatures. It was previously postulated that precipitation patterns and changes in the composition of forest landscapes were the main factors that would influence the fire activity trajectories below James Bay, as opposed to direct effects of warming on fire (Girardin et al 2009(Girardin et al , 2013a. The decoupling of temperature and CHAR krig trajectories in that particular area adds credibility to this thesis. Two periods of high sub-continental biomass burning and pĈ HAR (i.e. 4950-4750 BP and 1250-1050 BP) are found during some of the periods of warmest Holocene temperatures in Greenland (respectively about +2.4°C and +3.4°C above the 1951-1980 baseline; figure 4(e); Kobashi et al 2017). Both periods coincide with positive temperature anomalies that have been inferred from chironomid assemblages in eastern boreal North American lakes (∼+2.1°C and 2.2°C above baseline (Bajolle et al 2018, figure S11). Although the two high wildfire periods that have been identified herein had high temperatures in common, they developed under different precipitation and vegetation contexts. On one hand, the 4950-4750 BP episode occurred during a warm and dry period . To this may be added a third period, i.e. 3100-2850 BP and marked by low FRIs, which could fit into the warm anomaly of about 3400-2800 BP.
Our results covering past and modern times thus indicate that forest fires in the study area can be responsive to rapid warming at high latitudes. This is all the more important considering the changes that are currently underway in the North Atlantic basin. The pĈ NFD exhibits a drop in 1992-1993, consistent with high-latitude cooling caused by the eruption of Mount Pinatubo (Parker et al 1996), followed by an increase until 2005 that is consistent with Greenland warming (figure S14). For the most recent 11 years on record (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016), Greenland temperature has exhibited a slightly decreasing trend consistent with northern North Atlantic-wide cooling during that period (figure S15) and a slowdown of the Atlantic Ocean overturning circulation (Thornalley et al 2018). The proportion of forest that was affected by extreme biomass burning also declined during this period (figures 1 and S7). Clearly, a better understanding of the future decadal temperature variations in the Atlantic (e.g. Liu et al 2017) is of paramount importance in improving future boreal wildfire predictions in eastern boreal North America. Contextually, the spatial/temporal dimensions of warming and fire therein are likely to differ from other boreal regions of western North America where temperatures since the mid-20th century appear to be exceptionally warm in comparison with the last 13 000 years (

Conclusion
Our retrospective analysis, which was based upon a network of lake-sediment charcoal records, allowed us to identify two periods of widespread biomass burning across eastern North American boreal forests at 4950-4750 BP and 1250 to 1050 BP. Both periods coincided with positive temperature anomalies over Greenland (+2°C to +3°C above baseline) and are likely the best analogs for what could be anticipated in terms of impacts on fire in coming decades (midcentury) across much of eastern boreal North America as a consequence of human-caused climate changes (approximately +2°C to +4°C above the 1951-1980 baseline by 2040 for that particular area according to Chaste et al 2019). On the other hand, it is not to be ruled out that end-of-this-century conditions could exceed anything experienced in the Holocene (including the warm early-mid Holocene) and that the best analogs could likely come from more ancient periods (Burke et al 2018). We recommend that future studies aiming at clarifying the spatial/temporal dimensions of a high fire period in the past focus on very-high resolution charcoal and pollen collections targeting the sediment horizons spanning such events. Their detailed documentation could provide valuable insights regarding the effects of widespread burning on vegetation composition, boreal forest dynamics, and the carbon cycle. The capacity of boreal conifer tree species to maintain themselves, for boreal carbon stocks to recover, and for ecosystem services to be maintained in face of multiple fire disturbances could be better evaluated on the basis of these empirical data. This would also offer new opportunities to benchmark dynamic global vegetation model simulations. Such comprehensive analysis of past climate change impacts on fire activity, including past warm periods, is needed to inform future predictions of how changes in extreme wildfires could transform existing forest ecosystems (Sankey 2018).

Acknowledgments
This work was made possible thanks to financial support from the Strategic Grant and Discovery Programs of the Natural Sciences and Engineering Research Council of Canada (NSERC), and the Canadian Forest Service. We thank William J Parsons for the editing of the manuscript, and two anonymous reviewers for their comments on an earlier version.

Data availability statement
The data that support the findings of this study are available from the corresponding author upon reasonable request.