Climate seasonality and extremes influence net primary productivity across California’s grasslands, shrublands, and woodlands

Terrestrial vegetation is a substantial carbon sink and plays a foundational role in regional and global climate change mitigation strategies. The state of California, USA, commits to achieving carbon neutrality by 2045 in part by managing terrestrial ecosystems to sequester more than 80 MMT of CO2. We used a 35-year net primary productivity (NPP) remote sensing product with gridded climate, soil, topography, and vegetation data to evaluate spatiotemporal drivers of NPP variation and identify drivers of NPP response to extremes in water availability in California’s major grasslands, shrublands, and woodlands. We used generalized boosted models (GBMs) and linear mixed effects models (LMMs) to identify influential predictors of NPP and characterize their relationships with NPP across seven major vegetation cover types: annual grasslands, blue oak, chamise-redshank chaparral, coastal scrub, coastal oak woodland, mixed chaparral, and montane hardwood. Climate seasonality, specifically greater precipitation and warmer minimum temperatures in early spring and winter, was associated with greater NPP across space, particularly in chaparral, blue oak, and grassland systems. Maximum annual temperature and climatic water deficit (CWD) showed a negative relationship with NPP in most vegetation cover types, particularly chaparral and coastal scrub. We found a significant decrease in NPP over time in most vegetation types, appearing to coincide with the 2012–2016 California mega-drought. However, response to water availability extremes differed by vegetation type. In most vegetation types, especially grasslands, increases in NPP in extreme wet years were greater than declines in NPP in dry years. Our analysis characterizes several climate risks and conservation opportunities in using California’s natural lands to store carbon. Namely, shifts in climate seasonality and water availability extremes threaten these systems’ ability to fix carbon, yet hotspots of NPP resilience may exist and could be enhanced through conservation and restoration. Additional mechanistic work can help illuminate these opportunities and prioritize conservation decision making.


