Boreal forest vegetation and fuel conditions 12 years after the 2004 Taylor Complex fires in Alaska, USA

Background: Fire has historically been a primary control on succession and vegetation dynamics in boreal systems, although modern changing climate is potentially increasing fire size and frequency. Large, often remote fires necessitate large-scale estimates of fire effects and consequences, often using Landsat satellite-derived dNBR (differenced Normalized Burn Ratio) to estimate burn severity. However, few studies have examined long-term field measures of ecosystem condition in relation to dNBR severity classes in boreal Alaska, USA. The goals of this study were: 1) assess changes in dominant vegetation at plots resampled one and 12 years post fire; 2) use dNBR classes to characterize vegetation and downed woody fuels 12 years post fire; and 3) characterize the relationship between biophysical, topographic, and remotely sensed characteristics (e.g., moss and duff depth, canopy cover, elevation, aspect, dNBR) and understory species assemblages 12 years post fire. Results: Understory species richness doubled (from 39 to 73) between 2005 and 2016; some common species increased in cover over time (e.g., Ledum groenlandicum Oeder) while others decreased (e.g., Hylocomium splendens [Hedw.] Schimp.). In 2016, live and dead tree densities, tall shrub cover, and 1and 100-h woody fuels were significantly different among dNBR classes; moss and duff depth, canopy cover, and spruce seedling density were not. Elevation and aspect significantly influenced tall shrub cover, hardwood sapling density, and downed woody fuel loads. Understory plant communities differed between unburned and all burn classes, as well as between low and high dNBR severity. Ordination analysis showed that overstory (e.g., live tree density), understory (e.g., moss depth, woody fuel loading), and site (elevation, aspect, dNBR) significantly influences understory species assemblages. Conclusion: Remeasured sites (sampled one and 12 years post fire) showed recruitment of new understory species and differing, diverse responses to burning by several common plant species. In 2016, low-severity burned sites had generally the highest woody fuel loading, which may increase risk of repeated surface burning, although the reduction in live tree density would still result in decreased fire risk and behavior. Understory community composition correlated with multiple biotic and abiotic factors, including moss depth, canopy cover, elevation, aspect, and dNBR. Overall, our findings can improve landscape-level predictions of ecosystem condition following fire based on dNBR.


BA:
Basal Area (m 2 ha −1 ) BAF: Basal Area Factor DBH: Diameter of a tree measured at Breast Height (1.37 m) dNBR: differenced Normalized Burn Ratio NBR: Normalized Burn Ratio NMDS: Non-metric MultiDimensional Scaling perMANOVA: permutational Multivariate ANalysis Of VAriance SPH: tree Stems Per Hectare TRASP: TRansformed ASPect where 0 is northnortheast and 1 is south-southwest

