Preserving prairies: understanding temporal and spatial patterns of invasive annual bromes in the Northern Great Plains

Abstract

Two Eurasian invasive annual brome grasses, cheatgrass (Bromus tectorum) and Japanese brome (Bromus japonicus), are well known for their impact in steppe ecosystems of the western United States where these grasses have altered fire regimes, reduced native plant diversity and abundance, and degraded wildlife habitat. Annual bromes are also abundant in the grasslands of the Northern Great Plains (NGP), but their impact and ecology are not as well studied. It is unclear whether the lessons learned from the steppe will translate to the mixed‐grass prairie where native plant species are adapted to frequent fires and grazing. Developing a successful annual brome management strategy for National Park Service units and other NGP grasslands requires better understanding of (1) the impact of annual bromes on grassland condition; (2) the dynamics of these species through space and time; and (3) the relative importance of environmental factors within and outside managers' control for these spatiotemporal dynamics. Here, we use vegetation monitoring data collected from 1998 to 2015 in 295 sites to relate spatiotemporal variability of annual brome grasses to grassland composition, weather, physical environmental characteristics, and ecological processes (grazing and fire). Concern about the impact of these species in NGP grasslands is warranted, as we found a decline in native species richness with increasing annual brome cover. Annual brome cover generally increased over the time of monitoring but also displayed a 3‐ to 5‐yr cycle of reduction and resurgence. Relative cover of annual bromes in the monitored areas was best predicted by park unit, weather, extant plant community, slope grade, soil composition, and fire history. We found no evidence that grazing reduced annual brome cover, but this may be due to the relatively low grazing pressure in our study. By understanding the consequences and patterns of annual brome invasion, we will be better able to preserve and restore these grassland landscapes for future generations.

Introduction

During the last century, much of the prairie within the Northern Great Plains (NGP) has been plowed for cropland, planted with non‐native grasses to maximize livestock production, or otherwise developed, making it one of the most threatened ecosystems in the United States. Two decades ago, it was estimated that > 75% of the area of native mixed‐grass prairie had been lost since European settlement (Samson and Knopf ), and recent high soybean and corn prices have accelerated the rates of loss (Wright and Wimberly ). Native prairies remaining within the NGP are threatened by the invasion of exotic plant species. While many invasive plants threaten NGP prairies (e.g., Larson et al. ), two annual brome species—cheatgrass (Bromus tectorum L.) and Japanese brome (Bromus japonicus Thunb.)—have great potential to impair NGP native grasslands because they are abundant and widespread (Ogle et al. ).

Cheatgrass and Japanese brome are Eurasian, annual grasses that have been a part of the NGP landscape for more than a century, but their invasion in the region has accelerated since 1950 (Schachner et al. ). It is likely that cheatgrass spread into the NGP from other areas of the United States, but there may have also been direct introductions from European immigrants (Schachner et al. ). The history of Japanese brome invasion in the NGP is less well documented, but there are records of its presence in the early part of the 20th century (e.g., Wright and Wright ). The two species are functionally similar; both are cool‐season grasses that germinate in the fall, winter, or early spring and begin growth earlier in the spring than many of the native perennials. This gives them a head start in competing for moisture and nutrients with native species, thereby reducing native production and nutritional content of forage (Haferkamp et al. , ). Annual bromes' barbed seeds are a nuisance to grazing wildlife, and the species' early senescence (compared with native grasses) adversely affects wildlife forage availability later in the growing season (Ogle et al. ). In North America's western steppe ecosystems, where estimated presettlement fire return intervals averaged 196 yr (Balch et al. ), cheatgrass is well known for greatly increasing the frequency of fires, leading to widespread replacement of fire‐intolerant native sagebrush–bunchgrass steppe with fire‐prone annual grasslands (D'Antonio and Vitousek , Knapp , Balch et al. ). In contrast, the estimated fire return interval in native mixed‐grass prairies in the NGP is 9–12 yr (Guyette et al. ), and the lack of fire is often considered a disturbance. The effects of cheatgrass and Japanese brome invasion are not as well studied in the NGP as in steppe or other western U.S. ecosystems (Brooks et al. ). More information on the dynamics of annual brome invasion in the NGP is needed to protect and preserve the remaining native prairies.

The National Park Service (NPS) manages 11 park units within the NGP totaling 168,750 ha (Fig. ). While the NPS mission to “preserve ecological integrity and cultural and historical authenticity” (NPS ) includes protecting and restoring native prairie, land managers struggle to accomplish this because of the continual pressure of exotic invasive species, such as annual bromes, and fundamental changes in climate, fire, and grazing disturbance regimes that historically maintained native prairies. Successful management requires integrating science into decision‐making. This often entails translating results from small‐scale experiments into landscape‐scale responses, but this is complicated by the lack of relevant experiments, the complexity of ecological systems, high climate variability and its effects on ecosystem response, and logistical and financial constraints. Consequently, there is a lack of proven methods for restoration and long‐term invasive species control in the NGP.

Managing annual bromes in the NGP has proved particularly challenging for a number of reasons. First, annual bromes are already very abundant (Fig. ; USDA‐NRCS ). Second, annual brome populations respond strongly and quickly to climate fluctuations (Mack and Pyke ). Thus, the NGP's naturally high interannual climate variability (Borchert ) often causes large changes in annual brome abundance that are outside of a manager's control. Third, while a variety of management actions, including highly controlled grazing, herbicide application, and targeted prescribed fires, have shown promising short‐term results for controlling annual bromes in research‐scale plots (e.g., Whisenant and Uresk , Haferkamp et al. , Harmoney ), there is a lack of proven methods for long‐term, landscape‐scale annual brome control. Finally, land‐use history, topoedaphic factors, and recent grazing and fire history are thought to influence annual brome abundance across the landscape (Hulbert , Gundale et al. , Lovtang and Riegel , Reisner et al. , Bansal and Sheley ). This large variability and context dependency make it more difficult to determine the best approach to their control across large landscapes. Better understanding of the behavior of annual bromes in NGP prairies through space and time is critical for developing a successful management approach and preserving the ecological integrity of prairies within NPS units.

