The importance of seasonal resource selection when managing a threatened species: targeting conservation actions within critical habitat designations for the Gunnison sage-grouse

Additional keywords: Centrocercus minimus, Colorado, critical habitat, Gunnison sage-grouse, resource selection, species distribution.

The conservation of species-at-risk and the habitats they rely upon is a concern for resource managers. The ability to identify these priority habitats is an important and urgent first step in most conservation strategies (Johnson et al. 2004; Fedy et al. 2014). Identifying high quality habitat can be used to support conservation decisions regarding invasive species risk assessment, critical habitat designation, property acquisition and translocation of threatened or endangered species, but the use of these decisions in guiding conservation policies is still scarce (Guisan et al. 2013).

When species are listed under the US Endangered Species Act (ESA) 1973, there is a requirement to evaluate ‘critical habitat’ (Greenwald et al. 2012). Critical habitat is defined by the US Fish and Wildlife Service (USFWS) as the ‘specific areas within the geographic area, occupied by the species at the time it was listed, that contain the physical or biological features that are essential to the conservation of endangered and threatened species and that may need special management or protection’ (USFWS 2015). The designation of critical habitat requires the use of the best scientific data available (Kalen 2014; Murphy and Weiland 2016), but most critical habitat designations are identified based on known locations or general habitat features (Camaclang et al. 2015). Relatively few critical habitat designations are identified using habitat suitability models or knowledge of spatial structure, because it requires additional data (Camaclang et al. 2015). As a result, a critical habitat designation without using habitat suitability models, when the data is available, can lead to broad and sweeping geographic designations that have little utility in siting on-the-ground conservation efforts because they are neither ‘specific’ to the species nor occupied.

Gunnison sage-grouse (Centrocercus minimus; hereafter GUSG) are considered threatened under the ESA and critical habitat has been designated (USFWS 2014a). The distribution of GUSG has declined to an estimated 10.3% of its historical range (Schroeder et al. 2004). Studies have linked the grouse’s decline to habitat loss, degradation and fragmentation of sagebrush landscapes (Oyler-McCance et al. 2001; GSRSC 2005; Bukowski and Baker 2013). Current GUSG were recognised as a separate species from Greater sage-grouse (C. urophasianus), with the primary differences including size, plumage, courtship display and genetics (Young et al. 2000). GUSG distribution is limited to seven isolated populations in south-west Colorado and one in south-eastern Utah (Schroeder et al. 2004). Most of these populations are considered satellite except for one population, the Gunnison Basin (GB), which contains 85–90% of all GUSG rangewide, ~63% of occupied habitat and ~60% of leks (Fig. 1; GSRSC 2005; USFWS 2014b).

Fig. 1.  The Gunnison Basin population in south-west Colorado, USA along with the satellite populations (Crawford, Pinon Mesa, Cerro/Cimmarron/Sims, San Miguel, Dove Creek, and Poncha Pass) and the area in the Gunnison Basin covered by each dataset including the original data (2004–10), validation dataset A (2002), validation dataset B (2010–11), and active lek locations. Validation dataset B is the only dataset that includes the western portion of the Gunnison Basin.

The current critical habitat designation outlined by the USFWS uses the occupied habitat (called the species area mapping) defined by Colorado Parks and Wildlife (CPW) for the GUSG rangewide conservation plan (GSRSC 2005). Although this was the best available science in 2005, these maps were a rudimentary outline of the geographic range of GUSG based on vegetation cover types that GUSG could occupy, according to professional judgement and observation, and were mostly based on breeding season (i.e. lek locations) occupancy. There is a basic assumption in the range delineation that within the GUSG range all the plant communities are equally valuable; this is a potentially misleading oversimplification that ignores seasonal use differences and the landscape diversity that can influence GUSG movements and habitat use (Connelly et al. 2011). The GSRSC (2005) acknowledged these limitations and developed specific conservation strategies encouraging seasonal habitat mapping efforts to better identify areas to prioritise for protection, because the current occupied range may not reflect species behaviour or movement patterns (GSRSC 2005). Ultimately, the USFWS used this occupied range boundary as the basis for critical habitat designation and refined critical habitat further by excluding existing candidate conservation agreements with assurances, conservation easements and other management agreements with conservation measures applicable to GUSG. These exclusions reduced critical habitat by ~15% (USFWS 2014b), but many of these excluded areas may or may not be important GUSG habitat. Their exclusion from designation is one reason why the current critical habitat designation excludes biologically relevant habitat. In addition, it is recognised that to appropriately manage GUSG populations and habitat, managers must have a better understanding of seasonal habitat use to conserve the important segments of the landscape associated with distinct seasonal habitats throughout their annual cycle, regardless of land ownership (Connelly et al. 2000; Connelly et al. 2011; Fedy et al. 2014). The critical habitat designation by the USFWS acknowledges local scale seasonal habitat structure for GUSG, but does not (and cannot) use that information in designating critical habitat (USFWS 2014b), thereby ignoring the seasonal movement patterns of GUSG on the landscape.