Background
Wildfire has been a historically prevalent disturbance in the North American boreal forest (Viereck 1983;Ryan 2015), although in recent decades there has been an increase in annual area burned and frequency of large fires (Kasischke and Turetsky 2006;Kasischke et al. 2010;Calef et al. 2015). North American boreal forests characteristically burn in a stand-replacing fire regime (Duchesne and Hawkes 2000), in which crown fires are common and there is high tree mortality even in surface fires (Viereck and Schandelmeier 1980;Viereck 1983). Surface fires generally move easily into the crowns of spruce-dominated forests due to the distribution of fine branches and flammable lichen along the entire tree stem (Viereck 1983).
Unless fires burn into late summer under dry and hot conditions, fire generally reduces the forest floor organic layer without removing it down to bare mineral soil (Viereck 1983;Duchesne and Hawkes 2000). Thus, overstory consumption typical of high burn severity is not necessarily coupled to severe effects on understory vegetation (Fryer 2014). In this context, fires can produce highly variable effects on boreal forest overstory and understory as an agent of ecosystem change (Simard 1991;Lentile et al. 2006). Picea mariana (Mill.) Britton, Sterns & Poggenb. (black spruce) is a shade-tolerant species that is a dominant component of fire-prone Alaskan, USA, forests. Historically, P. mariana is generally successful in replacing itself after fire on moister sites and sites burned at low to moderate burn severity level when some organic layer is left unconsumed ). Germination success is variable but consistently higher soon after fire, with establishment from semi-serotinous cones beginning the first year after fire and continuing for up to 10 years (Viereck 1983;Duchesne and Hawkes 2000). The other dominant conifer in Alaska boreal forests, Picea glauca (Moench) Voss (white spruce), often decreases in dominance following fire except when fire coincides with a masting event (Peters et al. 2005). Hardwood trees in this region of the Alaskan boreal forest generally increase in density following fire and may dominate for years to decades before being replaced by spruce (Viereck 1983). Fires in young (<25 years old) stands or areas where organic layer consumption is high may reduce spruce regeneration to the point that hardwoods are never replaced as the dominant tree (Johnstone and Chapin 2006a; Johnstone et al. 2010). Although aboveground structures are easily killed even by low severity surface fires, Populus tremuloides Michx. (quaking aspen), Populus balsamifera L. (balsam poplar), and Betula papyrifera Marshall (paper birch) all resprout vigorously after high-severity burns. Betula papyrifera in particular can also seed in aggressively following fire, particularly in areas where the organic layer is consumed (Johnstone and Chapin 2006b).
Post-fire understory species succession in boreal forests generally proceeds from dominance of rapidly growing, vascular plants within the first few years following fire to increasing dominance of bryophytes (mosses and liverworts) and lichens as time since fire increases (Hart and Chen 2006). This succession can be measured by either resampling the same sites through time in successive years, or by sampling different sites along a known chronosequence of post-fire stand ages. Chronosequence sampling can alleviate some issues related to resampling sites and generally provides a longer view of time since fire, since it is often difficult to maintain resampling on the order of 100+ years necessary to capture the full range of succession in these slow-growing forests (e.g., Black and Bliss 1978;Morneau and Payette 1989;Chipman and Johnson 2002). However, resampling the same site has the benefit of tracking changes through time without the compounding factors associated with sampling different locations (e.g., differences in soil or microclimate). While generally a shorter timescale, these remeasurement studies allow the explicit relation of immediate post-fire plant communities to longer-term trends. Remeasurement studies have, however, shown generally similar trends in species composition through time as the chronosequence-based studies, while also highlighting the fact that time since disturbance is generally not enough to fully explain differences in observed understory species communities (Rydgren et al. 2004;Day et al. 2017).
Because live and decaying mosses and lichen create a deep organic layer that dominates the forest floor, particularly in later successional stages, fire effects on understory species assemblages in boreal forests are considered to be mostly driven by the depth of organic layer consumption (Wang and Kemball 2005;Gibson et al. 2016). Depth of this organic layer is a strong control on soil temperature and resulting permafrost layer melting or formation, which in turn affects site soil drainage and related factors (Viereck 1983;Yoshikawa et al. 2002;Kasischke and Johnstone 2005). Lower organic consumption by fire means that mosses (e.g., Pleurozium schreberi [Brid.] Mitt. [red-stemmed feathermoss], Hylocomium splendens [Hedw.] Schimp. [splendid feathermoss]) generally continue to dominate, whereas increased depth of burn results in a shift to higher forb cover, often Chamerion angustifolium L. Holub (fireweed) or grasses such as Calamagrostis canadensis (Michx.) P. Beauv. (bluejoint). Shrubs, including Salix spp. (willow), Vaccinium uliginosum L. (blueberry), Alnus Mill. spp. (alder), and Ledum groenlandicum Oeder (synonym Rhododendron groenlandicum, bog Labrador tea) frequently gain significant cover following fire and may dominate for up to 25 years post fire before succeeding to Picea mariana or P. glauca (Duchesne and Hawkes 2000).
These varied ecological effects following fire are dependent in part on burn severity, which is the degree of ecological change following fire (Lentile et al. 2006;Keeley 2009). In many forested systems, percent of overstory tree mortality is successfully used as a proxy for ecosystem burn severity, (i.e., both overstory and understory ecological effects; Lentile et al. 2006;Keeley 2009). Various indices have been developed to assess burn severity based on both surface and overstory effects (e.g., Composite Burn Index, CBI; Key and Benson 2006;Morgan et al. 2014). While field-based assessments of burn severity (e.g., CBI) are considered to be the preferred method for determining burn severity, these indices may not accurately capture ecosystem change due to collapsing of complex post-fire changes into simplified measures (Morgan et al. 2014).
Burn severity across larger scales or remote areas is commonly assessed using remote sensing (Lentile et al. 2006), with the most common method being the differenced Normalized Burn Ratio (dNBR; Eidenshink et al. 2007). However, particularly in boreal systems, some studies have found mixed success in the ability of remotely sensed indices to predict ground measures of burn severity and, particularly, that these indices fail to capture depth of organic layer burning Verbyla and Lord 2008;Alonzo et al. 2017). Other work, however, has found that dNBR is broadly effective at mapping burn severity in boreal systems as compared to field measures of CBI (Allen and Sorbel 2008;Boucher et al. 2016). Field-based methods for burn severity assessment (such as percent canopy mortality and CBI) have also been criticized as not characterizing ecologically significant burn severity in boreal systems, which further complicates validation of dNBR in this ecosystem Kasischke et al. 2008). While more ecosystem-specific methods for determining burn severity have been developed, such as estimation of organic layer consumption by looking at adventitious roots on Picea mariana Boby et al. 2010), there is still a need for remotely sensed burn severity estimation in order to characterize large-scale burn severity (Burton et al. 2008;Parks et al. 2018).
Although there is a significant body of research focused on post-fire effects in boreal forests, some gaps in knowledge still remain. While successional phases following wildfire have been well characterized in boreal ecosystems, studies explicitly considering differences in burn severity have usually defined severity using depth of organic layer consumption or similar (e.g., Mallik et al. 2010;Gibson et al. 2016), rather than directly using dNBR or a similar index as a measure of severity. Multiple studies have also focused on the relationship of remotely sensed indices to short-term, field-based measures of burn severity in boreal systems (e.g., Hoy et al. 2008;Murphy et al. 2008;Whitman et al. 2018), but few have focused on relating dNBR to field-based measures of long-term overstory and understory condition following wildfire. In addition, long-term measures of vegetation recovery should provide the more telling assessment of the ecological impact (i.e., severity) of the fire.
While acknowledging its shortcomings in boreal ecosystems, we proceeded to consider dNBR as our remotely sensed burn severity index for this study, for practicality and for lack of a clearly preferable alternative ). Since dNBR is the most commonly used burn severity index by researchers and managers, assessing how dNBR corresponds to post-fire ecological change is still warranted . We sampled 32 sites across four burn severity classes (unburned, low, moderate, and high severity) 12 years after the Taylor Complex fires, while also capitalizing on previously collected one-year post-fire data at a subset of our sites. The Taylor Complex (528 400 ha) was the largest fire complex in Alaska's record-breaking 2004 fire season, when 2.7 million ha burned statewide. The goal of this study was to survey long-term ecosystem condition following the Taylor Complex fires with three main objectives: 1) assessing changes in vegetation from one to 12 years post fire at resampled plots; 2) characterizing vegetation and downed woody fuels conditions 12 years post fire in relation to dNBR, elevation, and aspect; and 3) characterizing the influence of site biotic and abiotic characteristics (e.g., woody fuel loading, shrub cover, moss and duff depth, elevation, dNBR) on understory species assemblages.

