Restoring the fire–grazing interaction promotes tree–grass coexistence by controlling woody encroachment

Woody encroachment can convert grasslands and savannas to shrublands and woodlands, so understanding the processes which regulate woody encroachment is necessary to conserve or restore these ecosystems. We hypothesized that recreating the fire–grazing interaction would limit woody encroachment because focal grazing increases fuel accumulation on unburned areas and increases browsing on emergent woody plants in burned areas. This study was conducted in the Grand River Grasslands of Iowa and Missouri (USA) on 11 sites (15.4–35.0 ha). Each site was assigned to one treatment: patch-burn-graze (n = 4), with spatially discrete prescribed fires and free access by cattle (the fire–grazing interaction); graze-andburn (n = 4), with free access by cattle and one burn of the entire site every 3 yr; or burn-only (n = 3), with one site-wide burn every 3–5 yr and no grazing. The burn-only treatment increased woody encroachment fourfold compared to the graze-and-burn and patch-burn-graze treatments (130.2, standard error [SE] = 16.0; 20.9, SE = 12.0; and 46.3, SE = 10.8; plants/200 m). The patch-burn-graze treatment had 2– 3 cm more accumulated fuel and woody plants which were 12% shorter, on average, than the other treatments (comparing eight common species). The movement of large herbivores also appeared to decrease the frequency of woody species which spread vegetatively. Our work illustrates how the fire–grazing interaction may control woody encroachment and shows that cattle substitute, at least partially, for endemic


INTRODUCTION
Savannas are characterized by the coexistence of grasses and woody plants, where a discontinuous layer of trees and shrubs intermixes with a continuous grass layer (Van Langevelde et al. 2003). On grasslands, herbaceous plants dominate and woody plant establishment is minimal. Both grasslands and savannas can be converted to shrubland and woodland by uncontrolled woody encroachment (Archer et al. 1995, Van Auken 2000. Grasslands and savannas with ≥650 mm of precipitation are particularly unstable, as this climate allows for woodlands to develop unless there is periodic fire and herbivory to limit populations of resprouting woody plants (Higgins et al. 2000, Sankaren et al. 2005).
Long-lived woody plants store their reproductive potential, allowing herbaceous plants to dominate while woody plant recruitment is suppressed (Higgins et al. 2000). Periodic fire extends these intervals by consigning resprouting woody plants to a fire-trap, wherein established plants survive fire but do not mature within the fire-return interval. Woody plants eventually escape the fire-trap by growing to fire-resistant size (Higgins et al. 2007, Hoffmann et al. 2009). In addition, they can physiologically adapt after fire and herbivory by overcompensating (i.e., growing larger than before) and adopting a pole-like growth form to maximize height Keeley 2005, Ramula et al. 2019). Some woody species also bypass the fire-trap by reproducing vegetatively (Ratajczak et al. 2011, Scasta et al. 2014.
Increasing woody biomass suppresses grasses, causing reductions in fine fuel and less intense fires, which further promotes woody encroachment. Browsing inhibits this feedback loop (Van Langevelde et al. 2003), as consumption kills woody plants or decreases their size by directly reducing height and triggering stem sprouting, a defensive response at the cost of growing tall Bond 2003, Beschta andRipple 2007). Conversely, grazers decrease fire intensity by consuming fine fuel and creating discontinuous fuel beds (Thaxton and Platt 2006, Davies et al. 2016, Scasta et al. 2016). Yet, moderate grazing pressure, combined with periodic fire, reduces woody encroachment more than fire alone (Ratajczak et al. 2014, Venter et al. 2018. The mechanisms driving this pattern are unclear. Coupled fire and grazing are thought to maintain grasslands by interacting through a series of feedbacks which create a shifting mosaic in various stages of recovery. Instituting spatially explicit fires and moderate grazing pressure establishes focal grazing, wherein large herbivores concentrate on the most recently burned patch and uniformly eat the emergent vegetation Engle 2004, Duchardt et al. 2016). Over time, these intensively grazed areas are abandoned by large herbivores for newly burned patches (Scasta et al. 2016, Wilcox et al. 2018). If, instead, a single prescribed fire is applied to an entire site, large herbivores roam the whole area, mostly eating grass and leaving other plants ungrazed (Vinton et al. 1993). These effects on woody encroachment have never been tested. Furthermore, cattle (Bos taurus) are now the dominant large herbivore in Africa and North America , Hempson et al. 2017. Both continents evolved with diverse communities of large grazers and browsers, but cattle only browse about 10% of their diet (Plumb and Dodd 1993). The implications of this transformation are unknown.
We implemented the fire-grazing interaction with prescribed fire and cattle and compared it to two other regimes involving grazing and prescribed fire. We hypothesized that the fire-grazing interaction controls woody encroachment by establishing focal grazing. We predicted that focal grazing would increase browsing, as large herbivores are concentrated where woody plants are resprouting following fire, when they are most palatable, and would mitigate the fuel load reductions usually associated with grazing by reducing grazing pressure on the unburned patches.