Introduction
Terrestrial vegetation serves as a substantial carbon sink, absorbing more than 30% of anthropogenic CO 2 emissions over the last several decades (Ciais et al 2014, Walker et al 2021. Terrestrial ecosystem conservation and management, therefore, plays a fundamental role in global and regional strategies to mitigate the rate of climate change (Nabuurs et al 2008, Conant et al 2017, California Air Resources Board 2019. For example, legislation in California, USA, commits the state to reducing greenhouse gas emissions 40% below 1990 levels by 2030 and achieving carbon neutrality by 2045 (California Air Resources Board 2017). A central approach by the State to meet this goal is through doubling investments in land conservation and management, with the aim of using these investments to sequester more than 80 million metric tons of CO 2 by 2045 (California Air Resources Board 2019). However, California's diverse climate and topography yield highly variable patterns of vegetation productivity (Larsen et al 2015, Becchetti et al 2016, Carey et al 2020, posing a challenge for the State in understanding how to prioritize land management and conservation investments. A greater understating of the spatiotemporal variation in terrestrial ecosystem productivity and the environmental drivers of this variation are critical if we are to enhance conservation decision making and optimize the potential for terrestrial vegetation to store carbon. For example, knowledge of how ecosystem productivity varies among ecosystem type or along broad environmental gradients could factor into decisions about where land conservation or restoration investments are made. Understanding the potential for vegetation to store carbon fundamentally depends on understanding drivers of net primary productivity (NPP; g C m −2 yr −1 ) (Jackson and Prince 2016, Dangal et al 2020, which is a measure of carbon accumulated by plants through photosynthesis per unit area and unit time. In arid and semi-arid ecosystems in California, spatial variation in terrestrial vegetation NPP is strongly driven by climatic gradients in mean annual precipitation and temperature, with greater precipitation and warmer temperatures associated with higher NPP (Lieth 1973, Lehouerou et al 1988. However, climate seasonality could play an equally large role in driving spatiotemporal variation in NPP. For example, previous work in Eurasian temperate grasslands found that seasonal distribution of precipitation explained as much spatial variation in grassland NPP as mean annual precipitation, but the importance of precipitation seasonality diminished in grasslands with very low and very high mean annual precipitation (Guo et al 2012). Likewise, work on California grasslands found that both mean annual and monthly distribution of precipitation were primary drivers of the spatial distribution of NPP (Liu et al 2021). Understanding how climate seasonality impacts spatiotemporal variation in NPP in California is important as seasonal patterns of precipitation and temperature are expected to shift under a changing climate, including relatively wetter and warmer winter conditions and hotter and drier spring conditions (Pierce et al 2018).
Natural and anthropogenically driven changes in climate extremes may further contribute to spatiotemporal variation in NPP. For example, El Niño-Southern Oscillation, a major global driver of interannual climate variability significantly impacts NPP in the western United States (Dannenberg et al 2015) with the much wetter and warmer El Niño events providing more favorable conditions for NPP compared to the drier and colder La Niña events. In mid and high latitude systems, including the western United States, increasing temperatures associated with anthropogenic warming have increased annual NPP about 1.5% across much of the western United States (Nemani et al 2003). Continued warming and higher maximum temperatures, however, have increased the frequency and duration of drought across the globe, reducing global annual NPP about 1% (Zhao and Running 2010). How precipitation and temperatures extremes impact NPP in California is poorly understood. Warming temperatures have contributed to record breaking mega-droughts (Diffenbaugh et al 2015, Polade et al 2017 which negatively affect NPP. However, in some situations, warming temperatures could partly offset drought induced reductions in NPP by stimulating higher NPP during brief periods of water availability in otherwise dry years. Understanding how extreme climatic conditions impact spatiotemporal variation in NPP in California is important as future climate projections across California predict a greater frequency of climatic extremes, including a greater frequency of very dry and very wet years and higher maximum temperatures (Pierce et al 2018, Swain et al 2018.
Not all vegetation types in California may respond similarly to changes in climate means, seasonality, or extremes. California terrestrial vegetation represents a global biodiversity hotspot (Cowling et al 1996, Kt and Jetz 2007, Rundel et al 2016, and understanding drivers of NPP variation within and among a diverse range of vegetation cover types is essential for allowing the State to meet its carbon emissions reduction and carbon neutrality goals. Much of what we know about the potential for ecosystems to store carbon in California is focused on forests due to their high NPP and large standing biomass (Pan et al 2011, Gonzalez et al 2015. However, California contains over 17 million hectares of grasslands, woodlands, and shrublands that may show greater resilience to a drying climate compared to forests and serve as important carbon sinks (Dass et al 2018, Mayer andSilver 2022). Grasslands, woodlands, and shrublands in California span broad gradients of topography, soil, and climate conditions and show high spatial and temporal variation in productivity (Murphy 1970, Larsen et al 2015, Becchetti et al 2016. Of these vegetation types, grassland NPP dynamics are the most well studied due to their large extent and important role in forage production for California's livestock industry. California's grasslands consist primarily of introduced annual grasses where production is heavily influenced by precipitation amount and timing as well as minimum air temperatures (Bartolome et al 2007, Becchetti et al 2016, Liu et al 2021. Namely, high precipitation falling in the relatively warmer fall and spring enhances productivity while high precipitation in the relatively cooler winter has little effect on productivity (Becchetti et al 2016, Liu et al 2021. While we expect similar patterns of seasonality to be important drivers of NPP in a range of California's grasslands, woodlands, and shrublands, the relative importance and effect size of these factors for most vegetation types is unknown. To prioritize investments in terrestrial ecosystem conservation to sequester carbon, it is essential we understand how these drivers may differ among terrestrial ecosystem types.
This study had two objectives. First, we used a 35 year remotely sensed NPP product with gridded climate, soil, topography, and vegetation composition data to evaluate spatiotemporal drivers of NPP variation within and among California's major shrubland, grassland and woodland vegetation types. Second, we combined the same datasets with a multi-scalar drought index to evaluate drivers of NPP response to extreme variation in water availability within and among vegetation types. Collectively, the aim of this work is to provide greater insight into NPP drivers for a range of terrestrial ecosystems across California which may ultimately help enhance how land conservation and restoration can be used to meet the State's climate goals.

Study area
Our study focused on seven predominate rangeland vegetation cover types in California: annual grasslands, blue oak woodland and blue oak foothillpine (combined into one blue oak cover type), chamise-redshank chaparral, coastal oak woodland, coastal scrub, mixed chaparral, and montane hardwood, which collectively comprise over 14 million hectares, nearly half of the State's land area ( figure 1(a)). Intensive land use cover types, namely heavy agriculture (for example, pasture, cropland, and irrigated hayfield) and large urban areas, were masked and not included in our analysis. These lands appear as white pixels located primarily in the center of our study area (figure 1). Vegetation cover types are based on CALFIRE's FRAP vegetation (fveg15_1) dataset (Department of Foresty and Fire Protection 2015) within Major Land Resource Areas 14, 15, 17, 18, 19, and 20 (United States Department of Agriculture 2006). Our study area, extending from 32 • to 41 • N and 116 • to 123 • W, contains substantial variation in topography, ranging from −27 to 3029 m above sea level ( figure 1(a)). The region's Mediterranean climate with mild, wet winters and hot, dry summers (Becchetti et al 2016) results in a wide range of temperatures and precipitation amounts across our study region, including an annual minimum temperature range of −10 to 13 • C, annual maximum temperature range of 18 • C-40 • C (figure 1(c)), and mean annual precipitation range of 151-1624 mm ( figure 1(b)). Maximum precipitation and minimum temperatures typically occur in winter months (December-February). Spring (March-May) typically shows lower precipitation and higher temperatures than winter while maximum annual temperatures and lowest water availability occur in summer (June-September). Interannual variation in climate extremes are driven by complex local and global changes in energy budgets including effects of sea surface temperatures on atmospheric pressure over the northeastern Pacific Ocean (Swain et al 2017).

Response variables: NPP and NPP asymmetry
Our objectives were to understand the drivers of spatial and temporal variation in annual NPP and the drivers associated with spatial variation in NPP response to interannual water availability extremes. NPP measures the amount of carbon absorbed during photosynthesis minus carbon lost in respiration. While NPP is not a direct measure of longterm carbon storage, it is a primary factor involved in biomass accumulation and carbon sequestration, particularly in semi-arid ecosystems (Whittaker andLikens 1973, Keeling andPhillips 2007). We used annual cumulative NPP data from the Rangeland Analysis Platform (RAP); (Robinson et al 2019), available on Google Earth Engine (GEE; ee.ImageCollection('projects/rangeland-analysis-pla tform/npp-partitioned-v2')), as our estimate of vegetation NPP. RAP combines Landsat imagery at a 30 m pixel resolution with the MOD17 NPP algorithm to provide high quality estimates of annual NPP for the conterminous United States from 1986 through 2020. We calculated a second response variable, NPP asymmetry, to quantify vegetation NPP response to interannual water availability extremes. To do this, we calculated an asymmetry index (NPP AI) where NPP asymmetry is the difference between NPP gains during extreme wet years (Pulses, sensu Wu et al 2018) and NPP losses during extreme dry years (Declines, sensu Wu et al 2018). We used the 90th and 10th percentile of standardized precipitation evapotranspiration index (SPEI) (Vicente-Serrano et al 2010) to define extreme wet and dry years, respectively, at each location across our study area. We obtained SPEI data from gridMET (Abatzoglou 2013) on GEE (ee.ImageCollection('GRIDMET/ DROUGHT')). NPP Pulse, NPP Decline, and NPP AI were computed as: Figure 1. Spatial distribution of (a) elevation, (b) mean annual precipitation (mm), (c) and mean annual maximum temperature ( • C) across study area in California grasslands, shrublands, and woodlands. Large white areas in the center of maps are vegetation types not included in our analysis, primarily consisting of agriculture and urban areas.
where med(NPP) is the median NPP for the full measurement record (years 1986-2020), med(NPP p90 ) is the median NPP for wet years, defined as years that fell at or above the 90th percentile of SPEI, and med(NPP p10 ) is the median NPP for dry years, defined as years that fell at or below the 10th percentile of SPEI. High values of NPP pulse indicate a strong positive response of NPP to high water extremes, and high values of NPP Decline indicate a strong negative response of NPP to low water extremes. Therefore, positive values of NPP AI correspond to greater gains in productivity in wet years than losses in dry years, and negative values of NPP AI indicate greater losses in productivity during dry years than gains in wet years. to generate a list of potential abiotic and biotic variables to include in our models of spatial and temporal variation in NPP and spatial variation in NPP asymmetry (table 1). To assess the effects of climate over different temporal scales, we aggregated monthly precipitation, minimum and maximum temperature, and CWD from the basin characterization model (BCM) version 8 (Flint et al 2021) into seasonal and annual timesteps. We calculated long-term means and standard deviations for annual climate variables and means for seasonal climate variables. In California's Mediterranean climate, vegetation is adapted to seasonal summer dry periods (June-September) and most vegetation either severely reduces physiological activity, is summer dormant, or senesces during this dry period. As a consequence, NPP is minimal during the summer dry period (Xu and Baldocchi 2004). Because of this, seasonal timesteps were fall (October-November), representing the onset of the growing season with the first rains; winter (December-January), where water is abundant, but temperatures are low; early spring (February-March), which sees high precipitation and high temperatures; and mid-spring (April-May), where water is less abundant, but temperatures are high. To capture the effects of precipitation deviance from average conditions, we calculated annual precipitation variability following (Al-Yaari et al 2020). Lastly, we included average solar radiation from the BCM to represent long-term solar energy inputs.

Predictor variables
To assess the effects of physical terrain and geography on NPP and NPP asymmetry, our initial predictor list also included elevation, southness, Topographic Wetness Index (TWI; Beven and Kirkby 1979), and four soil variables. Elevation, southness, and TWI were calculated from a digital elevation model from the Shuttle Radar Topography Mission (Farr et al 2007) dataset on GEE (ee.Image('USGS/SRTMGL1_003')) at a 30 m resolution. Southness was calculated following (Ackerly et al 2020) and represents the amount of southfacing exposure at a site, ranging from −1 for a vertical north slope to 1 for a vertical south slope. TWI was calculated as upslope contributing area scaled by slope and represents topographically controlled soil wetness. Soil variables included percent silt, percent sand, percent clay, and pH, acquired from the National Resource Conservation Service's SSURGO  (2015) * Long term (1986-2020) mean was calculated and included as a predictor. * * Long term (1986-2020) mean and standard deviation were calculated and included as a predictors. and STATSGO national soil databases (Soil Survey Staff 2017). To examine impacts of broad scale geography, we included latitude and longitude.
We included vegetation type and functional group evenness as predictors of NPP to examine how vegetation structure relates to NPP and NPP asymmetry. We calculated functional group evenness as Pielou's evenness index (Pielou 1966) using percent cover of functional groups (herbs, shrubs, and trees) data acquired at 30 m resolution from RAP on GEE (ee.ImageCollection('projects/rangeland-analysisplatform/vegetation-cover-v2')). Functional group evenness quantifies the number of plant functional groups present (e.g. shrubs, grasses, trees) and the relative distribution of cover among functional groups. Sites with a higher number of functional groups and more even allocation of cover among groups would receive a higher evenness score than a site with the same number of functional groups but unequal distribution of cover among groups. Evenness of functional groups is a key component of functional diversity, where greater evenness may represent healthier ecosystem functioning and more complete utilization of available resources (Mason et al 2005). Species identity or proportion of native vs. nonnative species are not available at a statewide scale so were not included in this analysis.
Lastly, to control for disturbance effects on NPP, we included fire return interval in early models.
However, this variable was not included in final modeling efforts as it was uninformative.

Statistical modeling
We implemented GBMs to identify relevant predictors for NPP and NPP asymmetry from our initial list of candidate predictors (table 1). GBMs are a machine learning algorithm that combine the statistical power of regression trees (models that predict a response variable by making recursive binary splits in predictor variables) and boosting (a method of combining models to reduce bias and improve accuracy) (Elith et al 2008). GBMs capture the effects of many predictors and their interactions on complex processes and are less sensitive to spatial autocorrelation and missing values than generalized linear models, making them well suited for ecological applications (Elith et al 2008, Crase et al 2012. GBMs give an estimate of the relative influence of each predictor on the response variable, calculated as the percent reduction of squared error attributable to each variable in the final model (Greenwell et al 2020). From the GBM models, we identified a subset of predictors with high relative influence and utilized LMMs, built for each vegetation type, to characterize the nature of the relationships between this subset of predictors and temporal and spatial variation in NPP. For NPP asymmetry, our goal was to explore a range of potential variables and identify those with the greatest influence on NPP asymmetry across and within rangeland vegetation types in California. We did so by training GBMs with all candidate predictors (table 1) on the whole study area first, then training GBMs on each vegetation type individually.
All GBMs were built in R (Team 2022), Version 4.1.2, using package 'gbm' (Greenwell et al 2020). In each model, we used 10 000 trees, a shrinkage rate of 0.1, and five-fold cross validation. We limited the interaction depth to 2 to ensure results would be interpretable and applicable in the context of land management. GBMs were run on data at 270 m and 4 km resolution, allowing us to assess scaledependency of our results. We used nearest neighbor averaging to upscale spatial data where necessary. We trained models of NPP and NPP asymmetry across the entire study area on a random 20% of our full dataset, stratified evenly across the seven vegetation types.
Using the top predictors of NPP we identified with GBMs, we built LMMs to assess the relationships between annual NPP and environmental predictors. We then built models for each vegetation type with a subset of environmental predictors as fixed effects and location as the only random effect. Random intercepts were allowed, but slopes were fixed to reduce overfitting. Working from the list of variables with high relative influence from the GBMs, we first included spatial predictors to capture variation across space, then added temporal predictors until Akaike information criterion (AIC) stopped decreasing and log-likelihood stopped increasing by significant amounts. All spatial variables represented longterm means, and temporal variables were group mean centered to partition the spatial and temporal variation in our predictors. Lastly, interaction and quadratic terms were incorporated to establish a final model with an optimized AIC, log likelihood, and amount of variation explained for each vegetation type. Highly correlated and insignificant parameters were excluded (table S1). This process was repeated with data at the 270 m and 4 km scale.

Spatial patterns in NPP across the study region
Mean annual NPP (371.4 ± 176.5 g C m −2 yr −1 , mean ± SD) across the study region varied substantially among and within the major vegetation cover types (figures 2(a)-(d)). Highest mean annual NPP occurred in coastal oak woodlands (646.8 ± 249.0) and montane hardwood (576.9 ± 208.7) cover types while annual grasslands (267.6 ± 86.2) showed the lowest annual NPP. In general, mean annual NPP was highest on the west side of the coast range and transverse range and lowest in southern portions of the San Joaquin Valley. Our GBMs built for the entire study area suggested spatial variables related to vegetation structure and climate had the largest relative influence on mean annual NPP, with the most influential predictors being vegetation cover type (30% relative influence), mean maximum annual temperature (11%), mean functional group evenness (11%), mean early spring precipitation (7%), and southness (5%). Mean values for fall precipitation (4%), winter precipitation (3%), and maximum mid-spring temperature (3%) also emerged as top influences on NPP. In total, spatial predictors reduced the squared error in our GBM by 92.2%. Our analyses did not show strong differences between the 270 m and 4 km scale, so all values reported are from the 270 m analysis.

Drivers of NPP spatial variation within cover types
As expected, spatial variation in NPP in all cover types was strongly associated with spatial variation in temperature and precipitation (figure 3). Our LMM results, however, suggested the nature of the relationship of NPP with temperature and precipitation differed significantly across vegetation type. For, example, variation in early spring precipitation emerged as the largest factor associated with variation in in NPP in five of our seven cover types, but the positive relationship of early spring precipitation and NPP varied substantially and was notably greater in chamiseredshank chaparral than montane hardwood. Variation in coastal scrub and coastal oak woodland NPP appeared to be more strongly associated with variation in winter precipitation than the other systems. NPP of different systems also showed different relationships to temperature seasonality. In six of our seven cover types, greater mean annual maximum temperature was associated with a decrease in NPP, but in three systems (coastal scrub, mixed chaparral and montane hardwood) warmer early spring and winter minimum temperatures were associated with higher NPP. Overall, variation in predictors included in our LMMs for each vegetation type were associated with between 24.9% and 48.0% of NPP spatial variation. The model built for coastal scrub, a vegetation type characterized by low mean NPP and low mean precipitation, explained slightly more variation than any other model. In contrast, the model for montane hardwood, a vegetation type with high mean NPP and high mean precipitation, explained less variation than models for any other vegetation type. Due to spatial autocorrelation, LMMs tended to fit points with average conditions better than those with more extreme conditions (figure S1).

Drivers of NPP temporal variation within cover types
Temporal variation in climate predictors showed weak collinearity relative to long-term average climate, allowing a greater range of temporally varying predictors to be included in our models compared to our spatial predictors. Temporal variation in NPP within cover types was associated with complex interactions among our predictors (figures 4 and 5). Notably, NPP for all cover types showed a strong quadratic effect of year, with NPP increasing toward a maximum near the center of the study period then showing a negative slope in more recent years. In addition, in six of our seven vegetation types, the year effect on NPP interacted with climate predictors, and these interactions differed among vegetation types (figure 4). In the coastal scrub cover type, for example, NPP is maximized in Figure 3. Comparison of standardized parameter coefficients for spatial predictors (long term means) included in linear mixed effects models of net primary productivity (NPP; g C m −2 yr −1 ) for rangeland vegetation types in California. Standardized coefficients are calculated as the product of the unstandardized coefficient and the ratio between one standard deviation of NPP and one standard deviation of the spatial predictor. Error bars represent 95% confidence intervals. more recent years at below average fall precipitation amounts and in more distant years at above average fall precipitation amounts. In annual grasslands, high NPP was associated with lower-than-average winter minimum temperatures in more distant years and higher-than-average winter minimum temperatures in more recent years. In addition to these interactions with time, there also were several strong and reoccurring main effects of climate predictors associated with NPP. Namely, in five of our seven cover types, temporal variation in NPP was strongly associated with temporal variation in mid spring CWD, with coastal oak woodland temporal variation in NPP showing the lowest sensitivity to CWD and coastal scrub and chamise-redshank chaparral temporal variation in NPP showing comparable and high sensitivity to CWD (figure 5). Lastly, variation in NPP across all seven cover types was associated with quadratic effects of one or more climate variables, although the climate predictors, direction and concavity of the relationship differed across vegetation types (figure 5). Temporal variation in NPP explained in our models ranged between 4.65 and 15.1%.

Variation in NPP asymmetry across the study region
NPP asymmetry was positive over 61.8% of the study area but varied widely among and within vegetation types ( figure 6(a)). Mean NPP asymmetry was positive in all but one cover type, with annual grasslands showing the highest NPP asymmetry characterized by very high NPP pulse and high NPP decline. In contrast, coastal oak woodland and montane hardwood cover types showed positive NPP asymmetry characterized by low NPP pulse and very low NPP decline. Coastal scrub was the only cover type to show a negative NPP asymmetry with large NPP pulses in high water availability years exceeded by large NPP declines during low water availability years (figure 6(b)). We developed GBM models for each cover type to identify predictors influencing spatial variation in NPP asymmetry within each cover type. While some unique sets of predictors emerged within each cover type (table S2), five predictors (latitude, longitude, standard deviation of annual maximum temperature, mean functional group evenness, and standard deviation of function group evenness) consistently emerged within the top 12 for nearly all cover types (the one exception being longitude in the annual grassland model) (figure 7). Latitude and longitude were the two most influential predictors of NPP asymmetry in the blue oak, montane hardwood, and chamise-redshank chaparral cover types, while in all other cover types, standard deviation of functional group evenness was one of the top two predictors of NPP asymmetry. Longitude was a highly influential predictor of NPP asymmetry for every vegetation type except annual grasslands, which also has the largest geographic distribution of all vegetation types.  . Standardized parameter coefficients of highest order temporal predictors included in linear mixed effects models of NPP (g C m −2 yr −1 ) for rangeland vegetation types in California, excluding year parameters. Temporal predictors are group-mean centered, representing yearly divergence from long-term mean conditions at each point from 1986-2020 across the study area. Standardized coefficients are calculated as the product of the unstandardized coefficient and the ratio between one standard deviation of NPP and one standard deviation of the temporal predictor. Error bars represent 95% confidence intervals. temperature, and late-season maximum temperature. Importantly, though, we found that these relationships differ by cover type and not all systems respond strongly to the same seasonal climate patterns. For example, three vegetation types showed strong positive relationships between NPP and winter or early spring minimum temperature in addition to negative relationships with maximum temperature. Peng et al (2013) and Ma et al (2022) observed similar contrasting effects of minimum and maximum temperature on NPP in temperate grasslands. While our correlative approach cannot establish causal relationships, our results do suggest a possibility of increases in minimum winter temperatures due to climate change offsetting some reductions in NPP associated with a drying climate (Peng et al 2013, Ma et al 2022. Furthermore, our results suggest that future variation in precipitation seasonality could amplify overall negative effects of a drying climate on NPP with decreases in spring or winter precipitation being associated with significant NPP declines in all cover types included in our study. Importantly, subtle shifts in winter precipitation could have disproportionately larger effects on coastal oak woodland and coastal scrub ecosystems than other cover types. As our collective ability to forecast future changes in climate seasonality continues to advance (Doblas-Reyes et al 2013, Jones 2013, Urban 2019), these findings may provide some indication of how shifts in seasonal climate may impact the ability of grasslands, shrublands, and woodlands to help meet California's climate goals.
Across all vegetation cover types, NPP sensitivity to temporal variation in climate was lower than NPP sensitivity to spatial variation in climate. Standardized parameter coefficients, indicating strength of relationship with NPP, generally fell within −0.1 and 0.08 for temporal predictors and −0.35 and 0.45 for spatial predictors. The slope of the temporal relationship between climate and NPP is often much weaker than the slope of the spatial relationship between NPP and climate with one likely mechanism being interannual vegetation constraints limiting the magnitude of NPP response to temporal variation in climate (Lauenroth and Sala 1992, Sala et al 2012, Estiarte et al 2016, Knapp et al 2017. Nevertheless, several predictors not observed in the spatial analysis of NPP emerged as important predictors of NPP temporal variation. For example, an increase in annual or midspring CWD was associated with a temporal decrease in NPP in all cover types. Previous work has proposed legacy or lag effects of precipitation from previous years as a significant influence on NPP and a possible reason for discrepancies between spatialtemporal variation in climate and NPP (Bai et al 2008, Sala et al 2012, Knapp et al 2017. Our results suggest legacy effects of precipitation on current year water Figure 6. (a) Distribution of asymmetric response of net primary productivity (NPP; g C m −2 yr −1 ) to water availability anomalies across California rangelands, with outliers binned at 0.75 and −0.75, (b) spatial distribution of NPP asymmetry index (NPP AI), and (c) comparison of NPP AI, NPP pulse in wet years, and NPP decline in dry years among rangeland cover types in California. The red line in the histogram represents the overall mean NPP asymmetry index for the study area. Regions on the map marked 'x' and 'y' represent areas of potential NPP resilience to climate extremes. Large white areas in the center of the map are vegetation types not included in our analysis, primarily consisting of agriculture and urban areas. Error bars represent 95% confidence intervals. All mean values for NPP AI, NPP pulse, and NPP decline are significantly different from zero at the p < 0.05 level. NPP pulse and NPP decline for each vegetation type are significantly different from every other vegetation type at the p < 0.05 level. NPP AI for each vegetation type is significantly different from vegetation types marked with a different letter at the p < 0.05 level. NPP AI in vegetation types marked with the same letter are not significantly different at the p < 0.05 level. All p-values were adjusted with a Bonferroni adjustment for multiple comparisons. budget may be particularly important drivers of temporal variation of NPP in our shrub dominated systems as our chamise and scrub dominated systems showed the largest temporal reduction in NPP with an increasing CWD.
Across most vegetation types, we found a concave-down quadratic trend of annual NPP over time, reaching a maximum near the center of the study period (2005)(2006)(2007)(2008)(2009)(2010) then decreasing in recent years and showing complex interactions with climate. In line with previous work showing declines in NPP due to high temperatures and low precipitation (Chaplin-Kramer andGeorge 2013, Mu et al 2013), negative NPP trends in recent years align with a prolonged drought in California from 2012 to 2016 (Griffin and Anchukaitis 2014, Swain et al 2014, Mann and Gleick 2015. These findings provide insight into how increased frequency of extreme droughts may decrease vegetation productivity in California in the future. However, we also found evidence of changing relationships between NPP and climate over time. In several cases these interactions were intuitive. For example, in annual grasslands warmer winter minimum temperatures appeared to dampen the general decline in NPP through time. This aligns with our understanding of NPP in these systems being temperature limited in winter and the idea that warmer winter conditions could help dampen overall declines in NPP through time as California experiences prolonged drought (Griffin and Anchukaitis 2014, Swain et al 2014, Diffenbaugh et al 2015, Dong et al 2019. Similar relationships could be used to explain how increasing early spring precipitation in montane hardwoods and cooler spring temperatures in chamise-redshank chaparral helped buffer NPP declines through time. Not all interactive effects between climate and time on NPP, however, were intuitive. For example, both positive and negative deviations of mid spring temperatures from average were associated with a greater decline in NPP through time in blue oak systems and drier falls appeared to be associated with a smaller decrease in NPP through time in coastal scrub. The range of factors that drive climate and time may interact to alter NPP, including nonlinear responses of NPP to climate and the potential for NPP to be impacted by the synchronization of precipitation and temperature or possibly long term changes in vegetation health or structure (Knapp and Smith 2001, Huxman et al 2004, Bai et al 2008, Chaplin-Kramer and George 2013, Chen et al 2018, Al-Yaari et al 2020. While we are not able to assign a causal mechanism to these complex interactions, a deeper understanding of how interannual climate conditions differ through time and covary with drought and other climate extremes appears to be crucial to understand how California's vegetation may or may not support the State in meeting its climate goals through time. We observed a significant positive asymmetric response of NPP to water availability anomalies in 6 of 7 vegetation types, with a particularly strong positive NPP asymmetry in annual grasslands. Numerous studies investigating the asymmetric response of vegetation productivity to precipitation, particularly in grasslands, have found both positive and negative NPP asymmetries, with positive asymmetries often found in semi-arid ecosystems (Knapp and Smith 2001, Bai et al 2008, Wu et al 2011, Hoover et al 2014, Knapp et al 2017, Wilcox et al 2017, Al-Yaari et al 2020. Our findings suggest an overall positive NPP asymmetry in semi-arid ecosystems in California. Our use of SPEI to quantify anomalously wet and dry years is distinct from previous work which primarily used precipitation. The finding of positive asymmetry with respect to SPEI is important for understanding California's vegetation response to climate change, as water availability impacts on vegetation will be driven simultaneously by more variable precipitation and higher temperatures that increase rates of evapotranspiration (Cayan et al 2012, Diffenbaugh et al 2015, Polade et al 2017, Swain et al 2018. We also highlight vegetation types and regions which may be more resilient to changes in climate than others. Extreme positive response to climate variability in annual grasslands and low but stable response to climate variability in coastal oak woodland and montane hardwood cover types both have conservation value. Similarly, areas on the east side of the coast range (x in figure 6(b)) and along the southern edge of the Central Valley (y in figure 6(b)) show high NPP asymmetry. Collectively, these lands represent potential hotspots for NPP resilience in California which could be enhanced through conservation or restoration.
Our findings on NPP asymmetries also extend previous work by providing fine-scale, spatially explicit estimates of NPP asymmetry across a variety of California's ecosystems and identifying a suite of predictors associated with variation in NPP asymmetry within different vegetation cover types. Five variables consistently appeared as top predictors of NPP asymmetry in each of the seven vegetation types: latitude, longitude, functional group evenness mean and standard deviation, and annual maximum temperature standard deviation. These findings are consistent with a number of previous studies in other systems. For example, precipitation and interactions between precipitation and temperature have been widely noted as primary drivers of NPP asymmetry (Knapp and Smith 2001, Bai et al 2008, Hoover et al 2014, Wilcox et al 2017, Al-Yaari et al 2020, Liu et al 2021. Our predictors of latitude and longitude correspond to broad gradients in precipitation and temperature found across the state with higher latitude systems receiving more precipitation and lower temperatures than lower latitudes, and precipitation generally declining and temperature variation increasing with distance from moderating effects of the Pacific Ocean. Likewise, associations between functional group evenness and NPP asymmetry in our analysis is consistent with a broad body of research showing biodiversity as a stabilizing force on NPP across space and through time (Tilman et al 2001, 2006, Weigelt et al 2008, Wang et al 2019. Lastly, the influence of maximum temperature variability in models of NPP asymmetry across vegetation types may reflect an increasing role of temperature on ecosystem water stress. For example, Griffin and Anchukaitis (2014) showed that record-high temperatures in 2014 exacerbated drought by 36% under low but not unprecedented precipitation conditions, making 2014 the worst drought year in 1200 years. While these findings cannot be used to generate site specific conservation recommendations, these general patterns begin to highlight some of the key factors that may influence the spatial variation in ability of vegetation to store carbon in a future characterized by more extreme dry and wet years.

Conclusion
California's 14 million hectares of grasslands, shrublands and woodlands are fundamental to the State's overall climate strategy. Without the ability of these systems to fix and store carbon each year the state will not likely meet its sequential carbon emissions reduction and carbon neutrality goals. For the first time, we examine the abiotic and biotic drivers associated with spatial and temporal variation in NPP across these major vegetation cover types. Such insight could help guide conservation and restoration decision-making. This understanding also is essential for identifying how changing climate conditions may impact the NPP of these systems. We found that spatial and temporal NPP variation of most vegetation cover types is heavily influenced by precipitation and temperature seasonality. Thus, while mean climate impacts NPP, the seasonality of climate is uniformly important to the ability of these systems to fix carbon and forecasted changes in climate seasonality likely will have significant impacts on the amount of carbon these systems can fix. We also note a clear and largely uniform decrease in NPP through time, associated with the extensive 2012-2016 megadrought that California experienced. Our analysis of asymmetry of NPP response to extreme wet and extreme dry years shows that not all ecosystems within a vegetation type are equally able to respond to climate extremes. In some cases, ecosystems in drier and warmer locations appear to be in a substantial stronger position to maintain NPP during extreme dry years while also being able to rapidly capitalize on extreme wet years by upregulating NPP. Our findings contribute to a growing body of evidence that points to climate seasonality and extremes as a major factor driving ecosystem productivity in a variety of ecosystems across the globe. Our analysis here is not able to establish causal relationships but is able to characterize a range of climate risks and conservation opportunities that exist in California that may either help or hinder the State's progress towards their climate goals. This foundational understanding could help guide future mechanistic work, inform conservation decision making, and identify future threats to the ability of these systems to fix carbon.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: http:// rangeland.ntsg.umt.edu/data/rap/rap-vegetationnpp/v2/ (Robinson et al 2019).