Is the grass always greener? Land surface phenology reveals differences in peak and season‐long vegetation productivity responses to climate and management

Abstract Vegetation phenology—the seasonal timing and duration of vegetative phases—is controlled by spatiotemporally variable contributions of climatic and environmental factors plus additional potential influence from human management. We used land surface phenology derived from the Advanced Very High Resolution Radiometer and climate data to examine variability in vegetation productivity and phenological dates from 1989 to 2014 in the U.S. Northwestern Plains, a region with notable spatial heterogeneity in climate, vegetation, and land use. We first analyzed interannual trends in six phenological measures as a baseline. We then demonstrated how including annual‐resolution predictors can provide more nuanced insights into measures of phenology between plant communities and across the ecoregion. Across the study area, higher annual precipitation increased both peak and season‐long productivity. In contrast, higher mean annual temperatures tended to increase peak productivity but for the majority of the study area decreased season‐long productivity. Annual precipitation and temperature had strong explanatory power for productivity‐related phenology measures but predicted date‐based measures poorly. We found that relationships between climate and phenology varied across the region and among plant communities and that factors such as recovery from disturbance and anthropogenic management also contributed in certain regions. In sum, phenological measures did not respond ubiquitously nor covary in their responses. Nonclimatic dynamics can decouple phenology from climate; therefore, analyses including only interannual trends should not assume climate alone drives patterns. For example, models of areas exhibiting greening or browning should account for climate, anthropogenic influence, and natural disturbances. Investigating multiple aspects of phenology to describe growing‐season dynamics provides a richer understanding of spatiotemporal patterns that can be used for predicting ecosystem responses to future climates and land‐use change. Such understanding allows for clearer interpretation of results for conservation, wildlife, and land management.