In this study, we explore long‐term monitoring data with the goal of guiding research and adaptive management of annual bromes in native prairies, particularly in the unique management conditions of NPS units. We use a monitoring data set from 1998 to 2015 to explore the dynamics of invasive annual brome species in NPS units in the NGP. We ask three questions. First, what are the temporal and spatial trends in annual brome abundance from 1998 to 2015? Our hypothesis, based on casual field observations, is that annual bromes have been increasing in abundance within the NPS units despite large variations in abundance across space and time. Second, how is annual brome abundance related to grassland composition and native species richness? As has been shown in steppe ecosystems, we expect native species richness to be lower in invaded areas. Furthermore, more diverse prairies comprised of species consistent with an intact disturbance regime are expected to be more resistant to invasion. Finally, what are the relative impacts of climate, soils, physical environmental factors, fire, and grazing on brome abundance in NGP grasslands? We expect that patterns of annual brome abundance are complex within the NGP and are driven by many factors (Fig. ) including land‐use history, the native vegetation, nutrient availability, and fire frequency. We expect climate to be a strong driver of abundance and that years with wet autumns (when seeds are germinating) and springs (when bromes grow and set seed) increase abundance. Unlike in steppe ecosystems, we expect that fire will reduce brome abundance because it reduces the litter layer and consequently annual brome germination.

Methods

Study area

Data were collected from seven national park units in the NGP varying in size from ~337 to 40,000 ha (Fig. ): Fort Laramie National Historic Site (NHS) and Devils Tower National Monument (NM) in Wyoming; Scotts Bluff NM and Agate Fossil Beds NM in Nebraska; and Wind Cave National Park (NP), Jewel Cave NM, and Badlands NP in South Dakota. The region is dominated by short‐grass and northern mixed‐grass prairie, which grades into ponderosa pine forest in the Black Hills. Western wheatgrass (Pascopyrum smithii (Rydb.) Á. Löve), blue grama (Bouteloua gracilis (Willd. ex Kuth) Lag. Ex Griffiths), and needle and thread (Hesperostipa comata (Trin. & Rupr.) Barkworth) are the most common native grasses, but annual bromes and Kentucky bluegrass (Poa pratensis L.) are also common. In addition to prairie, some parks have ponderosa pine and grassland savannah (Devils Tower NM, Wind Cave NP, Jewel Cave NM) and badland plant communities (Scotts Bluff NM and Badlands NP). Despite its close proximity to the other parks in this study, we did not include monitoring data from Mount Rushmore National Memorial (Fig. ). Cheatgrass generally does not perform well in shade (Pierson et al. ) and annual bromes are rarely found within the heavily forested areas that cover Mount Rushmore.

The NGP has a continental climate, with hot summers and cold winters. At the seven park units included in this study, two‐thirds of annual precipitation falls in April–August, and mean annual precipitation ranges from ~40 to 50 cm (NOAA ). However, extreme variation in temperature and precipitation is typical in the region. For example, during the last 20 yr at Wind Cave NP, total annual precipitation has ranged between 33 and 73 cm. The native vegetation is adapted to this variation, and productivity responds strongly to decreases in summer precipitation (Yang et al. , Smart et al. ). Species richness and diversity in regional grasslands are also sensitive to temperature and precipitation fluctuations, but the response is complex and less predictable (Jonas et al. ).

In addition to climate, grazing and fire are drivers of vegetation patterns in these parks (e.g., Bachelet et al. ). Major grazers are American bison (Bison bison) and black‐tailed prairie dog (Cynomys ludovicianus) at Badlands and Wind Cave NPs, elk (Cervus elaphus) at Wind Cave NP, and domestic horses (Equus caballus) at Fort Laramie NHS. Pronghorn (Antilocapra americana), bighorn sheep (Ovis canadensis), mule deer (Odocoileus hemionus), and white‐tailed deer (Odocoileus virginianus) occur in the region but are generally much less abundant than the major grazers. Management control over the timing and intensity of grazing in park units is low compared with that in surrounding properties, as wildlife are free‐roaming within (and sometimes across) park boundaries.

Historically, fire was a common disturbance in NGP grasslands, with natural fire return intervals of 9–12 yr (Guyette et al. ). Natural fires have been suppressed for most of the last century, but the use of prescribed burning to mitigate the effects of the absence of natural fires has increased over time in these parks since its start at Wind Cave NP in 1973 (Wienk et al. ). As of 2015, there is a mosaic of recently burned and unburned areas in most of the parks in our analyses, but the history and extent of fire ranges widely among the parks. In general, however, prescribed fire is used more frequently in park units than in surrounding properties (Symstad and Leis, in press).

Vegetation, physical environment, and small‐scale disturbance data

We compiled vegetation composition data from 430 long‐term monitoring plots between 1998 and 2015 (inclusive), following methods established by NPS () and Symstad et al. (, b). The plots are maintained by two NP Service programs: the NGP Fire Effects program (NG‐FEP) and the NGP Inventory & Monitoring program (NGPN). The NG‐FEP established and began monitoring 136 plots from 1998 to 2010. During this time, data collection followed NPS national fire ecology program protocols (NPS ): In grassland sites, herb‐layer (≤ 2 m) vegetation cover and height data were collected using a point‐intercept method, with 100 points evenly distributed along a single 30‐m transect, and height measured at each point as the distance of the top‐most intercept from ground level. In forested sites, plots were 0.1 ha (20 × 50 m) in size and point‐intercept data (166 points) were collected along one of the two 50‐m sides. For each live tree with a diameter at breast height (dbh) > 15 cm located within the 0.1‐ha plot, the species and dbh were recorded. The densities of smaller trees (2.54 cm ≤ dbh ≤ 15 cm) were measured within a subset of the plot area.

