Plant functional groups and species contribute to ecological resilience a decade after woodland expansion treatments

Woody plant expansions are altering ecosystem structure and function, as well as fire regimes, around the globe. Tree-reduction treatments are widely implemented in expanding woodlands to reduce fuel loads, increase ecological resilience, and improve habitat, but few studies have measured treatment outcomes over long timescales or large geographic areas. The Sagebrush Treatment Evaluation Project (SageSTEP) evaluated the ecological effects of prescribed fire and cut-and-leave treatments in sagebrush communities experiencing tree expansion in North American cold desert shrublands. We used 10 yr of data from the SageSTEP network to test how treatments interacted with pre-treatment tree dominance, soil climate, and time since treatment to affect plant functional groups and dominant species. Non-sprouting shrub (Artemisia spp.), sprouting shrub, perennial graminoid, and annual grass responses depended on tree dominance and soil climate, and responses were related to the dominant species’ life-history traits. Sites with warm and dry soils showed increased perennial graminoid but reduced Artemisia shrub cover across the tree dominance gradient after prescribed burning, while sites with cool and moist soils showed favorable post-burn responses for both functional types, particularly at low to moderate tree dominance. Cut-and-leave treatments sustained or increased native perennial plant functional groups and experienced smaller increases in exotic annual plants in both soil climates across the tree dominance gradient. Both treatments reduced biocrust cover. Selecting appropriate tree-reduction treatments to achieve desired longterm outcomes requires consideration of dominant species, site environmental conditions, and the degree of woodland expansion. Careful selection of management treatments will reduce the likelihood of undesirable consequences to the ecosystem.


