The Effect of Seeding Treatments and Climate on Fire Regimes in Wyoming Sagebrush Steppe

Abstract

Wildfire size and frequency have increased in the western United States since the 1950s, but it is unclear how seeding treatments have altered fire regimes in arid steppe systems. We analyzed how the number of fires since 1955 and the fire return interval and frequency between 1995 and 2015 responded to seeding treatments, anthropogenic features, and abiotic landscape variables in Wyoming big sagebrush ecosystems. Arid sites had more fires than mesic sites and fire return intervals were shortest on locations first treated between 1975 and 2000. Sites drill seeded before the most recent fire had fewer, less frequent fires with longer fire return intervals (15–20 years) than aerially seeded sites (intervals of 5–8 years). The response of fire regime variables at unseeded sites fell between those of aerial and drill seeding. Increased moisture availability resulted in decreased fire frequency between 1994 and 2014 and the total number of fires since 1955 on sites with unseeded and aerially pre-fire seeding, but fire regimes did not change when drill seeded. Greater annual grass biomass likely contributed to frequent fires in the arid region. In Wyoming big sagebrush steppe, drill seeding treatments reduced wildfire risk relative to aerial seeded or unseeded sites.

1. Introduction

Fire regimes have changed considerably since the 1980s with increases in wildfire size and frequency across multiple ecosystems [1,2]. Climate change has led to shorter winters, decreased snowfall, shifts in precipitation regimes, and extended periods of drought, contributing to longer periods of fire risk [3,4,5]. Historical overgrazing and seed dispersal by livestock in arid lands and forests have contributed to annual grass invasions [6,7], and anthropogenic features, such as fuel breaks [7], roads, and utility lines [8], act as dispersal corridors. Global change has shifted many native plant communities to invasive-dominated communities, which promote fire [9,10]. Thus, climate change and annual grass cover increases have contributed to the increased numbers of large wildfires in the western United States [11,12], especially during periods of severe weather, such as dry lightning and strong winds [4,13]. Interactions among topography, climate change, fire history, invasive plants, and the difficulty in identifying vegetation with remote imagery provide uncertainty [14]; models predict increases in fire number, size, and fire season length [3,15,16,17] in the western United States. Fuel treatments may become insufficient to reduce future wildfire size, particularly during dry, windy conditions [4].

Sagebrush-dominated communities in the western United States once covered an estimated 38 million to 109 million ha [18]. Land use has converted 20% of these systems into private land, and grazing has impacted nearly all of the remaining sagebrush [19]. Invasive species also threaten the sagebrush steppe by altering the historical fire regimes [12,20,21]. In the United States’ Great Basin, fire regimes in sagebrush steppe ecosystems vary based on dominant sagebrush species. Within a sagebrush species, estimates of fire rotation or fire return vary due to the ranges of elevation found and regional climate variation [22]. For instance, one study estimated fire rotation for Artemisia arbuscula communities was 325–450 years [23] but could range from 93 to 187 in the Columbia Plateau and >1000 years in more xeric sites [22]. Artemisia tridentata subsp. vaseyana communities have a fire rotation of 70–217 years [23,24], but the lower elevation Artemisia tridentata subsp. wyomingensis ecosystems have fire rotations of 100-342 years [23,24]. Historical fires in Artemisia tridentata subsp. wyomingensis communities were larger, more continuous, stand-replacing fires than in Artemisia tridentata subsp. vaseyana communities [24].

Shifts in climate have favored a range expansion of the winter annual grass, Bromus tectorum (cheatgrass), an invasive species in the Intermountain West [25]. Consequently, B. tectorum has altered the historical fire regime in the Great Basin, leading to profound changes in fire size, frequency, and duration [12,21,26]. After a fire or other disturbance, B. tectorum takes advantage of the newly cleared space to recruit rapidly, increasing density and cover [10,27,28,29]. Even low B. tectorum cover will increase the chance of adjacent, non-invaded habitat burning in subsequent fires [26]. Models suggest climate change will increase the risk of B. tectorum invasion in Idaho, Montana, and Wyoming [25], making B. tectorum cover reduction by creating resistant plant communities an important goal of many fuels and post-fire rehabilitation treatments in the region [23,30].

Seeding treatments are commonly applied in the Great Basin to increase forage production, improve wildlife habitat, stabilize sites, and rehabilitate burned areas. Seed may be applied aerially from aircraft to avoid soil disturbance but may drop seeds in unfavorable habitat and decreases soil-seed contact limiting the likelihood of successful plant recruitment [31]. Drill seeding uses modified tractors to till the soil while simultaneously sowing seed [31]. Tilling soil may create microhabitats that promote desirable species recruitment, but the disturbance may also enhance B. tectorum recruitment [7,32]. Managers also combine aerial and drill seeding to enhance diversity and improve recruitment likelihood [31,33].

