Post-Fire Vegetation Response in a Repeatedly Burned Low-Elevation Sagebrush Steppe Protected Area Provides Insights About Resilience and Invasion Resistance

Sagebrush steppe ecosystems are threatened by human land-use legacies, biological invasions, and altered fire and climate dynamics. Steppe protected areas are therefore of heightened conservation importance but are few and vulnerable to the same impacts broadly affecting sagebrush steppe. To address this problem, sagebrush steppe conservation science is increasingly emphasizing a focus on resilience to fire and resistance to non-native annual grass invasion as a decision framework. It is well-established that the positive feedback loop between fire and annual grass invasion is the driving process of most contemporary steppe degradation. We use a newly developed ordinal zero-augmented beta regression model fit to large-sample vegetation monitoring data from John Day Fossil Beds National Monument, USA, spanning 7 years to evaluate fire responses of two native perennial foundation bunchgrasses and two non-native invasive annual grasses in a repeatedly burned, historically grazed, and inherently low-resilient protected area. We structured our model hierarchically to support inferences about variation among ecological site types and over time after also accounting for growing-season water deficit, fine-scale topographic variation, and burn severity. We use a state-and-transition conceptual diagram and abundances of plants listed in ecological site reference conditions to formalize our hypothesis of fire-accelerated transition to ecologically novel annual grassland. Notably, big sagebrush (Artemisia tridentata) and other woody species were entirely removed by fire. The two perennial grasses, bluebunch wheatgrass (Pseudoroegneria spicata) and Thurber's needlegrass (Achnatherum thurberianum) exhibited fire resiliency, with no apparent trend after fire. The two annual grasses, cheatgrass (Bromus tectorum) and medusahead (Taeniatherum caput-medusae), increased in response to burn severity, most notably medusahead. Surprisingly, we found no variation in grass cover among ecological sites, suggesting fire-driven homogenization as shrubs were removed and annual grasses became dominant. We found contrasting responses among all four grass species along gradients of topography and water deficit, informative to protected-area conservation strategies. The fine-grained influence of topography was particularly important to variation in cover among species and provides a foothold for conservation in low-resilient, aridic steppe. Broadly, our study demonstrates how to operationalize resilience and resistance concepts for protected areas by integrating empirical data with conceptual and statistical models.

Sagebrush steppe ecosystems are threatened by human land-use legacies, biological invasions, and altered fire and climate dynamics. Steppe protected areas are therefore of heightened conservation importance but are few and vulnerable to the same impacts broadly affecting sagebrush steppe. To address this problem, sagebrush steppe conservation science is increasingly emphasizing a focus on resilience to fire and resistance to non-native annual grass invasion as a decision framework. It is well-established that the positive feedback loop between fire and annual grass invasion is the driving process of most contemporary steppe degradation. We use a newly developed ordinal zero-augmented beta regression model fit to large-sample vegetation monitoring data from John Day Fossil Beds National Monument, USA, spanning 7 years to evaluate fire responses of two native perennial foundation bunchgrasses and two non-native invasive annual grasses in a repeatedly burned, historically grazed, and inherently low-resilient protected area. We structured our model hierarchically to support inferences about variation among ecological site types and over time after also accounting for growing-season water deficit, fine-scale topographic variation, and burn severity. We use a state-and-transition conceptual diagram and abundances of plants listed in ecological site reference conditions to formalize our hypothesis of fire-accelerated transition to ecologically novel annual grassland. Notably, big sagebrush (Artemisia tridentata) and other woody species were entirely removed by fire. The two perennial grasses, bluebunch wheatgrass (Pseudoroegneria spicata) and Thurber's needlegrass (Achnatherum thurberianum) exhibited fire resiliency, with no apparent trend after fire. The two annual grasses, cheatgrass (Bromus tectorum) and medusahead (Taeniatherum caput-medusae), increased in response to burn severity, most notably medusahead. Surprisingly, we found no variation in grass cover among ecological sites, suggesting fire-driven homogenization as shrubs were removed and annual grasses became dominant. We found contrasting responses among all four grass species along gradients of topography and water deficit, informative to protected-area conservation strategies. The fine-grained influence of topography was particularly important to variation in cover among species and provides a foothold for conservation in low-resilient, aridic steppe. Broadly, our study demonstrates how to operationalize resilience and resistance concepts for protected areas by integrating empirical data with conceptual and statistical models.