METHODS
We sampled 11 sites (15.4-35.0 ha) in the Grand River Grasslands of Ringgold County, Iowa, and Harrison County, Missouri, USA, a 62,000-ha area in the Central Forest-Grasslands Transition Zone, a savanna ecoregion ( Fig. 1; Appendix S1, Text S1). The average annual precipitation in the Grand River Grasslands is 925 mm, one-third falling from June to August, when the average temperature is 22.9°C.
In 2007, each site was divided into three patches of approximately equal area and assigned to one treatment: patch-burn-graze (n = 4); graze-and-burn (n = 4); or burn-only (n = 3; Fig. 1). In the patch-burn-graze treatment, one patch was burned annually on a rotating basis so that the entire site was burned over a three-year cycle with free access by cattle. The patch-burn-graze treatment mimicked the fire-grazing interaction. The graze-and-burn treatment featured one site-wide burn every three years with free access for cattle. The burn-only treatment had one site-wide burn every three to five years and no grazing. Pasture-wide burns were initiated in 2009. All burns on all sites occurred from mid-March to early April.
All grazed sites were fenced around their perimeters. Since 2012, the stocking rate of each site was adjusted annually (average of 2.5 animal unit month (AUM)/ha) on the basis of herbaceous plant biomass remaining at the end of the previous grazing season. The target was to establish moderate grazing pressure by leaving approximately 5000 kg dry mass/ha. Adaptive ❖ www.esajournals.org 3 February 2020 ❖ Volume 11(2) ❖ Article e02993 stocking ensured the grazing pressure relative to the amount of available forage was comparable among sites and among years (Appendix S1: Tables S1, S2).

Vegetation surveys
In 2007, permanent line transects were located to avoid wooded riparian corridors. One to three line transects were established on each patch (3-9 line transects per site; Fig. 2). Each line transect was 100-300 m in length, depending on the size of the patch, and the distance between each transect was ≥150 m. In July of each year, 30 0.5-m 2 quadrats per patch were placed 25 m to either side of each transect (Fig. 2). From 2012 to 2018, we measured litter depth (cm) and litter cover (%) along these transects. Litter was defined as horizontal dead plant material. Litter depth was measured to the nearest 0.25 cm at one random point in each quadrat. Litter cover was visually estimated and recorded as the midpoint of the following categories: 0%, 1-5%, 5-25%, 25-50%, 50-75%, 75-95%, and 95-100%.
In 2017 and 2018, we surveyed woody plants by subdividing the transects into 100-m subsections and establishing 100 9 2 m belt transects approximately overlapping the quadrats (Fig. 2). Within each belt transect, all woody ramets ≥0.5 m tall were measured once from late July to early August, recording species, maximum height (dm), and the number of live and dead stems arising from the root crown. Trees were counted as one ramet, and boles were counted as stems. Ramets (hereafter referred to as plants) are individual woody plants or clonal fragments that are rooted and may arise from vegetative nodes or seeds (Scasta et al. 2014). In 2018, we visually estimated the litter cover to the nearest 10% under the crown of each woody plant encountered on the belt transects (Twidwell et al. 2009).

