Does fire always accelerate shrub expansion in Arctic tundra? Examining a novel grass-dominated successional trajectory on the Seward Peninsula

ABSTRACT

Over the last century in the circumpolar north, notable terrestrial ecosystem changes include shrub expansion and an intensifying wildfire regime. Shrub invasion into tundra may be further accelerated by wildfire disturbance, which creates opportunities for establishment where recruitment is otherwise rare. The Seward Peninsula currently experiences more frequent and larger fires than other tundra regions in Alaska. There are areas of overlapping burn scars dating back to the 1950s. Using a chronosequence approach, we examined vegetation and ecosystem dynamics in tussock tundra. Increasing burn severity and fire frequency corresponded with an increase in grass cover and a decrease in shrub basal area. We used multivariate ordination analysis to create a single integrator variable of fire effect that accounted for time after fire, burn severity, and number of times burned. This fire effect was significantly associated with decreases in soil organic layer thickness and overall plant biomass. Unlike previous studies in Arctic Alaska tundra, we found that increases in fire frequency and severity did not increase shrub cover and biomass. Instead, intensifying fire disturbance, and particularly repeat fires, led to grass dominance. Our findings support the hypothesis that intensifying tundra fire regimes initiate alternative post-fire trajectories that are not shrub dominated and that are structurally and functionally quite different from sedge or shrub-dominated tundra.

KEYWORDS:

• Alaska
• nonmetric multidimensional scaling
• principal component analysis
• linear mixing model
• historical fire regime
• climate change
• repeat burning

Introduction

Climate change is affecting disturbance regimes in many regions of the world (Dale et al. 2001; Turner 2010). In high-latitude systems, a large body of research focuses on climate-induced changes in the fire regime of the boreal forest (e.g., Bergeron and Flannigan 1995; Soja et al. 2006; Chapin et al. 2010). In arctic tundra, however, wildfire has been less frequent or has occurred in remote localities where small short-lived fires or their fire scars were not observed and are therefore not part of the historical record (D. A. Walker and Walker 1991). The warming the Arctic has experienced over the last fifty years (Hinzman et al. 2005) is significant and coincident with an increase in the frequency and severity of wildfire (French et al. 2015; Masrur, Pretrov, and De Groote 2018). For example, the recent 2007 Anaktuvuk River fire was the largest and longest-burning wildfire known to occur on the North Slope of Alaska (Jones et al. 2015) in recorded history. Not only did the majority of the fire burn with moderate to high severity (Jones et al. 2009) but there also was no evidence of fire at this site for at least the previous 5,000 years (Hu et al. 2010; Chipman et al. 2015). Even when rare, tundra fires play an important role in releasing the carbon stored in tundra vegetation and soils to the atmosphere (Mack et al. 2011). Due to large carbon stores in tundra soils and their vulnerability to burning, increasing carbon emissions from tundra fires are predicted to feedback positively to climate warming (Genet et al. 2017).

Climate change is likewise expected to increase the extent and frequency of wildfire throughout Alaska (Melvin et al. [2017] and references within). The recent trend toward an increase in tundra fire frequency corresponds to high growing season temperature and low precipitation from 1979 to 2009 (Hu et al. 2010; Moritz et al. 2012). For example, across Alaska, when mean July temperatures exceed a threshold of 13.4°C and annual moisture availability falls short of 150 mm, there is a greater likelihood of fire occurrence (Young et al. 2017). In response to projected increases in temperature, modeling results show an increase in tundra area burning for northwestern and northern Arctic Alaska over the next century (Rupp et al. 2016), largely due to the positive effect of lower fuel moisture content on ignition and spread potential.

Over the last century, circumpolar changes in vegetation have included the expansion of shrubs and northward shifts of latitudinal treeline (Lloyd et al. 2003; Tape, Sturm, and Racine 2006; Myers-Smith, Hik, and Aerts 2017). The invasion into tundra by woody plants may be further facilitated by wildfire disturbance, which creates opportunities for establishment where recruitment is otherwise rare (Landhausser and Wein 1993; Joly, Chapin, and Klein 2010; Hewitt et al. 2017). Currently, fires in Arctic Alaska are limited by their ignition source (infrequent lightning) and the availability of favorable burn weather. However, the paleorecord shows that there is a strong historical link between fire regime and tundra vegetation in Alaska. The last time the Earth experienced dramatic climate warming, during the early Holocene (14,000–10,000 BP), wildfire occurrence corresponded with the presence and abundance of highly flammable birch (Betula) shrubs on the landscape (Higuera et al. 2008). These historical and contemporary patterns suggest positive feedbacks between increased shrub cover and increased tundra burning (Higuera et al. 2011). Changes in summer weather and the composition and abundance of fuels in tundra could contribute to further shifts in the fire regime that will directly impact post-fire successional trajectories. For example, in northwest Alaska, bryophytes, sedges, and graminoids—from both regrowth and seedling recruitment—dominated severely burned tussock–shrub and dwarf shrub tundra along a hillslope during the first decade postrecovery (Racine 1981). Then in the second and third decades, vegetation shifted to shrub dominance, mainly Rhododendron tomentosum spp. decumbens and the willow Salix pulchra, such that shrub cover was greater than before the fire (Racine et al. 2004). However, long-term post-fire tundra successional trajectories are relatively unknown beyond a few point localities, so it is unclear whether this observed trend is expected to be common or an anomaly (but see Racine, Dennis, and Patterson 1985; Racine, Johnson, and Viereck 1987). Understanding the consequences of a change in fire regime is complicated by the fact that there are relatively few large recent fires that are road accessible for easy study and monitoring post-fire over time.

Burn severity, or the change in plant biomass and soil organic material due to burning, can influence patterns of post-fire regeneration and reaccumulation of fuels by affecting mortality rates, post-fire seed and bud availability, and physical soil properties that affect plant establishment and growth (Johnstone et al. 2010; Hollingsworth et al. 2013). The magnitude of duff consumption affects both the availability of seed and bud banks (Schimmel and Granstrom 1996) and the physical properties of the post-fire soil (Dyrness and Norum 1983) that affect successful germination and establishment. In tundra ecosystems that accumulate thick duff layers (>50 cm), the degree to which these surface moss and organic soil layers are consumed during combustion is likely the most important control on plant regeneration.

