Samples and Data Collection

The climate model was developed from presence and absence points. Presence data, consisting of 131 occurrence points (Appendix 1), were derived principally from previous studies: McArthur and Sanderson (1999), Richardson et al. (2012), and Wilt et al. (1992). Techniques used to determine subspecies are described within each publication, but in nearly all cases flow cytometry or chromosome counts were used to confirm ploidy. The exception is Wilt et al. (1992), who use morphology and an assessment of phenolic compounds. Absence data was derived through several sources and contained 4464 points consisting of both target-group absences (TGA) and randomly selected background points. TGA, localities for other taxa that do not co-occur with the target species, have been used successfully in species distribution modeling (Mateo et al. 2010). A total of 3964 TGA were derived from previous studies (Richardson and Meyer 2012; Esque et al. unpubl. data), the USDA Forest Service, Forest Inventory and Analysis Program (FIA) (Bechtold and Patterson 2005), and the Consortium of California Herbaria (CCH) (data provided by the participants of the Consortium of California Herbaria [ucjeps.berkeley.edu/consortium]). The taxa used in TGA are listed in Appendix 2, along with their respective sources. A total of 500 background points were randomly selected from a group of three Level III Ecoregions (Omernick 1987). These ecoregions, Nebraska Sandhills, Northwestern Glaciated Plains, and Southwest Tablelands were chosen to provide additional absence points in areas lacking TGA to fill out the range of climatic variation.

Climatic Data

The geographic extent for both models and projections was set from 30° N to 55° N latitude and from 130° W to 100° W longitude to incorporate the entire range of possible sagebrush habitat. The baseline climatic data set was acquired from WorldClim (Hijmans et al. 2005), comprising 19 bioclimatic (BIOCLIM) variables (Appendix 3) for present conditions (mean 1950-mean 2000) at 30 arc-second resolution, which is roughly 1 km² at the equator. The BIOCLIM variables have been widely used in modeling work as variables that are biologically important for various species (Hijmans et al. 2005).

Bioclimatic Model

To model the climate-defined area of Wyoming big sagebrush, we estimated the likelihood that the climate was suitable across a large section of western North America. The estimate was derived from a climate profile, which is a multivariate description of the climatic niche. The climate profile was developed from bioclimatic models, that is, regressions of the presence and absence of a species on climate variables. The modeling techniques used here closely follow those of Rehfeldt et al. (2006) as explained in detail in Rehfeldt et al. (2009) and Crookston et al. (2010).

The Random Forests classification tree of Breiman (2001), implemented in R 3.02 (R Core Team 2013) by Liaw and Wiener (2012) in the package “RandomForest,” was used to predict the presence or absence of species from the climate variables. The Random Forests algorithm constructs a set of classification trees from an input data set and outputs statistics that reflect the likelihood that the climate at a location is suitable for the species (Rehfeldt et al. 2009). The trees in aggregate are called a forest. The climate profile was built on 12 forests, each with 100 trees (i.e., decision trees).

To create the trees, a majority of the data, usually around 64%, were used to create the model, and the remaining portion of the data set, the out-of-bag occurrence points, were used to test the model. The best-fitting model for each tree was built by comparing the out-of-bag error. Out-of-bag errors are comprised of rates of errors in commission (where the model predicts an occurrence when no plant is present), and errors in omission (where the model predicts an absence when a species is actually present). To make predictions about presence of a species, each tree in the forest provides one vote to the classification of an observation. Because classification errors approach a limit as the number of trees in the forest increase, collinearity and over-parameterization are inconsequential (Breiman 2001). The approach has been shown to be robust and has worked for widely distributed species (Ledig et al. 2010).

Assembling the presence-absence data for analysis requires satisfying Breiman's recommendation that presence data be in reasonable balance with absence data (Breiman 2001). Each forest would need one data set, and each data set was prepared within which presence and absence points represented 40% and 60% of the total, respectively. For each of the data sets, the amount of presences was fixed at 40% to limit the amount of out-of-bag errors, which increase when the number of presence points is less than 40% of total points used in the model (Rehfeldt et al. 2006).

All data sets contained all 131 presence points, each of which was weighted by a factor of two (each was included twice). This weighting assures that the resulting model is most robust for climates in which A. tridentata subsp. wyomingensis actually occurs (Rehfeldt et al. 2006; Ledig et al. 2010), and allows the number of absence points in the data set to be doubled, allowing for more complete sampling of the climatic variation. Each data set, therefore, included about 655 observations, with 262 observations with sagebrush, and about 393 observations without sagebrush.

