Evidence for fire in the Pliocene Arctic in response to amplified temperature

The mid-Pliocene is a valuable time interval for investigating equilibrium climate at current atmospheric CO₂ concentrations because atmospheric CO₂ concentrations are thought to have been comparable to the current day and yet the climate and distribution of ecosystems were quite different. One intriguing, but not fully understood, feature of the early to mid-Pliocene climate is the amplified Arctic temperature response and its impact on Arctic ecosystems. Only the most recent models appear to correctly estimate the degree of warming in the Pliocene Arctic and validation of the currently proposed feedbacks is limited by scarce terrestrial records of climate and environment. Here we reconstruct the summer temperature and fire regime from a subfossil fen-peat deposit on west–central Ellesmere Island, Canada, that has been chronologically constrained using cosmogenic nuclide burial dating to $\mathrm{3.9}+\mathrm{1.5}/-\mathrm{0.5}$ Ma.

The estimate for average mean summer temperature is 15.4±0.8 ^∘C using specific bacterial membrane lipids, i.e., branched glycerol dialkyl glycerol tetraethers. This is above the proposed threshold that predicts a substantial increase in wildfire in the modern high latitudes. Macro-charcoal was present in all samples from this Pliocene section with notably higher charcoal concentration in the upper part of the sequence. This change in charcoal was synchronous with a change in vegetation that included an increase in abundance of fire-promoting Pinus and Picea. Paleo-vegetation reconstructions are consistent with warm summer temperatures, relatively low summer precipitation and an incidence of fire comparable to fire-adapted boreal forests of North America and central Siberia.

To our knowledge, this site provides the northernmost evidence of fire during the Pliocene. It suggests that ecosystem productivity was greater than in the present day, providing fuel for wildfires, and that the climate was conducive to the ignition of fire during this period. The results reveal that interactions between paleo-vegetation and paleoclimate were mediated by fire in the High Arctic during the Pliocene, even though CO₂ concentrations were similar to modern values.

1 Introduction

Current rates of warming in the Canadian Arctic are now roughly triple the rate of global warming (Bush and Lemmen, 2019). Since 1850, global land surface temperatures have increased by approximately 1.0 ^∘C, whereas circum-Arctic land surface temperatures have increased by >2.0 ^∘C (Jones and Moberg, 2003; Francis and Skific, 2015). Such Arctic amplification of temperatures has also occurred during other warm climate anomalies in Earth's past. Paleoclimate records from the Arctic indicate that the change in Arctic summer temperatures during past global warm periods was 3–4 times larger than global temperature change (Miller et al., 2010). While earth system models (ESMs) have been able to provide fairly accurate predictions of the modern amplification of Arctic temperatures hitherto observed for some time (Marshall et al., 2014), they have only recently implemented mechanisms that simulate Arctic amplification of temperature for past warm periods such as the Pliocene (2.6–5.3) with a convincing pattern of seasonality (Zheng et al., 2019). The success of earlier models at capturing modern warming, contrasted with the additions needed to simulate the Pliocene Arctic temperatures, suggest that the array of fast and slow feedback mechanisms have not fully manifested themselves for the modern Arctic, and perhaps there are still further feedback mechanisms we are yet to understand and implement in climate models.

Figure 1 Global temperatures and atmospheric CO₂ concentration spanning the last 5 million years of Earth's history. Mean annual temperatures (MATs) are inferred from compiled δ¹⁸O foraminifera data (Lisiecki and Raymo, 2005) and plotted as anomalies from present values (a). Modern atmospheric CO₂ measurements (NOAA/ESRL) and ice core observations from EPICA (Luthi et al., 2008) are compared with proxy estimates (b; see Table S1) for the Pliocene epoch, indicated with beige shading. Smoothed curves have been fit to highlight trends in pCO₂ and temperature during the Pliocene.

The Pliocene is an intriguing climatic interval that offers important insights into climate feedbacks. Atmospheric CO₂ concentrations were, at times, as high as modern ones (Fig. 1), but generally show a decreasing trend throughout the Pliocene (Haywood et al., 2016; Pagani et al., 2010; Royer et al., 2007; Stap et al., 2016), Although CO₂ estimates from different methods do not converge, the modeled direct effects of these CO₂ discrepancies appear to be small (Feng et al., 2017). Of additional importance for comparability to the modern climate system, continental configurations were similar to present (Dowsett et al., 2016). While global mean annual temperatures (MATs) during the Pliocene were only ∼3 ^∘C warmer than in the present day, Arctic land surface MATs may have been as much as 15 to 22 ^∘C warmer (Ballantyne et al., 2010; Csank et al., 2011a, b; Fletcher et al., 2017). Further, Arctic sea surface temperatures may have been as much as 10 to 15 ^∘C warmer than modern ones (Robinson, 2009), and sea levels were approximately 25 m higher than present (Dowsett et al., 2016). As a result, the Arctic terrestrial environment was significantly different from today, with boreal ecosystems at much higher latitudes (Salzmann et al., 2008). These changes in vegetation due to climate may have also provided further important feedbacks to Arctic temperatures (e.g., Otto-Bliesner and Upchurch Jr., 1997).

To advance our understanding of Arctic ecosystem response and feedback to temperature amplification during past warm intervals in Earth's history, this investigation targets an exceptionally well-preserved Arctic sedimentary sequence to simultaneously reconstruct summer temperature, vegetation and fire from a single site.

2 Methods

2.1 Site description

To investigate the environment and climate of the Pliocene Arctic we focused on the Beaver Pond (BP) fossil site, located at 78^∘33^′ N (Fig. 2) on Ellesmere Island. The stratigraphic section located at ∼380 m above sea level (m a.s.l.) today includes unconsolidated bedded sands and gravels and rich organic layers including a fossil-rich peat layer, up to 2.4 m thick, with sticks gnawed by an extinct beaver (Dipoides spp.). The assemblage of fossil plants and animals at BP has been studied extensively to gain insight into the past climate and ecology of the Canadian High Arctic (Ballantyne et al., 2006; Csank et al., 2011a, b; Fletcher et al., 2017; Mitchell et al., 2016; Rybczynski et al., 2013; Tedford and Harington, 2003; Wang et al., 2017). Previous paleoenvironmental evidence suggests the main peat unit is a rich fen deposit with a neutral to alkaline pH, associated with open water (Mitchell et al., 2016), likely a lake edge fen or shallow lake fen, within a larch-dominated forest–tundra environment (Matthews and Fyles, 2000), not a low pH peat bog. While the larch species identified at the site, Larix groenlandia, is extinct (Matthews and Fyles, 2000), many other plant remains are Pliocene examples of taxa that are extant (Fletcher et al., 2017).

