Potential landscapes for conservation of the black-tailed prairie dog ecosystem

Aim: To identify potential landscapes for the conservation of the black-tailed prairie dog (BTPD) ecosystem, across their historical geographic range within the United States.

Location: Central Grasslands of the United States.

Methods: We used a structured decision analysis approach to identify landscapes with high conservation potential (HCP) for the BTPD ecosystem. Our analysis incorporated ecological, political, and social factors, along with changing climate and land use to maximize long-term conservation potential. We created scenarios that involved current and future projected suitable BTPD habitat, across the BTPD range within the United States. These were our RANGEWIDE scenarios. Additionally, because conservation policies and funding decisions are often made by political entities, we also identified STATE-LEVEL conservation priorities, under both present and projected future climate. Our STATE-LEVEL analysis sought conservation solutions within each of the states’ boundaries only, so do not consider a rangewide perspective.

Results: The landscapes we identified with HCP (top 30% range-wide) represented 22% of the historical distribution of black-tailed prairie dogs and remained strongholds under projected climate change. We provide a suite of HCP area scenarios to help inform different conservation and management interests, including those that consider projected climate change and jurisdictional (state-level) boundaries. STATE-LEVEL conservation priorities differed considerably from RANGEWIDE priorities, under both current and future climate scenarios. The largest difference was among the southern states (Arizona, New Mexico, and Texas), where climate change reduces the conservation priorities across this region more when viewed from a RANGEWIDE perspective than when viewed from a STATE-LEVEL perspective. Additionally, from a RANGEWIDE perspective, the eastern states have fewer areas with HCP compared to the western states within the BTPD range, but when viewed from a STATE-LEVEL perspective there are considerably more areas with HCP. We expected such differences because this question was aimed at understanding the HCP areas within each state, so the analysis was seeking conservation solutions within each of the states’ boundaries. Identifying STATE-LEVEL conservation priorities is important because funding sources and management priorities are often focused at the state-level, and not range-wide. This way, each state has information on conservation priorities within their own jurisdictional boundaries. We suggest each state focus conservation efforts for the BTPD ecosystem in those areas that remain priorities into the future at the STATE-LEVEL, while also considering those priorities identified within their state under the RANGEWIDE perspective.

Main Conclusions: Our findings highlight the large conservation potential for BTPDs and associated species, and the maps we generated can be incorporated into other large-scale, multi-species conservation planning efforts being developed for the Central Grasslands of North America.

Methods
Description of the data and file structure
Spatial Data Layers Used in Conservation Prioritization Analysis

We used the spatial conservation prioritization method and Zonation software (Moilanen et al. 2005) to evaluate how landscapes varied in their potential for prairie dog ecosystem conservation and restoration across the full range of the species in the United States. Our analysis included a total of 23 environmental input datasets for the full study area, based on the data sources described in Table 1. The most important layer we used to inform our analysis was the BTPD habitat suitability model, as it provided the basis for where, ecologically, the best places are to conserve and restore the BTPD ecosystem (Davidson et al. 2023). This habitat suitability model (HSM) was based on presence and absence data for BTPD occurrences across their geographic range within the United States (McDonald et al. 2015), and quantified how prairie dog occurrences related to climate, soils, topography, and land cover (see Davidson et al. 2023 for details). We also utilized HSMs for BTPDs under two future climate scenarios: 1) warm and wet and 2) hot and dry, to inform where the most ecologically suitable habitat will likely be located under a warming climate (Davidson et al. 2023).  

