Spatial data collection and processing

We downloaded publicly available boundary data for all BLM grazing allotments, BLM field offices, and EPA ecoregions from government websites ( Appendix A, Table A1 (mmc1.pdf)). We used geoprocessing tools in R and QGIS to remove duplicate allotment boundaries and removed allotments < 1 ha in area (R Core Team 2022). We assigned each allotment to one of nine ecoregions derived from the EPA ecoregion levels I–III: AZ/NM Highlands, Marine West Coast Forests, Eastern Cold Deserts, Forested Mountains, Northern Great Plains, Southern Great Plains, Mediterranean California, Warm Deserts, and Western Cold Deserts (see Fig. 1;  Table A2 (mmc1.pdf)). We excluded allotments from the Marine West Coast Forest ecoregion from further analyses due to small sample size and area (see Table 1).

We use BLM grazing allotments as a practical unit for summarizing local-scale rangeland trends in our analysis. While allotment boundaries are delineated by the BLM grazing program, the public lands within them are open to and managed for multiple uses. We do not evaluate the effects of livestock grazing or any other specific land use. Allotments often include private rangelands surrounding BLM parcels, and we used data from the US Geological Survey Protected Areas Database (PAD-US) to find the area of BLM and private lands within each allotment. Allotments averaged 4 491 ha in area and comprised 65% BLM, 24% private, and 9% other public lands (see Table 1). We retained private lands within allotment boundaries in our analysis so that trends are representative of these local administrative units. After processing, the allotment data included 21 011 grazing allotments totaling 94.37 10⁶ ha of rangeland (see Table 1).

We used estimates of vegetation cover and aboveground herbaceous production from the Rangeland Analysis Platform (RAP). The RAP uses machine learning to estimate cover of annual grasses and forbs (annuals), herbaceous perennial grasses and forbs (perennials), shrubs, trees, and bare ground, as well as aboveground production of annuals and perennials in kg ha^–1 on US rangelands from 30-m RS data (Allred et al. 2021; Jones et al. 2021). As described in Jones et al. (2018), RAP cover estimates correspond to first-hit cover as measured by the AIM and NRI protocols and species are grouped according to characteristics in the US Department of Agriculture (USDA) plants database (USDA, NRCS 2018). In the remainder of the paper, we refer to RAP cover and production categories (i.e., annuals, perennials, shrubs, trees, and bare ground) as vegetation functional groups. We masked out nonrangelands and areas of hay, alfalfa, and idle cropland within allotments using the National Land Cover Database (Homer et al. 2020) and the USDA National Agricultural Statistics Service Cropland Data Layer (USDA NASS 2016). After masking, we calculated allotment-level average cover and production using the RAP v.3.0 dataset for each yr from 1986 to 2020. Masking and allotment-level vegetation averages were calculated in Google Earth Engine (Gorelick et al. 2017).

Modeling long-term vegetation trends

We analyzed long-term trends in vegetation cover and production using separate linear mixed-effects models for each functional group. We fitted fixed effects for ecoregion and time (calendar year) and their interaction (year by ecoregion). Mixed-effects models were fitted with the lme4 package (Bates et al. 2015) using the following model formula in lme4 notation:

where y is yearly allotment-level cover or production, E is ecoregion, t is year, O is BLM field office, and A is allotment. We included random intercepts and slopes at the level of allotment and BLM field office nested within the ecoregion. For each model, we included 30 yr of data from 1991 to 2020. We log-transformed cover and production values before analysis. To avoid log-transforming values of zero, we excluded allotments that for any year had average cover of < 0.25%; in the herbaceous production models, we excluded allotments that for any year had production of < 0.25 kg ha^–1. This resulted in dropping several hundred allotments from most analyses and several thousand from the annual cover and production and tree cover models (see  Table A3-A9 (mmc1.pdf)). Adding a small increment to cover and production values before log-transformation and including all allotments in the analyses resulted in similar trend estimates. To avoid identifiability problems, we excluded a handful of allotments from field offices that had only a single allotment. After fitting models, we extracted trend estimates and 95% confidence intervals for each ecoregion using the “emmeans” package in R (Lenth 2022). We interpret trends as significant when 95% confidence intervals for an ecoregional trend do not overlap zero.

