INTRODUCTION

Controlling nonnative invasive grasses in historically fire-maintained grasslands is a major challenge for managers who want to maintain or increase native plant diversity because nonnative invasive grasses can outcompete and replace native grasses and forbs (e.g., Hoffmann et al. 2004; Reed et al. 2005; Gabbard and Fowler 2006; Brooks et al. 2010). Restoring the fire regime may control these invasive grasses, consequently improving ecosystem function, increasing species richness, and decreasing abundance of woody plant cover (Ditomaso et al. 2017; Starns et al. 2020; Nippert et al. 2021). However, if the invasive nonnative grass is fire-tolerant like the native grass species, the effects of fire are less predictable, especially if both species share life history traits (e.g., both are C₄, warm-season grasses) (Bahm et al. 2017; Kral et al. 2018; Palit and DeKeyser 2022). In such cases, a nonnative invasive grass may be favored by prescribed fires. It is not sufficient for prescribed fires to have negative effects on the target species; to be effective management tools for maintaining and increasing native biodiversity, the prescribed fires must differentially affect the target species.

Prescribed fire can be a difficult or expensive tool for many landowners. Other methods of removing aboveground biomass (via grazing, mowing, or clipping) are sometimes used to mimic the effects of fire, with mixed results. Mowing increased species richness in one study (Dee et al. 2016) and had no effect on the herbaceous community in another (Simmons et al. 2007). Targeted grazing can require extra fencing to ensure animals feed on targeted invasive species and comes with the risk that animals might spread seeds (Bailey et al. 2019). Clipping and fire treatments produced only slight differences in total annual biomass yield in Iowa grasslands (Ehrenreich and Aikman 1963), while clipping can also have no effect on the herbaceous community (Petersen and Drewa 2014) or can exacerbate negative effects of fire on a native grass species (Ansley et al. 2006; Ansley and Castellano 2007). It is important to learn more about the effects of these alternatives to burning so land managers who are not able to use fire can use the alternatives effectively to meet their goals.

Bothriochloa ischaemum (King Ranch bluestem; yellow bluestem), an Old World bluestem grass from central Asia, has invaded many grasslands in the southern Great Plains of the United States (Gabbard and Fowler 2006). B. ischaemum is a perennial warm-season bunchgrass that was introduced for domesticated livestock forage and erosion control in the early 20th century (Eck and Sims 1984; Wied et al. 2020). It is the most common nonnative invasive species on the Edwards Plateau, where it commonly forms dense stands that replace most native grass and forb species, dramatically reducing native species richness (Gabbard and Fowler 2006). Within the region, B. ischaemum has few known habitat limitations, except that it is never found under heavy shade of woody plant canopies (Gabbard and Fowler 2006).

B. ischaemum control has been attempted for several decades (Wied et al. 2020). Herbicides are generally ineffective (Robertson et al. 2013; Harmoney et al. 2017). B. ischaemum is tolerant of heavy grazing and mowing (Gabbard and Fowler 2006; Wied et al. 2020). Prescribed fires conducted in the winter (as is typical for this region) have had inconsistent effects on B. ischaemum; most studies report that these winter fires had no effect (Behr 2023) or in some cases increased abundances of this invasive species (Grace et al. 2001; Gabbard and Fowler 2006; Havill et al. 2015; Novak et al. 2021). Prescribed fires occurring in the summer or autumn have been more successful (Simmons et al. 2007; Ruckman et al. 2012; Twidwell et al. 2012; Whiting 2022). It has been suggested that some combination of summer and winter fires can be an effective restoration regime in this region, with summer fires shifting plant communities and winter fires maintaining the effects of the previous summer fire (Taylor 2007).

The differential effects of fire on B. ischaemum and native grasses have been investigated, but more work is needed. In general, it appears that fire management that has a negative effect on B. ischaemum does not have similarly negative effects on native grass species. Wonkka et al. (2017) found minimal native perennial bunchgrass mortality following winter fires but did not measure species responses to fires in other seasons. Other studies examining grass species' responses to summer fires have found that B. ischaemum and the dominant native perennial bunchgrass Schizachyrium scoparium responded differently to summer fires, with B. ischaemum showing greater decreases than S. scoparium, often with little to no damage to S. scoparium (Havill et al. 2015; Whiting 2022).