The fen-peat unit examined in this study was sampled in 2006 and 2010. The main sequence examined across the methods used in this study includes material from Unit II, the entire span of Unit III, and material from Unit IV sampled from Section A as per Mitchell et al. (2016; Fig. S1 in the Supplement; see Mitchell et al., 2016, Fig. 5), with a total sampled profile of 1.65 m. Unit III has been estimated to represent ∼ 20 000 years of deposition based on modern northern fen accumulation rates (Mitchell et al., 2016). The charcoal counts and measurements from this locality were based on 31 sample layers from the 2006 field campaign, while the temperature estimates from specific bacterial membrane lipids were taken from 22 of the sample layers collected in 2006 and an additional 12 samples collected in 2010. The same samples from the 2006 season were analyzed for mean summer temperature and char count where contents of the sample allowed. Pollen was tabulated from 10 samples from the 2006 sequence, located at different stratigraphic depths.

2.2 Geochronology

While direct dating of the peat was not possible, we were able to establish a burial age for fluvial sediments deposited approximately 4–5 m above and 30 m to the southwest of the peat. We used a method based on the ratio of isotopes produced in quartz by secondary cosmic rays. The cosmogenic nuclide burial dating approach measures the ratio of cosmogenic ²⁶Al ($t\mathrm{1}/\mathrm{2}=\mathrm{0.71}$ Ma) and ¹⁰Be ($t\mathrm{1}/\mathrm{2}=\mathrm{1.38}$ Ma) in quartz sand grains that were exposed on hillslopes and alluvium prior to final deposition at BP. Once the quartz grains are completely shielded from cosmic rays, the ratio of the pair will predictably decrease because ²⁶Al has double the radio-decay rate of ¹⁰Be. In 2008, four of the medium- to coarse-grained quartz samples were collected from a vertical profile of planar cross-bedded fluvial sands between 8.7 and 10.4 m below the overlying till surface. The samples were 5 cm thick, separated by an average of 62 cm, and should closely date the peat (the sandy braided stream beds represent on the order of ∼10⁴ years from the top of the peat to the highest sample). Quartz concentrates were extracted from the arkosic sediment using Frantz magnetic separation, heavy liquids and differential leaching with HF in ultrasonic baths. When sample aliquots reached aluminum concentrations <100 ppm (inductively coupled plasma–optical emission spectrometry, ICP-OES) as a proxy of feldspar abundance, the quartz concentrate was subjected to a series of HF digestion and rinsing steps to ensure that more than 30 % of the quartz had been dissolved to remove meteoric ¹⁰Be. Approximately 200 mg of Be extracted from a Homestake Gold Mine beryl-based carrier was added to 150 g of each quartz concentrate (no Al carrier was needed for these samples). Such large quartz masses were digested because of the uncertainty in the abundance of the faster decaying isotope. Following repeated perchloric-acid dry-downs to remove unreacted HF, pH-controlled precipitation, column chemistry ion chromatography to extract the Be and Al ions, precipitation in ultrapure ammonia gas, and calcination at temperatures above 1000 ^∘C in a Bunsen flame for 3 min, oxides were mixed with equal amounts of niobium and silver by volume. These were packed into stainless steel targets for measurement at Lawrence Livermore National Laboratory's accelerator mass spectrometer (AMS). Uncertainty estimates for ²⁶Al∕¹⁰Be were calculated as 1σ by combining AMS precision with geochemistry errors in quadrature. For a complete, detailed description of terrestrial cosmogenic nuclide (TCN) methods, see Rybczynski et al. (2013). The ages provided here are updated from Rybczynski et al. (2013) by using more recent production rate information and considering the potential for increasing exposure to deeply penetrating muons during the natural post-burial exhumation at BP.

Figure 2 Map of the Canadian Arctic Archipelago, highlighting the location of the Beaver Pond site (black star; 78^∘33^′ N; 82^∘25^′ W) and Eureka Climate Station (gray star; 80^∘13^′ N, 86^∘11^′ W – used for modern climate comparison) on west–central Ellesmere Island.

2.3 Paleotemperature reconstruction

Paleotemperature estimates were determined based on the distribution of fossilized, sedimentary membrane lipids known as branched glycerol dialkyl glycerol tetraethers (brGDGTs) that are well preserved in peat bogs, soils and lakes (Powers et al., 2004; Weijers et al., 2007c). These unique lipids are thought to be synthesized by a wide array of Acidobacteria within the soil (Sinninghe Damsté et al., 2011, 2014) and presumably other bacteria (Sinninghe Damsté et al., 2018) in soils and peat bogs but also in aquatic systems. Previously, it has been established that the degree of methyl branching (expressed in the methylation index of branched tetraethers; MBT) is correlated with mean annual air temperature (MAAT), and the relative amount of cyclopentane moieties (expressed in the cyclization index of branched tetraethers; CBT) has been shown to correlate with both soil pH and MAAT (Weijers et al., 2007b). Because of the relationship of the distribution of these fossilized membrane lipids with these environmental parameters, the distribution of these membrane lipids has been used for paleoclimate applications in different environments including coastal marine sediments (Bendle et al., 2010; Weijers et al., 2007a), peats (Ballantyne et al., 2010; Naafs et al., 2017), paleosols (Peterse et al., 2011; Zech et al., 2012) and lacustrine sediments (Loomis et al., 2012; Niemann et al., 2012; Pearson et al., 2011; Zink et al., 2010). In this study we reconstruct mean summer air temperature (MST), using a modified version of a calibration that was developed by Pearson et al. (2011) and is based on 90 core top lacustrine sediment samples from diverse climates and geographical areas.

Improved separation methods (Hopmans et al., 2016) have recently led to the separation and quantification of the 5- and 6-methyl brGDGT isomers that used to be treated as one since the 6-methyl isomers were co-eluting with the 5-methyl isomers (De Jonge et al., 2013). This has led to the definition of new indices and improved MAAT calibrations based on the global soil (De Jonge et al., 2014), peat (Naafs et al., 2017) and African lake (Russell et al., 2018) data sets.