Spatial heterogeneity in vegetation trends across scales

To quantify spatial variation in vegetation trends across scales, we fit separate linear-mixed effects models for cover and production of each functional group as described in the Modeling long-term vegetation trends section but with an additional term to capture variance in trends at the pixel scale:

where P identifies individual pixels as a grouping factor. Fitting this model to all pixels within all BLM allotments was prohibitively slow, so we used a subset of data. We randomly selected 2069 allotments (∼10% of total) and extracted yearly cover and production values from eight randomly selected 30-m pixels within each (n = 16 541). As in the allotment-level analysis, cover and production values were log-transformed before analysis, and we excluded pixels that had fewer than 20 yr of nonzero cover and production data. After fitting the models, we extracted group-level variance estimates for vegetation trends (slopes) at the field office, allotment, and pixel-scales.

Correlations between vegetation trends across scales

At each of the four spatial scales, we examined correlations between cover trends for annuals, herbaceous perennials, bare ground, and shrubs. Trend estimates for each ecoregion, field office, and allotment were drawn from models fitted in the Modeling long-term vegetation trends section. Trend estimates for pixels were drawn from the models fitted in the Spatial heterogeneity in vegetation trends across scales section. We used random effects slopes (Best Linear Unbiased Predictions) as trend estimates for individual field offices, allotments, and pixels. Since random effects were nested, (pixels in allotments, allotments in field offices and field offices in ecoregions), vegetation trend estimates at each scale represent the residual trend not captured by trends at larger spatial scales.

Supplemental analyses of annual invasion and woodland expansion

We performed three post-hoc analyses to further investigate trends in annual plants and trees. First, we quantified the number and proportion of allotments with average cover and production of annuals greater than that of perennials for each ecoregion for each decade. Second, we conducted a follow-up nonparametric test of tree cover trends among all allotments (including those below the 0.25% cover threshold used in the Modeling long-term vegetation trends section). We report the number and proportion of allotments in each ecoregion with positive and negative trends in tree cover as determined by the Mann-Kendall test for significant trends (P < 0.05). Thirdly, we used field data from the AIM program to identify which annual and tree taxa were most common in each ecoregion. This does not directly reveal which species are responsible for trends but can provide context for interpreting RS trends in overall annual and tree cover.

Status and trends of vegetation cover and production by ecoregion

Herbaceous perennial vegetation dominated allotments in the AZ/NM Highlands, Northern Great Plains, Southern Great Plains, Forested Mountains, and Western Cold Deserts (Fig. 3). Bare ground cover was dominant or codominant in Eastern Cold Deserts and Warm Deserts, while annuals dominated Mediterranean California. Average shrub cover of allotments was > 10% in all ecoregions except for the Northern Great Plains. Average tree cover was below 10% in all ecoregions except for the Forested Mountains (∼15%). Average herbaceous perennial cover and aboveground production were generally several times that of annuals in all ecoregions except Mediterranean California, where annuals were dominant (see Figs. 3 and 4).

Figure 3.

Mean vegetation cover of Bureau of Land Management allotments from 1991 to 2020. Each column shows a separate ecoregion and each row a separate functional group. Colored ribbons show interquartile range of allotment-level average cover. Trend lines from the mixed effects models are shown where significant. Y-axis is on a log scale.

We found positive trends in annual cover and aboveground production in all ecoregions (see Figs. 3 and 4). In contrast, there were declining trends in perennial cover in all ecoregions except the Northern Great Plains—with the most rapid declines in Mediterranean California and Warm Deserts (see Fig. 3). Perennial production trends were negative in the Eastern Cold Deserts, Forested Mountains, Mediterranean California, Western Cold Deserts, and Warm Deserts (see Fig. 4). Bare ground cover was declining in four ecoregions, most rapidly in the Northern Great Plains and Western Cold Deserts, but has increased slightly in Mediterranean California (see Fig. 3). Woody shrub cover declined in the Western Cold Deserts but increased in the Northern Great Plains and Mediterranean California (see Fig. 3). Since we modeled log-transformed cover and aboveground production, the trend coefficients (slopes) shown in Figure 5 can be interpreted as annual proportional rates of change. For example, a 0.028 rate of increase for annual grass cover in the Western Cold Deserts (see Fig. 5) indicates an increase of 2.8% per yr, equivalent to a doubling of cover in 25 yr.