Many prescribed fire experiments report only the effects of the treatments after 1 or 2 y. Previous studies of summer grassland fires in central Texas have not extended more than 3 y after treatment (Simmons et al. 2007; Reemts et al. 2019). Most of the other long-term studies of fire in grasslands in the region have been done in the quite different plant communities at Konza Prairie in Kansas (e.g., Towne and Kemp 2008) with the notable exceptions of Reemts et al. (2021) and Novak et al. (2021).

We took the opportunity in this study to return 12 y later to an experiment initially burned in 2009, and asked to what extent the treatments still differed and to what extent the initial plant community dominance condition of a plot still affected its composition. We studied both the short-term (2 y) and long-term (12 y) impact of a single summer fire (in 2009) on the invasive grass B. ischaemum and on the native plant community with a replicated and randomized experiment. We also asked whether a plot's initial vegetation (dominated with B. ischaemum or by native species) affected its species composition 2 and 12 y later. Finally, the experiment also included a clipping treatment to simulate mowing, to ask whether clipping can be an effective substitute for fire. If simply removing aboveground biomass is sufficient to tip the competitive balance toward native species, then clipping might achieve many of the same goals as prescribed fire without the obstacles that sometimes make using prescribed fire difficult or impossible, such as county burn bans (TX Loc Govt Code § 352.081 (2021), 2 OK Stat § 2-16-26 (2014)). In addition to measures of plant community composition (richness, cover), we examined possible relationships with soil depth, standing dead biomass, and litter.

Figure 1.

Map of randomized treatment plots within the study site. The experimental fire (2009) burned the whole study area (light-colored polygon) except for control plots (circles), which were not burned. Squares indicate plots that were clipped rather than burned. Maintenance fires (2014 and 2017) that were conducted by the land manager apart from this study burned the whole area including the control plots. Points are not to scale and never overlapped in the field.

METHODS

Study System

This study was conducted at a site in the Simons tract of the Balcones Canyonlands National Wildlife Refuge in Burnet County, Texas, USA (30.8536° N, -97.9744° W). At the beginning of the study, the dominant native species were short- and mid-height grasses, including Schizachyrium scoparium, Nassella leucotricha, Bothriochloa laguroides, and species of Bouteloua and Aristida. The nonnative invasive grass Bothriochloa ischaemum was also abundant at the site.

Experimental Design: In April 2009, we established 80 2 × 1 m randomly located plots (Figure 1). Half were designated as native-dominated (B. ischaemum cover ≤40%) and half were designated as invasive-dominated (B. ischaemum cover ≥60%). Within each of these categories, plots were originally randomly assigned to one of three treatments: burned, clipped, or control, creating a 2 × 3 factorial design.

Figure 2.