Sediment samples were freeze-dried and then ground and homogenized with a mortar and pestle. Next, using the Dionex^™ accelerated solvent extractor (ASE), 0.5–1.0 g of sediment was extracted with the solvent mixture of dichloromethane (DCM) : methanol (9:1, v∕v) at a temperature of 100 ^∘C and a pressure of 1500 psi (5 min each) with 60 % flush and 60 s purge. The Caliper Turbovap^® LV was utilized to concentrate the collected extract, which was then transferred using DCM and dried over anhydrous Na₂SO₄ before being concentrated again under a gentle stream of N₂ gas. To quantify the amount of GDGTs, 1 µg of an internal standard (C46 GDGT; Huguet et al., 2006) was added to the total lipid extract. Then, the total lipid extract was separated into three fractions using hexane : DCM (9:1, v:v) for the apolar fraction, hexane : DCM (1:1, v:v) for the ketone fraction and DCM : MeOH (1:1, v:v) for the polar fraction, using a column composed of Al₂O₃, which was activated for 2 h at 150 ^∘C. The polar fraction, which contained the GDGTs, was dried under a steady stream of N₂ gas and weighed before being redissolved in hexane : isopropanol (99:1, v:v) at a concentration of 10 mg mL⁻¹ and subsequently passed through a 0.45 µm polytetrafluoroethylene (PTFE) filter. Finally, the polar fractions were analyzed for GDGTs by ultra-high performance liquid chromatography – atmospheric pressure positive ion chemical ionization – mass spectrometry (UHPLC-APCI-MS) using the method described by Hopmans et al. (2016). The polar fractions of some samples were re-run on the UHPLC-APCI-MS multiple times, and the average fractional abundances of the brGDGTs was determined.

For the calculation of brGDGT-based proxies, the brGDGTs are specified by the Roman numerals as indicated in Fig. S2. The 6-methyl brGDGTs are distinguished from the 5-methyl brGDGTs by a prime. The novel indices, including MBT${}_{\mathrm{5}\mathrm{Me}}^{\prime }$, are based on just the 5-methyl brGDGTs and the CBT^′ that was used to calculate the pH (De Jonge et al., 2014):

$$\begin{array}{ll}{\displaystyle }{\mathrm{MBT}}_{\mathrm{5}\mathrm{Me}}^{\prime }& {\displaystyle }=\left(\right[\mathrm{Ia}]+[\mathrm{Ib}]+[\mathrm{Ic}\left]\right)/\left(\right[\mathrm{Ia}]+[\mathrm{Ib}]+[\mathrm{Ic}]+[\mathrm{IIa}],\\ \text{(1)}& {\displaystyle }& {\displaystyle }+\left[\mathrm{IIb}\right]+\left[\mathrm{IIc}\right]+\left[\mathrm{IIIa}\right]+\left[\mathrm{IIIb}\right]+\left[\mathrm{IIIc}\right]),{\displaystyle }{\mathrm{CBT}}^{\prime }& {\displaystyle }{=}^{-\mathrm{10}}\log \left[\right(\left[\mathrm{Ic}\right]+\left[{\mathrm{IIa}}^{\prime }\right]+\left[{\mathrm{IIb}}^{\prime }\right]+\left[{\mathrm{IIc}}^{\prime }\right]+\left[{\mathrm{IIIa}}^{\prime }\right],\\ \text{(2)}& {\displaystyle }& {\displaystyle }+\left[{\mathrm{IIIb}}^{\prime }\right]+\left[{\mathrm{IIIc}}^{\prime }\right])/(\left[\mathrm{Ia}\right]+\left[\mathrm{IIa}\right]+\left[\mathrm{IIIa}\right]\left)\right].\end{array}$$

The square brackets denote the fractional abundance of the brGDGT within the bracket relative to the total brGDGTs.

The distributions of aquatically produced brGDGTs in the lake calibration developed by Pearson et al. (2011) were used to determine MST. When this calibration is used, the fractional abundances of IIa and IIa^′ must be summed because these two isomers co-eluted under the chromatographic conditions used by Pearson et al. (2011):

$$\begin{array}{ll}{\displaystyle }\mathrm{MST}\left({}^{\circ }\mathrm{C}\right)& {\displaystyle }=\mathrm{20.9}+\mathrm{98.1}\times \left[\mathrm{Ib}\right]-\mathrm{12}\times \left(\right[\mathrm{IIa}]+[{\mathrm{IIa}}^{\prime }\left]\right)\\ \text{(3)}& {\displaystyle }& {\displaystyle }-\mathrm{20.5}\times \left[\mathrm{IIIa}\right]\mathrm{RMSE}=\mathrm{2.0}{}^{\circ }\mathrm{C}.\end{array}$$

MAAT and surface water pH were also calculated using a novel calibration created using sediments from East African lakes analyzed with the novel chromatography method and based upon MBT${}_{\mathrm{5}\mathrm{Me}}^{\prime }$ (Russell et al., 2018).

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{MAAT}=-\mathrm{1.2141}+\mathrm{32.4223}\cdot {\mathrm{MBT}}_{\mathrm{5}\mathrm{Me}}^{\prime }\\ \text{(4)}& {\displaystyle }& {\displaystyle }\mathrm{RMSE}\mathrm{of}\mathrm{2.44}{}^{\circ }\mathrm{C}{\displaystyle }& {\displaystyle }\mathrm{Surface}\mathrm{water}\mathrm{pH}=\mathrm{8.95}+\mathrm{2.65}\cdot {\mathrm{CBT}}^{\prime }\\ \text{(5)}& {\displaystyle }& {\displaystyle }\mathrm{RMSE}\mathrm{of}\mathrm{0.80}\end{array}$$

Analytical error (±0.38 ^∘C) was estimated as the average standard deviation of the duplicates run on 18 of the samples from throughout the section.

2.4 Vegetation and fire reconstruction

For charcoal, a total of thirty 2 cm³ samples were taken at 5 cm intervals from 380 to 381.45 m a.s.l. at the BP site, with an additional 2 cm³ sample collected at 381.65 m a.s.l. All samples were deflocculated using sodium hexametaphosphate and passed through 500, 250 and 125 µm nested mesh sieves. The residual sample caught on each sieve was then collected in a gridded petri dish and examined using a stereomicroscope at 20–40× magnification to obtain charcoal concentration (fragments cm⁻³). Charcoal area (mm² cm⁻³) was measured for each sample using specialized imaging software from Scion Corporation. For a detailed description of methods, see Brown and Power (2013).