Analyses
All analyses were performed with R Statistics (R Core Team 2019; Data S1). We estimated the effect of treatment on the woody plant density per belt transect with a generalized linear model with a Gaussian distribution and identity link function. The sampling unit was the belt transect, the independent variable was the number of woody plants of all species, and the dependent variables were, additively, treatment, yearssince-burn, mean elevation, mean slope, mean aspect, and site identity. The correlation between the number of plants per belt transect and site identity was low (|r| = 0.15), indicating that it is reasonable to treat each belt transect as an independent sample with respect to site identity (Dormann et al. 2013; also see Appendix S1: Text S1).
We calculated the mean slope, mean aspect, and mean elevation of each belt transect with the Spatial Analyst tool in ArcGIS version 10.5.1, using a 3-m resolution digital elevation model of Ringgold County and Harrison County from the National Elevation Dataset (USGS 2017).
To evaluate the effect of treatment on woody plant size, we compared the height and the number of live stems, on average and on a per-species basis, of eight species which were common in all three treatments (91% of all plants sampled). ❖ www.esajournals.org 4 February 2020 ❖ Volume 11(2) ❖ Article e02993 Data for the less common species were too sparse relative to the variance of the sample to estimate the effect of treatment. We compared height and the number of live stems, respectively, with generalized linear models with quasi-Poisson distributions and the log link functions. The plant was the sampling unit, and the dependent variables were, additively, species identity, years-sinceburn, and mean slope, aspect, and elevation of the belt transect where the plant occurred. We used the plant as the sampling unit because individuals independently express their phenotypes Bond 2003, He et al. 2017). The marginal effect of treatment was estimated by exponentiating the parameter coefficients (i.e., the odds ratios) and their 95% confidence intervals (CI). For analyzing height, plants >50 m tall were excluded as outliers (0.009% of observations). Then, we estimated the effect of treatment on the height and number of live stems of each woody plant species in the same way, except we used a separate model for each species instead of a covariate for species identity.
To compare the accumulated fuel among treatments, we estimated the effect of treatment on three metrics: patch-scale litter cover, patch-scale litter depth, and litter cover under the crown of each woody plant encountered on the belt transects. Patch-scale litter cover and litter depth were compared, respectively, with one-way ANOVAs with treatment as the dependent variable. The patch was the unit of replication, and the values were the mean quadrat measurements per patch. Third, we used a generalized linear model with a Gaussian distribution and identity link function to estimate the effect of treatment on the litter cover under the crown of each woody plant. Treatment and species identity were the additive dependent variables, and the plant was the unit of replication. For these three analyses, we only included patches or woody plants on patches in the year before prescribed fire was implemented, so all samples had a years-since-burn of two.
To evaluate the effect of treatment on woody plant community composition, we performed a PerMANOVA with the vegan package . We built the community matrix by calculating species frequency per site. Frequency accounts for uneven sampling effort by dividing density (plants/200 m 2 ) by the number of belt transects sampled. Then, the dissimilarity matrix was calculated with Bray-Curtis distances and visualized with rank-abundance curves. Species' traits (Appendix S1: Table S5) were acquired from regional floristic descriptions ) and the PLANTS Database (USDA 2019).

Woody plant community composition
Treatment accounted for 34% of the woody plant community dissimilarity (F 2,10 = 2.02, P = 0.04). Woody species with the ability to spread vegetatively always had the greatest frequency. The burn-only treatment was almost entirely composed of woody species which spread vegetatively. These species decreased in dominance in the graze-and-burn treatment, where long-lived woody species (i.e., ≥25 m ultimate height) increased in richness but maintained low frequencies (i.e., <1 plant per belt transect). The patch-burn-graze treatment was intermediate between the other treatments. Like the graze-and-burn treatment, long-lived woody species had higher richness and low frequency; similar to burn-only, woody species which spread vegetatively dominated ( Fig. 4; Appendix S2: Fig. S2).