The plant species used in aerial and drill seeding treatments have changed since the mid-twentieth century. In the 1940s, nonnative bunchgrasses, typically Agropyron desertorum or A. cristatum, were drill seeded in monoculture or low diversity mixes to increase cattle forage throughout the sagebrush steppe [18,34,35], but aerial seeding was used infrequently [36,37]. In the mid-1980s, the United States Bureau of Land Management (BLM) switched to using the “Hycrest” crested wheatgrass cultivar (Agropyron desertorum X A. cristatum) in rangeland seeding treatments [36,38], which created larger plants with more blades and greater fecundity than either parent species [39]. During the 1970s in the Great Basin, the number and size of fires increased and half the seeding treatments were conducted after wildfires as part of rehabilitation efforts [36]. After 1990, the BLM conducted drill and aerial seeding treatments predominately after wildfires in the Great Basin [36]. In the nineties, the most common species used in rangeland seeding treatments were either Hycrest or Agropyron fragile (Siberian wheatgrass) [36]. The proportion of seeding treatments using native grasses, forbs, and Artemisia tridentata was limited before 1990 [36] but increased considerably after 2000 [36,40] in response to Executive Order 13112. The order required federal land management to use native rather than nonnative species, when possible. After 2000, the BLM also increased the agency’s use of native grasses and forbs in seeding treatments. [41]. Before 1990, the number of drill seeding treatments in the Great Basin was four times greater than aerial seeding treatments [36]. In the 1990s, the number of aerial and drill seedings were roughly equal [36]. Since 2000, managers used aerial seeding twice as often as drill seeding treatments [36]. Since 1980, the proportion of treatments using Artemisia species has doubled, and treatments including native grasses have quadrupled in the Great Basin [36]. The use of native species has created habitats more suitable for sagebrush obligate species, such as Centrocercus species (sage grouse), Brachylagus idahoensis (pygmy rabbit), and Artemisiospiza species (sparrows) [42,43,44,45,46] compared to sites dominated by nonnative species [47].

Increasing wildfire number and size in the sagebrush steppe requires a better understanding of the impact of environmental variables and land management practices on natural recovery, annual grass invasion, and rehabilitation efforts after fire [48]. For example, fuel breaks, green strips, and prescribed fires create disturbances that can provide suitable habitat for B. tectorum to establish [7,49] and invade into an intact habitat. Disturbances along roads and power line corridors provide suitable habitat and corridors for B. tectorum dispersal [8]. Studies suggest that livestock grazing can facilitate B. tectorum invasion by transmitting seeds and clearing or damaging existing vegetation through preferential grazing [50,51]. Increasing B. tectorum would lead to more frequent annual grass fires resulting in further shifts in the fire regime [12,21]. Few studies have examined the effects of seeding treatment practices on fire regimes in rangelands [52,53]. Recent studies have examined the effect of post-fire seeding treatments in the sagebrush steppe at large scales [40,47], but those efforts focus on sites with recovery after a single fire. It is unclear what effect multiple fires and management actions over decades might have on invasive species, plant communities, and/or shifts in fire regime attributes.

We examined the relationships among fire history, seeding treatments, and environmental characteristics to test the effect of those predictor variables on the response of fire regime variables. The fire regime variables we used as response variables included (1) the number of fires between 1955 and 2015, (2) the mean fire return interval between 1995 and 2015 (time between fires), and (3) the fire frequency between 1995 and 2015 (number of fires per unit time) [54]. Though the fire return interval and fire frequency are correlated, there are distinctions. Fire return interval, or mean fire return interval, is the average number of years between fires [55], which indicates the time available for species to recruit and establish communities between fires [36,37,56,57]. Fire frequency is the recurrence of fire in a given area over time [55] and is used in modeling [12,58,59,60,61] because it reflects the probability that a fire will occur in any one year. Our objective was to understand how climate, environmental characteristics, proximity to roads and private land (potential dispersal corridors of B. tectorum and sources of ignitions), and seeding treatments affect the fire regime characteristics; namely the total number of fires since 1955, as well as fire return interval and fire frequency between 1995 and 2015. We hypothesized that seeded sites would have (1) fewer fires, (2) longer mean fire return intervals, and (3) less frequent fires.

2. Materials and Methods

2.1. Site Description

The study area was 209,000 ha, located in southern Idaho, U.S.A., south of the Snake River (42.5° N, −115.4° W, Figure 1). Elevation in this region ranges from 750 to 1800 m increasing in elevation southward. Precipitation and snowfall increased as elevation increased (Figure 1). Climate data from Remote Automated Weather Stations (RAWS) in the study region show an increase in temperatures from 1962 to 2012; the timing and amount of rain and snow have also changed (Figure 2). Snow at all elevations has decreased since 1962. Precipitation has increased during the spring months across all elevations and decreased during summer months at near 800 m and in June near 1600 m.

The dominant shrub in the study area was Artemisia tridentata subsp. wyomingensis Beetle and Young (Wyoming big sagebrush). Artemisia thinning and removal from 1940 to 1980 to increase cattle forage and wildfires 1980-present have impacted the once-dominant vegetation, though some isolated patches remain [62]. The most common grass species were Poa secunda J. Presl (Sandberg bluegrass, hereafter Poa), Agropyron cristatum (L.) Gaertn (crested wheatgrass, hereafter Agropyron), and Bromus tectorum L. (cheatgrass or downy brome, hereafter Bromus). Other native bunchgrasses included Pseudoroegneria spicata (Pursh) Á. Löve (bluebunch wheatgrass), AchnatherumAchnatherum thurberianum (Piper) Barkworth (Thurber’s needlegrass, hereafter Achnatherum), and Elymus elymoides (Raf.) Swezey (squirreltail, hereafter Elymus). Poa cover was >15% at 70% of our sites and co-dominated with either Agropyron or Bromus. The historical fire rotation for Artemisia tridentata subsp. wyomingensis is estimated between 100 and 342 years, depending on the method used to estimate fire return [23,24]. Fires have burned most of the study area at least once, while other regions have burned up to seven times in the past 65 years [63].

