Interannual climate variability mediates changes in carbon and nitrogen pools caused by annual grass invasion in a semiarid shrubland.

Exotic plant invasions alter ecosystem properties and threaten ecosystem functions globally. Interannual climate variability (ICV) influences both plant community composition (PCC) and soil properties, and interactions between ICV and PCC may influence nitrogen (N) and carbon (C) pools. We asked how ICV and non-native annual grass invasion covary to influence soil and plant N and C in a semiarid shrubland undergoing widespread ecosystem transformation due to invasions and altered fire regimes. We sampled four progressive stages of annual grass invasion at 20 sites across a large (25,000 km2 ) landscape for plant community composition, plant tissue N and C, and soil total N and C in 2013 and 2016, which followed 2 years of dry and wet conditions, respectively. Multivariate analyses and ANOVAs showed that in invasion stages where native shrub and perennial grass and forb communities were replaced by annual grass-dominated communities, the ecosystem lost more soil N and C in wet years. Path analysis showed that high water availability led to higher herbaceous cover in all invasion stages. In stages with native shrubs and perennial grasses, higher perennial grass cover was associated with increased soil C and N, while in annual-dominated stages, higher annual grass cover was associated with losses of soil C and N. Also, soil total C and C:N ratios were more homogeneous in annual-dominated invasion stages as indicated by within-site standard deviations. Loss of native shrubs and perennial grasses and forbs coupled with annual grass invasion may lead to long-term declines in soil N and C and hamper restoration efforts. Restoration strategies that use innovative techniques and novel species to address increasing temperatures and ICV and emphasize maintaining plant community structure-shrubs, grasses, and forbs-will allow sagebrush ecosystems to maintain C sequestration, soil fertility, and soil heterogeneity.

INTRODUCTION

Exotic plant invasions threaten biological diversity and ecosystem function (Buckley & Catford, 2016; Lonsdale, 1999) and have long‐term ecological consequences (Mack et al., 2000). Plants modify the soil physical and chemical environment in ways that feed back to vegetation community composition and structure, and non‐native, invasive plant species alter soil properties in ways that differ from native plants (D’Antonio et al., 2009; Sardans et al., 2017; Schimel et al., 1991, 1994). In addition, there is increasing recognition that “year effects” have been generally overlooked in ecological studies (Werner et al., 2020), and that interannual climatic variability (ICV) can lead to boom and bust cycles where plant establishment, plant growth, and soil nutrient cycling are high in wet years and very low in dry years (Hallett et al., 2019; Potts et al., 2019; Puritty et al., 2019). Understanding the differential effects of native and invasive plant species on soil and plant nitrogen (N) and carbon (C) and how these effects interact with interannual climatic variability may help explain the persistence of invasive species in water‐limited environments (Ehrenfeld, 2003; Ehrenfeld et al., 2001).

In the western United States, sagebrush (Artemisia L. [Asteraceae]) ecosystems are being invaded by a non‐native annual grass, Bromus tectorum L. (Poaceae) (Germino et al., 2015). In areas where native perennial grasses and forbs have been displaced, B. tectorum can become the dominant understory component of the A. tridentata community (Chambers et al., 2007, 2016). This dominance can modify the pools and fluxes of soil C and N (Hooker et al., 2008; Leffler et al., 2016), as well as increase fire probability by increasing the abundance and spatial continuity of fine fuels (Bradley et al., 2018; Chambers et al., 2019; Davies & Nafus, 2013). After a fire occurs, plant community composition often experiences several years of instability after which B. tectorum can dominate (Bradford et al., 2014; Davies et al., 2011; Mack, 1981; Mahood & Balch, 2019). Domination after a period of system instability implies a process of niche modification (Fukami, 2015; Grman & Suding, 2010) in which the invader alters the soil's physical and chemical environment in ways that facilitate its persistence (D’Antonio et al., 2009; Sardans et al., 2017; Schimel et al., 1991, 1994). This self‐reinforcing feedback is also likely to be influenced by external drivers like climate (Maestre et al., 2016).

Interannual climatic variability is the primary control on soil water content in dryland systems (Maestre et al., 2016). Declines in the mean or large increases in the variability of water availability may lead to nonlinear top‐down responses in ecosystem structure and function (Belnap et al., 2016; McCluney et al., 2012). These responses include increased invasion potential during periods when soil water content and nutrient availability exceed the uptake capacity of the native community (Davis et al., 2000). Years with increased aridity can result in decreased concentrations of C and N due to reduced inputs from litter decomposition (Crowther et al., 2016; Delgado‐Baquerizo et al., 2013), lower biomass production and soil inputs (Schwalm et al., 2010), and reduced cover of biological soil crusts (Elbert et al., 2012; Ferrenberg et al., 2015; Weber et al., 2016). Annual herbaceous plants are highly responsive to changes in temperature and moisture availability (Bansal & Sheley, 2016; Chambers et al., 2007). This makes them well suited to take advantage of episodic pulses of resource availability and their abundance can have profound effects on soil chemistry.

Soil nutrient cycling can differ between invaded and native ecosystems, in particular cycling of C and N (Nagy et al., 2021; Wilsey et al., 2020). Many invasive species have high growth rates and nutrient concentrations that can result in increased rates of litter decomposition, mineralization, and nitrification (Allison & Vitousek, 2005; Ehrenfeld et al., 2001; Liao et al., 2008; Sardans et al., 2017).

In A. tridentata ecosystems, pools of inorganic and total N are often higher in B. tectorum‐invaded sites compared to uninvaded sites (Booth et al., 2003; Gasch et al., 2013; Hooker et al., 2008; Johnson et al., 2011; J. B. Norton et al., 2004; Stark & Norton, 2015; but see Rau et al., 2011). This is attributed to greater biomass turnover, faster rates of soil N cycling, and redistribution of inorganic N in surface soils (Hooker et al., 2008). But single‐year measurements of high N content are not necessarily indicative of long‐term N storage, because the N cycle is highly sensitive to interannual climate variability (Austin et al., 2004). In years with low, but highly variable soil moisture, the transformation of N by microbes and uptake of N by plants can be asynchronous. This is because microbes respond faster to smaller amounts of moisture than plants (Dijkstra et al., 2012). Furthermore, higher evaporation and soil temperatures in annual grass monocultures (Turnbull et al., 2012) decrease the capacity of plants to activate and uptake labile N following small moisture pulses. Thus, inorganic N can accumulate in dry years due to the production of inorganic N by microbes coupled with a lack of plant uptake (Evans et al., 2001; Hanan et al., 2017; Norton et al., 2008). This unused, labile N may leave the system via leaching or nitrification (Austin et al., 2004; Hanan et al., 2017). When a system composed of deep‐rooted perennial plants is converted to shallow‐rooted herbaceous annuals, the depth that labile N must travel before it can no longer be accessed by the plant community decreases Kulmatiski et al. (2020).