DISCUSSION
Our experiment is the first to illustrate how the fire-grazing interaction promotes tree-grass coexistence by controlling woody encroachment. Focal grazing consigned resprouting plants to the fire-trap by increasing the fuel accumulation on unburned areas and increasing browsing on the most recently burned patch. In addition, the movement of large herbivores also appeared to decrease the frequency of woody species which spread vegetatively. Our data run contrary to the theory that grazing increases woody encroachment by decreasing fire intensity (Briggs et al. 2005), given that the burn-only treatment had far more woody plants despite having fuel beds that were both deeper and more continuous than in the graze-and-burn and patch-burn-graze treatments.   (Table 1). This conclusion is directly supported by Browning and Archer (2011), who tested the effects of cattle grazing on shrub encroachment in a semiarid grassland by mapping resprouting mesquite shrubs (Prosopis veluntina) inside and outside of exclosures over 74 yr of fire suppression. Contrary to their expectations, the absence of large herbivores promoted woody encroachment relative to grazed areas. Elsewhere, Ratajczak et al.
compared woody encroachment in a mesic grassland after 28 yr of treatment with either fire alone or fire and grazing by bison (Bison bison). At 4-yr fire-return intervals, the shrub cover increased twice as much in the ungrazed treatment compared to the grazed treatment. Our study adds to this growing line of evidence that moderate grazing by large herbivores controls woody encroachment and is the first, to our knowledge, to show evidence that attributes part of this effect to browsing by cattle Archer 2011, Ratajczak et al. 2014).
When plants are browsed, they can resprout more stems as a defensive response (Bond and Keeley 2005). Resprouting stems temporarily decreases resource allocation to the basal root structures and can delay reproduction or kill woody plants (Bowen andPate 1993, Vil a andTerradas 1995). Both fire-and-grazing treatments increased stem resprouting, so it is possible that combining prescribed fire and herbivory may limit woody encroachment by decreasing plant reproduction or increasing mortality. It is worth noting, however, that the patch-burn-graze treatment resulted in a proportional decrease in plant height whereas the graze-and-burn treatment did not. As a result, the graze-and-burn treatment had the largest plants of any of the three treatments, illustrating that plants in environments with herbivory may more quickly escape the fire-trap unless browsed by large herbivores.
Browsing intensity appeared to depend on fire behavior and, in response, how large herbivores changed their movement and foraging preferences. Coupling fire and grazing promotes intensive foraging on the most recently burned patch (Fuhlendorf and Engle 2004). Reflecting this behavior, nearly all woody plants in the patchburn-graze treatment were shorter than their conspecifics in the burn-only treatment and they decreased in height by about the same amount. In the graze-and-burn treatment, only three woody plant species were shorter. This pattern suggests that focal grazing establishes uniform browsing on the most recently burned patch, but otherwise, cattle only browse relatively palatable species (Hedtcke et al. 2009).
The continuous movement of large herbivores across the graze-and-burn sites also seemed to establish an environmental filter excluding woody species which spread vegetatively (Colgan and Asner 2014). These woody species dominated the burn-only and the patch-burn-graze treatments but were infrequent in the graze-andburn treatment; yet, there was no significant difference between the effects of the graze-and-burn and patch-burn-graze treatments on woody plant density. Thus, focal grazing may allow woody species which spread vegetatively to temporarily invade the unburned patches, but once large herbivores return, they extirpate most of these species (Van Uytvanck and Hoffmann 2009). We therefore propose that the presence and movement of large herbivores, at moderate stocking rates, reduce the frequency of woody plants which spread vegetatively.
Woody species which spread vegetatively facilitate transitions to shrubland and woodland by reproducing within the fire-return interval and forming dense "shrub islands" where fuel accumulation is reduced (Ratajczak et al. 2011). Gaps in the fuel load directly around plants confer fire resilience Clarke 2006, Twidwell et al. 2009), making it impossible to restore grasslands or savannas with fire alone (Van Langevelde et al. 2003, Bestelmeyer et al. 2018). This dynamic manifested in the burn-only treatment, where species which spread vegetatively reached the highest densities and fuel accumulation under the crowns of these plants was suppressed. There was no reduction in fuel load under the woody plants' crowns in either the graze-and-burn treatment or the patch-burngraze treatment, possibly because large herbivores prevented woody species which spread vegetatively from establishing or expanding to sufficient densities to alter fuel dynamics.
Controlling woody encroachment is necessary for grassland and savanna conservation and maintaining productive rangelands. To these ends, recreating the fire-grazing interaction with ❖ www.esajournals.org cattle may be an applicable approach commensurate with the breadth of the problem of woody encroachment ). Yet, this target may be more difficult to achieve in arid regions because desert rangelands are prone to stochastic droughts O'Connor 1999, Hempson et al. 2015). Droughts can kill woody plants and, so, reduce woody encroachment in the short term (Benigno et al. 2014, Pratt et al. 2014. However, the defoliation is locally unpredictable, which increases the risk of overgrazing and, in turn, facilitates the expansion of woody species which are unpalatable to cattle (e.g., honey mesquite Prosopis glandulosa; Hakkila et al. 1987, Bestelmeyer et al. 2018. To avoid these pitfalls, continuous grazing at low stocking rates or strategies, alternative to patch-burn grazing, which promote spatiotemporal heterogeneity, such as pastoral-nomadism or multi-paddock rotational grazing (M€ uller et al. 2007, Briske et al. 2008, Teague et al. 2011, scale the number of large herbivores to the amount of palatable forage and allow or encourage animals to move across the landscape (Oba et al. 2000, Teague et al. 2013).

