Piñon and juniper tree removal increases available soil water, driving understory response in a sage-steppe ecosystem

Abstract

Over the past century, piñon and juniper trees have encroached into sagebrush steppe lands of the interior United States, and managers have for many years removed trees to stimulate the favored understory. While consistent understory response to tree removal in these semiarid lands suggests that trees outcompete other plants for water, no studies have linked increased soil water to understory response after tree removal. We tested the hypothesis that tree removal at six sagebrush steppe sites increased soil water, leading to increased understory plant cover. Using a structural equation model, we found that before tree removal, trees suppressed shrubs (standardized coefficient [SC] = −0.87), perennial deep-rooted (SC = −0.50) and shallow-rooted bunchgrasses (SC = −0.36), but had no influence on cheatgrass. The model explained between 2% (cheatgrass) and 40% (shrubs) of pretreatment cover variation. Measurement of the same plots six years post-treatment showed that most cover variation was due directly to plant growth, with standardized coefficients between 0.51 (perennial shallow-rooted grasses) and 0.72 (cheatgrass). Competition between cheatgrass and perennial deep-rooted grasses was evident, with perennials having twice the influence on cheatgrass than vice-versa (SC = −0.24 vs. −0.11). Spring soil water (wet-degree days) increased significantly after tree removal, measured as cumulative over 6 years (SC = 0.30), and in the early Spring of year six (SC = 0.16). Treatment-induced increase of cumulative Spring wet degree-days explained variation in shrub cover at year 6 (SC = 0.12) and the increase of early Spring wet degree-days at year 6 led to increases in perennial deep-rooted grasses (SC = 0.24) and cheatgrass (SC = 0.23). We detected no influence of Spring wet degree-days on perennial shallow-rooted grasses. The post-treatment model explained between 34% (shallow-rooted perennial grasses) and 69% (deep-rooted perennial grasses) of variation in understory cover. Most variation was explained by re-measurement of the same populations, followed by treatment effects mediated through increased soil water availability, soil factors, and direct effects of the treatment itself. In conclusion, our model is consistent with the a priori hypothesis that additional wet degree-days due to tree removal is a significant mechanism behind observed increases in understory cover.

INTRODUCTION

Over the past century, woody plant expansion into semiarid grasslands and shrublands worldwide has altered ecosystem structure and function (Archer et al., 2017). The effects of both woody plant expansion and subsequent reduction on understory vegetation are highly dependent on hydrologic conditions that vary with climate (Dammeyer et al., 2016; Williams et al., 2017). Because semiarid lands have short growing periods when soil water is available and temperatures are favorable, expansion of woody plants tends to restrict available resources for the understory herbaceous plant community (Archer et al., 2011, 2017; Archer & Predick, 2014; Williams et al., 2017).

In the interior western United States, native piñon and juniper woodlands (Pinus monophylla; Juniperus occidentalis, and osteosperma) have expanded into sagebrush steppe communities (Artemisia), leading to suppression of understory vegetation (Abdallah et al., 2020; Miller & Tausch, 2001; Ray et al., 2019; Roundy, Miller, et al., 2014), altered hydrology (Ochoa et al., 2018; Pierson et al., 2010), altered fire regimes (Balch et al., 2012; Miller & Heyerdahl, 2008; Whisenant, 1990), and shifted patterns of biodiversity (Knick et al., 2014; McIver & Macke, 2014; Miller et al., 2011). In response to piñon and juniper expansion, land managers have undertaken large-scale tree removal in an effort to achieve a more desirable balance among vegetation functional groups (Miller et al., 2005, 2019). Since trees are likely the top competitors in this semiarid system, most researchers have hypothesized that tree removal causes a release of suppressed understory vegetation due largely to increased availability of soil water (Abdallah et al., 2020; Archer et al., 2011, 2017; Archer & Predick, 2014; Ochoa et al., 2018; Ray et al., 2019; Roundy et al., 2020; Roundy, Young, et al., 2014). Consistent with this hypothesis, several studies have measured parallel increases in both available soil water (Ray et al., 2019; Roundy et al., 2020; Roundy, Young, et al., 2014) and understory plant cover (Chambers et al., 2014; Freund et al., 2020; Miller et al., 2014; Ray et al., 2019; Roundy, Miller, et al., 2014) after piñon and juniper tree removal in sagebrush steppe systems. No studies, however, have definitively linked increased soil water to positive understory plant response after tree removal. Thus, while parallel responses of soil water and the plant understory to tree removal may suggest a relationship between the two variables, the definitive existence of such a relationship can only be demonstrated with a systems analysis, in which the relationship itself is tested.

Using a priori knowledge from previous studies of woodland-encroached sagebrush steppe systems, we developed a multivariate hypothesis that predicts how tree removal affects soil water and understory vegetation, and the relationship between them, given variation in pretreatment tree cover (Figure 1; Table 1). Our hypothesis starts with a pretreatment submodel that predicts how tree cover influences the understory vegetation, including shrubs and grasses (Table 1: Paths 1–5). Previous findings in this system have shown that the influence of trees on understory vegetation will be directly related to functional group rooting depth (Archer et al., 2011; Freund et al., 2020; Miller et al., 2000; Roundy, Miller, et al., 2014; Roundy, Young, et al., 2014), mediated somewhat by soil chemical factors (Miller et al., 2019; Rustad et al., 2020), and this relationship is indicated in our model by heavier dashed lines for deeper rooted functional groups, such as shrubs and perennial tall grasses. After tree removal, the model predicts that the cover of each understory functional group will depend in part on pretreatment tree and understory cover, due primarily to remeasurement of the same plots (Paths 6–9), and in part by the release of soil water as a consequence of tree removal (Paths 10–18) (Bates et al., 2000; Freund et al., 2020; Leffler & Ryel, 2012; Mollnau et al., 2014; Roundy et al., 2020; Roundy, Miller, et al., 2014; Roundy, Young, et al., 2014). Soil nutrient factors are predicted to play a modest role in determining understory cover, both pre- and posttreatment (Paths 5, 19) (Rau et al., 2014).