Absence points for the data sets were chosen in two steps. First, following the protocol of Rehfeldt et al. (2006), we defined an expanded climatic envelope as a 19-variable hypervolume corresponding to the climatic limits of distribution expanded by ±1 SD. Then, for each data set we randomly selected 40% of the points as absences from points that are within, and 20% of points were chosen randomly as absences from points outside, the climatic hypervolume detailed above. The number of forests was chosen by dividing the total number of absence points within the climate hypervolume described above by the number of presence points multiplied by two. Therefore, using 12 forests would assure that the probability would be high that all observations within the hypervolume would be used in at least one forest.

The final predictor variables used were culled from the 19 BIOCLIM variables through a variable reduction process following Rehfeldt et al. (2006) and Rehfeldt et al. (2009). Based upon the out-of-bag error, the predictors were eliminated using the mean decrease in accuracy to judge variable importance until only one variable remained. Then the top seven variables were chosen to use as predictive variables to create the climate profile.

Mapping

The climate profile from Random Forests analysis was mapped to the WorldClim climate grids. Each of the contemporary grid cells was evaluated for climatic suitability for sagebrush by the number of votes cast for the 100 trees in the 12 forests. A grid cell was considered to have suitable climate for sagebrush when the majority of the 1200 votes were cast in favor of the climate being suitable for sagebrush. This creates the bioclimatic model. Models were evaluated using the Area Under the Curve of Receiver-Operating Characteristic (AUC), a common measure for evaluating model fitness (Elith and Leathwick 2009).

The climatic data sets for the 2050s were acquired from WorldClim (Hijmans et al. 2005) and comprise the same bioclimatic variables as the contemporary data set. Climate surfaces for the 2050s (2040–2069) (Hijmans et al. 2005), derived from the IPCC (Intergovernmental Panel on Climate Change) 4th Assessment (IPCC 2007), were used to project the sagebrush bioclimate for this decade. To provide a consensus of 2050s projections, we used methods similar to Ledig et al. (2012) and Wang et al. (2012), where the outputs from five General Circulation Models (GCMs) are combined into an agreement map. GCMs included the Alb emission scenarios for the following five models: Canadian Center for Climate Modeling and Analysis (CCCMA CGCM3.1); Bjerkes Centre for Climate Research Norway (BCCR BCM2.0); Institute for Numerical Mathematics, Russia (INM-CM3.0); Commonwealth Scientific and Industrial Research Organization (CSIRO MK3.0); and the Center for Climate System Research (University of Tokyo), National Institute for Environmental Studies, and Frontier Research Center for Global Change (JAMSTEC), Japan (MIROC3.2 medres). Information on the GCMs and emission scenarios can be found elsewhere (IPCC 2007). Agreement mapping of the five GCM-scenario combinations were performed in R using the RandomForest package as above. The threshold used to calculate suitable area for the contemporary and each of the future models was 0.5. For the 2050s, the predicted presence of sagebrush-suitable climate is mapped only where more than two of the five GCMs showed agreement.

Ecoregional Assessment of Climate Niche Loss

For both contemporary and future (2050s) projections of the bioclimatic model, the total area (km²) predicted to have suitable climate was calculated. We also calculated the area where both contemporary and future models overlapped (stable), the area that is suitable in the contemporary model but not suitable in the future model (contracting), and the area that is not suitable in the contemporary model but is suitable in the future model (expanding).

Climate comparisons were made between two geographic regions that are predicted to have the greatest losses in climate niche: the Great Basin and the Great Plains. We defined the regions by a combination of Omernik's (1987) Level III Ecoregions. The Great Basin region is defined here by three of the ecoregions: Northern Basin and Range, Central Basin and Range, and Snake River Plain. The Great Plains region is here defined by several of the ecoregions: Middle Rockies, Southern Rockies, Northwestern Great Plains, Nebraska Sand Hills, High Plains, and Southwest Table Lands. The contemporary climate niche of Wyoming big sagebrush was split by whether the area is predicted to contract or remain stable by mid century. For these geographic areas, we compared Annual Dryness Index (ADI; mean annual precipitation / degree days >5 °C), as calculated in Rehfeldt et al. (2006), and summer-winter precipitation ratio (SWP; warmest quarter precipitation [PWQ, BIO 18] / coldest quarter precipitation [PCQ, BI019]). These values were extracted in each raster grid cell for the Great Basin and Great Plains. Two-way ANOVA and Tukey's HSD was used to test post-hoc mean differences between ecoregions and stable and contracting areas within each region.

Bioclimatic Model

The average model AUC was excellent at 0.979. To balance commission, omission, and out-of-bag error, we chose a seven variable model with an out-of-bag error of 2.65%. Commission and omission rates for this model were 4.51% and 0.03%, respectively (Table 1). The climatic data set comprised seven bioclimatic variables (Table 2). The seven-variable model was chosen as being reasonably parsimonious while providing a buffer against reliance on single variables. The most important variable for this model was the mean temperature for the warmest quarter. The second most important variable was the annual mean temperature, and the third most important variable was temperature annual range. Of the seven variables used in the model, six were related to temperature. Mean annual precipitation was the only precipitation related variable and was the sixth most important.

