Does burn severity affect plant community diversity and composition in mixed conifer forests of the United States Intermountain West one decade post fire?

Abstract

Background

Wildfire is an important ecological process in mixed conifer forests of the Intermountain West region of the USA. However, researchers and managers are concerned because climate warming has led to increased fire activity in recent decades. More area burned will result in larger land areas in early successional stages and will potentially limit tree establishment; therefore, evaluating long-term forest understory response to fire is important. We evaluate the impact of burn severity, overstory canopy cover, topography, and climate on understory plant community diversity and composition in seven wildfires, 9 to 12 years post fire, along a broad climate gradient from dry to moist mixed conifer forests in Idaho, Montana, and Washington, USA.

Results

Climate was the most important driver for species diversity and composition, but a burn severity gradient was detectable in the species data one decade post fire. A strong overlap in species composition between burn severity levels was documented, but dispersion was lower for high burn severity sites, indicating that those sites are still recovering. Local species richness and diversity had a nonlinear relationship with the burn severity index dNBR, with a maximum at low to moderate burn severity; the relationship was stronger in moist climates. Functional trait analysis revealed higher grass and forb cover in high-severity burns, higher cover of tree seedlings, residual and off-site colonizers in burned areas, and more shade-tolerant species in unburned areas. Of the 270 species recorded, 10% were introduced; however, only three were of noxious status and two were invasive annual grasses, generally occurring on dry sites.

Conclusions

The understory plant community was not fundamentally altered by these fires and fire contributed to increased species diversity both locally and regionally, suggesting that low to moderate burn severity fire is a treatment that contributes to long-term maintenance of a diverse and productive understory. Individual species traits were significant drivers of understory species assemblages and, as future change in climate and fire regimes leads to shifts in species composition, anticipation of consequences will be important. Although invasive species occurred at low cover levels, noxious weeds and invasive annual grasses will continue to be management challenges, particularly in dry regions of mixed conifer forests.

Resumen

Antecedentes

Los fuegos de vegetación representan procesos ecológicos importantes en los bosques mixtos de coníferas en la región intermontana del Oeste de los EEUU. Sin embargo, tanto los investigadores como los gestores de recursos están preocupados dado que el calentamiento global ha conducido a un incremento en la actividad de los incendios en décadas recientes. Más áreas quemadas van a resultar en mayores áreas en estadios sucesionales tempranos, y limitarán potencialmente el establecimiento de árboles. Es importante, por lo tanto, evaluar como el sotobosque responde al fuego en el largo plazo. Evaluamos el impacto de la severidad del fuego, la cobertura del dosel superior, la topografía, y el clima, en la diversidad y composición de la comunidad del sotobosque en siete incendios, 9 y 12 años después de ocurridos, a lo largo de un gradiente climático en bosques de coníferas, desde secos a húmedos, en Idaho, Montana y Washington, en los EEUU.

Resultados

El clima fue el conductor más importante en cuanto a diversidad y composición de especies, aunque un gradiente de severidad fue detectado en los datos de las especies una década después de los incendios. Una muy fuerte superposición fue documentada en la composición de especies entre niveles de severidad de los fuegos, aunque la dispersión fue menor en los sitios quemados a altas severidades, lo que indicó que esos sitios se están aún recuperando. La riqueza de especies locales y la diversidad tuvieron una relación no linear con el índice de severidad de quema dNBR, con un máximo en severidades bajas a moderadas; la relación fue más fuerte en climas más húmedos. El análisis de características funcionales reveló una más alta cobertura de pastos y malezas de hoja ancha en fuegos de alta intensidad, una más alta cobertura de plántulas de árboles, efectos de residuales, y colonización de exóticas en áreas quemadas, y más especies tolerantes a la sombra en áreas no quemadas. De las 270 especies registradas, 10% eran introducidas; sin embargo, solo tres tenían el estatus de nocivas y dos eran pastos anuales invasores, que generalmente aparecían en los sitios secos.

Conclusiones

La comunidad del sotobosque no fue fundamentalmente alterada por esos fuegos y el mismo fuego contribuyó a incrementar la diversidad de especies tanto a nivel local como regional, sugiriendo que los tratamientos de fuegos que queman de baja a moderada severidad son los que contribuyen al mantenimiento de un sotobosque diverso y productivo a largo plazo. Las características de las especies individuales fueron los conductores significativos de los ensambles de especies, y como el cambio climático y el régimen de fuegos llevan a variaciones en la composición de especies, es importante la anticipación de sus consecuencias. Aunque las especies invasoras aparecieron con bajos niveles de cobertura, las especies de malezas nocivas y los pastos anuales invasores continuarán siendo desafíos para el manejo, particularmente en regiones secas de bosques mixtos de coníferas.

Abbreviations

CC: overstory canopy cover

dNBR: Differenced Normalized Burn Ratio

fpp: length of frost-free period (days)

gsp: mean annual precipitation (map, mm), growing season precipitation (mm)

mat: mean annual temperature (10×°C)

mmax: mean maximum temperature of the warmest month (10×°C)

mmin: mean minimum temperature of the coldest month (10×°C)

MRPP: Multi-Response Permutation Procedure

NBR: Normalized Burn Ratio

NMS: Non-metric multidimensional scaling

NPMR: Non-parameteric multiplicative regression

NS: non-survivor (fire adaptation)

OC: off-site colonizer (fire adaptation)

RC: residual colonizer (fire adaptation)

smrp: summer precipitation (mm)

sprp: spring precipitation (mm)

SR: survivor rhizomes (fire adaptation)

SRCB: survivor taproot, caudex or bulb (fire adaptation)

Ƭ: Kendall’s tau coefficient

TRASP: transformed aspect where 0 is north-northeast and 1 is south-southwest

winp: winter precipitation (mm)

Background

Fire is an important and complex ecological process in conifer forests of the western United States (Agee 1993). Studies are projecting that larger fires in these and neighboring systems will occur more frequently (Barbero et al. 2015; Abatzoglou and Williams 2016; Bowman et al. 2017), impacting understory structure and cohorts of younger trees (Weiner et al. 2016; Smith et al. 2017; Sparks et al. 2017; Stevens-Rumann et al. 2018). Burn severity is usually defined via these and other understory and overstory characteristics that have changed due to fire and by the depth of soil combustion (Lentile et al. 2006; Lentile et al. 2007; Keeley 2009; Morgan et al. 2014). Burn severity is an important aspect of fire regimes, given its large influence on long-term forest structure (Goetz et al. 2007; Romme et al. 2011). Therefore, an improved characterization of the longer-term, second-order effects of fire on the understory in these conifer forests provides important descriptors of burn severity (Keeley 2009). This study sought to document understory characteristics one decade post fire in mixed conifer forests of the western US, across broad climate and burn severity gradients.

Estimates of burn severity across large spatial extents are usually achieved using the differenced Normalized Burn Ratio (dNBR; Key and Benson 2006) or relativized dNBR (RdNBR; Miller and Thode 2007), spectral indices applied to satellite sensor imagery. Spectral indices have been shown to correlate with the loss of biomass, duff consumption, and mortality of shrubs and trees induced immediately by fire, as well as with delayed mortality of trees within the first year post fire (Hudak et al. 2007; Lentile et al. 2007; Lentile et al. 2009; Strand et al. 2013; Miller et al. 2016; Weiner et al. 2016).

Relatively few studies have documented the relationships between burn severity and understory plant composition measures (such as forb, graminoid, and shrub cover) in mixed conifer forests of the western USA (Turner et al. 1999; Lentile et al. 2007; Perry et al. 2011; Burkle et al. 2015; Morgan et al. 2015; Romme et al. 2016; Willms et al. 2017). Forest understory plants are important contributors to the ecosystem because they provide most of the diversity in temperate forest plant communities (Roberts 2004; Gilliam 2007) and compose the early post-fire environment. Additionally, a diverse and productive post-fire plant community has been shown to contribute resistance to invasive species (e.g., McGlone et al. 2011) and constrain the abundance of invasives if established (Levine et al. 2004). The understory plant community further contributes to soil stability, thus preventing erosion (Roberts 2004); recycles nutrients (Gilliam 2007); provides cultural and traditional values (French 1965, Charnley et al. 2008, Wynecoop et al. 2019); and provides resources for wildlife and pollinators (Neill and Puettmann 2013).

