Black-tailed prairie dogs in Kansas

The historic extent of the BTPD in Kansas includes the portion of the state west of the Flint Hills and the tallgrass prairie ecosystem. Lantz (1903) estimated that the BTPD once occupied over two million acres in sixty-eight of the 105 Kansas counties. However, BTPD populations have declined dramatically over the past century. An aerial survey by the Kansas Department of Wildlife and Parks (KDWP) in 2001 estimated that the BTPD only occupied approximately 130,521 acres in the state (Pontius 2002). Based on a sampling of BTPD towns from the same aerial survey, Pontius (2002) estimated fewer than 5,000 total BTPD towns in western Kansas.

BTPD population declines have had adverse effects on other species, since the BTPD is a keystone species in the shortgrass prairie ecosystem in Kansas and the broader Great Plains region (Agnew, Uresk and Hansen 1986; Desmond, Savidge and Eskridge 2000; Kotliar et al. 1999; Miller, Ceballos and Reading 1994; Wuerthner 1997). Several endangered, threatened, or species of concern that may be found in parts of the Great Plains depend on the BTPD for survival. For example, BTPDs might provide critical food and habitat resources for the black-footed ferret (Mustela nigripes), ferruginous hawk (Buteo regalis), and swift fox (Vulpes velox). Monitoring and assessing BTPD populations may provide key information regarding the status of these other species.

Genetic Algorithm for Rule-set Prediction

Ecological niche modeling is based on the premise that a given species distribution is limited, in part, by the presence of suitable habitat on the landscape. This suitable habitat may be understood as the joint description of the environmental conditions that allow a species to satisfy its minimum requirements (Chase and Leibold 2003), and is defined by various ecological parameters. A benefit of such modeling is the ability to extend a limited number of known habitat locations for a species (or species occurrence data) to “predict” the geographic distribution of the species for a larger area. The Genetic Algorithm for Rule-Set Prediction (GARP) (Stockwell and Peters 1999) provides one such approach to ecological niche modeling. The GARP algorithm has been implemented in the Desktop GARP software of Scachetti-Pereira (2001). GARP is a genetic algorithm that is intelligent, iterative, and based on machine-learning. Ecological niche modeling with GARP requires the selection of two types of GIS data layers for input: 1) environmental variables that collectively characterize a species niche and 2) species occurrence data in the form of x and y geographic coordinates (latitude/longitude, UTM, etc.). A particular strength of GARP is the ability to generate accurate models with a limited number of species occurrence data points (Stockwell and Peterson 2002).

GARP describes relationships between species occurrence data and environmental variables through a number of rule types (atomic, range, negated range, and logistic regression [logit]), and, following internal testing and refinement, develops a set of rules to create a model (Stockwell and Peters 1999). GARP refines the rule-set until either the maximum number of iterations or model convergence is achieved, and outputs a user-defined number of models for each species. Each model is a raster map layer that predicts either absence or presence of the species for each pixel in the layer. An automatic “best subsets” option in GARP may be implemented to select the best model set from all models (see Anderson, Lew and Peterson (2003) for discussion of the best subsets selection procedures). Typically, a 10-best model set is selected from all GARP models (Anderson, Lew and Peterson 2003).

GARP has yielded successful predictive results in a variety of ecology, biodiversity, and related applications. For example, Peterson, Ball and Cohoon (2002) modeled ecological niches for twenty-five bird species in Mexico with GARP. Peterson et al. (2001) projected ecological niche models developed in GARP to predict the potential effects of global climate change on the distribution of the Cracidae bird species in Mexico. In addition to terrestrial species, GARP has been extended to aquatic environments for ecological modeling of fish species in both freshwater (McNyset 2005) and marine (Wiley et al. 2003) environments.

Environmental Variables

Generally, BTPDs prefer sites with flat or gentle slopes, deep soils with a low sand content, and vegetation characteristic of the shortgrass prairie. Environmental variables selected for ecological niche modeling were Normalized Difference Vegetation Index (NDVI) data layers for the spring, summer, and fall growing seasons (Fig. 1); slope (Fig. 2A); soil depth (Fig. 2B); and soil texture (Fig. 2C). Environmental variables were selected based on known characteristics of BTPD habitats, in addition to consultation with Kansas BTPD experts and a literature review of related studies for other regions of the Great Plains. Proctor, Beltz and Haskins (1998) used similar variables (vegetation, slope, soil depth, and soil texture) to model BTPD habitats in the northern Great Plains with GIS, and through a classification tree analysis and logistic regression found correlations between three of these variables (vegetation, slope, and soil texture) with BTPD habitat presence.

Figure 1.

Examples of Normalized Difference Vegetation Index (NDVI) environmental variable layers calculated from Moderate Resolution Imaging Spectroradiometer (MODIS) satellite imagery for (A) March 22, (B) July 12, and (C) September 30. The figure displays the NDVI layers for a portion of southwestern Kansas.