Table 1.

Confusion matrix of Wyoming big sagebrush (Artemisia tridentata subsp. wyomingensis) bioclimatic model showing the class error and number of observations classified by the Random Forests algorithm. Mean of the 12 Random Forests.

Mapped Projections

The contemporary climate niche predicts an area of nearly 108 million hectares (1,086,697 km², Figure 1A) for Wyoming big sagebrush. By midcentury, a 39% reduction is predicted in this climate niche, totaling 66 million hectares (Figure 1B). Only 32% of the contemporary climate niche is stable by the middle of the century, while 67% of the contemporary climate niche is predicted to be lost and 28% will be gained. Regions predicted to be most vulnerable to climate change extirpation include the trailing edge (i.e., the southern periphery of the subspecies), the western Great Plains, and lower elevations of the Columbia and Great Basin. Regions that retain or gain climate niche include western Wyoming and eastern Idaho, higher elevations in the Great Basin and the northern Great Plains.

Table 2.

Climate variables used to predict the climate niche of Wyoming big sagebrush (.Artemisia tridentata subsp. wyomingensis) bioclimatic model. Bioclimatic model variables are listed in order of importance.

Ecoregional Assessment of Climate Niche Loss

The range of climatic conditions affecting the predicted loss of Wyoming big sagebrush differed among spatial and temporal scales (Appendix 4). The SWP and ADI were significantly different between ecoregions (Great Basin versus Great Plains; P < 0.0001) and within ecoregions between predicted stable and contracting areas. For both ecoregions, the stable areas of the climate niche had a lower ADI than the contracting areas (P < 0.0001). ADI values and differences between stable and contracting areas were much greater in the Great Basin than Great Plains (P < 0.0001) (Figure 2A and C). While small, the difference between stable and contracting areas of the Great Basin was significantly different for SWP (P < 0.0001), and significantly large differences in SWP were observed between stable and contracting areas (P < 0.0001) within the Great Plains (Figure 2B and D).

Figure 1.

A) Mapped projection of the contemporary climate niche of Artemsia tridentata subsp. wyomingensis. Dark brown represents higher probability of occurrence (>0.75), whereas light brown represents lower probability (0.5 to 0.75). B) Mapped projection of the change in climate niche between contemporary and decade 2050. Dark purple represents areas that are predicted to have suitable climate for this subspecies in decade 2050 (i.e., stable or expanded), whereas light purple represents areas that are predicted to have unsuitable climates (i.e., contracted).

Model Development and Error

It is well known that big sagebrush subspecies are defined by climate and that Wyoming big sagebrush occupies the warmest and driest extent of this species distribution (Mahalovich and McArthur 2004). To be successful in the restoration of sagebrush ecosystems, it is imperative that subspecies are placed in the appropriate climate. Previous published bioclimatic models have used broader taxonomic hierarchies based on the constraints of GAP analysis data to define sagebrush climate niche (Bradley 2010; Schlaepher et al. 2012). In this study, our goals were to produce a management tool for contemporary and future restoration of Wyoming big sagebrush. Data were acquired from known occurrences that span much of the range of the subspecies and have been taxonomically identified to subspecies. Our strategy was to frame this subspecies' climate using targeted group absences using allopatric species in warmer and cooler climates (i.e., Coleogyne ramosissima Torr, and Cercocarpus ledifolius Nutt., respectively). A targeted group absence approach has been shown to be more accurate than pseudo-absences (Mateo et al. 2010). However, along the central and northern Great Plains, pseudoabsences were necessary because of the lack of suitable species to use as absence points.

As with any bioclimatic model, some modeling error can be expected. Sources of errors could come from the environment, including soil and small-scale topographic features. Ecological interactions (e.g., plant competition), disturbance, and land use histories could also be sources of error. Nevertheless, the resulting model generated low errors in the prediction of which much were due to commission. Geographically, we suspect the preponderance of commission errors occurs along the Great Plains and the boundary with the Chihuahuan Desert in New Mexico (Figure 1A). In these regions, Wyoming big sagebrush is less dominant. This is likely due to a change to increasing summer versus winter precipitation that favors grasslands (Ogle and Reynolds 2004). Here, soils and topography become more of an important influence on presence and absence, and thus predictive error based on a climate-only model.

Figure 2.