The number of allotments with annual cover and production greater than herbaceous perennials has increased over time in the AZ/NM Highlands, Eastern Cold Deserts, Forested Mountains, Mediterranean California, Western Cold Deserts, and Warm Deserts (see Tables A10, A11). And across all ecoregions, the proportion of allotments with annual cover surpassing that of perennials increased from 4.9% (n = 1 025) in the 1990s to 11% (n = 2 320) dur-ing the past decade. This shift to annual dominance was most dramatic in the Western Cold Deserts of Nevada and surrounding states, where more than a third of allotments (n = 1 277), representing some 18.2 million ha, now have more annual than perennial cover. Likewise, nearly one quarter of allotments in the Western Cold Deserts now have more production from annuals than perennials (see  Table A11 (mmc1.pdf)). The AIM data from 31 908 plots located within 6 744 allotments showed that exotic bromes such as cheatgrass (Bromus tectorum) and Japanese brome (Bromus japonicus/arvensis) were common in many ecoregions (see  Table A12 (mmc1.pdf)). For example, in the Western Cold Deserts, Northern Great Plains, and Eastern Cold Deserts, exotic annual bromes were found in 94%, 85%, and 65% of allotments, respectively. Native annuals were dominant in three ecoregions: the AZ/NM Highlands, Southern Great Plains, and Warm Deserts.

Figure 4.

Mean aboveground herbaceous production of Bureau of Land Management allotments from 1991 to 2020. Each panel shows a separate ecoregion. Colored ribbons show interquartile range of allotment-level average aboveground production for annuals and herbaceous perennials separately. Trend lines from the mixed effects models are shown where significant. Y-axis is on a log scale.

The average tree cover of allotments increased in the Eastern Cold Deserts, Forested Mountains, Northern Great Plains, and Western Cold Deserts (see Fig. 3). In the Northern Great Plains, tree cover increased at a rate of 2.2% per year, nearly doubling from 1.8% cover in the early 1990s to 3.5% cover in the past few years (see Fig. 5). Nonparametric tests showed that 52% of all allotments had increases in tree cover and only 5% had decreases (see  Table A13 (mmc1.pdf)). In the Northern Great Plains and Forested Mountains, tree cover increased in 65% and 60% of allotments, respectively. Notably, Mediterranean California was the only region where there were more declines than increases in tree cover. Field data from AIM showed that native conifers were the most common trees in many ecoregions (see  Table A14 (mmc1.pdf)). Junipers (e.g., Juniperus monosperma, J. osteosperma, J. occidentalis, J. scopulorum) were found on 30–50% of allotments in the AZ/NM Highlands, Eastern Cold Deserts, Forested Mountains, Southern Great Plains, and Western Cold Deserts and on nearly one fifth of allotments in the Northern Great Plains. Mesquite species (Prosopis glandulosa, P. veletuna) were also prevalent in the AZ/NM Highlands, Warm Deserts, and Southern Great Plains.

Spatial heterogeneity in vegetation trends at field office, allotment, and pixel-scales

Average vegetation trends reported for each ecoregion bely considerable variation in trends at field office, allotment, and pixel scales (Fig. 6). Except for annuals, trends varied more between pixels within the same allotment than between allotments in the same field office, or between field offices in the same ecoregion. Trend variance was especially pronounced for tree cover trends, and 59% of group-level variance occurred at the pixel scale. Trends for annual cover and production were unique in having proportionally more variance at the allotment scale and less at the pixel scale. Compared with cover trends, production trends showed a higher proportion of variance at the field office scale.

Correlations between vegetation trends across scales