Research in mixed conifer forests of the northwestern US indicates that climate factors and past fire occurrences are likely to impact forest regeneration in the future (Harvey et al. 2016; Kemp et al. 2016; Stevens-Rumann et al. 2018). Reduced natural tree regeneration has been documented in dry forests located at the edge of their climate tolerance zone (Donato et al. 2016) and in repeatedly burned mixed conifer forests (Stevens-Rumann and Morgan 2016; Stevens-Rumann et al. 2018). As climate- and fire-driven decreases in tree establishment, and conversion to non-forested systems, become greater management challenges in the next few decades, it will be increasingly important to understand the dynamics and resilience of understory vegetation in those ecosystems.

Plant mortality associated with fire depends on the heat received by the plant and the duration of the heat exposure. The plant’s growing points (buds or meristems) are particularly sensitive to heat and mortality depends on the amount of meristematic tissue injured (Brown and Smith 2000). Plants in fire-adapted ecosystems have developed a variety of strategies for surviving fire. Morgan et al. (2015) identified post-fire regeneration strategies in five categories: off-site colonizer; residual colonizer; survivor rhizomes; survivor taproot, caudex or bulb; and non-survivor. Off-site colonizers are plants that establish from a seed source outside the burned area, with seeds dispersed, for example, by wind or animals. Residual colonizers are plants that survive the burn (Stickney and Campbell 2000), for example, species with seeds that can survive deep in the soil for long periods of time (e.g., Ceanothus L.). Survivors can also be plants with tissue in the organic soil horizon or mineral soil that support dormant buds that generate new shoots on site following fire, for example rhizomes, caudex, or bulbs (Brown and Smith 2000; Stickney and Campbell 2000; Morgan et al. 2015).

Several measures of diversity have been developed for evaluating ecological differences and similarities between plant communities (McCune and Grace 2002; Legendre and Legendre 2012). Species richness refers to the number of species in a community, while diversity also accounts for the proportional abundance of those species, for example, Shannon’s diversity index (Shannon 1948). Diversity can be partitioned into alpha, beta, and gamma components, as originally proposed by Whittaker (1972). Alpha diversity reflects the local level diversity at sites within a habitat, while gamma diversity is the total species diversity in the region of interest. Beta diversity represents the dissimilarity in species composition between sampling units and can be calculated as gamma diversity divided by the mean alpha diversity (Whittaker 1972; Anderson et al. 2006; Legendre and Legendre 2012). Biotic cover or cover of live green vegetation is another measure commonly used to evaluate revegetation in post-fire environments (see, for example, Turner et al. 1999, Lentile et al. 2007, Edwards et al. 2015).

Several studies have applied these ecological measures to evaluate plant communities post fire. Generally, biotic cover and diversity are lower in high-severity burns (Turner et al. 1999; Lentile et al. 2007; Morgan et al. 2015), at least for a few years following a wildfire. Turner et al. (1999) found a greater number of forbs following crown fires compared to severe surface fires in conifer forests four years after the 1988 fires in Yellowstone National Park, USA. They further documented that, four years after fire, shrub cover was lower at increased burn severity. In a 25-year long-term study following the same Yellowstone fires, Romme et al. (2016) reported that species richness increased rapidly for the first five years post fire, and then levelled off to increase only slowly. Furthermore, Romme et al. (2016) found that most post-fire species had been present pre fire, surviving the fire via reproductive tissue below ground and rapidly sprouting post fire. Climate variables were the main control on the herbaceous understory immediately post fire; however, after 25 years, overstory canopy closure had become an important factor contributing to reduction of understory species, particularly those intolerant to shade.

In this study, we sought to improve our knowledge of long-term (decadal) post-fire understory vegetation response in mixed conifer forests of the Intermountain West region of the USA by characterizing plant species richness, diversity, cover, and community composition along broad burn-severity, topography, and climate gradients. To address this objective, field data were collected 9 to 12 years post fire in seven major wildfires that burned in mixed conifer forests in Montana, Idaho, and Washington, USA, during the time period 2003 to 2007.

One decade after wildfire, we centered our hypotheses around 1) species diversity, including alpha, beta, and gamma diversity (Whittaker 1972), and biotic green plant cover (hereafter, green cover); and 2) understory composition and associated environmental drivers.

Our diversity and green cover hypotheses were as follows:

• 1a) According to ecological theory (e.g., the Intermediate Disturbance Hypothesis; Huston 1979), moderate levels of disturbance are expected to increase local diversity. We therefore hypothesized that local species diversity (alpha diversity) is lower in areas burned at high severity compared to low to moderate burn severity about one decade post fire. We expected this relationship to be consistent along climate gradients.

• 1b) We expected beta diversity to be higher in areas burned at low severity because those areas contained both unburned and burned patches and therefore were more likely to support both species that survive the fire and post-fire colonizers.

• 1c) We expected wildfire to result in increased diversity across the sampled region (gamma diversity) because fire created patches with a variety of environmental characteristics suitable for species with different resource requirements (e.g., light).

• 1d) One decade after fire, we expected higher percentage of green cover in moderate- to high-severity burned areas compared to unburned because fire reduced overstory canopy cover, which allowed more light to reach the forest floor.

Our understory community composition and environmental gradient hypotheses were:

• 2a) We hypothesized that climate would be the main driver behind species composition, but that a burn-severity gradient would still be detectable one decade post fire.

• 2b) We expected statistically significant differences in species composition between unburned and high-severity burned areas, but not between unburned and low- to moderate-severity burned areas. We therefore expected unburned areas and areas burned at high severity to have a higher number of unique species compared to areas burned at low or moderate severity.

• 2c) We expected species functional traits (lifeform, lifecycle, origin, fire adaptation, and shade tolerance) to be significant drivers for plant community composition one decade post fire. We expected low- and moderate-severity burned areas to have a higher proportion of species that regenerate post fire via sprouting from underground reproductive tissue, and we expected a higher proportion of off-site and residual colonizers in high-severity burns. We expected shade tolerance to be an important factor determining species composition.

• 2d) Cover of noxious weeds and species with annual lifecycles would be highest at low elevation and in areas burned at high severity compared to unburned or other burn severity levels.

This study is unique because it documents longer-term changes in diversity and plant community composition after wildfires in mixed conifer forests across a broader climate gradient and larger region than most previous studies. It therefore contributes new knowledge to our understanding of how burn severity affects the understory component in these forests in the Intermountain West of the United States.

Methods

Study area

The seven wildfires incorporated in this study occurred in western Montana, central Idaho, and eastern Washington, USA (Fig. 1). The wildfires spanned a regional climate gradient ranging from dry to moist mixed conifer forest and included the Wedge Canyon and Robert fires west of Glacier National Park on the Flathead National Forest, the Black Mountain 2 and Cooney Ridge fires near Missoula on the Lolo National Forest, the School Fire in eastern Washington on the Umatilla National Forest, and the East Zone and Cascade fire complexes in central Idaho on the Payette and Boise national forests (Fig. 1, Table 1). The fires burned in 2003 to 2007 during July and August, which are the hottest and driest months of the year in the region (Pfister et al. 1977). Due to the varied topography, there are also steep temperature and precipitation gradients within each firescape, with some higher elevation areas receiving two to three times as much precipitation as some low elevation sites, with large aspect effects (Pfister et al. 1977). The fires exhibited all types of fire behavior, from crown fire to surface fire, and a range of fire intensities resulting in a patchy burn mosaic across the post-fire landscapes (Fig. 1, Table 2).