Vegetation was reconstructed using pollen and spores (herein pollen) at selected elevations chosen to capture upper and lower sections of the elevation profile and also elevations that correspond to changes in charcoal. The sample depths selected for pollen analyses were 380.3–380.4, 381.10–381.25 and 381.35–381.45 m a.s.l. Samples were processed using standard approaches (Moore et al., 1991), whereby 1 cm³ sediment subsamples were treated with 5 % KOH to remove humic acids and break up the samples. Carbonates were dissolved using 10 % HCl, whereas silicates and organics were removed by HF and acetolysis treatment, respectively. Pollen slides were made by homogenizing 35 µL of residue, measured using a single-channel pipette, with 15 µL of melted glycerin jelly. Slides were counted using a Leica DM4000 B LED compound microscope at 400–630× magnification. A reference collection and published keys (McAndrews et al., 1973; Moore et al., 1991) aided identification.

In addition to tabulating pollen and charcoal, a list of plant taxa derived from Beaver Pond was previously compiled in Fletcher et al. (2017). Extant species from this list were selected and their modern occurrences extracted from the Global Biodiversity Information Facility (GBIF.org, 2017). Observation data were grouped by 5^∘ latitude 5^∘ longitude grids cells, and the shared species count calculated using R (R Core Team, 2016). Modern fire frequency was mapped using the MODIS 6 Active Fire Product. The fire pixel detection count per day, within the same 5^∘ latitude 5^∘ longitude grids cells was tabulated over the 10 years of 2006–2015 and standardized by area of the cell. The modern climate maps were generated using data from WorldClim 1.4 (Hijmans et al., 2005). The values for the bioclimatic variables mean temperature of the warmest quarter (equivalent to MST) and precipitation of the warmest quarter (summer precipitation) were also averaged by grid cell. The shared species count, climate values and fire day detections were mapped to the northern polar stereographic projection in ArcMap 10.1.

3 Results

3.1 Geochronology

The burial dating results with ²⁶Al∕¹⁰Be in quartz sand at 10 m below modern depth provides four individual ages. From shallowest to deepest, the burial ages are $\mathrm{3.6}+\mathrm{1.5}/-\mathrm{0.5}$, $\mathrm{3.9}+\mathrm{3.7}/-\mathrm{0.5}$, $\mathrm{4.1}+\mathrm{5.8}/-\mathrm{0.4}$ and $\mathrm{4.0}+\mathrm{1.5}/-\mathrm{0.4}$ Ma (Table S2 in the Supplement), with an unweighted mean age of 3.9 Ma. The convoluted probability distribution function yields a maximum probability age of 4.5 Ma. Unfortunately, the positive tails of the probability distribution functions of two of the samples exceeds the radio-decay saturation limit of the burial age. Therefore, their probability distributions do not reflect the actual age probabilities and uncertainty. Given the positive tail in the probability distribution functions and the inability to convolve all samples, we recommend using the unweighted mean age, 3.9 Ma, with an uncertainty of $+\mathrm{1.5}/-\mathrm{0.5}$ Ma as indicated by the two samples with unsaturated limits. Despite the apparent upward younging of the individual burial ages, the 1σ uncertainties overlap, rendering the samples indistinguishable.

3.2 Paleotemperature estimates

3.2.1 Provenance of branched GDGTs

Previously, brGDGT-derived MAAT estimates ($-\mathrm{0.6}\pm \mathrm{5.0}$ ^∘C) from BP sediments were developed using the older chromatography methods that did not separate the 5- and 6-methyl brGDGTs and a soil calibration (Ballantyne et al., 2010). In marine and lacustrine sediments, bacterial brGDGTs were thought to originate predominantly from continental soil erosion arriving in the sediments through terrestrial runoff. More recent studies, however, have indicated that aquatically produced brGDGTs could be affecting the distribution of the sedimentary brGDGTs and thus the temperature estimates based upon them (Warden et al., 2016; Zell et al., 2013; Zhu et al., 2011). Since the discovery that sedimentary brGDGTs can have varying sources, different calibrations have been developed depending on the origin of the brGDGTs, i.e., soil calibration (De Jonge et al., 2014), peat calibration (Naafs et al., 2017) and aquatic calibrations (i.e., Foster et al., 2016; Pearson et al., 2011; Russell et al., 2018). Therefore, several studies have recommended that the potential sources of the sedimentary brGDGTs should be investigated before attempting to use brGDGTs for paleoclimate applications (De Jonge et al., 2015; Warden et al., 2016; Yang et al., 2013; Zell et al., 2013). In this study, we examine the distribution of brGDGTs in an attempt to determine their origin and consequently the most appropriate calibration to utilize in order to reconstruct temperatures from the BP sediments.

Branched GDGTs IIIa and IIIa^′ on average had the highest fractional abundance of the brGDGTs detected in the BP sediments (see Fig. S2 for structures; Table S4). A previous study that plotted the fractional abundances of the tetra-, penta- and hexamethylated brGDGTs together in a ternary diagram established that brGDGTs derived from soils lie within a distinct area of the plot (Sinninghe Damsté, 2016). To assess whether the brGDGTs in the BP deposit were predominantly derived from soils, we compared the fractional abundances of the tetra-, penta- and hexamethylated brGDGTs in the BP sediments to those from modern data sets in a ternary diagram (Fig. 3). Since the contribution of brGDGTs from either peat or aquatic production could affect the use of brGDGTs for paleoclimate application, in addition to comparing the samples to the global soil data set (De Jonge et al., 2014), peat and lacustrine sediment samples were added into the ternary plot to help elucidate the provenance of brGDGTs in the BP sediments. According to Sinninghe Damsté (2016), it is imperative to only compare samples in a ternary diagram like this where all of the data sets were analyzed with the novel methods that separate the 5- and 6-methyl brGDGTs since the improved separation can result in an increased quantification of hexamethylated brGDGTs. Recently, samples from East African lake sediments were analyzed using these new methods (Russell et al., 2018) and so these samples were included in the ternary plot for comparison (Fig. 3). Although the lakes from the East African data set are all from a tropical area, they vary widely in altitude and, thus, in MAAT. We separated them into three categories by MAAT (lakes >20 ^∘C, lakes between 10 and 20 ^∘C, and lakes <10 ^∘C). By comparing all the samples in the ternary plot, it was evident that the BP samples plotted closest to the lacustrine sediment samples from regions in East Africa with a MAAT <10 ^∘C, suggesting that the provenance of the majority of the brGDGTs from the BP sediments was not soil or peat but lacustrine aquatic production.

