Spatiotemporal patterns of rising annual plant abundance in grasslands of the Willamette Valley, Oregon (USA)

Plant communities are undergoing compositional changes that affect ecosystem function. These changes are not always uniform across the landscape due to heterogenous topographic and edaphic conditions. To predict areas most at risk of change, it is necessary to identify the landscape drivers affecting plant abundance. Annual plants are increasing across the western USA, largely driven by non-native annual invasions. Here, we quantified change in annual plant abundance and identified landscape factors contributing to that change over the past 35 years. We focused on Willamette Valley (Oregon) grasslands because they represent a new example in this phenomenon. To understand the spatiotemporal patterns of annual plant abundances between 1986 and 2020, we combined a remote-sensing vegetation cover dataset from the rangeland analysis platform with gridded soils data and topographic variables. We determined the rate of change in percent cover for each 30 × 30 m pixel and regressed cover against heat load, soil depth and sand content for > 5975 hectares to determine areas most sensitive to rising annual cover. We found a tendency toward increasing annual cover, with a median gain of + 15% cover among significantly increasing pixels. However, change was uneven across the landscape, with annual cover increasing markedly in areas with high heat load and shallower soils. We identified steep, south-facing slopes as being particularly sensitive to rising annual cover. Annual plant invasions may be lagging in this region compared to elsewhere in the western USA, but trends here suggest it may just be a matter of time.