INTRODUCTION
Sagebrush steppe ecosystems across western North America have been lost or degraded by cumulative impacts of human landuse legacies, biological invasions, and altered fire and climate dynamics (Davies et al., 2011). Steppe protected areas are therefore of heightened conservation importance but are few in number (Caicco et al., 1995;Noss et al., 1995;Stoms et al., 1998) and vulnerable to many of the same cumulative impacts broadly affecting sagebrush steppe. To address this problem, sagebrush steppe conservation science is increasingly emphasizing a focus on resilience to fire and resistance to non-native annual grass invasion as a decision framework (e.g., Chambers et al., 2019), now that it is well-established that the positive feedback loop between fire and annual grass invasion is the driving process of most contemporary steppe degradation (D'Antonio and Vitousek, 1992;Brooks et al., 2004;Chambers and Wisdom, 2009;Balch et al., 2013). Resilience to fire, which is defined as site capacity to return to pre-fire structure and function (e.g., native perennial grass cover is recovered over time), and the annual grass invasion resistance of sagebrush steppe ecosystems varies predictably along topographic-moisture gradients (Chambers et al., 2014a). The position of steppe landscapes along these gradients can be used to guide conservation decision-making (e.g., Chambers et al., 2014b). For example, high elevation (e.g., 1,500 m) steppe with mountain big sagebrush (Artemisia tridentata subsp. vaseyana) overstory tends to be much less invasible by annual grasses of Eurasian origin than low elevation (e.g., <1,000 m) steppe with Wyoming and basin big sagebrush (A. t. subsp. wyomingensis and A. t. subsp. tridentata) overstories (Chambers et al., 2014a). Low-elevation, low-productivity sites exhibit higher inertia and are much more difficult and expensive to restore (Suding et al., 2004;James et al., 2013), and when considered from a regional strategic perspective, high-elevation sites may therefore be considered as being of higher conservation priority (Chambers and Wisdom, 2009;Chambers et al., 2019). However, this confronts those low-elevation protected-areas that exhibit inherently low resilience to fire and resistance to invasion with an existential challenge requiring pro-active strategic landscape prioritization and triage. To fail to do so will likely result in the loss of low-elevation steppe protected areas as places of conservation value given current trajectories of degradation.
Like all protected areas, those in low-elevation sagebrush steppe are fixed in geographic space even as the ecological gradient space within which they occur (e.g., temperature and precipitation, fire regime) shifts over time toward novel conditions with less native biodiversity and lower conservation value. For example, National Park Service sites containing sagebrush steppe across the Columbia Basin and Snake River Plain have all experienced substantial fire-induced shifts toward novel annual grasslands in recent decades (Bangert and Huntly, 2010;Ellsworth and Kauffman, 2010;Powell et al., 2013;Rodhouse et al., 2014;Reed-Dustin et al., 2016;Esposito et al., 2019;Ellsworth et al., 2020;Nicolli et al., 2020), yet each of these parks are also charged with preserving native steppe as established by the 1916 National Park Service Organic Act. In order to meet this charge, these parks and others faced with the same kinds of challenges to protected-area relevancy need to develop granular, place-based ecological insights that are useful for park-scale preservation and restoration planning and contextualized within the resilience and resistance framework.
Inferences and predictions from statistical models fueled by field and remotely sensed observations can be an effective way to support park-specific evidence-based conservation decisionmaking, particularly when integrated into adaptive management supported by recurring monitoring. For example, Rodhouse et al. (2014) used field-based monitoring data and statistical regression to predict where native foundation bunchgrasses were most likely to occur along topographic gradients within the John Day Fossil Beds National Monument (JODA) in north-central Oregon, USA. These predictions were inferred to indicate apparent resilience and resistance and used to prioritize implementation of park burned-area emergency response and weed management plans (Hoh et al., 2015). More recently, the soil and ecological site attributes provided by the US Natural Resource Conservation Service (NRCS) web soil survey have been used to map current and future potential resilience and resistance across the sagebrush steppe biome (Maestas et al., 2016;Bradford et al., 2019), including in protected areas. For example, each of the separate management units of JODA are mapped almost entirely as low resilience currently and in the future (Figure 1A), a predicament also faced by other parks in the region. The mappable linkages between soil surveys and potential resilience and invasion resistance creates important opportunities to operationalize these concepts into a decision-making framework. However, it is evident that the generalized assessments by Bradford et al. (2019, e.g., Figure 1A) and others will need to be empirically evaluated and translated with more granularity and precision for park operations.
To support this need, we evaluate a 7-year (2011-2017) study of native and non-native grass species' responses to a 2011 wildfire in a repeatedly burned and inherently low resilience protected-area landscape in the John Day Fossil Beds National Monument. We use a newly developed class of zero-augmented ordinal regression (Irvine and Rodhouse, 2010;Irvine et al., 2016Irvine et al., , 2019 to model variation in quadrat-based observations of above-ground grass cover categories over time since fire and in relation to ecological site, topographic position, fire severity, and growing-season (October-May preceding each field survey) water deficit. The zero-augmented (hurdle-at-zero) model approach allowed us to explicitly evaluate patterns in occurrence (presence/absence) and in abundance (cover categories >0) simultaneously. We structured our model hierarchically so that observations could be grouped within ecological sites, an important step in contextualizing model-estimated trends over space and time within the resilience and resistance framework. We used the same ecological site descriptions developed through the NRCS web soil survey (NPS, 2012) as was used by Bradford et al. (2019) to develop their regional resilience and resistance map, although we retained the NRCS soil survey map unit polygons (see Figure 1B) rather than the coarse grid employed by Bradford et al. (2019). Importantly, we were able to utilize attributes of the ecological sites, including associated state-andtransition diagrams (sensu Stringham et al., 2003;Chambers et al., 2014c) and reference plant species lists as articulated hypotheses about fire-driven state transitions (Figure 2). Specifically, we hypothesized that the entire study area would exhibit accelerated transition to the annual grass state (state 3, Figure 2) because of the 2011 fire and earlier observations about shrub cover loss from previous burns, but with variation in that transition among ecological sites as a result of topographically-mediated effective soil moisture and site productivity that confers relative resilience and invasion resistance.
We used the models to examine cover patterns of two large-stature native foundational bunchgrass species and two invasive annual grasses of Eurasian origin. The two perennial species, bluebunch wheatgrass (Pseudoroegneria spicata; wheatgrass, hereafter) and Thurber's needlegrass (Achnatherum thurberianum; needlegrass, hereafter), comprise >75% of the total plant community biomass of the hypothetical historic reference state (state 1) in Figure 2 (NRCS citation), and ∼10% of the biomass in the invaded state (state 2 in Figure 2). The two annual grasses that we examined, cheatgrass (Bromus tectorum) and medusahead (Taeniatherum caput-medusae), are the dominant invasive species associated with both the invaded state and the annual grass state (states 2 and 3 in Figure 2). The relative abundances of these four species are therefore useful indicators of these hypothesized states described in Figure 2 and of resilience to fire and resistance to annual grass invasion for this study area, more generally. The inferences supported by our statistical survey design and our models provide insights about resilience and resistance that will support a more strategic approach to park steppe conservation than what has been pursued previously, specifically revealing the importance of topography in our study area as a strong mediator of resilience and resistance which can provide a conservation foothold. More generally, our approach demonstrates the importance of using empirical data and both conceptual and statistical models in FIGURE 2 | A simplified state-transition diagram for the aridic low-elevation sagebrush steppe ecological sites within the Clarno unit of the John Day Fossil Beds National Monument, modified from Natural Resources Conservation Service ecological site descriptions (NPS, 2012), Chambers et al. (2014c), andEllsworth et al. (2020). This diagram highlights three fundamental states with progressively steep threshold transitions into ecologically-novel conditions dominated by invasive annual grasses, from state 1 to state 2 and from state 2 to state 3. Note that other phases within these states dominated by big sagebrush and western juniper are possible (NPS, 2012;Chambers et al., 2014c;Ellsworth et al., 2020) but were not observed during our study due to past wildfires.
an integrated way to operationalize resilience and resistance concepts for place-based conservation decision-making in an increasingly at-risk North American biome.

Study Area
We evaluated post-fire response of native and non-native grasses over a 7-year period from 2011 to 2017 in the Clarno Unit of the John Day Fossil Beds National Monument, in north-central Oregon, USA (Figure 1). The Clarno Unit encompasses 5,800 ha of upland sagebrush steppe. Elevations range from 421 to 725 m. Clarno is topographically rugged, comprised of highly weathered volcanic ash-derived clay to stony clay-loam soils.
Thirty-year mean annual precipitation is ∼27 cm (ecological sites are attributed to 22-30 cm precipitation zones), and precipitation falls mostly as rain and ephemeral snow during October to May, which drives the growing season for the four cool-season grasses examined. Wildfires occurred in the study area in 1985in , 1994in , 1995 Figure 1.1), and we assumed that all quadrat locations experienced at least 2 burns within this time period, although lack of a published perimeter for the 1985 fire and variation in burn severity make it difficult to establish exact fire history. This fire history represents an ∼10-20-year fire frequency, which is substantially shorter than the multidecadal to centuries intervals established for low-elevation aridic Wyoming and basin big sagebrush ecosystems (e.g., Chambers and Wisdom, 2009;Bukowski and Baker, 2013).
Historically, Clarno was grazed by sheep and cattle for many decades prior to establishment as a National Park Service unit and effective fencing was constructed. Without quantitative information on grazing history (e.g., variation in timing and duration within the study area), however, it was impossible to account for land-use legacy in the study. Nonetheless, this legacy is an important aspect to the unfolding story of environmental change in the Monument and contextualizes our results.
Sagebrush steppe plant communities in Clarno were mapped prior to the 2011 wildfire by Erixson et al. (2011) following the National Vegetation Classification System to include ∼20% Wyoming-type sagebrush steppe (dominated by A. t. subsp. wyomingensis) and >50% as herbaceous, dominated by large-stature native perennial bunchgrasses (wheatgrass and needlegrass) or invasive annual grasses (cheatgrass and medusahead; Supplementary Figure 1.2). Western juniper (Juniperus occidentalis) trees were present in open savannah along several drainage bottoms. Five ecological sites were mapped in the study area by NRCS (NPS, 2012; Figure 1B), each of which share important abiotic attributes that pre-dispose the entire study area to inherently low resilience to fire and resistance to invasion, most importantly shared mesic soil temperature and aridic soil moisture regimes (Chambers et al., 2014a;Maestas et al., 2016;Bradford et al., 2019). The state 1 (Figure 2) reference conditions for each of the ecological sites list basin big sagebrush (A. t. subsp. tridentata) as the dominant shrub, although Wyoming big sagebrush was more common on these sites, as mapped by Erixson et al. (2011). Differences among ecological sites include two (ecological sites R010XB022OR and R010XB064OR) that were characterized as north-facing slopes dominated by wheatgrass more than needlegrass, and three (R010XB044OR, R010XB052OR, R010XB053OR) predominantly south-facing slopes dominated primarily by needlegrass. Notably, topographic heterogeneity measured from a 10-m digital elevation model was much greater than indicated by ecological site descriptions, and each ecological site included relatively wide gradients of slopes and aspects ( Supplementary Figures 1.3-1

Field Sampling Methods
We used vegetation monitoring data collected in the Clarno Unit of the John Day Fossil Beds National Monument during 2011-2017 as part of the National Park Service's vital signs monitoring program (Fancy et al., 2009). Field methods followed a protocol (Yeo et al., 2009) for ocular estimation of aboveground foliar cover categories of plant species found within 1-m 2 quadrats, following Daubenmire (1959) 7-level system (including a 0 cover class; 0% = 0, >0-5% = 1, >5-25% = 2, >25-50% = 3, >50-75% = 4, >75-95% = 5, and >95% = 6). Quadrats were randomly located in the field following a spatiallybalanced equal probability randomization (the Generalized Random Tessellation Stratified [GRTS] algorithm, Stevens and Olsen, 2004). We used the R library spsurvey (Kincaid and Olsen, 2019) to implement GRTS sampling. The GRTS design ensures spatial balance and representativeness of the sample within the area of inference (i.e., the sample frame) and has advantages over alternative approaches such as simple random and stratified random sampling that include more efficient design-unbiased population variance estimation, spatial balance, and reduction in autocorrelation, and flexibility for replacing sample units (Stevens and Olsen, 2004). Randomization was constrained within a sampling frame constructed in a GIS from vegetation maps that included all upland vegetated areas <30 degrees (58%) slope and that also excluded barren fossil-bearing ash beds that are sensitive to trampling (Supplementary Figure 1.2). Steep slopes >30 degrees presented safety risks to field observers and represented a relatively small amount of land removed from the sampling frame ( Supplementary Figure 1.2). New independent random samples were drawn each year. Quadrats were not permanently marked and never revisited (a "neverrevisit design" sensu Urquhart and Kincaid, 1999;MacDonald, 2003). The Monument discouraged use of permanent markers (e.g., steel rebar) but large independent samples each year offset the potential loss of power from not using permanent plots. Sample sizes ranged from 85 to 211 annually, totaling 1,338 quadrat observations (Table 1). Sample sizes were informed by a statistical power analysis developed specifically for the protocol, and examined via simulations (Irvine and Rodhouse, 2010;Irvine et al., 2016). Our annual samples covered the complete ranges, means, and standard deviations of the environmental gradients (slope, aspect, elevation) within the sample frame. In all years but 2012, our achieved sample sizes were twice the size recommended (∼n = 100) by Irvine and Rodhouse (2010). Surveys were conducted annually in early June, and species cover was visually assessed and binned into Daubenmire categories. Data are available on-line at: https://irma.nps.gov/DataStore/ Reference/Profile/2278547.

Statistical Analyses
We summarized plant species cover by plotting stacked bar charts with the proportion of quadrats in each cover class. The empirical proportion,p = x n , is a design-unbiased population estimator for categorical data (Agresti, 2010). We used a multi-level (hierarchical) ordinal zero-augmented beta regression (OZAB; Irvine et al., 2016) to examine the patterns of cover, in Daubenmire categories, along gradients of topography, fire severity, growing season water deficit, time (year), and ecological site. The OZAB model specifies a more flexible latent beta distribution to model plant cover (Damgaard and Irvine, 2019) rather than a latent logistic distribution as is done with the more familiar proportional-odds logistic regression (POLR; Agresti, 2010). An important distinction between OZAB and POLR is the ability to use the more flexible 2-parameter beta distributions to remove the restrictive assumption that the odds of moving into a higher plant cover category is the same between all categories. Adding a hurdle-at-zero model adds additional flexibility and the OZAB model provides an integrated framework to explore both presence/absence and abundance patterns and processes . Furthermore, the OZAB model directly links the cover classes to the partiallyobserved latent percent cover, allowing model parameters (e.g., environmental gradient coefficients) to be interpreted in terms of changes in mean percent cover rather than in terms of cumulative odds ratios, which are not intuitive . The OZAB model provides more precise parameter estimates than POLR or when using linear regression with category midpoints (Irvine and Rodhouse, 2010;Irvine et al., 2016). We reported the unconditional mean cover (and 95% credible intervals), as a proportion, among years and ecological sites, for each species.
Following, Irvine et al. (2016Irvine et al. ( , 2019, the observed cover categories were treated as a grouped continuous ordinal response from a multinomial distribution with J + 1 categories, [C i |π i ] = Multinomial(1, π 0i , . . . , π Ji ), for i = 1, . . . n observations (quadrats or plots) with category levels of 0, 1, . . . , J and Pr C i = j = π ji . The key advantage of the OZAB approach is specifying the latent variable (Y * i ) distribution as a mixture of a continuous beta distribution and a point mass at 0, as follows: where Y * i α, µ, φ is the zero augmented beta distribution. The abundance portion of the model was parameterized for regression by Ferrari and Cribari-Neto (2004) φ+1 . For the plant cover class data used in this study, µ is interpreted as the mean proportion cover conditional on presence. The unconditional mean cover is estimated by αµ.
The observed cover class category was defined as a discretized version of Y * i for observation i as follows: where θ represents the breaks between each cover class. The Daubenmire (1959) scale defines the vector as {θ 0 , . . . , θ 6 } = {0, 0.5, 0.25, 0.5, 0.75, 0.95, 1.0}. The multinomial probabilities are: The last two digits (in parentheses) of the ecological site alphanumeric code used by the Natural Resources Conservation Service is used in Figures 4, 5 for brevity.
Where F Beta denotes the cumulative distribution function for the partially observed conditional percent cover Y, which is beta distributed with parameters µ and φ.
The mean of proportion cover (µ) and probability of occurrence (α) were allowed to be influenced by explanatory variables (environmental gradient covariates) by using logit (µ i ) = X i β and logit (α i ) = W i γ , where X and W represent a matrix of covariates obtained for each plot . Continuous variable covariates were input as mean-centered and standardized quantities and categorical variable covariates were input as factors anchored at level 0 (e.g., unburned and study year 1) which facilitated computation and interpretation of µ and α. We specified diffuse normal priors on regression coefficients for γ ([γ ] = Normal (0, 0.01)). Diffuse normal priors on abundance regression parameters β [β] = Normal (0, 0.01) and a diffuse Gamma prior for the precision parameter ([φ] = Gamma (0.1, 0.1)) were also specified.

Environmental Gradients for Model Covariates
We used a 10-m digital elevation model (US Geological Survey National Elevation Dataset, http://ned.usgs.gov/) to obtain a measure of topography for each plot by calculating sin(slope) × cos(aspect), producing a variable that ranged from −1 to 1 (Supplementary Figures 1.3-1.5). Steep south-facing slopes trend toward −1 and steep north-facing slopes trend toward 1. This and similar calculations have been widely used to approximate insolation (Stage, 1976;McCune and Keon, 2002) which can indicate topographically mediated variation in effective soil moisture. We estimated burn severity for the 2011 fire by comparing two Landsat images immediately before (August 17, 2011) and after (September 9, 2011) the fire and estimating delta normalized burn ratio (dNBR; Supplementary Figure 1.6; Keeley, 2009). We used the US Geological Survey (USGS) Earth Explorer (earthexplorer.usgs.gov) to order and acquire Landsat 5 (based on USGS availability) normalized burn ratio scenes provided as part of the USGS suite of Landsat Surface Reflectance-Derived Spectral Indices products (see Masek et al., 2006). The two scenes subtracted to produce a dNBR raster were radiometrically corrected by USGS and ranked as level 1 precision and terrain corrected product (L1TP) quality. We binned severity values into five relative rankings categories, from unburned to high severity (Supplementary Figure 1.6) following approximately the thresholds described by Keeley (2009) and Key and Benson (2006), but tuned to our study area by an examination of the empirical distribution of observed dNBR values in our study and finding breakpoints that would ensure sufficient sample size within each category for inclusion as a 5-level factor model covariate. Including burn severity as a factor improved model computation and interpretation in terms of severity level description.
Water deficit during the growing season for cool-season grasses, which in our study region occurs during October-May preceding each sampling season (e.g., see ecological site descriptions, NPS, 2012) was used as an additional model covariate (Supplementary Figure 1.7). Water deficit is calculated as potential evapotranspiration minus actual evapotranspiration, and is a parsimonious, integrative variable for use in statistical models with fundamental influence on plant growth and vegetation, particularly in arid regions (Stephenson, 1990;Dilts et al., 2015;Thoma et al., 2020). Water deficit is derived from Thornthwaite-type water balance equations (Thornthwaite and Mather, 1955) that describe water availability to plants, accounting for latitude and solar loading, temperature and precipitation, and soil moisture (holding capacity). Water deficit is a more mechanistic description of water stress on plants in arid lands than using temperature and precipitation alone, or drought indices such as normalized difference water index that only measure leaf water content (Stephenson, 1990(Stephenson, , 1998Lutz et al., 2010;Dilts et al., 2015;Thoma et al., 2016Thoma et al., , 2020. For our study, water deficit was estimated by Thoma et al. (2020) for the National Park Service from Thornthwaite water balance equations calculated at daily time steps across the continental US at 1 km resolution, following methods developed by Thornthwaite and Mather (1955) but with modifications described by Lutz et al. (2010) and Thoma et al. (2020). The 1 km resolution was constrained by the required daily inputs of precipitation and temperature from the Daily Surface Weather and Climatological Summaries (Thornton et al., 2018;Daymet; https://daymet.ornl.gov/). October-May water deficit during our study was lowest (∼18 mm) in 2014 and highest (∼23 mm) in 2016 (Supplementary Figure 1.7).

Model Fitting
We used a Bayesian approach to fit OZAB models, using three chains and 5,000 iterations of Markov chain Monte Carlo (MCMC) algorithms within the JAGS software (Plummer, 2003), launched from R (R Core Team, 2019) using the rjags package (Plummer, 2019). Slight variations on the maximum likelihood estimations for β and γ were used to initialize the three chains, as well as the corresponding midpoint values on the proportion cover scale for the unobserved beta latent variable.
We developed a hierarchical (i.e., mixed-effects or partial pooling sensu Gelman and Hill, 2007) model structure, grouping quadrat cover observations by ecological site type and year, allowing for a unique coefficient to be estimated for each ecological site and year in order to reflect the survey design and support inferences over space (ecological site) and time. Covariate matrices for both occurrence and abundance included topography, burn severity, and water deficit. We added a quadratic parameter for topography to allow curvature in the modeled influence of topography on cover. Topography and burn severity varied over space at fine scales (10 m and 30 m, respectively) but not time (Supplementary Figures 1.5, 1.6). Water deficit varied coarsely (1 km) over space because of data resolution and because weather and climate are synoptic processes, and varied over time as a year effect, although deficit was not correlated with year ( Supplementary Figure 1.7). Full models were fit for both the abundance and presence portions of the OZAB model for each species separately. Unconditional means and 95% credible intervals were estimated from the posterior distributions of the three MCMC chains.
In order to evaluate model performance, we fit models with 90% of the dataset (1,251 plots), holding 10% of the dataset for model checking. We conducted posterior predictive checks for each model by assessing discrepancy between observed data and predicted data and by calculating Kendall's τ to measure the degree of concordance between observed cover classes and the predicted cover classes.

RESULTS
Woody overstory species were already rare in the study area prior to 2011 due to past wildfires and were almost entirely removed by the 2011 fire. This can be seen in the pairs of photographs in Supplementary Figure 1.8 that were taken at 2 locations within the study area in 1988 and 2012. In 2011, we estimated that frequency of occurrence (empirical proportion when cover >0%) of big sagebrush was only 6% (12 out of 215 quadrats; Supplementary Figure 1.9). This declined to 0% in 2012 after the 2011 fire and remained so through the duration of the study except that big sagebrush was encountered in 1 quadrat in 2014 (Supplementary Figure 1.9). Other shrubs were also removed entirely by fire, including gray rabbitbrush (2% in 2011 and 0% thereafter). The fire-resprouting sub-shrub broom snakeweed (Gutierrezia sarothrae) was the only woody species to persist on the landscape after fire, declining in 2012 but remaining rare throughout the study period (Supplementary Figure 1.9). Forbs were rare throughout the study and did not show response to fire. For example, locoweed (Astragalus spp.) and biscuitroot (Lomatium spp.) were the most abundant perennial forbs but did not vary in frequency of occurrence or in abundance over the seven years of study (Supplementary Figure 1.9).
The two perennial grasses wheatgrass and needlegrass and the two annual grasses cheatgrass and medusahead that were the focus of the study were the most frequently encountered and most abundant species throughout the study area over all 7 years of study (Figure 3). Cheatgrass was particularly ubiquitous, occurring in over 90% of quadrats during all years and in almost 100% of quadrats by the end of the study. Frequency of medusahead occurrence increased over the study from 34% in 2011 to 64% in 2017, with a concomitant increase in higher cover classes (Figure 3). All four of these species also exhibited resilience to the 2011 wildfire, maintaining or exceeding pre-burn cover by the end of the study (Figure 3). Notably, the two annual grasses appeared to increase over the study, in contrast with the two perennial bunchgrasses that exhibited no apparent trend. However, there was no apparent inflection in the medusahead increase following fire in 2012 as was evident for cheatgrass (Figure 3).
Model results for each of the grass species also revealed the apparent resilience to the 2011 fire among all four species, with strikingly low variation among ecological sites. Unconditional mean percent cover estimates, interpreted as average cover across each ecological site and year on unburned areas (burn severity = 0), for flat ground (topography and topography 2 = 0), and average water deficit conditions, showed strong effects of the fire on cheatgrass but no discernible effects on the other three species (Figure 4). Cheatgrass cover dropped precipitously the year after the burn but returned to pre-burn cover by 2015, 4 years after fire.
Burn severity had a strong positive influence on the model estimates of cover for the two annual grasses and a modest negative influence on cover for the perennial bunchgrasses. This influence was clear when comparing unconditional mean estimates for high burn severity (burn severity = 4; Figure 5 vs. Figure 4). This increase was the result of burn severity coefficient estimates that yielded ∼300-350% increases in both occurrence and abundance of the two annual grasses for unitlevel increases in burn severity (interpreting posterior means), with the strongest positive influence on medusahead. Coefficient estimates for the bunchgrasses were only slightly negative but included 0 within relatively narrow 95% credible intervals (Figure 6).
Topographic heterogeneity had very strong, non-linear, and contrasting influences on both the occurrence (hurdle-at-zero) and abundance (cover classes >0) model components estimated for all four grass species (Figure 6). When visualized as fitted lines along the topographic gradient, these influences show needlegrass occurrence strongly and needlegrass abundance modestly (due to rarity of high cover class observations for this species) associated with south facing slopes; cheatgrass abundance also associated with south-facing slopes but cheatgrass occurrence ∼100% along the entire gradient (although with increasing northward uncertainty shown in credible interval width because the few quadrats without cheatgrass occur only on north-facing slopes); medusahead associated with flat ground and shallow slopes of both aspects; and both the occurrence and abundance of wheatgrass strongly associated with north-facing slopes (Figure 7).
The two annual grasses were most strongly influenced by water deficit, but notably in opposite directions (Figure 6). Cheatgrass occurrence and abundance was negatively associated with dry high deficit years, whereas medusahead increased during high deficit years. Thurber's needlegrass was also modestly associated with water deficit, reflecting a higher drought tolerance than bluebunch wheatgrass, as was also evident from the examination of topographic patterns (south-slope abundance and occurrence curves; Figure 7).
Predictive performance was similar among OZAB models for each species, with 64.4-88.7% of predictions within 1 cover class (needlegrass: 88.5%, medusahead: 76.5%, cheatgrass: 64.4%, and wheatgrass: 73.7%, Supplementary Figures 1.10, 1.11, Supplementary Table 1.1). The discrepancies between the observed cover class data (10% hold-out data set) and the predicted cover class data from all MCMC iterations (observed classpredicted class) in each model were most often 0 (Supplementary Figure 1.10). Patterns of predicted cover classes were similar to observed cover classes, but the OZAB models predicted less range in cover classes for all species (Supplementary Figures 1.10, 1.11, Supplementary Table 1.1). The first and third quantiles of the posterior distribution of Kendall's τ based on predictions from each species' model were 25.6-34.0% for needlegrass, 19.0-26.0% for medusahead, 12.2-19.8% for cheatgrass, and 27.1-32.4% for wheatgrass (Supplementary Figures 1.10, 1.11, Supplementary Table 1.1).

DISCUSSION
We used empirical observations and statistical models of native perennial and non-native invasive annual grass cover to evaluate the hypothesis of unidirectional fire driven state-transition to ecologically novel annual grassland in a historically grazed, repeatedly burned, and inherently low resilience/low resistance sagebrush steppe protected area. Due to historical legacies of land use and previous fires and ecological site state transition diagrams, we anticipated that the entire park unit would accelerate to an annual grass dominated state following the studied fire that occurred in 2011, but with important spatial variation among ecological sites. We found substantial evidence in favor of this hypothesis that validated the generalized state and threshold-transition dynamics understood for this aridic system but with less spatial variation among ecological sites than anticipated and with more variation along the fine-grained environmental gradients examined. Variation among gradients underscored the importance of obtaining local-scale place-based insights to implement more effective strategic management and restoration actions. The topographic position of quadrats exerted the strongest influence on variation on grass cover and apparent resilience to fire and resistance to invasion FIGURE 4 | Unconditional means estimated for each ecological site from ordinal zero-augmented beta regression models fit to ordinal cover class observations made in 1,251 (90% of the dataset for training models) quadrats surveyed during 7 years of monitoring in Clarno, John Day Fossil Beds National Monument. Estimates were plotted here for model severity factor levels set to 0 burn severity (unburned), for flat ground (topography and topography 2 = 0), and for average water deficit conditions. Ecological sites are labeled for brevity with the last two numbers of the alphanumeric code used by the Natural Resources Conservation Service (see Table 1 for reference). among the two native foundational bunchgrass species examined (see photographs in Supplementary Figure 1.8). This insight provides a key conservation opportunity within an otherwise heavily invaded park. Steep north-facing slopes confer resilience and invasion resistance because of lower insolation and higher effective soil moisture and the Monument conservation strategy can be customized along topographic gradients (Chambers et al., 2014aRodhouse et al., 2014).
Variation in grass cover was minimal among the five ecological sites over time in our study. This was surprising and likely arises for three reasons. First, NRCS soil surveys are inherently imprecise for granular studies such as ours due to the logistical and financial resources available to identify and map soils and associated attributes at large scales (NRCS, 1995;Maestas et al., 2016). Because of this, there is a much higher amount of environmental heterogeneity included within mapped soil polygons than is indicated by ecological site descriptions. For example, variation in topography (slope and aspect) was as high within ecological sites as it was among ecological sites (Supplementary Figures 1.3-1.5). Steep, primarily north-facing slopes found in each of the ecological sites retained native perennial grasses more than on shallow slopes ( Supplementary Figure 1.8). A second reason for low variation in grass cover among ecological sites over time in our study area was the shared attributes of soil temperature and moisture regime across all sites and only subtle differences in soil texture, historical and degraded plant community composition, and other site attributes described by NRCS for the ecological sites in the study area. Similarly, a third reason for low variation among ecological sites was the homogenizing effect (sensu Foster et al., 2013) of the 2011 fire, with likely cumulative influences from past fires and grazing (Chambers and Wisdom, 2009;Davies et al., 2012;Bernards and Morris, 2016;Mahood and Balch, 2019;Wood et al., 2019;Ellsworth et al., 2020). This homogenization has occurred throughout the study area across ecological sites, as community composition of all ecological sites converged toward an annual grass state. The practical implications of this outcome are that ecological site descriptions provide important general information to guide expectations and conservation actions but are limited by scale and precision, particularly for small protected areas with narrow elevational gradients and with long-lasting legacies of past land use and disturbance histories. Homogenization across ecological sites within protected areas from repeated and widespread disturbances is particularly significant because it further reduces site resilience, including spatial resilience (sensu Chambers et al., 2019), and site invasion resistance as fewer contiguous stands of perennial vegetation remain after each successive disturbance event. We observed a substantial increase in the two invasive annual grasses over time following fire, and no apparent response to fire by the two native perennial bunchgrasses, confirming the anticipated fire resiliency of all four of these species summarized by Miller et al. (2013) and previously demonstrated for the perennial bunchgrasses within our study area by Reed-Dustin et al. (2016) and Ellsworth et al. (2020). However, the two annuals also increased with burn severity, most notably medusahead. The ubiquity of cheatgrass across the study area created a ceiling to the occurrence model component (hurdleat-zero), visualized in Figures 3, 7, such that the cheatgrass fire response resulted in an initial reduction in abundance (cover) immediately following fire and a rapid recovery to pre-burn abundance. In contrast, medusahead apparently was able to invade new, severely burned areas quickly (even during the first year after fire in 2012) and did not exhibit any apparent reduction in post-burn cover such that overall trend did not reflect a short-term reduction after fire as did cheatgrass. This contrast in response between the two annual grasses is striking, especially when the differences in response to water deficit are also considered.
Medusahead exhibited a positive association with mean water deficit during October-May prior to each survey, whereas cheatgrass exhibited a negative association with deficit. This finding is consistent with the different life-history traits of the two species. Cheatgrass depends more on early-season germination in late fall and during mild winter conditions and exhibits early senescence during dry spring conditions (Bradford and Lauenroth, 2006;Chambers et al., 2007;Mangla et al., 2011). Medusahead can withstand higher water deficit during this period in part because of later germination, moisture-retaining thatch build-up from previous growing seasons, and is able to utilize late season precipitation before senescence (Mangla et al., 2011;Nafus and Davies, 2014). The accumulation of medusahead thatch reinforces a positive feedback loop of soil conditioning, favorable soil moisture, and fine fuel build-up that reduces invasion resistance and restoration (e.g., native reseeding) potential (Davies, 2008;Mangla et al., 2011;Perkins and Nowak, 2013;Nafus and Davies, 2014). The conservation implications of these medusahead characteristics are notable, suggesting that medusahead has considerable invasibility in the kinds of settings represented by our study area, can withstand drought and fire, and can quickly invade severely burned areas before native perennial grasses and shrubs can recover (Nafus and Davies, 2014). The role of cheatgrass in exacerbating the fire cycle in sagebrush steppe is well-documented but much less so for medusahead (Miller et al., 2013), and our study insights for this species contribute to the growing understanding that this species can exhibit very different ecological dynamics than cheatgrass (Miller et al., 2013). In general, our study results exemplify the precarious conservation status of the John Day Fossil Beds National Monument and other similar low-elevation sagebrush steppe protected areas. At least two species of non-native invasive annual grasses, cheatgrass, and medusahead, and increasingly a third species, wiregrass (Jones et al., 2018;Nicolli et al., 2020; Ventenata dubia), are highly aggressive competitors across large portions of low-elevation aridic sagebrush steppe, including steppe protected areas, that are capable of driving thresholdtype state transitions to novel ecological conditions persisting vis a vis positive feedback loops (D'Antonio and Vitousek, 1992;Suding et al., 2004;Miller et al., 2013). We documented a complete removal of woody shrub cover that had begun prior to our study and was exacerbated by the 2011 fire ( Supplementary Figure 1.8). This loss of woody overstory, the low cover of deep-rooted perennial forbs, and the extent of annual grass invasion indicates that the majority of the Clarno study area has crossed the threshold into the annual grass state (state 3, Figure 2), with only steep and primarily northfacing slopes retaining substantial native perennial vegetation. The ability of land managers to reverse these kinds of state transitions is low both because aridic low-resilience sagebrush ecosystems exhibit low productivity and high inertia and because the required restoration inputs are technologically challenging and expensive (e.g., James et al., 2013). Therefore, these kinds of sites are often ranked low in priority, especially within the context of sage grouse (Centrocercus urophasianus) habitat management (e.g., Chambers et al., 2019). However, the slopes where tall-stature perennial bunchgrasses (specifically wheatgrass and needlegrass in our study area) persist with sufficient cover (e.g., ≥25%) does create a foothold for land managers when broader conservation objectives are pursued. Despite the species-specific fire resilience of wheatgrass and needlegrass exhibited in our study and elsewhere (Miller et al., 2013;Reed-Dustin et al., 2016;Ellsworth et al., 2020;e.g., Ellsworth and Kauffman 2010), the general resilience and invasion resistance of these ecological sites are so low that protection of remaining bunchgrass slopes from fire is reinforced as a core objective that is widely recommended (e.g., Chambers et al., 2019) and already outlined in Monument conservation plans (Rodhouse et al., 2014;Hoh et al., 2015;Shinderman et al., 2020). Protection of the extant contiguous perennial bunchgrass communities on both north-and south-facing slopes would contribute to diversity of bunchgrass community types (i.e., needlegrass and wheatgrass) although bunchgrass-dominant north slopes are more abundant (e.g., Figure 7) and are more likely to persist. However, over time, with sufficient operational capacity, it may also be possible to restore bunchgrass slopes outward from persistent bunchgrass stands, possibly even on hotter south aspects. However, without sufficient capacity, a triage strategy will be required so that protection of remaining quasi-intact bunchgrass slopes is prioritized (Chambers and Wisdom, 2009). The largest contiguous stands should be highest priority. We suggest that our strategy of integrating empirical data with conceptual and statistical models about state transitions is an effective sciencebased approach that could be replicated across the protected-area network to accelerate the operational application of resilience and resistance concepts to sagebrush steppe conservation, including integration into fire management operations.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
TR and KI designed the study. TR collected the data and wrote the manuscript. LB analyzed the data. KI and LB provided substantial contributions to the manuscript. All authors edited the manuscript and agreed to its content.