We further examined correlations between trends for different functional groups (Fig. 7). Annual cover trends were negatively correlated with bare ground trends at the ecoregion, field office, allotment, and pixel scales, with the strongest correlations at the ecoregion and field office scales (see Fig. 7A–D). Annual and herbaceous perennial trends were positively correlated at the field office, allotment, and pixel scales (see Fig. 7F–H). Annual trends were weakly negatively correlated with shrub trends at the allotment and pixel scales (see Fig. A1A). Bare ground trends were negatively correlated with perennial trends and positively correlated with shrub trends at subecoregion scales (Fig. A1B and C). The correlation between shrub and perennial trends was negative at the field office scale but switched to positive at the pixel scale (Fig. A1D).

Trends in vegetation cover and production by ecoregion

Our BLM-specific analyses reinforce the vexing challenge of managing annual plant invasions on western rangelands (Maestas et al. 2022; Tarbox et al., 2022). The continued increase in annual cover and loss of bare ground is evidence that rangelands are accumulating fine fuels that can further increase wildfire (Davies et al. 2012; Smith et al. 2022a). Exacerbating frequency and severity of wildfires can devastate ecosystems by removing keystone plants (e.g., big sagebrush Artemisia tridentata) and imperiling endemic wildlife such as greater-sage grouse (Centrocercus urophasianus; Coates et al. 2016). Even though exotic annuals began invading western rangelands over a century ago (Mack 1981), their continued spread shows little evidence of slowing (Smith et al. 2022b).

Figure 5.

Vegetation cover and aboveground production trends for Bureau of Land Management allotments from 1991 to 2020 estimated by mixed-effects models. Trend magnitudes are shown on the x-axis surrounded by 95% confidence intervals. Trends represent annual rate of change in log cover (left) or log aboveground production (right).

While annuals are not the dominant vegetation group in most ecoregions, a growing number of allotments now have annual cover and production exceeding that of perennials. Today, annuals are the dominant herbaceous component on 1 691 allotments encompassing 21 million ha of rangeland (Tables A10 and A11). The most dramatic rise in annuals was in the Western Cold Deserts of Nevada and parts of adjoining states where annual aboveground production more than tripled from roughly 90 kg ha^–1 in the early 1990s to > 300 kg ha^–1 over the past few years (see Fig. 4). This corresponds to an increasing trend of nearly 5% per year and highlights the growing amount of fine fuel available for destructive wildfires (see Fig. 5). In ecoregions where annual cover was increasing most rapidly, exotic bromes such as cheatgrass were the most common annual species in AIM field data ( Table A12 (mmc1.pdf))—strong evidence that these species are driving increases in annual cover. Climate change, specifically regional warming, may be facilitating the invasion of exotic annuals into higher elevations and cooler climates (Compagnoni and Adler 2014; Blumenthal et al. 2016; Smith et al. 2022b). However, hotter and drier conditions may be detrimental at lower elevations (Larson et al. 2017).

Figure 6.

Estimated variance in vegetation trends at three spatial scales. Text within bars shows percent of group-level variance found at each scale. Trend variances are estimated from mixed effects models fit to log-transformed cover or aboveground production.

Another more subtle yet concerning set of results were widespread loss of herbaceous perennial cover and production (see Fig. 5). Perennial grasses and forbs play a central role in forage production, carbon storage, biodiversity, and as indicators of overall rangeland health (West 2003; Pellant et al. 2005; Koteen et al. 2011). The epicenters of perennial loss were in Mediterranean California and the Warm Deserts (see Fig. 5). In Mediterranean California, perennial production has declined at a rate of 5.9% per year and is roughly one seventh what it was 30 yr ago (see Figs. 4 and 5). Our findings are consistent with a loss of cover and mortality of perennial herbs and bunchgrasses after prolonged drought in the Southwest and elsewhere (Munson et al. 2013; Brookshire and Weaver 2015; McIntosh et al. 2019; Winkler et al. 2019; Williams et al. 2022). The exception to perennial declines was in the Northern Great Plains, where a lengthened growing season and increased CO₂ concentrations could be offsetting the detrimental effects of heat or drought (Reeves et al. 2014; Hufkens et al. 2016; Brookshire et al. 2020).