Figure 2.

Environmental variable layers used in the ecological niche modeling: (A) Slope, (B) Soil depth, and (C) Soil texture. Data were not available for Doniphan and Johnson Counties (both outside the historic range of the BTPD) at the time of analysis.

NDVI, a measure of the amount and condition of biomass and vegetation, was selected to serve as a surrogate for land cover in the ecological niche modeling. Although the 250 m spatial resolution of the Moderate Resolution Imaging Spectroradiometer (MODIS) satellite sensor is coarser than the 30 m resolution of the Landsat Thematic Mapper (TM) sensor, a benefit of MODIS is improved temporal resolution. Other land cover mapping studies (e.g., Wardlow 2006; Wessels et al. 2004) have noted the advantages of the MODIS temporal resolution for discriminating among land cover types throughout the growing season. NDVI data layers were derived from biweekly MODIS satellite imagery composites for the year 2001 from March through October (fifteen dates in total) to encompass typical vegetation of the growing season in Kansas. NDVI was calculated for each of the fifteen MODIS images using the following equation: [(Near-Infrared Band − Red Band) / (Near-Infrared Band + Red Band)].

Soil depth and soil texture layers were derived from the Soil Survey Geographic (SSURGO) database, a detailed digital soil survey (1:24,000 scale) compiled on a county basis by the Natural Resources Conservation Service of the U.S. Department of Agriculture. Attribute data from the Soils Lite database for Kansas, a subset of commonly used SSURGO soil attributes, were joined to the SSURGO map unit polygons. The Soils Lite database includes measures of the average soil depth (in) to bedrock and an ordinal classification (1–21) of soil surface texture ranging from clay (1) to coarse sand (21) for each SSURGO map unit polygon. The SSURGO and Soils Lite databases were available for all Kansas counties in the study area.

Slope (%) was calculated from 30 m digital elevation models (DEMs) available from the U.S. Geological Survey's (USGS) National Elevation Dataset (NED). Striping, rows of cells with identical elevation values, and other data artifacts were present in the NED dataset in portions of Dickinson, Morton, Saline, and Sedgwick Counties. Corrected 10 m DEMs, available for each of the USGS topographic quadrangles with striping problems, were resampled to 30 m and replaced the quadrangles with striping in the NED dataset.

Species Occurrence Data

Species occurrence point data were selected from a database of known BTPD colonies provided by the KDWP. Latitude and longitude coordinates for each colony were visually estimated by observers in airplanes flying along north-south oriented transects across western Kansas. Although an advantage of the data collection method is a large sample size with coverage for much of the western half of the state, one limitation of the aerial survey data is the positional accuracy of each estimated coordinate. For this reason, the Kansas GAP Analysis Project land cover map (Egbert et al. 2001) was used to eliminate any species occurrence points that corresponded to land cover classes unsuitable for BTPD habitation (cropland, Conservation Reserve Program (CRP), forest/woodlands, marshland, urban, and water) as a result of positional error introduced during the coordinate estimation procedures. The remaining species occurrence points (455 total points) were randomly subset into training (80%, 364 points) and validation (20%, 91 points) datasets (Fig. 3). The training dataset was input into GARP for training and internal testing during the modeling process, whereas the validation dataset was withheld entirely for an independent validation of model results through Receiver Operating Characteristic (ROC) analysis. An independent validation data subset is preferable to data re-substitution for assessing model quality (Fielding and Bell 1997). Given the number of occurrence points available, an 80/20 training/testing subset was used as this satisfies the minimum required number of training points necessary to achieve maximal GARP model accuracy (Stockwell and Peterson 2002), while maintaining a testing subset of sufficient size to achieve statistical power in the ROC analysis (Hanley and McNeil 1982).

Figure 3.

(A) Training and (B) validation BTPD species occurrence data used in the ecological niche modeling.

Methods

Desktop GARP requires that all environmental layers be in raster data format at the same spatial resolution or pixel size. Therefore, the slope layer was resampled from 30 m to 250 m and the soil texture and soil depth to bedrock layers were converted from vector to raster format and sampled to 250 m to match the native spatial resolution of the NDVI layers derived from the MODIS satellite imagery.

A mask layer was created to exclude all geographic areas in Kansas east of the Flint Hills (i.e., beyond the historic range of the BTPD) in the resulting raster layers from the GARP modeling. In addition, cropland, CRP, forest/woodlands, marshland, urban, and water land cover classes were extracted from the Kansas GAP land cover map and added to the mask layer in order to exclude all land cover types unsuitable for BTPD habitation. CRP was included in the mask since CRP land is generally considered unsuitable for BTPD habitation (USDA Forest Service 2005).