Study area
We conducted this study in Picea mariana-dominated boreal forests along the Taylor Highway, primarily between the communities of Tok and Chicken in Interior Alaska. The climate is subarctic dry continental. Based on 30-year climate normals (1981 to 2010), the average daily January high temperature was −24°C with an average low of −34°C; in July, the average daily high was 22°C with an average low of 5°C. Average annual precipitation was 315 mm; the driest month was April with 6 mm of precipitation (generally snow), and the wettest month was July with 70 mm of rain (ACRC 2018). In this region, Picea mariana often co-dominates with Picea glauca. Other common tree species are Populus tremuloides, Populus balsamifera, and Betula papyrifera. This area is in the discontinuous permafrost zone, where the frozen permafrost layer is patchily distributed and may vary from being absent to several hundred meters thick (Jorgenson et al. 2008).

Landscape stratification and sampling design
Burn severity was indicated by the Monitoring Trends in Burn Severity (MTBS; Eidenshink et al. 2007) dNBR (differenced Normalized Burn Radio) product (http://www. mtbs.gov). The MTBS product was used because it is a free, validated source of burn severity maps for all large fires in the US and is commonly used by both researchers and managers for post-fire ecological studies, impact assessments, and to guide potential restoration. To create this product, continuous dNBR values are calculated as the pre-fire NBR (Normalized Burn Ratio) minus post-fire NBR, where NBR is a normalized ratio calculated from 30 m resolution Landsat bands 4 (near infrared band, NIR) and 7 (short wave infrared band, SWIR); (i.e., (NIR-SWIR)/(NIR+SWIR)). Continuous dNBR is then categorized as unburned, low, moderate, or high severity, or increased green vegetation based on thresholds determined by the MTBS analyst (Eidenshink et al. 2007). The NBR index is sensitive to changes in green vegetation, bare soil, and char, where areas with higher char and bare-soil cover fractions will have higher NBR values (Key and Benson 2006;Hudak et al. 2007). For image dates and proportion of burn severity in each severity class for each fire, see Additional file 1.
Burned sites were chosen from the Taylor Complex (http://www.mtbs.gov). Because these fires were part of the same complex, burned over the same time period, and all burned within the same vegetation types, they were treated as a single burn event when stratifying the sampling sites. Twelve sites were established in 2004 within low (n = 4), moderate (n = 5), and high (n = 3) burn severity conditions based on the full range of burn severity observed in the field Lentile et al. 2007), with some sites overlapping established pre-fire sites (Lewis et al. 2011). These 12 monumented sites were sampled in 2005 to collect one-year post-fire understory species composition and cover data , with each site consisting of five plots systematically placed within a 60 m diameter area (i.e., sampling 4 to 5 Landsat 30 m pixels, see Fig. 2). The 12 sites ended up falling within either low (n = 4), moderate (n = 5), or high (n = 3) burn severity classes, based on post-fire dNBR classified by MTBS (Additional file 1 for MTBS details). In 2016, in addition to remeasuring the 12 established sites, 20 new sites were added to achieve a balanced sampling design with regard to dNBR severity class (i.e., 32 total sites, with eight sites per unburned, low, moderate, and high severity classes). Because the original 12 sites did not include any unburned areas, all unburned sites in this study were established and sampled in 2016. Site selection followed a stratified random sampling design, with the strata being the four MTBS dNBR classes (unburned, low, moderate, high); high or low elevation (range 424 to 1529 m); and cool-wet or warm-dry aspect based on the topographic solar-radiation index (TRASP) transformation (Roberts and Cooper 1989) that varies from 0 to 1 (0 = cool, moist, 30°north-northeast, and 1 = warm, dry, 210°south-southwest). Elevation and TRASP were obtained from a 30 m resolution digital elevation model (https://viewer.nationalmap.gov/basic/). For accessibility, and to mitigate potential edge effects, sites were located within 0.8 km of, but at least 60 m from, the Taylor Highway.
Each site in 2005 and in 2016 consisted of five plots (Fig. 2). Field methods were generally described at the plot level, but measurements were summed or averaged to the site level for statistical analysis. The center of every plot was geolocated in 2005 and re-geolocated in 2016 with a Trimble GeoXT Global Positioning System (Trimble Inc., Sunnyvale, California, USA) with differential correction to 1 to 2 m precision. At all five plots per site, understory vegetation was identified to species (when possible) and the percent cover of each species was estimated within a 1 m 2 frame. Canopy closure (all trees and tall shrubs above 1.37 m) was estimated using a convex spherical densiometer. At the center plot, all standing trees were recorded if their base was located within a fixed 8 m radius. For each tree the species, status (dead, live), and diameter at breast height (DBH, 1.37 m) were recorded.
In 2016, additional variables were measured at all 32 sites (Fig. 1). Moss or litter and duff depth in millimeters was determined by cutting down through all organic layers until mineral soil was exposed, and fine downed woody fuel loading (1-, 10-, 100-h timelag classes, <1 cm, 1 to 2.5 cm, and 2.5 to 7 cm diameter, respectively) was estimated using a photoload guide (Keane and Dickinson 2007). Woody fuels were categorized by timelag classes, as is standard for many fire behavior and effects applications, which reflect estimated drying time of fuels of different sizes (Keane and Dickinson 2007). Spruce seedlings (<1.37 m height) and spruce and hardwood saplings (>1.37 m height and <10 cm DBH) were tallied within a 5.6 m radius of each plot center. Standing trees (>1.37 m height and >10 cm DBH) were recorded at each peripheral plot (four per site) within a variable radius using a 2 m 2 ha −1 basal area factor prism. Percent cover of tall shrubs and large downed woody fuel loading (1000-h timelag class, >7.62 cm diameter per Keane and Dickinson 2007) were ocularly estimated within a 5.6 m radius of the center plot at each site.

Data analysis
All analyses were run at the site level (n = 12 in 2005, n = 32 in 2016); plot-level (n = 60 in 2005, n = 160 in 2016) measurement data was aggregated to provide sitelevel measures and to preclude pseudo-replication due to spatial autocorrelation between plots within sites. The analysis for each objective, including any data transformation or subsetting, was done separately, as outlined below.

Vegetation change 2005 to 2016
These analyses were run primarily on the 12 sites sampled in both 2005 and 2016, although the unburned sites (sampled only in 2016) were used as well for some tests. Density of live trees, dead trees, and canopy cover in 2005 and 2016 was compared using a Wilcoxon signedrank test (v-statistic, V) in R (R Core Team 2017). To examine differences in understory species cover over time, we calculated total richness based on all species found on the sites. We then focused on changes in Fig. 2 Site layout (not to scale) showing central plot and four peripheral plots, along with the attributes measured within each plot. Note that shrub percent cover was measured only at the center plot, and that the variable-radius plot was based on a 2 m 2 ha − 1 Basal Area Factor prism percent cover of only common species, defined as any species that occurred in at least six of the 12 sites in both years. Narrowing our statistical testing to common species was done because it simultaneously highlighted species that are both ecologically and socially significant, and because it gave sufficient statistical power to run appropriate tests. For these common species, we tested for differences in cover between years using the Wilcoxon signed-rank test. We also tested for differences in cover among unburned sites, burned sites in 2005, and burned sites in 2016 using a Kruskal-Wallis non-parametric analysis of variance with post-hoc Dunn's test for pairwise comparisons (McKnight and Najab 2010). Test significance was defined when α < 0.10 for all tests. The Kruskal-Wallis was run on all severity levels combined because of the high variation inherent in plant cover data and because of the small number of sites (three) in high severity.

vegetation and fuels condition
We examined the relationship of field variables (tree and seedling density, canopy and shrub cover, downed woody fuel loading, moss and duff depth) to dNBR using Kruskal-Wallis analysis of variance to test for differences among categorical burn severity levels (unburned, low severity, moderate severity, high severity), with post-hoc Dunn's test, for all 32 sites sampled in 2016. All tests were run with significance defined at α < 0.10 with three degrees of freedom. To test the relationship between field variables and elevation and transformed aspect, we used Spearman's rho (ρ) rank correlations (McCune and Grace 2002). Spearman's ρ was also used to examine relationships between elevation and transformed aspect and understory plant community species evenness, species richness, and Shannon and Simpson indices of diversity (calculated in PC-ORD; McCune and Grace 2002; McCune and Mefford 2011). Correlations and P values were calculated using cor.test() in the R base package with significance defined at α < 0.10.

understory species community
Differences in understory species composition among burn severity levels were tested using a permutation multivariate analysis of variance (perMANOVA; Anderson 2001) in PC-ORD. The distance measure selected was Sorensen (Bray-Curtis; McCune and Grace 2002) and the number of randomizations was 4999. To test for differences in beta-diversity between severity levels, the homogeneity of group dispersions (or the difference in amount of within-group variance) was tested using the default methods with betadisper() (vegan package; Oksanen et al. 2018) to implement the multivariate test for homogeneity of dispersion described by Anderson (2006). This test can also serve as a proxy measure of beta-diversity; (i.e., the dissimilarity in species composition among sampling units for a given area; Anderson et al. 2006). The betadisper() test was run with three degrees of freedom and 999 permutations. An indicator species analysis (Monte Carlo method; Dufrêne and Legendre 1997) in PC-ORD determined the degree to which different species served as an indicator of a particular severity level based on both a species' consistent occurrence in, and relative abundance at, a given severity level.
The understory vegetation composition, and its association with various biotic and abiotic site characteristics (e.g., moss and duff depth, canopy cover, dNBR), were visualized using non-metric multidimensional scaling (NMDS) implemented with the metaMDS() wrapper from the vegan package (Oksanen et al. 2018), which, by default, implements the recommendations of Minchin (1987) in R. NMDS was chosen because it does not rely on an assumption of linear relationships among variables and allows examination of the trends in understory species in relation to multiple factors while reducing noise inherent in species diversity data (McCune and Grace 2002). The number of zeros in the data was reduced by excluding rare species that occurred at only a single site, resulting in a reduction from 118 to 74 species. A threeaxis solution with Bray-Curtis distance measure was chosen, with the default metaMSD() Wisconsin squareroot transformation (Minchin 1987) of the data to correct for the large range in cover values. The starting configuration was random and dimensionality was determined based on a scree plot of stress and dimensions (run in PC-ORD), whereby using four dimensions did not measurably improve the stress of the model (McCune and Grace 2002). Site characteristic vectors (e.g., dNBR, fuel loading) were then fit to the final ordination solution using the envfit() function (vegan package).

Vegetation change 2005 to 2016
There were no significant differences at the sites between 2005 and 2016 live tree density (mean = 60 ± 39 and 37 ± 25 trees ha −1 , respectively), dead tree density (mean = 280 ± 84 and 173 ± 36 trees ha −1 , respectively), or canopy cover (mean = 23 ± 3 and 20 ± 4%, respectively) (all P ≥ 0.181, V = 0, 18, 41.5, respectively). The overall or global understory species richness of the sites increased, from 39 species identified in 2005 to 73 species in 2016. Six understory species were common at the sites in both 2005 and 2016 (Additional file 2). When tested pairwise for differences in cover between years, only Ledum groenlandicum (V = 9, P = 0.037), Salix spp. (V = 1, P = 0.003), and Vaccinium uliginosum (V = 2, P = 0.012) had significant changes, all increasing in cover over time. When tested for differences in cover among unburned sites, 2005 burned sites, and 2016 burned sites (Fig. 3), Chamerion angustifolium had significantly higher cover at burned sites in both years than in unburned sites; Hylocomium splendens had significantly lower cover in 2016 than in 2005 or in unburned sites; Ledum groenlandicum had significantly lower cover in 2005 than in 2016 or in unburned sites; Salix L. spp. had significantly higher cover in 2016 than in unburned sites or in 2005; Vaccinium uliginosum had significantly higher cover in 2016 than in 2005 or in unburned sites; and Vaccinium vitis-idaea L. (lingonberry) had significantly lower cover in 2005 than in unburned sites.
Comparison of seedling and sapling density showed no significant difference among severity classes for separate spruce or hardwood species (data not shown) or for lumped total spruce or total hardwood seedling and sapling density by severity class (Fig. 6). Downed woody fuel loading in the 1-h fuels class (Fig. 7A) was significantly higher (P = 0.006) in low severity (0.052 ± 0.030 kg m −2 ) than in high severity (0.023 ± 0.011 kg m −2 ) burned areas, or in unburned (0.026 ± 0.009 kg m −2 ) areas (P = 0.05), and higher in moderate (0.052 ± 0.046 kg m −2 ) than in high severity (P = 0.06) burned areas. Similarly, 100-h fuels (Fig. 7C) were significantly (P = 0.038) higher in low severity burned areas (0.234 ± 0.218 kg m −2 ) than in unburned areas (0.033 ± 0.0478 kg m −2 ). There were no significant differences among severity classes for 10-h and 1000-h fuels ( Fig. 7B and D). Elevation (Fig. 8) was negatively correlated with canopy cover, hardwood sapling density, and 1000-h and 100-h fuel load, while positively Chamerion angustifolium, Hylocomium splendens, Ledum groenlandicum, Salix spp., Vaccinium uliginosum, and Vaccinium vitis-idaea. Boxplots show median, hinges (25th and 75th percentiles), and whiskers that extend to values no more than 1.5*IQR (inter-quartile range) from upper or lower hinge, respectively. Points outside this range are considered outliers and plotted separately. All Kruskal-Wallis tests were significant among unburned, all burned sites in 2005, and all burned sites in 2016 (df = 2, χ 2 ≥ 5.29, P ≤ 0.07); Dunn's pairwise test for differences are shown with letters to indicate significant pairwise differences (α < 0.1) correlated with tall shrub cover and duff depth. Transformed aspect (Fig. 9) was positively correlated with 1-, 10-, and 100-h fuel loads as well as with duff depth.

understory species community
Species richness, evenness, and diversity (Shannon's and Simpson's) were not significantly correlated with elevation, transformed aspect, or any of the field measurements (Additional file 3), except 100-h fuel load. Species richness and both Shannon and Simpson indices were negatively associated with 100-h fuel loading (all ρ ≥ |0.321|, P ≤ 0.07). The beta-dispersion analysis showed that, while low burn severity had the highest average distance to median (0.53), indicating higher beta-diversity, and high severity had the lowest (0.39), these differences were not significant based on a permutation test (data not shown; df = 3, F = 1.71, P = 0.18). Species richness of unburned areas (i.e., the total number of species found at all eight unburned sites) was 68, low severity burned area richness was 84, moderate severity burned area richness was 64, and high severity burned area richness was 58. There were 38 understory species (representing 32% of total species sampled) that were found on plots of all severity levels, including unburned, and 17 species were found at over 50% of the sites (Additional file 4). Polytrichum Hedw. moss and Vaccinium vitis-idaea were the most common, occurring at 97% of sites, and unidentified fungi of multiple species occurred at 94% of sites. In contrast, 13 species were found exclusively in unburned areas, 18 species exclusively in low severity burned sites, nine species in moderate severity burned sites, and seven in high severity burned sites. Many of these severity-specific species were found only at a single site, which may indicate that they simply represent rare species rather than species with strong severity-specific preferences.
Understory plant communities differed significantly with burn severity (perMANOVA, df = 3, F = 2.06, P = 0.002), where unburned sites had significantly different communities than low (P = 0.10), moderate (P = 0.001), and high severity (P < 0.001) burned sites. Low and high severity burned areas were also significantly different (P = 0.05). There were also several species of vascular plants, lichen, and moss that were significant indicators of unburned and high severity burned sites (Table 1).
These differences in understory community composition can be partially visualized in an NMDS ordination (3-dimensional [k = 3], type I stress [weak ties] equal to 0.134, 20 tries for two convergent solutions; Fig. 10). The first two axes contain the majority of variation (axis 1 R 2~0 .53, axis 2 R 2~0 .20) and are therefore the axes displayed. Elevation, dNBR, aspect, fuel loading (all size classes), moss and duff depth, and live tree density all influenced the distribution of understory species (Table 2). Certain variables, Kruskal-Wallis chi-square and P-value shown above each plot. Boxplots show median, hinges (25th and 75th percentiles), and whiskers that extend to values that are no more than 1.5*IQR (inter-quartile range) from upper or lower hinge, respectively. Points outside this range are outliers and plotted separately. Dunn's pairwise test significance are indicated by letters (α = 0.1) including 1-, 10-, and 100-h fuel loading and moss depth, were strongly associated with the first axis, while the remaining variables were not directly correlated with a single axis. The location of specific species in ordination space offers some insights on influences for species presence (Additional file 5); for example, shrub and hardwood species Betula neoalaskana Sarg., Populus tremuloides, and Populus balsamifera are all grouped to the far right along axis 1, indicating a correlation with high fuel loading and low moss and duff depth. However, the precise location of any individual species should be interpreted with caution given the high type I stress error in the final solution.

Overstory and understory 2005 to 2016
Neither density of live nor of dead trees changed significantly from one to 12 years post fire, indicating that the one-year post-fire sampling likely captured the majority of fire-related tree mortality. However, substantial but expected changes in cover of understory species occurred over time and a large number of species recruited between one and 12 years post fire: the total species  (2004), Alaska, USA. Kruskal-Wallis chi-square and P-value shown above each plot. Boxplots show median, hinges (25th and 75th percentiles), and whiskers that extend to values that are no more than 1.5*IQR (inter-quartile range) from upper or lower hinge, respectively. Points outside this range are outliers and plotted separately richness found on the 12 sites doubled. Bernhardt et al. (2011) also documented a sharp initial decline in the number of species post fire in this ecosystem, likely reflected in the low richness of our 2005 sampling, while our 2016 results highlight an influx of species into the burned sites over 11 years. This also reflects previous work showing that differing recruitment strategies (e.g., re-seeding from off-site, from seedbank, resprouting from underground structures) interacted with burn severity to result in both spatial and temporal changes in species richness with time since fire (Viereck 1983;Wang and Kemball 2005;Hollingsworth et al. 2013).
Individual species comparisons show a spectrum of responses. The increased cover of Vaccinium uliginosum, Ledum groenlandicum, and Salix spp. reflects trends found in other work, as does the decreased cover of Hylocomium splendens (Lecomte et al. 2005;Bernhardt et al. 2011). Bernhardt et al. (2011) found, in sites measured pre-and post-2004 wildfires in Interior Alaska, that Ledum groenlandicum, Vaccinium uliginosum, and V. vitis-idaea abundance significantly decreased one and two years post fire, which is reflected in the differences our study found between unburned and one-year results. However, our 12-year post-fire results show that, despite the initial decrease, L. groenlandicum, V. uliginosum, and Salix spp. have already recovered to or exceeded unburned average cover, likely due to their ability to resprout from underground rhizomes (Viereck 1983). On the other hand, cover of H. splendens at these sites continued to decrease from 2005 to 2016, reflecting continuing mortality of this climax-stage, shade-loving species (Tesky 1992). The cover of Chamerion angustifolium, an important off-site seeding species that may also resprout from rhizomes, generally decreased but remained above unburned levels. This is an unsurprising trend given that C. angustifolium is considered a diagnostic species for the first stage (1 to 20 years) of postfire recovery (Black and Bliss 1978).
Twelve year post-fire forest condition reflects dNBR class and site characteristics Density of live and of dead trees varied predictably with dNBR burn severity class, tracking with previous research that showed relatively strong relationships between overstory and remotely sensed indices Hudak et al. 2007). Spruce seedling density, on  (2004), Alaska, USA. Kruskal-Wallis chi-square and P-value shown above each plot. Boxplots show median, hinges (25th and 75th percentiles), and whiskers that extend to values that are no more than 1.5*IQR (inter-quartile range) from upper or lower hinge, respectively. Points outside this range are outliers and plotted separately. Dunn's pairwise tests are indicated by letters (α = 0.1) the other hand, was not significantly different among dNBR severity classes. Previous studies relating spruce regeneration to burn severity have alternately shown that Picea mariana is most successful on sites with some organic layer remaining (Johnstone and Chapin 2006b;Shenoy et al. 2011), that there is no relationship of seedling density to depth of organic layer burning (Gibson et al. 2016), or that regeneration success was highest in areas with exposed mineral soil and minimal organic cover (Mallik et al. 2010). The lack of trend in seedling density could reflect the lack of difference in moss and duff depth among dNBR classes at our sites, which may be traced back to the insensitivity of dNBR to organic layer consumption .
Given previous work relating field measures of surface burn severity to dNBR, we did not necessarily expect to see a strong trend in moss and duff depth with dNBR severity class (e.g., Hoy et al. 2008), and in fact there was no significant difference in organic layer depth among dNBR classes at our sites. Likely issues are the insensitivity of Landsat to height variation, and that the 30 m pixel resolution of Landsat is too coarse to capture the highly heterogeneous patterns of organic layer  (Alonzo et al. 2017). In addition, issues of topographic and canopy shadowing due to low solar elevation angles that affect dNBR and similar indices, even at lower latitudes, are compounded in high latitude areas such as Alaska . Overall, our results support the established understanding that Landsatbased dNBR reflects differences in overstory, but not necessarily understory, burn severity in this ecosystem.
On the other hand, downed woody fuel load did vary significantly among dNBR classes, with higher fuel loads in areas classified as low burn severity. In general, fire is expected to affect downed woody fuels by first consuming them during the active fire event (thus lowering fuel loading) but then by contributing fuels from mortality of trees and woody shrubs that become part of the surface fuel layer as they fall, as was evident upon revisiting the sites 11 years later. Higher fuel loads may raise the risk of more severe reburn effects but only at times when conditions are sufficiently dry to allow for fire to carry in the relatively high-moisture moss and duff ground layer ). However, we did not explicitly test the probability of reburn, and the canopy fuel load would still be much lower in these areas than in unburned stands (based on live tree  Johnston et al. 2015) without differentiating between burn severity levels, while others have focused only on finer fuel available for flaming consumption (e.g., Thompson et al. 2017) or on stands at least 30 years post fire (e.g., Hély et al. 2000), making comparison with our results difficult. Our fuel load estimates were similar to those measured in at least some similar stands of a study focused on time since fire (Johnston et al. 2015). However, our estimates of coarse woody debris (1000-h) were higher, likely because our sampling was concentrated in a timeframe when more dead trees had recently fallen to become surface fuel.
Species richness responded to both dNBR and woody fuel loading. Fuel load in the 100-h timelag class was negatively correlated with understory species richness and diversity indices. This contrasts with the findings of Table 1 Indicator species analysis Monte Carlo method (4999 permutations) results, where the observed indicator value (IV) is the percent of perfect indication based on relative abundance and relative frequency, the randomized IV is calculated from randomized groups, and P is the proportion of randomized trials with indicator value equal to or exceeding the observed indicator value. Only species significant at the α < 0.05 level are shown   Table 2 for correlations with environmental covariates and Additional file 5 for individual species Day et al. (2017) in a Pinus banksiana Lamb.-Picea mariana ecosystem, although Day et al. (2017) used only the presence or absence of woody debris rather than any measure of loading. Closer examination of the species composition at our low severity burned sites that may be driving this trend showed a higher understory cover of Betula neoalaskana and Populus tremuloides, which potentially increased fuel loading while shading out other understory species. This is also reflected in our NMDS results, in which those particular hardwood species were associated with higher woody fuel load. Total species richness was lowest for sites classified as high burn severity and highest for those classed as low burn severity, which mirrors previous work by Hollingsworth et al. (2013) in the area, as well as that of Lentile et al. (2007), who measured our same sites one year post fire, suggesting that the effects of burn severity on overall species richness lingers for at least a decade. Sites classified as moderate and high burn severity had different plant community composition than unburned areas 12 years after fire, with low and high burn severity areas also still differing from each other. The understory species community was influenced in part by surface factors, such as depth of moss and duff, as well as downed woody fuel loadings, which mirrors previous findings that bare ground and woody debris affect species composition (Day et al. 2017). Overstory variables like live tree density and basal area of dead trees, as well as factors like burn severity (dNBR), elevation, and aspect, were also significant drivers of understory vegetation composition. Whitman et al. (2018) found that similar factors, including overstory and understory tree density and basal area as well as soil moisture and heat-moisture index, were the dominant drivers of understory community dissimilarity, with burn severity having significant but secondary effects. Others have found that pre-fire forest type (Day et al. 2017) and burn severity (Hollingsworth et al. 2013) were primary drivers of species composition.
Several vascular and non-vascular plants were indicative of unburned areas. Empetrum nigrum L., Hylocomium splendens, and Pleurozium schreberi are known climax or late seral understory species that are maladapted to fire and may take decades to regain their pre-fire cover, making them an unsurprising indicator species for unburned sites (Lutz 1956;Johnson 1981;Tesky 1992). There is more limited information about post-fire response of Geocaulon lividum (Richardson) Fernald, but it has been shown to colonize early post fire and to increase in cover over time (Matthews 1994). We found no published information about the response of Flavocetraria cucullata (Bellardi) Karnefelt & A. Thell to fire. However, Peltigera aphthosa (L.) Willd. is generally killed by fire, is not a rapid recolonizer, and is often found in mature forests (Matthews 1993).
Chamerion angustifolium, a classic post-fire colonizer that often dominates high severity burns in Interior Alaskan boreal forests (Cormack 1953;Lutz 1956), was a significant indicator of high burn severity. We found no published data on the response of Vulpicida pinastri (Scop.) J.-E. Mattsson & M.J. Lai to fire specifically, although previous studies have theorized that certain lichens could serve as site pioneers Values in NMDS1 and NMDS2 (non-metric multidimensional scaling axis 1 and axis 2) are the direction cosines (i.e., the coordinates in ordination space for the head of the plotted vector arrow). R 2 correlations and P-values calculated based on permutation test Significance of P-value indicated as: *** < 0.001, ** < 0.01, * < 0.05, . < 0.1 (Ahti and Hepburn 1967), which could be why this species was a significant indicator of high severity burned sites. The soil crust lichen (unknown species) was unsurprisingly only found in high severity burned sites in this ecosystem, since these forests are characterized by a moss understory that often has no exposed soil. Morneau and Payette (1989) noted that a crustaceous lichen (Trapeliopsis granulosa [Hoffm.] Lumbsch), which was not specifically identified at our sites, was present and at maximum cover 14 years after fire in black spruce-dominated northern Quebec, Canada. Abiotic site factors, elevation and aspect, were significantly related to density of hardwoods, tall shrub cover, and downed woody fuels. Likely because the upper limit of elevation sampled (1081 m) was nearing treeline, there were negative correlations of elevation and hardwood sapling density, canopy cover, and 100-and 1000-h woody fuels. There was, however, a positive correlation between elevation and tall shrub cover, with the dominant shrubs being Salix, Alnus, and Betula species, reflecting the transition to a shrub-dominated landscape near treeline. Fine woody fuel loads (i.e., 1-, 10-, and 100-h classes) were higher on the warmer, south-southwest facing aspects than on the cooler north-northeast aspects. The warmer and dryer aspects may have had higher fine fuel loading because these areas are more likely to be dominated by Picea glauca or a mix of hardwood species, which have significantly higher average woody fuel loading than P. mariana stands that dominate cooler aspects (Pattison et al. 2018).