Figure 7.

Correlations between cover trends (x-axis) for annuals and bare ground (A–D) or between annuals and herbaceous perennials (E–H) at four spatial scales. Trends have been rescaled between 0 and 1 to facilitate visual comparison. Trends for ecoregions, field offices, and allotments are estimated from mixed-effects models fit to allotment-level average cover. Pixel-level trends are from the mixed-effect models fit cover at the individual pixel-scale. Blue lines and Pearson's r are shown for pairs of trends with a significant correlation at P < 0.05 level.

Despite a strong shift from perennial to annual dominance in many allotments (Tables A10–11), we did not find strong evidence that increases in annuals were directly responsible for declines in perennials. In general, areas with more rapid increases in annuals tended to show slower declines in perennials (positive correlations in Fig. 7E–H). This result is counter to the pattern we would expect if these groups strongly limited one another's abundances (Humphrey and Schupp 2004; Parkinson et al. 2013; Reisner et al. 2013; Chambers et al. 2016). To explain this, we hypothesize that perennials and annuals often respond similarly to variation in climate. For instance, if both perennials and annuals are water limited, they may show correlated trends across regions that have become drier or wetter over the past three decades (Munson et al. 2013; Larson et al. 2017). Positive correlations between trends were strongest at the field office scale, (see Fig. 7), a pattern that supports a role for climate variation rather than local interactions.

The proliferation of trees on BLM rangelands (see Fig. 3) is consistent with other studies showing the widespread expansion of trees into rangelands, both in the United States and worldwide (Nackley et al. 2017; Filippelli et al. 2020; Morford et al. 2022). In arid regions of North America, the expansion of native conifers into rangelands has been attributed to wildfire suppression and natural and anthropogenic climate change (Miller et al. 2008; Romme et al. 2009). Tree cover expansion can reduce plant species richness (Ratajczak et al. 2012) and habitat for wildlife that depend on treeless rangelands (Archer and Predick 2014). For example, even low abundances of tree cover may result in the extirpation of grouse (Lautenbach et al. 2017; Nackley et al. 2017) and other grassland birds (Brennan and Kuvlesky Jr. 2005). Given broad concern over the 50% reduction of grassland bird populations in the United States (Rosenberg et al. 2019), tracking further expansion of trees into rangelands will be key for future biodiversity conservation.

Figure 8.

Trends in annual cover (A–C) and tree cover (D–F) across three management scales. A, Field office-scale trends in annual cover across the western Cold Deserts. B, Allotment-scale trends in annual cover for allotments in the Jarbidge Field Office in Idaho (blue outline in A). C, Pixel-scale trends in annual cover within a single allotment (blue outline in B). D, Field office–scale trends in tree cover for field offices within the Western Cold Deserts. E, Allotment-scale trends in tree cover for allotments in the Jarbidge Field Office in Idaho (no trends were calculated for allotments in gray, where tree cover is ≈ 0%). F, Pixel-scale trends in tree cover within a single allotment.

The big picture afforded by RS revealed the vulnerability of rangelands to continued tree expansion in the Northern Great Plains (see Fig. 5). Despite having < 5% tree cover today, expansion has occurred on two thirds of allotments in the Northern Great Plains during our 30-yr analysis window ( Table A13 (mmc1.pdf)). Increases in tree cover represent both infilling of existing stands and recruitment of new trees into formerly treeless rangelands (Filippelli et al. 2020). The latter process may be the most consequential for conservation of grasslands as trees rapidly mature and become seed sources for continued invasion. Using RS-derived maps to focus management on recently invaded rangelands can be an efficient strategy for preventing tree invasion into the most vulnerable habitats (Twidwell et al. 2021). While we cannot determine which species are responsible for tree cover increases from the RS data alone, the AIM data showed that junipers were the most common trees in all four ecoregions with significant increases in tree cover ( Table A14 (mmc1.pdf)). Where junipers have recently invaded rangelands, targeted removal can have benefits for local shrub and grassland associated species (Olsen et al. 2021; Roberts et al. 2022), increase soil moisture availability (Roundy et al. 2014), and help maintain herbaceous productivity (Morford et al. 2022).