INTRODUCTION
Dryland regions are experiencing multiple stressors, including woody species expansion (Barger et al. 2011, Archer et al. 2017, Stevens et al. 2017, plant invasions, and climate warming, which lead to changes in wildfire frequency, severity, and extent (D'Antonio and Vitousek 1992, Brooks et al. 2004). Woody plant expansion is an ongoing process driving ecosystem change in a range of environments, including savannas of Africa, Australia, and South America (Lunt et al. 2010, Wigley et al. 2010, Silva and Anand 2011, and grasslands and shrublands of North America (Hibbard et al. 2001, Romme et al. 2009). These woodland expansions have consequences for hydrologic processes, nutrient cycling, fire regimes, and biological diversity (Hibbard et al. 2001, Parr et al. 2012, Ratajczak et al. 2012, Smit and Prins 2015, Pierson and Williams 2016. In North American drylands, the greatest changes in aboveground biomass from woody species expansion have been attributed to species of juniper and pine (Barger et al. 2011).
Woodland expansion by junipers (Juniperus occidentalis, Juniperus osteosperma) and pinyon pine (Pinus monophylla, Pinus edulis) is contributing to sagebrush ecosystem decline in the cold deserts of North America (Davies et al. 2011, Coates et al. 2017. Expansion of pinyon-juniper woodlands has been attributed to shifting climatic conditions, historic overgrazing, settlement-era tree harvest, and altered fire regimes (reviewed in Miller et al. 2019). In general, pinyon-juniper woodland expansion has resulted in increased risk of high severity fires due to progressive increases in woody fuels (Romme et al. 2009, Young et al. 2015. In addition, the extent and connectivity of sagebrush habitat has been reduced for many sagebrush-dependent wildlife species, including the Greater sage-grouse (Coates et al. 2017).
Tree-reduction treatments in pinyon-juniper woodlands, including mechanical treatments (such as cut-and-leave and mastication) and prescribed fire, have been widely implemented to reduce woody fuels and improve sagebrush-dependent wildlife habitat (Miller et al. 2019;Reinhardt et al. 2020). Tree-reduction treatment types differ in severity and can lead to long-term changes in ecosystem structure and function (Williams et al. 2018, Davies and. Effective ecosystem management requires identifying where on the landscape particular treatments achieve desired objectives while minimizing unanticipated consequences. Information on ecological resilience to disturbance (i.e., resilience) and resistance to invasive annual grasses (i.e., resistance) can be used to assess likely responses to treatments and prioritize areas for management (Maestas et al. 2016, Chambers et al. 2019a. In sagebrush shrublands, resilience and resistance change along environmental gradients and are strongly influenced by climate, topography, and soils , b, 2019a, b, Brooks et al. 2016. Soil climate (i.e., soil and temperature and moisture regimes) can be used to predict resilience and resistance following tree-reduction treatments, as sites in cool and moist regimes typically have greater resilience and resistance than those in warm and dry regimes , b, Roundy et al. 2018. Disturbances and management treatments can increase resource availability, promoting native species establishment and growth but also providing niche opportunities for invasive annual plants (Leffler andRyel 2012, Roundy et al. 2018). Native species' functional traits and life-history attributes are primary determinants of response to disturbance and treatments , Spasojevic et al. 2016) as well as resistance to invasive annual grasses within a given soil climate regime (Brooks et al. 2016). Effective competitors with exotic invasive plants are usually native or desirable plant species, such as perennial grasses and forbs that use similar resource pools and have resource use patterns that coincide with spatial or temporal aspects of the establishment and growth of the invasive plant (Leffler and Ryel 2012). To date, we lack a detailed examination of the long-term recovery of different plant functional groups after tree-reduction treatments and their relative contributions to resilience and resistance.
The responses of functional groups and species to management treatments can be predicted based on their life-history and physiological traits. In sagebrush shrublands, the dominant shrubs are Artemisia species that are fire-intolerant, obligate seeders and their reestablishment is highly dependent on environmental conditions (Shriver et al. 2018(Shriver et al. , 2019. Resprouting shrubs can regrow quickly and are typically more abundant on cooler and moister sites due to higher productivity (fuels) and more frequently historical fires (Pausas and Bradstock 2007, Davies et al. 2012, Spasojevic et al. 2016. Native perennial grasses and forbs vary in terms of response to fire and other disturbances, but typically increase in cover due to resprouting and seedling establishment (e.g., Ross et al. 2012, Bates et al. 2017b, Williams et al. 2017). Characterized by large seedbanks (Koniak andEverett 1982, Allen andNowak 2008), native annual forbs establish rapidly after fire, but are often outcompeted as perennnials establish and increase in abundance (Miller et al. 2013). In addition, lichen and moss biocrusts are highly sensitive to both fire and surface disturbances and can take decades to recover following management actions (Condon andPyke 2018, Root et al. 2020). Despite this general understanding of functional group and species responses to treatments, differences between climatic regimes and implications for longer-term resilience and resistance have not been investigated.
The long-term Sagebrush Treatment Evaluation Project (SageSTEP) provides a unique opportunity to evaluate the contributions of plant functional groups and species to resilience and resistance following tree-reduction treatments over both climatic and tree expansion gradients. SageSTEP was established in 2006 in part to evaluate the effectiveness of prescribed fire and mechanical treatments on sagebrush plant communities experiencing tree expansion (McIver and Brunson 2014). The multi-site experimental framework includes study locations that are distributed across an abiotic gradient of soil climate regimes and that encompass a tree expansion gradient. Ten or more years of post-treatment data now exist, enabling analyses of long-term plant community responses to treatments across the network.
We evaluated how dominant plant functional groups, species, and biocrusts responded to treereduction treatments over soil climate and tree dominance gradients. We predicted that (1) native perennial plant species have variable responses that can be attributed to their life-history traits, (2) native annual plants have ephemeral responses, particularly to prescribed fire, (3) exotic annual plants respond positively to treatment, especially prescribed fire, in warmer and drier site conditions, (4) and lichen and moss biocrusts respond negatively to both treatments. We relate the responses of the individual plant functional types and species to their adaptations to environmental conditions and disturbance and discuss the implications for prescribed fire and mechanical treatments.

Study area
We focused our analysis on 10 sites from the SageSTEP network (http://sagestep.org) in the Great Basin, USA, where sagebrush ecosystems are experiencing pinyon and/or juniper expansion (Table 1) Climatic and edaphic characteristics influence the ecological resilience of plant communities to disturbance and resistance to plant invasion (Hoopes and Hall 2002, Pollnac et al. 2012, Chambers et al. 2019a. Soil temperature and moisture regimes (Soil Survey Staff 1999) are based on annual and/or seasonal conditions of the belowground climate of a site and are a useful indicator of resilience and resistance (Chambers et al. , 2019a. Therefore, we grouped our sites by soil temperature and moisture regimes in analyses. Across the range of site variability for big sagebrush-dominated ecosystems experiencing tree expansion, mesic/aridicxeric soils are warmer and drier, and these sites typically support communities dominated by Wyoming big sagebrush (Artemisia tridentata ssp. wyomingensis). Sites with mesic/aridic-xeric soils include Bridge Creek, Greenville Bench, Marking Corral, Onaqui, Scipio, and South Ruby. Frigid/ xeric soils are cooler and moister, and sites with frigid/xeric soils are typically dominated by v www.esajournals.org mountain big sagebrush (A. tridentata ssp. vaseyana). Sites with frigid/xeric soils include Blue Mountain, Devine Ridge, Seven Mile, and Walker Butte. In situ measurements confirmed that sites classified as mesic/aridic-xeric had warmer soil temperatures in all seasons and were drier in spring and winter than sites classified as frigid/xeric (Roundy et al. 2020). Sites were fenced prior to treatment to exclude livestock. Further details on site history, treatments, and locations are presented in McIver and Brunson (2014).

Experimental design
At each of the 10 sites, there were three plots assigned the following treatments: (1) control, (2) prescribed fire, and (3) mechanical tree reduction by cut-and-leave, where trees >2 m tall were cut by chainsaw and left on site. Plots were 8-20 ha in size and contained between 12 and 22 sampling subplots of 30 9 33 m. Subplots within a plot covered a gradient of tree expansion from low to high tree cover. Bridge Creek, Marking Corral, Onaqui, and Walker Butte were treated in 2006, while Blue Mountain, Devine Ridge, Greenville Bench, Scipio, and Seven Mile were treated in 2007. South Ruby was treated in 2009, so data presented from this site does not include year 10. Because treatments were implemented in a staggered-start design, the calendar year varies across sites for each post-treatment year, allowing us to better isolate the effect of time since treatment from the effect of inter-annual weather variation associated with specific calendar years (Loughin 2006). Although every attempt was made to ensure that the functional group composition of the treatment plots was highly similar, when averaged across sites within a soil temperature and moisture regime, pretreatment functional group composition sometimes differed among treatments.

Data collection
Vegetation within each subplot was sampled along five transects using the line-point intercept method (Herrick et al. 2017). Vegetation measurements included canopy cover of shrub species, foliar cover of herbaceous species, and ground cover of lichen/moss biocrust; data were summarized by functional group. The shrub functional groups were Artemisia shrubs (all nonsprouting sagebrush species, calculated from the v www.esajournals.org sum of A. arbuscula, A. nova, A. rigida, and A. tridentata cover) and other shrubs (the sum of cover for all other shrub species, which are weak to strong re-sprouters). Herbaceous functional groups were perennial graminoids, annual grasses (represented predominantly by the exotic Bromus tectorum), perennial forbs, native annual forbs, and exotic annual forbs, which were calculated by summing the cover of all species within each group. Tree canopy cover was measured as the sum of aerial crown cover (A) for each tree >0.5 m tall in the subplot, using the formula A = p (D1 9 D2)/4, where D1 is the longest crown diameter and D2 is the perpendicular crown diameter. For pre-treatment measurements, tree cover was measured in the year of treatment implementation, whereas all other pretreatment vegetation measurements were collected in either the treatment year or the previous year. Because the full set of subplots was not sampled at every site in every year, we used data from five years that received full sampling: pretreatment, 1-yr post-treatment, 3-yr post-treatment, 6-yr post-treatment, and 10-yr post-treatment.

Statistical analyses
All statistical analyses were performed using R version 3.5.1 (R Core Team 2018). To determine whether annual weather conditions displayed temporal trends across sites that could affect vegetation responses by year-since-treatment (henceforth "year"), we calculated deviations from normal for each post-treatment year at each site by subtracting yearly values from 30-yr normals using 800-m gridded data (PRISM Climate Group 2019); we then ran an ANOVA with each weather variable as a dependent variable and year as the independent variable, followed by Tukey's honest difference tests to determine significant pairwise contrasts between the years. Weather variables were total precipitation (mm), minimum temperature (°C), mean temperature (°C), and maximum temperature (°C). A significant trend was found only for mean and minimum temperatures in year 10, which was cooler than other years (Appendix S1: Fig. S1).
We calculated the pre-treatment Tree Dominance Index (TDI) using the formula (tree cover/ (tree + shrub + perennial graminoid cover)), which was modified from Roundy et al. (2014a) to include all perennial graminoids. This calculation standardized tree dominance across sites with varying potential productivity by comparing tree cover to the total dominant perennial vegetation cover supported by a site. To determine the effects of treatments, year, and TDI on vegetation responses, we ran a linear mixed effects model for each response variable in each of the two soil temperature and moisture regimes separately. Treatment, year, TDI, and all two-way and three-way interactions among these variables were included as fixed effects, while site was included as a random effect in order to account for the wide variation in pre-existing vegetation cover and environmental conditions across the site network and to minimize the influence of autocorrelation among the subplots sampled at every site. We applied the Kenward-Roger approximation to ensure a conservative Type I error probability in the presence of unbalanced sample sizes and unequal variance among groups (Volaufova 2009). Models were fit using the lme4 package (Bates et al. 2015) and tested for significance using the lmerTest package (Kuznetsova et al. 2017). We refer to these as the temporal models henceforth to distinguish them from other analyses (Tables 2-4). We then examined pairwise differences with Tukey's honest difference tests using the estimated means between treatment-year combinations with the emmeans package (Lenth 2018); pairwise test results are presented in Appendix S1: Tables S3,  S4. Finally, to determine whether the soil temperature and moisture regime significantly affected treatment outcomes, we ran linear mixed effects models for year 10 cover and for the net change between pre-treatment and year 10 cover. In these models, treatment, TDI, soil temperature and moisture regime, and all two-way and threeway interactions among these variables were fixed effects, and site was a random effect. We refer to these as the soil climate models (Table 5).
We used the following vegetation cover groups as dependent variables in all analyses: four dominant functional groups (shrubs in the sagebrush genus Artemisia, other shrubs, perennial graminoids, annual grasses), and four sub-dominant functional groups (perennial forbs, native annual forbs, exotic annual forbs, and lichen/moss biocrust). We also analyzed the top three dominant shrub species and top three dominant perennial v www.esajournals.org graminoid species (by mean cover) for the respective soil temperature and moisture regime. We did not separate A. tridentata into subspecies because of uncertainty of field identification to this taxonomic level; evidence suggests local adaptation of populations to site conditions may be more  important than subspecies in explaining variation in A. tridentata performance (Lazarus et al. 2019). Dominant species that were not shared by both of the soil temperature and moisture regimes were not included in the models that tested for an interaction of soil temperature and moisture regime with treatment and TDI. Dependent variables were transformed with either logarithmic or square root transformations as necessary to meet linear model assumptions. In some cases, transformation was not sufficient to correct for severely zero-inflated and/or skewed data; therefore, these cases are depicted graphically in the figures but statistical results have not been presented in the tables. These included temporal models for perennial graminoids, native annual forbs, exotic annual forbs, lichen/moss biocrust, Purshia tridentata, and Achnatherum hymenoides at mesic/aridic-xeric sites; temporal models for exotic annual forbs, P. tridentata, Chrysothamnus viscidiflorus, and Poa secunda at frigid/xeric sites; soil climate models for year 10 cover of native annual forbs, C. viscidiflorus, and P. tridentata; and the soil climate model for the change in cover from pre-treatment to year 10 for P. tridentata.

Dominant species and treatment longevity
At mesic/aridic-xeric sites, the three shrub species with the highest mean cover across the study duration were A. tridentata (5.73%), C. viscidiflorus (1.81%), and P. tridentata (0.89%). The perennial graminoids with the highest cover were Pseudoroegneria spicata (10.16%), P. secunda (4.83%), and A. hymenoides (1.52%). The annual grass with mean cover >1% was Bromus tectorum (7.26%), and the exotic annual forbs with mean cover >1% were Alyssum desertorum (2.82%) and Ceratocephalum testiculatum (1.02%). No perennial forb or native annual forb species had cover >1% across these sites. At frigid/xeric sites, the shrub species with the highest mean cover were A. tridentata, P. tridentata, and C. viscidiflorus (6.34%, 2.53%, 1.54%, respectively). The perennial graminoids with the highest cover at these sites were Festuca idahoensis (10.70%), P. secunda (5.31%), and Elymus elymoides (2.87%). The annual grass with mean cover >1% was B. tectorum (3.65%). No forb species had cover >1% across these sites. Dominant species by treatment type in year 10 for each site are in Appendix S1: Tables S1, S2. In year 10 of the study, mean cover of pinyon and juniper trees in both the prescribed fire and cut-and-leave treatments remained below 2% (reduced from pre-treatment mean cover of 14.5% and 23.6%, respectively), while mean cover in the control plots had changed less, from 17.4% pre-treatment to 15.9% in year 10.

Functional group responses
Ten years after treatment, Artemisia shrub cover was generally lower in burned and higher in mechanical plots than in control and pre-treatment plots ( Fig. 1; Appendix S1: Table S3). The same effects were found for the differences between pre-treatment and year 10 cover. This treatment effect was generally consistent across the soil temperature and moisture regimes (Tables 2, 5). A significant interaction of treatment with TDI and year existed for both mesic/ aridic-xeric and frigid/xeric sites (Table 2). For both soil temperature and moisture regimes, the negative influence of burning on Artemisia shrub cover decreased with increasing TDI due to lower initial Artemisia shrub cover at higher TDI. The positive effect of the mechanical treatment on Artemisia shrub cover at both mesic/aridic-xeric and frigid/xeric sites was not influenced by TDI.
For other shrubs, both burning and mechanical treatment generally led to increased cover in year 10 compared to control plots and pre-treatment values, but significant interactions occurred among TDI, treatment, and year (Table 2), as well as among TDI, treatment, and soil temperature and moisture regime (Table 5). Burning had a positive effect on year 10 cover of other shrubs for both soil regimes, but this effect was greater at lower TDI in frigid/xeric sites. Other shrub cover also increased regardless of TDI in frigid/ xeric sites when mechanically treated (Fig. 1). In contrast, the mechanical treatment only had a positive effect on other shrub cover at higher TDI in mesic/aridic-xeric sites. Treatment effects on other shrub cover were generally not as apparent  Year 10 panels: response (predicted means with 95% confidence intervals) of dominant functional groups (Artemisia shrubs, other shrubs, perennial graminoids, and annual graminoids) to treatment (control, prescribed fire, and mechanical) and tree dominance index (TDI) in year 10. Change panels: difference in cover between year 10 and pre-treatment; values above the horizontal black dashes at 0 on the y-axis represent an increase in cover since pre-treatment and values below represent a decrease in cover since pre-treatment. Vertical gray dashes represent the first and third quartile of TDI across sites.
v www.esajournals.org 9 January 2021 v Volume 12(1) v Article e03325 in mesic/aridic-xeric sites until year 10 (Appendix S1: Table S3). At frigid/xeric sites, cover of other shrubs was lower in burned plots than in control plots in all but year 6, and cover of other shrubs was higher in mechanically treated plots than in control plots from year 6 onward (Appendix S1: Table S3). For perennial graminoid cover, treatment effects for both year 10 and the change in cover from pre-treatment to year 10 were influenced by soil temperature and moisture regime (Table 5). At frigid/xeric sites, there was a three-way interaction among treatment, TDI, and year, but the model for mesic/aridic-xeric sites could not be adequately fit (Table 2). Mesic/aridic-xeric sites showed increased perennial graminoid cover following both burn and mechanical treatments; this increase was greater with increasing TDI in mechanical treatments, but did not vary with TDI in prescribed burn treatments (Fig. 1). However, due to variability among sites there were no significant treatment effects within a given year (Appendix S1: Table S3). At frigid/xeric sites, mechanical treatment had a significant positive effect on perennial graminoid cover from year 3 onward ( Fig. 1; Appendix S1: Table S3). Both burning and mechanical treatment had greater differences between pre-treatment and year 10 cover than initial and year 10 cover in control plots (Fig. 1). This increase in perennial graminoid cover over time in the mechanical plots was greater at higher TDI, likely due to the low herbaceous cover present initially in locations with more mature woodland, and the lower cover of shrubs to compete with for resources in these locations. Site explained over 50% of the variation in perennial graminoid cover in the temporal model, which was the greatest effect of site in any of the functional group temporal models (Table 2). Although the effects of TDI, treatment, and year were significant at frigid/xeric sites, perennial graminoid cover was generally weakly characterized by these variables (28% of variance explained).
Annual grass cover was significantly greater in treated plots by year 10 compared to control and pre-treatment plots, but patterns over time were dependent upon TDI and soil temperature and moisture regime (Table 5; Appendix S1: Figs. S2, S3). At mesic/aridic-xeric sites, annual grass cover showed a large increase 10 yr after burning, an effect which did not vary with TDI. In contrast, the increase in annual grass cover at frigid/xeric sites after burning was smaller and showed a negative relationship with TDI. Mechanically treated plots had similar trends in both soil temperature and moisture regimes; the positive treatment effects were initially more pronounced at higher TDI, but then became apparent at lower TDI in year 10. Mesic/aridic-xeric sites showed a larger and earlier (year 3 vs. year 6) increase in annual grass cover following mechanical treatment than frigid/xeric sites. The change in annual grass cover from pre-treatment to year 10 differed significantly among all three treatments for both soil regimes. Annual grass cover increased beyond pre-treatment values only in treated plots. Burned plots had significantly greater increases from pre-treatment to year 10 than mechanical plots at high TDI in mesic/aridic-xeric sites and low TDI in frigid/ xeric sites (Fig. 1).
Perennial forb cover in year 10 was higher in treated plots than in controls, but temporal patterns in perennial forb cover were not well explained by our model parameters (Table 2). Significant interactions existed for soil temperature and moisture regime with both TDI and treatment for the change in perennial forb cover from pre-treatment to year 10 (Table 5). Mesic/ aridic-xeric sites showed modest increases in perennial forb cover 10 yr after treatment, primarily at higher TDI in burned plots (Fig. 2). Frigid/xeric sites showed a large increase in perennial forb cover following both burn and mechanical treatments regardless of TDI, but due to variability among sites, there were no significant treatment effects within a given year (Appendix S1: Table S3).
For native annual forbs, there was a significant interaction of soil temperature and moisture regime with TDI and treatment for the change in cover since pre-treatment (Table 5), but the model for year 10 cover could not be adequately fit. At mesic/aridic-xeric sites, native annual forb cover showed no to slight increases across the TDI gradient. At frigid/xeric sites, native annual forb cover in burned plots was greater than in both control and mechanical plots in all posttreatment years (Appendix S1: Table S3, Fig. S5). There was higher native annual forb cover in frigid/xeric burned plots in year 10 relative to prev www.esajournals.org   Fig. 2. Year 10 panels: response (predicted means with 95% confidence intervals) of sub-dominant functional groups (perennial forbs, native annual forbs, exotic annual forbs, lichen/moss biocrust) to treatment (control, prescribed fire, and mechanical) and tree dominance index (TDI) in year 10. Change panels: difference in cover between year 10 and pre-treatment; values above the horizontal black dashes at 0 on the y-axis represent an increase in cover since pre-treatment and values below represent a decrease in cover since pre-treatment. Vertical gray dashes represent the first and third quartile of TDI across sites. v www.esajournals.org treatment values, and this response increased with increasing TDI (Fig. 2). In mechanical plots, an increase in native annual forb cover since pretreatment was only apparent at high TDI (Fig. 2).

Frigid/Xeric (Cool/Moist) Sites
Exotic annual forb cover could not be fit to a temporal model, but soil climate models revealed significant two-way interactions between treatment and TDI, and between treatment and soil temperature and moisture regime, for the change in exotic annual forb cover from pre-treatment to year 10 (Table 5). At mesic/aridic-xeric sites, exotic annual forb cover increased in control plots primarily at low to moderate TDI; the reverse pattern occurred in mechanically treated plots, where increases in exotic annual forb cover were greater with increasing TDI. In burned plots at mesic/aridic-xeric sites, exotic annual forb cover increased across the entire TDI gradient (Fig. 2). At frigid/xeric sites, control plots had small increases in exotic annual forb cover at low TDI, and burned plots had increased cover at all TDI, but these responses were of a much smaller magnitude than at mesic/aridic-xeric sites (Fig. 2). Exotic annual forb cover did not increase at any TDI in mechanically treated plots at frigid/xeric sites (Fig. 2).
Lichen/moss biocrust cover was reduced by treatments at both mesic/aridic-xeric and frigid/ xeric sites. The model for mesic/aridic-xeric sites could not be adequately fit. Mesic/aridic-xeric sites did however exhibit a decrease in biocrust cover following both burn and mechanical treatments and this decrease was greatest for low TDI on burned plots (Fig. 2). At frigid/xeric sites, there was a three-way interaction among treatment, TDI, and year. Burned plots at frigid/xeric sites had lower lichen/moss cover than control plots in all post-treatment years, and burned plots had lower cover than mechanical plots in years 1, 3, and 10 (Appendix S1: Table S3, Fig. S3). Mechanical plots had significantly less lichen/moss than control plots, particularly at high TDI in year 3 onward (Appendix S1: Fig. S3). Changes in lichen/moss cover from pretreatment to year 10 differed between control plots and treated plots; cover increased modestly in control plots but decreased in treated plots. This decrease in lichen/moss cover since pretreatment was greatest at low to moderate TDI in burned plots and greatest at high TDI in mechanical plots.

Dominant species responses
The three dominant shrub species were the same for the frigid/xeric and mesic/aridic-xeric sites. An interaction among TDI, treatment, and soil temperature and moisture regime existed for the change in A. tridentata cover between pretreatment and year 10 (Table 5). In both regimes, A. tridentata cover was lower than pre-treatment in year 10 in burned plots at low to medium TDI, higher in mechanical plots across the full range of TDI, and did not change in the control plots (Fig. 3). However, frigid/xeric sites showed a greater positive response to burning at high TDI than did mesic/aridic-xeric sites. Similarly, A. tridentata cover was lower in burned plots than in control plots in all post-treatment years, although in years 6 and 10, this was primarily the case only at low to moderate TDI in frigid/xeric sites (Appendix S1: Table S4, Fig. S7). In mechanical plots, the change in A. tridentata cover was positive across the entire range of TDI, although this effect was apparent earlier in the post-treatment timeline at frigid/xeric compared to mesic/aridicxeric sites (year 3 vs. year 6; Appendix S1: Table S4, Figs. S6, S7).
Chrysothamnus viscidiflorus had higher cover in treated than control plots at mesic/aridic-xeric sites, but these differences were not apparent until year 6 and were primarily seen at low to moderate TDI (Appendix S1: Table S4, Fig. S6). At frigid/xeric sites, C. viscidiflorus cover increased in treated plots at low to moderate TDI, especially in burned plots. The change in C. viscidiflorus cover from pre-treatment to year 10 was significantly different between control and burned plots and between control and mechanical plots at mesic/aridic-xeric sites, where C. viscidiflorus cover increased modestly in both burn and mechanical plots ( Fig. 3; Appendix S1: Table S4). There was no change in P. tridentata cover in treated or control plots at mesic/aridic-xeric sites, where cover of this species was very low overall (Fig. 3). At frigid/xeric sites, P. tridentata cover increased only in mechanical plots (Fig. 3).
At mesic/aridic-xeric sites, P. spicata cover in both burned and mechanical plots was higher than in control plots from year 3 onward, primarily at moderate to high TDI (Appendix S1: Table S4, Fig. S8). In both treatment types, year 10 P. spicata cover was higher than pre-treatment

Frigid/Xeric (Cool/Moist) Sites
Difference in cover (%) Difference in cover (%) Fig. 3. Year 10 panels: response (predicted means with 95% confidence intervals) of dominant shrub species (Artemisia tridentata, Chrysothamnus viscidiflorus, Purshia tridentata) to treatment (control, prescribed fire, and mechanical) and tree dominance index (TDI) in year 10. Change panels: difference in cover between year 10 and pre-treatment; values above the horizontal black dashes at 0 on the y-axis represent an increase in cover since pre-treatment and values below represent a decrease in cover since pre-treatment. Vertical gray dashes represent the first and third quartile of TDI across sites. Year 10 panels: response (predicted means with 95% confidence intervals) of perennial graminoid species (Pseudoroegneria spicata, Poa secunda, Achnatherum hymenoides at mesic/aridic-xeric sites; Festuca idahoensis, Poa secunda, Elymus elymoides at frigid/xeric sites) to treatment (control, prescribed fire, and mechanical) and tree dominance index (TDI) in year 10. Change panels: difference in cover between year 10 and pre-treatment; values above the horizontal black dashes at 0 on the y-axis represent an increase in cover since pre-treatment and values below represent a decrease in cover since pre-treatment. Vertical gray dashes represent the first and third quartile of TDI across sites. v www.esajournals.org cover, and this increase was consistent across the gradient of TDI (Fig. 4). Cover of A. hymenoides increased in burned plots at low to moderate TDI (Fig. 4). At frigid/xeric sites, burning reduced F. idahoensis cover compared to mechanical or control plots (Fig. 4; Appendix S1: Table S4), although the model explained very little of the variation in cover (Table 4). The change in cover from pre-treatment to year 10 differed significantly between burned and control plots. In burned plots, F. idahoensis cover decreased at low TDI, but increased at high TDI, while in control plots F. idahoensis cover increased across the entire TDI gradient. In contrast to F. idahoensis, both burning and mechanical treatment led to higher cover of E. elymoides in year 10 (Appendix S1: Table S4), particularly at moderate to high TDI. Year 10 cover of E. elymoides in treated plots was also higher than pre-treatment cover, particularly in burned plots and at moderate to high TDI (Fig. 4). In contrast, there was no change in E. elymoides cover from pre-treatment to year 10 in the control plots.
Models testing interactions among TDI, treatment, and soil temperature and moisture regime showed significant interactions among these three parameters for year 10 P. secunda cover and the change in cover from pre-treatment to year 10 (Table 5). At frigid/xeric sites, burning led to slightly higher P. secunda cover in year 10 compared to control plots at low TDI, and slightly higher cover compared to pre-treatment at low to moderate TDI (Fig. 4). In contrast, mechanical and control plots showed similar levels of increasing cover in year 10. At mesic/aridic-xeric sites, both burning and mechanical treatments resulted in little to no change in P. secunda cover compared to control and to pre-treatment values (Fig. 4), with no effect of TDI (Table 3). In addition, treatment effects were not statistically significant in any given year (Appendix S1: Table S4), and control plots did not experience an increase in P. secunda cover from pre-treatment to year 10 as was seen in the frigid/xeric sites.
Notes: TDI, tree dominance index. Arrows and dashes indicate direction and relative magnitude of effect within the soil regime type: dash = minimal change, up arrow = increase in cover, down arrow = decrease in cover, two arrows = greater effect, three arrows = greatest effect.

DISCUSSION
Increasing resilience to wildfire and resistance to invasive annual plants by reducing woody expansion can only be achieved when the dominant functional groups and diverse native species recover. We found that despite variability in functional group and species identities and their relative abundances among sites, major differences existed in their responses depending on soil climatic regime, treatment type, and/or tree dominance ( Table 6). The differences in these responses can largely be explained by the adaptations of individual species comprising the functional group to environmental conditions and disturbance effects. Our results emphasize the need to consider not only soil climatic regimes and the degree of tree expansion but also different plant functional groups and species when evaluating sites for appropriate treatments.

Native perennial responses depended on soil climatic regime and species life-history
Obligate-seeding shrubs often decrease, while root-sprouting shrubs increase, as effective precipitation and productivity increase and fire-return intervals become shorter (Mermoz et al. 2005, Pausas and Bradstock 2007, Spasojevic et al. 2016. In sagebrush ecosystems, the keystone species are subspecies of Artemisia tridentata that are obligate seeders killed by fire. Ten years after treatment, Artemisia showed little recovery across the tree dominance gradient on mesic/aridic-xeric burned sites. This lack of Artemisia recovery was similar to other studies on warm and dry sites characterized by A. tridentata ssp. wyomingensis (Knudtson et al. 2014, Shriver et al. 2018 and in treated woodlands where burning led to greater long-term sagebrush losses than mechanical treatment . In contrast, Artemisia cover increased in burned plots in frigid/xeric sites at high tree dominance and increased in mechanical plots in both soil climatic regimes regardless of tree dominance. Prescribed burning was focused on treatment plots, and unburned patches adjacent to burned plots likely resulted in persistent seed sources for Artemisia and other species. Also, increased resource availability following tree reduction , Bates et al. 2017b, combined with lower cover of perennial graminoids and B. tectorum in higher tree dominance locations, likely resulted in competitive release of Artemisia . Finally, livestock exclusion after the treatments may have benefited Artemisia. Negative relationships between shrubs and perennial native grass cover in A. tridentata shrublands have been observed during grazing release (Anderson and Inouye 2001) in A. tridentata ssp. wyomingensis shrublands as well as after fire under moderate grazing in A. tridentata ssp. vaseyana shrublands (Harniss andMurray 1973, Hanna andFulgham 2015). Our results indicate a loss of resilience due to low establishment of Artemisia after fire on mesic/aridic-xeric sites, but an increase in resilience over time due to establishment and growth of Artemisia on frigid/xeric sites and on mechanical treatment plots. However, high tree dominance locations were characterized by a relatively high density of immature Artemisia shrubs, and thus neither plant community composition nor sagebrush stand structure reflected that of a mature community .
Resprouting shrubs (other shrubs) typically showed positive responses across the tree dominance gradient following both treatments. Chrysothamnus viscidiflorus, which produces large quantities of wind-dispersed seeds (Tirmenstein 1999, Miller et al. 2013, increased following treatment in both soil climatic regimes but was especially abundant at high tree dominance in mechanical plots in frigid/xeric sites. Purshia tridentata increased at low tree dominance in frigidxeric sites, reflecting resprouting of some ecotypes, as well as increased plant growth, seed productivity, and seedling establishment of this large-seeded and largely rodent-dispersed species (Vander Wall 1994, Zlatnik 1999a. Early seral resprouting shrub species with high seed production and long-distance dispersal mechanisms increased, providing rapid site stabilization following fire and on high tree dominance sites in both soil climatic regimes. However, effects on resilience were mixed, as C. viscidiflorus and similar wind-dispersed species in the Asteraceae family can preempt resources and prevent later seral species from establishing (Miller et al. 2013).
Perennial graminoids exhibited the most positive responses to the treatments, as observed v www.esajournals.org elsewhere (Bates et al. 2007, 2017b, Ross et al. 2012, Williams et al. 2017, Ernst-Brock et al. 2019. Most perennial graminoid species in sagebrush ecosystems are relatively fire-tolerant (Miller et al. 2013), and greater resource availability following prescribed burning and mechanical treatments likely promoted growth and seed production , Bates et al. 2017b. However, perennial grass responses to tree reduction at high pre-treatment tree dominance can be highly variable (Bates et al. 2013). We found major differences between particular species of perennial graminoids in the two soil climate regimes in response to treatments. At mesic/aridic-xeric sites, Pseudoroegneria spicata, a dominant, fire-tolerant bunchgrass (Wright et al. 1979, Zlatnik 1999b, increased in cover across the tree dominance gradient in response to both treatments. In contrast, Festuca idahoensis, the dominant bunchgrass at frigid/xeric sites, was significantly reduced by burning at low tree dominance where it was most abundant prior to treatment. Low tree dominance areas were characterized by relatively high cover of shrubs and grasses (woody and herbaceous fuels) and likely experienced contiguous burns of high enough severity to negatively affect F. idahoensis, a species with relatively low fire tolerance (Wright et al. 1979, Zouhar 2000, Boyd et al. 2015. Festuca idahoensis has the capacity to establish from seed after fire given seed survival or a seed source (Zouhar 2000), and slight increases at high tree dominance following burning may reflect seedling establishment. Elymus elymoides, an early seral, fire-tolerant bunchgrass (Simonin 2001, Miller et al. 2013) that is recommended for restoration seeding (Monsen and Stevens 2004), increased following both burning and mechanical treatments at frigid/xeric sites. This species produces large numbers of highly germinable, wind-dispersed seeds (Simonin 2001) and the strong increase in E. elymoides at high tree dominance may be due to the lack of competing species on these sites after treatment. Though smaller than the changes observed in treated plots, we also observed increases in perennial graminoid cover in control plots at frigid/xeric sites and attribute this primarily to the exclusion of livestock following the establishment of the study. Rest from grazing over intermediate time periods (5-6 yr) leads to increased cover of perennial graminoids (Davies et al. 2016), and the larger increases at frigid/xeric sites likely reflect these sites' greater potential for perennial grass productivity. Perennial graminoids were the major contributor to resilience over time, but differences in recovery potential existed that could be attributed to species life-history and physiological traits and pre-treatment tree dominance.
Forbs play critical roles in ecosystem functioning by providing food sources for wildlife and maintaining pollinator communities (Dumroese et al. 2016, Bates et al. 2017a. We found that perennial native forb cover increased in both treatments across the tree dominance gradient in both soil climate regimes, with the most pronounced increase on frigid/xeric sites. The most common and abundant perennial forb species were both fire and grazing tolerant, including Phlox spp., Linanthus pungens, Achillea millefolium (Miller et al. 2013), and legumes such as Lupinus spp. (Goergen and Chambers 2009). The increase in perennial forbs undoubtedly increased site resilience, but the dominant forb species and their responses to treatment likely reflected both a legacy of livestock grazing and fire history. High variability in perennial forb response was likely related to the high diversity of functional traits exhibited by this group (Barga et al. 2018) and differences in species composition and relative abundance that reflect variable site conditions and histories (Bates et al. 2017a).
Native annual species responded rapidly but persistence depended on tree dominance and soil climate regime Native annual forbs also provide important food sources for sagebrush obligate species like sage-grouse and help maintain pollinator communities (Dumroese et al. 2016). These species are often a major component of seed banks (Koniak andEverett 1982, Allen andNowak 2008) and can increase immediately after vegetation removal  or fire (Miller et al. 2013). Prescribed fire led to increases in native annual forb cover at high tree dominance in mesic/aridic-xeric sites, but this effect was ephemeral (Appendix S1: Fig. S4) as observed elsewhere (Miller et al. 2013). In contrast, at frigid/xeric sites prescribed burning led to a rapid flush of a different set of native annual v www.esajournals.org forbs in year one that generally persisted and became more abundant on high tree dominance sites through time (Appendix S1: Table S1, Fig. S5). In the absence of competition from residual perennials or other early seral species, we found that native annuals can persist for significant periods of time.

Treatments resulted in persistent increases in exotic invasive annuals
Higher invasive annual grass abundance in relatively warm and dry sites following fire and management treatments is common and can be explained largely by greater climatic suitability for establishment and growth (Brooks et al. 2016). A negative relationship between annual grass cover and perennial native grass cover often exists on mesic/aridic-xeric sites, especially where firetolerant perennial herbaceous species were removed or depleted prior to treatment (Chambers et al. 2007, Roundy et al. 2018. We found that despite increases in perennial grasses over time and an apparent negative relationship between annual and perennial grass cover on mesic/aridic-xeric sites, annual grasses continued to increase over the 10 yr of the study across the tree dominance gradient, especially on burned sites. The relationship of annual grass cover with perennial graminoid cover was less negative on frigid/xeric sites, likely due to more favorable conditions for cheatgrass seedling establishment in association with other species (Chambers et al. 2007, Goergen andChambers 2012). An increase in annual grasses on frigid/xeric mechanical plots may be due to more favorable conditions for seed entrapment and seedling establishment under downed trees (Monty et al. 2013).
Similar to invasive annual grasses, exotic annual forbs (e.g., Sisymbrium altissimum, Descurainia sophia, Alyssum desertorum, Ceratocephalum testiculatum, and Tragopogon dubius) increased following both treatments at mesic/aridic-xeric sites. Exotic annual forb increases were largest following prescribed burning and at high tree dominance, which may be attributed to high seed production, mechanisms for long-distance dispersal by wind or adhesion, and seed banking (Howard 2003). There was less cover and diversity of exotic annual forb species at frigid/xeric sites, presumably due to lower climatic suitability (Miller et al. 2013), and increases after treatment were much less pronounced than at mesic/aridic-xeric sites. High variability in exotic annual grasses and forbs among years reflects species' sensitivity to inter-annual differences in precipitation and temperature in drylands (West andYorks 2002, Munson et al. 2013) and other site characteristics. This sensitivity may explain the observed increases of exotic annual forbs in mesic/aridic-xeric control plots as well; this occurred primarily at lower TDI, where a lesser degree of resource capture by establishing conifers combined with livestock exclusion since the onset of the study may have further allowed exotic annuals to proliferate in climatically favorable years. Low resistance to invasive annuals on mesic/aridic-xeric sites decreased resilience by preempting resources for native species (Leffler and Ryel 2012), increasing fine fuels, and increasing the likelihood of more frequent and extensive fires (Brooks et al. 2004(Brooks et al. , 2016. Lichen and moss biocrusts were poorly adapted to disturbance from tree-reduction treatments Prescribed burning and mechanical treatments reduced the cover of lichen and moss biocrusts in both climate regimes. This negative effect was immediate and generally persisted through time, particularly in burned frigid/xeric sites, where mean biocrust cover was reduced to less than 1% with little to no recovery across the tree dominance gradient after 10 yr. The less severe response to burning in the mesic/aridic-xeric sites may be due to much higher initial cover of biocrusts and thus recovery potential. In mechanically treated plots at frigid/xeric sites, some biocrust recovery was evident in year 10, which was consistent with the lower end of timescales seen for biocrust recovery to physical disturbance (Weber et al. 2016). Loss of biocrust decreases site resilience to disturbance and resistance to invasive plants because biocrusts provide ground cover that resists the establishment of annual grasses (Condon andPyke 2018, Root et al. 2020), facilitates native perennial species that compete with annual grasses (Zhang et al. 2016, Condon and, and reduces ecosystem water loss (Chamizo et al. 2016a, b).

Synthesis and management implications
Although the wide variety in site histories and species across the SageSTEP network may have v www.esajournals.org precluded identifying how certain species and functional groups responded to treatments in this study, we nevertheless found strong network-wide relationships between woody expansion, soil climatic regime, and major plants of interest under common treatment scenarios. Prescribed fire in this system results in a more severe disturbance than mechanical treatments and has the potential to reduce both resilience and resistance. Perennial native grasses showed consistent increases over time across the tree dominance gradient for both soil climatic regimes. However, on mesic/aridic-xeric sites prescribed fire had negative effects on Artemisia shrubs, fire-sensitive perennial herbaceous species, and lichen and moss biocrusts. Also, large increases of exotic invasive annuals occurred (Table 6). Although lesser effects on functional group and species recovery were observed on frigid/xeric sites, Artemisia shrubs and the dominant perennial grass, F. idahoensis, decreased while both native and exotic invasive annual forbs increased. These 10-yr results indicate long-term decreases in resilience and resistance following prescribed fire on mesic/aridic-xeric sites. Changes in functional group composition relative to initial conditions on frigid-xeric sites indicate more positive effects on resilience and resistance. Focusing treatments on sites with low to moderate tree dominance can decrease site availability for opportunistic native and exotic invasive species that can alter successional trajectories and reduce resilience to future fires , Roundy et al. 2018. Low severity fires can minimize effects on species with low fire tolerances (Wright et al. 1979), and leaving unburned patches may help ensure long-term seed sources and increase the probability of establishment of Artemisia and other native species (Ricketts and Sandercock 2016).
Mechanical treatments have less severe effects than prescribed fire, resulting in more rapid shrubland community recovery and lesser effects on resilience and resistance. Ten years after mechanical treatment, tree cover remained low (less than 2%) and perennial grass and Artemisia shrub cover had increased across the tree dominance gradient on both soil climatic regimes (Table 6). Annual grasses also increased, but to a lesser degree than following prescribed fire. While it appears that overall resilience was increased, resistance was reduced, potentially increasing the risk of invasive grass-fueled fires in the future. Mechanical treatments that minimize surface disturbances may lessen the negative effects on lichen and moss crusts and reduce the loss of resistance to invasive annual grasses.
If treatments are desired in areas with higher tree cover or depleted understories, seeding a diverse mix of site-adapted, native species (rather than introduced species) after treatment may increase resistance to cheatgrass and maintain habitat values (Ott et al. 2019, Urza et al. 2019. Favoring existing genotypes that are better-adapted to future conditions, including tolerances to drought, pests, and disturbance, can increase a site's adaptive capacity to climate change (Butler et al. 2012, Finch et al. 2016). In addition, the flexibility to seed when weather conditions are appropriate and in more than a single year, especially on mesic/aridic-xeric sites, can help ensure successful seedings (Hardegree et al. 2018). Finally, appropriate livestock grazing management after treatment is a requisite for maintaining perennial herbaceous species (Bunting et al. 1998).
Long-term data across broad environmental gradients are essential for determining the effects of tree-reduction treatments on the functional groups and species that comprise plant communities experiencing woodland expansion. Vegetation dynamics of dryland areas with current woody species expansion are increasingly influenced by a variety of global change factors, including elevated CO 2 and temperature, altered fire regimes, and new species invasions (Archer et al. 2017, Tietjen et al. 2017, Miller et al. 2019. Additional research on the mechanisms driving the long-term responses of functional groups and species to these ongoing changes and the effects of management actions across environmental gradients will help to prioritize areas for treatment, select the most effective treatment method, and assure treatment effectiveness.

ACKNOWLEDGMENTS
We thank the SageSTEP team and their agency and university partners for development, implementation, and maintenance of the experimental network. Maggie Gray, Scott Shaff, and David Board provided information on site characteristics and data availability and structure. Jim McIver provided helpful comments on presentation of our results. This is Contribution Number 149 of the Sagebrush Steppe Treatment Evaluation Project (SageSTEP), funded by the U.S. Joint Fire Science Program, the National Interagency Fire Center, and the Bureau of Land Management.