In semiarid shrubland systems of the western United States, invasion by herbaceous species can interact with drought to reduce an ecosystem's capacity to sequester C (Esch et al., 2019). A. tridentata communities are a C sink during the growing season and this is amplified with increased summer precipitation (Germino et al., 2015; Huber et al., 2019). However, they may become a source with increases in dormant season precipitation (Huber et al., 2019). While plants are dormant, the accelerated microbial C mineralization triggered by precipitation generates inorganic C that is susceptible to leaching or atmospheric losses via respiration (Huber et al., 2019), without being offset by the C assimilation resulting from plant growth. A. tridentata typically grows in neutral to slightly alkaline soils (Hooker et al., 2008), and in alkaline soils, some of the inorganic C may be composed of carbonates which can also be lost via leaching, only to accumulate at the bottom of hillslopes (Johnson et al., 2014). B. tectorum monocultures sequester about half of the C of sagebrush communities (Germino et al., 2015), and can become a net source of C to the atmosphere with increasing summer drought (Germino et al., 2015; Prater et al., 2006). As with labile N, any labile C produced by episodic moisture pulses in the summer may be susceptible to leaching or atmospheric losses (Norton et al., 2004).

In sagebrush shrublands, soil beneath shrubs and bunch grasses have higher total C and N than soil in shrub interspaces (Bechtold & Inouye, 2007; Chambers, 2001; Ehrenfeld et al., 2005; Schlesinger & Pilmanis, 1998), and this may be altered with progressive invasion. Environmental heterogeneity creates niche diversity, which facilitates both invasion and coexistence between competitors via niche partitioning (Lundholm, 2009; Melbourne et al., 2007; Shea & Chesson, 2002; Stein et al., 2014; Tilman, 2004). The loss of shrubs from fire and subsequent invasion of annual grasses yields a homogenous vegetation structure. The persistence of this simplified structure may lead to a more homogeneous distribution of soil nutrients. Soil nutrients are lost in the deeper layers and concentrated near the surface (Nagy et al., 2021), where they are being cycled by shallow, annual roots, and over time may become more evenly distributed rather than being concentrated around shrubs. In order to establish in such an environment, the roots of native seedlings must compete directly with annual grasses for water and nutrients in the near‐surface soil (Cheng et al., 2004), which is typically not as fully occupied by deep‐rooted native perennials (Gibbens & Lenz, 2001).

Local variability in the abundance and chemical composition of plants may mediate the effect of climatic forces, but the idiosyncratic ways in which species respond to environmental change are not always linearly related to the effect they have on ecosystem processes (Felton et al., 2018; Suding et al., 2008). Climatic conditions affect plant tissue C:N ratios, which influence flammability, decomposition rates (Grootemaat et al., 2015), and soil inputs of C and N (Norton et al., 2004). In dry years, plants may allocate resources toward root or seed production at the expense of aboveground growth (Franklin et al., 2016; Grime, 1977). This leads to higher quality litter (lower C:N), which allows for faster N cycling in the ensuing years. In wet years, more aboveground growth may lead to higher C:N ratios, lower quality litter, and slower N cycling in the ensuing years. However, the relative abundance of forbs to grasses increases with growing season precipitation (Felton et al., 2018; Hallett et al., 2019). This increases litter quality and decomposition rates because forbs have higher leaf N content (Cornwell et al., 2008). Thus, while increased water availability may increase plant C:N ratios by species, the increased abundance of forbs may cause a decrease in community mean C:N ratios. The resulting lagged effect of climate on litter quality (Evans et al., 2001; Pilliod et al., 2017) feeds back to soil C and N cycles (Yelenik & Levine, 2010).

This broadscale, observational study was designed to evaluate how progressive stages of B. tectorum invasion alter soil and plant nutrients in A. tridentata ssp. wyomingensis ecosystems in the context of interannual climatic variability. In previous studies examining the mechanisms by which B. tectorum influences soil nutrient cycling, soil samples were typically collected from vegetation patches within single sites that were either invaded or not invaded by B. tectorum (Blank, 2008; Evans et al., 2001; Rimer & Evans, 2006; Sperry et al., 2006). Studies that sampled several locations (Norton et al., 2004) included a range of soil types. The only climatic context given was usually 30‐year normals. To our knowledge, no published studies have examined how progressive B. tectorum invasion influences soil and plant C and nutrients at the plant community or landscape scale. To evaluate how invasion and interannual climate variability affect soil and plant C and N, we examined four progressive stages of invasion. We sampled soil, litter, and plant tissue in 2013, which followed 2 years of dry conditions, and in 2016, which followed 2 years of wet conditions. We had three hypotheses: (1) Loss of Artemisia shrubs and increases in B. tectorum cover reduce soil total C and N pools and homogenize soil total C and N and other soil nutrients. (2) Higher water availability results in increased plant cover, greater plant tissue C:N ratios, and decreased soil total N, total C, and C:N ratios. (3) The primary drivers of soil and plant C and N differ for invasion stages with Artemisia shrubs and those dominated by B. tectorum.

METHODS

Study area