Boxplots illustrating the range of values for Annual Dryness Index (ADI) and summer-winter precipitation ratio (SWP) for stable and contracting areas of the contemporary climate niche of Wyoming big sagebrush (Artemisia tridentata subsp. wyomingensis). Increasing values of ADI indicate decreasing precipitation and / or increasing accumulating temperatures >5 °C. Increasing values of SWP indicate a higher ratio of summer relative to winter precipitation.

Contemporary and Future Projections

Comparisons of contemporary and future projections show considerable loss (39%) of Wyoming big sagebrush climate niche good correspondence to the previously published bioclimatic modeling of North American-biomes (Rehfeldt et al. 2012) and blackbrush (Richardson et al. 2014), an ecotonal species occurring between warm and cold deserts. Rehfeldt et al. (2012) showed that midcentury Mojave Desert climates would replace cold desert biomes in some areas of the Great and Columbia Basins. These areas that show biome turnover from cold deserts to warm deserts (i.e., Great Basin Scrub to Mojave) are also areas that show major reductions in Wyoming big sagebrush climate niche (this study), the most prevalence of cheatgrass (Bradley 2010), and areas estimated to have very low restoration potential (Wisdom et al. 2005). Moreover, midcentury projections of blackbrush climate niche show expansion into contemporary Wyoming big sagebrush climate niche in the Lahontan and Columbia Basins and Lower Snake River Plain (Richardson et al. 2014).

Understanding the biological association between climate niche loss and life history traits of the target species is an important aspect of ecological and adaptive processes. While we do not have direct data supporting which life histories traits are critical to success or failure of this subspecies, previous research provides opportunity to speculate. The success of big sagebrush seedling establishment has been shown to be dependent on the timing and amount of precipitation. Snowpack appears to be a critical component for big sagebrush seedling recruitment. Studies have shown that snowdrifts, either caused by other plants or fencing, can greatly increase the recruitment of seedlings (reviewed in Meyer 1994). Another component of climate that affects the distribution of big sagebrush is the seasonality of precipitation. As discussed above, predominant summer precipitation favors grasslands, whereas winter precipitation favors shrublands (Ogle and Reynolds 2004; Brooks and Chambers 2011). Changes in climate that affect the longevity of snowpack and the seasonality of precipitation could greatly impact big sagebrush geographic distribution.

In this study, we examined the changes in two regions that support the highest predicted midcentury loss of Wyoming big sagebrush, the western Great Basin and the northern and central Great Plains. Based on these analyses, the climate conditions that result in the loss of Wyoming big sagebrush are different between the two regions (Figure 2). An interaction between increasing summer temperatures and reduced precipitation (ADI) appears to be an important component to climate niche loss in the Great Basin. Greater aridity differences were observed between stable and contracting areas. Differences were considerably smaller between stable and contracting areas in the Great Plains (Figure 2A and C). In contrast, the seasonality of precipitation events from less winter to more summer is expected in the Great Plains, whereas relatively minimal change is expected in the Great Basin (Figure 2B and D). Grassland ecosystems would likely prevail in the western Great Plains based on these projections.

As historical plant migration rates have been estimated to be 10–30 km per century (McLachlan et al. 2005; Yansa 2006), it is likely that the net loss will actually be more than 39% as a large portion of the expanding area is more than 30 km from current localities. Therefore, the species may not be able to expand into the new suitable range in the short period of time (ca. 30 to 40 years). If sagebrush is unable to colonize the expanded areas of niche in the shortened window of climate change, assisted migration is one possible solution to the problem (Ying and Yanchuk 2006; Kramer and Havens 2009; Vitt et al. 2010). In Havens et al. (this issue), assisted migration is defined as “the purposeful movement of individuals or propagules of a species to facilitate or mimic natural range expansion or long distance gene flow within the current range, as a direct management response to climate change.” Successful assisted migration would need to ensure that the correct plant sources are transferred to the appropriate area. Such research in understanding the adaptive variation in Wyoming big sagebrush and other subspecies is ongoing.

Management Strategies and Planning

Restoration of Wyoming big sagebrush is a difficult and complex task. Restorationists will have to utilize a variety of management options and weigh a number of potential variables that can affect conservation and restoration outcomes (Chambers et al. 2013). Given the limited resources available, managers will have to focus on restoration sites that meet the most criteria for successful outcomes. Central among these criteria is an understanding of the impact of climate change. Our modeling focuses on the subspecies of the big sagebrush complex that is the most widespread and occupies the warmest and driest niche. The model suggests areas predicted to have an unsuitable climate niche in the upcoming decades (Figure 1B) would be poor choices for restoration of Wyoming big sagebrush; however, seed collected in these regions would be desirable for ex situ conservation or transfer to nearby suitable climates. Restoration should be focused on areas that are predicted to sustain Wyoming big sagebrush or areas of expansion.