Figure 3 A ternary plot illustrating the fractional abundances of the tetra- (Ia–c), penta (IIa–c and II^′a–c) and hexamethylated (IIIa–c and III^′a–c) brGDGTs. The global soil data set (open circles; De Jonge et al., 2014), the global peat samples (green circles; Naafs et al., 2017) and lake sediments from East Africa (black circles indicate samples from lakes >20 ^∘C, red circles indicate samples from lakes between 10 and 20 ^∘C and orange circles designate samples from lakes <10 ^∘C; Russell et al., 2018) are included for comparison with the Beaver Pond sediments (blue circles; this study).

The average estimated surface water pH for the BP sediments (8.6±0.2) calculated using Eq. (5) is within the 6–9 range typical of lakes and rivers (Mattson, 1999). This value is near the upper limit of rich fens characterized by the presence of S. scorpioides (Kooijman and Westhoff, 1995; Kooijman and Paulissen, 2006) and is higher than what would be expected for peat-bog sediments that are acidic (pH 3–6; Clymo, 1964) and which constitute most of the peats studied by Naafs et al. (2017). A predominant origin from lake aquatic production is in keeping with previous interpretation of the paleoenvironment of the BP site, which was at least at times covered by water as evidenced by freshwater diatoms, fish remains and gnawed beaver sticks in the sediment (Mitchell et al., 2016).

3.2.2 Aquatic temperature transfer function

Since there is evidence that the majority of the brGDGTs in the BP sediments are aquatically produced, an aquatic transfer function was used for reconstructing temperature. When we apply the African lake calibration (Eq. 4), the resulting estimated MAAT for BP is 7.1±1.0 ^∘C (mean ± standard deviation). This value is high compared to other previously published estimates from varying proxies, which have estimated MAAT in this region to be in the range of −5.5 to 0.8 ^∘C (Ballantyne et al., 2010, 2006; Csank et al., 2011a, b; Fletcher et al., 2017). A concern when applying this calibration is that it is based on lakes from an equatorial region that does not experience substantial seasonality, whereas the Pliocene Arctic BP site did experience substantial seasonality (Fletcher et al., 2017). Biological production (including brGDGT production) in BP was likely skewed towards summer and, therefore, summer temperature has a larger influence on the reconstructed MAAT. Unfortunately, no global lake calibration set using individually quantified 5- and 6-methyl brGDGTs is available yet. Therefore, to calculate MST (Eq. 3) we applied the aquatic transfer function developed by Pearson et al. (2011) by combining the individual fractional abundances of the 5- and 6-methyl brGDGTs. The Pearson et al. (2011) calibration was based on a global suite of lake sediments including samples from the Arctic, thus covering a greater range of seasonal variability. The resulting average estimated MST was 15.4±0.8 ^∘C (mean ± 1 standard deviation, n=34 samples), with temperatures ranging between 14.1 and 17.4 ^∘C (Fig. 4). This is in good agreement with recent estimates based on Climate Reconstruction Analysis using Coexistence Likelihood Estimation (CRACLE; Fletcher et al., 2017) that concluded that MSTs at BP during the Pliocene were approximately 13 to 15 ^∘C.

Figure 4 Reconstruction of mean summer temperature and fire for the Canadian High Arctic during the Pliocene. Mean summer air temperature reconstructed from a brGDGT-based proxy (blue; ±2σ) and relative 2010 data point in approximate relative position (purple; ±2σ). Charcoal counts reported as the number of fragments per volume (fragments cm⁻³) of peat (Orange ±2σ). Green boxes indicate relative depths of pollen sampling. Elevation of the deposit is reported as meters above sea level. (Data: Table S3.)

3.3 Vegetation and fire reconstruction

All sediment samples from BP contained charcoal (Fig. 4), indicating the consistent prevalence of biomass burning in the High Arctic during this time period. However, counts were variable throughout the section, with the middle and lower sections (mean 34 fragments cm⁻³) containing less charcoal compared to the upper section (mean 444 fragments cm⁻³). Overall, samples from BP contained on average 100.0±165 fragments cm⁻³ (mean ± 1σ), with charcoal area averaging 12.3±20.2 mm² cm⁻³. The variability of charcoal within any given sample was relatively low with a 1σ among charcoal area of approximately 2 mm² cm⁻³.

The three parts of the section analyzed for pollen (380.3–380.4, 381.10–381.25 and 381.35–381.45 m a.s.l.) reveal variations in vegetation (Figs. 4 and 5). Near the bottom of the section (380.3–380.4 m a.s.l.), Larix (26 %) and Betula (17 %) were the dominant trees. Alnus (6 %) and Salix (6 %) together with ericaceous pollen (4 %) were relatively high. In contrast, low numbers of Picea (3 %), Pinus (3 %) and fern spores were recorded. Additional wetland taxa like Myrica (5%) and Cyperaceae (6 %) were also noted. Overall, the non-arboreal (23 %) signal was well developed. Crumpled and/or ruptured inaperturate grains with surface sculpturing that varied from scabrate to verrucate were noted in the assemblage (12 %) but could not be definitively identified. It is possible that these grains represent Populus, Cupressaceae or additional Cyperaceae pollen. Between 381.10 and 381.25 m a.s.l., Larix (38 %) and Betula (21 %) increased in abundance, followed by ferns (7 %). Cyperaceae remained at similar levels (6 %), whereas Picea and Pinus decreased to 2 % and 1 %, respectively. Unidentified inaperturate types collectively averaged 14 %. Larix pollen (23 %) remained abundant near the top of the section (381.35–381.45 m a.s.l.), whereas Betula (2 %) decreased. Picea (16 %) Pinus (6 %) and ferns (23 %) increased in abundance. Of the ferns, trilete spores and cf. Botrychium were most abundant, followed by cf. Dryopteris. Inaperturate unknowns (10 %) were also observed. Other notables included Ericaceae (2 %) and Cyperaceae (2 %). While rare, Onagraceae grains were also observed (Fig. 5).

Figure 5 (a) Bar charts showing the relative pollen abundance in each portion of the section (error bars: 95 % confidence intervals; m a.s.l. – meters above sea level). (b) Pollen plate of select grains encountered in the BP section: (1) Pinus, (2) half a Picea grain, (3) Larix, (4) Betula, (5) Alnus, (6) Salix, (7) Myrica, (8) ericaceous grain, (9) Epilobium and (10) Cyperaceae; 50 µm scale: (1)–(3); 75 µm scale: (4)–(10).