Fig. 1

Locations of wildfires in Montana, Idaho, and Washington, USA, that burned in July and August between 2003 and 2007, and their burn severity index classifications. The East Zone, Cascade, and Robert fire footprints contain gaps in the data caused by the failure of the scan line corrector (SLC) onboard Landsat 7; the other fires were imaged by Landsat 5

Table 1 Information about the fires included in this analysis of wildfire sites in Montana, Idaho, and Washington, USA, that burned in July and August between 2003 and 2007. The fire size was obtained from the final MTBS fire perimeter. The range of elevations, mean annual temperatures, and mean annual precipitation for sites sampled within fires is reported. Climate data was obtained from Rehfeldt (2006) and downscaled to the site level according to Rehfeldt et al. (2015)

Table 2 The number of hectares burned at different severity levels within the fire perimeter of each fire site in Montana, Idaho, and Washington, USA, that burned in July and August between 2003 and 2007, based on remotely sensed data. Percent area is reported in parentheses. Percentages do not always add up to 100% because of the sensor failure leading to striping and no-data returns for Landsat 7 images of the Cascade, East Zone, and Robert fires. The number of sites in each category are reported below the percent area in each class. Increased greenness is a remote sensing class identifying areas with more green vegetation after the fire than before the fire. These areas constitute less than 2% of any fire and were therefore not sampled

The dry end of the sampled conifer forests was primarily composed of the Douglas-fir (Pseudotsuga menziesii [Mirb.] Franco) habitat type (Pfister et al. 1977), dominated by Douglas-fir, grand fir (Abies grandis [Douglas ex D. Don] Lindl.), and ponderosa pine (Pinus ponderosa Lawson & C. Lawson), with denser Douglas-fir stands occurring in areas that experienced greater fire suppression in the last few decades before the wildfires (Cooper et al. 1991). At higher elevation sites, lodgepole pine (Pinus contorta Douglas ex Loudon), subalpine fir (Abies lasiocarpa [Hook.] Nutt.), and Engelmann spruce (Picea englemannii Parry ex Engelm.) were part of the conifer species mix. Western larch (Larix occidentalis Nutt.) commonly occurred at many sites across all fires. Historically, the dryer regions were thought to have a fire return interval of 15 to 30 years, with the fire return interval increasing at higher elevation. For example, the subalpine fir habitat type is thought to have a mean fire return interval as short as 90 years at lower elevation and as long as 150 years at higher elevation on northern aspects (Arno 1980). According to a fire atlas going back to 1889 (Gibson and Morgan 2005), the sampled areas had not burned within the past 63 years, and most of them had no fire recorded in the fire atlas, except six sites in the Idaho dataset that burned 18 years prior to the 2007 fires. A list of plants observed on three or more sampled sites can be found in the Additional file 1.

Field data collection

Field data were collected at a total of 89 sites across the seven fires. The Montana fires were sampled in 2013 to 2015 (10 to 12 years post fire), the School Fire in Washington was sampled in 2015 (10 years post fire), and the East Zone and Cascade Fire complexes in Idaho were sampled in 2016 (9 years post fire). Sample sites were stratified by elevation, cosine transformed aspect (trasp), and burn severity (unburned, low, moderate, high). Classified burn severity maps used for site stratification were produced by the Monitoring Trends in Burn Severity (MTBS) project (Eidenshink et al. 2007). Aspect was cosine transformed to avoid the circular aspect variable ranging from 0 to 360 degrees, where both 0 and 360 are assigned to north (Roberts and Cooper 1989). For site stratification and sampling purposes, two wildfires that burned in close proximity within the same forest type at the same time were treated as a single wildfire event; these include Black Mountain 2 and Cooney Ridge fires, Robert and Wedge Canyon fires, and East Zone and Cascade fires. Areas treated following fire (e.g., salvage logged, seeded, planted, mulched, etc.) were avoided to focus the effort on understanding natural regeneration and plant community development post fire.