Although relatively unstudied in tundra, evidence from the neighboring boreal forest suggests that burn characteristics affect patterns of plant regeneration during this critical establishment period by exerting a strong influence on future vegetation characteristics and plant succession in two major ways. First, large-scale disturbances that expose mineral soil in tundra promote increased recruitment of shrub species (Lantz et al. 2009). Second, if post-fire tundra plant communities become dominated by shrubs, this will alter the ecosystem structure and increase the accumulation rate of aboveground fuels and duff and potentially future flammability.

Understanding the post-fire successional trajectories and ecosystem consequences of a changing tundra fire regime is critical to forecasting impacts on regional carbon balance, permafrost degradation, and shifts in species composition. Here we quantified vegetation composition and ecosystem characteristics (i.e., soil carbon [C] and nitrogen [N] stocks, aboveground biomass, shrub density, and foliar C and N) from multiple sites across the Seward Peninsula. Compared to elsewhere in Arctic Alaska, the Seward Peninsula in northwestern Alaska has experienced the greatest number of fires, cumulative area burned, and area reburned (Racine, Johnson, and Viereck 1987). In this particular tundra region, a mosaic of overlapping burns has been mapped since the 1950s (Alaska Interagency Coordination Center, unpublished data).

Our specific research questions in this study were as follows: (1) How do time after fire, burn severity, and times burned (collectively called “metrics to characterize fire regime”) affect the composition of species and plant functional types (e.g., sedges, grasses, shrubs)? We hypothesized that an intensifying fire regime—increases in burn severity and a shortened fire return interval—would favor shrub establishment and growth. To evaluate this objective, we described and compared species composition and plant functional types across plots along a gradient in fire history and regime. To better understand the underlying mechanisms driving post-fire successional pathways we asked: (2) How do metrics to characterize fire regime affect ecosystem function? We hypothesized that these metrics—time after fire (fire history), burn severity (fire history and regime), and times burned (fire regime)—would be coupled with an increasing shrub:tussock ratio, indicative of increasing shrub abundance in graminoid-dominated tussock tundra, leading to variation in nutrient cycling and carbon storage. To evaluate this objective, we correlated metrics to characterize fire regime, using a single synthetic fire effect variable, with above- and belowground ecosystem properties (specifically, foliar C and N, shrub:tussock characteristics, vegetation characteristics, soil organic depth, thaw depth, and soil C and N). Collectively, we used a combination of descriptive and correlation multivariate statistical approaches to disentangle the role of fire in shifting successional trajectories and present a conceptual model relying on results from both analyses and the literature to suggest mechanisms associated with variation in successional trajectories in western Alaska arctic tundra.

Methods

Study area

We studied tundra fire scars in two regions on the Seward Peninsula (Figure 1, Table 1). Our main study region was north of Nome near Quartz Creek and the Kougarok River, at the end of the Kougarok Road (Seward Peninsula West; Figure 1A). The Kougarok tundra fire field site has burned in five major fires in the decades since 1950 (Alaska Interagency Coordination Center, unpublished data). Lightning ignited four of these fires at the field site and one was human caused. The mosaic of overlapping fire scars allowed for the study of repeat fires. We also surveyed adjacent tundra with no record of burning west and south of the Kougarok tundra fire site. Our secondary study region was south of Kotzebue near the Buckland River (Seward Peninsula East; Figure 1B). Both regions have similar landscape characteristics and climate conditions (Lloyd et al. 2003), though the Kougarok tundra fire site has been more extensively studied due to its proximity to the road system (Liljedahl et al. 2007; Narita et al. 2015; Iwahana et al. 2016; Tsuyuzaki, Iwahana, and Saito 2017).

Table 1. Number of plots, fire scars sampled, and the years each plot burned at our two study regions on the Seward Peninsula, Alaska

Abbreviated sampling locality	Fire scars sampled (times burned)	Year(s) burned	Number of plots
Seward Peninsula West
Delome_71	Delome River (1)	1971	10
MP85_2(x2)	Delome River and Milepost 85 (2)	1971, 2002	5
MP85_2(X3)	Delome River, Garfield Creek, and Milepost 85 (3)	1971, 1997, 2002	5
Unburned	No evidence of recent burn	n/a	8
Seward Peninsula East
FishRr_11	Fish River (1)	2011	10
FishRr_Unk	Unknown fire (assumed 1)	Unknown; estimated 1982	5

Figure 1. Map of study sites on the Seward Peninsula, Alaska. Plots and fire scars are shown at (A) the end of the Kougarok road (Seward Peninsula West) and (B) Buckland River (Seward Peninsula East). The Minguk River fire scar (2010s), which burned in 2016 (Seward Peninsula West), is shown on the map, though sites were visited prior to this fire.

The central Seward Peninsula is characterized by continuous permafrost with a thickness of 15 to 30 m and a mean active layer thickness of 56 cm (Hinzman et al. 2003). Sloping hills with mixed shrub–tussock tundra and tussock tundra vegetation in the uplands are characteristic of the region. Three micrometeorological towers near the Kougarok field site recorded a mean annual temperature of −2.4°C, mean January temperature of −23.1°C, mean July temperature of +11°C, and mean summer rainfall (June–August) of 94 mm from 2000 to 2006 (Liljedahl et al. 2007).

Field methods

Vegetation sampling