According to the GBIF-based mapping exercise, the paleo-floral assemblage at BP most closely resembles modern vegetation found in northern North America, particularly on the eastern margin (e.g., New Hampshire, New Brunswick and Nova Scotia) and the western margin (Alaska, Washington, British Columbia and Alberta; Fig. 7a), and central Fennoscandia. Of these areas, the western coast of northern North America and eastern coast of southern Sweden have the most similarity to the reconstructed BP climate in terms of MST (Fig. 7b) and summer precipitation (Fig. 7c).

Figure 6 Examples of the feedbacks between temperature, vegetation and wildfire at the Beaver Pond site.

Figure 7 (a) Modern geographic distribution of observed occurrences of species common to the Beaver Pond species list. (b) Mean temperature of the warmest quarter (summer average) derived from WorldClim. (c) Mean precipitation of the warmest quarter (summer rain) derived from WorldClim. (d) Count of unique fire pixels detected per day over 10 years from MODIS 6 fire product, normalized by area of the latitude by longitude grid.

While high counts of active fire days are common in the western part of the North American boreal forest, they are not as common in the eastern part of the North American boreal forest (Fig. 7d), likely due to the differences in the precipitation regime. Low fire counts also typified Fennoscandia, likely due to historical severe fire suppression (Brown and Giesecke, 2014; Niklasson and Granström, 2004). Therefore, based on our reconstruction of the climate and ecology of the BP site, our results suggest that BP most closely resembled a boreal-type forest ecosystem shaped by fire, similar to those of Washington, British Columbia, the Northwest Territories, Yukon and Alaska (Fig. 7; but see Sect. 4.3).

4 Discussion

4.1 Geochronology

The plant and animal fossil assemblages observed at BP suggest a depositional age between 3 and 5 Ma (Matthews Jr. and Ovenden, 1990; Tedford and Harington, 2003). This biostratigraphic age was corroborated with an amino-acid racemization age ($>\mathrm{2.4}\pm \mathrm{0.5}$ Ma) and Sr-correlation age (2.8–5.1 Ma) on shells (Brigham-Grette and Carter, 1992) in biostratigraphically correlated sediments on Meighen Island, situated 375 km to the west–northwest. The previously calculated burial age of 3.4 Ma for the BP site is a minimum age because no postdepositional production of ²⁶Al or ¹⁰Be by muons was assumed. If the samples are considered to have been buried at only the current depth (ca. 10 m; see supplemental data), then the ages plot to the left and outside of the burial field, indicating that the burial depth was significantly deeper for most of the postdepositional history. The revised cosmogenic nuclide burial age is $\mathrm{3.9}+\mathrm{1.5}/-\mathrm{0.5}$ Ma. It is the best interpretation of burial age data based on improved production rate systematics (e.g., Lifton et al., 2014) and more reasonable estimates of erosion rate and ice cover since the mid-Pliocene (see Fig. S3; Table S5). As the stratigraphic position of the cosmogenic samples is very close to the BP peat layers, we interpret the age as representing the approximate time that the peat was deposited.

4.2 Fire, vegetation, temperature: a feedback triangle

Wildfire is a key driver of ecological processes in modern boreal forests (Flannigan et al., 2009; Ryan, 2002) and, although historically rare, has been becoming more frequent in the tundra in recent years (Mack et al., 2011). The modern increase in fire frequency is likely a consequence of atmospheric-CO₂-driven climate warming and feedbacks such as reduced sea ice extent (Hu et al., 2010) because the probability of fire is highest where temperature and moisture are conducive to growth and drying of fuels followed by conditions that favor ignition (Whitman et al., 2015). Young et al. (2017) confirmed the importance of summer warmth and moisture availability patterns in predicting fire across Alaska, highlighting a July temperature of ∼13.5 ^∘C as a key threshold for fire across Alaska.

The abundance of charcoal at BP demonstrates that climatic conditions were conducive to ignition and that sufficient biomass available for combustion existed across the landscape. brGDGTs-derived temperature estimates suggest mean summer temperatures at BP exceeded the ∼13.5 ^∘C threshold that drastically increases the chance of wildfire (Young et al., 2017). Indeed, the estimate of ∼15.4 ^∘C suggests mean reconstructed summer temperatures were ∼11 ^∘C higher than modern-day Eureka, Canada (∼4.1 ^∘C; Fig. 2), representing substantial additional amplification compared to the global average. Without the increased Arctic amplification of temperature that accompanies climate equilibrium with high CO₂, mean summer temperatures would be lower than the July temperature threshold that predicts increased wildfire. This is evidence that Pliocene Arctic amplification of temperatures was a direct feedback to increased wildfire activity. The increased extension of boreal forest into the Arctic was also possible due to Arctic amplification of temperatures. This biomass provided the fuel for combustion, and thus Pliocene Arctic amplification of temperatures is also an indirect feedback to wildfire (Fig. 6).

Conversely, the charcoal record at BP suggests that substantial biomass burning during the mid-Pliocene could have acted as a feedback mechanism amplifying or dampening seasonal warming (Fig. 6). Studies of the impact of wildfire on surface energy balance in present-day northern ecosystems have revealed the complexity of predicting wildfire's impact on climate. Ecosystems exhibit changing responses through time from the scale of years post-burning (Randerson et al., 2006; Bonan, 2008; French et al., 2016) to seasonal (Huang et al., 2014) and even diurnal differences post-deforestation that may impact net wildfire feedback to climate (Schultz et al., 2017). The radiative response to wildfire changes across latitudinal gradients (Jin et al., 2012) and between local and global scales (Ward et al., 2012; Liu et al., 2019). Additionally, the original vegetation type burned influences aspects of wildfire's impact on climate such as the original albedo (French et al., 2016), likely fire severity and intensity (Rogers et al., 2015) and time to pre-fire ecosystem recovery (French et al., 2016) or alternate ecosystem establishment (Johnstone et al., 2010b). The mechanisms that appear to have the largest effect include carbon release and sequestration (e.g., Harrison et al., 2018), changes in surface albedo (e.g., Huang et al., 2014), altered evapotranspiration (Liu et al., 2019), and aerosol effects both directly and also indirectly via cloud processes (e.g., Stone et al., 2008; Zhang et al., 2017). The potential role of wildfire as a feedback to climate in the mid-Pliocene Arctic is suggested by its prevalence through this > 20 000-year sequence, the impact of forest fire in modern ecosystems, and preliminary modeling of the complex direct impacts on the surface radiative budget (e.g., short-term black carbon deposition on snow and ice and long-term changes in albedo) and direct and indirect effects at the top of the atmosphere radiative budget (i.e., aerosol emissions; Feng et al., 2016). Further modeling experiments are needed to determine whether wildfire played a significant role in the magnitude and seasonal patterns of mid-Pliocene Arctic amplification of temperature.