FIGURE 1. Conceptual a priori model of woodland ecosystem describing anticipated principal factors explaining understory cover pretreatment and 6 years following clear-fell treatment at six warm/dry SageSTEP woodland sites. Dotted-line boxes represent variables hypothesized to influence understory cover. Line thickness indicates anticipated strength of relationship; dashed lines indicate negative relationship and solid lines indicate positive relationship. Rationale for each path of the overall hypothesis is described in Table 1.

TABLE 1 Paths and hypothesized mechanisms represented by a priori multivariate conceptual model (Figure 1), and prior studies upon which the model is based.

Path	Hypothesized mechanism
1, 2	Prior to treatment, trees suppress deeper rooted perennial vegetation, including shrubs and deeper rooted native perennial bunchgrasses
3, 4	Prior to treatment, trees have modest negative influence on shallow-rooted perennial grasses (Poa secunda; Sandberg's bluegrass) and annual grasses
5	Prior to treatment soil depth, phosphorus, nitrogen, and magnesium positively influence understory cover
6, 7, 8, 9	After treatment, cover of each functional group will be dependent on the cover of each group prior to treatment
10	Each year after treatment, tree removal adds additional spring soil water (cumulative soil water)
11	After treatment, tree removal provides additional spring soil water during the year at which soil water is measured
12	After treatment, tree removal provides additional early spring soil water during the year at soil water is measured
13, 14	After treatment, cumulative spring soil water boosts shrub and deep-rooted perennial grass cover during the year at which plant cover is measured
15	After treatment, additional spring soil water boosts deep-rooted perennial grass growth during the year at which plant cover is measured
16, 17	After treatment, additional early spring soil water modestly boosts cheatgrass and shallow-rooted perennial grass cover during the year at which plant cover is measured
18	After treatment, cumulative spring soil water modestly boosts shallow-rooted perennial grass cover during the year at which plant cover is measured
19	The influence of soil factors (depth, P, N, and Mg) will have little relative influence on understory cover measured after treatment, compared with the pretreatment model

Note: Prior studies upon which the model is based: Paths 1–4 (Archer et al., 2011; Freund et al., 2020; Miller et al., 2000; Roundy, Miller, et al., 2014; Roundy, Young, et al., 2014); Paths 5 and 19 (Miller et al., 2019; Rustad et al., 2020); Paths 6–9 (Chambers et al., 2020; Freund et al., 2020); and Paths 10–18 (Bates et al., 2000; Leffler & Ryel, 2012; Mollnau et al., 2014; Roundy et al., 2020; Roundy, Miller, et al., 2014; Roundy, Young, et al., 2014).

We used structural equation modeling (SEM) to test the multivariate hypothesis that an increase in available soil water is an important driver behind the frequently observed increase in understory vegetation after piñon and juniper tree removal. SEM is useful in this context because it allows for the specification of relationships among predictors as a multivariate network hypothesis. In our model hypothesis, for example, the effect of tree removal on the understory is mediated through the path represented by changes in available soil water (Figure 1). Importantly, because the presumption of this kind of mediation is based on a priori knowledge from other sagebrush steppe studies of woodland encroachment (Bates et al., 2000; Freund et al., 2020; Leffler & Ryel, 2012; Mollnau et al., 2014; Roundy et al., 2020), the test of this multivariate hypothesis can contribute to the development of a causal understanding of this system (Cox & Wermuth, 1996). The primary objective of the current study, therefore, is to use SEM to test for the existence of a mediating path (i.e., soil water) that largely explains the often-observed recovery of understory vegetation after piñon and juniper tree removal in sagebrush steppe systems. To test this hypothesis, we use experimental data from several sagebrush steppe sites encroached by piñon and juniper, which are part of the SageSTEP study (Sagebrush Steppe Treatment Evaluation Project).

STUDY SITES

The study area comprises six widely separated SageSTEP sites located across a 1500-km gradient from north-central Oregon to south-central Utah (Table 2; Figure 2). We used data for vegetation (Freund et al., 2020), and soil water and temperature (Roundy et al., 2020), collected both pre- and 6 years posttreatment from untreated control plots and from plots in which all piñon and juniper trees had been removed by clear-felling.

TABLE 2 Study site information.

Site location and treatment year by tree species	Site detailsa	Dominant shrub and grass species
Western Juniper
BC: Bridge Creek, OR, 2006	30 cm, Pacific Maritime; 800–900 m; 25% slope, NW-facing; loam soil from basalt, 43 cm depth, 43% sand	Basin big sagebrush (Artemisia tridentata); bluebunch wheatgrass (Pseudoroegneria spicata); Sandberg bluegrass (Poa secunda); Idaho fescue (Festuca idahoensis); cheatgrass (Bromus tectorum)
Pinyon–Utah Juniper
SR: South Ruby, NV, 2009	38 cm, Pacific/Continental; 2100–2200 m; 8%–30% slope; all aspects; loamy soil from carbonate sedimentary rocks, 52 cm depth, 32% sand	Bitterbrush (Purshia tridentata); Wyoming big sagebrush (Artemisia tridentata); rabbitbrush (Chrysothamnus viscidiflorus); bluebunch wheatgrass (Pseudoroegneria spicata); Great Basin Wild Rye (Leymus cinereus); Indian ricegrass (Achnatherum hymenoides); cheatgrass (Bromus tectorum)
MC: Marking Corral, NV, 2006	35 cm, Pacific/Continental; 2300–2400 m; 6%–20% slope; all aspects; loamy-skeletal soils from andesite and quartz alluvium, 51 cm depth, 58% sand	Wyoming big sagebrush (Artemisia tridentata); bluebunch wheatgrass (Pseudoroegneria spicata); Thurber needlegrass (Achnatherum thurberianum); Sandberg bluegrass (Poa secunda); cheatgrass (Bromus tectorum)
Utah Juniper
ON: Onaqui Mt., UT, 2006	35 cm, Continental/Pacific; 1700–2100 m; 2%–30% slope; E-facing; loamy-skeletal soils from carbonate sedimentary rocks, 46 cm depth, 30% sand	Wyoming big sagebrush (Artemisia tridentata); broom Snakeweed (Gutierrezia sarothrae); bluebunch wheatgrass (Pseudoroegneria spicata); Indian ricegrass (Achnatherum hymenoides); Sandberg bluegrass (Poa secunda); cheatgrass (Bromus tectorum)
SC: Scipio, UT, 2007	37 cm, Continental/Pacific; 1700–1800 m; 2%–28% slope; W-facing; loamy-skeletal soils from mixed sedimentary rocks, 46 cm depth, 39% sand	Wyoming big sagebrush (Artemisia tridentata); rabbitbrush (Chrysothamnus viscidiflorus); Stansbury cliffrose (Purshia stansburiana); bluebunch wheatgrass (Pseudoroegneria spicata); Sandberg bluegrass (Poa secunda); cheatgrass (Bromus tectorum)
GB: Greenville Bench, UT, 2007	34 cm, Continental/Pacific; 1750–1850 m; 2%–28% slope; N-facing; gravelly to cobbly loam soils from volcanic alluvium and colluvium, 44 cm depth, 36% sand	Wyoming big sagebrush (Artemisia tridentata); rabbitbrush (Chrysothamnus viscidiflorus); squirreltail (Elymus elymoides); needle and thread (Hesperostipa comata); cheatgrass (Bromus tectorum)

a Details include annual precipitation, elevation, slope, aspect, and soil type and composition, respectively.

FIGURE 2. Location of six study sites (encircled) in Great Basin and surrounding sagebrush-steppe lands, relative to hydrographic Great Basin, interior western United States (background map made by U.S. Forest Service).

All six study sites have been encroached by piñon or juniper trees in the past 80 years. Five sites are dominated by Wyoming big sagebrush (Artemisia tridentata ssp. wyomingensis), while one site (Bridge Creek) has primarily basin big sagebrush (A. tridentata ssp. tridentata). While the influence of summer monsoons from the Gulf of Mexico tends to increase toward the southeast, the great majority of precipitation falls between October and May at all six sites. The sites vary widely in a number of characteristics relevant to the current analysis. Bridge Creek (BC) is a steep, low-elevation, west-facing western juniper site located in central Oregon. Among the six sites, BC had more than twice the cover of deep-rooted perennial grasses pretreatment. Soils are sandy loams derived from basalt and other volcanic rocks with nearly twice the concentration of magnesium as the other sites. South Ruby (SR) is located approximately 1000 km southeast of BC, and is a high-elevation central Nevada piñon and juniper site with measurement plots arranged across a wide range of slopes and aspects. SR had relatively low grass cover and intermediate shrub cover pretreatment. Soils are loamy, derived from carbonate sedimentary rocks, with relatively high organic carbon and total nitrogen. Marking Corral (MC), located about 80 km southeast of SR, is a high-elevation central Nevada piñon and juniper site, with variable slopes and aspects. Pretreatment functional group cover pattern was very similar to SR. Soils are loamy-skeletal, derived from andesite and quartz alluvium. Onaqui (ON) is a Utah juniper site located approximately 250 km northeast of MC and is an intermediate elevation, variable slope, east-facing site. Prior to treatment, ON was characterized by low cheatgrass cover, relatively high perennial grass cover, and intermediate shrub cover. Soils are loamy-skeletal derived from carbonate sedimentary rocks, with relatively low sand and bicarbonate concentrations. Scipio (SC) is a Utah juniper site located 120 km southeast of ON and is an intermediate elevation, variable slope, west-facing site. SC had very high cover of cheatgrass and shrubs, and relatively low cover of perennial grasses pretreatment. Soils are loamy-skeletal, derived from mixed sedimentary rocks, with low sand and bicarbonate concentrations. Greenville Bench (GR) is located about 150 km southwest of SC and is an intermediate elevation, variable slope, north-facing site. Prior to treatment, GR had the lowest cover of perennial grasses, but intermediate cover of cheatgrass and shrubs. Soils are gravelly to cobbly loams and derived from volcanic alluvium and colluvium.

TREATMENTS AND MEASUREMENTS

Each of the six study sites was a statistical block, with one untreated 10–25 ha control core plot and one similar-sized clear-felled plot, where all trees of any size were cut down and left on site. The clear-fell treatment, intended to remove trees to release the understory plant community, was applied in the fall in 2006 (BC, MC, and ON), 2007 (GR and SC), or 2009 (SR) (Table 2). Because treatments were applied in different years, we expected that some variation in response would be due to interannual weather variation. The clear-fell treatment reduced tree canopy cover to less than 1% in all treated subplots, and while some tree regrowth did occur, average tree cover remained less than 2% by posttreatment year 6. With the exception of the ungrazed BC site, all plots were fenced before the spring of each treatment year to exclude cattle grazing.

Prior to treatment, 15 measurement subplots (0.1 ha) were established within each plot (for a total of 90 vegetation subplots each for the pre- and posttreatment models), which spanned a gradient defined by the relative dominance of trees. Each subplot was assigned to one of three tree dominance “phases,” determined primarily by pretreatment tree cover, as well as local factors such as soils, aspect, and slope: Phase I: dominance of shrubs and herbaceous plants with small trees present; Phase II: codominance of trees with shrubs and herbaceous plants; and Phase III: dominance of trees with lower cover of shrubs and perennial herbaceous vegetation. All plant functional group variables were measured using the line-point intercept method, which generated percent cover estimates for each 0.1-ha subplot (Freund et al., 2020; McIver et al., 2010). The same 180 subplots were measured in this way both pretreatment (year 0) and 6 years posttreatment (year 6). For the present analysis, we used cover estimates for trees, shrubs, perennial deep-rooted grasses, perennial shallow-rooted grasses (Sandberg's bluegrass: Poa secunda), and annual grasses (primarily cheatgrass: Bromus tectorum) (Table 2).

A total of 35 soil water and temperature measurement stations were deployed, one within each of the three tree “phases” within each plot at each of the six sites (the BC control did not have available Phase III subplots, so we were only able to pair vegetation and soil water measurement for 35 stations) (Table 3). We were therefore able to link soil water and temperature measurements from each station to vegetation measurements taken from corresponding subplots having similar pretreatment tree dominance conditions (Appendix S1: Figure S1; Table 3). Each of the 35 soil water and temperature stations was equipped with a Campbell Scientific, Inc. (Logan, UT) CR10X or CR1000 micrologger and multiplexer that deployed 16 soil temperature and soil water matric potential sensors (Roundy, Young, et al., 2014). Thermocouples that measured temperature and gypsum blocks that measured soil water matric potential were buried at 1–3, 13–15, 18–20, and 28–30 cm depths at the eastside dripline of tree and shrub canopies and associated litter mounds and on two interspaces between shrubs and trees (4 depths × 4 microsites = 16 sensors). Volumetric water content was measured in one interspace for each station using a Campbell Scientific CS616 water content reflectometry probe, 30 cm long and slanted from 10 to 20 cm deep. This probe uses time-domain reflectometry (TDR) to measure soil water content. We used a standard calibration curve and thermocouple temperature at the 13–15 cm depth to calculate volumetric water content from TDR output (Campbell Scientific, Inc., 2012–2016). Microloggers were programmed to read sensors every 60 s and to store hourly averages. Thermocouples were made by twisting and soldering copper and constantan wires together on one end (T-type thermocouple wire, Omega Engineering, Norwalk, CT). Resistance data from gypsum blocks (model GB-1, Delmhorst Instrument Company, Towaco, NJ) were converted to water potential using standard calibration curves (Campbell Scientific, 1983). Although several derived variables were calculated for all four seasons at each microsite and at each depth, we focused on spring wet degree-days (sum of hourly temperatures >0°C for each hour of soil matric potential >−1.5 MPa, divided by 24 h/day), and averaged measurements taken across depths and microsites to calculate each variate used in the analysis. We considered that wet degree-days calculated this way was our best representation of favorable soil temperatures and soil water availability for plant growth because degree-days above freezing (>0°C) are commonly used to predict plant phenology and growth, while matric potentials >−1.5 MPa are necessary for nutrients in soil water to flow through soil pores to roots (Leffler & Ryel, 2012). We omitted the 1–3 cm soil depth for wet degree-day calculations to best represent rooting depth for the upper profile-rooting zone where nutrients are most available and where soil water fluctuations are less than in the surface soil (Cline et al., 2018; Leffler & Ryel, 2012). Finally, to aid interpretation, we include precipitation measurements for posttreatment year 6, under the assumption that herbaceous cover in particular would be sensitive to precipitation of the previous wet season. Precipitation was measured at 1–1.5 m height on one station at each site (untreated Phase III) with an electronic tipping bucket rain gage (model TE 525 WS-L, Texas Electronics, Dallas, TX) and removable precipitation adapter for snowfall (model CS-705, Campbell Scientific).

TABLE 3 Distribution of soil water/temperature stations as per pretreatment tree phase in untreated control and clear-fell plots at six study sites, and number of vegetation measurement subplots per station plus pretreatment tree cover range (in percentage) linked to each soil water/temperature station.

Site	Tree species	Treatment	Phase	No. subplots	Pretreatment tree cover range (%)
Bridge Creek	Western Juniper	Control	I	12	1.6–19.1
			II	3	13.7–20.8
		Clear-fell	I	9	2.5–10.9
			II	4	16.9–25.1
			III	2	22.8–29.1
Marking Corral	Pinyon–Juniper	Control	I	3	2.1–10
			II	5	18.3–38.8
			III	6	45.8–62.2
		Clear-fell	I	4	4.7–11.6
			II	6	13.6–27.6
			III	7	24.3–54.5
South Ruby	Pinyon–Juniper	Control	I	4	0.1–12.3
			II	4	10.4–30.1
			III	7	37.2–75.7
		Clear-fell	I	2	5.9–9.6
			II	4	13.7–24.6
			III	9	21.4–68.6
Greenville Bench	Utah Juniper	Control	I	4	5.8–12.5
			II	1	16.1
			III	10	22.6–54.1
		Clear-fell	I	4	4.2–10.5
			II	3	19.4–26.2
			III	6	24.1–44.0
Onaqui	Utah Juniper	Control	I	6	3.2–11.5
			II	8	10.0–18.4
			III	4	23.9–30.8
		Clear-fell	I	3	1.6–11.8
			II	10	11.7–26.4
			III	3	22.2–31.4
Scipio	Utah Juniper	Control	I	4	5.2–10.9
			II	5	13.6–26.9
			III	5	37.3–47.0
		Clear-fell	I	2	8.0–8.8
			II	5	15.7–25.7
			III	6	27.8–42.5

Soil physical and chemical characteristics were analyzed from samples taken before treatment from inter-canopy areas on the northeast corner of one vegetation subplot per phase per treatment plot. Soils were sampled between 1 June and 31 September 2006. Soil cores were 7.62 cm in diameter and taken in 15-cm increments to a depth of 90 cm or until an impenetrable obstruction was encountered. We used mean values for the entire depth profile for the current analysis. Standard procedures for determining soil physical and chemical variables were used, as described in Rau et al. (2014). We analyzed a subset of soil variables that were thought to potentially influence the relationship between soil water and vegetation, in the context of treatment (Rau et al., 2014). Variables analyzed included soil depth, % sand, % silt, % clay, pH, ppm Bray phosphorus, ppm Mg, ppm HCO₃, and % nitrogen.

STATISTICAL ANALYSIS

We used SEM (Grace, 2006) to test the hypothesis that tree removal by clear-felling increased available soil water over 6 years of time, which in turn increased understory cover (Figure 1). We used data from the 180 vegetation subplots and 35 soil water stations to populate our model. The pretreatment model evaluated the cover of four plant functional groups (shrubs, perennial deep-rooted grasses, perennial shallow-rooted grasses, and cheatgrass) as a function of pretreatment tree cover, and within the context of key soil factors that were expected to influence the relationships. These included soil depth and percent sand, considered to be potential indicators of water-holding capacity, and total N, Bray, and bicarbonate-P, and Mg⁺, potential indicators of nutrient availability. Distributional properties of all variables were examined, along with the linearity of relationships, and standard transformations (log and square-root) were applied as needed. Because the soil variables were only sampled pretreatment, they served as a control to make sure posttreatment differences in vegetation cover were not mainly the result of pretreatment soil differences among sites.

We used the piecewise SEM R package (Lefcheck et al., 2016) and a weight-of-evidence model comparison approach (Grace, 2020) to evaluate six alternative models and chose a final pretreatment model for interpretation (Appendix S1: Table S1). Variables in the pretreatment model included tree cover, understory cover (shrubs, deep- and shallow-rooted perennial grasses, and cheatgrass), and soil factors (Figure 1; Appendix S1: Tables S1 and S2). The posttreatment model used the same variables as for pretreatment, but added the measurements of understory cover taken in year 6, as well as the soil water variables expected to be critical in explaining the variation of year 6 understory cover. Posttreatment variables in the analysis included: (1) early (March and April), late (May and June), and total (March–June) spring wet degree-days, both cumulative years 2 through 6 (insufficient data from year 1) and year 6 only; (2) maximum early, late, and total spring volumetric water content at 10–20 cm soil depth, year 6 when it would be most related to year-6 plant growth; and (3) early, late, and total spring wet degree-days, averaged over each of the six posttreatment years. We used a weight-of-evidence approach, along with corrected Akaike information criterion (AIC_c) to select the best fitting of five year-6 models (Grace, 2020). Starting with our initial model, this process involved the use of d-separation tests and coefficient p values to reach a model supported by these metrics, while accumulating a set of models that were subject to multi-model inference using information criteria. The result of this two-stage modeling process was a final inferential model that included the percent variation explained by the model for each of the response variables and quantified the strength of direct relationships between variables, using range-standardized coefficients (Grace et al., 2018). Standardized coefficients were preferable for two reasons: (1) The set of variables we analyzed included percentages, concentrations, and counts, and standardization aligned the scales of all variables analyzed to allow for comparison in response; and (2) range standardization allows for easier interpretation of results, for example, a standardized coefficient of 0.50 between two variables would indicate that the dependent variable would be expected to change 50% across the full range of the independent variable.

RESULTS

The final inferential structural equation model (Figure 3) confirmed that data were consistent with the overall hypothesis of the a priori conceptual model (Figure 1). Specifically, results indicate that removal of trees by clear-felling was associated with an increase in wet degree-days, which in turn led to increases in understory cover. The effects of tree cover in the pretreatment system were as expected, with their strength increasing from no effect for the shallow-rooted cheatgrass to modest effect for the shallow-rooted perennial grass P. secunda, to intermediate effect for deep-rooted perennial grasses, and finally, to a strong influence on the deepest rooted shrubs (Figure 3). Prior to treatment, soil depth was negatively associated with both shallow (standardized coefficient [SC] = −0.13) and deep-rooted perennial grasses (SC = −0.08), phosphorus was negatively associated with cheatgrass cover (SC = −0.13), and the composite variable (Bray-P and bicarbonate-P) was significantly associated with shrub cover (SC = 0.21) (Appendix S1: Table S2; Figure 3). Six years after treatment, remeasurement showed that cover values of all four functional groups were strongly determined by their respective pretreatment covers (SC between 0.51 and 0.72). The clear-fell treatment also had a significant effect on cumulative total spring wet degree-days (SC = 0.29) and on early spring wet degree-days (SC = 0.16) measured at year 6 (Figure 3). In turn, cumulative spring wet degree-days led to increased shrub cover at year 6 (SC = 0.12), while early spring wet degree-days led to increases in both deep-rooted perennial grass cover (SC = 0.24) and to increases in cheatgrass cover (SC = 0.23). There were also direct treatment effects—not mediated through changes in soil water—on all four functional groups, positive for shrubs, deep-rooted perennial grasses and cheatgrass, and negative for shallow-rooted perennial grasses (Figure 3). In the posttreatment model, soil concentration of phosphorus was negatively related to both shrubs (SC = −0.10) and perennial shallow-rooted grasses (SC = −0.20), but the pretreatment soil depth effect disappeared. We detected two composite soil variables in the posttreatment model, explaining significant variation in cheatgrass cover (magnesium, depth = 0.27) and in early spring wet degree-days (phosphorus, magnesium = 0.42). Finally, we detected a modest competitive effect of pretreatment deep-rooted perennial grass cover on cheatgrass cover at year 6 (SC = 0.24), accompanied by a weaker, reciprocal effect of cheatgrass on perennial grasses (SC = 0.11) (Figure 3).

FIGURE 3. Final inferential model indicating data fit to hypothesized conceptual model of tree removal effect on soil water and understory plant growth 6 years after clear-fell treatment at six warm/dry SageSTEP woodland sites. Relative magnitudes of standardized path coefficients are indicated by line thickness; dashed lines indicate negative relationship and solid lines indicate positive relationship. Proportion of variation of endogenous variables explained by the model is indicated by r2 values above each variable. Soil boxes represent unique soil variable found to explain significant variation in vegetation variables (see Appendix S1: Table S2). Cheat, cheatgrass; CSWDD6, cumulative spring wet degree-days, years 2–6; DPG, deep-rooted perennial grasses; ESWDD6, early spring wet degree-days, year 6; SPG, shallow-rooted perennial grasses.

DISCUSSION

Results support the hypothesis that additional soil water, made available through the removal of trees, is an important mechanism behind observed increases in understory sagebrush steppe vegetation measured at posttreatment year 6. While separate univariate analyses reported closely parallel responses to tree removal in soil water (Roundy et al., 2020; Roundy, Young, et al., 2014) and in the understory (Freund et al., 2020; Williams et al., 2017), SEM has confirmed a link between the two responses and takes us one step closer to understanding the cause and effect in this system. While our analyses clearly support the overall hypothesis linking soil water to vegetation, some details of the confirmed model require interpretation. The following discussion will focus on five such details of the final model: (1) The link between cumulative spring wet degree-days and posttreatment shrub cover and the lack of such a link leading to deep-rooted perennial grasses; (2) the link between early spring wet degree-days at year 6 and posttreatment cheatgrass and deep-rooted perennial grass cover, and the lack of such a link leading to shallow-rooted perennial grasses; (3) the relative strength of coefficients in the model, particularly from treatment to wet degree-days to plant cover, compared with the direct link between treatment and plant cover; (4) the influence of soil factors on functional groups in both the pre- and posttreatment models; and (5) the presence of unanticipated relationships detected by modeling, such as competition. Interpretation of these details should bolster understanding of how additional soil water, released by tree removal, contributes to the growth of the sagebrush steppe understory.

Cumulative spring water and shrubs, perennial grasses

The significant increase in shrub cover 6 years after tree removal was due to growth of existing individuals, and to recruitment of new seedlings (Chambers et al., 2020). In both cases, it makes sense that cumulative spring soil water, released by tree removal, explains some of the shrub cover increase measured at year 6. Cumulative spring wet degree-days is an estimate of the degree-days when soil water is available for plant growth each year after treatment (Roundy et al., 2020). Because shrub individuals grow during each spring when there is available water, and because annual growth will add cumulatively to the woody structure of each plant, we would expect that cumulative spring soil water would be linked to shrub cover (Dahl, 1995; Kramer & Kozlowski, 1979; Roundy & Ruyle, 1989; Turley et al., 2003). That shrub cover is not linked to soil water measured only at year 6 can be explained by assuming that any growth in shrubs that occurred in the spring of year 6 would represent only a fraction of the total cover increase we measured, which is a consequence of growth through six successive springs following treatment. On the other hand, deep-rooted perennial grass cover at year 6 was linked to wet degree-days on year 6 but not to cumulative wet degree-days. Perennial grasses would be expected to add or lose tillers each spring and summer depending on water availability (Dahl, 1995; Roundy et al., 1979; Tomlinson & O'Connor, 2004; West, 1988). With additional soil water availability from tree removal, individual plants would produce more tillers and leaf growth than the same individuals had pretreatment. Perennial grass cover after tree removal would be expected to respond to additional soil water availability each year, but especially during the year of measurement.

Early spring water and grasses

The link between year-6 early spring wet degree-days and cheatgrass makes sense, because cheatgrass is a shallow-rooted, winter or early spring annual and would be expected to respond to the amount of water available in the early spring within each year of growth. The lack of a similar link to shallow-rooted perennial grasses (P. secunda) may be associated with a decline in its density during the 6 years since treatment (Williams et al., 2017). This decrease could be associated with increases in tall grass cover or shrubs, especially in clear-fell plots. However, competitive links of this kind were not detected in our model. The highly variable response of P. secunda cover to the clear-fell treatment among sites (Appendix S1: Figure S2) and the lack of a significant link to other variables included in our model suggest that it is responding to other factors. Interestingly, we found a negative association between pretreatment phosphorous and posttreatment cover of P. secunda (Figure 3), supporting Aanderud et al. (2017), who reported a negative association between phosphorous and P. secunda cover after tree mastication. P. secunda variability in response to temperature (Monaco et al., 2005), tree cover, and added resources among sites (Bates et al., 2000; Williams et al., 2017) may be associated with its ecotypic variation (Solomon, 2019), and could explain the lack of response to early spring water in our model.

The link between early rather than late spring wet degree-days and deep-rooted perennial grass cover was associated with more consistent additional wet degree-days from tree reduction across the sites in early spring compared with late spring for year 6 (Appendix S1: Figure S3). For the Utah sites, winter and early spring precipitation was relatively high, but late spring precipitation was very limited for year 6 (Appendix S1: Figure S4). This resulted in very limited additional wet degree-days from tree reduction in late spring, probably because soils were dry for both control and clear-fell treatments for most of late spring (Appendix S1: Figure S3). The MC site had a similar precipitation pattern as the Utah sites, but had greater additional wet degree-days in late spring from tree cutting than the Utah sites. MC is a higher elevation and cooler site than the Utah sites, so early spring soil moisture may have carried over from early to late spring on this site more than on the Utah sites. This may have resulted in more additional wet degree-days due to treatment in late than early spring on that site. Both BC and SR had moderate precipitation in winter and early spring. Late spring precipitation was high for BC and was especially high for SR (Appendix S1: Figure S4). This precipitation pattern resulted in high additional wet degree-days for the clear-fell treatment for BC for both early and late spring (Appendix S1: Figure S3). For SR, additional wet degree-days were higher in early than late spring, probably because both control and clear-fell treatment soils were wet from the unusually high precipitation in late spring. It is intuitive that additional soil water from tree reduction for a season is associated with seasonal precipitation patterns for a site and year. However, the timing of precipitation during cooler or warmer seasons may greatly affect soil water recharge or output through evapotranspiration (Williams et al., 2017). Roundy et al. (2020) found that across all SageSTEP woodland sites additional spring soil water from tree reduction was greater on drier than wetter years. For year 6, deep-rooted perennial grass cover was higher on the clear-fell than control treatment. This associated better with wet degree-days in early spring because it was more consistently higher on the clear-fell than control treatment in early rather than in late spring. This interpretation illustrates the challenges in identifying cause and effect relationships in studies that are entirely field based. In this analysis, the challenge was particularly acute because we evaluated only a single posttreatment year, for sites that were widely distributed geographically, and had been treated in three successive years (i.e., 2006, 2007, and 2009).

Strength of coefficients

In the posttreatment model, the highest coefficients leading to measured cover at year 6 are those that originate from the same functional group variables measured pretreatment. This pattern has been observed in several previous analyses (Chambers et al., 2014; Davies & Johnson, 2011; Davies & Sheley, 2007; Miller et al., 2014; Roundy, Miller, et al., 2014) and is largely a consequence of propagule pressure in the system. Thus, one of the best predictors of posttreatment functional group cover is pretreatment cover, regardless of whether or not a plot has been treated. Yet our model indicates that tree removal also strongly influences posttreatment cover, particularly for shrubs, deep-rooted perennial grasses, and cheatgrass. For these three functional groups, tree removal influenced year 6 cover both directly, and indirectly through soil water. While the indirect effects through soil water were expected and are consistent with our conceptual understanding of this system (Roundy et al., 2020), the significant direct effects require explanation. We believe that the significant effects of tree removal, unexplained by soil water changes, is due primarily to the difference in the granularity of data between soil water and vegetation measurements. Recall that vegetation cover estimates, derived from line-point data, were taken both pre- and posttreatment at 180 0.1-ha subplots, 30 for each of the six sites. Soil water estimates, on the other hand, were taken from a total of 35 stations, only six for each of the six sites (five for BC). Thus, each soil water station was coupled to an average of five vegetation measurement subplots, classified by pretreatment tree phase (Table 3; Appendix S1: Figure S1). Yet, functional group cover estimates from subplots within each phase were highly variable, and it is this variation that most likely contributed to the significant direct effects from treatment to vegetation cover, measured at year 6.

Soil factors

Soil-vegetation paths in the final model are mostly related to among-site patterns of variation, and the granularity of soil core data, compared with vegetation data (soil sampled in only one third of the vegetation subplots). Thus, although paths may indicate cause and effect, they are more likely due to other variables associated with vegetation cover that we did not measure or include in the model. For example, soil paths 2 and 3 show a weak negative relationship between soil depth and pretreatment deep and shallow-rooted perennial bunchgrass cover (Figure 3). Yet for individual sites, there were both positive and negative associations of shallow and deep-rooted perennial grass cover with soil depth. Another example of this is soil path 7, which indicates a relatively strong relationship between the composite variable comprised of magnesium and soil depth, and the cover of cheatgrass at year 6 (Figure 3). Inspection of mean soil features measured at the site level indicates, however, that sites with higher cheatgrass cover (SC and GB) have comparatively shallower depth and lower concentrations of magnesium. Other examples of soil effects that can be explained by among-site patterns of variation include soil paths 1, 4, 5, and 8. All of these significant paths identified in our model could have been due to other soil factors (i.e., Ca, K, soil texture, etc.) that were not included in the model, which varied significantly among sites, and interacted with the factors we selected for analysis. As well, because soil depth, texture, and nutrients were only sampled once for each plot/phase combination (at the same general locations as for soil water), it is also possible that the specific locations we selected for sampling soils inadvertently biased the associations we could detect in our analyses. Thus, while it is possible that our results may point to valid relationships within the system—for example, that root distribution may be more important than depth for sagebrush steppe grasses (Nippert & Holdo, 2015) or that depth by nutrient interactions may play a role in plant growth (Hooker et al., 2008; Leffler & Ryel, 2012; Ryel et al., 2008; Schenk & Jackson, 2002)—we believe it is more likely that our results reflect methodological challenges in testing a complex hypothesis with field data. Finally, recall that each site in our analysis served as a treatment block, and that we had only six warm/dry woodland sites available for analysis. Thus, while the sites served as statistical “replicates,” in reality, all analyses so far have demonstrated that each SageSTEP site is unique (Freund et al., 2020; Roundy et al., 2018; Roundy, Miller, et al., 2014; Roundy, Young, et al., 2014), comprised of an idiosyncratic set of features, many of which contribute in unknown ways to the response variables we measure. We therefore believe that, for the most part, identified linkages between soil factors and vegetation in our model are at best a strong reflection of the challenges in testing system hypotheses in field studies.

Unanticipated relationships

Some additional pathways that we detected require interpretation. We did not expect shrubs to have a significant positive influence on soil water. Yet the cover of shrubs, both pre- and posttreatment, explained modest but significant variation in early spring wet degree-days at year 6. A possible explanation for this is that the two-layered root system of big sagebrush is competitive for soil water (Cook & Lewis, 1963; Sturges, 1973), but trees and shrubs may intercept snow and incident solar radiation, increasing soil temperatures by insolation and increasing soil water by earlier and prolonged snowmelt (Kormos et al., 2017; Myers-Smith & Hik, 2013; Sturges, 1977).

Finally, our model clearly indicates that pretreatment deep-rooted perennial grass cover negatively influenced cheatgrass cover at year 6, while there was also a more modest effect of cheatgrass pretreatment on perennial posttreatment. Many other researchers have found these negative relationships (Blank & Morgan, 2012; Bybee et al., 2016; Chambers et al., 2007, 2014; Roundy et al., 2018), and thus maintaining perennial grass cover has become a key management strategy for supporting resistance to cheatgrass dominance in the sagebrush system (Chambers et al., 2017; Pyke et al., 2017). While these initially omitted linkages make sense based on our understanding of this system, their detection illustrates the value of SEM for uncovering patterns in the data that can then be subjected to further study (Grace, 2006).

CONCLUSION

The identification of mechanisms is a key enterprise in science, because it leads to a better understanding of cause and effect (Grace & Irvine, 2020). In the tree-encroached sagebrush steppe system, prior work had assumed that piñon and juniper trees were the top competitors for the limited resources in this system. Thus, an experiment that removed trees would be expected to produce a competitive release of the understory, and that is exactly what has been found in every univariate analysis thus far conducted (Freund et al., 2020; Ochoa et al., 2018; Ray et al., 2019). It has also been suspected that the primary mechanism behind competitive release in this semiarid system would be the increase in soil water available for growth in the spring and early summer of each year following tree removal, and univariate analysis has also confirmed this (Abdallah et al., 2020; Roundy et al., 2020). Yet the univariate analyses of response to tree removal, measured separately for the understory and for soil water, are inadequate for quantitative evaluation of mechanisms, because they cannot specifically confirm the mechanistic link between the two types of measured responses.

Univariate analyses by themselves cannot therefore definitively demonstrate how this ecosystem functions, cannot fully capture their multivariate complexion, and thus fall short in the key enterprise of understanding mechanisms. This paper has demonstrated that SEM, through its capacity to test a multivariate hypothesis consistent with a prior understanding of relationship, takes us one step closer to identification of mechanisms in this system, which in turn leads to better predictive ability. Simply put, understanding the mechanisms that operate in a multivariate world requires multivariate analyses, applied in such a way that system hypotheses, constructed from insights that stem from prior work, can be tested as a whole, with specific directional pathways that reflect concepts of cause and effect. As Bollen and Pearl (2013) put it, “Researchers do not derive causal relations from a structural equation model. Rather, the model represents and relies upon the causal assumptions of the researcher … the credibility of the structural equation model depends on the credibility of the causal assumptions in each application.” In our application, detecting significant linkages between tree removal, water availability, and understory plants is credible, not because it arises from a data-based modeling exercise, but rather because that modeling exercise rests on a foundation of credible assumptions from prior work. The overall conclusions of this analysis are robust, because they arise from a multivariate model that itself was constructed from a substantial body of knowledge collected over the years about the sagebrush steppe ecosystem. More than anything else, the present model generally confirms expectations implicit in previous analyses and is a testament to the insight of the researchers who have preceded it.

ACKNOWLEDGMENTS

This is Contribution Number 146 of the Sagebrush Steppe Treatment Evaluation Project (SageSTEP), funded by the U.S. Joint Fire Science Program, the Bureau of Land Management (Washington Office), the National Interagency Fire Center (NIFC), and the Great Northern Land Conservation Cooperative (USFWS). Logistical support, for treatment implementation and ongoing site maintenance, has been provided by the Bureau of Land Management and the U.S. Forest Service (Region 5). Supported by the USGS Climate R&D and Ecosystems Programs. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

DATA AVAILABILITY STATEMENT

The dataset used in the analysis for this paper and accompanying metadata and variable definitions (McIver et al., 2022) are available from Zenodo: https://zenodo.org/record/6591345#.YpdsgpPMJp8.