To characterize vegetation composition and biomass, we sampled 1 m² plots within fire scars in both Seward Peninsula study regions and adjacent unburned tundra near the Kougarok tundra fire field site in July 2012. Plots were arranged within each fire complex—area of fires, with multiple fires—so that we sampled vegetation that burned once, twice, or thrice and that had no evidence of recent fire for a total of forty-three vegetation plots (Table 1). Our sample size was limited by the difficulty of traversing the tundra landscape by foot, even though one of our study regions is road accessible. At each plot the following eight measurements were made: (1) full vascular and nonvascular species lists and their percentage cover (plant nomenclature follows the Arctic Vegetation Archive’s Pan-Arctic Species List; D. A. Walker et al. 2018), (2) height and absolute cover of plant functional types (shrubs [deciduous shrubs and evergreen shrubs], graminoids [grasses and sedges], forbs, bryophytes, and lichens), standing dead and woody debris, litter, bare soil, and standing water; (3) shrub basal diameter for all shrubs with >1 cm basal diameter; (4) aboveground biomass clipped in a 25-cm subplot for all vascular <1 cm basal diameter and nonvascular plants; (5) number of tussocks and their diameter for Eriophorum vaginatum, the sedge that is the foundation of the sampled tundra communities. Tussock diameters were measured twice in perpendicular directions across each tussock and then averaged; (6) foliar C and N content from two common upland shrubs (Salix spp., either S. pulchra or S. glauca if S. pulchra was not present and Betula nana); and the tussock forming sedge E. vaginatum; (7) abiotic factors including latitude and longitude coordinates, topographic position, elevation, parent material, slope (%) measured with an inclinometer, and scalar site moisture (Johnstone, Hollingsworth, and Chapin 2008); and (8) mid-season synchronous thaw depth, soil organic layer (SOL) depth, and mineral soil pH at the four corners of each 1 m² plot (see below).

Soil sampling

Thaw depth was measured by inserting a metal rod into the soil until it hit ice, marking the surface of the green moss on the rod, removing it, and measuring the distance from the tip to the mark with a meter stick. In some cases rock was encountered above the permafrost table, preventing an accurate measurement of thaw depth and this was recorded. SOL depth was measured by slicing a square pit with a serrated knife, removing a monolith of organic soil, exposing the surface of the mineral soil, and measuring the distance from the surface of the green moss to the mineral soil on two sides of the pit (averaged to yield one SOL depth measurement/point). If ice was encountered before mineral soil, this was reflected to account for incomplete measurements of SOL depth. Mineral soil pH was measured in situ with an Oakton waterproof pHTestr 30, and a soil sample was collected for determination of bulk density and C and N concentrations.

At one randomly selected corner in each plot, we harvested a second 10 × 10 cm monolith of organic soil from the exposed side of the pit. The monolith extended from the surface of the green moss to the surface of the mineral soil (roughly 5–30 cm, depending on location). We clipped aboveground vegetation from the surface of the monolith, measured its dimensions, wrapped it in multiple layers of tinfoil to preserve structure, placed it in a cooler with blue ice, and returned it to the lab at the University of Alaska Fairbanks, where these samples were frozen and shipped to the University of Florida for additional soil analyses (C and N). We also collected a 10 cm deep × 5.5 cm diameter mineral soil core from the pit for analyses of soil pH to compare with in situ measurements of soil pH.

Metrics used to characterize fire regime

We surveyed various burn metrics in the field for each plot. Substrate and vegetation burn severity composite scores for each plot were collected using the Alaska Interagency Fire Effects Task Group protocol (Alaska Interagency Fire Effects Task Group 2007), optimized for the tundra region. Composite scores vary from low to moderate to high, using a suite of measurements and indicators in the field (Supplementary Table 1). This is a relative scalar that depends on time after fire, so we use it to interpret evidence of fire. For example, the post-fire appearance of tussocks indicates a burn severity of “low” if tussocks have moderate quantities of senesced leaf litter and less than 20 percent of the basal area consumed and a severity of “high” if tussocks have no leaf litter remaining and greater than 60 percent of their basal area consumed. This protocol, however, makes it difficult to interpret burn severity with increasing time after fire.