The USFWS designates ~84% of the GB as critical habitat, but there is no differentiation of spatial or temporal variability in quality and use. For over a decade, fine scale models based on detailed patterns of animal use have shown potential to identify crucial habitat not obvious at larger scales (Klar et al. 2008). These more detailed seasonal models are needed for identifying restoration and rehabilitation of areas important to the recovery and viability of GUSG (Wisdom et al. 2011). For example, previous models (Aldridge et al. 2012) have identified nesting habitat for GUSG and suggest that ~50% of the GB is important nesting habitat – less than the current designated critical habitat (USFWS 2014b). In an effort to define seasonal habitat selection patterns, we collected radio-telemetry data from over 200 GUSG from 2004 to 2010, along with two independent validation datasets that could be used to model GUSG habitat selection in a more robust fashion specific to the GB during two seasons (breeding, 1 April–15 July; and summer, 16 July–30 September) and covering those 6 months of the year.

We had two main objectives for this study. First, we developed habitat selection models for the breeding and summer seasons using GUSG radio-telemetry data collected in the GB population from 2004 to 2010. Second, we compared the specific seasonal models with the critical habitat designated during the ESA threatened listing process. Because of the broad generalisation of critical habitat designation by the USFWS, we believe that our analysis will assist the USFWS in the development of a recovery plan for GUSG to more specifically focus on conservation actions in the most important areas within the GB.

Our study area was the Gunnison Basin population of GUSG in Gunnison and Saguache counties, Colorado, USA. This area comprises 85–90% of all GUSG (GSRSC 2005; Fig. 1). The Gunnison Basin is a 200 km² intermontane basin ranging in elevation from 2300 to 2900 m (Hupp and Braun 1989). Mountainous terrain areas border the north, east and south-east sides of this population and are not commonly used by GUSG.