Introduction
Global environmental change is contributing to plant invasions around the world (Dukes and Mooney 1999), with novel climate and resource conditions presenting opportunities for species to colonize new habitat and increase in abundance (Burke and Grime 1996;Thuiller et al. 2008;Jauni et al. 2015).However, the effects of global change for plant invasions often play out differently depending on geographic Vol:.( 1234567890) factors on the landscape.At a landscape scale, heterogeneity in topographic and edaphic conditions that affect hydrology, nutrient quality, or composition of the existing biotic community can dictate where and how rapidly introduced species take advantage of opportunities (Theoharides and Dukes 2007).To predict possible new invasion fronts, it is necessary to identify the landscape drivers leading to changes in plant abundances.
The historic and ongoing expansion of annual species across the western United States is one of the largest plant invasions in the world (Young and Evans 1970;Jackson 1985;Seabloom et al. 2003;Norton et al. 2007;Germino et al. 2016;Bishop et al. 2019;Smith et al. 2022).In this region, nonnative annuals (predominantly grasses such as Bromus spp., Taeniatherum caput-medusae, and Vulpia spp.) were first introduced through livestock grazing during Euro-American settlement, dating back to the mid-19th century (Jackson 1985;Knapp 1996).Since then, rising annual plant dominance has displaced much of the native perennial vegetation and resulted in ecological consequences including altered soil nutrient and moisture regimes (Melgoza et al. 1990;Dyer and Rice 1999;Prevéy and Seastedt 2014), increasing fire frequency (Fusco et al. 2019;Davies et al. 2021), and declines in soil carbon (Koteen et al. 2011) and native plant diversity (Davies 2011;Brooks et al. 2016), among others.Understanding how the landscape shapes annual plant dominance is thus important for predicting areas vulnerable to future invasion and association losses in ecosystem function and biodiversity.
To date, much of the focus on annual plant invasions in the western US has been concentrated in California grasslands, the Great Basin, and other arid and semi-arid ecoregions, with less attention paid to Mediterranean grasslands of the western Pacific Northwest (spanning from southwestern Oregon to southwestern British Columbia, Canada).Previously dominated by native perennial bunchgrasses (e.g., Festuca roemeri and Deschampsia cespitosa) and a diversity of forbs, this region now tends to be dominated by non-native perennial grass species (e.g., Schedonorus arundinaceus and Arrhenatherum elatius) (Sinclair et al. 2006), where the cooler, moister conditions likely allowed these grasslands to largely resist a shift toward annual dominance during Euro-American settlement (Everard et al. 2010).
However, rising temperatures and increasing precipitation variability are amplifying the seasonality of the Mediterranean climate system in the western Pacific Northwest, causing winters to become warmer and summers to become longer, hotter and drier (Mote and Salathé 2010;Dalton and Fleishman 2021), conferring an advantage to annuals at the expense of perennial species.Indeed, two recent climate manipulation experiments in the western Pacific Northwest found that warming plots by + 2.5 °C increased the relative abundance of non-native annuals, which benefitted from reduced temperature limitation earlier in the season (Pfeifer-Meister et al. 2016;Reed et al. 2021b).As a result, western Pacific Northwest grasslands may experience a new invasion wave as climatic conditions become increasingly favorable for annual plant dominance, similar to what has already occurred in California grasslands.
The interplay of climate change with landscapelevel determinants of heat and water stress, such as slope, aspect, soil depth and texture, may create "hot spots" of initial annual plant invasions in the western Pacific Northwest.Indeed, despite large regional differences in climate and land-use histories, these spatial factors have been important predictors of annual plant abundance in both California and the Great Basin.In California, for example, non-native annual grasses have largely displaced native perennial vegetation from the hot and dry interior Central Valley and surrounding foothills (Jackson 1985;Clary 2012).Perennial grasses (both native and non-native) tend to only persist toward the cooler, moister coast and their abundance is negatively correlated with slope (Clary 2012), likely due to reduced water availability associated with greater runoff and shallower soils.Annual plants also tend to have a competitive advantage in shallower soils due to their concentration of root growth and soil-water utilization in the upper soil profile (Holmes and Rice 1996).Similarly, in the Great Basin, areas with higher heat load (steep, south-facing slopes) and coarser soil texture (i.e., areas with greater water stress) have been associated with higher non-native annual grass abundance (Chambers et al. 2007;Condon et al. 2011;Reisner et al. 2013).
A lesson from California and the Great Basin, however, is that once some annual species have sufficient seed accrual, their populations may expand beyond these initial, environmentally determined Vol.: (0123456789) hot spots.For example, patches of non-native annual grasses (Bromus diandrus and Avena fatua) in California produce such high propagule pressure that they can rapidly take over neighboring perennial-dominated patches through seed rain alone following disturbance (Larios et al. 2013).Likewise, the invasive annual forb yellow starthistle (Centaurea solstitialis) is able to overcome factors that initially limit its success in California serpentine communities through long-term buildup of propagules (Carr et al. 2019).These mass effect phenomena allow areas initially resistant to invasion to be overwhelmed by invader propagule pressure, leading to invader dominance despite environmental conditions previously acting as invasion filters (D' Antonio et al. 2001;Von Holle and Simberloff 2005;Simberloff 2009).Thus, as invading annuals become increasingly abundant and build up propagule loads, the strength of landscape drivers that initially led to their heterogenous abundance may appear to reduce over time as they expand and take over more of the landscape.
In anticipation of ongoing climate change in the Pacific Northwest, it is critical to understand how annual plants have changed in abundance to date and what areas of the landscape are most sensitive to future change.Understanding the topographic and edaphic factors that may be driving this phenomenon will be critical to facilitating the conservation and restoration of these ecosystems in the future.Moreover, documenting a shift in annual plant abundance at a landscape scale provides an opportunity to test shifts between landscape-level determinants of initial annual establishment and subsequent mass-effect patterns of spread.Here, we combine a large remote sensing vegetation cover dataset from the Rangeland Analysis Platform (https:// range lands.app/) with soils data from the USDA's Gridded Soil Survey Geographic Database (gSSURGO) and topographic variables from digital elevation models to understand the spatiotemporal patterns of annual plant abundance in western Pacific Northwest grasslands, concentrating on the Willamette Valley in Oregon.We address the following questions: 1. Have grasslands in the western Pacific Northwest exhibited an increase in percent cover of annual species over the past 35 years? 2. Is there evidence of a spatiotemporal shift in annual cover than can be related to certain land-scape factors?What topographic and edaphic conditions are most sensitive to rising annual cover?
We hypothesize that western Pacific Northwest grasslands have shifted toward greater annual percent cover since 1986, concurrent with ongoing climate change.Furthermore, we expect that annual species cover would have been historically associated with steep, south-facing slopes and/or shallower, sandier soils where heat load is highest and available water storage is low.We hypothesize that annual cover fidelity to these conditions has decreased through time, as climate change and a buildup of propagules has allowed for increasing colonization of alternative microsites.

Materials and methods
Our study takes advantage of recently developed raster datasets from the Rangeland Analysis Platform (Allred et al. 2021;Jones et al. 2018; described in greater detail below) to understand the spatiotemporal change in annual herbaceous cover in this region.We investigated the importance of different landscape drivers leading to that change using a set of variables with known relevance to annual plant invasions in the western US, namely heat load (a composite of slope and aspect), soil depth and soil texture (sand content).We focused our analyses on a set of conserved sites in the Willamette Valley, Oregon with known natural grassland habitat.Whereas many grasslands in California and the Great Basin have already shifted toward annual dominance (> 50% relative cover), which has been reported on extensively, grasslands in the western Pacific Northwest may represent a new example in this phenomenon as climate change makes conditions increasingly favorable for annual plant species.We focused on the Willamette Valley in particular (as opposed to western Washington and/or British Columbia grasslands) for two primary reasons: because, in general, grasslands of the Willamette Valley ecoregion are taxonomically and functionally representative of others across the Puget Lowland (Washington) and Georgia Basin (British Columbia) ecoregions (Dunwiddie and Bakker 2011), and because its geographic location toward the southern end of the region (where it is relatively warmer Vol:.( 1234567890) and drier) make it particularly susceptible to rising annual dominance (Pfeifer-Meister et al. 2016;Reed et al. 2021b).The Willamette Valley also features the most extensive grassland habitat across the Willamette Valley-Puget Lowland-Georgia Basin ecoregions (~ 800,000 hectares in the Willamette Valley (Christy and Alverson 2011).
We focused on conserved sites with natural grassland habitat since much of the land classified as "grasslands" in the Willamette Valley are managed agricultural fields of grass seed production, which we sought to avoid.We selected 16 sites (Fig. 1a) identified by the Willamette Valley Oak and Prairie Cooperative as having permanent conservation status, containing a considerable component of prairie and savanna habitat, and being at least 40 hectares in size (WVOPC, 2020).Additionally, we selected these sites out of a set of candidates based on the criteria of having relatively wide distributions of soil depth, texture, slope, and aspect, as well as inspection of satellite imagery to verify the presence of grassland habitat.All sites exist within a relatively narrow elevational range between 45-644 m a.s.l., typical of land across the Willamette Valley.Representative photos for 11 of the 16 sites can be found in the Supporting Information  to provide examples of the dominant plant community compositions within these grasslands.

Vegetation cover data
The Rangeland Analysis Platform (RAP) is a thoroughly vetted online resource of fractional vegetation cover estimates for the United States, produced by learned multitask models correlating on-the-ground data to satellite-derived imagery through regression tree approaches (Allred et al. 2021).It provides yearly vegetation percent cover estimates from 1986 to present for the contiguous United States at 30 m pixel resolution (Jones et al. 2018;Allred et al. 2021).In the most recent release (Version 3.0), these estimates are provided for six different coverage classes: annual forbs and grasses (hereafter "annuals"), perennial forbs and grasses (hereafter "perennials"), shrubs, trees, litter and bare ground.
These data were produced through an observationbased machine learning model that related nearly 75,000 line-point intercept plots in the field to the historical Landsat Collection 2 satellite imagery record.Field vegetation data comes from the Natural Resource Conservation Service's Natural Resources Inventory program, the Bureau of Land Management's Assessment, Inventory and Monitoring program, and the National Park Service's Northern Colorado Plateau Network program and was collected between 2004 and 2019 in uncultivated, undeveloped rangelands across 17 western US states.The algorithm is able to differentiate between annuals and perennials by using Landsat 5 TM, 7 ETM + and 8 OLI surface reflectance measurements alongside Landsat vegetation indices calculated by normalized difference vegetation index and normalized burn ratio two as predictors of fractional cover (Allred et al. 2021).Generally, perennial and annual herbaceous vegetation conveniently separate out by these measurements, giving distinct signatures that differentiate them.The RAP algorithm predicts functional group percent cover estimates with high accuracy, with annual forbs and grasses having a mean absolute error of 7.0% (Allred et al. 2021).This absolute error refers to the difference in magnitude between the prediction of the RAP model with the true value of that observation (i.e., the on-the-ground cover measurement for any given pixel).
It is important to note that the RAP dataset is unable to distinguish between native and non-native vegetation or between grasses and forbs, within the annuals functional group.This means that all of the data presented here is only descriptive of all annual plants as a whole and is not empirical evidence for a rise in any particular subset of annual plants, such as non-native annuals or annual grasses.However, when we interpret our results throughout the discussion, we point to an abundance of evidence from other published studies to suggest that the results we observe are likely driven by non-native annuals in particular.To provide additional supporting evidence for this suggestion, we include ranked abundance data of herbaceous species observed within twenty 100 m 2 relevé plots collected from two study sites in 2021 (Howard Buford Recreation Area and South Eugene Meadows) as a way to ground-truth the dominant annual and perennial species within these grasslands.These data, along with a brief description of their collection methods, can be found in the Supporting Information (Supplemental Methods, Figure S13, Tables S1 and S2).
Using Google Earth Engine, we downloaded RAP data files (coordinate system WGS 1984) of the Willamette Valley for annuals and perennials for each year between 1986 and 2020.To eliminate forested and other non-grassland areas, we masked RAP data to grassland/herbaceous, pasture/hay, and shrub/scrub classifications in the National Land Cover Database 2019 Land Cover dataset (Dewitz, 2021) using the raster package (Hijmans, 2020) in R version 4.0.2(R Core Team, 2020), and clipped data to polygons of our 16 study sites.Prior to this, we had previously resampled the land cover dataset from 10 to 30 m resolution using nearest neighbor assignment and reprojected to WGS 1984 using ArcMap 10.5.This resampling technique uses the cell center from the input raster (i.e., the 10 m scale) that is closest to the processing cell center (i.e., the 30 m scale) and assigns it as the output value for the processing cell.It is important to note that resampling can add uncertainty to an analysis because one is interpolating data and potentially losing information by going from a finer to a coarser scale.However, this step was necessary because all raster datasets in the analysis needed to be at the same spatial resolution.
For each year in the 1986-2020 RAP dataset, we calculated new rasters for total herbaceous cover by summing the annual and perennial layers.We then filtered to cells with ≥ 20% average total herbaceous cover across the 35 years to avoid areas which may have skewed estimates due to low herbaceous cover in general (e.g., grassland borders near forests).Upon masking and filtering, we had 66,391 pixels of data, accounting for > 5975 hectares.

Topographic and edaphic data
For soil data (depth and texture), we obtained 2019 gSSURGO data (10 m resolution) for the state of Oregon through the USDA Natural Resources Conservation Service web portal (Soil Survey Staff).Within this geodatabase, each soil map unit corresponds to one or more soil polygons represented in raster form.Each map unit contains several soil component records ('chorizon' table records), each representing a soil horizon and each having a percent attribute indicating their relative proportional composition within the map unit.For soil depth, we used the 'brockdepmin' variable in the 'muaggatt' (map unit-aggregated) table when available.For map unit records with no data for 'brockdepmin', we calculated a minimum depth by choosing the deepest non-rock horizon of each component (i.e., the max 'hzdept_r' from the 'chorizon' table, indicating the distance from the top of the soil to the upper boundary of bedrock), and then averaging these numbers weighted by the component percent.For percent sand, we used the attribute for the topmost horizon in the dominant components of the mapunit ('sandtotal_r' from the 'chorizon' table where 'hzdept_r' = 0).We then converted these layers to WGS 1984 and resampled from 10 to 30 m resolution using ArcMap 10.5, again using the nearest neighbor assignment.
Aside from soil depth and sand content, we also considered but excluded other relevant soil predictors including silt and clay content, organic matter, pH and available water storage to 100 cm depth.For texture, we decided to only use sand content since silt content explained little variation and sand and clay content were largely interchangeable as they were inversely related.The pH and organic matter data from the gSSURGO database were at such a coarse scale that they both had very little variation across our study sites, making inference difficult.We decided to exclude available water storage because this was generally a byproduct of heat load and soil depth (i.e., areas with high heat load and/or shallow soils tend to have have low water storage while areas with low heat load and/or deeper soils have higher water storage).
For topographic variables (slope, aspect and heat load), we downloaded a 30 m resolution digital elevation model file from OpenTopography (NASA, 2013).To calculate slope and aspect, we used the slope and aspect tools from the spatial analyst toolkit in Arc-Map 10.5.Heat load is a unitless index of the amount of heat a ground surface receives through solar radiation based on its slope, aspect and latitude.We calculated heat load in R using the McCune and Keon (2002) method with south-facing slopes (180°) at a maximum.Heat load is a useful composite variable to represent the inherent interaction between slope and aspect on solar radiation, given that steep, south-facing slopes receive the greatest radiation while steep, north-facing slopes receive the least radiation (in the northern hemisphere).
In our analyses, we considered using elevation as an additional topographic variable and slope and aspect explicitly instead of heat load (using northness as a linear index of aspect calculated using the cosine of aspect (in radians, θ), where − 1 indicates due south and + 1 indicates due north).However, we excluded these variables because elevation varied minimally across our study sites and heat load combines the interaction of slope and northness in a biologically meaningful way with just a single composite variable.

Analyses
To address our first question of whether western Pacific Northwest grasslands have exhibited an increase in annual cover over the past 35 years, we used our processed RAP data in R to conduct pixel by pixel linear regressions of annual percent cover against year.We then generated a new raster layer for the rate of change using the coefficients from these models.For our second question regarding the topographic and edaphic conditions most sensitive to rising annual cover, we first averaged yearly annual cover estimates to half-decade intervals to reduce some noise generated by strong year to year fluctuations in vegetation cover.We then modeled this average annual cover variable (using a logit-transformation to bound estimates by 0-100%) against heat load, soil depth, percent sand and two-way interactions of each of these variables with time (a continuous representation of half-decades, where 1 = 1986-1990, 2 = 1991-1995 and so on), using centered and scaled values for all predictors.Collinearity among heat load, soil depth, and percent sand was small (r ≤ 0.41; Figure S1).Using the 'quantreg' package in R (Koenker, 2020), we fit quantile regressions to characterize spatiotemporal trends in the 10th, 50th, and 90th percentiles of annual cover.We used quantile regressions since we did not expect that all pixels would be changing to the same degree and we wanted to understand the drivers behind those differences.In addition to estimating the central tendency, quantile regression allowed us to identify and characterize trends for the most rapidly annualizing pixels (90th percentile) as well those lagging at the slowest rate or even losing annual cover (10th percentile).

Results
We found that grasslands across the western Pacific Northwest have generally increased in annual percent cover since 1986 (Fig. 1).Across the 66,391 pixels in our analysis, 16,903 (25.5%) exhibited a significant increase in annual cover (p < 0.05), 8,088 (12.2%) exhibited a significant decline and 41,400 (62.4%) exhibited no significant change (Fig. 2).Among cells that significantly annualized (i.e., gained annual cover), the median rate of change was + 0.429% cover year −1 , resulting in an increase of 15% annual cover across the 35 years.Among cells that exhibited a significant decline in annual cover, the median rate of change was − 0.231% cover year −1 , resulting in an 8% annual cover decline across the 35 years.Across all pixels, the total median rate of change was + 0.052% cover year −1 (for an overall 1.8% annual cover gain across the 35 years; Fig. 2).Individual sites differed idiosyncratically, with South Eugene Meadows exhibiting the greatest median rate of change across all pixels (+ 0.209% cover year −1 ; Figure S2K) and Rattlesnake Butte exhibiting the lowest median rate of change across all pixels (− 0.04% cover year −1 ; Figure S2G).
As expected, high annual cover was more prevalent toward steep, southern aspects (positive heat load effect) and shallower soils (negative depth effect) for all three quantiles (Fig. 3; Table 1).With percent sand, the 90th quantile of annual cover had a positive relationship whereas the 50th and 10th quantiles had negative relationships (Fig. 3).
Fig. 2 Distribution of coefficient values in the pixel by pixel regressions of annual forb and grass (AFG) percent cover against year.Grey bins denote pixels that did not exhibit a significant change in annual cover (p > 0.05), purple bins exhibited a significant decline while orange bins exhibited a significant increase.The solid black vertical line displays the population median value (0.052% cover year −1 ), while the purple and orange vertical lines display the median value among pixels in which annuals significantly declined and increased, respectively Fig. 3 Standardized coefficient effect sizes from quantile regressions of average annual forb and grass percent cover (logit-transformed) against landscape predictor variables and their interactions with time.All predictor variables were centered and scaled before modeling.Solid symbols denote a significant effect while the hollow symbol indicates a nonsignificant effect (Table 1) Vol:.( 1234567890) Counter to our hypothesis, we found that annual fidelity toward areas with high heat load was amplified through time (positive time × heat load interaction) for all three quantiles (Fig. 4a).Based on our model results, the 10th, 50th, and 90th quantiles of annual percent cover at highest heat load rose, respectively, from ~ 4%, 9%, and 20% in 1986-1990 to ~ 6%, 14%, and 40% in 2016-2020.Contrarily, at lowest heat load, the 10th, 50th, and 90th quantiles declined, respectively, from ~ 6%, 9%, and 9% in 1986-1990 to ~ 1%, 4%, and 5% in 2016-2020.Additionally, fidelity to shallower soils was amplified through time for the 90th and 50th quantiles (albeit minimally for the latter) but decreased through time for the 10th quantile, again with a relatively minimal effect size (Fig. 4b).At the 90th quantile, annual percent cover rose from ~ 20-34% in the shallowest soils and from ~ 13-18% in the deepest soils across the 35 years.The relationships between annual cover and sand changed minimally through time (Fig. 4c).
Patterns of annual cover against landscape factors (including slope and aspect) were largely consistent across sites (Figures S3-S7), although some factors had too little variance within sites to generate much inference (e.g., soil depth at Chip Ross Park (CRP); Figure S6A).Maps of each landscape factor can be found in the supporting information, Figures S8-S12.

Discussion
Using a recently developed remote-sensing data product (Jones et al. 2018;Allred et al. 2021), we provide a spatially explicit, fine resolution analysis to show the change in annual percent cover in grasslands of the western Pacific Northwest over the past 35 years.Overall, there has been a slight trend toward increasing annual cover, but this has been highly variable across the landscape.In particular, rising annual cover has occurred to the strongest degree in areas with high heat load as well as (but to a lesser extent) in areas with shallower soils.Conversely, areas with lowest heat load tended to be those that exhibited a decline in annual percent cover through time.Despite dramatic year to year fluctuations in annual herbaceous cover, probably linked to interannual precipitation variability (Bradley and Mustard 2005), we detected a significant positive trend for approximately one quarter of all pixels.At the moment, this remains only a subset of the landscape, suggesting that environmental filtering, rather than mass effects (when propagule buildup allows for colonization despite environmental filtering), is still largely shaping where annuals are found in abundance.
Our analysis reveals that steep, south-facing slopes are likely the most sensitive areas on the landscape to rising annual cover, and that this sensitivity has been increasing through time since the late 1980s.The fact that annual abundance is greatest in areas with high heat load is expected and consistent with studies from the Great Basin on cheatgrass (Bromus tectorum) and other non-native annual grasses (Chambers et al. 2007;Condon et al. 2011;Reisner et al. 2013).However, this result is counter to a remote sensing study from the Colorado Plateau which found cheatgrass hotspots associated with flatter slopes (Bishop et al. 2019).Nonetheless, the pattern we observed may be revealing for consequences on native perennial vegetation, as one study from California found that native grasses (five of six being perennial) were least abundant while non-natives (13 of 15 being annual) were most abundant on warm, south-facing slopes of (Gelbard and Harrison 2003).As annuals continue to rise within high heat load areas, they may be displacing perennial vegetation.We found that the strength of the relationship between annual cover and heat load continued to increase over time.This suggests that this region is in the relatively early stages of invasion, such that annual species have been overcoming dispersal limitation over time or that environmental change is increasingly selecting for annual species or a combination of the two.This is counter to our hypothesis, which predicted a decrease in the strength of environmental filtering due to annual seed accrual and mass effects.It seems that this expectation may have been ahead of schedule: it is possible that climatic conditions in low heat load regions are still not adequate for this expansion and/or that the more favorable high heat load regions are simply not yet at capacity of annual abundance.Most pixels still have relatively low annual percent cover even in highest heat load areas (median of only 14% cover in highest heat load areas by 2016-2020), so there may still be time before propagule buildup, spillover and expansion takes place.
To our knowledge, no studies in other regions of the western US suggest a threshold value for annual cover that represents a turning point in the mass effects/seed accrual response.Any threshold would likely vary dramatically based on dominant species and site conditions.With that said, data in Larios et al. (2013) from California suggest that when the introduced annual grasses Bromus diandrus or Avena fatua reach a relative abundance of approximately 25% in a native-dominated patch, this is enough  val ('86-90 = 1986-1990; '91-95 = 1991-1995; and so on).Black lines show model predictions back-transformed from the logit scale.Top line is the 90th quantile, middle line is 50th quantile, and bottom line is 10th quantile.Background datapoints are colored by soil depth in a and heat load in b and c to trigger a dominance shift to introduced annual dominance (> 50% relative abundance) the following year due to heavy seed rain.Both of these nonnative annual grasses are common representatives of this functional group within our study region, so it is possible that as little as 25% annual cover could trigger a mass effect response in some cases.Additionally, across the Great Basin, areas identified using RAP data as having "annual grass dominance" have a median annual cover of just 29% (Smith et al. 2022), suggesting that it does not require extremely high cover for an area to be considered annual dominant.Thus, when considered at this scale, our study sites may be behind in the invasion process, but perhaps not immensely so.
Perhaps more revealing is the observation that the lowest heat load areas (i.e., steep, north-facing slopes) are generally losing annual percent cover through time (Fig. 4a).This suggests that these areas of the landscape have been highly resistant to annual invasion thus far or that enough propagules have simply not been introduced yet.Alternatively, this also could be a sign that native annual species are being further squeezed out from these parts of the landscape by perennial vegetation.Like in California grasslands (Gelbard and Harrison 2003), these cooler, moister areas may be serving as refugia for and fortification of perennial vegetation.Indeed, lower heat load has been suggested to increase perennial grass cover in Great Basin rangelands (Davies et al. 2007), and greater perennial cover increases resistance to annuals (Chambers et al. 2007).Eventually, propagule pressure from neighboring high heat load areas may be enough to overcome the strong annual plant filters in low heat load areas (D 'Antonio et al. 2001;Von Holle and Simberloff 2005;Simberloff 2009), but it appears that this time has not arrived yet.
While heat load was the primary environmental driver, soil depth was also a substantial control on patterns of increasing annual cover.Areas with deeper soils are not as sensitive to rising annual percent cover as are areas with shallower soils, yet even the deepest soils are still exhibiting a small rise in annual abundance, particularly at the 90th quantile (Fig. 4b).This finding of greater sensitivity in shallower soils agrees with a recent study from eastern Oregon on the invasive annual grass Ventenata dubia, which attributed invasion susceptibility in shallow-soil sites to lower resident cover and biomass (Tortorelli et al. 2022).This is also consistent with another study from eastern Oregon in 2009 which found higher non-native annual grass abundance in xeric (shallow soil) relative to mesic (deep soil) sites (Endress et al. 2020).However, those particular sites in the Endress et al. (2020) study were resampled in 2016 and were found to have had a significant increase in non-native annuals (particularly V. dubia) in the mesic but not xeric sites over that seven year time span (Averett et al. 2020).Likewise, a study in California grasslands suggested that deeper soils are more invasible by the non-native annual grass Lolium multiflorum (Kolb et al. 2002), which they attributed to higher nitrogen or water availability.In these latter two examples, it may be that results were a product of a later stage of the invasion or due to the specific environmental tolerances of these two particular species.If shallower soils were invaded first, then studies that occur after the primary invasion may find expansion toward deeper soils and draw misleading conclusions.Alternatively, these results may be contextualized by regional differences or differences in species environmental tolerances.Given that California is hotter and drier than the western Pacific Northwest, soil depth may be less of a mediator in those circumstances.
As for soil texture, whether an area had high or low sand content did not seem to matter in terms of its susceptibility to rising annual cover, especially for the central tendency (50th quantile; Fig. 4c).In the Great Basin, some invading annuals like cheatgrass (B.tectorum) are associated with sandier soils (Reisner et al. 2013) while others like medusahead (T.caput-medusae) seem to prefer clayey soils (Young and Evans 1970;Applestein and Germino 2022).In the semi-arid Colorado Plateau, cheatgrass hotspots have been found to be associated with higher clay contents (Bishop et al. 2019).Thus, differing soil texture niches among species within the annual functional group may be negating a general trend in response to sand content.Alternatively, a lack of a strong trend could also be explained by the relatively low variability in sand content across pixels in the gSSURGO database, as evident by the sparsity of datapoints with > 40% sand (Fig. 4c).This low variability was also true for clay and silt content, which likewise did not exhibit trends with annual percent cover (data not shown).
Despite clear patterns in annual percent cover through time and in relation to some landscape Vol.: (0123456789) factors, the quantile regressions explained relatively little variation overall (pseudo-R 2 values of 0.14, 0.02 and 0.02 for the 90th, 50th and 10th quantiles, respectively).This indicates that there are still many other unexplained factors at play which affect where and how strongly annuals occur.Most notably, our analysis does not consider the disturbance or management history of the landscape, despite these being critical drivers of annual invasions in the western US (Bradford and Lauenroth 2006;Brooks et al. 2016;Williamson et al. 2020).Ideally, variables such as time since conversion of land from native prairie habitat to agricultural usage and vice/versa or the frequency and type of disturbances (e.g., prescribed fire, mowing, herbicide use, grazing, etc.) might be informative and relate to annual cover through time.For example, areas with relatively low heat load and deep soils (e.g., parts of Finley National Wildlife Refuge; Figures S8, S11) may still be highly invaded by annuals if they had been heavily disturbed and had an abundance of propagules introduced.Conversely, areas with relatively high heat load and shallow soils (e.g., Native Oaks Ridge; Figures S8, S11) may have low annual cover if they have not been disturbed and are still dominated by relatively drought-tolerant native perennial bunchgrasses such as Danthonia californica, or that conditions are so harsh that total herbaceous cover is low in general.Indeed, although we found general trends across sites, the large site by site idiosyncrasies in annual cover may reflect unique land-use histories and seed availability.The primary reason we left management and land use history out of model development was because these data do not exist in raster form at this resolution.These other factors are important in the context of our results, and local knowledge of site conditions should always be considered before interpreting remotely sensed data products such as this for on the ground management (Allred et al. 2022).
While remote-sensing data provides a powerful way to understand changes over large spatial gradients and longer timespans, there are a few limitations associated with this type of analysis.In particular, we are unable to differentiate between native and non-native annuals or between grasses and forbs.We make the assumption that the rise in annual cover observed here is most likely amongst non-native annuals (and grasses in particular) as this is the overwhelmingly strong pattern observed in other western US regions (Norton et al. 2007;Germino et al. 2016), and because other evidence from western Pacific Northwest prairies suggests native annuals have declined precipitously since the 19th century (Dunwiddie et al. 2014) and are poised to continue declining with climate change (Pfeifer-Meister et al. 2013;Reed et al. 2021a, b), while non-native annual grasses are poised to increase (Pfeifer-Meister et al. 2016;Reed et al. 2021b).We recognize the possibility that the rise in annual cover observed here could indicate a rebound in native annual species, and that our data as presented give no indication to the validity of our assumption.However, no evidence exists to our knowledge to suggest that this is a possibility whereas all the available evidence points to the contrary (that non-native annuals are the ones on the rise).Indeed, in addition to published studies, local land managers have also been observing an alarming increase in non-native annual grasses such as Vulpia myuros, Ventenata dubia, and Taeniatherum caput-medusae across the Willamette Valley in recent years (P.Reed, personal communications).
Additionally, there is some inherent error within the satellite-derived RAP dataset (Allred et al. 2021), and the edaphic factors from gSSURGO seem to be relatively coarse estimates across the landscape.Finer scale estimates of soil conditions and inclusion of other variables such as nutrient availability (which is not currently available in the gSSURGO database) could improve predictive performance.
Even with these caveats, our results provide important implications for local grassland management and restoration.In particular, Willamette Valley sites with highly sensitive landscape factors (e.g., steep, south-facing slopes and shallow soils) but intact perennial communities may be poor areas to target for restoration, as there could be high risk in converting these areas to annual dominance following disturbance (Brambila et al. 2022).Indeed, postfire recovery in a nearby forested ecosystem saw non-native annual grasses taking over in warm, dry areas (Dodson and Root 2015), and others in the Great Basin have suggested avoiding disturbance in xeric sites to prevent the spread of cheatgrass dominance (Condon et al. 2011).Conversely, there may be lower risk (for now) associated with disturbing perennial sites that have less sensitivity (e.g., northfacing slopes), and thus these areas should be higher priorities for restoration.It is important for land managers in this region to be aware of the increasing risk of annual dominance to prevent unintended consequences to ecosystem functions.
Overall, this study is a novel example of using a remote-sensing data product (the Rangeland Analysis Platform) in conjunction with other large geospatial datasets (the USDA's gSSURGO database) to understand the landscape drivers contributing to a change in plant functional group abundance through time.For grasslands in the western Pacific Northwest, the change in annual percent cover between 1986 and 2020 has not been uniform across the landscape.Increases in abundance have outweighed decreases and have occurred more often in areas with high heat load and shallow soils, while decreases have generally occurred in areas with low heat load.Western Pacific Northwest grasslands may be lagging behind the explosive annual invasions that have already resulted in their widespread dominance in California or ongoing rapid expansion in the Great Basin, but trends here suggest it may just be a matter of time.

Fig. 1 a
Fig. 1 a Map of the Willamette Valley in Oregon, depicting the rate of change in annual forb and grass (AFG) percent cover between 1986-2020 across the 16 study sites.Orange indicates an increase in annual cover; purple indicates a decline.Insets B-Q provide zoomed-in views of each site.b Baskett Slough National Wildlife Refuge; c Ankeny National Wildlife Refuge; d Kingston Prairie Preserve; e Chip Ross

Fig. 4
Fig. 4 Results from quantile regressions showing the effects of a time x heat load (with soil depth and sand held at their means), b time x soil depth (with heat load and sand held at their means), and c time x percent sand (with soil depth and heat load held at their means) on annual forb and grass (AFG) percent cover.Each individual plot shows a half-decade inter-

Table 1
Standardized coefficient estimates of quantile regression modelsPredictor variables were mean-centered and scaled by their standard deviation Vol.: (0123456789)