An increase in atmospheric convection has been simulated in response to diminished sea ice during warmer intervals (Abbot and Tziperman, 2008), but this study did not confirm whether this increase in atmospheric convection was sufficient to cause lightning ignitions. An alternative ignition source for combustion of biomass on Ellesmere Island during the Pliocene is coal seam fires, which have been documented to be burning at this time (Estrada et al., 2009). However, given the interaction of summer warmth and ignition by lightning within the same climate range as posited for BP, we consider lightning the most likely source of ignition for Pliocene fires in the High Arctic.

Fire return intervals cannot be calculated from the BP charcoal counts due to the absence of a satisfactory age–depth model and discontinuous sampling. As strong interactions are observed between fire regime and ecosystem assemblage in the boreal forest (Brown and Giesecke, 2014; Kasischke and Turetsky, 2006) and in response to climate, comparison with modern fire regimes for areas with shared species compositions and climates may inform a potential range of mean fire return intervals (MFRIs).

Matthews and Fyles (2000) indicated that the Pliocene BP environment was characterized by an open larch-dominated forest–tundra environment. The modern area with the most species in common with BP is central northern Alaska (Fig. 7a). The area over which shared species were calculated is largely tundra but includes the ecotone between tundra and boreal forest. Other zones that share many species with BP are continuous with Alaska down the western coast of North America to the region around the border of Canada and the United States, the eastern coast of North America in the region around the border of Canada and the United States (∼50^∘ N), and central Fennoscandia. Of these zones, the MSTs of Alaskan tundra sites (6–9 ^∘C) are less similar to BP (15.4 ^∘C) than ∼50^∘ N on both western and eastern coastal North American sites and central Fennoscandia (12–18 ^∘C, Fig. 7b). The eastern coast of North America has higher rainfall during the summer (>270 mm) than the west coast and Alaska (Fig. 7c), which correlates to the timing of western fires. The low summer precipitation for much of the west (<200 mm) is consistent with previously published summer precipitation estimates for BP (∼190 mm). As a result, the fire regime of the west coast ∼50^∘ N may be a better analogue for BP than the east coast of North America. In central Fennoscandia there is also a west vs. east coastal variation in summer precipitation, with the western, Nordic part of the region experiencing higher summer precipitation (252–>288 mm) than the more similar eastern, Swedish part of the region (∼198 mm).

Investigation of the modern fire detection data (Fig. 7d) suggests that the two regions most climatically similar to BP – ∼50^∘ N western North America and central Sweden – have radically different fire regimes. It is likely this is caused by historical fire suppression in Sweden that limits the utility of modern data for comparison with this study (Brown and Giesecke, 2014; Niklasson and Granström, 2004). To understand the fire regimes, as shaped by climate and species composition rather than human impacts, we considered both the modern and recent Holocene reconstructions for these regions (Table 1). This shows that (a) within any region variation arises from the complex spatial patterning of fire across landscapes and (b) that the regions most similar to BP (∼50^∘ N western North American and eastern Fennoscandian reconstructions for the recent Holocene) have shorter fire return intervals than the cooler Alaskan tundra or wetter summer ∼50^∘ N region of the eastern North American coast.

Table 1 Modern and recent Holocene fire return interval reconstructions for the candidate analogous regions considered in this study.

^a The reciprocal converted from burn rate (%) (see Van Wagner et al., 2006). ^b Estimates likely effected in some areas by human activity. In such instances Recent Holocene is preferred. ^c Fire cycle. ^d “Recent” here refers to records (distinct sections) that begin after the end of the Holocene Climate Optimum and end near the present.

While the shared species for Siberia appear low, the total number of observations for Siberia in the modern biodiversity database used is likewise low – and the latter is a potential cause of the former. Given the similar climate to BP on the Central Siberian Plateau and some key aspects of the floras in Siberia such as the dominance of larch, we considered the fire regime of the larch forests of Siberia. Kharuk et al. (2016, 2011) studied MFRIs across Siberia, from 64 to 71^∘ N, the northern limit of larch stands. They found an average MFRI across that range of 110 years, with MFRI increasing from 80 years in the southern latitudes to ∼300 in the north (Table 1). Based on the similarity of the climate variables, the more southerly MFRIs (∼80 years) may be a better analogue. Key differences between boreal fires in North America compared to Russia are a higher fire frequency with more burned area in Russia but a much lower crown fire and a difference in timing of disturbance, with spring fires prevailing in Russia compared to mid-summer fires in western Canada (de Groot et al., 2013; Rogers et al., 2015).

The pollen-based vegetation reconstruction derived in this study indicates that open Larix–Betula parkland persisted in the basal (380.3–380.4 m a.s.l.) parts of the sequence. Ground cover was additionally dominated by shrub birch, ericaceous heath and ferns. While the regional climate may have been somewhat dry, the record suggests that, locally, a moist fen environment dominated by Cyperaceae existed near the sampling location. Shrubs including Alnus and Salix likely occupied the wetland margins.

The corresponding relatively low concentration of charcoal in this stratigraphic interval may reflect lower-severity fires or higher sedimentation rates. We consider the former more likely due to the depositional environment of Unit III from Mitchell et al. (2016), a lake edge fen peat in a beaver pond or small lake, without evidence of high sediment influx overwhelming peat production. We posit that a surface fire regime, somewhat like that in southern central Siberia existed. This premise is also supported by the fire ecology characteristics of the dominant vegetation. Larix does not support crown fires due to leaf moisture content (de Groot et al., 2013) and self-pruning (Kobayashi et al., 2007). The persistence and success of larch in modern-day Siberia appears to be driven by its high growth rate (Jacquelyn et al., 2017), tolerance of frequent surface fire due to thick lower bark (Kobayashi et al., 2007) and tolerance of spring drought due to its deciduous habit (Berg and Chapin III, 1994). Arboreal Betula are very intolerant of fire and easily girdled. However, they are quick to resprout and are often found in areas with short fire return intervals. Like Larix, arboreal Betula have high moisture content of their foliage and are not prone to crown fires. Betula nana L., an extant dwarf birch, is a fire endurer that resprouts from underground rhizomes or roots (Racine et al., 1987) thus regenerating quickly following lower-severity fires (de Groot et al., 1997). The vegetation and fire regime characteristics are similar further up the sequence at 381.10–381.25 m a.s.l., with the exception that ferns increased in abundance while heath decreased.

