Forecasting climate change impacts on plant populations over large spatial extents

Plant population models are powerful tools for predicting climate change impacts in one location, but are difficult to apply at landscape scales. We overcome this limitation by taking advantage of two recent advances: remotely sensed, speciesspecific estimates of plant cover and statistical models developed for spatiotemporal dynamics of animal populations. Using computationally efficient model reparameterizations, we fit a spatiotemporal population model to a 28year time series of sagebrush (Arte­ misia spp.) percent cover over a 2.5 × 5 km landscape in southwestern Wyoming while formally accounting for spatial autocorrelation. We include interannual variation in precipitation and temperature as covariates in the model to investigate how climate affects the cover of sagebrush. We then use the model to forecast the future abundance of sagebrush at the landscape scale under projected climate change, generating spatially explicit estimates of sagebrush population trajectories that have, until now, been impossible to produce at this scale. Our broadscale and longterm predictions are rooted in smallscale and shortterm population dynamics and provide an alternative to predictions offered by species distribution models that do not include population dynamics. Our approach, which combines several existing techniques in a novel way, demonstrates the use of remote sensing data to model population responses to environmental change that play out at spatial scales far greater than the traditional field study plot.


IntroductIon
Forecasting the impacts of climate change on plant populations and communities is a central challenge for ecology (Clark et al. 2001, Petchey et al. 2015. Population models are ideally suited for meeting such a challenge because they provide a way to link climate drivers directly to population dynamics (Hare et al. 2010, Adler et al. 2012, Ross et al. 2015, Shriver 2016. However, inference from population models is typically limited to small spatial extents because the data required are difficult to collect across broad species ranges. Almost every study of plant population v www.esajournals.org TREDEnniCK ET AL. dynamics relies on demographic observations recorded at the meter-to-submeter scale (see, e.g., Salguero-Gómez et al. 2015). Local-scale demographic data make building population projection models an easy task (Ellner and Rees 2006, Rees and Ellner 2009, Adler et al. 2012), but it is very difficult to extrapolate small-scale studies to large spatial extents with any certainty because the data likely only represent a small subset of parameter space and environmental conditions . The real challenge is not to simply make population forecasts, but to do so at spatial scales relevant to policy and management decisions .
The ideal tool would be a broadscale, dynamic population model (Schurr et al. 2012, Merow et al. 2014), but developing useful models at this scale has been limited by the availability of time series data at large spatial extents and statistical methods for fitting high-dimensional spatial models. Fortunately, new advances in remote sensing and statistics now allow us to overcome both of these limitations. First, new remote sensing (RS) methods are now producing accurate time series of species-specific plant cover at landscape scales. These data can be fit with dynamic population models which include yearly fluctuations in climate as covariates. Such RS time series have revolutionized models of how climate affects ecosystem-level processes (e.g., Running et al. 2004) and have been used to detect longterm trends in plant population abundance (e.g., Homer et al. 2015), but they have yet to be used to drive a dynamic population model. Second, animal population modelers have developed dimension reduction and reparameterization techniques to efficiently fit high-dimensional spatiotemporal models (see Conn et al. 2015 for a review). These new statistical methods have yet to be applied to RS-derived plant population data at broad scales.
Large-scale, spatially explicit population models based on RS data could offer a valuable new way to investigate the effects of large-scale environmental changes playing out at landscape and regional scales. Most current assessments of how plant and animal populations will respond to climate change rely on species distribution models (SDMs). SDMs rely on static associations between contemporary climate and a species' distribution or, more rarely, abundance to project future distribution or abundance (Elith and Leathwick 2009) and they are easily applied at landscape to continental scales (e.g., Maiorano et al. 2013, Clark et al. 2014. However, the shortterm and small-scale population dynamics that actually drive the large-scale distributions of species are not represented in most SDMs. Because SDMs typically rely on occurrence data, their projections of habitat suitability or probability of occurrence provide little information on the future states of populations in the core of their range-areas where a species exists now and is expected to persist in the future (Ehrlén and Morris 2015). Furthermore, because they lack short-term dynamics, SDMs usually cannot produce any estimate of the rate at which local populations will increase or decrease in the near term and instead project a future equilibrium species distribution that may or may not ever be reached. Direct validation of such predictions is extremely rare (Roberts and Hamann 2012). Large-scale dynamic population models could overcome these limitations. They would produce spatially explicit estimates of species abundance within the species range (Ehrlén and Morris 2015), have the potential to model expansion in abundance outside the range when coupled with dynamic models of dispersal, and would provide testable predictions of how populations should respond to short-term climate perturbations. These shortterm predictions also would give modelers the opportunity to repeatedly validate and refine their models (Luo et al. 2011).
Sagebrush (Artemisia spp.) ecosystems offer an ideal testing ground for new spatially explicit population models derived from RS data. Sagebrush species are widely distributed (Kuchler 1964), they are sensitive to climate (Perfors et al. 2003, Miglia et al. 2005, Poore et al. 2009, Dalgleish et al. 2011, Xian et al. 2012, Apodaca 2013, Schlaepfer et al. 2014a, b, Harte et al. 2015, Homer et al. 2015, new landscape-and regionalscale time series of sagebrush cover are now being produced from aerial imagery , and forecasts of future sagebrush ecosystems are in high demand due to the precarious conservation status of the greater sagegrouse (Centrocercus urophasianus) (Arnett and Riley 2015). SDMs typically predict that much of the area occupied by sagebrush ecosystems v www.esajournals.org TREDEnniCK ET AL. today will become unsuitable for sagebrush due to climate change, resulting in a dramatic loss in the extent of sagebrush habitat by the end of this century (Shafer et al. 2001, neilson et al. 2005, Bradley 2010, Schlaepfer et al. 2012, Still and Richardson 2015. Ecohydrology models supply a possible mechanism for sagebrush losses predicted by SDMs: Climate warming could lead to earlier snowmelt, increased evaporation, and ultimately less recharge of deeper soil layers in the spring (Schlaepfer et al. 2012(Schlaepfer et al. , 2014a. in warmer parts of its range, increased temperature could be especially detrimental to sagebrush as it depends on water from deeper soil to survive and grow in this arid region (Pechanec et al. 1937, Schlaepfer et al. 2011, Germino and Reinhardt 2014. in contrast, at higher elevations and in colder regions, warming and earlier snowmelt could lengthen the growing season and increase sagebrush occurrence (Schlaepfer et al. 2012(Schlaepfer et al. , 2014a. Direct observations of individual plants and experimental plots tend to agree with these models: Growth tends to respond negatively to spring and summer temperatures (Miglia et al. 2005, Poore et al. 2009, Apodaca 2013) except at higher elevations where earlier snowmelt may allow for a longer growing season (Perfors et al. 2003, Harte et al. 2015. A large-scale, spatially explicit population model for sagebrush driven by interannual climate variability would provide a valuable new tool for assessing how sagebrush could respond to climate change in the future. Building on recent technological advances in spatial statistics (Latimer et al. 2009, Conn et al. 2015 and anticipating ever-increasing availability of RS data (He et al. 2015), we demonstrate how large-scale plant population models could be used to predict population impacts of climate change. As a proof of concept, we use a process model motivated by Gompertz densitydependent population growth and a remotely sensed time series of sagebrush cover from Wyoming (Homer et al. , 2015. We account for spatial autocorrelation with dimension reduction techniques (Latimer et al. 2009, Conn et al. 2015 and produce spatially explicit estimates of sagebrush percent cover. Unlike most SDMs, our approach models the dynamics of plant abundance through time, and thus is a population model, in the same spirit that models of animal counts through time are population models. The modeling framework we propose can be applied to any spatially explicit time series of plant cover or density, but its application to remotely sensed data products offers the greatest potential to combine the information of population models (e.g., population status and temporal dynamics) and the spatial extent of SDMs.

Data
Remotely sensed time series.-To demonstrate our modeling approach, we use a subset of a remotely sensed time series of sagebrush (Arte misia spp.) canopy cover in Wyoming . As part of a separate study, Homer et al. (2012) estimated sagebrush percent cover using a regression tree to relate ground reflectances retrieved by three sources of optical imagery (QuickBird, Landsat, and AWiFS) to 1780 field observations of sagebrush cover distributed across Wyoming. The regression tree model was further validated using another 297 field observations. For Wyoming sagebrush, the model achieved an R 2 = 0.65 and an out-of-sample RMSE of 5.46% . To hind-cast sagebrush cover, the regression tree model was applied to historical remote sensing images to generate yearly predictions of sagebrush cover for all of Wyoming for the years 1984-2011. This resulted in an annual time series of sagebrush cover at 30 m resolution from 1984 to 2011 (Appendix S2: Fig. S1). in this remote sensing product, values represent the percentage of a 30 × 30 m pixel covered by sagebrush. in our study, we focused on a 5070 × 2430 m subset totaling 13,689 30 × 30 m pixels each year (Fig. 1). Thus, the subset of the remote sensing product we use contains 369,603 observations spanning 27 year-to-year transitions (27 yr × 13,689 pixels).
Climate covariates.-Our approach models interannual changes in plant cover as a function of seasonal climate variables. We used daily historic weather data for the center of our study site from the nASA Daymet data set (available online: http://daymet.ornl.gov/). The Daymet weather data are interpolated between coarse observation units and capture some spatial variation. We relied on weather data for the centroid of our study area. We calculated five climate variables from the Daymet data for the v www.esajournals.org TREDEnniCK ET AL. time period coinciding with our remotely sensed data .
We narrowed our focus to climate covariates that we know are important for sagebrush and that could be calculated from general circulation model projections. The five climate variables in our population model are as follows: (1) cumulative, "water year" precipitation for year t − 2 (lagPpt), (2) year t − 1 fall-through-summer precipitation (ppt1), (3) year t fall-through-summer precipitation (ppt2), (4) year t − 1 average spring temperature (TmeanSpr1), and (5) year t average spring temperature (TmeanSpr2), where t − 1 to t is the transition of interest. We selected these variables a priori based on previous studies (see Introduction), although not all emerge as important predictors in our model.

Additive spatiotemporal model for sagebrush cover
We use a descriptive model for sagebrush cover that includes additive spatial and temporal effects similar to that described by Conn et al. (2015). interannual change in percent cover represents the integrated outcome of recruitment, survival, growth, and retrogression (shrinkage) of individual plants from year to year. We model observed integer percent cover (y) in cell i at time t as conditionally Poisson: where μ i,t is the expected percent cover of pixel i in year t: Our model of percent cover change includes a density-dependent effect of log-transformed cover in the previous year (y i,t − 1 ), climate effects (x t ), and a spatial random effect (η) for each pixel i. Climate effects were standardized [(x i −x)∕σ(x)] to improve convergence during the model fitting stage and to allow for easier prior specification. The intercept, β 0,t , was allowed to vary through time; these random year effects recognize that all observations from a particular year share the same climate covariates and thus are not independent. We used a Poisson likelihood because integer percent cover values in the sagebrush data product can be considered a form of count data. We also evaluated a negative binomial model, but found little evidence for overdispersion beyond what our model was already accommodating via the spatial random effects (η). There was no evidence of zero inflation in our data, but see below (Accommodating zeros) for how we handled the small number of zero percent cover observations. We assume that the remotely sensed estimates of percent cover are "true" and free of error. This need not be the case, and if measurement or sampling error is known, then it could be included in our Bayesian model as a "sampling model" .
The spatial random effect (η) accounts for spatial autocorrelation among pixels that occur near each other in space. Thus, η acts as an offset on the intercept (β 0,t ), creating a spatial field that defines how pixels differ from the mean, on average, in space (e.g., areas of perennially low or high cover, relative to average cover). Fitting the model with a spatial random effect (η) is computationally demanding for large data sets such as ours. The computational demand is due to the required calculations of the spatial covariance matrices, which increase as a cubic function of the number of locations (Wikle 2010). Key to our approach is a dimension reduction strategy that greatly reduces the number of parameters needed to be estimated to account for spatial variation by reducing the size of the spatial covariance matrices that need to be inverted at each Markov chain Monte Carlo (MCMC) iteration. Fitting models that appropriately account for spatial autocorrelation over large spatial extents would not be feasible without these modern techniques. Our dimension reduction strategy expresses the high-dimensional spatial random effect, η, as the product of an expansion matrix, K, and a smaller parameter vector, α (e.g., Hooten et al. 2003, Hooten and Wikle 2007, Conn et al. 2015. We can then approximate the spatial effect as: in this case, α is a m × 1 vector of reduced spatial random effects, and K is a S × m matrix that maps the reduced effects to the full S-dimensional space, where S is the total number of observed locations. Thus, we are able to reduce the effective number of parameters from S to m. The last remaining obstacle is to parameterize the matrix of basis functions, K. We use kernel convolution (Barry et al. 1996, Higdon 1998) to interpolate the spatial random effect between m "knots" that are nonrandomly distributed across the space of our study area. This means we are modeling spatial random effects at the knot level, and we use K to interpolate those effects between knots. We use an exponential kernel density to define the distance-decay function around the knots (w), such that the entries of K are: where and d s,m is the Euclidean distance between the centroid of sample cell s and the location of knot m, and σ is the kernel bandwidth. it is possible, through exhaustive model selection and fitting, to determine the optimal form of the kernel and to estimate optimal values for σ (Higdon 2002, Hooten and. However, given the relative size of our data set and computational limitations, we defined kernels around 231 knots (Appendix S3: Fig. S2) whose nearest-neighbor distances are approximately equal to the range of spatial dependence in residuals from a simple generalized linear model (GLM) fit without climate covariates and the spatial random effect (~500 m; Appendix S3). An infinite number of knots would result in an exact representation of the spatial process and covariance model. Computationally, using an infinite number of knots is not possible; thus, the use of dimension reduction techniques serves as an approximation, where the accuracy increases with the number of knots. Given the trade-off between knot number and computation time, we chose to base our knot number on the spatial dependence as described above.
The Bayesian posterior distribution of our spatiotemporal model can be expressed as:

Accommodating zeros
Our process model (in Eq. 2) includes a log transformation of the observations (log(y t − 1 )). Thus, our model does not accommodate zeros. Fortunately, we had very few instances where pixels had 0% cover at time t − 1 (n = 47, which is 0.01% of the data set). Thus, we excluded those pixels from the model fitting process. However, when simulating the process, we needed to include possible transitions from zero to nonzero percent cover. We fit an intercept-only logistic model to estimate the probability of a pixel going from zero to nonzero cover: where y is a vector of 0s and 1s corresponding to whether a pixel was colonized (>0% cover) or not (remains at 0% cover) and μ i is the expected probability of colonization as a function of the mean probability of colonization (b 0 ). We fit this simple model using the "glm" command in R (R Core Team 2013). For data sets in which zeros are more common and the colonization process more important, the same spatial statistical approach we used for our cover change model could be applied and covariates such as cover of neighboring cells could be included.

Fitting the model
We fit the spatiotemporal model in R (R Core Team 2013) using the "no-U-Turn" Hamiltonian Monte Carlo sampler in Stan (Stan Development Team 2014a) and the RStan package (Stan Development Team 2014b). We obtained posterior distributions of all model parameters from three MCMC chains comprised of 1000 iterations each, after discarding an initial 1000 iterations as burn-in. Short chains of samples are a hallmark of the Stan algorithm, which is extremely efficient. Compared to other samplers, fewer iterations are required to achieve convergence. Each chain was initialized with unique parameter values, and the model was fit in parallel using the Utah State University High-Performance Computing facility. Model fitting required 5 d on a four-node central processing unit with 2× AMD Opteron Processor 4386 @3.10 Ghz, 64 GB of RAM per node, 16 cores per node, and each chain launched in parallel on separate cores. We assessed convergence visually and calculated scale-reduction factors (Appendix S4; R < 1.1 for all parameters) Rubin 1992, Gelman andHill 2009).

Simulating the process
We performed four sets of simulations to (1) compare observed and simulated equilibrium cover, (2) compare observed and simulated yearand location-specific cover, (3) forecast future equilibrium population states under projected climate change, and (4) make temporally explicit forecasts of sagebrush cover starting the final year of our observations and ending in year 2098. Using the posterior distribution of model parameters, we simulated a matrix of pixels equal to the size of the study area (13,689 pixels or matrix elements). For simulations (1) and (3), we initialized all pixels with arbitrarily low cover (1%) and then projected the model forward by randomly drawing climate covariates from the observed climate time series (for 1) or a perturbed climate time series (for 3). We ran equilibrium simulations (1 and 3) for 2000 time steps and then compared the output across simulations, after discarding an initial 100 time steps. To calculate average future equilibrium sagebrush cover, we ran simulation (3) for each global circulation model (GCM) and representative concentration pathways (RCP) scenario TREDEnniCK ET AL. separately and then averaged the results over GCMs. For simulation (2), we initialized each pixel with its actual percent cover value for time t and cell s and projected the model forward one time step and compared the one-step-ahead forecast with the observed value. For simulation (4), we initialized each pixel with the final observed value in 2011 and then projected the model forward based on GCM yearly weather projections. We ran these simulations for each GCM and RCP scenario combination separately and then aggregated the results over the GCMs by calculating the mean and the 90th percentiles for each RCP scenario.
We used the posterior mean of each parameter for all simulations except for (4) where we ran 50 simulations with unique sets of parameters from the chains. Random year effects were included in simulations by randomly drawing a posterior mean year effect (β 0,t ) for each iteration (simulations 1 and 3), using the posterior mean year effect for a specific year (simulation 2), or by a drawing a future-year random effect from the posterior mean and standard deviation of the mean intercept (simulation 4, e.g., β 0,T ∼ normal(β 0 ,σ 2 β 0 ) for some future year T). Our simulation approach provides a reasonable and computationally efficient approximation to the true posterior predictive mean when used in these scenarios with our data.
We required future projections of climate for our study area to conduct the equilibrium and temporally explicit forecasts described above. Thus, we used the most recent climate projections from the intergovernmental Panel on Climate Change (iPCC), the Coupled Model intercomparison Project 5 (CMiP5; available online: http://cmip-pcmdi. llnl.gov/cmip5/). The CMiP5 provides projections from a suite of GCMs; we used projections from 18 GCMs (Appendix S1: Table S1) that produced weather projections for three RCPs: RCP 4.5, RCP 6.0, and RCP 8.5 (described online: http://tntcat. iiasa.ac.at/RcpDb/). The three RCPs correspond to stabilization of radiative forcing before 2100, after 2100, and ongoing increase in greenhouse gas emissions, respectively.
To simulate equilibrium sagebrush cover under projected future climate, we applied average projected changes in precipitation and temperature to the observed climate time series.
For each GCM and RCP scenario combination, we calculated average precipitation and temperature over the 1950-2000 time period and the 2050-2098 time period. We then calculated the absolute change in temperature between the two time periods (ΔT) and the proportional change in precipitation between the two time periods (ΔP) for each GCM and RCP scenario combination. Lastly, we applied ΔT and ΔP to the observed 28-year climate time series to generate a future climate time series for each GCM and RCP scenario combination. These generated climate time series were used to simulate equilibrium sagebrush cover. We simulated equilibrium cover separately for each GCM and RCP scenario combination before averaging the results, but we show the average projected climate changes across all models in Table 1.
For the temporally explicit forecasts, we used yearly GCM projections from 2012 to 2098 to simulate the process starting from the end point of the remotely sensed sagebrush cover data (ends in 2011). We aggregated daily GCM output for each GCM and RCP scenario into the seasonal climate covariates used to fit our model. These yearly climate time series were not aggregated further because we ran simulations for each GCM and RCP scenario, rather than one simulation per RCP scenario averaged over GCMs. The key assumption of our forecasting approach is that the historical correlations between weather and sagebrush cover change will continue to hold in the future.

results
Averaging across all GCMs, precipitation and temperature in our study area are projected to increase; the magnitude of increase depends on the RCP scenario (Table 1). Trajectories of our climate covariates from GCM projections show similar trends (Fig. 2).
All parameters in our model converged on stable posterior distributions (Appendix S4). Only the lagPpt climate covariate can be considered important based on a 90% credible interval, and it had a positive effect on sagebrush percent cover change (Fig. 3). in other words, if the year 2000 water year was wetter than average, sagebrush cover would increase from the 2001 to the 2002 growing season. Other climate effects strongly overlapped zero, but their posterior means were positive, except for fall-through-spring precipitation the first year of a cover transition (t − 1), whose posterior mean was negative (Fig. 3). The posterior mean for the spatial random effect, η, captured the overall spatial structure of the observed data (Appendix S5: Fig. S1). This indicates our choice of knot placement and dimension reduction strategy was adequate for describing permanent spatial variation in the data.
When we simulated the pixel-based population model based on observed climate, it was able to reproduce the spatial pattern of observed percent cover, averaged over time (Fig. 4A, B). Our model shows a tendency to underpredict perennially low-percent cover pixels (Fig. 4C), but does a better job at predicting high-cover pixels. Point predictions are most confident, although slightly biased, in low-percent cover pixels (Fig. 4D). The model is also able to adequately reproduce observed dynamics when  we make one-step-ahead predictions based on observed climate and cover in the previous year for each pixel. When we made these in-sample, one-step-ahead forecasts, the model achieved an RMSE = 4.31, in units of percent cover. The Pearson correlation coefficient between observations and predictions was 0.62.
When we apply the fitted model to iPCC climate change scenarios, the model predicts gains in sagebrush percent cover, on average (Figs. 5 and 6A). The spatial effect remains strong enough in low-cover regions to counteract the positive   effect of projected precipitation increases (Fig. 5). Thus, our model predicts an increase in the heterogeneity of sagebrush cover because projected cover increases are smaller in low-cover pixels than in high-cover pixels ( Fig. 5; Appendix S6: Fig. S1). For the temporally explicit forecasts, we show spatially averaged values and the associated uncertainty due to variability in GCM projections, variability in model parameters, and uncertainty in our process model (Fig. 6A). Based on our model and GCM projections, we forecast an average increase in sagebrush cover at our study area, but a decrease is not outside the realm of possibility (shaded regions in Fig. 6A). The generally increasing trend reflects the positive effect of precipitation on sagebrush cover change estimated for our study area (Fig. 3). We also show how our model is capable of near-term forecasts in Fig. 6B.

dIscussIon
Despite the need to forecast population responses to climate change over large spatial extents, as demonstrated by the wide application of SDMs (e.g., Clark et al. 2014), landscape-scale population models for plant species remain more concept than reality (Schurr et al. 2012, Merow et al. 2014). We introduced a new approach that uses methods from the dynamic spatiotemporal modeling literature (e.g., Conn et al. 2015) to fit a population model to remotely sensed estimated of plant percent cover. As a proof of concept, we applied our approach to a remotely sensed data product of sagebrush percent cover from 1984 to 2011 in Wyoming . We first discuss our results specific to sagebrush ecology and response to climate and then discuss the more general implications and limitations of our proposed approach.

Sagebrush response to climate and climate change
The climate effects we estimated, based on cover data at 30 m spatial resolution, are consistent with individual-level responses of sagebrush to climate-related variables. Research on individual plants has shown that wetter winters are correlated with greater stem growth in sagebrush (Poore et al. 2009, Apodaca 2013) and that warmer spring temperatures may enhance sagebrush growth in cold climates by advancing the date of snowmelt and increasing the length of the growing season (Perfors et al. 2003, Harte et al. 2015. in agreement with those individual-level responses, posterior means for all precipitation and temperature effects in our model were positive, except for the effect of fall-through-spring precipitation in the first year of a cover transition (ppt1 ; Fig. 3). The cumulative amount of precipitation the year before a cover transition (pptLag in our model) emerged as the strongest predictor of sagebrush cover change (Fig. 3). However, mean estimates for the climate effects are relatively weak (Fig. 3).
Such small effects could indicate that sagebrush are not very sensitive to interannual climate variability, that our model is poorly specified, or that climate responses are difficult to detect using coarse-scale data. Given findings from previous research demonstrating the importance of precipitation and temperature to sagebrush growth (Pechanec et al. 1937, Schlaepfer et al. 2011, Germino and Reinhardt 2014 and regeneration (Schlaepfer et al. 2014b), it is unlikely that sagebrush are insensitive to climate. We used aggregated climate covariates that may not completely capture the climate dependence of sagebrush cover change. However, the covariates we chose closely match the climate-related variables that have been shown to drive sagebrush growth, survival, and regeneration (e.g., Dalgleish et al. 2011, Schlaepfer et al. 2014b). More likely, aggregated estimates of plant abundance, such as percent cover, mask interannual variability at the level of the individual plant and make it more difficult to detect the drivers of interannual variability. Additionally, we chose not to downscale the Daymet weather data, meaning that in a given year all pixels shared the same climate, which limits our statistical power. nonetheless, our model was capable of detecting climate effects that agree with our knowledge of sagebrush ecology and allowed us to make forecasts of future sagebrush abundance.
Under projected climate, we forecast modest increases in sagebrush cover for all RCP scenarios in the long term (Figs. 5 and 6A). Our forecasts reflect both the estimated effect size for each climate covariate and the amount of change in those covariates projected by the GCMs. Cumulative precipitation the year before a given year-to-year transition was the strongest standardized effect (Fig. 3), but precipitation is projected to increase only moderately (Table 1, Fig. 2) and the negative effect of fall-through-spring precipitation in the first year of a cover transition (ppt1) had an offsetting effect. in contrast, mean spring temperature had a weak positive effect on sagebrush cover changes, but the projected temperature increase is large (Table 1, Fig. 2).
An interesting consequence of explicitly modeling the effect of space (through η) is the forecasted increase in spatial heterogeneity (Appendix S6: Fig. S1). Our model projects little change in lowcover pixels but substantial increases in the cover of high-cover pixels (Fig. 5). Had we not explicitly accounted for spatial dependence in our model, we would have missed this result. We were unable to attribute the spatial structure apparent in the data (Fig. 4A) and approximated by our model (η; Appendix S5: Fig. S1) to slope, aspect, elevation, or coarse soil type (results not shown). The lack of correlation between η and landscape factors leads us to conclude that the spatial structure in our data set emerges from some combination of fine-scale microhabitat associations and legacy effects of disturbance.
While we forecast an increase in sagebrush cover at our study area, SDM studies typically project dramatic declines in climate suitability for sagebrush with warming (Shafer et al. 2001, neilson et al. 2005, Bradley 2010, Schlaepfer et al. 2012, Still and Richardson 2015. There are many potential explanations for this apparent contrast, ranging from the type of model used to the particular climate covariates considered, but the location of our study area in a cold portion of sagebrush's geographic distribution may be the best. The response of plant species to weather varies along climatic gradients (e.g., Clark et al. 2011, Vanderwel et al. 2013, and sagebrush are especially sensitive to the timing of snowmelt because their growth depends on recharge of deep soil water (Schlaepfer et al. 2012(Schlaepfer et al. , 2014a. in warmer parts of the sagebrush range, earlier snowmelt is detrimental to growth and survival (Pechanec et al. 1937, Schlaepfer et al. 2011, Germino and Reinhardt 2014. in colder regions, earlier snowmelt due to temperature increases can lengthen the growing season and increase sagebrush occurrence and cover (Schlaepfer v www.esajournals.org TREDEnniCK ET AL. et al. 2012, 2014a. The average annual temperature across the sagebrush steppe biome is 6.9°C (SD = 1.6; Schlaepfer et al. 2011), whereas average temperature at our study area from 1980 to 2013 was 4.6°C (calculated from Daymet estimates). Our study area lies at the cold extreme of the sagebrush range; thus, the weak positive response to temperature that we estimated (Fig. 3) and carried through to our forecasts (Figs. 5 and 6A) likely represents the positive effect of earlier snowmelt, and thus higher moisture availability early in the growing season.
A previous analysis of a different subset of the remote sensing data set we used also came to a different conclusion, projecting future sagebrush decline (Homer et al. 2015). The discrepancy between the results of Homer et al. (2015) and ours primarily reflects a difference in the climate projections used for projecting future changes rather than differences in our inference about responses to historical variation in weather. Homer et al. (2015) used downscaled weather projections from a single model from the iPCC 4, whereas we used native-resolution weather projections from a suite of models from the iPCC 5. Consistent with our study, Homer et al. (2015) found a generally positive relationship between pixel-level sagebrush cover and precipitation, but the future climate scenario they chose resulted in a mean decrease in precipitation, causing a predicted decline in sagebrush cover. A second difference is that Homer et al. (2015) relied on regressions of decadal trends in sagebrush cover against decadal trends in climate at the level of individual pixels. Our current approach is fundamentally different in that we specifically model the impact of interannual variation in weather on year-to-year changes in sagebrush cover using a dynamic population model. Thus, our model takes advantage of the additional information contained within short-term responses to climate fluctuations. Lastly, the location of Homer et al.'s (2015) study area is, on average, at a lower elevation than our current study area. The geographic difference results in different historical and projected climate, and, as discussed above, sagebrush may respond differently to warming depending on geographic location.
We projected sagebrush cover to the end of this century, but an important feature of our approach is that it can also produce short-term forecasts (Fig. 6B). For example, we could forecast the effects of a multiyear regional drought on sagebrush cover (Debinski et al. 2010). Validating spatial population models against short-term predictions would give ecological forecasters a way to assess and improve the performance of their models, which would greatly increase our confidence in long-term forecasts. This cycle of prediction, validation, and refinement is missing from most currently available population-level forecasts of the effects of climate change.

A landscape-scale plant population modeling approach: opportunities and limitations
Our approach for modeling plant populations overcomes two major hurdles for spatially explicit population models. First, we used moderate-resolution, remotely sensed estimates of sagebrush percent cover as a response variable, enabling us to fit a dynamic population model over a large spatial extent. Speciesspecific estimates of plant abundance are becoming commonplace as remote sensing technology develops (e.g., Baldeck andAsner 2014, Colgan and, and in a few years, several remotely sensed time series may be available. Second, borrowing from new methods in spatiotemporal modeling of animal abundance (e.g., Conn et al. 2015), we fit the model using a dimension reduction strategy that accounted for spatial autocorrelation within a feasible computational time. Accounting for spatial autocorrelation allows for statistically rigorous inference on the effects of interannual climate on sagebrush cover change in our study region. The spatial covariance structure also provided a way to obtain spatially explicit predictions at a resolution below that of the climate covariates (i.e., within the study region; Figs. 4 and 5). Our approach is amenable to any spatially explicit time series of plant abundance, but we see remote sensing data sets offering the largest opportunity for landscape-scale population models. Furthermore, it would be straightforward to include additional covariates related to disturbance (e.g., fire) or biotic interactions. Thus, we see our method as a first step toward coupling the mechanistic power of dynamic population models with the spatial extent of SDMs. The spatially and temporally explicit forecasts made possible by our approach should v www.esajournals.org TREDEnniCK ET AL. be especially relevant to land management decisions based on near-term forecasts.
Several a priori modeling decisions determined the spatial extent and resolution of our results. We retained the native spatial resolution of the remote sensing data (30 × 30 m). This constrained the extent that we could reasonably model because of the computational challenges in estimating spatial random effects. Even with our dimension reduction technique, modeling a larger area at this resolution would require a greater number of spatial knots, and computation time would increase substantially (Wikle 2010). To model a larger spatial extent, we could aggregate the original remote sensing time series data to a coarser spatial resolution. This would allow us to model a much greater spatial extent with a similar number of knots and a similar computation time. While a coarser scale model would lose some fine-scale detail, it could be applied to a much larger area, potentially gaining some strength in estimating climate effects by spanning a greater range of climate variation. However, gains made by incorporating greater regional variability by modeling at a coarser resolution could be offset by the loss of information inherent when aggregating plant responses into larger pixels.
Our spatial extent and resolution also affected our use of climate covariates. We did not downscale Daymet data to match the spatial resolution of the sagebrush data, meaning that in each year all pixels share the same climate covariates. This is a potential limitation of our study, and could explain the weak effect of climate covariates that we observed (Fig. 3). We also did not allow different portions of our study area to respond to climate in different ways. Doing so would require spatially varying climate effects and a substantial increase in computational time. However, in future applications, it will be important to allow climate effects to vary over space to better capture reality. Conn et al. (2015) provide examples of how such spatiotemporal interactions can be included in abundance models. We might expect climate effects to interact with spatial covariates such as soil type, slope, and aspect. in our relatively small study area, we did not observe important effects of these factors, but it is possible to include such abiotic data layers as predictors when fitting models at larger spatial extents where variability may be greater.
The uncertainty associated with our forecasts highlights several opportunities to improve our approach. First, parameter uncertainty could be reduced by regulating the variance of the posterior distributions of climate covariates via ridge regression (e.g., Gerber et al. 2015). Second, uncertainty associated with climate projections could be reduced by identifying GCMs that perform exceptionally well for a particular study location (e.g., Rupp et al. 2013). Such considerations will be important when forecasting in support of particular management objectives. However, knowledge of uncertainty is itself important knowledge for management (Bradshaw and Borchers 2000). Deciding that no actions should be taken based on the data at hand is itself a management decision.

conclusIon
We introduced a new approach to fitting and simulating population models at large spatial extents with plant population data derived from state-of-the-art remote sensing. We used the model to forecast future abundances of sagebrush in Wyoming and found that at our relatively cold site sagebrush should be expected to increase in cover. As more species-level remote sensing data sets become available and computing power increases, this approach will be applicable to a wider number of species and even larger spatial extents. Future modeling could include the effects of nonclimate driversincluding the effects of species interactions and disturbance. For sagebrush, including fire and competition with nonnative annual grasses in the model may be especially important for a complete assessment of the effects of climate change (Bradford and Lauenroth 2006). Fortu nately, our spatiotemporal modeling framework could easily be extended to model additional species and dynamic processes as the data become available. The approach we have developed here fills an important gap in spatial scales between SDMs and local-scale demographic population models.

acknowledgMents
This work is the outcome of a distributed graduate seminar led by PBA and supported by a national Science Foundation CAREER award (DEB-1054040).
David T. iles, Eric LaMalfa, and Rebecca Mann participated in project conception as part of the distributed graduate seminar and provided comments that improved the manuscript. ATT was supported by an nSF Postdoctoral Research Fellowship in Biology (DBi-1400370), and AK was supported by an nSF Graduate Research Fellowship. Additional support came from the Utah Agricultural Experiment Station, Utah State University, and this article is approved as journal paper number 8856. We are grateful to Debra K. Meyer at USGS EROS for extracting the data set used in this study and to David Koons and two anonymous reviewers for comments that improved the manuscript. Compute, storage, and other resources from the Division of Research Computing in the Office of Research and Graduate Studies at Utah State University are gratefully acknowledged. We acknowledge the World Climate Research Programme's Working Group on Coupled Modelling, which is responsible for CMiP, and we thank the climate modeling groups (listed in Appendix S1: Table S1) for producing and making available their model output. For CMiP, the U.S. Department of Energy's Program for Climate Model Diagnosis and intercomparison provides coordinating support and led development of software infrastructure in partnership with the Global Organization for Earth System Science Portals. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government. lIterature cIted