Seeding treatment data was acquired from the Land Treatment Digital Library [38] and the Bureau of Land Management (BLM) [64,65]. Seeding treatments included aerial, drill, combined seeding, and unseeded treatments. Seeded species varied by treatment and year (see Section 1) [38,65]. As the number and size of wildfires increase in the 1980s [63], the BLM stopped thinning Artemisia stands and focused both aerial and drill seeding treatments on post-fire rehabilitation [38,64]. Seeded species for drill seeding treatments before 2000 were primarily low diversity seed mixes dominated by nonnative Agropyron species. After 2000, the proportion of seeding treatments using native grasses exceeded the use of nonnative grasses. Seeding treatments rarely included Artemisia tridentata prior to 1990, but 70% of seeding treatments included Artemisia after 2000. In our study area, aerial seeding was used primarily to seed shrubs species, although a small percentage of aerial seed mixes also included grasses and forbs when seeded on Wilderness Management Areas or rocky sites.

2.2. Data Extraction

We used ArcMap 10.3 software (ESRI) to generate 4000 random points across our study site with a minimum of 1 km distance between points. We eliminated points on features BLM land managers would never treat such as roads and highways [66], water features (e.g.- rivers, ponds, reservoirs) [67], and land owned not owned by the agency [68] by reviewing data extracted using the Intersection tool. We then removed points not located in pixels designated as “Wyoming Big Sagebrush-Wheatgrass” in “Group Name” category of LANDFIRE Environmental Site Potential Biophysical Settings raster [69]. The points were then verified using 1 m resolution orthorectified aerial imagery with a 6 m horizontal accuracy [70]. The result was 573 points (Figure 1a) used to extract data from polygon and raster layers. The response variables (fire number between 1955 and 2015, mean fire return interval 1995–2015, and mean fire frequency from 1995 to 2015) were derived from BLM fire polygon shapefiles [63]. Prescribed fires in the sagebrush steppe are much less common than wildfire, but we included prescribed fires when calculating the number of fires, if they occurred. The increased number and size of wildfires in the sagebrush steppe makes understanding the impact of environmental variables, seeding treatments, and potential sources of anthropogenic ignition sources on natural recovery and rehabilitation efforts after a fire imperative [48]. To address these concerns, we used information on seeding treatment history [64,65,71], climate [72], elevation [73], and proximity to anthropogenic features on the landscape [66,68] as explanatory variables. For a complete list of the 99 explanatory variables, see Table S1. Data extraction included unburned and burned sites. Unburned locations were typically unseeded or drill seeded, while burned areas could be unseeded or seeded via aerial, drill, or both methods.

2.3. Field Collection

Using 67 field sites from a concurrent field study conducted in 2014 and 2015, we collected standing herbaceous and woody biomass along the environmental gradient (Figure 1b). There were 11 unburned sites and 56 sites that burned one, two, three, or six times. We divided sites into the most recent treatment type (drill seeded, aerially seeded, or unseeded sites) nested within the number of fires that had occurred since 1950. The project did not include locations with both aerial and drill seeding conducted at the same time. Biomass was collected from six randomly placed 1 m² quadrats within a 180 m² location and classified into one of four functional groups: perennial bunchgrass, annual grass, forb, and shrubs. Biomass was air-dried in a closet for over a year.

2.4. Statistical Analysis

According to historical fire data from the BLM [63], there were few fires in our study area between 1900 and 1950. We defined fire number as the number of times a site burned in between 1955 and 2015. To determine how past vegetation management, fire history, and environmental conditions shaped the current fire regime characteristics, we defined fire return interval as the average number of years between fires from 1995 to 2015. We assigned sites that had burned only once during the 20-year period a fire return interval of 20 years. We labeled places that did not burn between 1995 and 2015 a fire return interval of 25 years. Fire frequency was determined by dividing the number of fires on a site over the twenty years from 1995 to 2015 by twenty.

We used a Nonparametric Multiplicative Regression (NPMR) in HyperNiche 2.3 [74] to determine how climate, seeding treatments, fire history, and proximity to anthropogenic features predicted the total number of fires since 1955, as well as the fire return interval and fire frequency between 1994 and 2014 (Table S1). The NPMR. evaluates how explanatory variables interact in nonlinear and multiplicative ways to alter the dependent variable and allows for the detection of nonlinear but significant relationships among explanatory and dependent variables. While Generalized Additive Models, Generalized Linear Models, and Multiple Linear Regression can analyze the relationship between response and explanatory variables, they can only handle a few potential explanatory variables. The accuracy of those analyses is limited to linear relationships. NPMR analyzes relationships between more variables and yields models that reflect nonlinear relationships [75].

Our NPMR used a quantitative local mean Gaussian weighting model. We assessed model fit using cross-validated R² (xR², Equation (1)), which a measure of the relationship between the residual sum of squares (RSS) and the total sum of squares (TSS) such that:

$$\mathrm{x}{\mathrm{R}}^2=1-\frac{\mathrm{R}\mathrm{S}\mathrm{S}}{\mathrm{T}\mathrm{S}\mathrm{S}}$$

(1)

The xR² uses a ‘‘leave-one-out’’ cross-validation and does not require withholding data for validation purposes [76]. To control overfitting, we set the model improvement criteria to a data-to-predictor ratio of ten, and a minimum of a 2% improvement was required to increase the number of variables in a model by an additional variable. Bootstrap resampling (each dataset resampled with replacement 100 times to generate 100 new datasets, each with n – 1 plots) was used to assess model stability against the inclusion of particular plots in a given analysis by providing an average xR² (±SE) between the final model and 100 resampled datasets. We selected the models with the highest xR² value, average neighborhood size (>10), and an improvement criterion of 2% [76].