Fire history reflects the fire regime, and we often derive metrics from fire histories for the purpose of characterizing fire regime. For our analyses, we later grouped plots according to quantitative, categorical, and scalar metrics to characterize fire regime: (1) time after fire (fire history characteristic): number of years since a given plot burned (or <30 years [N = 23] and ≥31 years classes [N = 20]); (2) times burned (fire frequency or a fire regime characteristic): number of times a given plot burned (none [N = 8], once [N = 25], or more than once [twice and thrice; N = 9]); and (3) burn severity (both a fire history and fire regime characteristic): composite score developed for this study from 0 to 8 for a given plot, related to both the substrate (mineral soil exposed and charred layer in soil) and vegetation burn severity (abundance of charred stems or culms) and whether there was sign of past fire on shrubs or tussock (low [0–1, N = 10], moderate [2–4; N = 17], or high [5–8; N = 16]. Putatively unburned plots (burn severity = 0) were grouped with those with the lowest burn severity (burn severity = 1), because there were only two burned plots classified as low severity; both occurred at the 1971 Delome fire scar, having burned once. There are no plots in our study that are both thirty-one years or more after fire and showing high burn severity, due to the methodological issue of decreasing ability to detect severity with time after fire.

Laboratory methods

Species composition

Voucher species for all vascular and nonvascular species were verified in the laboratory. All vascular species were verified at the University of Alaska Fairbanks Museum of the North Herbarium. Lichen and bryophyte specimens were sent to the Komarov Botanical Institute at the Russian Academy of Sciences and verified by Misha Zhurbenko and Olga Afonina, respectively.

Biomass estimates

All clipped aboveground biomass samples were dried at 25°C for 48 hours and weighed (dry [fresh] weight; hereafter called plot biomass). We created allometric equations (Berner et al. 2015) to estimate shrub biomass for shrubs >1 cm² in basal diameter by collecting multiple shrubs of various species and basal diameters outside of plots and weighing them (hereafter called shrub biomass). We used a published allometric equation to estimate tussock biomass (Mack et al. 2011), using measurements of tussock diameters.

Soil analyses

A soil monolith was collected at each plot, comprised of 5- to 10-cm increment soil samples. In the lab, each monolith was thawed and sliced into 5-cm depth intervals with a serrated electric carving knife, with the last sample of variable depth depending on the location of the organic/mineral interface. Samples were homogenized by hand, and coarse organic materials (>2.5-cm twigs and roots) and rocks were removed. Coarse and fine organic fractions were weighed wet, dried at 70°C for 48 hours to determine dry matter content, and then ground on a Wiley mill with a 40-mm sieve. Carbon and N content were measured on a Costech Elemental Analyzer (EA, Costech Analytical, Los Angeles, CA), calibrated with the standard NIST peach leaves (SRM 1547, National Institute of Standards and Technology, Gaithersburg, MD). Rocks were weighed and the volume was measured by displacement in water. The volume of each monolith layer was calculated as depth multiplied by area minus the volume of rocks. Bulk density and C and N pools were calculated for both fine and coarse organic fractions. Pools were calculated for each depth increment by multiplying bulk density by concentration and depth and extrapolating to 1 m². Soil pH was measured on dried soil samples. Distilled deionized water was added to create a slurry and pH was measured using an Oakton pH 2700. Soil pH was calculated on a mean profile basis by multiplying the log-transformed pH in each depth increment by the proportional soil mass of that increment, summing the products, and back-transforming the sum.

Foliar C and N

Vegetation samples from E. vaginatum, B. nana, and Salix spp. (S. pulchra or S. glauca) from each plot were pooled and dried at 25°C for 48 hours and ground to a fine powder. Foliar C and N were measured on a Costech EA, as described above.

Statistical methods

All statistical analyses were performed using the R statistical package 3.3.2 (R Core Team 2016) unless otherwise specified.

Effect of fire on species composition and plant functional types

We evaluated diversity at three scales: alpha diversity (the average species richness per plot), beta diversity (species turnover, ratio of the total number of species across all sites to the average number of species per plot), and gamma diversity (regional diversity, or the total number of species across all plots; McCune and Grace 2002). We used analysis of variance (ANOVA) with Tukey’s post hoc tests to test for significant differences in alpha diversity between burn severity classes (defined in Metrics Used to Characterize Fire Regime) and an independent t-test to test for significant differences in time after fire classes (defined in Metrics Used to Characterize Fire Regime). Assumptions of ANOVA were met, including normality (normal probability plot) and homogeneity of variance (Levene’s test).

To examine the abundance of plant functional types in relationship to metrics used to characterize fire regime, we related the mean plot absolute cover by time after fire, burn severity, and times burned classes. Due to overlapping canopy layers, the sum of absolute cover can exceed >100 percent. We did not include rushes (Juncaceae) in the graminoid class or seedless vascular plants (horsetails and lycopods) in the resulting figures, because there were few taxa of these functional types within the sampled plots. We discarded one quadrat from all analyses because it was a multivariate outlier (Mahalanobis distance; McCune and Grace 2002). The discarded quadrat was composed of >70 percent forb (Petasites frigidus) cover, unlike any other sampled plot, because the maximum cover by this species elsewhere was <20 percent.

To quantify the differences in species composition between plots and correlate overall species composition to metrics to characterize fire regime and the local environment, we used nonmetric multidimensional scaling (NMDS; Mather 1976; Kruskal and Wish 1978), complemented with Kendall correlation coefficients (τ; PC-ORD v7.0; McCune and Mefford 2017). NMDS is an ordination technique that is well suited for species compositional data because it does not assume linear or orthogonal relationships between the axes (McCune and Grace 2002). Singletons were removed from the species composition data (i.e., species that only occurred at one plot). A preliminary analysis, using a Monte Carlo test of significance for a six-dimensional solution stepping down to a one-dimensional solution, was used to determine the appropriate dimensionality (number of axes) for ordination. The preliminary run included 500 iterations with fifty runs of real data, fifty runs of randomized data, and an instability criterion of 0.000001. Based on preliminary results, we chose a three-dimensional solution for the final analysis. All analyses used the Sorenson distance measure and random starting configurations. High Kendall correlations (r² > 0.25) between our main matrix (species vs. plot) and an environmental matrix (environmental variables vs. plot) were displayed using biplot overlays that indicate the strength and direction of the correlations with the environment. Finally, we calculated a post hoc proportion of variance represented by each axis (the coefficient of determination or r²) by calculating the Sorenson distance measure between pairwise distances of objects in ordination spaces and distances in the original matrix (McCune and Grace 2002).

To better visualize changes in vegetation in response to metrics to characterize fire regime in our ordinations, we post hoc grouped our vegetation plots into four distinct plant communities that primarily differ in percentage cover of graminoids vs. shrubs following the classification of Viereck et al. (1992). Plots with at least 25 percent cover of erect shrubs were classified as either shrub tundra (>75 percent shrub cover; N = 9) or mixed shrub–tussock tundra (N = 11). The remaining plots are all mesic graminoid herbaceous and classified as either tussock tundra (N = 14) or graminoid–forb tundra (N = 8). Graminoid–forb tundra, however, is not a vegetation type identified by Viereck et al. (1992) but is most similar to their bluejoint meadow, which is dominated by Calamagrostis canadensis and is described as often being pure stands or having other grasses and forbs present but not codominant. Bluejoint meadows are most common south of latitudinal treeline where their development is generally initiated by disturbance such as fire and land clearing. Bluejoint communities on the Seward Peninsula have been reported and described as small and largely restricted to disturbed sites such as villages and recently drained lake basins (Racine and Anderson 1979). We chose to designate these plots as graminoid–forb tundra because they differ in being codominant by three different grasses (Calamagrostis lapponica, Calamagrostis neglecta, and Arctagrostis latifolia), none of which are the dominant species of Viereck et al.’s (1992) bluejoint meadow. This vegetation type is likely a northwestern variant of the bluejoint meadow, but because we did not sample with the intent to classify plant communities we instead refer to this vegetation type as a less descriptive graminoid–forb tundra that warrants further description and research.

Effect of fire on ecosystem properties

Our multivariate descriptive ordination analyses revealed a complex and often correlated relationship between our metrics used to characterize fire regime (time after fire, burn severity, and times burned) and vegetation. Metrics to characterize fire regime were highly correlated with one another (Spearman’s ρ = −0.60 times burned and time after fire, Spearman’s ρ = 0.80 times burned and burn severity score, and Spearman’s ρ = −0.81 time after fire and burn severity score). Conceptually, burn severity represents the amount of biomass (particularly duff) consumed in the last fire and times burned (i.e., fire frequency) represents cumulative biomass consumed over multiple fires. However, due to the chronosequence approach of our study (as opposed to studying plots over time that have burned multiple times) the effect is similar. Times burned creates site conditions that mimic a more severe fire. Alternatively, time after fire allows biomass to reaccumulate, thus creating site conditions that mimic a low-severity fire.

We see these differences in site and vegetation characteristics that are tied to their recent fire histories as reflective of an intensifying fire regime (more severe and frequent fires) and therefore integrated these three variables into a synthetic index that represents an overall effect of fire through the use of principal component analysis (PCA) hereafter referred to as “fire effect.” PCA is a method that can distill a suite of correlated variables into one or a few variables and address issues of multicollinearity among variables (e.g., Hamann et al. 1998; Graham 2003). By synthesizing fire effect, we could more explicitly quantify correlations between our metrics used to characterize fire regime and specific ecosystem properties than in our initial descriptive ordination analysis (San-Miguel et al. 2020). This analysis did not incorporate species composition; instead, the analysis examined the relationship among our metrics used to characterize fire regime as a single fire effect variable and ecosystem properties. These data met the assumptions of normality. We ran a distance-based PCA (PC-ORD v7; McCune and Mefford 2017), and our cross-product matrix contained correlation coefficients between our three variables.

To test the effect of fire on above- and belowground ecosystem properties, we used linear mixed effects models in the nlme package (Pinheiro et al. 2017). Because we tested multiple response variables with the same experimental units in relation to fire effect, we corrected p values by controlling for false discovery rates with the Benjamini-Hochberg procedure (Benjamini and Hochberg 1995). We exported the first PCA axis and used it as our comprehensive fire effect variable. With separate models, we tested whether the fixed factor fire effect (PCA axis 1, the independent variable) was important in explaining the variability in belowground ecosystem properties (SOL depth, thaw depth, profile weighted C:N, C concentration, N concentration, C pool, N pool, fine pool, coarse woody debris pool, and surface soil C:N, the dependent variables) with quadrat nested in burn scar as a random factor to account for variability in observations based on environmental heterogeneity. When testing the relationship between fire and SOL depth, we excluded sites (n = 5) where we hit ice before reaching the bottom of the organic horizon (n = 38 vs. n = 43 plots).

We used the same model structure with fire effect as a fixed factor and quadrat nested in site as a random factor when testing fire effects on aboveground ecosystem properties, specifically related to a shift in the shrub:tussock ratio. We tested relationships between our synthetic fire effect variable and total species richness, bryophyte and deciduous shrub richness, shrub basal area and density, plot biomass, foliar C:N, % N, and % C of E. vaginatum, B. nana, and Salix spp. We excluded variables that did not meet the assumption of homoscedasticity (i.e., grass richness, shrub richness, and metrics of tussock abundance) and variables that were highly correlated with other aboveground variables.

Results

Effect of fire on species composition and plant functional types

Across all forty-three plots, alpha diversity (average number of species/plot ± SE) was 16.95 ± 0.75, beta diversity was 6.01 (total number of species/alpha diversity), and gamma diversity was 102 (total number of species across all plots). Species richness varied significantly among burn severity classes (ANOVA F = 6.97, p = .003), with richness generally declining as burn severity increased. Moderate-severity plots had the greatest species richness per plot (19.53 ± 1.09), significantly greater than that of high-severity plots (13.94 ± 1.07; Tukey’s post hoc p = .002). There was no significant difference, however, between species richness in high-severity and low-severity plots (17.40 ± 1.28; Tukey’s post hoc p = .129) or between moderate-severity and low-severity plots (Tukey’s post hoc p = .439).

The NMDS ordination had a final stress of 12.96, indicating that the ordination was a satisfactory representation of species composition patterns (McCune and Grace 2002). Three axes captured approximately 86 percent of the variance in species composition across all plots, and the final ordination solution was highly stable with an instability of 0.00001 after seventy-eight iterations (Figure 2). Axis 1 represented 35 percent of the variation in community composition and was not directly correlated with metrics to characterize fire regime. Instead, this axis was directly correlated with variables related to a gradient in composition from tussocks to shrubs (positively correlated with shrub basal area:tussock basal area and negatively correlated with tussock basal area and sedge percentage cover; Supplemental Table 2). The shift in tussock to shrub abundance is further demonstrated by the vegetation classes plots are grouped in the ordination (Figure 2). Although there is strong overlap in species composition along Axis 1 between the tussock tundra vegetation plots and the mixed shrub–tussock tundra vegetation plots, there is little overlap in species composition between tussock tundra plots and shrub tundra plots (Figure 2a).

Figure 2. NMDS ordination of all sites. Three axes captured approximately 86 percent of the variance in species composition across all plots. Ordinations were correlated with environmental and vegetation characteristics and grouped by vegetation class. Strong Kendall correlations (r² > 0.25) with the ordination axes are indicated by biplot vectors where plots are grouped by vegetation class showing (a) Axis 1 vs. Axis 2 and (b) Axis 2 vs. Axis 3.

The relationship between metrics to characterize fire regime is complex; therefore, we evaluated Axis 2 and Axis 3 together. Axis 2 represented 32 percent of the variation in species composition and was most directly linked to the effect of repeat burning (positively correlated with organic matter depth and site moisture and negatively correlated with grass cover, times burned, and elevation). Axis 3 represented 19 percent of the variation in species composition and was the axis most directly linked to a burn severity gradient (positively correlated with burn severity score, bare soil cover, and thaw depth and negatively correlated with time after fire, shrub basal area, and shrub/total biomass; Figure 2b). These two axes together contribute to the distinct graminoid–forb tundra vegetation plots in ordination space (Figure 2). Axis 2, however, also had strong correlations with variables generally indicative of shortened fire return interval, namely, times a plot burned and SOL depth. Axis 2 more accurately represented both a complex environmental and fire regime gradient, likely driven by the coupling of landscape position with site moisture, and likelihood of repeat burning and loss of SOL (Figure 2, Supplemental Table 2).

We were interested in evaluating whether there were noticeable shifts in plant functional types and, in particular, shrubs related to metrics used to characterize fire regime (Figure 3). Our results indicate that total shrub cover was aproximately 15 percent to 20 percent greater in moderate-severity (56.8 percent) or low-severity (57.9%) plots than in high-severity plots (40.0%). We found that total graminoid cover was greater (~15 percent) in high-severity plots (48.2 percent) compared to low-severity (34.6 percent) or moderate-severity plots (33.2 percent). Because percentage cover of sedges did not vary among severity classes (~30 percent for all severities), this difference in graminoids can be accounted for by shifts in the percentage cover of grasses. The most notable shift in plant functional types related to metrics to characterize fire regime was observed for graminoid cover, predominantly for the grasses Arctagrostis latifolia and Calamagrostis lapponica (Figures 2 and 3). The older burn scars had <1 percent grass cover, and the moderate and low burn (including unburned) severity tundra plots showed only 2 percent and 1 percent cover, respectively. In contrast, the younger burn scars and those with high burn severity had a mean grass cover of >15 percent. We also observed a trend toward lower bryophyte and lichen cover in high-severity plots.

Figure 3. Mean absolute cover of plant functional type, with plots grouped by time after fire, burn severity, and times burned classes.

The most notable shift in plant functional types related to metrics to characterize fire regime was observed for graminoid cover, predominantly for the grasses Arctagrostis latifolia and Calamagrostis lapponica (Figures 2 and 3). The older burn scars and sites that had not burned had <1 percent grass cover, and moderate and low burn (including unburned) severity tundra plots and once burned plots showed only 2 percent, 1 percent, and 3 percent cover, respectively. In contrast, the younger burn scars and those with high burn severity had a mean grass cover of >15 percent. Moreover, we observed nearly 30 percent grass cover in fire scars where repeat burning occurred. We also observed a trend toward lower bryophyte and lichen cover in high-severity plots.

Effect of fire on ecosystem properties

Quantifying fire effect

Axis 1 of the PCA described approximately 84 percent of the variability and is highly significant (p = .009). The eigenvalue rankings on Axis 1 include burn severity score = 0.96552, time after fire = −0.91067, and times burned = 0.88776.

Belowground ecosystem properties

With increasing fire effect, SOL depth was significantly reduced (χ² = 14.35, p < .001), concurrent with an increase in thaw depth (χ² = 21.07, p < .001; Figures 4a and 4b, Supplementary Table 3a). The C pool was also reduced with increasing fire effect (χ² = 5.22, p = .04), and the N pool showed a marginal reduction (χ² = 3.47, p = .06; Figures 4c and 4f, Supplementary Table 3a). The profile-weighted soil C:N was not sensitive to variation in fire effect, but the surface soil C:N increased with greater fire effect (χ² = 4.12, p = .04; Figure 4e, Supplementary Table 3a). The fine fuels pool was reduced with increasing fire effect (χ² = 5.44, p = .02), whereas the coarse woody debris fuels pool was not responsive to fire (Figure 4d, Supplementary Table 3a).

Figure 4. The response of belowground ecosystem characteristics to increased fire effects. Fire effect is a composite variable based on time after fire, number of times burned, and burn severity. Solid lines indicate significant relationships (p < .05), and dashed lines indicate marginally significant relationships (0.1 > p > .05). (c), (d), (f) represent profile-weighted pools. Surf = surficial (upper 5 cm of soil).

Aboveground ecosystem properties

Plot-level biomass marginally declined with a stronger fire effect (χ² = 4.69, p = .06; Supplementary Table 3b), whereas grass cover increased (χ² = 8.10, p = .03; Figure 5b, Supplementary Table 3b). Foliar C:N of Salix sp. was positively related to fire effect (χ² = 9.27, p = .03), apparently due to an increase in %C (χ² = 9.70, p = .03), concurrent with an increase in foliar %C for E. vaginatum (χ² = 13.78, p = .03) and a decrease in %N (χ² = 8.01, p = .03; Figures 5c, 5d, 5f–5h, Supplementary Table 3b). Although there was not a shift in the foliar C:N of Betula nana (Supplementary Table 3b), we did see a decline in %N with increased fire effect (χ² = 13.68, p = .03; Figure 5e). We did not observe an effect of fire on metrics representing shrub cover, basal area, or biomass (Supplementary Table 3b); however, total plot biomass likely reflects some shifts in shrub biomass due to the large percentage of shrub biomass in total biomass. These results show that the shrub plant functional type is less responsive to fire effects than anticipated. We could not conduct analyses on tussock metrics versus fire effect due to heteroscedasticity.

Figure 5. Response of aboveground vegetation characteristics to increased fire effects. Fire effect is a composite variable based on time after fire, burn severity, and times burned. Solid lines indicate significant relationships (p < .05).

Discussion

Post-fire successional trajectories and ecosystem properties are driven by the complex fire regime of the Seward Peninsula. In this region, not only is fire frequency increasing and burn severity intensifying (Racine, Johnson, and Viereck 1987; Rupp, Starfield, and Chapin 2000; Hu et al. 2015) but in the last fifty years reburning of tundra is a more common phenomenon (Liljedahl et al. 2007). The strong historical link between fire regime and tundra vegetation (Higuera et al. 2011), coupled with the observed climate-related increases in shrub dominance in Arctic Alaska—particularly in recently burned areas (Racine et al. 2004; Tape, Sturm, and Racine 2006; Lantz, Gergel, and Henry 2010; Myers-Smith et al. 2011)—led us to hypothesize that changes in metrics to characterize fire regime in northwestern Arctic Alaska would be tied to shifts from graminoid tussock tundra to a mixed shrub–tussock tundra or shrub tundra. We predicted that burned areas would show an increase in shrub abundance and a decrease in tussock abundance relative to unburned tundra and that as more fire-related shrub recruitment occurred, more biomass would accumulate.

Our ordination results suggest a strong correlation among our metrics used to characterize fire regime, though a large amount of the variation in sampled communities was attributed to a complex tussock to shrub gradient. Nevertheless, this is not a causal analysis, and we cannot conclude that tussock density and abundance are not related to fire regime. Tussocks dominate both fire-disturbed and undisturbed tundra sites and occur in all stages of succession, regardless of time after fire (Jandt and Meyers 2000); however, E. vaginatum cover, biomass, and density on burned sites was found to exceed that of unburned tundra even twenty-five years post-fire (Fetcher et al. 1984; Racine, Johnson, and Viereck 1987; Jandt and Meyers 2000). Despite the strong correlations we observed in overall plot species composition and variation in the tussock-to-shrub ratio (i.e., tussock basal area, shrub basal area:tussock basal area), these variables were not directly correlated with metrics to characterize fire regime in the descriptive ordination analysis. However, we saw significant negative correlations between fire effect and plot biomass in our quantitative mixed model analysis, which is likely a reflection of a decrease in shrub biomass, given the strong correlation between these two variables. Our results suggest a potential sensitivity of shrub abundance to fire and a potential resilience of E. vaginatum to fire, but these results in no way support the positive relationship between fire and increases in shrub abundance found in other tundra ecosystems (e.g., Lantz, Gergel, and Henry 2010).

Shrub and graminoid (attributed to grasses) plant functional types showed the greatest response to differences in burn severity. Regardless of time after fire, total shrub cover was highest in moderate- and low-severity plots, where shrubs survive belowground and resprout post-fire, and graminoid cover was greatest in high-severity plots, where shrubs have greater mortality. We expected high-severity fires that expose mineral soil in the tundra to promote increased shrub recruitment and, hence, cover of shrub species over time (Lantz, Gergel, and Henry 2010). In contrast, we found that high-severity fires promoted recruitment and perpetuation of grasses. This was also evident from our mixed model analysis that showed a positive relationship between fire effect and grass cover, whereas we observed no fire effect on our various shrub metrics. This finding is particularly of interest because we observed a shift from a complete absence of grasses in unburned tundra to a range of 7 percent to 79 percent cover following fire, most notably in areas where high-severity and repeat fires occurred. This graminoid–forb vegetation type was completely distinct in species composition.

Similar to the tallgrass prairie in central North America (Ratajczak, Nippert, and Ocheltree 2014), for graminoid–forb tundra to remain grass-dominated requires periodic disruption by wildfire to prevent shrub and tussock dominance. The observed “grassification” of tussock tundra has been described as a short-lived, early successional post-fire phenomenon on the Seward Peninsula and North Slope (Raynolds, Walker, and Maier 2005; Racine et al. 2010; Jones et al. 2013). This is a distinct variant of tussock tundra (Eriophorum vaginatum–Sphagnum spp. subtype Chamerion angustifolium) that Raynolds, Walker, and Maier (2005) reported includes the grasses Arctagrostis latifolia and Calamagrostis canadensis as indicator species. Our study did not document Calamagrostis canadensis but rather two other Calamagrostis species, Calamagrostis lapponica and C. neglecta, and for now we are describing a less descriptive graminoid–forb tundra that is more stable than an early successional variant of tussock tundra post-fire. Seeding by rhizomatous grasses is a locally important revegetation mechanism following a disturbance in the tundra (Forbes, Ebersole, and Strandberg 2001), including wildfire (Bliss and Wein 1972; Racine et al. 2004, 2010). However, we observed grasses in burn scars that were up to forty-six years old and therefore report for the first time a post-fire graminoid–forb-rich tundra plant community with the potential to persist and be self-perpetuating. At the time of our field study, the Kougarok fire complex experienced three major wildfires (Table 1; 1971, 1997, 2002). Since then, two more wildfires occurred within the study area, in 2015 (shown in Figure 1) and in 2019 (not shown).

In other tundra regions, post-fire tundra succession and vegetation development is strongly linked to active layer thickness (Fetcher et al. 1984; Bret-Harte et al. 2013). At the Kougarok fire complex, the 2002 fire was reported as moderate to severe with 50 percent of the 14-cm organic layer consumed by fire. Removal of the SOL precipitated a subsequent increase in soil temperature, earlier freezing, and deeper thaw that persisted five years post-fire (Liljedahl et al. 2007), as well as the development of thermokarst subsidence (Iwahana et al. 2016). These belowground shifts in ecosystem properties were strongly linked to vegetation structure, where deep thaw corresponded to graminoid-rich areas and shallower thaw corresponded to shrub-rich areas (Narita et al. 2015). Our results indicate that fire effect was positively correlated with thaw depth and negatively correlated with SOL (as anticipated), and these trends also corresponded to a decrease in overall plot biomass and an increase in grass cover. The increased grass dominance associated with increased thaw depth has been observed in other tundra ecosystems (e.g., Natali, Schuur, and Rubin 2012), where grasses have faster turnover and are better able to access nutrients released by thawing. In addition, arctic grasses are observed to generally be more deeply rooted than deciduous shrubs (Iversen et al. 2015).

An increased fire effect resulted in foliar traits and soil characteristics that suggest plant-available N limitations with increasing fire severity. With an increase in fire effect, we found a significant increase in foliar C:N of Salix species, concurrent with marginal increases in C concentration and a decline in N concentration for E. vaginatum. Betula nana foliage also showed a negative relationship between fire effect and N concentration. At the same time, we saw a marginal decline in the total soil N pool, and surface soil (top 5 cm) C:N increased as fire effect increased. Together, these patterns could indicate reduced plant N availability with increased fire effect. Fire is a “pulse” disturbance, where surface organic N is rapidly volatilized (Alexander and Mack 2016), with a potentially long time to reaccumulation (Mack et al. 2011).

Our results also suggest a relationship between increased fire effect and reduced forage quality and quantity of S. pulchra and S. glauca, which are browse species for caribou, moose, and ptarmigan (Christie et al. 2014). An increase in shrubs across the Arctic has resulted in changes to the quality and quantity of forage for northern ungulates (Tape, Sturm, and Racine 2006); 45 percent to 50 percent of spring and summer forage for caribou in Alaska consists of deciduous shrubs, including willow and birch species (Thompson and Barboza 2014). Our observations of species-level variation in foliar C:N in response to fire suggest that, regardless of shifts in abundance, fire is uniquely affecting shrub foliar traits that could have cascading effects on habitat (Mekonnen, Riley, and Grant 2018).

The climate and local flora of the Seward Peninsula are linked to the North Slope tundra, where underlying parent material and cold temperatures contribute to the lack of trees (D. A. Walker et al. 2018). However, the complex fire regime of the Seward Peninsula and the lasting effect this has on vegetation communities are more similar to the boreal forest, where fire is an integral component of the landscape (Chapin et al. 2010). For example, the strong negative relationship we found between fire effect and SOL depth is much more similar to the patterns of burn severity and SOL in the boreal forest (Boby et al. 2010) than other arctic tundra regions where post-fire SOL has been measured (e.g., Mack et al. 2011). It is likely that, due to the repeat burning that we observed, there was less variability in pre-fire SOL depths and therefore residual or post-fire SOL is a better proxy for burn severity than in other less fire-prone tundra ecosystems. In boreal forest systems, successional trajectory is strongly linked to burn severity (Johnstone, Hollingsworth, and Chapin 2008; Hollingsworth et al. 2013), where low-severity fires promote relay succession and a quick return to pre-fire vegetation and soil conditions, whereas high-severity fires can trigger a shift in successional trajectory and the recruitment of novel species, both native (Johnstone et al. 2010) and invasive (X. J. Walker et al. 2017). The effect of burn severity on vegetation and recruitment in Alaska tundra is documented but limited (Bret-Harte et al. 2013; Hewitt et al. 2016). However, our study area experienced repeat burning, thus shortening the time after fire, which can be analogous to a shortened fire return interval and undoubtedly exacerbates fire effects on vegetation succession. There is evidence that at the northern limit of black spruce distribution in the boreal forest of the Yukon Territory, repeat burning and short (fifteen-year) fire return intervals limited tree germination due to the buildup of grass litter (Brown and Johnstone 2012). Although our research area contained no trees, we also observed an increase in grass and subsequent litter as well as a decrease in tussocks, mosses, and lichens, with an increase in burn severity and decrease in time after fire. Together, those two observations likely limit the germination of shrubs, as has been observed in tundra systems that burned at high severity (Barrett et al. 2012). Our results suggest that repeat fires or shortened fire return intervals, as well as an increase in burn severity, shift tussock tundra toward a more grass–forb-dominated tundra with significantly fewer mosses and lichens.

Tundra post-fire succession is not yet fully understood, and the few long-term studies published document a quick return to pre-fire tussock tundra (Bret-Harte et al. 2013) or transition to shrub tundra (Racine 1981; Barrett et al. 2012; Jones et al. 2013).

However, based on our observations, we put forth an alternative post-fire successional trajectory, in which grass dominance persists and is self-perpetuating post-fire, with repeat burning (Figure 6). In tussock tundra, high moss accrual results in deep soil organic layers, fostering shallow active layers and low microbial decomposition due to cold soil temperatures (Hobbie et al. 2000). These conditions in turn favor low-severity fires that return tussock tundra to its pre-fire state, similar to the return to a pre-fire black spruce state in the cold, wet sites of boreal Alaska (Johnstone et al. 2010). Under the historical tundra fire regime, this internal feedback has likely been a source of ecosystem resilience due to cool, moist soils that are poorly drained and underlain by permafrost. However, changes in vegetation via the direct effects of climate change, coupled with shifts in tundra fire regimes, have the potential to result in alternative successional trajectories. As shrub cover increases in tussock tundra as a direct result of changes in climate (Myers-Smith et al. 2011), there is a shift in biomass and fuels. In the historic past, this shift from graminoids to shrubs is linked to changes in fire regime (Higuera et al. 2008). After high-severity fires that burn deeply into the SOL, we hypothesized that shrub species establishment is favored and leads to a switch to an alternative plant successional trajectory dominated or codominated by shrubs. However, changing fire regimes—with shifts toward shorter fire return intervals, increases in the area burned and reburned, and higher severity fires—may in fact have the opposite effect: decreased woody plant recovery and development of novel, grass-dominated graminoid–forb-rich tundra. Because our study did not include older (more than thirty-one years) high burn severity plots, our observation of a potential alternative graminoid–forb-rich tundra successional trajectory warrants further study.

Figure 6. Post-fire successional trajectories in Alaska arctic tundra. Arrows 1–4 depict trajectories among vegetation types and arrows 5–7 depict fire regimes that maintain existing vegetation. Colors represent observed (black), hypothesized (gray), and novel (orange) pathways. The trajectories among vegetation types include (1) tussock tundra to shrub tundra with a long fire-free interval and climate warming through a shift in dominance of existing shrubs and subsequent shading, accumulation of deciduous leaf litter, and loss of the understory (Epstein et al. 2004); (2) tussock tundra to shrub tundra after high-severity fire through combustion of vegetation and soil organic layers, thereby exposing mineral soil and facilitating shrub recruitment (e.g., Landhausser and Wein 1993; Racine et al. 2004; Jones et al. 2013); and (3), (4) tussock tundra (this study) and shrub tundra to a novel graminoid–forb tundra after high severity fire and repeat burning. (5), (6) Fire regimes that maintain existing vegetation include low-severity fire maintains tussock tundra (e.g., Bret-Harte et al. 2013; this study) and shrub tundra and their vegetation and ecosystem properties (7) and a short fire return interval (this study) maintains graminoid–forb tundra with its distinct vegetation and ecosystem properties.

We observed strong and significant shifts in vegetation and soil properties related to metrics to characterize fire regime. In contrast to our hypotheses, however, we found that increases in burn severity were associated with decreases in shrub abundance, and plot biomass was negatively correlated with our integrated variable of fire effect. Our results challenge the prevailing hypothesis that more fire activity will lead to higher plant productivity in tundra. Rather, we observed increased grass abundance with an intensifying fire regime. Incorporating multiple metrics to characterize fire regime into our analyses enabled a more thorough understanding of the complex relationship between fire and vegetation on the Seward Peninsula, indicating the continued need for research that addresses the effects of fire on successional trajectories in tundra ecosystems.

Data availability

All data associated with this project are available at: http://www.lter.uaf.edu/data/data-detail/id/752.

Acknowledgments

We thank four anonymous reviewers for their constructive feedback that greatly improved the article. We also thank Kirsten Barrett for field support, Angie Floyd and Jamie Hollingsworth for creating Figure 1, and Carolyn Rosner for creating Figure 6. Bonanza Creek LTER (jointly funded by NSF and PNW Research Station, USDA Forest Service) provided logistical support, and the National Park Service and Bureau of Land Management provided housing in Kotzebue, Alaska.

Disclosure statement

The authors declare no conflict of interest.

Supplementary material

Supplemental material for this article can be accessed on the publisher’s website.

Additional information

Funding

This research was supported by grants from the USGS Alaska Science Center (RW0195) and the National Science Foundation (DEB-1026415).