Predictive Models of Selective Cattle Use of Large, Burned Landscapes in Semiarid Sagebrush-Steppe

ABSTRACT The fire-exotic annual grass cycle is a severe threat to shrub-steppe rangelands, and a greater understanding of how livestock grazing relates to the problem is needed to guide effective management interventions. Grazing effects vary throughout shrub-steppe rangelands because livestock are selective in their use within pastures. Thus, knowing where cattle are located and concentrate their use in a postfire landscape is important for enhancing plant community resiliency to disturbance and resistance to exotic annual grass invasion. We asked how the distribution and intensity of cattle use varied across 113 000 ha of recently burned, environmentally varied shrub-steppe. Generalized linear mixed effects models were used to determine the relationship of cattle dung (presence/absence and counts), which was recorded during the third to fifth postfire year (after grazing deferment) on 1166 (531-m2) plots, to water sources, burn severity, grass cover, and topographic predictors. Our distribution and intensity of use models revealed similar relationships between cattle use and landscape predictors. Cattle use was greater in areas that were flatter and closer to water and that had moderate burn severity and less heat load and ruggedness. Slope had the strongest effect on cattle use of the predictors. The probability of cattle being present decreased by 10% for every 5° increase in slope until slope exceeded 15°, and then the effect of slope weakened. Despite moderate slopes (X̄ = 14°), cattle use was greater in areas of moderate burn severity, presumably because these areas provided greater perennial grass production. While there was much unexplained variation, these models suggest that cooler climate, water access, topographic factors, and burn severity affect maneuverability to create greater livestock use of certain areas within grazing pastures. Restoration investment planning or assessments and expectations of restoration success could be improved by considering that these livestock hotspots may recover differently from the surrounding landscape.