Unlike traditional regression, NPMR does not fit coefficients in an equation. Instead, NPMR reports tolerances—the standard deviations used in the Gaussian smoothing. High tolerance values, relative to the predictor range, indicate a greater distance among points targeted for estimation. We report the average neighborhood size and sensitivity for each model. Neighborhood size is the average number of sample units contributing to the estimate of occupancy at each point on the modeled surface. Sensitivity indicates the relative importance of each quantitative predictor in the model. A sensitivity of 1 means that, on average, the percent change in the value of a predictor will result in a similar percentage change in the estimate of the response variable. In contrast, a sensitivity of 0 indicates the predictor does not affect the response variable. Since tolerance and sensitivity indicate how much the dependent variable changes in response to a predictor variable, sensitivity can only be calculated for continuous variables.

We examined the effect of 99 explanatory variables (Table S1) on fire regime characteristics (fire number, fire return interval, and fire frequency). We evaluated the models with the 100 greatest xR² values, and we selected two to three models based on the predictor variables that appeared most frequently. We analyzed the sensitivity of the predictor variables and selected the model with the greatest sensitivity.

In the Great Basin, fire frequency and return interval are affected by the types and amount of available fuels, particularly annual grasses such as B. tectorum. We analyzed fuels composition (i.e., the types of fuels) and amount by collecting biomass in four functional groups (perennial bunchgrass, annual grass, forb, and shrubs). We analyzed the composition of functional group biomass using a Multi-Response Permutation Procedure (MRPP), a multivariate approach that analyzes the relationship among dependent variables (the biomass of four functional groups) and determines if there are differences among groups while allowing for an unbalanced design [77]. The MRPP compares the observed data to a randomized data. The MRPP measures the degree of separation among predetermined groups using the test statistic, T, which is the difference between the observed and expected weight mean within-group distance divided by the variance of the expected mean within group distances. A greater negative value of T indicates greater separation among groups. In MRPP, the A-statistic is the effect size of within group similarity compared to the expected similarity based on the randomized data. When A > 0 groups are more homogenous than chance by chance; A < 0 indicates groups are more heterogeneous than expected. The MRPP includes a pairwise comparison of individual groups but does not correct for multiple comparisons automatically. We used the Benjamini and Yekutieli false discovery rate (FDR) method described by Narum [78] to establish a conservative threshold for significance for the multiple comparison test (α = 0.0115). We conducted three MRPP analyses using predetermined groups based on (1) seeding treatment (aerial, drill, or unseeded), (2) four 250 m elevation ranges: low (770–1020 m), medium low (1021–1270 m), medium high (1271–1520 m), and high (1521–1780 m) and (3) seeding treatment nested within elevation bands [79]. Though the NPMR. data includes sites that received both aerial and drill seeding in the same treatment, the 67 sites used for acquiring biomass data used in the MRPP did not include such treatments. There were 11 unburned sites and 56 sites that burned one, two, three, or six times. Sites were divided into the most recent treatment type (drill seeded, aerially seeded, or unseeded sites) nested within the number of fires that had occurred since 1950. If one variable’s values are orders of magnitude greater it can skew the analysis, yielding results based primarily on that one variable. In our data, woody biomass was much greater than the biomass of the other three groups. To normalize data, biomass was square root transformed. We used a Jaccard distance measure to further decrease the influence of shrub biomass.

We used a Nonmetric Multidimensional Scaling (NMS) ordination with a Jaccard distance matrix to visualize relationships among MRPP groups (e.g., elevation, seeding treatment) within the ordination space determined by biomass. Drill seeded sites at 770–1020 m and 1521–1780 m had only one representative site each and were dropped from the analysis. Stress was determined by comparing observed data to randomized data 50 times using a stability criterion of 0.0001. We used an analysis of variance in MyStat 12 (Systat Software Inc) to evaluate if the biomass in one or more plant functional groups were different among elevation and treatment groups. Shapiro–Wilks test determined if the data was normally distributed and a Levene’s test was used to test for homoscedasticity. When biomass variables did not meet the assumptions for a parametric ANOVA, we used the Kruskal–Wallis Test.

3. Results

The number of fires since 1955 ranged from unburned to seven fires (

Table 1. Summary statistics of the fire regime response variables. The total fire number was the sum of the fires between 1955 and 2015. The mean fire return interval and frequency were based on fire data between 1995 and 2015. Sites that did not burn between 1995 and 2015 were designated as a fire return interval of more than 25 years and had a fire frequency of 0.

Fire Regime Variable	Category	Count (n)
Total Fire (n)	Unburned	164
	1	114
	2	154
	3	63
	4	42
	5	24
	6–7	12
Return Interval (years)	<5	6
	5–10	143
	11–20	188
	>20	236
Fire Frequency	0.00	236
	0.05	1880
	0.10	81
	0.15	43
	≥0.20	26

). The fire return intervals ranged from 3 to 20 years, with 217 of the 573 locations not burning after 1995 (Table 1). Fire frequency between 1995 and 2014 varied from 0 to 0.3 with 150 sites burning with a frequency of 0.1 or greater (Table 1). Treatment history varied among locations with the total number of aerial seedings on a site ranging from 0 to 5 treatments and drill seedings ranging from 0 to 3 treatments (Table S2). The treatment history included 260 sites (45.3%) with some history of aerial seeding, but only 83 sites (14.5%) had two or more aerial seeding treatments since 1950. Sites were more likely to have a history of drill seeding treatments in our dataset. There were 309 sites (54%) with at least one drill seeding and 65 locations (11%) had more than one drill seeding treatment. Though there was a tendency for the number of aerial and drill seedings to increase with the number of fires (Table S3), the weighted average of aerial seeding treatments did not exceed two until there were more than six fires, and the weighted average of drill seeding treatments did not exceed one until the seventh fire (Figure S1).