Conclusions
Species richness at remeasured sites increased by almost double the number of understory species over the 11 years between measurements. Cover of late-successional moss species (Hylocomium splendens) decreased over time, while cover of resprouting and early-successional shrub species (Salix spp., Vaccinium uliginosum, Ledum groenlandicum) increased. Twelve years post fire, remotely sensed burn severity classes (Landsat-derived dNBR) were related to long-term ecological condition in many of the ways that we expected based on previous research using ground-based severity estimates (e.g., Hollingsworth et al. 2013;Gibson et al. 2016;Day et al. 2017). Density of live trees decreased and shrub cover increased at higher burn severity classes, while downed woody fuel loading tended to peak in low and moderate burn severity classes. Understory plant community composition was still significantly different at all burn severity levels as compared to unburned, and at low burn severity as compared to high burn severity. As expected based on previous studies, plant composition was influenced by a variety of factors in addition to burn severity, including downed woody fuel loading and moss depth. Elevation and aspect also influenced understory species composition, in addition to influencing other factors such as hardwood sapling density, woody fuel loading, and shrub cover.
As wildfires increase in both size and number in boreal ecosystems, the need to estimate burn severity and its effects on a large scale will only increase in importance (Kasischke and Turetsky 2006;Mann et al. 2012;Calef et al. 2015;Young et al. 2017). While dNBR and the MTBS product have significant flaws, particularly in their ability to estimate depth of surface burning in boreal systems, these methods remain one of the primary sources for managers and others seeking to estimate burn severity. Efforts to improve upon the calculation of dNBR (Parks et al. 2018), and, for boreal systems specifically, to use LiDAR (Alonzo et al. 2017) or hyperspectral imagery (Lewis et al. 2011) to estimate organic layer consumption post fire offer promising avenues for applications of burn severity estimation that require more accurate data on surface burning (e.g., estimates of ecosystem carbon). However, the results of this study can help inform broad-scale estimates of post-fire ecosystem condition such as those created by LANDFIRE (Landscape Fire and Resource Management Planning Tools Project, https://www.landfire.gov/).

Additional files
Additional file 1: Monitoring Trends in Burn Severity (MTBS) metadata information for the burn severity classification of each fire from the Taylor Complex (2004), Alaska, USA, used in this analysis; all information is from the MTBS fire bundle download (www.mtbs.gov). Note that percentages of burn severity classes will not sum to 100 due to non-processing mask (scan line error, cloud cover, water, etc.) included in the product. (PDF 137 kb) Additional file 2: Common species (occurring in 50% or more of the 12 remeasured sites) encountered in 2005 and 2016 in sites sampled across the Taylor Complex fires (2004), Alaska, USA. Species are presented in alphabetical order first as the species common in both years and then those common in just one year (below the line). See Additional file 4 for common names and growth types. (PDF 162 kb) Additional file 3: Spearman rank correlation rho (ρ) and P-value for relationship between site species richness, evenness, and Shannon and Simpson indices of diversity on sites sampled in 2016 across the Taylor Complex fires (2004) review comments. The Tetlin Native Corporation provided access to some of the field sites. The authors thank two anonymous reviewers for their constructive and insightful comments.