In order to achieve more regional variation in the resulting ecological niche models, the environmental variables and species occurrence datasets for the state were divided into two zones, one north and one south of 38° 42′N latitude. Preliminary results by Stockwell and Peterson (2002) have demonstrated that such regional subdivision may improve GARP model results. A north-south division of the state was preferable to an east-west split in order to ensure an equivalent number of species occurrence data points in each zone.

Jack-knifing was performed in Desktop GARP separately for the north and south zones to determine if any environmental layers should be excluded from the ecological niche modeling due to high correlations with omission, i.e., an error corresponding to an absence prediction where species occurrence data actually occur, or a “false negative”. Environmental layers with a correlation (r) to omission of greater than 0.10 were removed from the ecological niche modeling. One NDVI layer (October 16) for the south zone (r = 0.11), and one NDVI layer (September 14) for the north zone (r = 0.13) were removed due to correlations with omission that exceeded 0.10.

GARP was run separately for the north and south zones with a convergence level of 0.001 and 10000 maximum iterations for both runs. The “best subsets” option in Desktop GARP was enabled (omission threshold = 5%, commission threshold = 50% of distribution, 20 total models under hard omission threshold) to select automatically the 10-best model set for each of the two zones. The 10-best model sets for both zones were merged together to create the final BTPD habitat suitability map.

Model Validation with Receiver Operating Characteristics (ROC) Analysis and Average Omission Error

Two methods of model validation commonly used in ecological niche modeling, Receiver Operating Characteristics (ROC) analysis (e.g., McNyset 2005; Wiley et al. 2003) and average omission error calculation (e.g., McNyset 2005), were used to assess the overall accuracy of the 10-best model set. ROC analysis compares model results with the species occurrence validation data withheld from the GARP modeling. The area under the curve (AUC) statistic is a measure of the overall accuracy of the model set; an AUC statistic of 1.0 indicates total agreement between the model set and the species occurrence test data whereas a model set developed randomly yields an AUC of 0.5. Average omission error is a measure of the average omission for the 10-best models. The BTPD model set for Kansas yielded an AUC = 0.701 and an average omission = 7%, indicating good overall model set predictions.

Model Validation with an Additional Independent Dataset: An additional dataset of thirty-seven active BTPD colonies in the Cimarron National Grassland in Morton County was available for an independent validation of the 10-best model set (Fig. 5). These colonies were on publicly held portions of the Cimarron National Grassland, and were identified by visually and audibly locating BTPDs in the area. Perimeters of these BTPD colonies were mapped as polygon features with a differential Global Positioning System (GPS) receiver and then post-processed with data from a base station in Elkhart, KS. The use of differential GPS for mapping these colonies resulted in locations that were considerably more accurate than the statewide BTPD species occurrence dataset.

Figure 5.

BTPD colony perimeters delimited with GPS overlaid on GARP 10-best model set for Morton County, Kansas.

The centroid of each BTPD colony polygon was computed and used to calculate the AUC and average omission for Morton County (AUC = 0.628, average omission = 6%). In addition, the 10-best model set values for each pixel in each of the polygons were summed and converted to percentages of total area for each BTPD colony (Table 1). The results indicate a high degree of accuracy in the model predictions. For example, each of the 10-best models predicted presence (i.e., values of 10) for at least 75% of the total area of the colony for thirty of the thirty-seven (81.1%) colonies. Likewise, each of the 10-best models predicted presence for 100% of the total colony area for seventeen of the thirty-seven (45.9%) colonies.

Table 1.

Calculated area (in square meters) of masked and 10-best model values for each BTPD colony in Morton County, Kansas as well as percentage of total colony area.

Interestingly, seventeen of the thirty-seven (45.9%) colonies fell partially within the masked land use/land cover types with percentages ranging from 2.4% to 99.7% of the total colony area. To evaluate this artifact, colonies were overlaid on top of the GAP land cover map used to generate the mask. Perhaps since the landscape is highly fragmented in western Kansas, several BTPD colony boundaries slightly extend beyond prairie land cover types into cropland or other land use/land cover types that are deemed unsuitable for BTPD habitat. These particular colonies contain relatively small percentages of the masked area.

Colonies with higher percentages of the masked area typically fell on land classified as CRP by the GAP land cover map. While land enrolled in CRP is planted to grassland species, it is generally considered to be unsuitable for BTPD habitat (USDA Forest Service 2005). However, it is possible that BTPDs have successfully established colonies on CRP in western Kansas and would suggest that CRP should not have been masked a priori to the GARP modeling. Perhaps a more likely explanation is that native prairie in these areas were misclassified as CRP in the GAP land cover map used for the mask. While general land cover classes (herbaceous, scrubland, forest, cropland, urban, and water) were mapped at an overall accuracy of 88%, accuracy levels for specific grassland alliance-level communities were substantially lower (Egbert et al. 2001). Field verification of the GAP land cover types that correspond to BTPD colonies in Morton County is necessary to verify these possible explanations.