All field sites were located at least 60 m but within 250 m of a road to avoid impacts related to roads and to optimize travel time. Each field site consisted of five 1 m² plots: a central plot surrounded by one plot upslope and one downslope from the center plot and a plot in each direction opposite one another on the slope contour. All outer plots were located 30 m from the center plot, such that each site was composed of five plots spanning four or five 30 m Landsat pixels. A minimum of 150 points were logged at each plot location with a Trimble GeoXT or 7X (https://www.trimble.com) global positioning system (GPS) unit with sub-meter precision at differential correction. At each of the five plots, understory cover of vascular plants (below 1.37 m in height) by species were ocularly estimated within a 1 m² quadrat, then averaged to the site level for analysis. Cover of moss (Bryophyta) was also recorded because it was occurring across the study area; moss was not recorded to the species level. Overstory canopy cover was recorded with a convex spherical densiometer in four cardinal directions; the four measurements were then averaged for each plot.

Image processing

The differenced Normalized Burn Ratio (dNBR; Key and Benson 2006) across the spatial extent of each wildfire was calculated from Landsat 5 TM or Landsat 7 ETM+ satellite sensor data (https://earthexplorer.usgs.gov). The dNBR was calculated as pre-fire NBR minus one year post-fire NBR; dNBR burn severity class thresholds followed Key and Benson (2006): unburned was classified as <100 dNBR values, low burn severity class from 100 to 269, moderate burn severity class from 270 to 439, and high burn severity class as >440. Images used in the analysis are listed in Table 3. For multivariate regression and correlation analyses, we used the average of the dNBR values of the four to five pixels where field data were collected.

Table 3 Landsat scenes and dates for images used to compute the dNBR index and burn severity classification for wildfire sites in Montana, Idaho, and Washington, USA, that burned in July and August between 2003 and 2007

Climate data

Climate variables included 30-year normals (1981 to 2010) for each site scaled to 30 m resolution. The downscaled data were obtained by applying a spline model (Hutchinson 2000) to ~1 km climate data developed for the western United States (Rehfeldt 2006; Rehfeldt et al. 2015) using a 30 m digital elevation model from the Shuttle Radar Topographic Mission (https://www2.jpl.nasa.gov/srtm). We evaluated the 20 climate variables for collinearity and removed variables that were correlated with another variable above 0.90 determined by the Pearson product moment correlation coefficient. The remaining nine downscaled climate variables included length of frost-free period (ffp, days), mean annual temperature (mat, 10×°C), mean maximum temperature of the warmest month (mmax, 10×°C), mean minimum temperature of the coldest month (mmin, 10×°C), mean annual precipitation (map, mm), growing season precipitation (gsp, mm), summer precipitation (smrp, mm), spring precipitation (sprp, mm), and winter precipitation (winp, mm). Soil properties, parent material, or other environmental gradients influencing plant communities were not considered in this study.

Data analysis

Species richness was computed as the number of species recorded within the five 1 m² plots at each site and included trees, shrubs, grasses, forbs, and moss. Species alpha diversity was calculated using Shannon’s diversity index (Shannon 1948, Legandre and Legandre 2012) for each site. Beta diversity was calculated for each burn severity level according to the formula

$$ {D}_{\upbeta}=\frac{D_{\upgamma}}{D_{\upalpha}} $$

(1)

(Whittaker 1972; Anderson et al. 2006), where D_β is the beta diversity; D_γ is the gamma diversity or the diversity of sites combined within each burn severity level; and D_α is the alpha diversity, which is the mean diversity of sites by burn severity level. We followed the suggestion by Whittaker (1972) and used true diversity (i.e., the exponential of Shannon’s index) for the calculation. Total green cover was calculated as the sum of cover of all plants and is therefore in some instances above 100% due to overlap of species within plots. Differences in species richness, alpha diversity, and green cover between burn severity levels (unburned, low, moderate, high) were tested with analysis of variance (ANOVA) with a post-hoc Tukey test for pairwise comparisons in the SYSTAT statistical software, version 13 (Systat Software Inc., San Jose, California, USA).

We used non-parametric multiplicative regression (NPMR) in HyperNiche 2.30 (McCune and Mefford 2009) to determine how overstory canopy cover, dNBR, and topographic and climatic variables interacted in non-linear, multiplicative ways (McCune 2006) to influence site-level response variables of understory species richness, alpha diversity, and green cover. For each NPMR analysis, we used local mean Gaussian model functions to identify combinations of predictor variables that maximized model fit and minimized overfitting. We controlled for overfitting by defining a minimum average neighborhood size of 5% of the samples (n = 89), a minimum data-to-predictor ratio of 10, and an improvement in fit criteria of 5%. Model fit was assessed using cross-validated R² (xR²). For each analysis, we evaluated the 200 models with highest xR² and determined the most commonly occurring predictor variables in the models. To further control for parsimony, we used an increase in xR² of 0.02 as a threshold to justify selection of a model with additional predictor variables. The model that included the most commonly occurring predictors, had the highest xR², and had the fewest predictors was selected as the best model. To test model stability, we ran a bootstrap resampling with 100 replacements to generate 100 datasets; we reported an average fit (bootstrap xR²) and standard error between the final model and 100 resampled datasets using the same model parameters. We also reported the average neighborhood size (N), which is the average number of sites contributing to the model estimate at each point on the modeled surface. We reported tolerance and sensitivity values for predictor variables in each selected model. The tolerance value indicates how broadly a given point estimate is influenced by the surrounding sample space. High tolerance values, relative to the range of the predictor, indicate that data points with a greater distance contribute to the estimate of the response variable. Sensitivity ranges from 0 to 1 and indicates the relative importance of each predictor in the model. A sensitivity of 1 indicates that a change in the predictor variable results in a change in the response variable of equal magnitude, while a sensitivity of 0 means that changing the value of the predictor has no effect on the response variable (McCune and Mefford 2009).

To analyze the large species cover dataset, multivariate ordination analysis was used to explore patterns in plant community composition, including gradients of dNBR, overstory canopy cover, elevation, cosine transformed aspect, and climate using the PC-ORD software (McCune and Grace 2002; McCune and Mefford 2011). Specifically, we used Nonmetric Multidimensional Scaling (NMS) with the Sørensen distance measure, a non-parametric ordination method designed to reveal gradients in species data and correlations with environmental variables. This method was chosen because it works well with non-normal data, does not assume a linear relationship, and avoids the zero truncation problem that many times occurs with community data (McCune and Grace 2002). We used the autopilot mode with the slow and through setting for the NMS runs, which assists in choosing the best solution. We removed species occurring two or fewer times in the dataset because they contributed little information and created challenges in calculating meaningful distances in the ordination (Peck 2010). Following ordination, axes were correlated with the species data and environmental variables using Kendall’s tau (Ƭ) coefficient (Kendall 1938). Kendall’s Ƭ ranges from −1 to 1, where 0 indicates lack of correlation. Explanatory variables included burn severity class, dNBR, elevation, cosine transformed aspect, overstory canopy cover, and climate variables. We tested the difference in understory species composition by fire and burn severity level using Multi-Response Permutation Procedure (MRPP) analysis with pairwise comparisons (McCune and Grace 2002). The MRPP generates two statistics, the P- and the A-values. The P-value represents the probability of a type I error under the null hypothesis of no difference between samples. The A-value is a measure of agreement between groups, where A = 1 for complete within-group homogeneity, A = 0 when the heterogeneity within groups is equal to the expectation, and A < 0 if there is less agreement within groups than expected by chance. As a measure of dispersion, we computed the average distance to the centroid (in ordination space) for each burn severity class (Anderson 2006) and differences in dispersion between burn severity classes were tested with ANOVA.

To identify potential indicator species within burn severity classes, we used the analysis method developed by Dufrêne and Legendre (1997) in PC-ORD (McCune and Grace 2002). The method evaluates the relative abundance and frequency of each species within groups (i.e., burn severity index levels in our analysis). To evaluate the probability that a species occurred according to expectation or exceeded expectation within a specific group, we ran 4999 Monte Carlo randomization tests in PC-ORD. We used α = 0.05 as the statistical significance threshold.

To evaluate differences in species functional traits, we classified the species by five traits including lifeform (tree, shrub, grass, forb, moss), lifecycle (annual, perennial), origin (native, introduced, noxious), fire adaptation (NS = non-survivor, OC = off-site colonizer, SR = survivor rhizomes, RC = residual colonizer, SRCB = survivor taproot, caudex or bulb; Morgan et al. 2015), and shade tolerance (high, moderate, low) (see Additional file 1). Species capable of regenerating from both seed and rhizomes (SR), or from taproot, caudex, or bulb (SRCB) were included in the SR or SRCB category. Differences in species cover by functional traits within burn severity groups were tested with ANOVA.

Finally, we used NPMR in HyperNiche 2.30 (McCune and Mefford 2009) to determine how overstory canopy cover, dNBR, and topographic and climatic variables influenced the presence of invasive annual grasses and species with a noxious status. We analyzed the influence of environmental drivers on cover of the introduced annual flammable cheatgrass (Bromus tectorum L.) and noxious weeds (USDA 2018) including spotted knapweed (Centaurea stoebe L.), Canada thistle (Cirsium arvense [L.] Scop.), and bull thistle (Cirsium vulgare [Savi] Ten.).

Results

Species diversity (hypothesis 1a–c)

Overall, 270 species were recorded at the 89 field sites, with 115, 167, 190, and 130 species recorded on unburned, low-, moderate-, and high-severity burn sites, respectively. Means and standard errors within burn severity group are reported in Table 4. Of those species, 26 were observed only on unburned sites, 26 only in low-severity burn sites, 35 only in moderate-severity burns, and 16 species were observed only in high-severity burns, whereas 62 species were recorded across all burn severity levels including unburned. It should be noted that, among the species that occurred in only one burn severity category, 84% occurred only once in the dataset, indicating that they were relatively rare species, in general. Overall, species richness was different between burn severity levels (df = 3, F = 3.693, P = 0.015), with pairwise comparisons confirming higher species richness at moderate burn severity compared to unburned (P = 0.016; Table 4); no other pairwise comparison was significant (P > 0.05). The overall ANOVA test for difference in alpha diversity between burn severity levels was significant (df = 3, F = 2.978, P = 0.036), with pairwise comparisons confirming higher alpha diversity at low burn severity compared to unburned (P = 0.044; Table 4).

Table 4 Summary of species richness and diversity measures by burn severity level for wildfire sites in Montana, Idaho, and Washington, USA, that burned in July and August between 2003 and 2007. Species richness denotes the number of species recorded at each severity level. The number of species unique to each severity class is also noted. Species diversity (Shannon’s diversity index) average is the average alpha diversity for sites within each severity class, SE denotes the standard error. Gamma diversity within a severity class is the Shannon’s index calculated for all sites within that severity class. Beta diversity represents the dissimilarity in species composition between sampling units. Letters (A, B, C) indicate statistically significant difference (P < 0.05) between burn severity classes for species diversity and richness. There was no difference in green cover between severity classes

Non-parametric multiplicative regression (NPMR) identified the predictor variables most influential on species richness and species diversity (alpha level). Richness and alpha diversity were most sensitive to climate variables; however, overstory canopy cover and dNBR, which both relate to burn severity, were also identified as predictors (Fig. 2a through f, Table 5). Species richness showed a non-linear relationship with dNBR with the highest richness occurring at intermediate levels of dNBR, indicating low to moderate burn severity (Fig. 2a, Table 5). Species richness generally increased with mean annual temperature as long as sufficient precipitation was available (Fig. 2b). Similar to species richness, alpha diversity also showed a non-linear relationship with dNBR such that higher alpha diversity was observed in low to moderate burns compared to unburned and high dNBR; the non-linear relationship was more pronounced at higher precipitation (Fig. 2c). Species alpha diversity further showed a non-linear relationship with both overstory canopy cover and dNBR, diversity being highest at intermediate levels of overstory canopy cover and dNBR (Fig. 2d, Table 5). Species alpha diversity generally increased with the length of the frost-free period as long as sufficient precipitation was available (graph not shown). Gamma diversity and beta diversity were both higher for sites that burned at low and moderate severity compared to high-severity burn sites or sites that did not burn in these fires (Table 4).

Fig. 2

Species richness (a and b), species alpha diversity (c and d), and understory green cover (e and f) graphed against predictor variables identified in NPMR analysis for wildfires in Montana, Idaho, and Washington, USA, that burned in July and August between 2003 and 2007. Both species richness and alpha diversity show a non-linear relationship with the burn severity index dNBR, with the highest richness and diversity at intermediate levels of burn severity (a through c). Abbreviations for predictor variables are: ffp = length of frost-free period (days), cc = canopy cover (%), sprp = spring precipitation (mm), dnbr = burn severity index, mat = mean annual temperature (10×°C), map = mean annual precipitation (mm). See Table 5 for model evaluations. The model surface is displayed in a red-to-black gradient: lighter colors indicate higher values and darker colors indicate lower values on the 3-dimensional surface. The vertical lines at the corners of some graphs are included to assist with reading the value on the y-axis and with seeing the 3-dimensional relationship. Areas colored in gray did not have enough data for model prediction

Table 5 Best predictor variables for species richness, alpha diversity, and green cover across climate, burn severity, and topographic gradients for wildfire sites in Montana, Idaho, and Washington, USA, that burned in July and August between 2003 and 2007. Abbreviations for predictor variables are: ffp = length of frost-free period, cc = canopy cover, sprp = spring precipitation, dNBR = burn severity index, mat = mean annual temperature, map = mean annual precipitation. xR² is the cross-validated R². Bootstrap xR² is the average fit of 100 models run with the selected predictors and associated standard error (SE). N is the average number of sites contributing to the model estimate at each point on the modeled surface. Tolerance is a smoothing parameter that indicates how broadly a given point estimate is influenced by the surrounding sample space. The tolerance value is reported in the original scale of the predictor variable. A low percent tolerance means that the response variable is sensitive to a change in the predictor variable. Sensitivity ranges from 0 to 1 and indicates the relative importance of each predictor in the model

Understory green cover (hypothesis 1d)

No difference was observed between burn severity levels for green cover (df = 3, F = 2.207, P = 0.093). Predictor variables most strongly influencing understory green cover were elevation, mean annual temperature, and the length of the frost-free period, with green cover decreasing with shorter frost-free period and higher elevation (Fig. 2e-f, Table 5). Neither the dNBR indicator of burn severity nor overstory canopy cover were selected predictors in the NPMR analysis, suggesting that one decade after fire, burn severity no longer influenced understory green cover in these wildfires.

Burn severity and species composition (hypothesis 2a and 2b)

After removing species with two or less occurrences from the 89 sites across all fires, 143 species remained in the dataset prepared for ordination analysis (see species list in Additional file 1). Species composition significantly differed between some of the sampled fires according to an MRPP test (P < 0.001, A = 0.069). A multiple comparisons test suggested significant differences in species assemblage between fires that were geographically farther apart. The Montana fires were different from Idaho fires (P < 0.001, A = 0.036) and the Washington fire (P < 0.001, A = 0.061), and the Idaho fires were also different than the Washington fire (P < 0.001, A = 0.050). We therefore organized the species dataset by state and ran one ordination for the fires in each state.

Similarities between the ordinations of fires within the three states were observed. Ordinations for all three states resulted in solutions for which the ordination axis with the highest explanatory power was related to climate variables and the second axis was related to burn severity variables. Strong overlap was observed between burn severity levels but high-severity burn sites were generally different from unburned sites. The convex hull polygon for high burn severity represented smaller region in species space compared to other burn severity levels, indicating lower dispersion in the high-severity burn class.

The ordination for the Montana fires (41 sites) yielded a three-axis solution that explained 71.1% of the variation in the data and a stress factor of 16.6. The first ordination axis explained 40.8% of the variation in the plant community composition data, while the second and third axis explained 14.8% and 15.5%, respectively. Vectors of explanatory variables for the axes are displayed in Fig. 3a and b. Correlation (Kendall’s Ƭ) with the explanatory variables (Table 6) showed that climate variables related to precipitation and temperature were the most important variables explaining axis 1, for which a high value on axis 1 corresponded to lower mean annual temperatures and high precipitation annually and during all seasons. Axis 2 was positively correlated with overstory canopy cover and negatively correlated with dNBR, representing the burn-severity gradient. Axis 3 was most strongly correlated with dNBR, overstory canopy cover, and elevation (Fig. 3a and b, Table 6). Strong overlap was observed between the burn severity levels, and it was also noted that the polygon for high burn severity occupied a smaller region in species space compared to low and moderate burn severity and areas not burned in these fires. Analysis of the distance to centroid within burn severity groups confirmed differences in dispersion between burn severity classes (df = 3, F = 2.945, P = 0.045; Fig. 4a). Pairwise comparisons did not reveal differences at P < 0.05, but it can be noted that lower dispersion was observed for high burn severity compared to unburned (P = 0.052) and moderate burn severity (P = 0.060). MRPP analysis indicated significant differences between burn severity groups (P < 0.001, A = 0.0371), although the low A value suggests that the differences may not be ecologically significant (McCune and Grace 2002). A pairwise comparisons test suggested significant differences only between unburned and high burn severity (P < 0.001, A = 0.0258).

Fig. 3

Ordination for fires that burned in Montana (a and b), Idaho (c and d), and Washington (e and f), USA, in July and August between 2003 and 2007. Convex hulls for each burn severity group and vectors with |Ƭ| > 0.200 are shown as black lines. Axis 1 for Montana and Washington fires and Axis 2 for Idaho fires represent a temperature–precipitation gradient, while Axis 2 for Montana and Washington and Axis 1 for Idaho represent a burn severity gradient. The axes represent ordination scores and do not have units. Abbreviations for predictor variables are: elevation = elevation, dNBR = burn severity index, trasp = transformed aspect, cc = canopy cover, ffp = length of frost-free period, gsp = growing season precipitation, map = mean annual precipitation, mat = mean annual temperature, mmax = mean maximum temperature of the warmest month, mmin = mean minimum temperature of the coldest month, smrp = summer precipitation, sprp=spring precipitation, and winp = winter precipitation

Table 6 Correlation of the explanatory variables to the ordination axes derived from species data by state (Kendall’s tau statistic, Ƭ) for wildfire sites in Montana, Idaho, and Washington, USA, that burned in July and August between 2003 and 2007. Explanatory variables are: elevation, index of burn severity (dNBR), transformed aspect (trasp), overstory canopy cover (cc), length of frost-free period (ffp), mean annual temperature (mat), mean maximum temperature of the warmest month (mmax), mean minimum temperature of the coldest month (mmin), mean annual precipitation (map), growing season precipitation (gsp), summer precipitation (smrp), spring precipitation (sprp), and winter precipitation (winp). Correlations in bold (Ƭ > 0.200) are shown as vectors in Fig. 3

Fig. 4

Box plots of the distance to ordination centroid within burn severity classes in Montana fires (a), Idaho fires (b) and the Washington fire (c), USA, that burned in July and August between 2003 and 2007. The center line represents the sample median, the box represents the first and third quartile of the sample, the whiskers represent the minimum and maximum after outliers have been removed. Outliers are defined as data points located farther away than >1.5 times the distance between the first and third quartile (the length of the box). Outliers are noted with asterisks (*). In each graph, the box plots represent, from left to right, unburned, low-, moderate-, and high-severity burn sites. P-values for ANOVA tests for differences by burn severity class are listed. Letters (A, B, C) indicate statistically significant difference (P < 0.05) between burn severity classes. Pairwise comparisons between high burn severity and unburned (P = 0.052) and moderate burn severity and unburned (P = 0.060) for Montana fires are shown in (a)

The ordination for the Idaho fires (32 sites) yielded a three-axis solution that explained 71.0% of the variation in the data and a stress factor of 15.1. The first ordination axis explained 24.5% of the variation in the plant community composition data while the second and third axes explained 28.0% and 18.6%, respectively. Vectors of explanatory variables for the axes are displayed in Fig. 3c and d. Correlation (Kendall’s Ƭ) with the explanatory variables (Table 6) showed that climate variables related to precipitation and temperature were the most important variables explaining axis 2, the axis with the highest explanatory power, for which a high value on axis 2 corresponded to high elevation, lower temperatures, a shorter growing period, and high precipitation annually and during the growing season. Axis 1 was positively correlated with transformed aspect and negatively correlated with overstory canopy cover. Axis 3 was most strongly correlated with summer precipitation (smrp; Fig. 3c and d, Table 6). Similar to the Montana fires, strong overlap could be observed between the burn severity levels, and the polygon for high burn severity covered a smaller region in species space. Analysis of the distance to centroid confirmed no difference in dispersion between burn severity classes (df = 3, F = 0.496, P = 0.688; Fig. 4b). MRPP analysis did not show any significant differences between burn severity groups (P = 0.496, A = −0.001).

The ordination for the Washington fire (16 sites) yielded a three-axis solution that explained 82.4% of the variation in the data and a stress factor of 9.9. The first ordination axis explained 39.0% of the variation in the plant community composition data while the second and third axes explained 29.3% and 14.1%, respectively. Vectors of explanatory variables for the axes are displayed in Fig. 3e and f). Similar to the Montana and Idaho fires, the correlation (Kendall’s Ƭ; Table 6) showed that climate variables related to precipitation and temperature were the most important variables explaining axis 1, for which a high value on axis 1 corresponded to lower mean annual temperatures and high precipitation annually and during all seasons. Axis 2 was positively correlated with dNBR and negatively correlated with overstory canopy cover, representing the burn-severity gradient. Axis 3 was most strongly correlated with overstory canopy cover and dNBR (Fig. 3e and f, Table 6). Also for the Washington fire, overlap could be observed between the burn severity levels, and the polygon for high burn severity represented a smaller region in species space. Note that there were not enough low burn severity sites to display a polygon. Analysis of the distance to centroid within burn severity groups confirmed differences in dispersion between burn severity classes (df = 3, F = 3.825, P = 0.039; Fig. 4c). Post-hoc pairwise comparisons revealed lower dispersion for high burn severity sites compared to moderate burn severity (P = 0.039). MRPP analysis did not suggest any differences in composition between burn severity groups (P = 0.2355, A = 0.020), although an ad-hoc pairwise comparisons test suggested significant differences between unburned and high burn severity sites (P = 0.008, A = 0.112).

Indicator species analysis computed indicator values (IV) and identified species indicative (P < 0.05) of the burn severity levels in the large species dataset. Wild strawberry (Fragaria virginiana Duchesne; IV = 28.5, P = 0.024), hawkweed (Hieracium L. spp.; IV = 25.4, P = 0.024), Sitka valerian (Valeriana sitchensis Bong.; IV = 16.7, P = 0.018), and starry Solomon’s seal (Smilacina stellata [L.] Link; IV = 15.3, P = 0.045) were identified as indicator species for low burn severity sites. Western larch (IV = 21.0, P = 0.023) was identified for moderate burn severity sites. Fireweed (Chamerion angustifolium [L.] Holub; IV = 39.2, P = 0.006), common dandelion (Taraxacum officiale F.H. Wigg.; IV = 24.9, P = 0.017), and snowbrush (Ceanothus velutinus Douglas ex Hook.; IV = 15.0, P = 0.026) were identified as indicators for high burn severity sites. Indicator species identified in areas that did not burn in the sampled fires were sidebells wintergreen (Orthilia secunda [L.] House; IV = 23.5, P < 0.001), rattlesnake plantain (Goodyera oblongifolia Raf.; IV = 19.7, P = 0.013), and fernleaf biscuitroot (Lomatium dissectum [Nutt.] Mathias & Constance; IV = 16.2, P = 0.015). Species with canopy cover among the top 10 in a burn severity class and state are listed in Additional file 2; note the many species present in all burn severity levels 10 years post fire.

Functional traits analysis (hypothesis 2c)

Functional trait analysis revealed a few differences between burn severity groups by lifeform, lifecycle, origin, fire adaptation, and shade tolerance (Fig. 5) 9 to 12 years post fire. For the combination of lifeform, lifecycle, and origin, perennial grass and annual forb cover increased with burn severity in Montana and Idaho but was variable in Washington (Fig. 5a, b and c). Cover of tree seedlings increased with burn severity in Montana (Fig. 5a). The proportion of residual and off-site colonizers increased with burn severity in Montana (Fig. 5d) and differences in rhizome survivors were significant in Washington (Fig. 5f). The cover of shade-intolerant species increased with burn severity level in Montana (Fig. 5g) and the cover of moderately shade-tolerant species was higher at low and moderate burn severity in Washington (Fig. 5i).

Fig. 5

Percent understory canopy cover by species traits and burn severity class for sites in Montana, Idaho, and Washington, USA, that burned in July and August between 2003 and 2007. In each graph, the bars represent, from left to right, unburned (U), low- (L), moderate- (M), and high- (H) severity burn sites. Canopy cover of species with a combination of life form, lifecycle, and origin are shown for Montana (a), Idaho (b), and Washington (c). The cover of species with different fire adaptive traits are shown for Montana (d), Idaho (e), and Washington (f); and cover of species in shade tolerance classes are shown for Montana (g), Idaho (h), and Washington (i). P-values for ANOVA tests for differences by burn severity class are listed. Letters (A, B, C) indicate statistically significant difference (P < 0.05) between burn severity classes

Introduced and noxious species (hypothesis 2d)

Introduced species were present at low cover levels within the studied fires. Introduced annual grasses were observed on 13 of the lower-elevation sites (900 to 1500 m) with mean cover values ranging from 0.1 to 1.2% across burn severity classes, with one low burn severity site having as high as 12% cheatgrass cover. Annual grass species were dominated by cheatgrass but also included field brome (Bromus arvensis L.) and rattlesnake brome (Bromus briziformis Fisch. & C.A. Mey.). Annual grasses were absent on higher-elevation sites (>1500 m), except for one low burn severity site with 0.2% cheatgrass. Evaluation of cheatgrass cover in relationship to dNBR, canopy cover, topographic and climate variables with NPMR (Table 7) showed that areas at risk of increased cheatgrass cover are characterized by spring precipitation in the 100 to120 mm range (Fig. 6a) and a mean minimum temperature of the coldest month above −10 °C (Fig. 6b). Cheatgrass was more likely to occur in areas classified at dNBR <600 (Fig. 6a and b), corresponding to low and moderate burn severity rather than high.

Table 7 Best predictor variables for invasive species occurring in more than 10% of the wildfire sites in Montana, Idaho, and Washington, USA, that burned in July and August between 2003 and 2007. Abbreviations for predictor variables are: sprp = spring precipitation, dNBR = burn severity index, mmin = mean minimum temperature of the coldest month, map = mean annual precipitation, and winp = winter precipitation. xR² is the cross-validated R². Bootstrap xR² is the average fit of 100 models run with the selected predictors and associated standard error (SE). N is the average number of sites contributing to the model estimate at each point on the modeled surface. The tolerance value indicates how broadly a given point estimate is influenced by the surrounding sample space. Sensitivity ranges from 0 to 1 and indicates the relative importance of each predictor in the model

Fig. 6

Bromus tectorum cover (a and b), Centaurea stoebe cover (c and d), and Cirsium spp. cover (e) graphed against predictor variables identified in NPMR analysis for sites in Montana, Idaho, and Washington, USA, that burned in July and August between 2003 and 2007. Abbreviations for predictor variables are: sprp = spring precipitation (mm), winp = winter precipitation (mm), dnbr = burn burn severity index, mmin = mean minimum temperature of the coldest month (10×°C), map = mean annual precipitation (mm). See Table 7 for model evaluations. The model surface is displayed in a red-to-black gradient: lighter colors indicate higher values and darker colors indicate lower values on the 3-dimensional surface. The vertical lines at the corners of some graphs are included to assist with reading the value on the y-axis and with seeing the 3-dimensional relationship. Areas colored in gray did not have enough data for model prediction

The relatively recently introduced wiregrass (Ventenata dubia [Leers] Coss.; Wallace et al. 2015) occurred on three sites on the School Fire in Washington at cover levels ranging from 0.2 to 6%. Because of the low sample size, we did not relate V. dubia occurrences to environmental variables. Other observed introduced perennial grasses, introduced as pasture grasses to North America, were Kentucky bluegrass (Poa pratensis L.), orchardgrass (Dactylis glomerata L.), and timothy (Phleum pratense L.).

Noxious invasive species observed on sampled sites included spotted knapweed, which occurred on eight sites in Montana fires across burn severity levels, with cover ranging from 0.2 to 9.2%. Introduced thistle species were observed on 12 low- to high-severity burn sites at cover levels at or below 6% and included Canada thistle and bull thistle. Spotted knapweed and thistle species were also evaluated against environmental variables with NPMR (Table 7, Fig. 6c through e). Although the models were not very strong, most likely due to low presence, results suggested that both of these noxious invasives are more likely to occur in the warm and dry region of the mixed conifer forest in areas that burned at high severity. Other observed introduced annual forbs were mouse-ear chickweed (Cerastium L. spp.), yellow lucerne (Medicago falcata L.), and starwort (Stellaria L.spp.), and introduced perennial forbs included bird’s-foot trefoil (Lotus corniculatus L.), common sheep sorrel (Rumex acetosella L.), common dandelion (Taraxacum officiale), common mullen (Verbascum thapsus L.), common plantain (Plantago major L.), clover (Trifolium L. spp.), oxeye daisy (Leucanthemum vulgare Lam.), prickly lettuce (Lactuca serriola L.), yellow salsify (Tragopogon dubius Scop.), and veronica (Veronica L. spp.). Introduced shrubs that occurred on three sites were red raspberry (Rubus idaeus L.) and cutleaf blackberry (Rubus laciniatus Willd.).

Discussion

Species diversity along a burn-severity gradient (hypothesis 1a through c)

Although climate variables were the most important predictors for species richness and alpha diversity, burn severity was a contributing factor. In agreement with our first hypothesis (1a), we observed lower local richness and alpha diversity one decade post fire in areas burned at high severity. Analysis of richness and alpha diversity along gradients of predictor variables revealed non-linear relationships between the burn severity index (dNBR) and species richness and diversity, for which the highest diversity was observed at a low to moderate burn severity. We interpreted this observation in light of the intermediate disturbance hypothesis (IDH; Huston 1979, Huston 2014). Huston (1979) stated that local species diversity is optimized at intermediate levels along disturbance gradients. Several researchers have since applied the IDH to various aspects of disturbance regimes (e.g., intensity, extent, and duration; see review by Shea et al. 2004). Our research suggested that the IDH may also apply to the burn-severity gradient within a fire. In the case of burn severity levels observed one decade post fire, the mechanism for the observed humpback-shaped relationship is likely short-term prevention of competitive exclusion rather than long-term coexistence as suggested in the original description of the IDH (Connell 1978; Huston 1979; Huston 2014).

For species richness, this humpback-shaped relationship with the burn severity index spanned the entire moisture gradient (Fig. 2a); however, for alpha diversity, the relationship was more pronounced in the moist region of the mixed conifer forest (Fig. 2c). We had not expected differences in the diversity–severity relationship along the climate gradient. An explanation could be that there was a larger proportion of fire-adapted species at the drier end of the climate gradient where fires occur more frequently. Given that moist mixed conifer forest burns at a lower frequency than dry forest, the accumulation of fuels and the conditions required to sustain a wildfire may result in higher plant mortality and loss of seedbank, contributing to a slower post-fire recovery, particularly in areas burned at high severity.

In agreement with our hypotheses 1b and c, both beta and gamma diversity were higher at low and moderate burn severity compared to high burn severity areas or areas that did not burn in the sampled fires. Burkle et al. (2015) also documented higher local species richness, beta diversity, and species dispersion at mixed-severity burns compared to high-severity burns, confirming that fire disturbance contributes to increased diversity in fire-adapted systems. Furthermore, patches of variable burn severity create post-fire environments that are suitable for establishment of different species, thereby increasing the overall richness and diversity at a landscape perspective (i.e., gamma diversity). Although we did not study the spatial pattern of patches burned at different severity levels, our results support the hypothesis proposed by Halofsky et al. (2011), stating that mixed burn severity levels within fires create landscapes where early- and late- successional plant communities are present, creating a landscape that potentially contributes to diversity also in wildlife species. One recent example demonstrating that diversity in burn severity can cascade into diversity of avian species is given by Tingley et al. (2016), with greater effects observed 10 years post fire compared to one year post fire.

Green cover along a burn-severity gradient (hypothesis 1d)

According to hypothesis 1d, we expected that lower overstory canopy cover in moderate-to high- severity burns would result in a higher percentage of green cover on the forest floor compared to unburned or low burn severity areas. However, we found that, one decade post fire, understory green cover was more related to elevation and climate than to burn severity measures such as dNBR or as indicated by overstory canopy cover. This result does not support our hypothesis 1d. It is possible that other variables that we did not asses introduced variability in understory green cover. For example, we did not record soil characteristics or geology, and we did not have a good measure of forest floor consumption and soil combustion following fire (e.g., Wang and Kemball 2005), which could introduce fine-scale variability in burn severity with consequential impacts on understory green cover. Our results suggest that one decade is enough time for understory plants to return to pre-fire cover levels, which is in agreement with other studies stating that understory green plant cover is back to pre-fire cover levels one decade after fire. Following the Yellowstone fires in 1988, Turner et al. (1999) found that understory biotic cover in areas burned in surface fires was not detectably different from adjacent unburned areas four years later; while understory biotic cover in areas burned in a crown fire was only half of that in unburned areas. Edwards et al. (2015) did not find significant differences in understory green cover between burn severity classes in sub-boreal conifer forest in British Columbia, Canada, five to six years after an experimental burn. A remote sensing study by Lewis et al. (2017), on some of the same Montana fires as in this study, suggested that, 10 years post fire, green cover had increased to approximately 60% regardless of initial condition, but that this took about twice as long on high burn severity sites. Remote sensing studies have also shown that Net Primary Productivity (NPP) in North American boreal forests has returned to pre-fire levels nine years following fire (Hicke et al. 2003).

Climate, burn severity, and species composition (hypothesis 2a and b)

In agreement with our hypothesis 2a, understory species ordination and correlation with explanatory variables revealed a strong environmental temperature–precipitation gradient along the axis with the highest explanatory power (axis 1 for Montana and Washington and axis 2 for Idaho). This is not surprising since the fires in both Montana and Idaho spanned a large elevation gradient. Information about how species assemblages change along temperature, precipitation, and burn-severity gradients will be particularly important in the future given anticipated changes in climate. Regional climate predictions forecast increased temperature and water pressure deficit across the western continental US, combined with increased area burned (Barbero et al. 2015; Abatzoglou and Williams 2016). The axis with second highest explanatory power (axis 2 for the Montana and Washington fires and axis 1 for the Idaho fire) was correlated with dNBR for the Montana and Washington fires, which we interpret as the burn-severity gradient. Similar to Romme et al.’s (2016) study in Yellowstone, we found that the burn-severity gradient was detectable in the understory species composition 10 years post fire. Romme et al. (2016) further found that, 25 years post fire, the burn-severity gradient was no longer detectable. Given the state of post-fire recovery observed in our data, we will likely be able to draw the same conclusion.

According to our hypotheses 2b, we suggested that differences in species composition would be detectable between unburned and high burn severity areas, but not between unburned and low to moderate burn severity areas. A common trend for all three ordinations, and confirmed with MRPP tests, was the strong overlap in species space between unburned areas and the different burn severity levels, indicating similar species composition across burn severity levels. Species composition was statistically different between high-severity burns and unburned areas in Montana and Washington; however, the MRPP A-value was low, suggesting weak ecological relevance. Our results are comparable to those reported by Romme et al. (2016), who concluded that the 1988 Yellowstone fires did not fundamentally alter the understory plant community, and that most post-fire species had been present pre fire. One decade post fire, we observed a smaller polygon in species space and significantly lower dispersion around the centroid for sites that burned at high severity compared to the other burn severity levels, similar to observations by Burkle et al. (2015). In our study, indicator species confirmed that few species were unique to any burn severity level, particularly moderate burn severity. We conclude that our hypothesis 2b was supported by our findings; however, the differences in species composition on high burn severity sites compared to other burn severity levels or unburned areas lie predominantly in lower dispersion and proportions of species rather than different species. For both the Montana and Washington fires, the high burn severity polygon was located toward the moist end of the climate gradient, which can be expected given a longer fire return interval with a mixed or stand-replacing fire regime (Arno 1980) in higher-elevation mixed conifer forests.

Species functional traits (hypothesis 2c)

One decade post fire, differences in the functional traits analyzed (lifeform, lifecycle, origin, fire adaptation, and shade tolerance) were still detectable in some of the fires. As expected according to our hypothesis 2c, higher proportion of off-site and residual colonizers were documented in high-severity burns, and also identified via indicator species analysis (e.g., fireweed and dandelion [off-site colonizers] and snowbrush ceanothus [residual colonizer]). Snowbrush ceanothus seed can remain viable in the soil for a very long time and germinate following fire (Agee 1993). Weiner et al. (2016) also documented snowbrush ceanothus cover in high-severity burns in western juniper (Juniperus occidentalis Hook.) woodlands, and Dodge (2018) documented ceanothus in high-severity burns in ponderosa pine. We demonstrated that species functional traits and responses to fire harbor important information that can help forest stewards anticipate consequences of changes in climate and fire regimes; Johnstone et al. (2016) called this ecological memory.

In reviewing the list of most common species in each burn severity class (Additional file 2), we noted that a large proportion of species classified as non-survivors were mosses. Moss was present at all burn severity levels in all fires, often at higher cover in burned areas, even severely burned (Additional file 2). Additional research about the role of mosses in the post-fire environment of mixed conifer forest in the US is warranted. In future studies, we recommend that mosses are identified to the species level because species differ in their fire adaptation. For example, moss species likely to be found within the study region range from fire adapted (e.g., Ceratodon purpureus [Hedw.] Brid. and Polytrichum juniperinum Hedw.), to surface rhizome adapted (e.g., Lycopodium annotinum L.), to mosses that do not survive fire or recolonize quickly post fire (e.g., Pleurozium schreberi [Brid.] Mitt. and Selaginella densa Rydb.).

Although this study was not designed to analyze tree regeneration, it is worth noting that higher cover of tree seedlings was recorded in burned areas compared to unburned. Tree regeneration is of concern particularly in the warm and dry climate region of the mixed conifer forest due to a warming climate (Harvey et al. 2016; Kemp et al. 2016), and due to frequently repeated fires (Stevens-Rumann et al. 2018).

Noxious weeds and annual grasses (hypothesis 2d)

Noxious weeds and annual invasive grasses were not widespread in sampled areas; however, their presence deserves attention. We hypothesized (2d) that noxious weeds and annual grasses would be highest at low elevation and in areas burned at high severity. Noxious weeds and annual grasses indeed had a higher presence in the lower elevation dry mixed conifer forest; however, cheatgrass was more likely found at lower burn severity while noxious weeds were positively associated with higher burn severity. Cheatgrass was present in all fires except in the colder regions of Montana (Wedge Canyon and Robert fires). Most occurrences and those with highest cover (up to 12%) occurred on the School Fire. Cheatgrass, well known for its flammable properties and ability to carry fire quickly, has recently been shown to alter fire regimes even at cover levels below 5% (Bradley et al. 2018). It has contributed to an increase in fire frequency in sagebrush steppe communities throughout the western US, leading to loss of diversity and native species (Svejcar et al. 2018). Wiregrass, observed on the School Fire, is becoming increasingly invasive in perennial grass systems such as prairies, pastures, lands in conservation reserve programs, and open forests (Wallace et al. 2015; Fryer 2017). Similar to cheatgrass, wiregrass is a non-native, invasive, annual grass that reduces productivity on the lands it invades (Fryer 2017). It is different from cheatgrass because it can inhabit moister and higher elevation sites; it has the ability to expand into and replace established, seemingly resilient, grasslands; and it has low nutritional value (Fryer 2017). Spotted knapweed was found at lower elevation in Montana in the Robert, Black Mountain 2, and Cooney Ridge fires. It can outcompete native vegetation via its deep taproot, its high seed production and quick seed spread, and its low palatability to grazing animals (Zouhar 2011).

Conclusions

While climate variables were most important for influencing understory species diversity and composition among variables tested, fire was an important ecological process for maintaining species diversity both locally and regionally in the mixed conifer forest of the US Intermountain West. Within the fires that we sampled, the mixed conifer forests appeared resilient even to high-severity fire. However, it is taking longer for the understory in high burn severity areas to recover; evidence of fire in the form of lower understory diversity and dispersion was still apparent after one decade. We did not find evidence for divergence in plant community assemblages or differences in successional trajectory for areas that burned with high, moderate, or low severity in the fires studied. Although we documented some differences in plant community composition in high burn severity compared to low burn severity or unburned areas, those species appear to be a subset of the plant species assemblage in the study area. Given time, plant species in the unburned, low, and moderate burn severity sites will likely continue to expand into areas burned at high severity. We suggest that the areas we sampled were in “safe operating space,” as defined by the forest resilience framework proposed by Johnstone et al. (2016), with regard to the forest understory. Areas with a potential to develop “resilience debt” in understory components were those with relatively high cover of invasive species, particularly cheatgrass and spotted knapweed, but also wiregrass. “Resilience debt” is broadly defined by Johnstone et al. (2016) as loss of an ecosystem’s capacity to recover due to misalignment between disturbance regimes or post-disturbance conditions and community adaptation to fire; for example, changes in ecological communities that affect the response to future disturbance.

Forests at the lower-elevation montane ecotone that, according to previous research, are susceptible to loss of tree canopy due to drought or frequent reburns (Donato et al. 2016; Stevens-Rumann and Morgan 2016; Stevens-Rumann et al. 2018) are likely also the most susceptible to invasive community-altering forbs and grasses. Maintaining a resilient understory will be particularly important in areas where loss of tree canopy can be expected, to manage the land for functional and productive plant communities in a changing climate. Understanding individual species functional traits and responses to climate and fire will be essential for anticipating plant community adaptation to a changing climate and fire regimes.