The NPMR model that best predicted the number of fires between 1955 and 2015 consisted of the first year a site was burned, minimum spring vapor pressure deficit (VPD), total number of seeding treatments, and most recent seeding treatment method (

Table 2. The Nonparametric Multiplicative Regression best-fit models for each fire regime variable. The explanatory variables that compose the models include the first year a site burned (1st fire year), the 30-year mean minimum vapor pressure deficit from March-May (Spring VPD min), the seeding treatment before the most recent fire (Pre-fire seeding treat), the total number of times a site was aerially seeded (Aerial seeding (n)), the total number of times a site was drill seeded (Drill seeding (n)) and 30-year mean spring precipitation from March-May (Spring precipitation).

Response Variable	Eval xR²	Model Fit xR² (±SE)	Avg Size	Predictor Variable	Sensitivity	Tolerance
Fires (n) *	0.8	0.83	31.2	1st fire year	0.09	9.45	(15%)
1950–2014		(±0.0)		Spring VPD min	0.05	0.26	(25%)
				Treatment (n)	0.05	1.5	(25%)
				Pre-fire seeding treat †	-	-	-
Fire Return Interval *	0.78	0.81	43.32	Aerial seeding (n)	0.22	0.5	(10%)
1994–2014		(±0.0)		1st year treated	0.15	10.35	(15%)
				Spring VPD min	0.02	0.46	(45%)
				Pre-fire seeding treat †	-	-	-
Fire Frequency *	0.76	0.77	33.29	Drill seeding (n)	0.07	0.75	(25%)
1994–2014		(±0.0)		Spring precipitation	0.05	6.28	(25%)
				Aerial seeding (n)	0.03	0.25	(5%)
				Pre-fire seeding treat †	-	-	-

Eval xR2 = The correlation from original model analysis. Model Fit xR2 = Correlation based on the bootstrap analysis. Avg. Size = The number of sample units contributing to the estimate of occupancy. * p < 0.01, SE = Standard error. † Categorical variable. Sensitivity and tolerance can only be calculated with continuous variables.

). The number of fires was unaffected by vapor pressure deficit (VPD) for sites that first burned before 1990 (Figure 3a). For sites that burned for the first time after 1990, arid sites (greater spring VPD) tended to have more fires than more mesic sites. At sites that burned for the first time in 2014, the most arid sites burned 2.5 times on average while more mesic sites burned only once (Figure 3a). The number of fires increased with the number of seeding treatments regardless of the first year a site burned (Figure 3b). Fire number increased with greater numbers of treatments, particularly in areas with more arid spring weather; six treatments contributed to two fires at the most mesic sites and nearly four fires at more arid sites (Figure 3c).

The seeding treatment prior to the most recent fire (pre-fire) method had a significant effect on the modeled number of fires since 1955 (Figure 4). Sites that were drill seeded or aerial and drill seeded prior to the most recent fire tended to only burn twice; unseeded sites burned three to four times and aerial seeded sites burned three to five times (Figure 4a). Increasing aridity resulted in increased fire number for locations that were aerially seeded or unseeded prior to the most recent fire, but fire number at sites that were drill seeded, aerial, and drill seeded, or had no prior history of seeding treatments did not increase with aridity (Figure 4a). Fire number decreased as the year a site burned became more recent for aerially seeded treatments, but the number of fires for other pre-fire treatments were relatively unaffected by the year a site first burned (Figure 4b). Regardless of the number of historic treatments on a site, aerially seeded sites burned more times than other treatments (Figure 4c). On sites with one treatment, if the site was aerially seeded or left unseeded after the penultimate fire, there were three fires compared to two when drill seeded (Figure 4c). On sites with one treatment, drill seeded sites burned twice, while aerially seeded and unseeded sites burned three times (Figure 4c).

The NPMR model that best predicted the fire return interval between 1995 and 2015 consisted of four variables: the number of aerial seeding treatments, the first year a site was treated, the minimum spring VPD, and the pre-fire seeding method (Table 2). When the number of aerial seeding treatments increased, the fire return interval decreased across all levels of spring VPD but was reduced at high VPD (Figure 5a). On sites without aerial seeding treatments, the fire return interval was shortest on sites that first burned between 1975 and 1990 (Figure 5b). Adding aerial seeding treatments resulted in decreased fire return intervals, particularly at locations first treated before 1970 (Figure 5b). The fire return interval was shorter when on more arid sites regardless of the first time a site was treated (Figure 5c). Sites first treated in 1975–1995 had the shortest fire return interval on the most arid sites, while the shortest fire return interval on the more mesic sites were when the site was treated between 1985 and 1995 (Figure 5c).