Prior to 2010, NG‐FEP plot locations were located randomly within major vegetation types within areas planned for prescribed burning (burn units) in the near future. In 2011, the NG‐FEP adopted the data collection protocols of the NGPN (next paragraph), which was developed in conjunction with the NG‐FEP, but continued to establish new plots only within new burn units. Thus, the number and distribution of NG‐FEP plots varies widely among parks.

The NGPN plots are 0.1 ha (20 × 50 m) in size and are randomly located within a park, and the number of plots is approximately proportional to the park's size. Plot locations were chosen based on a spatially balanced probability design (Generalized Random Tessellation Stratified [GRTS], Stevens and Olsen , ). The GRTS design supports inference from plots to the entire park and is spatially balanced. Moreover, the annual subset of plots visited is also spatially balanced. Data on ground cover (e.g., litter, bare soil, rock), herb‐layer (≤ 2 m) height, and plant cover by species were collected on two parallel 50‐m transects using a point‐intercept method. Ground cover type was collected at every point, whether or not the point had plant cover. Tree density was measured as described above. The type and approximate area of common disturbances, including rodent mounds, animal trails, grazing, and fire, within a 27‐m radius circle centered in each plot was recorded. The 27‐m radius encompassed the grid cell used for the GRTS sampling regime. The size of each disturbance (m²) ranged from 0 (not present) to 2290 (the whole plot area). Total disturbance was calculated as the sum of all individual disturbances, so the value can be > 2290 m². More detailed methods can be found in the NGPN Vegetation Monitoring Protocol (Symstad et al. , b). A total of 294 plots from 2011 to 2015 were surveyed by the NG‐FEP and NGPN using these methods. At NG‐FEP and NGPN plots, slope grade and aspect of both the transect on which data were collected and the hill on which the plot occurred were measured, and Universal Transverse Mercator (UTM) coordinates and elevation of the center point were recorded in the field using a GPS. Most of the 430 plots were visited more than once between 1998 and 2015, yielding a total of 1162 separate visits, but revisit schedules and frequencies varied. On average, plots were revisited 2.8 times. Because plots were randomly located within parks or burn units and were not established with any a priori knowledge of brome abundance, and because most plots were only sampled two to three times, we have no reason to expect a spatial bias in the data over time with respect to annual bromes.

The absolute herbaceous cover of plants was derived from the total number of intercepts along a transect. This value could total more than 100, because multiple layers of vegetation were recorded at each location along the transect, or less than 100 when there were large amounts of bare ground. Absolute herbaceous cover was used as a potential explanatory variable because it serves as a proxy for site productivity. Our response variable in analyses was the relative cover of annual bromes, which was calculated by dividing the total number of cheatgrass and Japanese brome intercepts by the total number of vegetation intercepts on the transect(s) in a plot. We choose to use relative, instead of absolute, cover because it represents the relative influence of the annual bromes (or other plant groups) in the plant community regardless of fluctuations in weather or differences in soil that influence total plant cover. Results were similar using absolute cover (data not shown). Plant life forms, nativity, and family were based on definitions from the USDA Plants Database (USDA‐NRCS ). Plant species richness was calculated for each plot using the total number of species intercepted along the transects. Average height was calculated over all points on the transects. Relative cover of native forbs and sedges (all Carex spp.) was calculated for each plot. Along with other metrics described above, these were used as a proxy to assess grassland condition, because forbs contribute positively to species richness and some evidence from the NGP suggests that Carex cover is sensitive to disturbance in these systems (Wienk et al. ).

Data for explanatory factors

