Potential natural vegetation and NPP responses to future climates in the U.S. Great Plains

. Asymmetric climate projections throughout the U.S. Great Plains may intensify the existing latitudinal temperature gradient and magnify the longitudinal precipitation gradient. These potential changes present a unique challenge to understanding the ecological consequences of future climates in the region. Here we investigate how climate change may affect the spatio-temporal patterns of potential natural vegetation types (PVT) and net primary production (NPP) throughout the 21st century with the global dynamic vegetation model MC2. Simulations were driven by projected climate variables from ﬁ ve global climate models under representative concentration pathway (RCP) 8.5. MC2 simulated C3 and C4 grassland, shrubland, forest, and woodland (shrubland + forest) PVTs, and total NPP for each PVT. The largest increases in woodland and grassland NPP occurred in the Northern Plains (17.5% and 4.7%), followed by the Central Plains (10.6% and 0.0%), while NPP in the Southern Plains remained unchanged compared to historic means (1981 – 2010). A shift from grassland to woodland in the Northern and Central Plains further affected regional NPP; regional woodland NPP increased 72% and 26% in the Northern and Central Plains, respectively, while regional grassland NPP decreased 18% and 12%, respectively. The most pronounced shift in PVT was associated with increasing, rather than decreasing, mean annual precipitation in the Northern Plains where grassland contracted in response to westward expansion of woodland. C3 grassland was gradually replaced by C4 grassland in the Northern Plains by 2080, and only a trace remained at centuries end. C3 grassland decreased to a trace amount ca. 2060 in the Central Plains, while C4 grassland increased slightly. The relative stability of PVTs in the Southern Plains suggests that species and functional trait diversity may buffer grassland responses to future climates by providing the capacity for species reordering. The asymmetric response of simulated vegetation and NPP to 21st century climate change suggests that the provision of ecosystem services — beef cattle production, carbon sequestration, and grassland bird habitat — will be modi ﬁ ed in distinct ways along a latitudinal gradient throughout the Great Plains.