In the upper part of the sequence (381.35–381.45 m a.s.l.), where charcoal was abundant, the Larix–Betula parkland was replaced by a mixed boreal forest assemblage with a fern understory. Canopy cover was more closed compared to the preceding intervals. The forest was dominated by Larix and Picea, with lesser amounts of Pinus. While Betula remained part of the forest, it decreased in abundance possibly due to increased competition with the conifers. Based on exploratory CRACLE analyses of climate preferences using GBIF occurrence data (GBIF.org, 2018a, b, c, d) of the dominant taxa (Larix–Betula vs. Larix–Picea–Pinus), the expansion of conifers could indicate slightly warmer summers (MST∼15.8 ^∘C vs. 17.1 ^∘C). This result differs from the stable MST estimated by bacterial tetraethers, although within reported error, and the small change is certainly within the climate distributions of both communities. The CRACLE analyses also suggest that slightly drier conditions may have prevailed during the three wettest months (249–285 mm vs. 192–219 mm). While the interaction between climate, vegetation and fire is complex, small changes in MST and precipitation could have directly altered both the vegetation and fire regime, which in turn further promoted fire-adapted taxa. In addition to regional climatic factors, community change at the site may have been further influenced by local hydrological conditions, such as channel migration, pond infilling and ecosystem engineering by beaver (Dipoides sp.).

The high charcoal content of the upper portion (∼ Unit IV) of the sequence has three potential explanations: reworking of previously deposited charcoal, decreased sedimentation or increased wildfire production of charcoal. The first explanation seems unlikely because there is no difference in the shape of the macro-charcoal between the upper and lower portions of the sequence. A change in the dimensions of the charcoal would be expected if it had undergone additional physical breakdown from reworking (see Fig. S4). The second, decreased sedimentation, may occur if the deposition is a result of infrequent, episodic flooding intermixed with long periods during which charcoal was deposited. The recorded sedimentology does not support this explanation but, due to the complexity of flooding processes, also does not disprove this explanation. The third explanation, that increased charcoal reflects increased wildfire, is supported by the change in plant composition and suggests that frequent, mixed-severity fires may have persisted at this time. While Larix is associated with surface fire, Picea and Pinus are adapted to higher-intensity crown fires. A crown fire regime may have established as conifers expanded, altering fuel loads and flammability. For example, black spruce sheds highly flammable needles, its lower branches can act as fuel ladders facilitating crown fires (Kasischke et al., 2008), and black spruce was previously tentatively identified at BP (Fletcher et al., 2017). While it has thin bark and shallow roots maladapted to survive fire (Auclair, 1985; Brown, 2008; Kasischke et al., 2008), it releases large numbers of seeds from semi-serotinous cones, leading to rapid reestablishment (Côté et al., 2003). The documentation of Onagraceae pollen at the top of the sequence could potentially reflect post-fire succession. For example, the species Epilobium angustifolium L. is an early seral colonizer of disturbed (i.e., burned) sites, pollinated by insects.

It appears that the Larix–Betula parkland dominated intervals correspond to the peat- and sand-stratigraphic Units II and III described by Mitchell et al. (2016), whereas the mixed boreal forest in the upper part of the sequence is contemporaneous with Unit IV, described as peat and peaty sand, coarsening upwards. Thus, while vegetation and fire regimes seemingly changed through time at this Arctic site, temperatures appear more stable or at least to have no apparent trend within analytical and reconstruction uncertainty. Thus, it is suggested that the fire regime at BP was primarily regulated by regional climate and vegetation and perhaps additionally by changing local hydrological conditions. Regarding climate, MST remained high enough ($>\sim \mathrm{13.5}$ ^∘C) throughout the sequence to allow for fire disturbance, and the pollen suggests that temperatures may have marginally increased in the upper part of the sequence. Alternatively, other climate variables, such as the precipitation regime, or local hydrological change may have initiated the change in vegetation community. Up-sequence changes in vegetation undoubtedly influenced fine fuel loads (e.g., surface layer needles, mosses and twigs) and flammability. Indeed, the fire ecological characteristics of the vegetation are consistent with a regional surface fire regime yielding to a crown fire regime.

Betula and Alnus, which occurred earlier in the depositional sequence, are favored by beaver in foraging (Busher, 1996; Haarberg and Rosell, 2006; Jenkins, 1979). Moreover, the presence of sticks cut by beaver in Unit III reveals that beavers were indeed at the site, moistening the local land surface. The lack of beaver-cut sticks and changes in sediment in Unit IV may indicate that the beavers abandoned the site, possibly in response to changes in vegetation (i.e., increased conifers and decreased Betula) limiting preferred forage or due to lateral channel migration, as evidenced by the coarsening upward sequence described by Mitchell et al. (2016). As a result, the local land surface may have become somewhat drier, contemporaneous with the change towards Larix–Picea–Pinus forest and a mixed-severity fire regime.

5 Conclusion

The novel temperature estimates presented here confirm that Ellesmere Island summer temperatures were considerably warmer (15.4±0.8 ^∘C) during the likely > 20 000-year mid-Pliocene interval ($\mathrm{3.9}+\mathrm{1.5}/-\mathrm{0.5}$ Ma) investigated, compared to the modern Arctic. The ∼11 ^∘C higher than present-day summer temperatures at Beaver Pond support an increasing effect of Arctic amplification of temperatures when CO₂ reaches and exceeds modern levels. Our reconstruction of the paleo-vegetation and ecology of this unique site on Ellesmere Island suggests an assemblage similar to forests of the western margins of North America and eastern Fennoscandia. The evidence of recurrent fire and concurrent changes in taxonomic composition are indicators that fire played an active role in mid-Pliocene Arctic forests, shaping the environment as it does in the boreal forest today. Evidence from fire in the modern boreal forest suggests that fire may have had direct and indirect impacts on Earth's radiative budget at high latitudes during the Pliocene, acting as a feedback to Pliocene climate. The net impact of the component process remains unknown and modeling experiments are needed to quantitatively investigate the effects of the kind of fire regime presented here, on the Pliocene High Arctic. Collectively, these reconstructions provide new insights into the paleoclimatology and paleoecology of the Canadian High Arctic, ∼3.9 Ma.