We built an initial conceptual model for annual brome abundance in NGP parks based on available monitoring data and factors thought to affect annual brome such as climate (Chambers et al. , Bradley ), soils (Belnap et al. , Bansal and Sheley ), and grazing (Mack and Pyke , Harmoney ). This yielded 34 potential explanatory variables for brome abundance (Fig. ). We extracted total monthly precipitation, mean minimum temperature, and mean maximum temperature for 1997–2014 by plot location using the 800 m scale PRISM data (PRISM Climate Group, Oregon State University, http://prism.oregonstate.edu). Seasonal total precipitation and seasonal average minimum and maximum temperatures calculated from the monthly data were used as climatic explanatory variables. Soil composition (percent sand, silt, and clay), pH, and cation‐exchange capacity were determined by a spatial join of soil maps for each park (NRCS ) and the center of all plot locations using ArcGIS software 10.2.1 (ESRI, Redlands, California, USA). Fire history maps were compiled for each park and cross‐referenced with plot locations. For each plot visit, we determined the number of years since it burned and the number of recorded fires. The length of the fire history record varied by park, but most began in the 1980s. For plots where no burns were recorded, we calculated the difference between the year of data collection and the oldest fire recorded in the park. This is likely an underestimate of the true time since it burned because fires were infrequent prior to the 1980s. Annual nitrogen (N) deposition values were based on data from the National Atmospheric Deposition Program's total N deposition maps (http://nadp.sws.uiuc.edu/committees/tdep/tdepmaps). Grids were downloaded and overlaid with park boundaries to get one value per park per year for 2000–2013. The total number of bison culled during roundups from 1997 to 2015 for Wind Cave NP and Badlands NP was provided by the NPS. Potential explanatory variables derived from monitoring data (see Vegetation, physical environment, and small‐scale disturbance data) were tree basal area, herbaceous cover, native species richness, vegetation height, ground‐level litter cover and bare cover, transect slope grade, hill slope grade, slope aspect, and elevation. Aspect was adjusted based on equations in McCune and Keon () to account for the similarity between 360 and 0 degrees.

Statistical analysis

We included only plots with tree basal area < 15 m²/ha (e.g., grasslands or open savannah; this excluded 53 forested plots) and with annual bromes present (> 0% cover) during at least one visit in 1998–2015 in our analyses (295 of the 430 monitored plots and 823 plot visits). Some visits with no recorded annual brome cover remained in the data set because of the large interannual variation in brome abundance in our analyses. We chose not to include plots where annual bromes had never been recorded because a large number of zero values for the response variable in the data set make models mathematically difficult. Moreover, the lack of annual bromes in plots where they were never recorded may occur because of a lack of propagule pressure. Our data set is not well suited for determining propagule pressure, and therefore, we focused on understanding spatio‐temporal variability in areas where we knew the species occurred. However, analyses including plots in which annual bromes were never recorded yielded similar results to those without these plots (data not shown).

All statistical analyses were completed in R version 3.2.2 (2015‐08‐14, R Foundation for Statistical Computing, Vienna, Austria). To address our first question, we used linear mixed models to determine trend over time in relative percent cover of annual bromes across all parks; year, park, and plot were included as nested random factors. To address our second question, we used linear mixed models to test for the relationship between brome relative cover and seven measures of grassland condition: native species richness, native forb cover, sedge cover, bare ground cover, litter cover, total herbaceous plant cover, and average canopy height. These models were used to explore the patterns in our observational data set rather than to infer causation. As in the above models, year, park, and plot were included as nested random factors. Relative brome abundance was log‐transformed prior to analysis to meet assumptions of heteroscedasticity and normality. Best‐fit models were chosen based on Akaike's information criterion (AIC) and the most parsimonious model with the lowest AIC score was chosen when the difference between models was not significant (P > 0.05). We report the conditional R², which describes the proportion of the variance explained by both the fixed and random factors. For our third question, our goal was to explore the relationship between multiple factors and brome abundance. However, the final data set included one response variable (relative cover of brome at 823 site visits), 34 potential explanatory variables, and missing information for some of these variables for some site visits. Grazing disturbance data were commonly missing, as were recent estimates of N deposition (2014–2015) and climate parameters (2015) and ground cover information (e.g., litter cover) in pre‐2010 plots. Thus, we used two analytical approaches for this question. In the first, we used regression and conditional inference trees from ctree in the party package (Hothorn et al. ) for exploratory analyses using all site visits and all explanatory variables. In the second, we evaluated the fit of linear mixed models of relative brome abundance (lme in nlme package; Pinheiro et al. ) and the explanatory variables using a step function (stepAIC) to test for the best‐fit model based on AIC. The stepAIC function starts with the largest model and uses a stepwise algorithm to add and remove all variables and calculates AIC for each resulting model. Analysis of variance yields the final model with the lowest AIC. We used these two approaches because they have different strengths. Using the regression trees, we were able to include the full sample size and a large number of explanatory factors. The mixed model approach was limited in size and number of potential factors, but was able to account for the repeated measures at plots.

In the second approach, the initial regression model included 19 of the 34 explanatory variables as fixed factors: park, litter cover, herbaceous cover, native species richness, tree basal area, slope aspect, hill slope grade, transect slope grade, years since fire, number of fires, soil composition (% silt and clay), and seven climate variables: fall, winter, spring, and summer precipitation and average maximum temperature in fall, winter, and spring (see Appendix S1 for mean and range of each factor). Plot was included as a random factor nested in year and park. Potential explanatory variables in this approach were reduced from the 34 shown in Fig.  to 19 for two reasons. First, this model evaluation method requires consistent sample sizes across models, so we eliminated variables with mostly missing information to maintain the greatest sample size possible. For example, evidence of grazing was assessed in few plots and grazing evidence was rare where it was assessed. This refinement yielded 418 site visits for which all information was available (e.g., no missing cells). Second, where variables were highly correlated (R² > 0.5; e.g., vegetation height and herbaceous canopy cover or minimum and maximum temperatures), we retained only the variables that appeared in the regression tree analysis and/or resulted in the model with the lowest AIC. We also used variance inflation factors (vif.lme in R) to help detect multicollinearity and eliminated factors causing it (e.g. annual N deposition). As park was consistently an important factor, we repeated the model building procedures using the groupings of parks that appeared in the regression tree analysis: short‐grass prairie parks (Scotts Bluff NM and Fort Laramie NHS), northern mixed‐grass parks (Agate Fossil Beds NM and Badlands NP), and Black Hills parks (Wind Cave NP, Jewel Cave NM, and Devils Tower NM).

Results

Trends in annual brome abundance

Over the 18‐yr monitoring period, the average relative percent cover of annual bromes throughout NGP parks varied significantly across years (F_1,16 = 17.8, P < 0.001), ranging from < 1% in 1998 to 27% ± 3% cover in 2010, and increased significantly over the monitoring period (Fig. ; conditional R² = 0.93). This increase was not steady, however, in that lows occurred approximately every three years. This irregular cycle did not follow precipitation in the region, although a severe drought in 2012 likely contributed to relatively low cover in that year and the following year (Fig. ). Fort Laramie NHS has high annual brome cover and was only visited in the last few years of data collection (2011–2015); because of the relatively few visits in this park, the increasing trend was still significant when Fort Laramie NHS was excluded from the analysis (F_1,16 = 16.3, P < 0.001).

Furthermore, relative brome cover varied greatly among the seven parks (F_6,68 = 15.3, P < 0.001; Fig. a), but the increase over time was similar in all the parks (time × park interaction F_6,68 = 0.9, P = 0.505). When parks were analyzed in separate models, three parks showed a significant increase in annual brome cover over time (Fig. b): Scotts Bluff NM (F_1,15 = 19.1, P < 0.001, conditional R² = 0.91), Wind Cave NP (F_1,15 = 9.5, P = 0.007, conditional R² = 0.90), and Agate Fossil Beds NM (F_1,7 = 8.0, P = 0.025, conditional R² = 0.92).

Relationship between brome abundance and grassland condition

There was a negative relationship between the relative cover of annual bromes and native species richness (F_1,551 = 36.3, P < 0.0001) (Fig. ). Other metrics of condition such as native forb cover (F_1,551 = 20.3, P < 0.001) and bare ground cover (F_1,428 = 10.8, P = 0.001) also showed a negative relationship with relative cover of annual bromes. Conversely, as average height (F_1,551 = 5.7, P = 0.017) and litter cover (F_1,428 = 17.7, P < 0.001) increased, there was an increase in relative brome abundance. There was no significant relationship between Carex cover and annual brome cover (F_1,551 = 2.3, P = 0.1303).

Regression trees and linear models describing brome abundance

In the regression tree approach, we found that the factors related to annual brome abundance varied significantly among parks. Short‐grass prairie parks (Fort Laramie NHS and Scotts Bluff NM) differed from the rest in that they had litter cover as a significant node (Fig. ). At these parks, a less continuous litter layer (≤ 80% of ground cover) was associated with lower cover of annual bromes (P < 0.001). Annual brome cover was very high (> 60% cover) in sites with high (> 80%) litter cover and in plots with low native species richness (≤ 4 native species per transect).

Badlands NP (where all plots were < 930 m in elevation) and the few plots in Wind Cave NP that were < 1107 m in elevation shared similar patterns of annual brome abundance (Fig. ). Annual brome abundance was lowest when the average maximum winter temperatures were > 5.2°C (1.8°C warmer than average) and herbaceous canopy cover was moderate or lower (≤ 140%, average herbaceous cover across all sites was 135% [Appendix S1: Table S1]). When average winter maximum temperature was ≤ 3.4°C, annual brome cover tended to be low in level areas (slope grade < 3%) with low‐to‐moderate herbaceous cover (≤ 132%).

In the Black Hills parks (Devils Tower NM, Jewel Cave NM, and Wind Cave NP), patterns of annual brome abundance were best predicted by herbaceous cover (Fig. ). The lowest annual brome cover in these parks occurred in areas of low‐to‐average herbaceous cover. In Agate Fossil Beds NM, steeper transect slopes (> 2%) tended to have higher annual brome cover.

The best‐fit linear regression model for all parks combined included 11 explanatory factors (Table ; conditional R² = 0.94). The final model included park; herbaceous cover; native species richness; litter cover; fall and winter precipitation; soil % silt and % clay content; transect slope grade; hill slope grade; and years since fire. Park, herbaceous cover, native species richness, and hill slope grade were more significant than the other factors (Table , P < 0.001). Higher herbaceous and litter cover, lower native richness, wetter winter and fall, lower silt and clay content, steeper slopes, and more years since fire were associated with higher relative cover of annual bromes (Table ). In the short‐grass prairie parks, higher relative cover of annual bromes was most associated with high herbaceous and litter cover, wetter springs, and sites with lower soil silt content and southern aspects. Park was not a significant factor in this analysis or the other ecoregions, suggesting that parks in the same ecoregions share similar patterns of brome abundance. In northern mixed‐grass prairie parks, higher herbaceous cover, lower native species richness, wetter winter, steeper transect slope grade, and less frequent fires were the most significant factors associated with higher annual brome relative cover. In Black Hills parks, where annual brome cover was relatively low, the best‐fit model included herbaceous cover, native species richness, winter and spring maximum temperature, fall, winter, and spring precipitation, hill slope grade, and aspect (Table ). The most significant factors were herbaceous cover and spring precipitation. The direction of the relationships within the park groups generally followed those seen in the larger model. However, there was some variation between short‐grass prairie and Black Hills parks. In the short‐grass prairie, warmer and wetter springs were associated with higher brome cover. In the Black Hills parks, we saw the opposite pattern where cooler and drier springs were associated with increases in annual brome cover.

Explanatory factors from the best‐fit linear mixed models predicting annual brome relative cover in Northern Great Plains grasslands

Factor	Northern Great Plains	Short‐grass prairie	Northern mixed‐grass prairie	Black Hills
Slope	df	F	Slope	df	F	Slope	df	F	Slope	df	F
Park		6, 31	15.1***	…	…	…	…	…	…	…	…	…
Plant community
Herbaceous cover	+	1, 363	45.5***	+	1, 65	31.7***	+	1, 187	37.4***	+	1, 92	23.0***
Native richness	−	1, 363	13.7***	−	1, 65	2.6	−	1, 187	11.6***	−	1,92	0.7
Litter cover	+	1, 363	4.5*	+	1, 65	21.3***	…	…	…	…	…	…
Tree basal area	…	…	…	NA	NA	NA	NA	NA	NA	…	…	…
Temperature
Fall maximum	…	…	…	+	1, 65	1.2	…	…	…	…	…	…
Winter maximum	…	…	…	…	…	…	…	…	…	+	1, 92	4.2*
Spring maximum	…	…	…	+	1, 65	1.1	…	…	…	−	1, 92	6.8*
Precipitation
Fall	+	1, 363	3.0****	+	1, 65	0.2	…	…	…	+	1, 92	1.7
Winter	+	1, 363	4.7*	…	…	…	+	1, 187	10.0**	+	1, 92	0.1
Spring	…	…	…	+	1, 65	8.6**	…	…	…	−	1, 92	7.1**
Summer	…	…	…	−	1, 65	< 0.1	…	…	…	…	…	…
Soil texture
% Silt	−	1, 363	7.5**	−	1, 65	15.9***	…	…	–	…	…	…
% Clay	−	1, 363	5.2*	…	…	…	…	…	…	…	…	…
Physical environment
Transect slope grade	+	1, 363	9.6**	−	1, 65	1.1	+	1, 187	30.6***			–
Hill slope grade	+	1, 363	15.7***	+	1, 65	5.9*	…	…	…	+	1, 92	6.8*
Slope aspect	…	…	…	+	1, 65	9.7**	…	…	…	−	1, 92	2.1
Elevation	NA	NA	NA	−	1, 65	6.9*				…	…	…
Fire
Years since fire	+	1, 363	6.6*	+	1, 65	2.5	…	…	…	…	…	…
No. fires	…	…	…	+	1, 65	3.1****	−	1, 187	11.6***	…	…	…
Model fit
Conditional R2			0.94			0.99			0.99			0.99
Akaike's information criterion (AIC)			1216			255			579			342
Plot visits			418			91			206			121

Notes

Models were selected using AICs. Fixed effects with positive slopes are indicated by (+) and effects with negative slopes are indicated by (−). F values, degrees of freedom, and significance (***P < 0.001, **P < 0.01, *P < 0.05, ****P < 0.1) are provided for all explanatory factors. Factors with “NA” were not included in the model. The model selection process included some factors that were not individually significant in the best‐fit model. The Northern Great Plains data included all park units; the short‐grass prairie data include Scotts Bluff National Monument (NM) and Fort Laramie National Historic Site; northern mixed‐grass prairie data include Agate Fossil Beds NM, Badlands National Park (NP); and the Black Hills data include Wind Cave NP, Devils Tower NM, and Jewel Cave NM.

Discussion

Implications of increasing annual brome abundance for NGP grasslands

Over the last 18 yr, annual bromes comprised up to 30% of herbaceous cover in prairies of seven national park units in the NGP (Figs.  and ). Across these same parks, this degree of annual brome invasion is associated with an approximately 20% decrease in native species richness (Fig. ). This is consistent with the large body of literature and experiments that have shown a negative effect of annual bromes on native ecosystems in this region and elsewhere (e.g., Melgoza et al. , Knapp ). Because our data set is observational, rather than experimental, it is impossible to differentiate between cause and effect in this relationship, and responses could go both ways. Annual bromes may preferentially establish in areas of low native species richness and, conversely, negatively affect the recruitment and persistence of native species. In competition experiments, Japanese brome has been shown to reduce tiller and leaf growth of native perennial grass seedlings found in NGP grasslands (Romo and Eddleman ). Numerous studies have demonstrated that cheatgrass can be an aggressive competitor with perennial grasses (Harris , Melgoza et al. , Lucrecia and Johnson ). While our data show patterns rather than mechanisms, annual bromes are clearly a concern for land managers in the NGP trying to maintain native prairie species richness and diversity. Thus, better understanding of the behavior of annual bromes in the NGP through space and time is critical for protecting and preserving the remaining prairies in this region.

The spatial and temporal distributions of monitoring plots in our data set are not ideal for detecting trends over the 1998–2015 period. For example, the high spatial variability of annual brome abundance in the data set (Fig. ) might suggest that low brome abundance in 1998 is due to low sample size in that year (15 sites). However, the increasing trend over all parks is still significant when that year is removed from the analysis (F_1,15 = 12.6, P = 0.003). In addition, low sample size cannot explain the drop in cover seen in over 100 sites in 2012. The large temporal and spatial variation exhibited in our data set emphasizes that land managers should be cautious when using 1‐ to 2‐yr data sets to assess the condition and extent of brome invasions. Rather, it is clear that continued long‐term monitoring is critical to better understand the extent and trends of brome invasions in the NGP.

Although not ideal, our data set is one of the few that quantifies the trend in annual brome abundance at a regional scale, and it reveals a more than doubling of the influence of these species on the grasslands of national park units over the past 18 yr (Fig. ), supporting general observations of accelerating invasion in this region since the 1950s (Schachner et al. ). While the overall trend was significant across all parks, park‐specific analyses show that not all parks experienced a significant increase over the sampling period. Two parks with no trend and the lowest degree of invasion (annual bromes < 5% relative cover), Devils Tower NM and Jewel Cave NM, may have experienced less propagule pressure than other parks due to their locations in generally more forested areas. The other two parks with no trend (Badlands NP and Fort Laramie NHS) have a substantially greater degree of invasion. Although the lack of trend at Fort Laramie NHS may be due to its short monitoring period (5 yr), both of these parks have historical records of agricultural use. Determining the land‐use history of any given point on the landscape is time‐consuming and not always possible due to lack of records, but this type of information may shed light on the variation in annual brome trends and the factors that influence them among parks.

On a regional scale, relative cover of annual bromes appears to follow a 3‐ to 5‐yr cycle (Fig. ). Although a large decline in their abundance corresponded with a severe drought in 2012, our analyses show that other factors in addition to precipitation contribute to annual brome behavior. More information is needed to determine whether the cyclic pattern is driven by biological factors, such as a lag effect of climate or productivity, or by management cycles, such as prescribed fire implementation. A complex interaction of temperature, precipitation, and ecological factors determining annual brome emergence and growth, as demonstrated in the Great Basin (Mack and Pyke , Chambers et al. , Bansal and Sheley ), is also evident in the NGP parks in our analyses.

Both the regression tree and linear model approaches showed that relative cover of annual bromes is most associated with park unit, weather, the extant plant community, soil composition, and slope grade and aspect (Table , Fig. ). While our initial hope was to find a general model that could predict annual brome abundance across the NGP, ecological complexity and distinct land‐use histories made this difficult. Park units did, however, group together by ecoregion, and annual brome behavior in Black Hills parks, short‐grass prairie parks, and northern mixed‐grass prairie units was similar to one another within those regions (Fig. ). These groupings may be due to similar vegetation, climate, or other factors we did not explore. Despite differences across the ecoregions, the direction of the relationships between brome abundance and environmental factors other than weather were fairly consistent (e.g., herbaceous cover was positively associated with annual brome cover in all ecoregions). However, the variable response of annual brome cover to seasonal temperature and precipitation indicates that as land managers build programs to manage annual bromes across the region, it will be important to recognize the potential differences and similarities among park units.

Native plant communities resist invasion?

In addition to the negative effect of annual bromes on native species richness described in other studies cited in the previous section, the negative relationship between the two (Fig. ) is consistent with the hypothesis that intact native plant communities can resist invasion of annual bromes. This hypothesis is also supported by the fact that, when we used native herbaceous canopy cover instead of total herbaceous cover as a potential explanatory variable (data not shown), it replaced native transect richness as a significantly positive factor in the regression models. While our data cannot elucidate cause and effect, it is consistent with experimental evidence from a ponderosa pine–bunchgrass ecosystem in Arizona that suggests that established native grasses can resist invasion from cheatgrass even when propagule pressure is high (McGlone et al. ). However, land managers must be cautious in relying solely on native plant resistance because despite the overall trend, there were numerous sites that had high relative cover of annual bromes (> 20%) despite high native species richness (Fig. ).

One of the most consistent results across all park units and models was that two other aspects of the extant plant community, high total herbaceous cover and high litter cover, were positively associated with relative cover of annual bromes (Fig. , Table ). Although high absolute cover of annual bromes contributes to a high total cover of herbaceous plants, there is no clear causative reason for high relative cover of annual bromes to yield higher total cover. In contrast, high herbaceous cover is likely correlated with sites that have high available moisture and nutrients, which favor brome establishment. A thick litter layer can promote the establishment and spread of annual bromes, because it can provide a matrix to catch brome seeds and store moisture, allowing them to germinate and grow more effectively (Whisenant and Uresk , Bansal et al. , Jones et al. ). Litter and brome abundance may create a positive feedback cycle because of bromes' relatively low litter quality, which may slow decomposition rates (Evans et al. , Ogle et al. ). Thus, management strategies that reduce the depth and extent of litter may help control annual bromes, but the effectiveness may vary with the timing of that reduction, as well as climate (Bansal et al. ). Moreover, our results suggest that a reduction in the total herbaceous cover may reduce annual brome establishment. Grazing and fire are two management tools that strongly influence litter and herbaceous cover in NGP grasslands.

Grazing and fire

Surprisingly, large grazers and small mammal disturbance (e.g., prairie dogs) did not explain any variation in relative brome cover in our data set (Fig. ). We ascribe this to two factors. First, the density of grazers in park units, even those units with bison herds, is generally lower than most NGP lands grazed by cattle. For example, Badland NP manages for a herd of 700 bison, which corresponds to a low consumption rate of approximately 12% of plant productivity (Licht ). Second, the grazers in NGP national park units are free‐roaming, so the timing and intensity of their grazing are highly variable from one location to another. Thus, the grazing timing and intensity shown to reduce annual brome cover in experimental work (Harmoney ) are probably not occurring over a large portion of national park units in the NGP. Second, prairie dogs are also significant grazers and are present in Wind Cave NP, Scotts Bluff NM, Badlands NP, and Devils Tower NM, but due to the random allocation of NGPN plots within parks and avoidance of prairie dog towns for NG‐FEP plots, they were only present in a few plots within our data set. Continuous clipping of vegetation by prairie dogs may provide little opportunity for annual brome species to produce seeds, or the large amounts of bare ground and associated soil desiccation in prairie dog towns may not be conducive to seed germination and seedling survival, explaining the relatively low cover of annual bromes occurring in prairie dog towns in this region (Fahnestock and Detling ). On the other hand, areas with pocket gopher activity (i.e., mounds of loose, bare earth surrounded by undisturbed vegetation) seem to have high brome cover in some NGP parks (A. Symstad, unpublished manuscript). More directed experiments are needed to better understand the relationship between grazing and annual brome abundance in prairie systems. Likewise, as more monitoring data about these disturbances are collected, clearer patterns of their effects on annual brome abundance may emerge.

Our initial hypothesis was that fire was a driver of annual brome abundance in NGP grasslands because it removes the litter layer that promotes brome establishment and spread. Furthermore, spring or fall burns may burn seedlings when native warm‐season perennial grasses are dormant. This strategy notably contrasts with that in sagebrush habitats, where cheatgrass promotes a reduced fire return interval and causes declines in native shrubs (D'Antonio and Vitousek , Balch et al. ). Our data set did not show strong support for our hypothesis. The time since an area was burned was a relatively weak predictor of brome abundance and the number of fires in a site was only significant in northern mixed‐grass prairie parks (Table ). However, the direction was as expected, and there tended to be an increase in annual brome abundance as the time since a plot had burned increased and a decrease with more frequent burns. Future strategies that incorporate prescribed burning may be beneficial to restoring native prairies in the NGP, but more research is needed to understand how the timing, frequency, and severity of fires influence annual brome abundance.

The relationship between precipitation, temperature, and physical environment and brome abundance

As winter annuals, cheatgrass and Japanese brome are known to respond strongly to precipitation and temperature (Mack and Pyke , Melgoza et al. , Chambers et al. ). From the linear mixed model, we found that seasonal precipitation was an important factor driving annual brome abundance, but the direction of the response varied across ecoregions. At Scotts Bluff NM and Fort Laramie NHS, warmer and wetter springs and falls were associated with increased annual brome abundance. This pattern is consistent with a nearby experimental study showing a strong positive response of cheatgrass to warming (Blumenthal et al. ). It is also similar to the pattern seen in Washington, where a dry autumn can cause large amounts of seedling mortality and reduce brome abundance in the following growing season (Mack and Pyke ). Warmer winter temperatures in the Black Hills parks also corresponded with increased cover of annual bromes, which is consistent with patterns seen further west (e.g., Lovtang and Riegel , Bansal and Sheley ). Further research is needed to understand why cooler and drier spring conditions were associated with increased annual brome abundance in the Black Hills. The regression tree model also highlights the regional differences. In that analysis, the cover of annual bromes in Scotts Bluff NM and Fort Laramie NHS responded to winter precipitation, while in the other regions cooler winter temperatures, rather than precipitation, were associated with higher annual brome cover. These differences highlight the fact that annual bromes can exhibit a wide range of responses to climate and suggest that land managers should be cautious when making decisions based on climate responses seen in other areas.

Transect and hill slope grades were the facets of the physical environment with the strongest predictive power for annual brome cover in the NGP (Table , Fig. ). Across all parks in the NGP, annual bromes had a tendency to have the lowest cover in flat areas (< 2% grade). It may be that these areas are too hot and dry to maintain annual brome. Consistent with other studies, slope aspect and elevation showed weak predictive power for annual brome abundance in our data set. Studies have shown conflicting evidence of increasing cheatgrass abundance on north–south gradients (Bradley and Mustard , Banks and Baker ), and differences in abundance associated with slope aspect are most likely the result of complex interactions between solar radiation, climate, and resource availability (Hinds , Rickard and Warren ). In its introduced range, cheatgrass occurs in a broad span of elevations, but its greatest densities are typically between 600 and 1800 m (Stewart and Hull ), and it is limited by increased snow depth (Griffith and Loik ) and low temperatures at high elevations (Chambers et al. ). It is unlikely that elevation alone would limit annual bromes in the region covered by our data set (the highest plot was 1717 m above sea level), but it may interact with tree density, resource availability, and fire history.

We found in the linear mixed model analysis that across the NGP park units and within ecoregions, soil texture derived from soil maps did influence annual brome abundance. In general, the relationship between soil type and annual brome invasion is not well studied, but in some areas cheatgrass is commonly associated with deep, sandy soils and is less likely to be found in fine‐textured soils (Beatley ). However, silt and clay soils may result in increased water availability, which tends to favor annual bromes (Miller et al. , Bansal et al. ). Our data set suggests annual bromes in the NGP are less abundant in soils with high silt content, and this is particularly true in the short‐grass prairie (Table ). Studies have demonstrated that annual brome abundance is influenced by soil chemistry and is favored by high potassium and phosphorus availability (Belnap et al. ). Also, annual brome invasion can result in increased rates of nutrient cycling in soils and ultimately a loss of soil organic matter, further reinforcing their dominance in infested sites (Norton et al. ). Although we were not able to include annual N deposition in the linear mixed model because of its strong correlation with precipitation and park unit, the parks with higher rates of N deposition (Fort Laramie NHS and Scotts Bluff NM) had the highest annual brome cover. Annual bromes are N limited in some sites (e.g., Evans et al. , Concilio et al. ) and this may be the case in the NGP. In a recent N addition experiment in Wind Cave and Badlands NPs, annual brome cover responded positively to low and moderate levels of N addition (Symstad et al. ). Our data set relied on large‐scale soil maps for explanatory variables related to soil texture and chemistry and modeled estimates of N deposition because we did not have field measures. The potential importance of these factors in our models of annual brome abundance suggests that future monitoring should incorporate more direct measures of N loading, soil chemistry, and texture.

Conclusion

Annual exotic brome grasses are widespread in NGP grasslands, including the national park units in the region. The reductions in native species richness and forage associated with their invasion are a large concern not only for NPS land managers tasked with preserving and protecting native prairies, but also for private landowners depending on these grasslands for their livelihood. The increase in annual brome relative cover over the last two decades demonstrated in our data emphasizes the need for increased attention to these species in the NGP, but managers are struggling with limitations in available data and for proven methods for controlling annual bromes. The unique data set explored here has many limitations, but it provides valuable insight into the behavior of annual bromes over space and time in the NGP. We found that annual brome behavior in response to climatic variability in this region is consistent with expected behavior based on work in other regions. This climatic variability will have to be accounted for when designing and evaluating the effectiveness of landscape‐scale management approaches to control and reduce annual brome prevalence. Our results suggest that promising avenues to pursue are approaches that prevent excessive litter buildup and promote native plant species. Fire and grazing are two tools land managers can use to do this, but determining appropriate regimes that accomplish annual brome control at a landscape scale will require a significant, focused effort that combines replicated experiments, monitoring, and hypothesis‐driven adaptive management with creative management techniques. This is especially true in NPS units where management is traditionally light when compared with management on other public and private lands.

Acknowledgments

We thank the following NPS employees for helping us extract climate, soils, wildlife, and fire history data: J. Gross, W. Monahan, C. Boever, D. Licht, K. Hansen, S. Rockwood, L. LaFleur, and A. Jarding. We thank K. Paintner‐Green for providing valuable discussions. Comments from S. Bansal, T. Rodhouse, C. McIntyre, J. Jonas, and two anonymous reviewers improved this manuscript. This work was funded by the National Park Service, with support for AJS from the U.S. Geological Survey. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.