Introduction
Changes in fire regimes during the past 50 yr are an increasing problem in semiarid and arid ecosystems and one that is expected to worsen under predicted climate change ( Abatzoglou and Kolden 2011 ;McKenzie et al. 2004 ). Interactions among grazing, large fires ( > 40 0 0 0 ha), and exotic annual grasses have degraded native shrub-steppe plant communities across large portions of the western United States (reviewed in Pyke et al. 2016 ). Exotic annual grasses such as cheatgrass (Bromus tectorum) can occupy interspaces between native plants and form highly combustible and continuous fuel beds, resulting in a positive feedback loop with fire (reviewed in Germino et al. 2016 ). Many of the foundational shrub species (i.e., Artemesia spp.) in this ecosystem are fire intolerant, and decades or even centuries may elapse before plant communities are restored to their prefire conditions ( Baker 2006 ). The slow growth of sagebrush and community shifts to exotic annual grasslands have been linked to population declines in sagebrush obligate wildlife species such as the greater sage-grouse ( Centrocercus urophasianus, Coates et al. 2016 ). Moreover, the quality and quantity of forage for cattle decreases as communities shift from native perennial bunchgrass to exotic annual grass ( Germino et al. 2016 ). Given the magnitude and severity of this transformation on natural processes and ecosystem services, broad-scale effort s to rehttps://doi.org/10.1016/j.rama.2022.07.007 1550-7424/© 2022 The US Geological Survey. Published by Elsevier Ltd on behalf of The Society for Range Management. All rights reserved. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ) store native plant communities are being implemented, especially after wildfire ( Pilliod et al. 2017 ), but there is considerable uncertainty as to how the selective use of cattle grazing after fire affects restoration outcomes.
Cattle use is ubiquitous but highly variable in intensity, timing, frequency, and spatial distribution across shrub-steppe rangelands ( Adler et al. 2011 ), resulting in a broad range of effects on vegetation, which can be positive, negative, or neutral ( Jones 20 0 0 ). Where plant-community integrity is high (e.g., < 10% exotic annual grass cover) and climate and other site conditions are favorable ( Chambers et al. 2019 ), moderate grazing intensity before fire may increase resistance to exotic annual grasses, at least in small areas (e.g., 65 ha, Bates and Davies 2014 ;Davies et al. 2016 ). There are also examples where light to moderate grazing had neutral to positive effects on site resistance to annual grass invasion ( Davies et al. 2021 ). However, inappropriate historic and current grazing practices in some areas have depleted native perennials or biological soil crusts and thereby contributed to increases in exotic annual grasses ( Reisner et al. 2013 ;Condon and Pyke 2018 ;Root et al. 2020 ;Williamson et al. 2020 ). Recovering perennial herbaceous plants may be more vulnerable to grazing pressures because their scarcity could lead to greater utilization and their incompletely developed state may render them with greater physiological vulnerability to grazing. Hotspots of grazing use and impacts after wildfire are thus critical to know for promoting resistance to invasion by exotic annual grasses and recovery of resilient perennials after wildfire, in both planning management actions and assessment.
Postfire livestock management strategies that minimize stress to recovering perennial herbs and avoid increases in exotic annual grasses are critical for breaking the exotic annual grass-fire cycle . However, guidelines for postfire grazing management and deployment of restoration investments generally assume implicitly that grazing intensity is uniform within pastures, even though cattle use and grazing impacts are likely unevenly distributed within most pastures. In addition to biotic factors (e.g., forage quality and plant community type), abiotic factors such as water sources, topographic position, and hill slope have a strong influence on how cattle are distributed and how their intensity of use varies within permitted grazing areas ( Holechek 1988 ;Bailey 2005 ;Gersie et al. 2019 ;Raynor et al. 2021 ). Furthermore, cattle are known to shift their distribution when resources are moved intentionally (e.g., salt and minerals) or altered by fire ( Clark et al. 2014( Clark et al. , 2016. Thus, there are significant opportunities to adjust grazing management practices and assist in restoration planning to improve the condition and resilience of shrub-steppe communities that are threatened by exotic annual grasses, particularly after fire ( Burkhardt and Sanders 2012 ;Boyd et al. 2014 ).
Public grazing records provide information on the season of use and intensity (animal unit months [AUMs]) that are useful for investigating grazing impacts at broad spatial scales, such as the large pastures on federal rangelands in the western United States. However, fair assessments of grazing impacts require some ability to measure or predict variable utilization of the available grazing area by livestock. Methods that have been used to quantify relationships between biotic and abiotic factors and cattle use within pastures include, for example, Global Positioning System (GPS) devices, direct observation, plant utilization, and fecal deposits (i.e., dung) ( Gillen et al. 1984 ;Hart et al. 1993 ;Ganskopp 2001 ;Porath et al. 2002 ). Dung has been used to determine relative habitat use, intensity of use, and behavior of cattle ( Porath et al. 2002 ;Sullivan et al. 2012 ;Cole and North 2014 ). This noninvasive method for characterizing cattle use can be conducted with relatively little additional costs during postfire monitoring. Our primary objective was to describe and quantify some of the resources within grazing pastures that influence the distribution of cattle and their intensity of use after fire. To accomplish this, we collected data on cattle dung (i.e., presence/absence and counts) during an intensive restoration and monitoring effort within a postfire shrubsteppe landscape and then modeled the data using generalized linear mixed effects models (GLMMs) with the intent of determining which landscape features best explained variation in cattle use during the spring growing season. Next, we assessed how well our top models fit the data and their accuracy in predicting cattle use across pastures.

Study area
Data were collected between 2018 and 2020 in the Soda Wildfire perimeter, which burned ∼113 0 0 0 ha and 109 grazing pastures along the border of southern Idaho and Oregon in the western United States during the summer of 2015 ( Fig. 1 ). Approximately 60% of the study was mapped as low burn severity, 30% moderate burn severity, < 0.4% high severity, and 10% unburned. The fire-affected area was characterized as arid and mountainous shrub-steppe and comprised primarily sagebrush (e.g., Artemisia tridentata ssp . wyomingensis ), perennial bunchgrasses and forbs, and exotic annual grasses such as cheatgrass (Bromus tectorum). The mean monthly precipitation was 28.9 mm (range: 10.3-45.9 mm), and the mean monthly temperature was 10 °C (range: 6.5-13.5 °C) across the study area ( Fig. 2 ). Mean elevation was 1286 m (range: 729-1947 m), and mean slope was 11 °(range: 0-60 °). The Bureau of Land Management (BLM) implemented various combinations of herbicide and seeding treatments after the fire totaling > 10 0 0 0 0 ha. Pastures were rested from livestock grazing for ∼2 to 3 yr after the fire to allow plants to recover and increase resistance to exotic annual grasses and resiliency to disturbances ( US DOI 2007 ). Cattle grazing mostly occurred during the spring, summer, and fall at a stocking rate of approximately 0.25-0.50 AUM per hectare. Except for these temporary periods of rest and reductions in AUMs in some pastures after grazing resumed, the grazing season of use and intensity was consistent before and after fire.

Data collection
Cattle dung was measured in permanent 13-m radius field plots (531 m 2 ) established in 2016 across the study area. Plot location coordinates were generated via a stratified-random method (i.e., by pasture) at 1 plot per 54.5 ha or denser ( n = > 20 0 0 plots total) and were moved if they overlapped roads, had > 20% trail area within an 18-m radius, or fell within 0.40 km of water troughs or ponds (see Fig. 1 , Applestein et al 2018 ). The number of cattle dung were counted in plots between 28 April and 20 August each year after cattle were in pastures for at least 15 d. The current years' cattle dung was identified on the basis of its dryness using a "kick test" where cattle dung, when kicked, was either displaced (i.e., previous year) or not displaced (i.e., current year). For this analysis we used data that were collected from 1166 plots for which cattle were in pastures between 1 January and 30 June (i.e., spring grazing).
The percent cover of perennial bunchgrasses and exotic annual grasses was determined from 6-m 2 rectangular aerial photographs at plot center captured from nadir at 2-m height (with Nikon Coolpix AW130, 16 megapixel, Nikon Inc., Melville, NY). We analyzed species cover using grid-point intercept (GPI) software (Sam-plePoint v. 1.43, 100 points image −1 ; Booth et al. 2006 ).
We attained several spatial variables (i.e., topographic, burn severity, water sources, and grass cover) that have been reported or suspected to influence cattle use. Specifically, we used GIS (Ar-cGIS; ESRI Corporation, Redlands, CA) to extract values of elevation, slope, aspect, heat load, ruggedness, burn severity from 30-m pixels, distance to water from 10-m pixels, and percent cover of perennial bunchgrass and exotic annual grass from 32 m × 32 m pixels for each monitoring plot. The US Geological Survey's digital elevation models (30-m pixels) were used to calculate values for elevation (m), slope ( °), and aspect ( °). A burn severity map, which represented low, moderate, and high burn severity and unburned areas, was attained for the Soda Wildfire through the Monitoring Trends in Burn Severity interagency program ( www.mtbs.gov ).
Burn severity categories primarily represent significant differences in vegetation biomass, which are derived from Differenced Normalized Burn Ratio images via Landsat data ( Eidenshink et al. 2007 ). The BLM's range improvement GIS layer was used to determine active anthropogenic water sources. We also used a spatially continuous model to identify which streams had the highest probability of annual streamflow permanence (PROSPER, Jaeger et al. 2019 ). After exploring the models and consultation with local experts, we decided to include streams that had a probability of > 0.30 of having water for our analysis. We combined the anthropogenic water and streams layers and then used Euclidean distance in ArcGIS PRO (because it allowed us to confine our measurements within pastures) to create a 10-m pixel raster, which represented the distance from anthropogenic water and streams within each pasture. Heat load (30-m pixels), which is a unitless value ranging from 0 (low incident radiation) to 1 (high incident radiation), was derived from values of aspect, slope, and latitude ( Welty and Jeffries 2018 ). Vector ruggedness measure (30-m pixels), another unitless value ranging between 0 (no terrain variation) and 1 (complete terrain variation), was calculated from variability in slope and aspect ( Welty and Jeffries 2018 ). Lastly, we used the near distance tool in ArcGIS to calculate the Euclidian distance between monitoring plots and roads and pasture fences.
For each functional group (perennial bunchgrass and exotic annual grass) and year, we conducted an empirical Bayesian kriging regression (combines ordinary least-squares regression and kriging) in ArcGIS Pro at 32 × 32 m pixel resolution with the following covariates used as predictors: mean annual precipitation (PRISM), heat load (calculated from US Geological Survey digital elevation models), number of herbicide treatments, and number of drill treatments. In the empirical Bayesian kriging routine in ArcGIS Pro, predictor variables are transformed into their principal components before the least-squares regression. The dependent variable (perennial bunchgrass or exotic annual grass) was log empirically transformed (with a constant -0.01 added to all values to remove zeros), the semivariogram model type was K-Bessel, the minimum neighbors was set at 5, and the maximum neighbors was set at 15. The resulting maps were constrained to be between 0% and 100% cover  ).

Cattle distribution model
GLMMs were used to model the distribution of cattle and their intensity of use within the Soda Wildfire. For our "distribution model," we used a binomial distribution with a logit link to model the binary response of cattle dung being either present or absent in plots. For our "intensity model" we initially developed a GLMM with a Poisson distribution to model the amount of cattle dung in plots. However, our count data were overly dispersed (dispersion ratio 2.7; P < 0.001; check_overdispersion function; performance package; Lüdecke et al. 2021 ); therefore, we used a negative binomial model instead of the Poisson model. We modeled both re- sponse variables as a function of the 10 explanatory variables (i.e., predictors) described earlier, plus two additional variables to assess nonlinear effects of elevation and aspect and included random effects of survey year and monitoring plot identity using the "glmmTMB" function and glmmTMB package ( Brooks et al. 2017 ) in Rversion 4.1.2 ( R Core Team 2021 ). Burn severity was a categorical variable and unburned was the reference category. We also removed the high burn severity category from the analyses because it represented a trivial portion of the study area ( < 0.4%). Correlation among predictors did not exceed our a priori threshold ( > 0.60); therefore, we used all 11 variables in the analysis. We selected the best predictive models by developing full GLMMs with all 12 predictors for each response variable (i.e., two separate models). Then we removed the predictor with the highest coefficient P value and compared the nested models (full model vs. reduced model with predictor removed) using Akaike Information Criteria with a bias correction term for small sample size (AIC c ; Burnham and Anderson 2002 ;Burnham et al. 2011 ). The model with the lowest AIC c was retained, and the process was repeated until removing a predictor resulted in a higher AIC c value and all predictors were statistically significant at P < 0.05.
We assessed goodness of fit for our final models with marginal (predictors) and conditional (random effects) R ² values ( Nakagawa and Schielzeth 2013 ). Predictive accuracy of our final models was assessed using k-fold cross validation using the train function in the Caret package ( Kuhn 2021 ). We randomly partitioned the data where 80% of the data was used to train the model and 20% of the data was withheld to test the model. We used repeated cross validation with nine folds and repeated the process 100 times. Predictive accuracy measures from k-fold cross validation for the distribution model were accuracy and Cohen's kappa, and for the intensity model were root mean squared error and the coefficient of determination ( R ²). Models with high Cohen's kappa, accuracy, and R ²; and low root mean squared error values were considered to have high predictive accuracy.

Results
From 2018 to 2020 we observed cattle dung on 511 of the 1166 plots. Cattle selected areas that had gentler slopes and ruggedness, were closer to water sources, had lower heat loads, and had moderate burn severity ( Fig. 3 , Table 1 ). Slope had the greatest effect on the probability of cattle presence of all the predictors ( Fig. 4 ). The relationship between cattle presence and slope was moderated as slope increased. For example, the probability of cattle being present decreased by 10% for every 5 °increase in slope until slope exceeded 15 °, and then the effect of slope weakened  ( P = < 0.001, see Fig. 4 ). Cattle avoided areas farther from water sources. The probability of cattle being present decreased ∼3% for every 500-m increase in distance from water sources ( P = 0.008, see Fig. 4 ). Cattle also exhibited negative associations with both heat load and ruggedness (see Fig. 4 ). The probability of cattle being present decreased ∼15% for every 0.5 unit increase in heat load ( P = 0.024) and ∼9% decrease for every 0.02 unit increase in ruggedness ( P = 0.029). Cattle selected areas that were classified as moderate burned severity more than unburned or low burned severity ( P = 0.031), whereas cattle did not differentiate between unburned and low burned severity ( P = 0.659, Fig. 5 ).
The intensity of cattle use based on abundance of dung shared the same relationships as cattle presence ( Fig. 6 , see Table 1 ). There were strong nonlinear associations between cattle use intensity and slope (see Fig. 6 ). The effect of slope was greatest until 20 °, at which point the effect weakened and the predicted number of dung was reduced (see Fig. 5 ). For example, the predicted numbers of dung amounts per plot were 4.53, 1.57, 0.54, and 0.19 on 0 °, 10 °, 20 °, and 30 °slopes, respectively. Complete model selection results for the distribution and intensity models are in Table S1 (available online at […]).
The variation in cattle use explained by both models was low. The distribution model explained 32% (conditional R 2 ) of the variation in the probability of cattle use with 11% of the variation explained by the predictors (marginal R 2 ). For the intensity model, the predictors explained 19% of the variation and the random effects of survey year and monitoring plot identity explained another 46% of the variation in cattle use. The accuracy and error for predictions of cattle use across the study area were not strong, with an accuracy of 0.61 and Cohen's kappa of 0.20 for the distribution model, and a root mean square error of 12.61 and r 2 of 0.03 for the intensity model.

Discussion
Heterogeneity in fire-grazing interactions and their relationship to plant community transitions is well known in rangelands such as tall-grass steppe ( Fuhlendorf et al. 2009 ), but there are fewer studies in shrub-steppe and the interactions might be expected to differ strongly. In grass-dominated systems, which have relatively high herbaceous production, fire-grazing interactions (i.e., pyric herbivory) are essential to ecosystem structure and function and can reduce invasions of non-native species ( Starns et al. 2019 ;Sherrill et al. 2022 ). Moreover, maintaining ecological processes in grassland systems through disturbance-driven heterogeneity is a common application for managing cattle and wildlife ( Allred et al. 2011 ;Hovick et al. 2014 ;Wilcox et al. 2021 ). However, shrubsteppe has relatively less herbaceous production than grasslands, and frequent disturbances reduce resistance to exotic plant invasions and thereby facilitate shifts from native perennial to exotic annual grass −dominated communities ( Bradford and Lauenroth, 2006 ;Condon and Pyke 2018 ). Previous studies have revealed cattle to be selective of landscape features in shrub-steppe and related rangelands, and here we revealed selection in the postfire environment, demonstrating that our models may represent cattlehabitat relationships at broader spatial scales. Our findings suggest resource selection by cattle in a postfire landscape deviates slightly from unburned rangeland, perhaps owing to increased production of perennial bunchgrasses in burned areas, but the influence of topographic features and water on cattle use strongly persists. Still, both the distribution and intensity model did not have strong explanatory power and predictive accuracy, indicating mechanisms at finer spatial scales (e.g., unmeasured variation in forage quality) may be driving variation. More direct methods for quantifying cattle use (e.g., GPS) are needed.
Topography has a strong influence on the distribution and intensity of cattle use across a wide range of ecosystems, including shrub-steppe rangelands ( Mueggler 1965 ;Clark et al. 2014Clark et al. , 2016Liao et al. 2018 ;Raynor et al. 2021 ). The relatively strong effect of hill slope and ruggedness on cattle distribution and intensity reported here is consistent with these previous studies. Cattle use generally declines with increasing hill slope, yet there are instances where cattle will use moderate slopes more frequently, such as recently burned areas ( Clark et al. 2014( Clark et al. , 2016. However, cattle tend to avoid slopes exceeding 30 ° ( Gillen et al. 1984 ;Bailey 2005 ) and where topographic variation is high, therefore limiting their distribution and grazing capacity in steep rugged terrain ( Bailey et al. 2015 ;Gersie et al. 2019 ;Raynor et al. 2021 ).
Cattle are dependent on water and concentrate around these resources within grazing pastures, especially in semiarid and arid environments ( Ganskopp 2001 ;Bailey 2005 ;Pinchak et al. 1991 ). In the western United States, much of the utilization by cattle occurs within 1.6 km of water sources and areas beyond 3.2 km are considered ungrazable in many contexts ( Valentine 1947 ; Holechek  1988 ). Consistent with these patterns, we observed an effect of water where the probability of cattle use declined with increasing distance from water. In the western portion of our study area, Clark et al. (2014Clark et al. ( , 2016 examined habitat use of cattle during the spring growing season for 5 yr after prescribed fire and reported a negative relationship between water and cattle use during the first year after fire, a positive relationship during yr 4 and, a nonsignificant effect of water the other 3 yr. These authors speculate that fire may have played a role in these conflicting use patterns by shifting preferred feeding sites within pastures ( Clark et al. 2014( Clark et al. , 2016. Heat load had a moderate, negative effect on both the distribution of cattle and their intensity of use. Heat load values repre-sent potential direct incident radiation and are derived from values of aspect, slope, and latitude ( McCune and Dylan 2002 ). High heat loads are characterized as southwest-facing steep slopes at low latitudes, whereas low heat loads, which cattle in our study area preferred, are described as northeast low-angled slopes at high latitudes ( McCune and Dylan 2002 ). We are not aware of this predictor being used in any previous studies involving cattle use, but its usefulness in our models might be mirrored in other, topographically varied landscapes.
Heterogeneity in forage quantity and quality influences where cattle graze and how much time they spend grazing in an area (see reviews Senft et al. 1987 ;Bailey 2005 ). Increases in the quantity and quality of forage after fire ( Hobbs and Spowart 1984 ;Bates et al. 2009 ) can shift the distribution and increase utilization of burned areas by large herbivores ( Boindini et al. 1999 ;Van Dyke and Darragh 2007 ). In our study area, cattle use was greater in areas where the burn severity was moderate. These areas were also characterized as having steeper slopes ( ∼4 °) and greater cover of perennial bunchgrasses ( ∼5%) than unburned areas or areas of low burn severity. Although cattle use was negatively associated with slope across the study area, cattle deviated from this relationship and used moderately burned areas that were on average 14 °, perhaps because they provided greater forage. Despite decreases in forage quality 1-2 yr post fire ( Hobbs and Spowart 1984 ), forage production often remains high for > 7 yr ( Wambolt et al. 2001 ;Bates et al. 2020 ). Our findings provide further evidence from Clark et al. (2014Clark et al. ( , 2016 that fire can indirectly influence cattle use for at least 5 yr in shrub-steppe rangelands; however, the mechanisms driving these patterns require further investigation. The low amount of variation explained in both response variables and weak prediction of our final models reflect the general challenge of predicting cattle use. These challenges arise because the amount of usable space for cattle is determined by many biotic and abiotic factors and these factors are influential at multiple spatial scales ( Senft et al. 1987 ). Cattle use is a hierarchal process: at broad spatial scales, landscape features such as the location of water, hill slope, and forage patchiness contribute to where cattle will feed (i.e., feeding area selection), whereas at finer spatial scales (i.e., within feeding area selection), their use is driven by characteristics of forage quality and quantity ( Bailey et al. 2015 ). Furthermore, reproduction, life stage, breed type, and cognitive and social behavior can all contribute to selection processes ( Sahu et al. 2020 ). Thus, our models may represent how a limited number of abiotic factors influence the selection of feeding areas by cattle (albeit poorly), but to fully understand the mechanisms driving selection processes and their impact of postfire vegetation recovery, other biotic factors should be examined ( Bailey 2005 ).

Implications
Grazing management is an integral component of multiple use natural resource management; therefore, knowing how cattle use varies across large pastures can assist with natural resource planning (e.g., fuels, restoration). Here, we demonstrate how resource heterogeneity contributes to the selective use of cattle within pastures. Restoration planning could be improved with consideration of these combined effects on cattle use during their postfire planning effort s. For inst ance, the models presented here could be used to develop maps that could specify areas that may need longer periods of rest (high use areas) or areas that may be grazed earlier (low use areas). Furthermore, an understanding of how grazing impacts are distributed across fire affected landscapes is needed to fairly evaluate restoration outcomes. Our distribution and intensity models were reasonable yet provided a modest gain in explaining the variability in cattle use of a large, recently burned landscape based on evidence of cattle dung. Recent advances in DNA metabarcoding and GPS telemetry tracking could provide additional knowledge for understanding variation in grazing impacts. Diet composition revealed by fecal DNA metabarcoding would help elucidate relationships between restoration plant species and preferred forage of cattle. GPS collars would provide robust sample sizes of individual cattle locations at fine temporal resolutions ( ∼10 min) to comprehensively represent population-level relationships with both fine-scale (e.g., forage quality) and landscapescale predictors. These methods and others are likely necessary to achieve higher explanatory power and predictive accuracy that is required to comprehensively understand grazing-plant interactions after fire.

Declaration of Competing Interest
All contributing authors have approved this manuscript and have agreed to Rangeland Ecology and Management's submission policies. The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. This manuscript has not been previously published and is not currently in review by another journal.