As the number of aerial seeding treatments in a site’s treatment history increased, fire return intervals tended to decrease but varied according to seeding treatment type preceding the most recent fire (Figure 6a). The fire return interval was 20 years with drill seeding regardless of whether or not there was aerial seeding before the drill seeding (Figure 6a). For sites with both aerial and drill seeded prior to the most recent fire, the fire return interval decreased from 15 to 16 years with one aerial seeding treatment to 14 years with a history of two aerial seeding treatments (Figure 6a). The fire return interval decreased from 16 to 19 years when there was no aerial seeding to as little as 7 years when there was a history of two aerial seeding treatments (Figure 6a). When the most recent seeding treatment was aerial seeding, the fire return interval decreased from as much as 10 to as little as 5 years as the number of aerial seeding treatments increased (Figure 6a). On sites with a history of one aerial seeding, sites that were aerially seeded after the penultimate fire had the shortest fire return intervals (Figure 6a). Drill seeding in addition to aerial seeding increased the fire return interval, although fire return intervals were greater on sites that were only drill seeded (Figure 6a). For aerial seeded sites, the fire return interval increased as the first treatment year became more recent; however, fire return intervals decreased slightly when aerial and drill seeding were combined as first treatment year became more recent (Figure 6b). The fire return interval was unchanged for sites with drill seeding treatments regardless of when the site was first treated (Figure 6b) or across the spectrum of average minimum spring vapor pressure deficit (Figure 6c). Fire return interval decreased as the minimum spring VPD increased on unseeded or aerially seeded sites after the penultimate fire (Figure 6c). Fire return interval on drill only or aerial and drill seeding were unaffected by the increase in spring VPD (Figure 6c).

The NPMR model that best predicted fire frequency from 1995 to 2015 consisted of four variables: the number of drill seeding treatments, spring precipitation, the number of aerial seeding treatments, and the pre-fire seeding method (Table 2). Fire frequency decreased slightly as precipitation increased (Figure 7a,b). Fire frequency increased with the number of aerial seeding treatments with the greatest increase at high precipitation from a frequency of 0.02 (once every 50 years) to 0.18 (once every 5.5 years, Figure 7a). Fire frequency increased slightly as the number of drill seeding treatments increased with the greatest increase in fire frequency at low spring precipitation (Figure 7b). Fire frequency did not change as the number of drill seeded treatments increased when aerial seeding treatments were one or fewer (Figure 7c). Fire frequency decreased with increasing drill seedings when a site had three aerial seeding treatments (Figure 7c).

Sites varied in fire frequency based on the seeding treatment used after the prior fire. Fire frequency increased with the number of aerial seeding treatments, but the frequency was less when sites were drilled alone or in combination with aerial seeding treatments (Figure 8a). Fire frequency was greatest on sites aerially seeded after the penultimate fire regardless of whether there was one or two drill seeding treatments previously (Figure 8b). Fire frequency was similar on drill seeded or unburned and unseeded sites prior to the most recent fire and was not affected by spring precipitation (Figure 8c). Fire frequency was greater for aerial and drill seeded sites than sites with only drill seeding and showed a slight decrease as precipitation increased (Figure 8c). Variation in fire frequency for unseeded and aerially seeded sites decreased as precipitation increased (Figure 8c). Fire frequency for aerially seeded sites decreased as precipitation increased (Figure 8c).

Multi-response permutation procedure analysis showed there was no significant difference in biomass when grouped by most recent seeding treatment method (T = −1.10, A = 0.016, p = 0.13). Biomass components were significantly different along the elevational gradient (p = 0.007), but within group homogeneity was low (T = −3.17, A = 0.058), suggesting a high degree of variability. When sites were nested by treatment type within elevation, there was significant differences among groups and within group similarity increased (T = −3.08, A = 0.109, p = 0.004), but there was still considerable variation within treatment and elevation groups. The best fit NMS ordination had two axes that explained much of the variation (R² = 0.9). The final stress was 11.7. Axis 1 was strongly correlated with shrub biomass and, to a lesser extent, perennial bunchgrass while axis 2 was strongly correlated with annual grass biomass and perennial bunchgrass (Table S4, Figure S2). The poor differentiation and within group variability with sites were grouped by elevation is demonstrated by the considerable overlap in ordination space (Figure 9).

The MRPP pairwise-comparisons, using the conservative α = 0.0115, showed biomass composition was significantly different between aerial seeded sites above 1521 m and four other groups: aerially seeded sites below 1020 m (T = −3.51, A = 0.211, p < 0.004); unseeded locations below 1020 m (T = −6.48, A = 0.247, p < 0.0002); aerially seeded sites between 1021 and 1270 m (T = −3.02, A = 0.096, p < 0.01) or 1271–1520 m (T = −3.92, A = 0.194, p < 0.003, see Table S5 for a complete list). Plant biomass composition at sites below 1020 m with no seeding treatment was the most variable, while aerial seeded sites between 1021 and 1270 m were the least variable, dominated by grasses (Figure 10a). Two plant functional groups were not normally distributed (annual grass W = 0.94, p = 0.005; shrub W = 0.93, p = 0.002) and the shrub biomass was heteroscedastic (W = 6.19, p < 0.001) requiring Kruskal–Wallis tests. Annual grass biomass was significantly different among seeding treatments and elevation groups (H₍₉₎ = 21.32, p = 0.01). Biomass among seeding treatments for bunchgrasses and shrubs tended to be different (bunchgrass F_(9,55) = 1.86, p= 0.08; shrub H₍₉₎ = 15.91, p = 0.07), but there was no significant difference in forb biomass among treatments (F_(9,55) = 1.29, p = 0.26). Bunchgrass and shrub biomass made up the largest biomass portion at sites above 1521 m, with annual grasses only found on aerially treated sites (Figure 10b). On aerial seeded sites below 1520 m, annual grass biomass comprised a more substantial proportion of biomass with little or no shrub biomass present (Figure 10b).