Shrub cover was increasing in the Northern Great Plains and Mediterranean California and declining in the Western Cold Deserts (see Fig. 3). Shrub losses in the Western Cold Deserts likely reflect losses of sagebrush to increased wildfire and highlight the threat facing sagebrush-dependent wildlife (Coates et al. 2016). This region also saw both rapid increases in annuals and declines in bare ground (see Fig. 5), trends that are correlated with one another (see Fig. 7) and loss of shrubs (Fig. A1). Taken together, our results support the idea that exotic annual grasses are reducing bare ground, fueling more frequent and severe wildfires, and ultimately driving losses of cover of sagebrush and other shrubs (Chambers et al. 2014; Coates et al. 2016). Notably, we did not observe significant increases in shrubs in the AZ/NM Highlands and Warm Deserts, where previous studies raised concerns of shrub expansion (Eldridge et al. 2011). Satellite-based results do not disprove earlier cases of shrub expansion in the Southwest; however, they may indicate that the process of shrub invasion has stopped in recent decades or that it is a local phenomenon rather than ecoregion wide.

Quantifying spatial heterogeneity in vegetation trends

Heterogeneity and scale present grand challenges in monitoring rangelands across space and time (Fuhlendorf et al. 2017; Sayre 2017). Our findings show that trends varied greatly at all subecoregion scales: between field offices, allotments, and especially pixels within the same allotment (see Fig. 6). For rangeland managers, a critical implication of this result is that trends measured at a single field plot are unlikely to represent trends across an entire allotment—let alone a field office or ecoregion (West 2003; West and Wu 2003). Fine-scale spatial heterogeneity may be especially problematic for quantifying trends in tree, shrub, and perennial cover, all of which had the majority of trend variance at the pixel scale (see Fig. 6). On the other end of the spectrum, annual cover and production trends showed more variation at the field office and allotment scales than at the pixel scale (see Fig. 6). This may indicate the importance of larger regional-scale factors such as climate in determining the direction and rate of annual trends. In general, the large proportion of total trend variance at the 30-m pixel scale cautions against extrapolating plot-level trends to larger spatial scales. Failure to consider the bigger picture provided by RS may result in poorly targeted management actions that fail to add up to beneficial outcomes at larger scales. The wall-to-wall coverage provided by RS data allows today's rangeland scientists to pinpoint hotspots of vegetation change and develop a spatially targeted strategy for defending intact habitat from annual invasion and woodland expansion (Doherty et al. 2022; Maestas et al. 2022; Tarbox et al. 2022).

In this analysis, we summarize vegetation trends at the ecoregion scale; however, BLM managers will frequently want to examine trends at much finer spatial scales. In Figure 8, we illustrate how vegetation trends derived from RS could be visualized at the appropriate scale for management. For instance, trends aggregated at the scale of BLM field offices highlight subregions where invasions by annuals and expansion of trees are most rapid (see Fig. 8A, 8D). Allotment-scale trends within a single field office (Fig. 8B, 8E) help locate priority areas for reviewing land use authorizations and identify core areas of intact rangelands for conservation (Maestas et al. 2022). At an even finer scale, trends measured at the pixel scale (Fig. 8C, 8F) make it clear which areas within an allotment have experienced increases in annuals and trees and which areas have experienced no change or even declines. Despite their utility, the large size of RS datasets can make them difficult to work with in desktop GIS software (Tarbox et al. 2022). To facilitate adoption of RS data in rangeland management, stakeholders may wish to prioritize the development of online dashboards for displaying trends at field office, allotment, and pixel scales (see Fig. 8). Examples of web apps that display spatial data at a variety of spatial and temporal scales are the U.S. Forest Service's Rangeland Production Monitoring Service (US Forest Service, RMRS 2018) and the US Drought Monitor (University of Nebraska, Lincoln 2022).