INTRODUCTION
The U.S. Great Plains are dominated by a continental climate with pronounced latitudinal temperature and longitudinal precipitation gradients (K€ oppen 1936(K€ oppen , Shafer et al. 2014). These gradients have a prominent effect on ecological patterns and processes, and human activity, including land use and population density (Conant et al. 2018, Seager et al. 2018b). The potential for asymmetric patterns of climate change to modify existing climatic gradients throughout the region presents a unique challenge to understanding the ecological consequences of future climates (Epstein et al. 1998, Bradford et al. 2006.
The latitudinal gradient in drought severity may further increase given projections for increasing mean annual precipitation (MAP) in the Northern Plains, and increasingly variable MAP in the Southern Plains , combined with higher mean annual temperatures (MAT; Vose et al. 2017). MAT is approximately 10°C higher in the Southern than the Northern Plains (PRISM Climate Group), and it is projected to increase similarly throughout the region (5°C in late century). Warming-induced increases in evapotranspiration have been the major contributor to increased drying throughout the region (Kukal andIrmak 2016, Wehner et al. 2017). Increasing climate variability and extreme events are anticipated to accompany the projected changes in MAT and MAP (Conant et al. 2018, Kloesel et al. 2018. Several investigations have suggested an increasing probability of severe drought in the Central and Southern Plains as the century progresses (Christian et al. 2015, Cook et al. 2015, Seager et al. 2018a. Collectively, climate change may have substantial ecological impacts on net primary production (NPP), species range distributions, plant community composition, including woody plant cover, and the provision of diverse ecosystem services (Polley et al. 2013, Conant et al. 2018, Kloesel et al. 2018. Wildfire frequency and extent have increased in the region during the past several decades (Barger et al. 2011, Donovan et al. 2017, Wilcox et al. 2018, and projections for greater interannual precipitation variability in the Southern and Central Plains suggest that this trend is likely to continue (Weatherly andRosenbaum 2017, Stambaugh et al. 2018). The consequences of climate change will both directly and indirectly impact human livelihoods and rural economies by altering the economic viability of dryland cropping and rangeland beef cattle production (Shafer et al. 2014). These climatic effects may further manifest as shifts in land ownership, use, and cover, including the relative proportion of grasslands and croplands, and modify the provision of ecosystem services throughout the Great Plains (Conant et al. 2018, Kloesel et al. 2018, Seager et al. 2018a).
Here we investigate how climate change may modify the spatio-temporal patterns of potential natural vegetation types (PVT) and net primary production with a global dynamic vegetation model, MC2 (Bachelet et al. 2015). MC2 simulations were driven by climate variables from five global climate models under representative concentration pathway (RCP) 8.5. MC2 produced spatio-temporal simulations for C3 and C4 grassland, shrubland, forest, woodland (shrubland + forest) PVTs, and total NPP for each PVT from 2015 to 2099.
This assessment will provide a useful regional reference from which to address the potential ecological consequences of future climate variability and change by addressing the following questions: How will the relative proportions of grassland and woodland PVTs be altered by future climates? How will C3 and C4 grassland PVTs respond to climate change throughout the century? How will changes in temperature and precipitation affect the temporal and spatial trends of NPP at local and regional scales? Finally, we draw inferences for the future provision of important regional ecosystem services-rangeland cattle production, C sequestration, and obligate grassland bird habitat-based on changes in climate variables, PVTs, and associated NPP.

Site description
The U.S. Great Plains encompass approximately 1.3M km 2 in portions of 10 states in the central USA. The study domain for this investigation was bound by the 95th and 105th longitudes west, and the 49th and 30th latitudes north (Fig. 1). The region was divided along the 36th and 42nd latitude into Northern, Central, and Southern Plains.
Native vegetation consists of three major grassland ecoregions that follow a pronounced east-west precipitation gradient: shortgrass prairie in the semiarid west, mixed-grass prairie in the central region, and tallgrass prairie in the mesic eastern region (K€ uchler 1964). However, the vast majority of tallgrass and more mesic mixed-grass prairies have been converted to cropland early last century with the remaining native grassland occurring primarily west of the v www.esajournals.org 100th meridian (Yu and Lu 2018). Grasslands and cropland, in similar proportions, currently account for 89% of the total land cover in the Great Plains (Drummond et al. 2012).
The region's climate is broadly continental, but it ranges from cold semiarid in the northwest to humid subtropical in the southeast (K€ oppen 1936). Annual precipitation peaks bimodally (spring and fall) in the Southern Plains and unimodally (summer) in the Northern Plains. Large temperature fluctuations exist between winter and summer seasons (Shafer et al. 2014). Mean annual temperature (MAT) from 1981 to 2010 in the Northern, Central, and Southern Plains were 6.6°, 11.5°, and 16.9°C; mean annual precipitation (MAP) during the same period were 511, 623, and 683 mm, respectively (not shown; RISM Climate Group).

Climate data
Downscaled climate projections for MAT and MAP were used from five climate models within the Coupled Model Intercomparison Project phase 5 (CMIP5; Taylor et al. 2012). Projections were downscaled by Abatzoglou (2013) to a 4-km grid using the MACAv2-METDATA downscaling algorithm. Five models were selected from 41 available CMIP5 models based on availability of downscaled versions and a literature assessment of the accuracy of retrospective model projections to observed climate averages and extremes. This assessment was based on concordance with the non-downscaled version of each model for winter and summer total precipitation and average temperature, the number of high precipitation days, and the number of hot days for central North America (Sheffield et al. 2013); the 95th percentile of precipitation, growing season length, and daily average minimum and maximum temperatures globally (Sillmann et al. 2013), as well as model responses to changes in ENSO circulation and multi-year trends of temperature and precipitation in the south-central USA (D. Rosendahl, personal communication). Models with a seasonal precipitation bias >25% or seasonal temperature bias >4°C were excluded, and the remaining 7 models were ranked based on their accuracy in the aforementioned criteria.
Five of these climate models were selected (Table 1) based on their respective differences with the historical climate data to produce a broad range of future climate projections (dry/ wet/warm/cool/close to ensemble average) and, subsequently, a diverse range of MC2 simulations. Climate projections were based on representative concentration pathway 8.5 (RCP 8.5), the highest of a series of radiative forcing scenarios (Riahi et al. 2011, van Vuuren et al. 2011, Hayhoe et al. 2017. Simulations for both RCP 4.5 and 8.5 were initially assessed, but the outcomes were similar so we chose to present only RCP 8.5, which appears to most appropriately represent current emission trends (IPCC 2014, Sanford et al. 2014. PRISM climate observations for MAT and MAP from 1981 to 2010 were used as a 30year reference climatology, a widely accepted temporal reference for climate data analysis (WMO 2011). All climate data were acquired as gridded datasets with a 4-km spatial resolution.

MC2 dynamic global vegetation model
The dynamic global vegetation model MC2 (Bachelet et al. 2018) and its predecessor, MC1, (Bachelet et al. 2001, Bachelet 2015 have been used widely to investigate ecosystem responses to climate change, including vegetation dynamics and range shifts (Daly et al. 2000, Neilson et al. 2005, Bachelet et al. 2015, consequences of various fire regimes and atmospheric CO 2 concentrations (Lenihan et al. 2008, Sheehan et al. 2019, forest function (Turner et al. 2015, Kim et al. 2017, and carbon sequestration (Bachelet et al. 2018). The model contains three modules that simulate biogeography, biogeochemistry, and wildfire interactions (see Bachelet et al. 2015 for further detail).
The biogeography module simulates shifts in potential vegetation types based on climate and biomass thresholds (Bachelet et al. 2015 The model simulates dynamics of plant lifeforms (functional groups), rather than individual species. The module uses environmental gradients of minimum monthly temperature and growing season precipitation to simulate the relative dominance of woody lifeforms. The relative dominance of C3 vs. C4 grasses, including forbs, sedges, and other herbaceous vegetation, is simulated by calculating the potential production of pure C3 and pure C4 vegetation types using soil temperature. The model relies on thresholds of carbon pool values to distinguish between forest, shrubland, and grassland PVTs. The MC2 model simulates potential natural vegetation, not actual vegetation cover, which is influenced by socio-political drivers in addition to biophysical variables. MC2 and the CMIP5 climate models used here do not account for feedbacks between land cover and climate as is the case for Earth system models (e.g., Community Earth System Model version 1, CESM1) that possess the capacity to combine climate, land cover, ice, and carbon cycle components (Hurrell et al. 2013).
The biogeochemistry module is a modified version of the CENTURY model (Parton et al. 1994) that simulates carbon and nitrogen cycles, including NPP and ecosystem carbon balance, decomposition, and soil respiration. NPP is determined by temperature, soil water availability, soil nitrogen, and atmospheric CO 2 (Bachelet et al. 2015). Projected increases in annual atmospheric CO 2 concentrations were prescribed by RCP 8.5 (Riahi et al. 2011) for all five climate models. This RCP projects an increase from current concentrations of 415-936 ppm in 2100 (Hayhoe et al. 2017). The CO 2 concentration informed simulation of NPP and vegetation types, which influenced competition for light, water, and nutrients between herbaceous and woody vegetation (Bachelet et al. 2015. The fire module simulates fire occurrence, area burned, and fire impacts including mortality, consumption of aboveground biomass, and nitrogen volatilization. Wildfire occurrence is simulated in two ways: natural fire presence and fire suppression, the latter of which assumes that fires below a certain threshold can be extinguished while fires above the threshold cannot (Sheehan et al. 2019). Fire occurrence is simulated as discrete events in response to calculated ignition probabilities. The module runs on daily time steps by using a randomly distributed set of daily precipitation values derived from monthly precipitation values. Leaf moisture content is calculated as a function of the ratio of available soil Note: Changes in temperature and precipitation represent differences between the projected mid-century mean (2040-2069) and the mean historic reference (1971PRISM Climate Group).
v www.esajournals.org water to potential evapotranspiration, and it is used to determine live fuel moisture content and fire behavior. To reflect a realistic geographic extent of a fire under assumed ignitions, the fire module limits the area burned with an algorithm based on fire return interval and years since last fire. In the fire module, each vegetation type is assigned both a maximum and minimum fire return interval (Bachelet et al. 2015, Sheehan et al. 2019. Model simulations with an activated fire mode are presented here to reflect the recent increase in fire frequency and magnitude throughout the Great Plains (Donovan et al. 2017). The fire presence mode simulated increased grassland area and reduced woodland area in all three regions, compared to the fire suppression mode, as anticipated. These differences were greatest in the Northern Plains and smallest in the Southern Plains (Appendix S1: Fig. S1a-c). Differences in total annual NPP (g C/m 2 ) between the two fire modes were similar among the three regions and less pronounced than differences in potential vegetation type.
The MC2 model runs in three successive phases. In the initialization phase, the biogeography module generates a map of PVTs for the average climate between 1895 and 1924. This map is used by the biogeochemistry module to calculate initial values for carbon and nitrogen pools associated with each vegetation type and their prescribed fire return intervals. The initialization phase ends when the resistant soil carbon pool size changes by less than 1% among successive years. Spinup, the second phase, is run iteratively using detrended historical climate data  to allow for readjustments of vegetation type and carbon pool sizes in response to interannual variability and simulated wildfires. The spinup phase ends when the net biological production (net ecosystem production minus carbon consumed by wild fire) reaches a value near zero. In the third, transient phase, the model is run with time series of historical and future climate variables.
NPP simulations of the MC2 model were compared to Landsat-based NPP data provided by the Rangeland Production Monitoring Service (RPMS) from 1984 to 2010 (Reeves et al. 2020). Simulations of annual aboveground primary production (ANPP; g/m 2 ) for grasslands were derived from total g C/m 2 model output via a 5:1 root-to-shoot ratio (Jackson et al. 1996, Mokany et al. 2006) and a 1:2 carbon-to-biomass ratio (Zhou et al. 2018). The spatial domain of this analysis coincided with the three dominant grassland types identified by RPMS in the Great Plains-northern mixed-grass prairie (NMGP), central and southern mixed-grass prairie (SMGP), and shortgrass steppe (SGS)-that were assessed in this investigation (Reeves et al. 2020; Appendix S2: Fig. S1). MC2 and RPMS data were spatially masked to include only natural grassland and preclude pastures, agriculture, urban, water, barren, and other non-vegetated surfaces.

Data analysis
MC2 was used to develop both historic  and future (2015-2099) simulations of PVTs and NPP. Future and historic MC2 simulations were acquired from Dominique Bachelet (Conservation Biology Institute, www.consbio.org, personal communication) in annual time steps from 2015 to 2099 and 1981 to 2010, respectively, as gridded datasets with a 4km spatial resolution fitted to the climate data. Future MC2 simulations were driven by aforementioned downscaled CMIP5 climate projections, and historic MC2 simulations were driven by PRISM climate observations (PRISM Climate Group).
Twenty individual PVTs were simulated by MC2 throughout the entire region. These PVTs were consolidated into three broad categories that conformed to K€ uchler's (1964) vegetation classification of the Great Plains (Table 2). The grassland PVT was divided into C3 and C4 grasslands; the shrubland PVT consisted of C3 and C4 shrubland (i.e., photosynthetic designation refers to associated grassland species); and forest PVT was comprised of temperate and subtropical woodlands and forests. The analysis focused on spatial and temporal shifts of C3 and C4 grassland, shrubland, and forests throughout the 21st century. Shrubland and forest PVTs were further aggregated into a woodland category for a portion of the analyses to emphasize grasslandwoodland dynamics. Each grid cell in MC2 was assigned a single PVT and corresponding NPP value. The proportional area of an individual PVT was defined as the percentage of all grid cells occupied by the PVT within a region in a v www.esajournals.org given year. In the case of woodland PVT, grid cells occupied by either shrubland or forest PVTs were aggregated. Similarly, proportional total NPP of a PVT was defined as the percentage of total NPP in a region assigned to a specific PVT.
Simulations of the proportional area of PVTs and NPP were calculated annually and compared to projected historical mean values . Projected NPP is expressed as relative change to eliminate the need to partition total NPP (g C/m 2 ) simulated by the MC2 model into above and belowground NPP for each PVT. Historical mean values for PVTs and NPP were simulated by MC2 with observed climate data (PRISM Climate Group). MC2 variables represent unweighted ensemble averages of all five MC2 simulations, each driven by individual climate model projections (Table 2). However, transitions from C3 to C4 grassland PVT (Fig. 5) and NPP and PVT simulations of grassland and woodland (Fig. 6) were presented in annual and decadal time steps as individual MC2 simulations and their unweighted ensemble average. Simulated changes in grassland and woodland NPP per m 2 were combined with projected spatio-temporal changes for both PVTs to develop regional NPP estimates. All simulations of PVTs and NPP were developed under conditions of an activated wildfire module.

Projected temperature and precipitation
Projected MAP for 20-year increments beginning in 2020 were 6.9-12.4% (35-63 mm) higher in the Northern Plains and 1.6-6.8% (10-43 mm) higher in the Central Plains, compared to their respective historic means (Fig. 2a, b). Projected MAP in the Southern Plains remained 1.6% (11 mm) above historic means until 2039, but then decreased 0.5-4.8% (3-33 mm) below historic values after 2039 (Fig. 2c). MAT in each of the three regions are projected to rise between 4.9°and 5.5°C by the end of the century (Fig. 2d-f).

Projected spatio-temporal shifts in PVTs
Grassland was dominant west of the 100th meridian throughout the Great Plains until approximately 2050 (Fig. 3), while shrubland and forest dominated the eastern half of the region. However, in the Northern Plains grassland grid cells rarely occupied over 60% of a given longitudinal transect (0.5 degrees width) with the remaining grid cell representing woodland throughout the century (Fig. 3a). In contrast, in the Central and Southern Plains grassland grid cells occupied up to 90% of a given longitudinal transect, primarily in the western half of both regions (Fig. 3b, c). Grassland shifted westward and decreased in area in the Northern Plains, and to a lesser extent in the Central Plains, c. 2050 as shrubland and forest PVTs increased in the eastern portions of both regions. In the Southern Plains, the proportion of PVTs remained constant throughout the century, with the minor exception of eastward grassland expansion near the end of the century.
MC2 simulations of a grassland shift toward the northwest in the Northern Plains and west in the Central Plains began to diverge among the five climate model projections (wet/dry/warm/ cool/average) as the century progressed (Fig. 4). Agreement among MC2 simulations was greatest early in the 21st century where 94% of total grassland area was projected by all five MC2 simulations. Model agreement decreased to 78% by mid-century and further decreased to 70% by the end of the century. Although each climate projection and associated MC2 simulation are assumed to have a similar probability of occurrence, greater agreement among simulations is suggestive of a more probable outcome.
Increasing divergence among MC2 simulations throughout the century, particularly in the Northern Plains (Fig. 4a, c), may have resulted from disproportionate changes in MAT and MAP among individual climate projections. The difference in MAT between the warm and cool climate models increased from 1.9°C between 2015 and 2049 to 4.3°C between 2050 and 2099 in the Northern Plains, much more than during the same periods in the Southern Plains (1.6°-2.2°C, not shown). At the same time, differences in MAP projections of the wet and dry climate models increased from 108 to 127 mm during these respective periods in the Northern Plains, but decreased from 215 to 174 mm in the Southern Plains (not shown).
C4 grassland gradually replaced C3 grassland as the century progressed (Fig. 5a-c). In the Northern Plains, C4 grassland exceeded C3 grassland c. 2030. C4 grassland occupied an average of 31% of the Northern Plains between 2027 and 2099, while C3 grassland gradually decreased to a trace amount by the end of the century. In the Central Plains, C4 grassland occupied between v www.esajournals.org 7 October 2020 v Volume 11(10) v Article e03264 41% and 47% of the region, while the proportion of C3 grassland gradually decreased from 10% in 2015 to a trace in 2099. In the Southern Plains, C4 grassland occupied 45% of the region between 2015 and 2085 and increased to just over 50% in 2099, while C3 grassland existed in only trace amounts throughout the century. Individual simulations for C4 grassland area showed greater variation in the Northern than Southern Plains throughout the century.

Simulated NPP per m 2 and per region
Simulated NPP per m 2 increased from the Southern to the Northern Plains for both grasslands and woodlands in all three regions compared to historic means (Fig. 6a-c). In the Northern Plains, mean projected grassland and woodland NPP per m 2 , from 2015 to 2099, were 4.7% and 17.5% higher, respectively, than historic means (Fig. 6a). In the Central Plains, grassland NPP per m 2 remained similar to the historic mean while woodland NPP per m 2 was 10.6% greater (Fig. 6b). In the Southern Plains, projected NPP per m 2 remained similar to historic means for both grassland and woodland PVT (Fig. 6c).
MC2 simulations underestimated ANPP production for NMGP and SMGP by 23.2% and 17.9%, respectively, and overestimated production for SGS by 18.5% relative to Landsat-based RPMS NPP data (Appendix S3: Fig. S1). Standard deviations for MC2 simulations were 19.8% and 4.3% smaller for NMGP and SMGP, respectively, and 35.1% greater for SGS, compared to RPMS. RPMS and MC2 were strongly correlated in all regions (0.59-0.70), and 49.5%, 35.2%, and 34.8% of the variance (R 2 ) between MC2 and RPMS for NMGP, SMGP, and SGS, respectively, were explained (Appendix S3: Fig. S1). These correspondent values indicate that the ANPP simulations of the MC2 model, which are based on potential vegetation types, are relatively robust compared to those derived from Landsat data, especially considering uncertainty in root-to-shoot ratios.

Proportion of grassland and woodland PVTs
Spatio-temporal changes in the proportions of grassland and woodland are expressed as decadal averages of future simulations for each of the five climate projections compared to their respective historic values (Fig. 6d-f). Differences between the proportions of grassland and woodland were greatest in the Northern Plains where woodland increased 45% and grassland decreased 35% (Fig. 6d). However, differences among individual MC2 simulations were also greatest in the Northern Plains which resulted in a more diffuse boundary between grassland and woodland. In the Central Plains, the grassland PVT decreased 13%, with a similar increase in woodland (Fig. 6e). Differences among individual MC2 simulations were smaller than in the Northern Plains.
In the Southern Plains, simulated changes for both grassland and woodland area were minimal; but, unlike the Northern and Central Plains, grassland area increased slightly near the end of the century while woodland decreased (Fig. 6f).

Total NPP per region
The product of NPP per m 2 and the proportion of PVT per region provides an estimate of regional NPP for each PVT throughout the century (Fig. 6g-i). In the Northern Plains, grassland NPP decreased 18%, while woodland NPP increased 72% above the historic mean (Fig. 6g). In the Central Plains, changes in NPP between grassland and woodland PVTs were less pronounced, and individual model projections deviated less than that in the Northern Plains (Fig. 6h). Regional grassland NPP decreased 12%, while total woodland NPP increased about 26% compared to historic means. Regional NPP exhibited the least change in the Southern Plains, grassland NPP remained unchanged and woodland NPP increased only 3% compared to historic values (Fig. 6i).

DISCUSSION
Inferences were drawn for each of the initial research questions based on the MC2 model simulations. First, NPP per m 2 progressively v www.esajournals.org increased from the Southern to the Northern Plains and mean increases were consistently greater for woodland than for grassland. Second, major shifts in PVTs were confined to the Northern Plains, where grassland shifted westward in conjunction with woodland expansion. Third, the proportional area of C4 grassland approached that of the C3 grassland in the 2020s in the Northern Plains, and C3 grassland gradually decreased to only a trace by the end of the century. Collectively, these results indicate that the provision of ecosystem services may be modified in distinct ways along a latitudinal gradient throughout the Great Plains.

Climate projections
MAT and MAP projections from the five selected GCMs are consistent with previous climate projections for the Great Plains, including the National Climate Assessment (NCA; Shafer et al. 2014) and the 6th Assessment Report of the Intergovernmental Panel on Climate Change (IPCC 2014). MAT increased incrementally by 2.7°and 5.2°C in mid-century and late century, respectively, in our projections (Fig. 2), compared to increases of 3.2°-6.6°C late in the century as reported by the NCA . Precipitation in the Northern Plains is projected to increase 10-30% during the growing season while winter precipitation is projected to decrease 10-20% ). The MC2 ensemble average projected a mean increase in MAP of 8.9% from 2015 to 2099 compared to the historic mean for the Northern Plains. MAP in the Southern Plains is projected to fluctuate between À10 and +10% in any season ). Our climate projections indicated that MAP in the Southern Plains will increase 1.4% between 2015 and 2049, but then decrease 2.7% after 2050, in addition to increased interannual variability in the second half of the century (Fig. 2).
Climate projections often underestimate extreme climate events (Huang and Gao 2017) which may have greater ecological effects than those driven by more incremental changes in MAT and MAP. Large and frequent variation in projected MAP below those of mean historic values in the Central and Southern Plains corroborate previous projections of future severe droughts (Cook et al. 2015, Seager et al. 2018a. Consequently, eastward shifts in the semiarid zone (100th meridian) throughout the Great Plains as suggested by Seager et al. (2018a), and a greater probability of severe drought after midcentury in the Central and Southern Plains (Cook et al. 2015) may have greater impacts on the spatial distribution of PVTs.

NPP simulations
The MC2 ensemble average indicated that grassland NPP per m 2 was 4.7% and 0.6% greater in the Northern and Central Plains, respectively, throughout the century than the historic reference. These values are considerably lower than those recorded at the Prairie Heating and Carbon Dioxide Enrichment (PHACE) experiment in southeastern Wyoming. Combined warming and CO 2 enrichment as projected by RCP 8.5 increased ANPP of mixed-grass prairie by 25% (Mueller et al. 2016), 33% (Morgan et al. 2011), and 38% (Augustine et al. 2018). Aboveground NPP was observed to increase 41% in response to CO 2 enrichment (700 ppm) in open-top chambers experiments in the shortgrass steppe of eastern Colorado (Morgan et al. 2004). Free-air CO 2 enrichment at the BioCON experiment in Minnesota increased ANPP of grassland species 33%, which is similar to the Wyoming PHACE experiment, but the authors cautioned that soil water and nutrient limitations may prevent this optimal value from consistently being realized (Reich et al. 2014). However, MC2 simulations run with individual climate models projected NPP increases as high as 31% (Fig. 5a-c), approaching those observed in experimental investigations late in the century (Morgan et al. 2011, Reich et al. 2014, Mueller et al. 2016, even though the MC2 ensemble average only projected increases of 3.8% and 5.4% for 2015-2049 and 2050-2099, respectively, for the Northern Plains. Variation among individual MC2 simulations was likely driven by differences in MAT and MAP projections, but other variables may have been involved. Cumulative April-July precipitation and actual evapotranspiration (AET), as well as the ratio of cumulative April-July AET to potential evapotranspiration rates, are more strongly correlated to interannual changes in v www.esajournals.org ANPP in the Great Plains than is MAP, which was used in these model simulations (Chen et al. 2019). Further, CO 2 fertilization has been observed to have the greatest impact on grassland ANPP when it coincides with spring precipitation, which is largely the case throughout the Great Plains (Hovenden et al. 2019).

PVT simulations
The most pronounced shifts in PVTs occurred in the Northern Plains where proportional grassland area decreased in response to the westward expansion of woodland. This shift in PVT was associated with increasing, rather than decreasing, MAP (Barger et al. 2011, Scholtz et al. 2017. Increasing MAP and atmospheric CO 2 concentrations may have sufficiently offset increasing evapotranspiration associated with higher MAT to allow greater percolation of deep soil water. Woody plants may have gained a competitive advantage over grasses based on their ability to access this deep soil water, that is, dual layer soil hypothesis (Noy-Meir 1973, Tietjen et al. 2017. The replacement of C3 grassland by C4 grassland in the Northern and Central Plains prior to mid-century is widely anticipated, and it is primarily assumed to result from the competitive advantage that increasing MAT conveys to the C4 photosynthetic pathway (Epstein et al. 1997, 1998, von Fischer et al. 2008). However, this interpretation merits caution, because well-established short-term responses of C3 and C4 species to elevated atmospheric CO 2 may not necessarily predict long-term outcomes (Mueller et al. 2016, Reich et al. 2018. Future climate had a much smaller influence on the spatio-temporal variability of PVTs in the Central and Southern Plains even though higher MAT and more variable MAP were projected in these regions. Although not simulated by MC2, we speculate that grassland stability in the Southern and Central Plains may be a consequence of the large number of species that are uniquely adapted to and distributed along a longitudinal precipitation gradient, that is, short, mid-species, and tallgrass species from west to east (K€ uchler 1964(K€ uchler , Epstein et al. 1996(K€ uchler , 1998. This diversity of species and functional traits may buffer grassland responses to future climates by providing the capacity for rapid composition change in response to dynamic climatic conditions (Dı́az andCabido 2001, Knapp et al. 2015). For example, more drought-adapted midgrasses shifted eastward into mesic tallgrass prairie during the 1930's drought (Weaver 1943). Further, interannual fluctuations in grassland NPP may be stabilized by increased species evenness, as dominant species decrease and subordinate species increase, in response to warming and elevated atmospheric CO 2 (Zelikova et al. 2014).
The absence of shrubland expansion in the Southern Plains was unexpected given the extent of woodland encroachment that had occurred throughout the region in the previous century (Barger et al. 2011). Three possible explanations, either independently or in combination, may account for this outcome. First, shrubland may have approached a carrying capacity, established by biotic or abiotic processes, that has limited further encroachment in the region (Huang et al. 2018). Second, warmer, drier future climates, and especially extreme events, may contribute to shrub mortality (Twidwell et al. 2014). Third, it is possible that MC2 simulations may have underestimated the area of this PVT based on the climate variables associated with future climates.

Implications for ecosystem services
Climate-induced changes in PVT and NPP may both directly and indirectly modify the provision of ecosystem services. The severe drought of 2011-2014 throughout the Great Plains contributed to adverse agricultural and ecological outcomes (Rippey 2015, Moore et al. 2016). Here we develop implications for three major ecosystems services within the region-rangeland beef cattle production, C sequestration, and obligate grassland bird habitat-based on the projected climate variables and simulated PVTs and NPP derived in this investigation.
Carbon cycling and sequestration are strongly influenced by NPP and the composition of plant functional groups (Sleeter et al. 2013, Ahlstr€ om et al. 2015. Increasing woodland cover is often associated with greater carbon sequestration, but various disturbances may affect the magnitude of this response (Barger et al. 2011, Eldridge et al. 2011. Previous MC2 simulations have indicated that carbon sequestration may increase in the central and western USA throughout the 21st century, primarily as a consequence of an v www.esajournals.org increase in woodland cover (Bachelet et al. 2015).
However, carbon sequestration is strongly mediated by MAP, in addition to the composition of plant functional groups. For example, the Great Plains was a net carbon sink throughout the period 2000-2008 with an annual C sequestration rate of 0.3-47.7 g CÁm À2 Áyr À1 (Zhang et al. 2011). However, 4 yr of successive drought substantially decreased the rate of C sequestration, especially in the semiarid western region. Similar drought-induced decreases in C sequestration have been observed in the shortgrass steppe of eastern Colorado . Therefore, increasing future woodland cover and NPP in response to increasing MAT and MAP in the Northern Plains are likely to increase carbon sequestration in that region. In contrast, C sequestration in the Southern Plains may decrease in response to increasing MAT and surface evapotranspiration, and decreasing MAP, which will collectively decrease soil water availability (Zhang et al. 2011. Increases in both MAT and MAP, with minimal shifts in PVTs, may counteract one another to maintain a stable carbon sink in the Central Plains.
Previous assessments have concluded that beef cattle production will increase in the Northern Plains in future climates, but likely decrease in the Southern Plains, based on the respective changes in climatic variables supporting grassland NPP (Polley et al. 2013, Briske et al. 2015). However, MC2 simulations of woodland encroachment into grasslands and greater increases in woodland than grassland NPP seriously challenge this interpretation. Increasing woodland area and NPP is anticipated to decrease beef cattle production by reducing the amount and availability of forage to support cattle. It has been estimated that a 1% increase in woodland area may reduce cattle production by 2.5% on productive sites in the USA (Anad on et al. 2014). An increase in woodland expansion in the Southern and Central Plains has contributed to the emergence of numerous prescribed burn cooperatives that promote the use of fire to manage woody plants and maintain grasslands (Twidwell et al. 2014). Projected decreases in C3 grassland and NPP in the Northern Plains may further reduce the availability of high-quality forage to support cattle production early in the growing season. However, an extended growing season for C4 grassland may partially compensate for the decrease in C3 grassland production (Reyes-Fox et al. 2014, Hufkens et al. 2016. Woodland expansion will further reduce habitat for obligate grassland bird species in the Northern Plains and to a lesser extent in the Central Plains (Thompson et al. 2014). Woodland encroachment is known to decrease both the richness and abundance of grassland bird species (Scholtz et al. 2017). This will exacerbate recent habitat losses resulting from grassland conversion to croplands (Lark et al. 2015), including a 25% decrease in Conservation Reserve Program (CRP) acreage since 2007 (Morefield et al. 2016).
The most pronounced effects of climate change in the Great Plains may be an intensification of the existing latitudinal temperature gradient and a magnification of the longitudinal precipitation gradient. These asymmetric projections of climate change produced asymmetric responses of simulated PVTs and NPP. NPP increased from the Southern to the Northern Plains throughout the century as expected, but woodland NPP exceed grassland NPP in all three regions. Unexpectedly, the greatest shift in vegetation type in future climates was associated with increasing, rather than decreasing MAP in the Northern Plains, where woodland replaced grassland in the central portion of the region. Simulations of vegetation type and NPP to future climates suggest that the provision of ecosystem services will be modified in distinct ways along a latitudinal gradient. These simulations may provide a reference from which to develop and test hypotheses to further investigate the asymmetric consequences of climate change throughout the Great Plains.

ACKNOWLEDGMENTS
This work was supported by USDA-NIFA grant 2016-67003-24970 to DDB. The insightful comments of two anonymous reviewers greatly improved the manuscript. Linda Joyce constructively commented on an earlier version of this manuscript. IT support for data storage and processing was generously provided by the South Central Climate Adaptation Science Center. Dominique Bachelet and Ken Ferschweiler graciously provided insights into several aspects of MC2 model function.