4. Discussion

Our research adds to a growing body of work that examines what shapes post-fire communities in the sagebrush steppe [54]. Artemisia tridentata subsp. wyomingensis communities exist in dry climates and require ≥ 80 years to recover to pre-fire levels of cover variability [23,24,80]. We found sites burned 2–7 times between 1955 and 2015 and the mean fire return intervals being as little as 3–5 years between 1995 and 2015 for aerial seeded sites. This makes the fire regime far outside the normal range of variability [23,24]. Pre-fire seeding treatment and site aridity consistently predicted fire regime characteristics. The number of seeding treatments predicted fire regime characteristics. The recovery rate in A. tridentata subsp. wyomingensis communities after a fire may increase with greater precipitation [81], but our results suggest a more nuanced effect. The effect of site moisture on fire regime characteristics was dependent on the type of seeding treatment used after a fire. Greater spring moisture had no effect on fire regime after drill seeding treatments; however, increased moisture—in the form of precipitation or lower vapor pressure deficit—resulted in fewer fires, lower fire frequency, and longer fire return intervals on aerial seeded sites. In our study area, precipitation increased, and vapor pressure deficit decreased, as elevation increased and resulted in cooler, more mesic summer conditions and longer, colder winters. Despite the elevation and climatic gradients, there was little differentiation in plant biomass when grouped by seeding treatment and elevation. When seeding treatments were nested in elevation, aerial treatments at high elevations showed differences with other treatments. Additionally, there was a difference in the amount of annual grass biomass along the elevation gradient.

The interaction between seeding treatment method and moisture availability suggests climatic influence over native bunchgrass and Bromus establishment likely contributed to the difference in fire regime characteristics. As moisture availability increased, fire number and fire frequency decreased at aerially seeded and unseeded sites but not at drill seeded sites. Established perennial bunchgrasses can inhibit Bromus plant growth [82,83], cover and density [84], but using aerial seeding is much less reliable at ensuring plant establishment than drill seeding [31]. Though the increase May rain since 1987 at low elevation sites may promote bunchgrass germination, rain between June and September has decreased since 1987, likely limiting the successful establishment of bunchgrasses [85,86,87] and opening habitat for Bromus establishment. In addition, mild winters with little snow at the low elevations and precipitation in May promote Bromus survival and fecundity [88,89]. Though native Vulpia species were present at the biomass collection sites, annual grass biomass was primarily composed of Bromus tectorum. Greater annual grass biomass at low elevations shown here demonstrates Bromus populations can take advantage of favorable growing conditions to increase their populations into the available space left open at lower elevation sites. Nonnative grasses can alter fire regime characteristics in a variety of native plant communities by increasing fire risk [26], creating continuous fuel beds [10], and/or increasing the likelihood of ignition [90] leading to increased fire frequency [12] and area burned [12]. That is particularly true of Bromus tectorum [20,21,26,91]. At high elevation sites, summer moisture appears adequate to promote bunchgrass establishment; when combined with snow accumulation and snow resident time, these climatic conditions likely decreased Bromus survival. The result is less annual biomass and fewer continuous fuels, and a shorter fire season that limit fires.

Fire regime characteristics were often determined by complex, nonlinear relationships between multiple variables. The effect of the year a site first burned, climate, and the total number of seeding treatments interacted to affect the number of fires. On sites that first burned prior to 1970, the number of fires did not change as the climate became drier, but the number of fires increased when the site first burned in or after 2000. This lack of change in the response of fire number before 1970 suggests a connection to changes in Bromus cover over time. Sites with that first burned prior to 1970 may have more time to recover from the fire. This is important because Bromus cover tends to higher in the years immediately following fire [27,92,93] and may decrease over time [94,95,96]. Bromus outcompetes bunchgrass seedlings [97,98,99], but some native species (e.g.-Poa secunda and Elymus multisetus) are able to establish or persist with Bromus [100,101]. If bunchgrasses establish, mature bunchgrasses can inhibit Bromus growth and reproduction [51,83,84], but the shift from a Bromus dominant community to a bunchgrass community could take decades [40,96]. In our study area, Bromus decreased as native bunchgrasses cover or density increased (unpublished data). Given the literature and our observations, it is likely bunchgrass communities on sites burned prior to 1970 had more mature bunchgrasses than recently burned sites that reduced Bromus cover consistent numbers of fires on sites burned before 1970 than more recently burn sites regardless of climate or the number of seeding treatments.

Fire return intervals responded differently from fire number to the interaction between time a site was first treated, climate, and the number of aerial treatments. The change in fire return interval had a markedly different relationship with time since first treated when the first seeding treatment occurred between 1970 to 1995 compared to the pre-1970 or post-1995 periods. That relationship suggests that changes in species composition used by the Bureau of Land Management in seed mixes for seeding treatments contributed to changes in fire return intervals. Agropyron desertorum, the dominant species used prior to 1970 [38,102], inhibits Bromus growth and seed production [82,84]. Bromus can recruit and survive in the presence of A. cristatum [99,103] and remain co-dominant in stands with A. cristatum and Poa secunda [104]. The increased use of A. cristatum after 1970 and Hycrest in the 1980s–1990s may have allowed Bromus to establish and even small amounts of Bromus cover can increase fire frequency [12,20,21]. The use of native grasses after 2000 has created a fire regime similar to the one created by treatments using A. desertorum.