The study area encompassed 25,000 km2 Northern and Central Basin and Range Ecoregions (hereafter, Great Basin), in the north‐central part of the state of Nevada in the western United States. Twenty sites with the same soil type were selected for study after reviewing the Bureau of Land Management (BLM) fire and soil maps and consulting with land managers from the BLM Winnemucca Field Office. All sites had similar slopes (0%–5%) and elevations (1297–1607 m) and were located at least 1.5 km apart to ensure statistical independence. Mean temperatures ranged from 19°C in July to −1°C in January. Precipitation ranged from 254 to 304 mm with most arriving as snow in fall and winter. All sites had similar soil textures and minimal, if any carbonate content (Table S1; O’Geen et al., 2017), and represented the same ecological site type: Loamy 10–12 precipitation zone/A. tridentata ssp. wyomingensis (Wyoming big sagebrush) and Achnatherum thurberianum (Piper) Barkworth (Poaceae, Thurber's needlegrass). All sites had the same grazing regime, which was grazing from summer through fall. The grazing allotments where our study sites were located were allotted an average 10 ± 6 billed animal unit months per square kilometer. In other words, one cow with her calf, or five sheep plus their lambs, for 10 months per square kilometer (NRCS, 2003).

The sites had four stages of B. tectorum invasion: (1) six sites with a mature sagebrush canopy, 0%–5% aerial cover of B. tectorum, and no fires within the last 30 years (I. intact sagebrush); (2) five sagebrush sites with 5%–20% B. tectorum cover and no recent fires (II. invaded sagebrush); (3) five cheatgrass‐dominated sites with recent (<5 years before 2013) fires (III. cheatgrass‐dominated); and (4) four cheatgrass die‐off sites, where B. tectorum failed to emerge, establish, or reproduce in areas that supported dense B. tectorum stands in previous years (Baughman & Meyer, 2013) (IV. cheatgrass die‐off). These later sites were a year or two post‐die‐off in 2013, had some B. tectorum plants (~5% cover), and were generally dominated by nonnative annual forbs such as Erodium cicutarium (L.) L’Hér. ex Aiton (Geraniaceae), Sisymbrium altissimum L. (Brassicaceae), and Lepidium perfoliatum L. (Brassicaceae). Five sites (three cheatgrass‐dominated, two intact sagebrush) that were sampled in 2013 were not resampled in 2016 due to access issues.

Field data collection

At each site, we established six, 50‐m transects that were parallel and spaced 20 m apart. Transects were measured in adjacent pairs, yielding three replicates per site. We measured shrub cover and collected composite samples of soils and litter biomass along these transect pairs before active growth began in early‐mid April 2013 and in late March 2016. The transect pairs were resampled during peak production in early‐mid June 2013 and in mid‐late June 2016. During the second visit, we estimated the cover of all plant species and ground cover and collected composite samples of herbaceous vegetation for nutrient analysis. A list of all species encountered is in Table S2.

We measured the cover of all shrub species at each site using the line‐intercept method in 2013 (Elzinga et al., 1998). The length of each shrub that intercepted the transect line was measured, and then, the lengths were summed to calculate total shrub distance. Percent shrub cover was determined by dividing the total shrub distance by the total transect distance.

In April 2013 and March 2016, litter and soil samples were collected systematically every 5 m along each transect for a total of 11 samples per transect (n = 22 per transect pair, n = 66 per site). Litter biomass was removed from a 15 × 15 cm patch at each sampling location along the transect, and then, soil samples were collected from the same locations with a punch auger to a depth of 0–10 cm. Soil samples from each transect pair were placed in a bucket and homogenized, returned to the lab, and sieved to <2 mm. Soil samples (n = 3 subsamples per site) were air‐dried and sent to the Soil, Water and Forage Analytical Laboratory (SWFAL) at Oklahoma State University where they were analyzed for percent total C and total N using an LECO TruSpec CN Analyzer (LECO Corp.), and NO3‐N and NH4‐N using KCl extraction and quantification with a Lachat® autoanalyzer. Soil Mg and Ca were measured using Mehlich 3 extraction and quantification with a Spectro ICP, and SO4‐S using calcium phosphate monobasic extraction and quantification with a Spectro ICP. One bulk density sample was collected systematically at each transect pair, for a total of three samples per plot, with a 6 cm tall × 5 cm diameter metal cylinder. Samples were returned to the lab, dried at 105°C, and weighed after the mass stabilized. Samples were then sieved to <2 mm. The weight and volume of any gravel removed were subtracted from the raw bulk density measurements. We used the mean of the three adjusted bulk density measurements at each plot to convert concentrations measured in the lab to area‐based content. Litter biomass was oven‐dried at 60°C, ground in a Udy mill, and then sent to SWFAL where it was analyzed for total C and N using an LECO TruSpec CN Analyzer (LECO Corp.).

In June, a 0.1 m2 quadrat was placed along the transect adjacent to the litter and soil sampling location and used to visually estimate basal and aerial cover of all species, litter, bare ground, rocks, cryptograms, and dung. Standard protocols were used and frequent calibrations among estimators were conducted (Elzinga et al., 1998). We used these estimates to summarize species cover by functional groups: native shrubs (SH), perennial native grasses (PNG), perennial native forbs (PNF), annual native forbs (ANF), annual invasive grasses (AIG), and annual invasive forbs (AIF). After cover was estimated, aboveground biomass of B. tectorum, Poa secunda J. Presl (Poaceae), and all other herbaceous plants were harvested from the same quadrats. Samples were homogenized for the three different categories for every transect pair (n = 3 subsamples per site) as was done for litter and soils. If either B. tectorum or P. secunda were not collected from within the quadrats, an additional sample was taken nearby from within the site. Homogenized vegetation subsamples were returned to the lab where they were oven‐dried at 60°C, ground, and weighed. Vegetation subsamples were then sent to SWFAL for analysis of total C and N using an LECO TruSpec CN Analyzer (LECO Corp.).

Climate data

We used minimum temperature (T min), actual evapotranspiration (AET), precipitation, and climatic water deficit (CWD) derived from gridMET (Abatzoglou, 2013; Figure 1) to evaluate how climatic variability interacted with plant–soil relationships. We chose these variables because they closely reflect physical conditions important for plant growth and survival. Minimum temperature is important because of the effects of freezing, while AET and CWD represent the supply and demand components of atmospheric water balance (Dobrowski et al., 2013). We calculated precipitation for the water year preceding sampling because of the effects of antecedent precipitation on annual grass biomass and cover (Pilliod et al., 2017).

FIGURE 1

Climatic conditions from January 2011 to January 2017 at the study sites. Dashed vertical lines are the dates of soil and litter sampling, dotted vertical lines are the dates of vegetation cover and biomass sampling. (a–c) The monthly anomaly from 30‐year normals of actual evapotranspiration, climatic water deficit, and minimum temperature, with each line representing one study site. In (d), bars show the monthly precipitation anomaly at the Winnemucca, Nevada airport. The red lines represent yearly averages

Analyzing changes in community composition

Because vegetation biomass and cover can fluctuate in years with different moisture inputs, we examined the observed differences in the understory plant communities from 2013 to 2016 using ordinations. We grouped species cover data in two ways: (1) origin and life form (AIG, AIF, PNG, ANF, and PNG); and (2) life form (annuals, perennials, grasses, and forbs). We then created an ordination with non‐metric multidimensional scaling (NMS) using the origin and life form groupings (Minchin, 1987). We used the plant functional group cover values as arguments for the “envfit” function in the vegan library (Oksanen et al., 2019). This function calculates how each variable relates to the two NMS ordination axes and was analogous to a regression model with the response variable being predicted by the two axes. The strength (R 2) and significance of the relationship was calculated using a permutation test with 9999 replications. To examine how interannual climate, soil, litter and plant nutrients correlated with the ordination axes, we used these variables as arguments to the “envfit” function. To examine how the composition at the different invasion stages changed between the two sample years, we calculated the yearly means of the first two ordination axes for each invasion stage.

Analyzing variation in soil, litter, and plant tissue characteristics between invasion stages and year of sampling

We used linear mixed models in the R package “lmer” (Bates et al., 2015) to determine if soil and plant nutrients differed based on year of sampling and stage of invasion. Invasion stage, year, and the interaction between the two were fixed effects and site (containing three replicated transects) was a random effect. Significance of the fixed effects was determined using a type 3 ANOVA with a Kenward–Roger approximation for degrees of freedom. Mean comparisons were performed using Tukey's honest significant difference tests for multiple comparisons and considered significant at the 95% confidence level (p < .05).

To determine if B. tectorum invasion or interannual climatic variability influenced the plot‐scale heterogeneity of soil resources, we calculated the standard deviation of each soil characteristic for the three transect pairs for each site. We conducted separate Kruskal–Wallis tests on the standard deviations using: (1) all four invasion stages; (2) the two sampling years; and (3) the shrub dominant invasion stages (I and II) versus the herbaceous‐dominated invasion stages (III and IV). We evaluated invasion stages with different levels of shrub dominance based on the hypothesis that loss of shrub cover would correspond with loss of islands of fertility. We used a Bonferroni adjustment when necessary.

Analyzing interactions between functional group cover and interannual climate

To gain an understanding of how interannual climatic variability interacted with plant functional group cover to affect soil nutrients, we built path models using the R package “lavaan” (Rosseel, 2012). We structured our path models such that climatic variables would affect soil nutrients both directly and indirectly through their effects on plant functional group cover. We built two path models with soil total N and C as response variables. One model was for the shrub‐dominated invasion stages (I and II) and the other was for the herbaceous‐dominated stages (III and IV).

We used two groups of exogenous variables and one group of endogenous variables in the path models. Exogenous variables can predict both the endogenous and response variables, but are not predicted by other variables in the model. The first group of exogenous variables accounted for interannual climate variability. We used mean T min and AET, the standard deviation of CWD () for the 180 days preceding sampling, and antecedent precipitation (P ant), calculated as the sum of precipitation for the water year that ended the year prior to sampling. The second group of exogenous variables accounted for site‐level characteristics that are relatively static from year to year. These were biological soil crust (BSC) cover and shrub cover. Endogenous variables predict the response and are predicted by the exogenous variables. These accounted for site‐level characteristics that vary from year to year, likely as a result of climatic variability. These were the cover of the different plant functional groups and litter.

We built our models bottom‐up, beginning with simple, exploratory path models, and then increased complexity. We used a modification index function from the “lavaan” package to check for paths that could be added to the model and improve fit. We removed statistically insignificant variables unless their removal resulted in residual correlations rising above 0.1. We bootstrapped the coefficients of the standardized solution for each model with 5000 random draws and considered the coefficients significant if the 95% confidence intervals did not cross zero. Models were built to maximize (>0.9) two measures of good fit, the Comparative Fit Index (CFI) and Tucker–Lewis Index (TLI), and to minimize (<0.1) two measures of error, the Root Mean Square Error of Approximation (RMSEA) and Standardized Root Mean Square Residual (SRMR).

It is possible that an exogenous climate variable may not have a significant direct effect on the response variable, but may have a direct effect on an endogenous functional group predictor, which then has an indirect effect on the response variable. To examine this phenomenon, we calculated the indirect and total effects of the exogenous climate variables on the response variables for each path model. Indirect effects were the effects of the exogenous predictors on the endogenous plant functional group cover predictors, and then on the response variables (soil total N and C). Total effects were the indirect effects plus the direct effects of the exogenous predictors on the final response variables.

All statistical analyses were done in R (R Core Team, 2020). Data and code are publicly available on Github at https://www.github.com/admahood/invasion_climate_soil.

RESULTS

Changes in community composition

The NMDS ordination that evaluated community composition fit well (stress = 0.108; non‐metric fit = 0.988; linear fit = 0.946), and sites were clustered according to invasion stage (Figure 2a). The NMDS axes represented gradients of perennial to annual cover and forb to grass cover (Figure 2b). NMDS axis 1 was positively correlated with annual plant cover and negatively correlated with perennial plant cover (R 2 = .95, Table S3). Axis 2 was positively correlated with forbs and negatively correlated with grasses (R 2 = 0.87, Table S3). Interannual compositional change was greater for invasion stages III and IV than stages I and II and followed the AIG vector (Figure 2a). This reflected greater herbaceous cover, particularly AIG cover, in the wet year at most sites. Total soil C and N were positively related to NMDS axis 2 (corresponding to forb cover) and negatively correlated with axis 1, which corresponded to perennial cover (Figure 2c). B. tectorum, P. secunda and other plant C, litter N, and P. secunda C:N ratio were positively correlated with axis 1, which corresponded to annual versus perennial cover, and negatively correlated with axis 2, which corresponded to grass versus forb cover (Figure 2c). Litter C and other plant N concentrations were positively correlated with axes 1 and 2, corresponding to annual forb cover (Figure 2c).

FIGURE 2

Non‐metric multidimensional scaling ordination illustrating the differences in understory composition, and soil, plant, and litter N and C between 2013 and 2016. Graphs illustrate (a) the mean change from 2013 to 2016 by invasion stage; (b) significant correlations (p < .05) of plant functional group cover with the ordination; and (c) significant correlations of soil, plant, and litter nutrient concentrations with the ordination. The arrows in (a) represent the direction and magnitude of change. The length of each arrow in (b) and (c) corresponds to the R 2 of the correlation. AIF, annual invasive forb cover; AIG, annual invasive grass cover; PNF, perennial native forb cover; PNG, perennial native grass cover

Soil and litter characteristics

Soil C, N, C:N ratios, and secondary soil nutrients were all lower in 2016 than in 2013 (Table 1). Few meaningful patterns existed among invasion stages for the soil variables (Table 1), but there were several variables with significant interactions between the years and invasion stages. In the intact sagebrush stage (I), total N content was similar in 2013 and 2016; however, in the other stages, total N content was higher in 2013 than 2016 (Figure 3a). A similar but weaker pattern occurred for soil total C (p = .039; Figure 3b). Soil SO4 was lower for all invasion stages in 2016 except in stage IV, where it was similar for both years (Figure 3d). Soil Ca was lower for all invasion stages in 2016 except in stage I, where it was similar for both years (Figure 3e). Litter N was higher in 2016 than in 2013 for stages I, II, and IV but did not differ in stage III (Figure 3g). Litter C and C:N were lower in 2016 than 2013 for stages I and II, but did not differ in stage III, while in stage IV, C:N was lower in 2016 and C was similar in 2 years (Table 1; Figure 3h,i).

TABLE 1

Type III ANOVA results for soil total C and N, soil C:N ratio, soil nutrients, and litter and plant %C, %N, and C to N ratios

Response variable	Year		Invasion stage		Year × Invasion stage
	F	p(>F)	F	p(>F)	F	p(>F)
Soil Total N	32.6	<.001	0.8	.516	6	.001
Soil Total C	105.2	<.001	0.3	.856	3.1	.031
Soil C:N	117.9	<.001	0.6	.594	0.5	.716
Soil NO3	765.7	<.001	1.5	.256	4.5	.005
Soil NH4	104	<.001	1.9	.172	2.3	.085
Soil NH4 + NO3	363.6	<.001	1.8	.182	2.3	.085
Litter N	77.4	<.001	0.8	.519	6.6	<.001
Litter C	7.2	.008	4.5	.015	2.5	.066
Litter C:N	48.3	<.001	2.3	.104	5.4	.002
Bromus C	128.9	<.001	7.4	.002	4.9	.003
Bromus C:N	173.8	<.001	0.7	.563	2.1	.101
Bromus N	163.8	<.001	1.4	.264	3.3	.023
Poa C	3.1	.081	0.6	.628	0.4	.726
Poa C:N	122	<.001	1.6	.241	1.7	.177
Poa N	248.3	<.001	3.7	.037	11.8	<.001
Other C	20.8	<.001	5.6	.006	4.5	.005
Other C:N	60.4	<.001	6.2	.004	2.2	.088
Other N	59.2	<.001	5.9	.004	3	.035

ANOVAs are based on mixed models with year and invasion stage as fixed effects, and site as the random effect. Denominator degrees of freedom (not shown) and F‐statistics were calculated using the Kenward–Roger method.

FIGURE 3

Effects of year and invasion stage on soil total C and N, soil C:N ratio (a–c), litter %C, %N, and C:N ratio (d–f), and soil total mineral N, nitrate, and ammonium (g–i). Open circles represent individual measurements. Diamonds and error bars represent the marginal means and 95% confidence intervals estimated from the linear mixed models. Year was always significant. When the interaction between invasion stage and year is significant, colored lowercase letters above the bars illustrate differences between years within each invasion stage based on the Tukey Honest Significant Differences post hoc test. Gray capital letters illustrate significant differences among invasion stages when the interaction was not significant

Analyses of standard deviations of soil total C and N, litter %C, %N, and C:N ratios showed differences among invasion stage and year. When sites were grouped into the four invasion stages, within‐site standard deviations differed only between stages I and IV for soil total C (Table S4). However, when sites were grouped by amount of A. tridentata cover (stages I and II vs. III and IV), within‐site standard deviations in soil total C and soil C:N ratio were lower for stages III and IV than for stages I and II (Table 2). When the data were grouped by year, standard deviations of soil total N were higher in 2013, the dry year, while variation in litter N was higher in 2016, the wet year (Table 2).

TABLE 2

Comparisons of stages I–II and stages III–IV and of the wet (2016) and dry (2013) years based on standard deviations calculated from three replicate samples from each site, averaged within each category

	Stages I–II	Stages III–IV	2016 (Wet)	2013 (Dry)
Soil total N (kg/ha)	96.25	85.37	59.73b	116.33a
Soil total C (kg/ha)	1275.23a	847.25b	939.78	1184.09
Soil C:N	1.33a	0.73b	0.86	1.20
Litter N (%)	0.12	0.09	0.14a	0.08b
Litter C (%)	2.79	2.54	3.22	2.23
Litter C:N	5.91	6.60	6.68	5.87

Letters indicate significant differences between grouped stages and years based on Kruskal–Wallis tests.

Plant tissue characteristics

Carbon concentrations of P. secunda, B. tectorum, and other plants were higher in the wet year (2016) than the dry year (2013), and nitrogen concentrations were lower in the wet year. Because N concentrations were lower in the wet year, C:N ratios were higher (Figure 4; Table 1). B. tectorum had lower C concentrations in stage I than all other stages in the dry year (2013; Figure 4). P. secunda had higher N concentrations in stage III than II in 2013, but N concentrations did not vary by invasion stage in 2016. Both species had higher C concentration in stage III than stages I and II and higher N concentrations in stage IV than in stages I and II. As a result, C:N ratios were lower in stage IV than in stages I and II. The contrast in nutrient concentrations for other plants between the shrub‐dominated versus herbaceous‐dominated invasion stages is likely due to the difference in species composition. The other plant species present in stages I and II were generally annual and perennial native forbs, while those in stages III and IV were mostly annual invasive forbs (Table S5).

FIGURE 4

Effects of year and invasion stage on plant tissue C and N concentrations. Open circles represent individual measurements. Diamonds and error bars represent the marginal means and 95% confidence intervals estimated from the linear mixed models. Year was always significant. When the interaction between invasion stage and year is significant, colored lowercase letters above the bars illustrate differences between years within each invasion stage based on the Tukey Honest Significant Differences post hoc test. Gray capital letters illustrate significant differences among invasion stages when the interaction was not significant

Path models for soil total nitrogen and carbon

Path models fit well with values approaching 1 for measures of fit (CFI, TLI) and values approaching 0 for measures of error (RMSEA, SRMR, Table S6). Covariance matrices for each path model are in Tables S7 and S8. The models for soil total C and N were quite different when split between the shrub‐dominated invasion stages (I and II) and those dominated by invasive annuals (III and IV). The model for stages I and II had many direct climate effects on soil total C and N, with indirect effects through the litter C:N ratio and perennial native grass cover. In contrast, in the model for the later invasion stages, soil total C and N were strongly influenced by annual invasive grass cover.

The model for invasion stages I and II explained 84% of the variation for soil total C and 61% for soil total N (Figure 5a). Litter C:N, PNG, shrub cover, and were positively related to soil total C. P ant was negatively related to soil total C, but positively influenced PNG and negatively influenced litter C:N. The indirect and total effects of P ant on soil total C mediated through litter C:N were negative, while the effects of P ant mediated through PNG were neutral (Table 3). Shrub cover and PNG had positive effects on soil total N, while P ant and AET had negative effects. The indirect effects of P ant and AET mediated through PNG were both neutral.

FIGURE 5

Path models predicting soil total carbon and nitrogen for (a) invasion stages I and II, and (b) invasion stages III and IV. Paths are based on bootstrapped standardized solutions (R = 5000); those with bootstrapped confidence intervals crossing zero have been omitted. The numbers along each path are the median standardized coefficients, which also correspond to the line widths. The numbers below each endogenous or response variable are the median bootstrapped R 2 values. Gray and black lines indicate positive relationships; orange and red lines indicate negative relationships. , standard deviation of climatic water deficit; AET, actual evapotranspiration; AIF, annual invasive forb cover; AIG, annual invasive grass cover; BSC, biological soil crust cover; Csoil, soil total carbon; Litter, litter cover; Nsoil, soil total nitrogen; NF, native forb cover; P ant, antecedent precipitation; PNG, perennial native grass cover; shrubs, shrub cover; T min, minimum temperature

TABLE 3

Indirect and total effects of climatic variables on soil total nitrogen and carbon

Exogenous variable	Effect type	Mediator(s)	Soil C	Soil N	Invasion stages
Actual evapotranspiration	Indirect (single pathway)	PNG	0.19	0.25	I and II
	Total	PNG	−0.06	−0.26	I and II
Antecedent precipitation	Indirect (single pathway)	LCN	−0.3*		I and II
	Indirect (single pathway)	PNG	0.26	0.34	I and II
	Total	LCN	−0.74*		I and II
	Total	PNG	−0.17	−0.27	I and II
Actual evapotranspiration	Indirect (multiple pathways)	AIF AIG	−0.098	−0.079	III and IV
	Indirect (single pathway)	AIG	−0.31	−0.27	III and IV
	Total	AIF AIG	−0.17	−0.14	III and IV
	Total	AIG	0.081	0.13	III and IV
Antecedent precipitation	Indirect (single pathway)	AIG	−1.2*	−1	III and IV
	Total	AIG	−0.7*	−0.5	III and IV
Minimum temperature	Indirect (single pathway)	AIG	0.55*	0.48	III and IV
	Total	AIG	NA	0.54*	III and IV

Values are the median of the bootstrapped (n = 5000) coefficients, and asterisk indicates that the 95% confidence intervals did not cross zero. Median coefficients for the direct effects are displayed in Figure 5.

The model for invasion stages III and IV explained 62% of the variation for soil total C and 59% of the variation for soil total N (Figure 5b). Both soil total C and soil total N were positively influenced by and negatively influenced by AIG. The indirect and total effects of P ant mediated through AIG were negative for total soil C and neutral for soil total N. The indirect effects of T min mediated through AIG were positive for soil total C and neutral for soil total N (Table 3). The total effect of T min mediated through AIG was positive.

DISCUSSION

Here, we demonstrate that the well‐documented sensitivity of herbaceous cover to inter‐annual climatic variability (Felton et al., 2021; Hallett et al., 2019; Pilliod et al., 2017; Witwicki et al., 2016) mediates the effects of progressive invasion of exotic annual grasses and changes in plant functional group 3 composition on soil total C and N, and plant tissue C and N. Our NMS ordination indicated that differences in community composition were affecting soil total N and total C pools, while our linear mixed models indicated that year effects were a primary driver. Using path models to incorporate the effects of interannual climatic variability, we showed that different variables were determining soil total C and soil total N in the different invasion stages. Furthermore, climate‐driven differences in herbaceous plant cover drove soil total C and soil total N in different directions depending on the composition of the herbaceous vegetation, with annual invaders resulting in losses of soil total C and soil total N, and native perennials resulting in increases (Figure 5).

To our knowledge, this is the first study in this system where multiple sites were sampled across a broad geographic extent, and both the stage of invasion and interannual climatic variability were taken into account. Prior studies on the relationship between B. tectorum invasion and soil total N and C pools are conflicting, with reports of both higher (Booth et al., 2003; Gasch et al., 2013; Hooker et al., 2008) and lower (Rau et al., 2011; Stark & Norton, 2015) levels after invasion. Failure to adequately account for spatial and interannual climate variability within an ecosystem is a common shortcoming in ecological research (Werner et al., 2020) and helps explain the discrepancies in prior studies.

Annual grass invasion and loss of shrubs reduce and homogenize soil total N and C

The path model for invasion stages III and IV indicated direct positive relationships between climatic variables and AIG cover, and direct negative relationships between AIG cover and soil total N and soil total C (Figure 5b). It is likely that the observed interactions between invasion stage and year effects (Table 1) are the result of climate amplifying the negative effect of annual grass productivity on soil total N and soil total C in stages III and IV. Our finding that increased annual grass cover was associated with soil N and soil total C losses in invasion stages III and IV in the wetter year is consistent with an emerging consensus that the rate of nutrient cycling is likely higher in B. tectorum‐invaded environments in the Great Basin, and this is amplified in higher moisture years. One mechanism for faster nutrient cycling under B. tectorum in favorable years may be increased root exudation of N into the soil by B. tectorum (Morris et al., 2016), which initiates a self‐reinforcing feedback that facilitates B. tectorum growth (Blank et al., 2016; O’Conner et al., 2015). Increased nutrient cycling can lead to labile N and C removal from near‐surface soil layers via percolation, respiration, denitrification, and runoff (Austin et al., 2004; Norton et al., 2004). In drier years, water from smaller precipitation events that does not provide sufficient moisture for plants to activate is still used by soil microbes to mineralize N and C (Dijkstra et al., 2012; Evans & Burke, 2013; Fierer et al., 2003; Fierer & Schimel, 2002).

The loss of biological soil crust (BSC) due to soil disturbance and fire often accompanies B. tectorum invasion, and further reduces the capacity of invaded sites to retain soil N and C. At our study sites, BSC cover was lower at the invaded sagebrush stage (II) than in the uninvaded stage, and essentially absent at cheatgrass‐dominated (III) and cheatgrass die‐off (IV) stages (Table S5). In addition to providing N input (Belnap et al., 2016; Condon & Pyke, 2018; Weber et al., 2016), BSC reduces the rate of soil drying, and thus gives plants more time to activate and uptake newly available resources (Austin et al., 2004). The loss of BSC as a source of N input, combined with increased nutrient cycling leading to N and C exiting the system, may add up to significant yearly losses in near‐surface soils (Austin et al., 2004).

The lack of shrub cover and uniform distribution of annual herbaceous plants at our sites in invasion stages III and IV is common post‐fire and leads to the loss of “islands of fertility” (Austin et al., 2004; Bechtold & Inouye, 2007; Doescher et al., 1984; Evans et al., 2001; Germino et al., 2018). We found statistically lower within‐site variability of soil total C and soil C:N ratios in invasion stages that had lost shrub cover (Table 2). Lower resource heterogeneity reduces the availability and diversity of niches that species can occupy, making it harder for native species to reestablish (Melbourne et al., 2007; Tilman, 2004). Thus, a reduction in within‐site variability and homogenization of soil total C, as we observed here, may contribute to the continued long‐term dominance of invasive annual plants.

Soil, plant, and litter C and N concentrations are strongly influenced by interannual climate variability

Interannual climate variability affected every variable we measured: soil total C, soil total N, and the C:N ratios of litter and soil, and the tissues of P. secunda, B. tectorum, and other plants (Table 1). The large interannual differences in soil total C (Figure 3b) were unexpected, because soil organic C is usually considered to be fairly stable (but see Lehmann & Kleber, 2015; Morgan et al., 2016). In sagebrush ecosystems, soil water content is the primary driver of soil C mineralization (Norton et al., 2012). When soil is wet, C cycling rates can be an order of magnitude greater than when soil is dry (Saetre & Stark, 2005). Increased litter inputs associated with B. tectorum invasion are associated with higher labile C in the top 10 cm and higher C mineralization rates (Norton et al., 2012). In addition to losses from higher rates of soil respiration, this increased nutrient cycling may be leading to significant C losses through vertical movement of labile C from the top 0–10 cm sampled here to deeper soil depths (Nagy et al., 2021).

The cover and tissue nutrient content of the herbaceous plants were strongly influenced by interannual climatic variability (Figure 4). The resulting effect on litter chemistry may have regional scale consequences for fire risk and nutrient cycling. Plant biomass produced in dry periods becomes litter with low C:N ratios (Evans & Burke, 2013), which may persist for 2–3 years (Pilliod et al., 2017). This can lead to increased rates of nutrient cycling and, thus, loss of mineralized N through leaching and gaseous loss of CO2, nitrification of NH4, and denitrification of NO3 when soil moisture is available (Austin et al., 2004; Evans & Burke, 2013). In contrast, biomass produced in wetter periods becomes litter with higher C:N ratios. In the wet period, litter C:N was higher, B. tectorum and P. secunda C:N ratios were universally higher, and other plant C:N ratios were higher in invasion stages III and IV. The layer of dead and cured fine fuel created by the pulse of growth in wet years will potentially persist longer than the fuel produced in dry years, because litter with higher C:N ratios may decompose more slowly (Grootemaat et al., 2015; Pausas et al., 2017). Consequently, in wet years, there is not only an increase in fine fuel biomass, but those fine fuels may be more flammable and take longer to decompose.

Drivers of soil C and N differ with invasion stage

Changes in species composition and simplification of plant community structure to a single herbaceous cover layer results in different pathways through which C and N travel. In all stages, the strongest direct predictors of C and N were abundance of the dominant grass type, that is, PNG in stages I and II, and AIG in stages III and IV (Figure 5). In stages I and II, PNG had positive relationships with soil total N and soil total C, and the negative effect of higher water availability on soil total N and C was dampened by its positive influence on PNG cover. In these stages, sites with higher cover of shrubs, BSC, NF, and PNG were all associated with higher soil total N and C. This agrees with prior work indicating that sagebrush communities sequester C (Prater et al., 2006), but that their ability to sequester C declines as the cover of key ecosystem components declines. In stages III and IV, the simplified species composition yielded a simple path model centered on a single pathway through AIG. In these sites, P ant led to higher AIG which led to lower soil total N and soil total C. In contrast to Artemisia‐dominated stages I and II, the negative effect of higher water availability on soil total N and soil total C was amplified by its effect on the abundance of AIG, agreeing with Prater et al. (2006) who found B. tectorum monocultures to be a net C source.

Litter chemistry and abundance are known to change after invasion, and this has important consequences for soil nutrient cycling (Norton et al., 2008). Because the vast majority of litter is dead plant material from the previous growing season, litter C:N ratios are reflective of the conditions of the preceding growing season. Here, litter C:N ratios operated differently between the soil total C and soil total N path models. In stages I and II, litter C:N ratio was positively related to soil total C. Antecedent precipitation drove decreases in both soil total C and litter C:N ratios, and higher litter C:N ratios amplified the effect of P ant on soil total C (Table 3). Huber et al. (2019) also found that increases in winter rainfall can turn soil under sagebrush plants into a net C source. In stages III and IV, litter C:N ratio was not a significant part of the model, perhaps because the litter layer was composed mainly B. tectorum detritus, and the C:N ratio did not differ between years.

Implications for management

Our results indicate that reductions in aboveground carbon storage resulting from the loss of shrubs (Fusco et al., 2019) are accompanied by belowground carbon losses. These losses may increase over time with ongoing climate change and progressive conversion to annual grass dominance. In the Great Basin, temperature has increased by 1°C above historic averages and is projected to increase by 2–6°C by 2100 (Bradford et al., 2020; Melillo et al., 2014). Future projections of precipitation are mixed, but models agree that there will be higher interannual variability (Bradford et al., 2020). These trends are resulting in high climate suitability for B. tectorum and other annual invaders (McMahon et al., 2021) as well as increased risk of large‐scale wildfire across much of the area (Brown et al., 2021; Gao et al., 2021). Sagebrush shrublands are currently experiencing among the largest wildfires in the western United States (Brooks et al., 2015; Dennison et al., 2014), loss of the dominant shrub species (fire‐intolerant A. tridentata), and active conversion to grasslands dominated by invasive annuals. This conversion will likely result in progressive decreases in carbon storage if our ability to restore these ecosystems does not improve (Arkle et al., 2014; Knutson et al., 2014).

Restoring sagebrush ecosystems and maintaining their capacity to store and retain carbon and soil nutrients are becoming more challenging as temperatures and interannual climate variability increase. After wildfire in many sagebrush ecosystems, there is a short temporal window outside of which restoration is not likely to succeed. In the first 1–3 years after fire, B. tectorum typically has relatively low density and cover and competition with seeded species is low (Chambers et al., 2014; Urza et al., 2017, 2019). The soil resource islands beneath shrubs increase post‐fire establishment of seeded perennials (Germino et al., 2018) and can persist for several years after fire (Bechtold & Inouye, 2007). Our results indicate that after conversion to annual grass dominance, these resource islands eventually disappear. The homogenization and depletion of soil nutrient pools may negatively affect establishment processes by reducing both resource availability (Davis et al., 2000) and niche diversity (Fukami, 2015; Melbourne et al., 2007). Strategies for overcoming the post‐fire establishment bottleneck include optimizing the timing of restoration by anticipating cooler and wetter years and taking advantage of them when they occur (Bradford et al., 2018; Hardegree et al., 2018; Shriver et al., 2018). Restoration techniques for addressing interannual climate variability include seed priming and coating to shorten germination times and maximize water availability for seedling establishment (Pedrini et al., 2020). Various site preparation techniques also increase seedling establishment (Farrell et al., 2021; Herriman et al., 2016), and methods used in more arid systems for increasing soil water availability, such as soil imprinting (Doherty et al., 2020), may be beneficial.

Selecting restoration species that are adapted not only to increased temperatures and aridity but also interannual climate variability around these trends may increase restoration success (Bradley St. Clair et al., 2013; Gibson et al., 2016). Current approaches use climate‐informed seed transfer zones to minimize the possibility of introducing maladapted genotypes (Bower et al., 2014; Gibson et al., 2019). It may be necessary to modify these approaches to use current climate to determine suitability but use future climate scenarios to help determine the species and zone of application (Richardson & Chaney, 2018; Shryock et al., 2018).

In many areas, sagebrush restoration, regardless of the technique, is unlikely to succeed because of prolonged annual grass dominance or temperature increases that have already occurred (Germino et al., 2018; Shriver et al., 2018, 2019). Our results indicate that after invasive annuals have displaced perennials, the system is more sensitive to interannual climate variability and more susceptible to soil C and N losses. Thus, it may be prudent for land management in these areas to shift the focus from restoring particular species to restoring plant community structure—deep‐rooted shrubs, grasses, and forbs (Balazs et al., 2020; Bradley St. Clair et al., 2013; Laughlin, 2014). Restoration with native species other than sagebrush may not have the same habitat benefits as successful sagebrush restoration. However, if native species can be established that compete effectively with invasive annuals, tolerate fire, and persist in future climate conditions (Butterfield et al., 2017), the system may be able to maintain C sequestration, soil fertility, and soil heterogeneity in the face of interannual climate variability.

Implications for future research

Future research should investigate the effects of annual grass invasion on ecosystem function in different stages of invasion at larger spatial extents and with enough repeated measurements to capture long‐term trends in the context of interannual variability (Cusser et al., 2021). Functional components include stocks of above‐ and belowground biomass, soil total C and N storage, and soil water storage. Annual or subannual data collected at temporal extents long enough to capture the full range of variability are necessary to disentangle long‐term trends of these functional components from the confounding effects of interannual variability (Werner et al., 2020), time since invasion, time since disturbance, and soil characteristics (Austin et al., 2004; Noy‐Meir, 1973). Some prior studies had well‐defined invasion stages (Rau et al., 2011), or a well‐defined time of introduction (Stark & Norton, 2015), but to our knowledge, no study has adequately captured the range of interannual variability in this system. Time since initial introduction may not always be feasible to ascertain, but annual grass invasions initiated during the satellite record may be detected from seasonal differences in NDVI (Boyte et al., 2019; Bradley & Mustard, 2006). Time since the loss of woody plant cover may be determined from satellite‐derived datasets of fire occurrence (Balch et al., 2020; Hawbaker et al., 2017) and functional group fractional cover (Jones et al., 2018). Understanding the inter‐ and intra‐annual variability in each component of above‐ and belowground C and N pools and fluxes over decadal temporal extents will allow for the determination of rates of change at each site, rather than relying on space for time substitutions.

CONCLUSIONS

Episodic differences in temperature and water availability have predictable effects on B. tectorum invasion and soil nutrients, but have seldom been adequately analyzed. Our results indicate that annual grass invasion changes mechanisms and pathways of nutrient cycling and carbon storage, especially after the loss of shrubs. While soil water content is perhaps the single most important factor in nutrient cycling in dryland ecosystems, our research illuminates some nuance. We found that the strength and direction of the effects of water availability on soil total C and N are dependent on plant functional group composition, because long‐lived, deep‐rooted plants tend to store nutrients belowground, while short‐lived, shallow‐rooted plants tend to extract nutrients. Thus, while changes in the mean and variation of water availability was the primary control, the effects of water availability on soil C and N were amplified by invasive annuals and dampened by native perennials. These changes in mechanisms and pathways of nutrient cycling and carbon storage that result from the conversion of perennial to annual dominance will likely lead to the long‐term depletion of soil resources and a decline in resource heterogeneity.

Supplementary Material

Click here for additional data file.