| INTRODUC TI ON
Ecological information across space and time can enhance decisionmaking by allowing researchers and practitioners to address problems within appropriate geographies, derive landscape-scale indicators of change, and prioritize in situ monitoring efforts (Dale et al., 2000;Lausch et al., 2018;Marvin et al., 2016). For example, adequate understanding of the timing, magnitude, and duration of seasonal events (phenology) of an ecological community can inform management and restoration efforts and help managers estimate how sensitive a community may be to ongoing changes in climate.
Many aspects of phenology can be measured by satellite-based remote sensing, termed "land surface phenology," a combined measurement of variation across multiple plant species and individuals as well as background surfaces such as soils (Hanes et al., 2014).
Studies of land surface phenology have identified numerous changes across the globe in the seasonal cycle of vegetation growth due to climate, including changes to start of spring, end of season, and season-long productivity Morisette et al., 2009). Greening, measured through either increased growing-season production or higher peak productivity, is evident over considerable portions of the Earth's terrestrial surface Zhu et al., 2016). However, parts of the terrestrial surface have experienced no change or are browning-a decrease in growing-season production and/or lower peak productivity (de Jong et al., 2011;Meng et al., 2020;Zhu et al., 2016). Trends in greening, browning, and other aspects of phenological change often challenge the predictions of land surface models, suggesting that our understanding of these processes must be improved.
Phenological measurements, be they dates or magnitudes, are not always strictly increasing or decreasing; many factors play into their variability and patterns (de Jong et al., 2011;Li et al., 2019;Piao et al., 2019). Although analyses of an individual phenological measure through time (interannual trend) provide important information, the amount of the land surface that exhibits statistically significant trends in these phenological measurements reflects how constraining climate factors are changing through time and how ecosystems respond to these changes. Interannaul trends may depend heavily on the selected starting and ending years of the time series, which may mask more subtle changes in the progression of vegetation greenness in an area (e.g., Yang et al., 2021;Yuan et al., 2019). Including climatic drivers in models of phenology helps elucidate the mechanisms behind observed changes across years. In addition, there are important differences in the implications of climate change for peak productivity versus season-long productivity. Unfortunately, only peak or season-long productivity is typically investigated, which can obscure comparisons of the mechanisms underlying different aspects of phenological changes . Precipitation and temperature have been identified as important drivers of seasonal variation in productivity; however, these relationships can vary across ecosystems and vegetation types (Fu et al., 2017;Maurer et al., 2020;Reed et al., 2019;Yang et al., 2021). Furthermore, other natural and anthropogenic factors may confound this spatial dynamic, such as disturbance, land-use and land-cover change (including invasive species), CO 2 enrichment, and changes to water use and soil water-holding capacity (Nemani et al., 2003;Piao et al., 2019;. The U.S. Northwestern Plains (NWP; Figure 1) provide a meaningful geography to expand from prior studies of phenological change due to large expanses of intact native vegetation (Auch et al., 2011) interspersed by more intensively managed agricultural areas. Climatic drivers of phenology vary across gradients in the NWP, which exhibits strong west-to-east increases in precipitation and south-to-north decreases in temperature ( Figure A1; Epstein et al., 1996). Interannual variability in precipitation is high in the NWP, leading to variable productivity (Petrie et al., 2016); however, productivity of higher elevation areas, which are typically forested, has varied less over time and has responded more strongly to temperature than precipitation (Fu et al., 2017;Potter, 2020).
Collectively, these factors lead to variability in spatiotemporal patterns of phenology across the NWP. Here, interannual trends of increasing productivity outweigh browning trends; about 20% of the NWP has exhibited a significant increase in peak productivity (as measured by maximum annual NDVI) over the past few decades (Brookshire et al., 2020), and further increases are anticipated (Hufkens et al., 2016). Additionally, growing-season length has also increased, primarily due to later end-of-season dates rather than earlier start-of-season dates (Ren et al., 2020). How these different measurements of phenology combine to determine changes in vegetation seasonality and productivity across different ecosystems in the NWP remains less clear.
Consequently, we aimed to examine the spatiotemporal variability in dates and magnitudes of phenological events in the NWP and how prevailing climatic conditions may influence and/or govern this variability. Ecological models, in our case phenological ones, can be K E Y W O R D S disturbance, grasslands, productivity, remote sensing, shrub-steppe, trade-offs hierarchical, complex, and challenging to parameterize (Morisette et al., 2009;Newman et al., 2019). Therefore, to explore patterns at various scales, for example, temporal, regional, and community level, we analyzed spatiotemporal patterns in six phenological metrics across seven plant communities and an entire ecoregion.
Specifically, our objectives were to first evaluate interannual trends and then compare those findings to a model incorporating annual precipitation and temperature across space and vegetation community types of the NWP (Figure 1). We then used these findings to examine the predictive ability of annual weather variables for phenological measures and evaluate their role in patterns, such as variation of greening versus browning between different measures of productivity.

| Study area
We defined the NWP study area using a two-part process. First, we selected the Northwestern Great and Glaciated Plains ecoregions (https://www.epa.gov/eco-resea rch/ecore gions; Omernik, 1987) and then included semi-arid areas within the upper Missouri River basin, the transitional areas from the Great Basin (parts of the Middle Rocky Mountains) and Wyoming Basins (Bighorn Basin part of this ecoregion, Figure 1). We chose this boundary to include similar vegetation types of semi-arid grasslands and shrublands with intermixed agricultural and forested areas and to cover a wide range of climatic regimes. Vegetation communities range from barren and shrubland systems in the west to mixed-grass and tall-grass prairie to the east.
Forested areas are found in higher elevations, such as the western edge of the NWP, in mountain islands like the Black Hills, and adjacent to rivers. Areas converted to agriculture primarily occur in the northern and eastern portions of the study region and in irrigated areas along waterways. The NWP receives from 33 to 56 cm (13 to 22 in) of precipitation per year and mean annual temperature ranges from 3.8 to 10.5℃ (39 to 51℉), although precipitation can exceed 160 cm (63 in) per year, and annual mean temperatures average below 0℃ (32℉) at higher elevations ( Figure A1; USDA, 2006).

| Data sources and processing
We obtained phenology data derived from the Advanced Very High Resolution Radiometer (AVHRR) satellite processed by the USGS for the years 1989-2014 (USGS EROS, 2015). These data utilized composites of daily, 1 km pixel values to calculate the normalized difference vegetation index (NDVI), from which phenological measures F I G U R E 1 The Northwestern Plains study area, which comprises the Northwestern Great and Glaciated Plains Level 3 Ecoregions (Omernik, 1987), plus areas containing semi-arid vegetation in the upper Missouri River basin. These additions include portions of the Middle Rocky Mountains and Wyoming Basins Ecoregions. Land cover data from the Multi-Resolution Land Characteristics Consortium (NLCD 2016, available at www.mrlc.gov) were derived (Eidenshink, 2006;Reed et al., 1994). The general approach applied by the USGS Remote Sensing Phenology Center for these data is to create smoothed NDVI composites (Eidenshink, 2006) using a weighted least squares approach (Swets, 1999) and employ a delayed moving-average algorithm to identify dates of phenological events and their magnitudes (Reed et al., 1994). Due to degradation in the orbit at the end of life for different platforms carrying the AVHRR, solar zenith angles (SZAs) have changed through time and are associated with trends in phenological measures, degrading the temporal consistency (Ji & Brown, 2017) and precluding analyses beyond 2014. Although much of the study area had limited area with significant correlations between SZA and phenological measures, we chose to remove problematic years (1992, 1993, 1994, 1999, and 2002) to minimize issues from orbital degradation (Ji & Brown, 2017). Our comparison of results including versus excluding the five problematic years found limited large-scale effects on the results.
We used six phenological measures ( Figure 2) similar to those used in other studies of semi-arid systems (Maurer et al., 2020;Ren et al., 2020;Yang et al., 1998). These included four date-related measures (start, end, and length of season, and day of peak productivity (NDVI)) and two measures related to productivity and the shape of the phenological curve (maximum NDVI and time-integrated NDVI).
While maximum NDVI (hereafter peak productivity) is a surrogate for vegetation productivity at a single point in time, the peak of the growing season, time-integrated NDVI (hereafter season-long productivity) is a surrogate for vegetative productivity across the entire growing season.
We acquired climate data from the Parameter-elevation Relationships on Independent Slopes Model (PRISM) version AN81d at 30-arcseconds (~800 m, for our study domain) spatial resolution (Daly et al., 2008) for each calendar year of our phenology dataset . We chose annual values ( Figure A1) to provide a parsimonious overview of the complex mix of phenological measures and vegetative communities. In effect, annual measures are a relative index of the average climate conditions that each "pixel" experiences over the study period. We recognize the decreased fit and sensitivity that may occur by combining seasons, in addition to a reduction in identifying the mechanistic understanding of seasonal factors influencing phenology. By reducing the complexities, this approach allowed us to assess whether relatively simple models provide meaningful predictions that can connect with annual climatic variability.  (Wickham et al., 2014, Yang et al., 2018. There are sixteen classes, including developed, agricultural, and natural land-cover classes. We aggregated these 30-m pixels to assign the primary (modal) land-cover class within each 1-km AVHRR pixel for 2001, 2004, 2006, 2008, 2011, and 2013. We then identified AVHRR pixels where the modal land-cover class remained the same for every year and used those in analyses by land-cover type. Values for land surface phenology measures over heterogeneous areas can be driven by a subset of the area (e.g., a subdominant land-cover class) (Zhang et al., 2017). Therefore, we also used the change index NLCD layer to calculate the per AVHRR pixel (1 × 1 km resolution) percentage of Landsat pixels (30 × 30 m resolution) that changed between 2001 and 2016. This index was used to measure any subdominant land-cover changes in the study area ( Figure A2) to assess whether such changes may be influencing trends. Our analysis considered the land-cover classes of barren (i.e., vegetation generally minimal or constituting a minority of the pixel, e.g., South Dakota badlands), deciduous forest, evergreen forest, shrubland, grassland, agriculture (a combination of pasture, intensively managed under an annual or perennial cycle for feed, grazing, crops, and/or woody agriculture), and wetlands (combined woody and herbaceous wetlands) (Table 1).

| Data analysis
We conducted two sets of analyses to demonstrate how the spatial variation in phenology measures and observations, such as greening, can be enriched with the inclusion of climate variables that are likely drivers of potential trends. First, we modeled interannual trends for each of our six phenology measures ( Figure 2). For each AVHRR pixel, we calculated the Sen's slope to account for outliers, skewed data, and heteroscedasticity and used the Mann-Kendall test to assess significance of the interannual trend. For each phenology measure, we removed any record with more than three total missing values and then interpolated any gaps of one value. Any remaining records with missing data (gaps of two or more years) were removed.
As Sen's slope requires a continuous record, we included all AVHRR years in this analysis, but removed years with orbital degradation from the Mann-Kendall test of the significance of the trend. We used the land-cover data to examine the distribution of Sen's slopes between cover classes.
Second, to identify associations of climate factors in addition to time-based trends in phenological measures, we fit linear models including the explanatory variables of annual mean temperature, annual precipitation, and time. Our model for each AVHRR pixel was: where Y i is the ith observation of the phenological measure (dependent variable), B 0 is the intercept, B 1… B k are the slope coefficients for independent variables, Year i is the ith year of the observation, P i is the total annual precipitation in the ith year, T i is the annual mean temperature in the ith year, and ε i is the error term for the ith observation.
This allows us to expand beyond simple trend analyses to disentangle associations of time, precipitation, and temperature. We included time (Year) as a variable to represent temporally dependent linear effects that are not accounted for by temperature and precipitation measures, such as changes to land management or use, disturbance recovery, other climatic factors including increasing atmospheric CO 2 concentrations, or sensor effects (e.g., Auch et al., 2011;Brookshire et al., 2020;Nemani et al., 2003;Stoy et al., 2018;Zhu et al., 2016). To display the relative contributions of each predictor across space, we calculated their partial correlation. Then, we took the absolute value of these and normalized values to cover the range of 0-255 so each independent variable could be used as an input band into a RGB image.

| Phenology trends
Notable patterns emerged from the interannual trends of phenological measures (Figure 3; Table A1). For example, the majority of trends for the start and end of season increased, suggesting both are occurring later, but less than 15% of AVHRR pixels exhibited significant trends. The majority of interannual trends toward earlier start-of-season days were found in the southern portion of the study domain, whereas the northeastern portion of the study area was characterized by later end of season and some later days of peak productivity ( Figure 3). Season-long productivity had significant trends for 23% of AVHRR pixels in the study area, the majority of which (84%) increased. Similarly, 29% of pixels had significant trends for peak productivity of which 95% increased. The larger increases in season-long productivity and peak productivity were in the northeastern NWP, whereas the majority of decreases occurred in the southwestern portion of the study area ( Figure 3).
Significant trends in phenological measures varied by plant community type and responses were not synchronous, and trends are likely only associated with major land-cover changes in a small number of cases (Tables 1 and A1; Figures A3 and A4). Some notable exceptions between land-cover classes for phenological measures include shrub communities having a balanced split between trends of earlier (48%) and later (52%) spring start compared to other communities wherein the preponderance of significant AVHRR pixels (1) moved toward later spring onset (>70%). Evergreen and deciduous forests and barren areas tended toward an earlier end of season, while other communities tended toward later end-of-season dates.
Significant trends in peak productivity (95%) and season-long productivity (96%) largely increased; however, barren communities and to some extent evergreen forests differed from other communities with less dominant splits between increasing and decreasing peak and season-long productivity. Between 2001 and 2016, 94% of the study area (as measured in changes to Landsat pixels) remained in the same land-cover class ( Figure A2) and the vast majority (92%) of the AVHRR pixels in the study area had limited, less than 20%, subpixel change (Table 1). AVHRR pixels with significant trends in phenological measures were not more likely to have increased subpixel cover class changes ( Figure A4).

| Associations of climate and time to phenology measures
Coefficients from our multiple regression analyses of phenology (Equation 1) varied across space, between each phenology measure, and across community types. Annual values of precipitation and temperature ( Figure A1) in combination with time performed well in accounting for variation in peak productivity and season-long productivity but had limited explanatory power for date-based measures ( Figure A5). After accounting for annual climate, many areas maintained strong associations with time (a surrogate for other potential effects not covered by precipitation and temperature) and in similar areas as those with significant trends in the interannual trend analyses (Figures A6 and A7). In general, peak productivity and season-long productivity increased as annual precipitation increased (Figures 4, 5, A8, and A9). However, higher mean annual temperature tended to increase peak productivity yet decrease season-long productivity (Figures 4, 5, A10, and A11). Higher mean annual temperatures were predominantly associated with earlier start-of-season dates, but changes to end of season and length of season were much more mixed between earlier/shorter and later/ longer ( Figures A10 and A11). In analyses investigating the relative contributions of all predictors through a partial correlation analysis, spatial patterns emerged for each of the six phenology measures ( Figure A15). In the northeastern portion of the study area, time-a surrogate for multiple possible factors-was strongly associated with end-ofseason dates, while areas farther west tended to have stronger associations with temperature ( Figure A15). Temperature often had the largest partial correlation with date-based phenological measures. However, for productivity-based measures, precipitation had larger correlations with peak productivity, except in agricultural areas to the NE and forested areas in the SW (Figure 4).
For season-long productivity, temperature often in combination with annual precipitation had the larger partial correlations for most of the study area, with the same patterns for time as were seen for peak productivity ( Figure 5). Within land-cover types, the most prominent shifts in correlations due to precipitation for peak productivity versus precipitation and temperature for season-long productivity occurred in the shrub, pasture/crop, and wetland communities (Figures 4 and 5).   2001, 2004, 2006, 2008, 2011, and 2013. The "Other" land-cover type includes AVHRR pixels where the modal landcover type changed over the time period. Change classes based on the NLCD 2016 Land Cover Change index, as the proportion of 30-m pixels coded under any change category within each AVHRR pixel.

| D ISCUSS I ON
We identified three main findings from this research: (1) Over our study period, we found that only a small proportion of our study area had significant temporal trends in any of our six phenological measures; (2) annual climatic variables had robust explanatory power for productivity-based phenological measures; and (3)

F I G U R E 5
The combination of partial correlation coefficients for different drivers on season-long vegetation productivity. Scatter plots (left column and bottom row) and RGB image (upper right) illustrate the partial correlations of annual precipitation (AP), mean annual temperature (MAT), and time (Yr) on productivity as measured by AVHRR time-integrated NDVI from 1989 to 2014. The thick solid line represents the study area boundary, and the dashed lines are U.S. state boundaries. The scatter plots identify the absolute values of partial correlations from 2000 randomly selected pixels from the study for each vegetation community. In the map, color represents the relative association (absolute values) of the correlation of each factor. The brightness of pixels is relative to the combined partial correlations where darker colors have smaller and brighter colors have larger summed partial correlations. Hatching indicates areas where at least one partial correlation coefficient of the three variables considered was significant (p < .05) Second, annual precipitation and temperature had strong explanatory power for productivity-related phenological measures but tended to predict date-based measures poorly. By including relatively simple explanatory variables, we were able to identify changes through time and across communities and to identify divergent influences of temperature on peak versus season-long productivity. Therefore, the definition of greening and implications of both measures should be clearly understood and incorporated in phenological studies . Furthermore, illustrating the consequences of greening and/or browning under contemporary or future scenarios should include assessment of their patterns across varied plant communities and assess the trade-offs/exchanges of higher spring/peak productivity for lower season-long production (Butterfield et al., 2020;Hu et al., 2010;Wang et al., 2011).  (Gu et al., 2003;Hao et al., 2019;Piao et al., 2019;Turner, 2010;Zhang, Liu, et al., 2019), all of which complicate temperature-and precipitation-based explanations for phenological phenomena. Accordingly, specific drivers or surrogates such as time should be quantitatively included along with climate.

| The role of climate in phenological variability
Increasing precipitation had a clear association with increased productivity across large areas of the NWP (Figure 5), such as in low elevation, dry, warm areas like the Powder River basin ( Figure 6).
We found that grasslands and shrublands had larger positive responses in both peak and season-long productivity to precipitation relative to other communities ( Figure A13). Precipitation is a strong driver of grassland productivity in the NWP (Petrie et al., 2016;Yang et al., 1998), and in other similar biomes, globally (Meng et al., 2020;Wu et al., 2020). Furthermore, responses by a community to water limitations, and therefore productivity, vary between communities (Maurer et al., 2020;Ponce-Campos et al., 2013;Webb et al., 1978) and even between grassland types (Konings et al., 2017;Yang et al., 1998). Conversely increasing temperature mostly had a divergent association, increasing peak productivity but lowering season-long productivity. Implications of warming temperatures on productivity can vary between vegetation functional groups (Livensperger et al., 2016), have cascading ecosystem effects (Beard et al., 2019;Rehnus et al., 2020;Renner & Zohner, 2018), and include consequences to human society (Bezerra et al., 2019;Stevenson et al., 2015).
Overall, both precipitation and temperature were needed to fully explain phenological measures (e.g., Gao et al., 2020). Spatial and ecological differences in limiting factors-such as temperature restrictions on the length of the season in water-sufficient areas versus water-limited systems with ample growing days-lead to spatial variation Cowles et al., 2018;Piao et al., 2019;Reed et al., 2019;Yang et al., 2021). The precipitation and temperature controls on productivity are stronger in arid and semi-arid regions such as much of the NWP (see Figure A1), compared to more mesic areas (Maurer et al., 2020). There are important trade-offs that may come from shifts in the growing season (Figure 2; see also Fu et al., 2017;Butterfield et al., 2020). Specifically, whereas warm, wet springs may increase early-season productivity and peak productivity, this may be offset by lower soil-water recharge and decreased production during dry, hot summers (Butterfield et al., 2020;Hu et al., 2010;Lian et al., 2020;Wang et al., 2011). Conversely, a later start of spring from delayed snow melt can lead to increased, sustained, summer productivity (Hu et al., 2010;Potter, 2020). These relationships and potential trade-offs may impose conflicting selection pressures on other members of the ecosystem due to asynchronous changes between start of season and peak versus season-long productivity.
Further investigations of the role of vapor pressure deficit (VPD), or other metrics of atmospheric moisture and water balance, may help explain seasonal vegetative processes, beyond those informed by temperature and precipitation. In North American grasslands, VPD during July and August was a strong predictor of productivity (Konings et al., 2017). Higher VPD leads to decreased gross primary productivity across global ecosystems (Zhang, Ficklin, et al., 2019), and VPD even was a strong predictor in small mammal distribution models Smith et al., 2019). In our study area, northeastern Montana and much of North and South Dakota showed decreasing VPD between 1979 and 2013 in the spring and summer, while VPD in central Montana and Wyoming increased (Bromley et al., 2020;Ficklin & Novick, 2017). Although VPD is a promising explanatory variable for vegetation responses to climate, it is unclear how to best aggregate VPD across time to account for its nonlinear impacts on productivity (Lasslop et al., 2010).

| Nonclimate associations
We included time in our models as a surrogate for multiple potential linear factors other than temperature and precipitation that can influence phenological measures, and in the NWP time strongly governed response in specific areas (Figures 4-6). Fire and insect disturbances are common in forests of the NWP and may have affected productivity (Pfeifer et al., 2011;Turner, 2010). For example, the Greater Yellowstone Ecosystem experienced stand-replacing fires across large areas in 1988 that were followed by reestablishment of forest communities such as those dominate by lodgepole pine (Pinus contorta) (Turner, 2010). We found stronger associations between time and season-long productivity within fire perimeters compared to outside, illustrating postdisturbance succession such as from patchy pioneer plant species to dense early-stage evergreen forests ( Figure 6).

| Opportunities for application and addition study
We used annual values for temperature and precipitation to test the improvement from interannual trends with only a small increase to the dimensionality of our models. Although this approach had strong explanatory power for productivity measures, it was more limited for date-based phenological measures. Phenology is responsive to seasonal conditions prior to and during the growing season that may diverge from annual patterns (e.g., Bianchi et al., 2019;Fu et al., 2017;Ren et al., 2020), and although there is often strong correlation between annual and seasonal values, climate trends and species response to climate may not vary homogeneously (Bromley et al., 2020). Temperature was the clearer control on season duration: Warmer temperatures can advance spring (Allstadt et al., 2015;Fu et al., 2017;Ren et al., 2020) and extend the growing season through longer snow-free periods and later freeze events (Potter, 2020).
Later spring starts and spatial differences in the direction of spring movement are not unprecedented in the Northern Hemisphere (e.g., Li et al., 2019;Ren et al., 2020). For example, spring is trending earlier in western and central Europe but later in eastern Europe, attributable to differing weather patterns (Ahas et al., 2002). Shifts to later end-of-season dates (as found in much of our study area) are found in central Asia and China, and movement toward earlier end-of-season dates is found in more southern areas of the Great Plains, western and central Asia, and parts of Mongolia (Ren et al., 2020). Overall, more complex models addressing seasonal time periods (Ren et al., 2020), additive growing degree days (de Beurs & Henebry, 2007), and/or ecological memory  have promise for date-based phenological measures.
In our study, we identified that trends toward higher peak productivity did not always translate to higher season-long production. Drivers of peak and season-long productivity can differ in general and across vegetation communities (Fu et al., 2017(Fu et al., , 2019Xia et al., 2015). Mixed responses of greening versus browning in grassland, shrub, steppe, and forested areas are also found in regions of Patagonia (Bianchi et al., 2019) and across vegetation communities in Mongolia (Meng et al., 2020). We demonstrated that greening and browning studies should consider peak and season-long productivity differently, given that their relationships to climate drivers differed markedly. Earlier or higher growth and/or longer growing-season peaks can be offset by limited production in a dry summer (Livensperger et al., 2016;e.g., Hu et al., 2010). These shifts, trade-offs, and asynchronies in vegetation phenology have feedbacks and consequences to ecological communities, including consumer-resource dynamics, species' behaviors, and other ecosystem processes (Beard et al., 2019;Morisette et al., 2009).
Understanding the implications of future climate change is important for predicting potential ecological and social-ecological impacts (e.g., Allstadt et al., 2015;Beard et al., 2019;Beever et al., 2019;Epstein et al., 2021;Hufkens et al., 2016). A collective understanding and consideration of multiple phenological measures and their drivers is important for effective management actions and continued development of a multiscale ecological perspective for management (Bezerra et al., 2019;Morisette et al., 2009).

CO N FLI C T O F I NTE R E S T
None declared.  Note: 1-Increasing and decreasing values are calculated from only significant pixels, with total trends being the number of significant pixels divided by the total number pixels of the community type. For each community type, the percentage increasing (decreasing) is the number of significant pixels with a positive (negative) slope divided by the total number of significant pixels. Therefore, increasing slopes correspond to transition dates later in the year, longer seasons, or higher productivity as applicable to each phenology measure.

A D D ITI O N A L S PATI OTE M P O R A L A N D CO M M U N IT Y PAT TE R N S A N D FI N D I N G S
These figures and tables provide additional information on temperature and precipitation patterns, potential contributions of land-cover change, differences between communities, and distributions of slopes from the multiple linear regression analyses. In addition, results for all phenological measures are included from the partial correlation analyses showing the spatial differences between predictors.

F I G U R E A 4
Stacked proportional bar plot comparison of study area wide change of land-cover types (ForestDecid = Deciduous Forest, ForestEverg = Evergreen Forest, see text for description of each land-cover class) within AVHRR pixels to those with significant interannual trend of peak productivity (MAXN). Land cover derived from the Landsat-based (30 m) NLCD Change Index record. MAXN used as the example comparison as it has the largest area of significant interannual trend and overlap with the majority of other significant interannual trends for the other five phenological measures studied F I G U R E A 5 R 2 values for multiple linear regression models of the change from mean annual temperature, annual precipitation, and time for phenological measures in the Northwestern Plains. The study area is depicted with a dashed black line and U.S. states with gray solid lines F I G U R E A 6 Slopes for change through time for phenological measures calculated from the AVHRR satellite record in the Northwestern Plains after accounting for annual precipitation and mean annual temperature. The study area is depicted with a dashed black line and U.S. states with gray solid lines. Note different axes for each plot F I G U R E A 7 Histograms for the change over time for phenological measures calculated from the AVHRR satellite record for the Northwestern Plains after accounting for annual precipitation and mean annual temperature. Note different axes for each plot F I G U R E A 8 Slopes for the effect of annual precipitation on phenological measures calculated from the AVHRR satellite record in the Northwestern Plains after accounting for mean annual temperature and time. The study area is depicted with a dashed black line and U.S. states with gray solid lines. Note different axes for each plot

F I G U R E A 9
Histograms of the change from annual precipitation for phenological measures calculated from the AVHRR satellite record for Northwestern Plains communities after accounting for time and mean annual temperature. Note different axes for each plot F I G U R E A 1 0 Slopes for the effect of mean annual temperature for phenological measures calculated from the AVHRR satellite record in the Northwestern Plains after accounting for annual precipitation and remaining time-based trend. The study area is depicted with a dashed black line and U.S. states with gray solid lines. Note different axes for each plot

F I G U R E A 11
Histograms of the change due to annual mean temperature on phenological calculated from the AVHRR satellite record measures for Northwestern Plains communities after accounting for time and annual precipitation. Note different axes for each plot F I G U R E A 1 2 Boxplots for the change over time in phenological measures calculated from the AVHRR satellite record for major Northwestern Plains vegetation types (ForestDecid = Deciduous Forest, ForestEverg = Evergreen Forest, see text for description of each land-cover class) after accounting for annual precipitation and mean annual temperature. Note different axes for each plot F I G U R E A 1 3 Boxplots for the change from annual precipitation on phenological measures calculated from the AVHRR satellite record for major Northwestern Plains vegetation types (ForestDecid = Deciduous Forest, ForestEverg = Evergreen Forest, see text for description of each land-cover class) after accounting for time and mean annual temperature. Note different axes for each plot F I G U R E A 1 4 Boxplots for the change from mean annual temperature on phenological measures calculated from the AVHRR satellite record for major Northwestern Plains vegetation types (ForestDecid = Deciduous Forest, ForestEverg = Evergreen Forest, see text for description of each land-cover class) after accounting for time and annual precipitation. Note different axes for each plot F I G U R E A 1 5 Combined partial correlations for time (Yr), annual precipitation (AP), and mean annual temperature (MAT) in phenological measures for the Northwestern Plains. Color represents the relative correlation of each factor in explaining the measures between 1989 and 2014. Brightness of the color is relative to the absolute value of the combined correlations where darker colors are smaller and brighter colors larger combined partial correlations. Dashed lines are U.S. state boundaries