The effect of the number of treatments on fire regime characteristics varied among treatment strategies. Sites with recent fires often received post-fire rehabilitation seeding treatments [38]; therefore, the number of treatments should increase as the number of fires increases. We found increasing the total number of seeding treatments on sites that first burned between 1965 and 1975 only increased the number of fires slightly, while fire numbers increased from one to three on sites that burned for the first time between 2005 and 2014. The number of fires and fire frequency increased as the number of treatments increased on sites where the pre-fire seeding treatment was aerial or unseeded but were unaffected when the most recent treatment included drill seeding. Aerial seeding grass species result in only low seeded species cover in forested systems [105]. In the sagebrush steppe, aerial seeding may increase nonnative perennial bunchgrasses at high elevations [40] but requires herbicide to reduce annual grass cover prior to seeding at low elevations [106]. One drawback of aerial seeding is the high potential of low seed-to-soil contact, resulting in low seeded species [31]. Great Basin native bunchgrass and forb recovery is not enhanced by aerial seeding after a fire [107]. Native species are inhibited when nonnative grasses, like A. cristatum, are included in the seed mix [40]. Aerial seeding does not increase Artemisia tridentata cover or density after a fire [33] even after two decades [40]. Drill seeding also had longer fire return intervals and lower fire frequency than other treatments in the last 20 years and mitigated the effect of aerial treatments. The number of fires on sites with drill and aerial treatments used in tandem was similar to those that were only drill seeded and had an intermediate fire return interval and fire frequency. Further, on sites that had aerial seedings in the past, treatments with both aerial and drill seedings before the most recent fire had in less frequent fires.

In our study area, the majority of recent aerial seed treatments included only A. tridentata seeded into recently burned, fire-resilient grasslands, such as Agropyron cristatum, A. fragile, or E. wawawaiensis [38]. The typical method of aerially seeding A. tridentata subsp. wyomingensis in the winter without drill seeding is unlikely to promote shrub establishment after fire [33]. In addition, the A. tridentata subsp. wyomingensis aerial seeding treatments in our study area were used on sites with established A. cristatum stands [38,64]. Agropyron cristatum outcompetes the seedlings of all A. tridentata subspecies [108], making it unlikely that A. tridentata subsp. wyomingensis will establish when aerially seeded into drill seeding treatments that used A. cristatum. Agropyron cristatum inhibits native forbs and grasses, leading to sites with low diversity or evenness [109,110,111] even when native species initially establish at the same rate as A. cristatum [112]. Therefore, sites seeded with both A. cristatum and native forbs and grasses in the 1980s and 1990s are unlikely to retain a diverse native species component. Since ‘Secar’ E. wawawaiensis is less palatable to cattle than other bunchgrasses used in recent post-fire rehabilitation seeding efforts [113,114] and avoided by goats [115], herbivore pressure on native forbs and other bunchgrasses likely increase on sites with E. wawawaiensis. This may shift the plant community to an unpalatable or invasive-dominated ecosystem [116], where fuels accumulate over time and create continuous fuel beds that promote more frequent fires. Whether the seeding included A. cristatum or E. wawawaiensis, the long-term result is a low diversity grassland that may be more susceptible to invasion by small-seeded, winter annual species [117], like Bromus.

Many factors may affect fire regimes, including average climate and interannual variation in weather [54,118], available fuels [119,120,121], historical fire management practices [52,53], and interactions between climate, fuels, and management [14]. There is a lack of understanding how seeding treatments affect fire regimes across arid landscapes [52,53,122]. Though some research has looked at post-fire seeding treatments after one fire [33,40,84,94,123], this study included sites with multiple wildfires and treatment histories over >50 years. Our results suggest that spring moisture is an important determinant of the fire regime, but the effect of spring moisture varies among seeding techniques. Spring moisture likely increased B. tectorum fecundity. Greater precipitation in May, could have stimulated germination of perennial bunchgrass, but the recent shifts in decreased precipitation between June and September would have increased water stress and mortality for seeded species. The consequence would be reduced bunchgrass establishment and less competition for B. tectorum in subsequent years, leading to increased B. tectorum establishment. As precipitation shifts and temperatures increase in the arid western United States, it may be necessary to provide additional summer water in the first year after seeding treatments to increase post-fire rehabilitation success.

Drill seeded sites had fewer fires over the 60-year study period and in recent years had longer fire return intervals and less frequent burns than aerial seeded sites. This trend did not differ along the moisture gradient and fire regimes on recently drilled sites had similar fire regimes to sites with older treatments. Aerially seeded sites had more frequent fires and were dominated by A. cristatum, E. wawawaiensis, or Bromus. The accumulation of fine fuels from Bromus combined with unconsumed, senesced A. cristatum and E. wawawaiensis likely contributed to frequent fires. Drill seeding helped mitigate the frequent and greater number of fires associated with aerial seeding treatments. If the goal is to minimize the number and frequency of fires [124], drill seeding larger areas with native bunchgrasses and forbs may help achieve that goal on drier sites historically dominated by A. tridentata subsp. wyomingensis. Additional treatments, such as the reduction of A. cristatum and E. wawawaiensis, may better establish diverse communities that are more palatable to livestock and wildlife, more resilient to Bromus invasion, and less prone to fire.