CONCLUSIONS
Recreating the fire-grazing interaction mitigates the fuel reductions that usually result from grazing and establishes uniform browsing on the most recently burned patch. High fire intensity and browsing act together to maintain resprouting woody plants in the fire-trap (Higgins et al. 2000, Van Langevelde et al. 2003. Furthermore, grazing correlates with more long-lived woody plant species, which typify savannas (Higgins et al. 2000), and excludes woody plants which spread vegetatively, a life-history trait that accelerates transitions to shrubland (Ratajczak et al. 2011). By creating these dynamics, cattle substitute, at least partially, for endemic large herbivores as long as the stocking rate is moderate.

ACKNOWLEDGMENTS
This material is based upon work that is supported by the National Institute of Food and Agriculture, U.S. Department of Agriculture, under award number ILLU-875-918, and by the Iowa Ornithologist's Union and crowdfunders via Experiment.com/ grasslandbirds. We thank TJ Benson, G Spyreas, JJ Coon, SB Maresh Nelson, SJ Halsey, B Clifford, and two anonymous reviewers for comments on the manuscript; and B Ewing, J Donohue, and A Lewanski for their expertise collecting field data. JFC collected and analyzed the data; JFC and JRM wrote the manuscript; JRM, DMD, WHS, and JFC designed the study; and JRM, DMD, and WHS procured most of the funding.

Study region
The Central Forest-Grasslands Transition Zone is an ecoregion which supports North American savannas, characterized as an ecotone between the eastern deciduous forests and the tallgrass and mixed-grass prairies of the Great Plains. The landcover currently is a mosaic of row-crop agriculture and grassland, shrubland, savanna, and woodland habitats (Olson et al. 2001).
Situated within this ecoregion is the Grand River Grasslands, which has been identified as the best-known opportunity to restore a functional tallgrass prairie in the Midwest because 15% of the land is protected and the > 80% is grass ). The topography is rolling hills, with loess soils along ridgetops and glacial till on hill slopes. The predominant land cover within 1 km of our sites was grassland (58%), comprising pastures, hayfields, and native prairie; woodlands (22%); and row crops (18%;

Study design
In 2014-2018, the patch-burn-graze and graze-and-burn sites received one stocking treatment: season-long or early-intensive. In the season-long treatment, sites were stocked from early April to early September (6 months) at a moderate stocking rate. In the early-intensive treatment, the sites were stocked from early April to early July (3 months) at heavier stocking rates. As a result, instantaneous grazing pressure (AU/T DM of plant biomass) is 2-4 times greater in the earlyintensive treatment than the season-long treatment (Maresh Nelson et al. 2019). These treatments tested the response of Schedonorus arundinaceus, the invasive grass tall fescue. In addition, glyphosate herbicide was systematically applied to the graze-and-burn sites to reduce S. arundinaceus and various adaptive management practices were used to reduce other invasive plants (Table S1; Table S2).

Vegetation surveys
The number of belt transects sampled per patch in each year (Table S3) are provided.

Analyses
To model plant density due to treatment, we used generalized linear models rather than mixed models to account for the effects of site identity (Table S4). During preliminary analysis, we tested two types of mixed models which were built with the 'lme' function in the 'nlme' package (Pinheiro et al. 2019). When pasture was incorporated as a repeated measure, the maximum likelihood estimation did not converge on a solution. When pasture was treated as a random variable, the least-squared mean number of plants per transect due to each treatment was the same as the means estimated by the generalized linear model, but the mixed model was imprecise (CI = 27.4-232.9 plants/200 m 2 in the burn-only treatment; −55.8-120.8 in the grazeand-burn treatment; and −52.9-123.1 in patch-burn-graze). Giant confidence intervals signal that simplifying the model is appropriate (Bolker et al. 2008).
When least-squared means due to treatment were calculated, we held all other variables at their average values, and used Tukey's test to compute the significance of the pairwise differences between the means.
There is a small discrepancy in the degrees of freedom between the models for the number of live stems and height because a few samples had missing data (e.g. the plant's height was measured, but due to observer error, its stems were not counted).
The species traits used to evaluate community composition are listed (Table S5). Patch-burn-graze Season-long 1 glyphosate could introduce bias by reducing woody plant density; but note there was no significant difference between density in the patch-burn-graze and graze-and-burn treatments (P = 0.32) and that species which spread vegetatively can do so within the same year as disturbance, giving them time to recover from herbicide 2 disking and herbicide could introduce bias by killing woody plants and reducing litter accumulation; but note, the burn-only treatment still had the highest woody plant density and highest litter accumulation 3 woody plants were not sampled here in 2018 4 potential bias accounted for by the year-since-burn covariate included in the models rancher did not clarify if they stocked 19 or 20 fall cows 2 site was stocked with the same number of cattle as previous year; but the stocking rate cannot be calculated because of unknown turn-in and removal dates and unplanned fence maintenance Table S1: Tukey's test results of the significance of the pairwise differences (P) between the least-squared means due to treatment for woody plant density (plants/200m 2 ).

Treatment
Burn-only Graze-and-burn Patch-burn-graze Burn-only -P < 0.0001 P = 0.0001 Graze-and-burn --P = 0.32 Patch-burn-graze --- Table S2: Tukey's test results of the significance of the pairwise differences (P) between the least-squared means due to treatment for patch-scale litter depth (cm).

Treatment
Burn-only Graze-and-burn Patch-burn-graze Burn-only -P < 0.0001 P = 0.0008 Graze-and-burn --P = 0.06 Patch-burn-graze --- Table S3: Tukey's test results of the significance of the pairwise differences (P) between the least-squared means due to treatment for patch-scale litter cover (%).

Treatment
Burn-only Graze-and-burn Patch-burn-graze Burn-only -P = 0.08 P = 0.12 Graze-and-burn --P = 0.96 Patch-burn-graze --- Table S4: Tukey's test results of the significance of the pairwise differences (P) between the least-squared means due to treatment for the litter cover (%) under the crown of each plant.