However, the goal of our analysis was to not only determine potential landscapes for conservation based on local habitat suitability, but also to examine how the distribution and connectivity of native grassland habitat at broad spatial scales, the distribution of threats to prairie dog habitat (such as development and conversion to cropland), and the political and social landscape collectively influence opportunities to conserve and restore the BTPD ecosystem (Table 1; Fig. S1). We used the 2016 National Land Cover Database (NLCD) to inform on the location, extent, and connectivity of favorable habitat (grassland/shrubland), versus unfavorable habitat (forest/woodland and emergent wetland) for prairie dogs (USGS 2019a). We also created a landscape fragmentation layer by mapping the degree of rangeland fragmentation across the historical BTPD range. To do this, we followed the methods of Augustine et al. (2021), except that we used the 2016 NLCD as the source data layer rather than a combination of the 2011 NLCD and USDA Cropland Data Layers. Briefly, every pixel was classified as either (1) rangeland, which we defined as grassland, shrubland, and improved pasture/hay cover types, (2) a fragmenting land cover type, which we defined as forest, cropland, or developed lands, or (3) neutral land cover types which were not rangeland, but also did not fragment adjacent rangelands. In the final fragmentation map, we set all pixels mapped as either a fragmenting or a neutral land cover type to a value of zero, and then calculated the distance to the nearest fragmenting land cover type for each rangeland pixel (e.g., Figure 3 of Augustine et al. 2021). Additionally, we incorporated spatial data on land use: oil and gas well locations, distance to transmission lines, wind turbine count, and road density (Homeland Security Infrastructure Program (HSIP) 2020; United States Census Bureau 2020; Federal Aviation Administration 2021; Welldatabase 2021). These land use data layers provide information on anthropogenic activity that reflect the presence of humans and habitat quality. Areas that have higher levels of human activity may be less favorable for the BTPD ecosystem because of the increased potential for shooting of prairie dogs, impacts on associated species through behavioral modification, and habitat degradation. We also included spatial layers on projected habitat loss. The tillage risk layer (Olimb & Robinson 2019) informs where habitat is most likely to be lost to cropland. Further, we included scenarios of overall landcover change projected into the future (Sohl et al. 2018), with a focus on areas that would retain the greatest amount of favorable grassland habitat. We then obtained PAD-US (USGS 2019b), National Conservation Easement Database (NCED; (Ducks Unlimited & The Trust for Public Land 2021)), and other private conservation land data to determine the landownership of identified HCP areas (Table 1). We also obtained data from Carlson et al. 2022 (Wildlife Model in Figure 1) to relate HCP areas to plague risk.

We also included social and political spatial data in our analysis. We collated percent of Conservation Reserve Program (CRP) grasslands per county and the League of Conservation Voters Conservation Score Card (LCVCSC) to reflect political and social support for the environment (on a per county basis) (USDA Farm Service Agency 2020; League of Conservation Voters 2022). We also included data from a novel survey of wildlife governance preferences delivered to Canadian, Mexican, and American residents (Anonymized et al. 2023(a), 2023(b)) to determine the probability that a region would support increases in prairie dog populations or support federal or private incentives for prairie dog conservation. Census tract level estimates were generated using a Bayesian multi-level regression with post stratification wherein the demographics of survey respondents are used to map the probability to census geographies based on the demographic composition of the Census tracts (Anonymized et al. 2023(b); Gelman 2007; Hanretty 2020). Finally, we created a spatial layer of the count of Land and Water Conservation Fund (LWCF) projects  (The Wilderness Society 2015) to reflect a regions’ institutional capacity to actualize conservation.

Data Preparation

The data layers were integrated into a nested hexagon framework (NHF). A NHF grid is based around a 1 km2  hexagon unit that is aggregated up by units of 7 to generate coarser scale cells of 7 km2 (cogs), 49 km2 (wheels), and 343 km2 (rings), allowing for cross-scale multidisciplinary analysis while obscuring precise sensitive location data.  

A total of 31 data layers representing point, polygon, and raster formats were processed and summarized into the NHF for consideration in the Zonation analysis (Table S1). While the exact process used to integrate the data layers into the NHF and subsequently into raster files for the Zonation analysis was slightly different for each data layer, the general process was the same. All GIS data processing was done using ESRI ArcMap 10.7 software. Input data layers were intersected with the NHF and the data layers were summarized per NHF hexagon cell using Zonal Statistics, Tabulate Area, or other similar geoprocessing tools to generate a summary of the source layer data per hexagon. Examples of the resulting tabular summaries conveyed the area of each landcover class per hexagon cell (later converted to a percent), the mean tillage risk, majority landscape condition, the sum of the meters of road or number of wells within a cell, or the presence of wind turbines within each 1 km2 hexagon cell.