A timeline of the study with key event dates (arrows) and Keetch-Byram Drought Index (KBDI) (solid line) from 2009 to 2021 ( https://mrcc.purdue.edu/CLIMATE, weather station ID: USC00411250). The experimental burn took place in July 2009 (widest arrow) when KBDI was >600 (dashed line). KBDI values of 600–800 are associated with severe drought and extreme fire behavior. The winter maintenance burns (medium black arrows) were not conducted as a part of this experiment but did burn the entire study site (including controls) in 2014 and 2017. Both took place when KBDI was <200 (moist soil, no fire danger). Surveys (small black arrows) occurred before the experimental burn (2009) and twice after the experimental burn (2011 and 2021).

To impose the burn treatment, a prescribed fire was conducted in July 2009. Flame lengths were about 1 m in a predominantly backing fire. The Keetch-Byram Drought Index (KBDI) was above 600 (severe drought) on the date of the experimental fire (Figure 2,   Supplemental Table S1 ()). The entire field was burned except for plots assigned to remain unburned (control and clipped treatments), which had mown, wetted strips around them.

The experimental design had to be altered because some of the control plots were burned by accident. The design was split into two: one 2 × 2 factorial design testing burn vs. unburned × initial vegetation (native- or invasive-dominated) and a second 3 × 1 design comparing the burned, clipped, and unburned treatments in only invasive-dominated plots. The native-dominated, clipped treatment was dropped from the design. The two designs shared the invasive-dominated unburned plots and the invasive-dominated burned plots (Figure 3). These two designs had a joint total of 80 plots (counting each plot once even if it was part of both designs).

Clipping was done within a week after the fire in 2009. Plants were clipped to the base and the resulting litter was evenly spread across the plot.

Additional prescribed burns were conducted by the Refuge in January of 2014 and 2017 (  Supplemental Table S1 ()). During these maintenance burns, the entire study area, including all plots, were deliberately burned. The Keetch-Byram Drought Index (KBDI) was very low (<200: no drought) during the two subsequent January burns (Figure 2). No other management or major disturbances occurred at the study site between 2009 and 2021.

Vegetation Surveys: Cover was visually estimated in each plot in the 2009 and 2011 surveys by Andruk. Cover was estimated in each plot by Behr in 2021 using a presence–absence quadrat method (Critchley and Poulton 1998). The percent cover of each species was recorded in spring 2009 (pre-treatment), spring 2011 (2 y post-treatment) and spring 2021 (12 y post-treatment). Survey data from 2011 from 6 of the 80 plots was not available. Only 74 plots (5 treatment–initial vegetation combinations: burned–native N = 26, unburned–native N = 12, burned–invasive N = 18, unburned–invasive N = 10, clipped–invasive N = 8) with data from all three years were used in the analyses that compared the three time points (2009, 2011, and 2021).

Figure 3.

Conceptual diagram showing the 2 × 2 factorial design (surrounded by the double line) and the 1 × 3 clipping treatment design (dashed line). In the 2 × 2 factorial design, the burning treatment and whether the plots were native or invasive dominated at the beginning of the experiment were contrasted with each other. In the 1 × 3 clipping treatment design, the burned–invasive, unburned–invasive, and the clipped–invasive plots were contrasted with each other. The burned–invasive and unburned–invasive plots were used in both designs.

Soil depth, cover of bare ground, litter, and standing dead vegetation were recorded in all 80 plots in the 2021 survey. Soil depth was measured by inserting a metal probe into the ground until it would not go deeper. This probe length in the ground was recorded as the soil depth at each plot (Fuhlendorf and Smeins 1998). Soil depths at these plots ranged from 1 to 30 cm. Litter was defined as dead vegetation lying flat on the ground that was not attached to any plant (Bai et al. 2001). Standing dead vegetation was defined as dead vegetation that was still attached at the base and was standing upright rather than lying flat on the ground (Bai et al. 2001). This distinction between litter and standing dead vegetation was used because the vertical structure of the dead biomass can potentially affect light penetration, soil moisture, fire behavior, and other environmental factors (Knapp 1985; Bai et al. 2001).

Statistical Analyses: Univariate: These analyses each used data from the 74 plots for which we had data from all three surveys. Six response variables were analyzed in separate analyses. The first response variable was species richness, i.e., the number of species present in a plot in a given year. The second response variable was Bothriochloa ischaemum cover. The other four response variables were richness and cover of native forbs and native grasses.

We constructed generalized linear mixed models with either negative binomial or Poisson distribution (these distributions are used for discrete count data) (Bolker 2008; Bolker et al. 2009). Our initial models each included four predictor variables. Because, as described above, the experiment did not include all six combinations of 2009 pre-treatment vegetation (native or B. ischaemum dominance) and treatment (burn, clip, control), we constructed models with a combination treatment–vegetation term of five levels (the five combinations of treatment and 2009 vegetation) as a predictor variable. Survey year (2011 or 2021) was a second predictor variable, and the interaction of survey year and treatment–vegetation was a third predictor variable. Each of the initial models included the 2009 pre-treatment value of the response variable as a covariate because pre-treatment values tend to explain much of post-treatment variation (Gelman and Hill 2007). To account for repeated measures of the same plots, the residual variation was modeled as a repeated measure, with plot nested within treatment as its subject and an R matrix with compound symmetry (Littell et al. 1996; Bolker 2015). AICc (Akaike information criterion, corrected for small sample sizes) values were used to determine whether to retain the interaction between post-treatment survey year and the treatment–initial vegetation combination as a term in the final model of each response variable (Arnold 2010).

After determining that the effect of the treatment–initial vegetation combination term was significant for the given response variable, we did a series of contrast tests to increase the precision of our results (Day and Quinn 1989). If the slopes of the response variable versus the covariate differed significantly among the treatment–initial vegetation combinations, the contrasts compared slopes (Gelman and Hill 2007). If the slopes of the response variable versus the covariate did not differ significantly among treatment–initial vegetation combinations, the interaction term was dropped from the model, and the contrasts compared least-squares means. Because of the missing native-dominated clipped treatment combination, we chose contrasts based on conceptualizing the experiment as having two different experimental designs that had some overlapping terms (Figure 3). First, we used contrasts to test the main effects and interaction term of a 2 × 2 design with burn/unburn and initial vegetation as its two factors (  Supplemental Table S2 ()). Second, we used additional contrasts to compare all three treatments in a one-way design that included invasive-dominated plots only. Each contrast was a function of least squares mean values from the best statistical model for that response variable and was tested with a Scheffé test.

In separate analyses using only the community data collected in 2021, the same six response variables (richness, B. ischaemum cover, cover and richness of native grasses and native forbs) were analyzed with models with the following predictor variables: pre-treatment (2009) value of the response variable as a covariate, soil depth, cover of bare ground, cover of litter, and cover of standing dead biomass. Generalized linear mixed models were constructed with either negative binomial or Poisson distribution, as appropriate (Bolker et al. 2009). AICc values were used to identify the best model for each response variable (Arnold 2010).

SAS 9.4 (SAS Institute 2023) PROC GLIMMIX was used for the univariate analyses.

Statistical Analysis: Multivariate: A separate multivariate analysis was done using nonmetric multidimensional scaling (NMDS) to relate separation in plot characteristics and treatments to differences in plant community composition using the vegan package in R 4.2.1 (Oksanen et al. 2023). We evaluated whether composition differed among treatments with PERMANOVA using the adonis function in the vegan package. To assess between-treatment differences, we calculated Bonferroni-corrected P values using the pairwise.adonis function in the pairwiseAdonis package (Martinez Arbizu 2020). Soil depth was used as an environmental variable in this analysis for all three years because that variable presumably had similar values in 2009 and 2011. Bare ground, standing dead, and litter cover were only used with the 2021 vegetation data as those likely varied from year to year.

RESULTS

As expected, community structure in 2011 resembled that of the pre-treatment community (2009) in many ways. Burned invasive-dominated plots had a weaker relationship between 2009 and 2011 in species richness and richness of native forb species (Figures 4a and 4c). The burned invasive-dominated plots had roughly the same number of native forb species and total native species in both years.

Native grass cover in 2011 was significantly greater in the burned invasive-dominated plots than in any others, including the corresponding control plots (unmanipulated invasive-dominated plots) and even the unmanipulated native-dominated plots (Figure 5a). Native forb cover was similar in all plots except the unmanipulated invasive-dominated plots, where it was significantly lower (Figure 5c).

Figure 4.

Species richness, native grass and forb richness responses to treatment were stronger in 2011 than in 2021. The best model results for each richness response variable are shown here. Axis labels are in log scale. Each line has its own slope in these graphs. Results from slope contrast tests for these variables are in   Supplemental Tables S6, S7, and S8 ().

Native grass and forb richness in 2011 in the clipped plots had little or no relationship to 2009 (pre-treatment) grass or forb richness in those plots (i.e., flat lines in Figures 4c and 4e). Native forb cover was significantly higher in clipped plots than in the corresponding control plots in 2011 (Figure 5c). Native grass cover was also higher in the clipped plots than in the corresponding control plots in 2011, but not significantly so (Figure 5a).

Twelve years after the fire, in 2021, the relationships were weaker than they had been in 2011 but some effects could still be detected. By 2021, the relationships between richness or cover with their pre-treatment values were much weaker and mostly not significant (i.e., the dashed lines of Figures 4 and 5 are close to horizontal; this created significant contrasts between the 2009 and 2011 slopes [  Supplemental Tables S6–S9 ()]). However, the relationship between native grass richness in 2009 and 2021 was significantly different between the burned native-dominated plots and the unburned native-dominated plots; the latter had a negative slope (Figure 4d).

Native grass cover remained significantly greater in the burned invasive-dominated plots than in the unburned invasive-dominated plots (Figure 5b), although there was no detectable difference between burned native-dominated and unburned native-dominated plots in 2021. Native grass cover also remained significantly greater in the unburned native-dominated plots than in the unburned invasive-dominated plots (Figure 5b), indicating a degree of continuity in community composition even in the absence of a summer fire.

Figure 5.

Native grass and forb cover responses to treatment were stronger in 2011 than in 2021. Native grass cover was significantly greater in the burned invasive-dominated plots than the unburned invasive-dominated plots in 2011 and this relationship remained significant in 2021. The best model results for each native species cover response variable are shown here. Axis labels are in log scale. Results contrast tests for these variables are in   Supplemental Tables S10 and S11 ().

B. ischaemum cover was high in the invasive-dominated plots in 2009 and had decreased significantly (from ∼80% to ∼30%) in the burned invasive-dominated plots by 2011 (Figure 6, right side). B. ischaemum cover in the native-dominated plots was low by design in 2009, and remained low in 2011 in these plots, regardless of whether they were burned (Figure 6, left side).

Unexpectedly, by 2021, B. ischaemum cover in the invasive-dominated plots had declined in the unburned plots almost to the level of the burned plots, in which B. ischaemum cover remained low (Figure 6, right side). B. ischaemum cover remained low in 2021 in all the plots that were initially dominated by native species (Figure 6, left side).

Clipping had little or no effect on B. ischaemum cover in 2011, again comparing clipped plots with the corresponding control plots (Figure 6). Burning was significantly more effective than clipping in reducing B. ischaemum cover in 2011 (  Supplemental Table S4 ()). None of the effects of the 2009 clipping treatment were still detectable in 2021.

Multivariate Analysis Results

Prior to treatment, burned–native and unburned–native communities were similar to each other, as were burned–invasive and unburned–invasive communities (Figure 8a). In 2011, the burned–invasive communities were significantly different from the unburned–invasive communities while native-dominated communities remained similar to each other (Figure 8c). In 2021, burned–invasive and unburned–invasive plots returned to their pre-treatment level of similarity (Figure 8e). Communities in the clipped–invasive and the unclipped–invasive plots were never significantly different from each other in any year (Figures 8b, 8d, and 8f). Results of pairwise PERMANOVA tests between the treatments in each year are reported in   Supplemental Table S5 ().

Figure 6.

Cover of the invasive grass B. ischaemum was significantly reduced in plots that were burned and had been dominated by the invasive prior to treatment. This reduction was still detectable 12 y after the summer fire (right side). Plots that were native dominated before the experiment remained native dominated (left side). Symbols indicate treatment: burned (triangles), unburned (circles), and clipped (squares, right side only). Bars are standard errors. The results from contrast tests for this variable are in   Supplemental Table S9 ().

Soil depth was a significant environmental variable in the NMDS analysis in 2009 and 2021, but not in 2011. Soil depth was associated with the native communities more strongly in 2021 than in 2009.

DISCUSSION

A single summer fire in 2009 had strong effects on this central Texas grassland plant community, as did whether a plot was dominated by native species or by the invasive grass B. ischaemum before the fire. The effects of the 2009 fire and initial community dominance on richness were similar to their effects on grass and forb cover, and largely explicable by the negative effects of B. ischaemum on native species. Two years after the 2009 fire, the effects of the fire on the richness of native species were positive. Twelve years after the fire, only the positive effect on native grass richness was still detectable. Additionally, soil depth appears to be a relevant environmental factor affecting the spread of B. ischaemum, with shallower soil appearing to exclude the species.

In this study, a single summer fire successfully reduced the invasive B. ischaemum without harming the native grasses to the same extent and this effect was still detectable 12 y later. Similar reductions in B. ischaemum cover following a summer fire have been reported by Reemts et al. (2021) and Whiting (2022), although Twidwell et al. (2012) found slight increases in the species following a summer fire. The benefit to native grasses in this study supports the findings of other studies that found summer fires did not damage native grass species (Limb et al. 2011; Wonkka et al. 2017), even if they did cause reductions in invasive grass species (Havill et al. 2015; Whiting 2022). The lasting positive impact of the summer fire and resulting reduction of B. ischaemum on native grass species provides encouraging evidence that this treatment can decrease this invasive grass without damaging the native grass species.

The reductions in B. ischaemum cover in the unburned plots over the course of the experiment remain unexplained. A similar unexplained decrease in B. ischaemum cover in unburned plots was also observed in a central Texas grassland by Reemts et al. (2021). They speculated that the species may create such a dense litter layer and reduce available nitrogen in the soil to such low levels that it impedes its own growth (Reemts et al. 2021). It may be that the two cool-season fires of the entire site during the interim, in 2014 and 2017, slowly reduced B. ischaemum cover in the unburned plots as well as maintaining the effects of the 2009 fire, but this is not likely since cool-season fires are not known to cause such reductions in this species (Grace et al. 2001; Gabbard and Fowler 2006; Havill et al. 2015; Novak et al. 2021; Behr 2023). Further study is necessary to determine the mechanisms behind this unexplained decrease in B. ischaemum cover in the unburned plots over time.

B. ischaemum was uncommon in shallower soils on our site, suggesting a previously undocumented habitat limitation for the species. Whereas B. ischaemum (a relatively deep-rooted species; Wied et al. 2020) can probably survive in all but the shallowest soils in this region, it may be outcompeted by natives in shallower soils. The obvious intermingling of B. ischaemumdominated and native-dominated patches in the site precludes explanations based on a shortage of B. ischaemum seeds to colonize native-dominated patches. To our knowledge, this is the first study linking soil depth with the ability of B. ischaemum to invade, become established, and exclude native species, although other studies have found a similar relationship between soil depth and other invasive plant species (Morghan and Rice 2006; Xavier and D'Antonio 2016; Li et al. 2018). Native grasses in this location might be better adapted to shallow soils, as was found in a Saskatchewan study in which the native grass Festuca hallii preferentially rooted in shallower soils whereas the invasive grass Bromus inermis rooted in deeper soil (Li et al. 2018). It is also possible that the native species in our study would also thrive in deeper soils, but they have adaptations that allow them to persist longer than B. ischaemum in shallower, drier soils. In one California study, shallow soil and the dry conditions associated with it were associated with reductions in an invasive thistle's performance and increased performance of a native grass species (Morghan and Rice 2006). Both species performed better in deeper soils, but the native species was better able to withstand the arid conditions associated with shallow soil (Morghan and Rice 2006). Native grasslands with patches of shallow soil may be able to resist invasion by and establishment of B. ischaemum, at least for a while.

The positive relationship between B. ischaemum and standing dead cover reflects this species' ability to produce large amounts of persistent standing dead biomass, unlike most native grass and forb species. Standing dead cover was negatively associated with native grass and forb cover, likely as a result of the stresses imposed by this persistent standing dead B. ischaemum biomass.

The clipping treatment was not effective for reducing B. ischaemum cover, particularly in the long term. Other studies have found similar results for biomass removal treatments like mowing or clipping (Ehrenreich and Aikman 1963; Ansley et al. 2006; Ansley and Castellano 2007; Petersen and Drewa 2014). Dee et al. (2016) found that mowing was effective at increasing species richness, but these increases quickly disappeared once regular mowing treatment ceased. It is possible that repeated clipping might have reduced B. ischaemum over time, but it is doubtful that even repeated clipping would produce the lasting reductions in B. ischaemum observed in this study.

Figure 7.

B. ischaemum cover was positively associated with, and native species richness and cover were negatively associated with, standing dead biomass (a and b) and soil depth (c and d). This suggests a previously undocumented limitation to deeper soils for B. ischaemum.

Summer prescribed fire can have dramatic effects on plant communities in the southern Great Plains. A single summer fire reduced the cover of the invasive, nonnative grass B. ischaemum 2 y later and thereby increased the cover and richness of native species. Some of the effects of this single summer fire, but not of clipping, persisted for 12 y after the initial fire. The results of this study have also increased our understanding of the ecology of B. ischaemum with the finding that the grass does not seem to be able to persist well in shallow soils. The continued existence of native-dominated communities at this site over the 12 y of this study may have been facilitated by patches of shallower soil, which apparently excluded, or at least slowed, the local invasion of B. ischaemum. We conclude that summer fire is likely to be an effective management tool for the control of B. ischaemum in this region, while single clipping or mowing treatments are not. Winter fires may help maintain the effects of a summer fire, but summer fires are needed to cause lasting shifts in community composition.

Figure 8.

Nonmetric multidimensional scaling (NMDS) of plant community compositions in each year and treatment, based on Bray-Curtis dissimilarity (k = 2 for all). Soil depth was used as an environmental variable for all three years since that probably did not change between 2009 and 2021. Soil depth, bare ground, standing dead cover, and litter cover were used as environmental variables in 2021. See   Supplemental Table S5 () for within-year differences between treatment types.