We captured GUSG from March to early May from 2004 to 2010, using spotlighting techniques (Giesen et al. 1982; Wakkinen et al. 1992). We fit GUSG with a 16 or 17 g necklace-style VHF radio-transmitter (model A4050 by Advanced Telemetry systems, Isanti, MN, USA or model R12B by Holohil Systems, Ontario, Canada) equipped with a 4-h mortality sensor, and also fit each bird with a uniquely numbered leg band (National Band and Tag Co., Newport, KY USA). The transmitter was <2% of an average GUSG female (1270 g, s.d. 90 g) or male (2110 g, s.d. 190 g). Trapping and handling protocols were approved by the Colorado Parks and Wildlife Animal Care and Use Committee (permit # 02-2005).

Following release, we located radio-marked individuals on the ground using hand-held Yagi antennas once every 1–3 days (from date of capture through September) to monitor status (dead or alive) and movement patterns. Monthly flights were conducted during the winter to assess survival, but due to logistical restraints, on the ground telemetry locations were not collected. We used triangulation to estimate locations using maximum likelihood estimates generated in program LOCATE II (Nams VO (1990) Locate II. Pacer, Truro: Nova Scotia, Canada).

Locations were assigned to one of two seasons: breeding or summer. The cutoff dates were based on the Gunnison sage-grouse Conservation Plan, which guides GUSG management in Colorado (GSRSC 2005). We used a resource selection framework comparing used vs available locations, so we randomly generated a sample of ‘available’ locations within the same geographical extent of the GB (Stokland et al. 2011). We conducted a sensitivity analysis of the available sample size for each season (Northrup et al. 2013) and found that coefficients converged at n = 9000 available samples in both seasons. Average daily movement distances were estimated across all the birds in both seasons and used to buffer all presence and available locations within that season (180.5 m in the breeding season and 223.0 m in the summer season). Due to terrain constraints and access issues, there were a few telemetry errors that were larger than this average daily movement. Therefore, we used this buffer to remove those locations and summarised all habitat variables at a biologically relevant scale.

We classified vegetation type into eight biologically meaningful categories (available as Supplementary Material 1 on the journal website) using the Colorado basinwide vegetation layer (i.e. sagebrush (Artemisia spp.), forest, etc.). This land cover layer was constructed in 2005 from 25-m resolution landsat imagery as part of the Colorado Vegetation Classification Project administered by CPW in collaboration with the Bureau of Land Management and the USDA Forest Service. We believe that the vegetation cover type was consistent and reflected current conditions in the GB as there were few relevant environmental perturbations (e.g. wild or managed fire) before or during our study. From the eight categories, we excluded those that consisted of <0.1% of the cover within the buffers (e.g. urban, riparian and water). Thus, the vegetation categories we used in model development were irrigated agriculture (pastureland), sagebrush, grassland, bare and forest, which have all been shown to be important predictors of habitat use in previous studies on both Greater sage-grouse and GUSG (Rice et al. 2013; Stanley et al. 2015; Walker et al. 2016). In addition, these vegetation categories have been cited as influential in multiple seasons for GUSG habitat and are often the focus of management actions in Colorado (GSRSC 2005).

We obtained elevation data from the USA Geological Survey (USGS) digital elevation model and used the national hydrography dataset to measure density of water bodies (Dzialak et al. 2012). We used the National Wetlands Inventory layer from USFWS to measure the distance to wetlands. Development variables included roads, easements and address points (hereafter referred to as residential) collected by Gunnison County. Using these layers, we calculated distance to highways (paved roads only), distance to residential, distance to easements, residential density and road density (includes paved and unpaved roads). We also measured the distance to sagebrush, because observations of GRSG feeding on the edges of vegetation cover near the ecotone with sagebrush have been recorded, particularly in the summer season (Connelly et al. 2011). A list of the variables summarised by seasonal buffer is in Table 1.

Table 1.  Variables used in Gunnison sage-grouse habitat models with the mean value and standard errors for presence and available buffers in the breeding (1 April–15 July) and summer (16 July–30 September) seasons in Gunnison Basin, Colorado, USA (2004–10)

We first calculated Pearson correlation coefficients (r) for all the habitat variables (r > 0.70; McGarigal et al. 2000) to remove highly correlated variables that may cause multicollinearity issues in the model and removed the variable with a larger P-value. Table 1 indicates the final variables used in each season after removing correlated variables. We fit our used and available location data with a binomial generalised linear model using the package ‘lme4’ in program R (R Core Team 2013; Bates et al. 2013). We used a random intercept for each individual grouse within each season to account for unbalanced sampling among animals (Gillies et al. 2006).

We constructed a set of alternative models from all linear combinations of the habitat variables in each season (McAlpine et al. 2008). We generated predictions from each of the best-fitting models within the 95% confidence set and averaged the predictions into a final model for each season to strengthen inference (Burnham and Anderson 2002). The unstandardised coefficients for each of the seasonal models are equivalent to selection ratios (Manly et al. 2002) and exp(β_i) can be interpreted directly as the odds ratios. We also used standardised coefficients (averaged across the same set of models) to assess the relative effects of different covariates measured at different scales (Schielzeth 2010). To do this we centred and scaled all variables in each season (mean = 0, s.d. = 1;McAlpine et al. 2008).

The model weights in the 95% set were recalculated to sum to 1 and each model within the set was used to create a model prediction surface (Anderson 2008). Each prediction surface within the model set was multiplied by its weight and then added together to produce a final model averaged prediction surface (Aldridge et al. 2012). This process allows all plausible models in a set to be used in multimodel inferences for spatial predictions (Anderson 2008).

We used the final model averaged prediction surface for each season to create a prediction surface in ArcMap 10.1 based on the associated habitat variables (ArcGIS 10.1; Environmental Research Systems Institute, Redlands, CA). We applied the logistic equation in the following form to create the relative probability of GUSG presence across the Gunnison Basin:

where x_n are covariates at location i for bird j with fixed effect regression coefficients β_n, and β₀ is the mean intercept. The logistic function was used to create a relative probability of presence surface with values between 0 and 1 across the GB for each season (1 = high, 0 = low).

In order to compare our models with the critical habitat, we completed an error matrix for the breeding model, the summer model and the combination of both seasonal models to critical habitat. We used prediction values of >50% relative probability of presence to define occupied habitat and all relative probability of presence values of <49% as unoccupied habitat. We used overall accuracy and the true skill statistic (TSS) as measures of comparison accuracy (Allouche et al. 2006). The TSS ranges from –1 to 1, where 1 equals perfect agreement and values of 0 or less indicate a performance no better than random (Allouche et al. 2006).

The gold standard for any model is to test model predictions against independently collected data (Wiens et al. 2008). We used two independent telemetry datasets for validation of the breeding and summer seasonal models as well as lek location data in the breeding season. We captured and marked GUSG for both independent datasets using the same techniques described earlier. The first dataset was collected in 2002 (validation set A) and the second dataset was collected during 2010–2011 (validation set B). In both studies, we circled each bird at a <30 m radius and manually corrected the Universal Transverse Mercator (UTM) location to reduce observation error. In addition, because it was important to test predictions in areas where data had not been collected for the original model development, and to verify that the model results could extrapolate across the GB, validation set B was collected in the western portion of GB where our seasonal models were extrapolated (Fig. 1). We separated each validation dataset into seasons based on the original models (breeding: 1 April–July 15; summer: 16 July–30 September). For the validation datasets, we calculated the proportion of locations or leks for which the model predicted >50% relative probability of presence for each seasonal prediction model (Sawyer et al. 2007).

We collected 7643 locations between April 2004 and September 2010, from 210 radio-marked GUSG. There were 188 GUSG and 3922 locations for the breeding season and 171 GUSG and 3721 locations for the summer season. Movement buffers for GUSG in the GB tended to be smaller during the breeding season (180.5 m) than during the summer (223.0 m). GUSG locations in both the breeding and summer seasons occurred in areas with higher proportions of sagebrush, higher water density, proximity to sagebrush, lower road density, lower residential density and further from residential than the available locations (Table 1).

Our breeding model suggests that GUSG locations were 21 times more likely to be within sagebrush and 16 times more likely to be within grassland than available locations (Table 2, Fig. 2a). Additionally, GUSG locations were five times more likely to be further from wetlands than available locations, but two times as likely to be in areas of high water density. We also detected smaller effects where GUSG were closer to easements, further from residential development and highways, occupying slightly lower elevations and within lower residential densities (Table 2).

Table 2.  Top model standardised and unstandardised coefficients for Gunnison sage-grouse resource selection in Gunnison Basin, Colorado, USA (2004–10) in the breeding season (1 April–15 July) including the standard error (s.e.), lower (LCI) and upper (UCI) 95% confidence intervals and odds ratios

Fig. 2.  (a) Relative probability of Gunnison sage-grouse use in Gunnison Basin, Colorado, USA (2004–12) during the breeding season (1 April–15 July). (b) Relative probability of Gunnison sage-grouse use in the Gunnison Basin, Colorado, USA (2004–12) during the summer season (16 July–30 September).

Our summer model indicated that GUSG locations were 20 times as likely to be in sagebrush, 57 times as likely to be in irrigated agriculture, and 12 times as likely to be in bare cover types compared with available locations, indicating a diverse use of habitat types (Table 3, Fig. 2b). Smaller effects included GUSG being further from wetlands and residential development, closer to highways and in areas with low road density and high water density during the summer season.

Table 3.  Top model standardised and unstandardised coefficients for Gunnison sage-grouse resource selection in Gunnison Basin, Colorado, USA (2004–10) in the summer season (16 July–30 September) including the standard error (s.e.), lower (LCI) and upper (UCI) 95% confidence intervals and odds ratios

Both of our models validated well with the independent datasets with slightly better results for the breeding model. Our validation set A had 507 locations (breeding, n = 262; summer, n = 245) and validation set B had 1600 locations (breeding, n = 762; summer, n = 838). The percentages of locations with a relative probability >0.50 for the breeding model were 83.0% for validation set A, 84.0% for validation set B and 85% for leks. The percentages of locations with a relative probability >0.50 for the summer model were 75.0% for validation set A and 79% for validation set B.

The current critical habitat designation without exclusions was approximately the same as the rangewide map provided by Colorado Parks and Wildlife, but was reduced by ~15% with exclusions removed (Table 4). For comparison, we estimated the area of occupied habitat for the breeding and summer models by summing the predictive surface with values >0.50 relative probability of presence threshold, assuming that values above 0.50 represent the habitat most likely to support GUSG (Table 4). Although the critical habitat is not divided by seasons, the area encompassed by the breeding and summer models was 27.9% and 17.9% smaller, respectively, compared with critical habitat (Table 4). Because we created two models based on habitat use during the breeding and summer seasons, whereas critical habitat represents overall habitat (regardless of seasonal use), direct comparison between the current critical habitat and our models is difficult. However, we intersected the breeding and summer seasonal models and used this as a conservative way to look at overall potential differences compared with critical habitat. By intersecting the breeding and summer models, the area was 4.8% smaller compared with critical habitat (Table 4), however, this approach devalues the benefit and importance of using the seasonally specific models, given the pattern and configuration of the occupied and unoccupied habitat is quite different to the critical habitat. Overall, the occupied and unoccupied habitat in the critical habitat map agreed with the combined breeding and summer models 68.5% of the time (Supplementary Material 2). Therefore, 31.5% of the cells in the critical habitat designation did not align with our combined models (Supplementary Material 2). The breeding and summer models agreed with the critical habitat in 56.4% and 58.4% of the cells, respectively. We found that our combined model predicted unoccupied habitat in areas the critical habitat predicted occupied in 20.9% of the cells (Supplementary Material 2). In contrast, we found our combined model predicted occupied habitat where the critical habitat predicated unoccupied habitat in 10.6% of the cells. Therefore, our models identified areas in the unoccupied critical habitat designation that were predicted to be relatively better quality habitat and areas in the occupied critical habitat predicted to be relatively poorer in quality (Fig. 3). This illustrates that although the overall amount of habitat was similar between the critical habitat designation and our models, the pattern and configuration of that habitat were very different.

Table 4.  Approximate area in hectares estimated by the seasonal models using the 50% relatively probability cutoff compared with both the Colorado Parks and Wildlife (CPW) occupied range map for Gunnison sage-grouse and the critical habitat designated by USFWS in 2014 (raster values have been rounded so totals may not be exactly the same)

Fig. 3.  Occupied and potential critical habitat designated by US Fish and Wildlife Service in 2014 compared with the breeding (1 April–15 July), summer (16 July–30 September) and combined seasonal model predictions (occupied habitat is based on >0.50 relative probability cutoff value).

Our breeding season model confirms similar habitat patterns found in the only other landscape modelling effort for GUSG that predicted GUSG use in higher proportions of sagebrush, further from roads and lower residential densities (Aldridge et al. 2012). Our analysis provides the first predictive summer season model for GUSG. Our results show that GUSG are dependent on large patches of sagebrush during both the breeding and summer seasons, although they tend to have more diverse habitat selection during the summer when they utilise numerous habitats, including irrigated agriculture and bare ground. Use of smaller sagebrush openings in a diverse suite of cover types has been found within Colorado with GRSG (Hausleitner 2003; Thompson 2012). In addition, use of agriculture is not uncommon during the summer season and has been documented in the literature (Connelly et al. 2011) and by many residents in the GB (Knapp et al. 2013).

The USFWS (2014b) recognises that local scale habitat structure and quality vary spatially, with some areas providing habitat for one or more seasons, but to date, critical habitat has not been defined on a seasonal basis. Although the USFWS (2014b) identified fine scale habitat structure attributes, they acknowledged that the data to delineate seasonal habitats at a landscape level did not exist at the time of designation (USFWS 2014b). We suggest that our seasonal models provide an analysis of a long-term dataset that the USFWS can use in addition to the estimated critical habitat within the GB to address these information gaps while establishing recovery goals and siting conservation actions. Designation of critical habitat does not lend itself to seasonal models, because any given site is usually classified as either critical habitat or not. We are not proposing to change how critical habitat is designated, but rather to include more data-driven models when they are available as a second tier of designation. In addition, the spatial arrangement and relative probability values provide a more accurate, data-driven assessment of the environmental variables that influence GUSG habitat use within the GB. This information provides a much stronger ability to identify on the ground conservation measures for species conservation when the information is available and can assist in targeted recovery of the species.

There are multiple reasons why multi-scale spatially explicit predictive surfaces could help conservation-reliant species such as GUSG (Scott et al. 2010). First, the critical habitat designation does not explicitly address the environmental variables important to GUSG. The incorporation of data-driven relationships with environmental variables allows management agencies to direct their management efforts on these variables in a recovery effort. For example, the summer model indicates that multiple habitat types are important to GUSG; this might guide management actions. If restoration and rehabilitation of habitat is a key component to GUSG recovery (Wisdom et al. 2011), management needs to identify finer scale priorities for habitat management. Accomplishing meaningful habitat protection requires more than just a large spatial area, but rather broad-scale land use planning (Kalen 2014), which our models can support.

Second, our analysis provides the most current and best science available on GUSG resource selection in the GB. We agree with the recommendation by Murphy and Weiland (2016) to use contemporary data that reflect spatial and temporal patterns in resource use by the species as one way to meet the best available science directive in the ESA (section 4(b) (2)) as well as ensuring data is current, quantitative and documented (Clark et al. 2002). This analysis attempts to fulfil these recommendations because it directly incorporates telemetry data over a 6-year period throughout two seasons of the GUSG lifecycle.

Finally, there is value in evaluating critical habitat at multiple spatial scales (Shirk et al. 2014) because species may respond differently at larger or smaller scales (Cunningham et al. 2014). Consideration of scale is necessary for deciding how habitat data should be applied in resource management (Boyce et al. 2003). Thus, we suggest incorporating both the general critical habitat outlined by the USFWS and also our finer scale models to focus and prioritise habitat restoration and recovery efforts. We caution managers that drawing conclusions about habitat selection based on observations at any one scale may misconstrue the importance of variables driving system behaviour overall (Doherty et al. 2010). Using a hierarchical approach to conservation planning is not novel (Johnson et al. 2004), but we are unaware of any use of the application for critical habitat designation, nor the use in policy actions related to ESA in the USA.

Our models provide important insights into habitat use patterns of GUSG in the breeding and summer seasons, but information illustrating winter habitat use patterns is still lacking. We recommend allocating future resources and research on collecting data during the winter season. The GB population has been the focus of much research (Aldridge et al. 2012; Davis et al. 2014; Davis et al. 2015), so we caution applying our models to the other seven satellite populations (<10–15% rangewide population; Fig. 1), where data is more scarce and local seasonal use patterns may differ. More monitoring and research could be directed towards the satellite populations to improve these data gaps.

Gunnison Basin stakeholders are keenly aware of the seasonal habitat use patterns of GUSG (Knapp et al. 2013) and have the perception, real or assumed, about the restrictions that can be imposed by critical habitat designation on the public and private sectors (Salzman 1990). Because of these concerns, stakeholders assisted us in selecting variables for our model development; they are invested and engaged in the recovery of GUSG. Gunnison County has filed an intent to sue the federal government for numerous issues, such as failing to use the best available science and issues related to the critical habitat designation (Gunnison County 2014), including failure to demonstrate that critical habitat designation is essential to conserving the species, designating habitat that was not beneficial to the species and designating critical habitat based on flawed and unsubstantiated conjecture regarding historic range (Gunnison County 2014). Our analysis provides an opportunity for the USFWS to use a more detailed, data-driven analysis when identifying where to target conservation actions and inform recovery goals. We do not believe that our models can replace the current critical habitat designation because the intent of critical habitat incorporates economic valuation of the landscape and pre-existing conservation easements by law. However, by utilising these models in the recovery planning process, there’s an opportunity to do so with a secondary product that has stakeholder support, which could lead to more successful, collaborative, strategically sited conservation efforts. We recommend that the USFWS take an adaptive approach to conserving this population by using our models as new information combined with stakeholder input to guide conservation efforts going forward. Using a hierarchical approach to GUSG recovery – by including our models to strategically identify seasonal areas to concentrate habitat improvements, conservation and restoration within the GB – may assist in gaining local support in GUSG conservation and recovery to accelerate the delisting process. This hierarchical approach for GUSG could then be used by USFWS on other species that have more detailed data available to inform critical habitat designations.

The authors declare no conflicts of interest.