Within the attribute table of the hexagon feature class, a series of new attribute fields were created to convey the newly summarized data (e.g., % grassland, number of wells). Using the unique hexagon ID’s, the data tables of the summarized information were joined with the feature class attribute table, and the summarized data was copied into the newly created hexagon attribute fields using the “calculate field” process. Due to the number of hexagons (over 2 million record rows) being calculated, this process often took several days so researchers later began using a python script to “update cursor” that proved much faster than join/calculate field process. The resulting attribute table of the NHF one-kilometer cells provided a summary of the datasets integrated, all pre-summarized to the same framework for compatibility and easy use (Table S1). Some source data layers like percent of CRP and the political voting data were originally in coarse (county/voting district) spatial resolutions. As a result of summarizing these datasets to the hexagons, the results displayed a false level of spatial precision regarding the data values conveyed. In cases where coarse data was summarized and displayed at a higher spatial resolution, many individual hexagons share the same value that originally represented the district/county as a whole, not a specific hexagon. 

The hexagon feature class data was exported to a series of raster layers using the ArcMap Feature to Raster function to accommodate the conservation prioritization software requirements that all input data be in a raster format. Output raster layers were specified to have a 90 m resolution, were snapped to the same 90 m pixels as the ensemble habitat suitability models, and the raster values were derived from the values in each of the feature class attribute fields representing the 1 km2 hexagon summarized data. The intersect, calculate field, and convert to raster processes were done in batches using the 5x5 degree NHF tile or by regional groupings of 7 tiles for the northern half of the range and 9 tiles for the southern half of the range for efficient processing. After each tile was converted to a raster layer, they were mosaiced together to create a series of range-wide raster layers, and then clipped to the BTPD range boundary (Fig. S2).   

Prioritization Analysis

 We used Zonation, an approach and software for spatial conservation prioritization, to select HCP areas for the conservation of the prairie dog ecosystem. Zonation produces a hierarchical spatial priority ranking of the study region, accounting for complementarity by considering local representation of the biodiversity features (species, ecosystem types, etc.; Moilanen et al. 2005). Zonation iteratively removes cells whose removal causes the smallest loss in feature representation across the overall remaining region until no cells are left in the region. The hierarchical conservation rank of the region is based on the order of cell removal, which is recorded and can be used later to select any given top fraction (e.g., best 25%) of the region. We used the additive benefit function (ABF) removal rule, which is based on the sum of the features representation in each cell, favoring places containing high habitat quality for a large number of biodiversity features.

The relative weighting of data layers is an important component of the Zonation algorithm and impacts the order in which cells are removed from the prioritization landscape. Cells that contain a high-weight feature are kept longer in the analysis than cells with only low-weight features. Features with a negative weight are considered undesirable. Consequently, they are found among the cells with low conservation priority and removed from the landscape early in the analysis. To identify those areas with the highest potential for prairie dog ecosystem conservation, we used a weight of 10 for spatial layers describing habitat suitability for BTPDs, a weight of 1 for landscape-scale land use/land cover features that have a positive influence on conservation potential and a weight of 1 for social environment layers with a positive influence on conservation potential. The spatial layers were considered as features in the analysis with positive values (i.e., higher values indicated favorable places for BTPD conservation). Because suitable habitat is ultimately the most important variable for conservation, we assigned habitat suitability features with the highest weighting among all positive features. We also considered land use in the selection of priorities, aiming to avoid places with high intensity of anthropogenic activities and potential conservation conflicts. Those layers within the landscape-scale land use/land cover and risk of habitat loss categories that negatively affect conservation potential were given negative weights (-4). These areas consequently had low values of conservation priority and were removed from the study region early in the analysis. Details on each feature used can be found in Table 1. Areas with low habitat suitability or high sandy soil (>90%) were masked out of the analysis using an area mask file, where cells with value “1” were included in the analysis, while cells with value “0” were excluded (Table 1).

We used Zonation to evaluate conservation potential under various scenarios. First, we evaluated HCP areas across the geographic range of BTPDs using suitable habitat under the current climate. Next, we created scenarios that involved current and future projected suitable BTPD habitat, across the BTPD range within the United States. To do this, we used the interaction function that induces connectivity of suitable sites for the interacting features to account for distribution shifts due to climate change. Additionally, because conservation policies and funding decisions are often made by political entities, we also identified conservation priorities within each state, under both present and projected future climate. For this, we used the Administrative Units (ADMU) function in Zonation to also select state priorities in the final conservation ranking (Moilanen & Arponen 2011).