The U.S. Fish and Wildlife Service (USFWS) is responsible for reviewing the biological status of hundreds of species to determine federal status designations under the Endangered Species Act. The longleaf pine Pinus palustris ecological system supports many priority at-risk species designated for review, including five species of herpetofauna: gopher tortoise Gopherus polyphemus, southern hognose snake Heterodon simus, Florida pine snake Pituophis melanoleucus mugitus, gopher frog Lithobates (Rana) capito, and striped newt Notophthalmus perstriatus. To inform status decisions and conservation planning, we developed habitat suitability models to 1) identify habitat features that best predict species presence and 2) estimate the amount and distribution of suitable habitat across each species' range under current conditions. We incorporated expert judgment from federal, state, and other partners to capture variation in ecological settings across species' ranges, prioritize predictor variables to test in models, mitigate data limitations by informing the selection of pseudoabsence points, qualitatively evaluate model estimates, and improve the likelihood that experts will trust and use model predictions for conservation. Soil characteristics, land cover, and fire interval strongly influenced habitat suitability for all species. Suitable habitat was distributed on known species strongholds, as well as private lands without known species records. Between 4.7% (gopher frog) and 14.6% (gopher tortoise) of the area in a species' range was classified as suitable habitat, and between 28.1% (southern hognose snake) and 47.5% (gopher frog) of suitable habitat was located in patches larger than 1 km2 (100 ha) on publicly owned lands. By overlaying predictions for each species, we identified areas of suitable habitat for multiple species on protected and unprotected lands. These results have direct applications to management and conservation planning: partners can tailor site-level management based on attributes associated with high habitat suitability for species of concern; allocate survey effort in areas with suitable habitat but no known species records; and identify priority areas for management, land acquisitions, or other strategies based on the distribution of species records, suitable habitat, and land protection status. These results can aid regional partners in implementing effective conservation strategies and inform status designation decisions of the USFWS.
Currently, the U.S. Fish and Wildlife Service (USFWS) is responsible for reviewing the biological status of hundreds of at-risk species to determine if listing is warranted (hereafter, listing decisions) under the U.S. Endangered Species Act (ESA 1973, as amended; USFWS 2020a). Each of these at-risk species require formal status assessments to characterize risk of extinction and inform listing decisions (USFWS 2016; Smith et al. 2018), and many species are also the focus of conservation planning efforts by federal, state, and other partners to recover or sustain species before they become rarer (Sutherland and deMaynadier 2012; Pickens et al. 2017). Estimates of demographic rates, population abundance, and risk of extinction are valuable criteria that inform listing and conservation decisions (Beissinger and Westphal 1998; McGowan et al. 2017); however, acquiring these estimates remains challenging for most species including rare or cryptic species on which conservation efforts often focus. In light of these challenges, status assessments and conservation planning efforts frequently rely on information about the spatial distribution of species and their habitats, which can be more easily estimated for data-limited species (e.g., Engler et al. 2004) and linked to projections of abundance or extinction risk for species if or when additional data are available (e.g., Larson et al. 2004; McGowan et al. 2017).
Habitat suitability models (HSMs: also called species distribution models) are now commonly used tools for estimating the distribution of species, their habitats, and threats (Guisan and Zimmermann 2000; Franklin 2010). In short, HSMs use measures of environmental and landscape attributes (e.g., soil characteristics, canopy cover, fragmentation, rainfall) in places where a species was known to occur over some time scale to project where similar conditions occur throughout the species' range (Boyce et al. 2002; Engler et al. 2004; Elith and Leathwick 2009; Franklin 2010). Known species locations (presence data) can be collected from records in natural history collections, via systematic surveys (e.g., Enge et al. 2014), or opportunistically (e.g., through volunteer-reported, or “citizen science”, data: Bradter et al. 2018). The projected habitat distributions can then be used to understand species–habitat relationships, predict where potentially suitable habitat is likely to occur, and prioritize areas for surveys for new populations, habitat management, parcel acquisition or designation, species translocation, and other applications related to conservation planning (Larson et al. 2004; Elith and Leathwick 2009; Franklin 2010; Guisan et al. 2013; Villero et al. 2017). Robust HSMs can also be used to project future habitat distributions under different environmental change scenarios (Thuiller 2003; Guisan et al. 2013) and inform long-term planning.
Developing HSMs that inform listing and conservation decisions at large spatial scales (i.e., a species' range) requires addressing at least two issues. First, ecological and genetic variability can exist in wide-ranging species, so researchers can improve predictions of HSMs developed across a species' range if models adequately account for this regional variation (Murphy and Lovett-Doust 2007). If researchers know or expect ecological or genetic variability, they can identify geographic units that capture this variability, and build HSMs to estimate different species–habitat relationships among these units (e.g., May et al. 2011). Second, spatially extensive information about species presences is often available from state-maintained and other databases (e.g., natural heritage programs: Groves et al. 1995), but absence information is not usually collected and available across a species' range. Multiple approaches have been developed to predict habitat suitability using presence and pseudoabsence (“background”) data (e.g., generalized linear models or machine-learning models, such as MaxEnt: Guisan and Zimmermann 2000; Phillips et al. 2006, 2009). However, the accuracy of HSM predictions generally increases with the quality of absence information provided (Brotons et al. 2004) and reduces the effects of spatial or detection biases (Gu and Swihart 2004). For example, Bradter et al. (2018) used records when experts reported nontarget, but not target, species as “inferred absences,” and found that supplying HSMs with inferred absences improved model accuracy compared to using random background points.
When no absence information is available, studies have used expert judgement to inform different model components, such as the set of habitat predictors evaluated in HSMs and prior estimates for species–habitat relationships, to supplement presence-only data and improve model predictions (Murray et al. 2009; O'Leary et al. 2009; Arfan et al. 2018; Reside et al. 2019). We note that some applications using expert judgment have found its inclusion did not substantially improve HSM predictive ability (e.g., Pearce et al. 2001; Seoane et al. 2005; Charney 2012); however, it is reasonable that adopting recommended best practices for expert elicitation can increase the accuracy of expert judgment and, in turn, the predictions of HSMs that it informs (Martin et al. 2012; Addison et al. 2013; Arfan et al. 2018; Reside et al. 2019). Moreover, involving the same experts during the model development phase who will ultimately be the end-users of model outputs can ensure those outputs are trusted by experts and likely to be relevant and used in decision making (Addison et al. 2013; Guisan et al. 2013).
We integrated expert opinion into range-wide HSMs for at-risk species designated for status assessments by the USFWS. Our focal species, five sympatric herpetofauna associated with the longleaf pine Pinus palustris ecosystem, were a priority for study by the USFWS, our principal partner. Our goals were to 1) identify habitat features that best predict species presence and 2) estimate the amount and distribution of potential suitable habitat across each species' range under current conditions at the time of the analysis (2018). We incorporated expert judgment from federal, state, and other partners with experience in research and management of these species to capture variation in ecological settings across species' ranges, prioritize predictors to test in models, mitigate data limitations by informing the selection of pseudoabsence points, and qualitatively evaluate model estimates. We also intended the collaborative involvement of experts in the modeling process to improve the likelihood that experts will trust and use model predictions for conservation planning (Villero et al. 2017). The results of this study will inform upcoming status assessments and spatially explicit conservation planning for these five focal species, and our broader approach—using expert judgment to improve the accuracy of HSMs—is applicable for other data-limited species.
Study system, species, and extent
The longleaf pine ecosystem supports many priority at-risk species designated for status reviews, including our five focal species of herpetofauna: the gopher tortoise Gopherus polyphemus, southern hognose snake Heterodon simus, Florida pine snake Pituophis melanoleucus mugitus, gopher frog Lithobates (Rana) capito, and striped newt Notophthalmus perstriatus. Collectively, these species span several ecoregions across the southeastern coastal plain of the United States from Louisiana to North Carolina. They are associated with uplands characterized as xeric longleaf pine savannahs, scrub, flatwoods, hammocks, coastal dunes, and sandhill habitats that typically have well-drained, sandy soil, low canopy cover, and adequate herbaceous ground cover (Auffenberg and Franz 1982; Tuberville et al. 2000; Smith et al. 2006; Roznik et al. 2009; Miller et al. 2012; Beane et al. 2014; Farmer et al. 2017). Frequent fire has long been established as the key driver of the natural disturbance regimes that influence population and community dynamics in longleaf pine and related systems (Platt et al. 1988; Glitzenstein et al. 1995; Van Lear et al. 2005). The gopher frog and striped newt additionally use isolated, ephemeral wetlands with open canopies (Greenberg 2001; Humphries and Sisson 2012; Enge et al. 2014; Farmer et al. 2017). Structural and functional features of these systems have been affected to varying degree by extensive habitat loss, fragmentation, and degradation from land use change and fire exclusion (Outcalt and Sheffield 1996). Consequently, regional or range-wide population declines have been observed or suspected over the past century due to these and other stressors for the gopher tortoise (Auffenberg and Franz 1982; Diemer 1986; Hermann et al. 2002; McCoy et al. 2006; Smith et al. 2006), southern hognose snake (Tuberville et al. 2000), Florida pine snake (Franz 1992; Miller et al. 2012), gopher frog (Semlitsch et al. 1995; Enge et al. 2014), and striped newt (Farmer et al. 2017). Although previous efforts have assessed habitat suitability at local scales (e.g., a single protected area, multiple counties) for the gopher tortoise (Baskaran et al. 2006; Kowal et al. 2014; Johnson et al. 2017) and at the range-wide scale for the gopher frog (Johnson et al. 2017) and striped newt (May et al. 2011), range-wide HSMs are needed for all five species that incorporate recent occurrence data and account for potential regional variation in species–habitat relationships.
Our study extent encompassed the combined ranges of the five focal species (range maps downloaded from the U.S. Geological Survey (USGS) Gap Analysis Project; USGS 2001) in the southeastern United States (Figure 1). The extent included four U.S. Environmental Protection Agency (EPA) Level III ecoregions—the Middle Atlantic Coastal Plain, Southeastern Plains, Southern Coastal Plains, and Southern Florida Coastal Plain—which were composed of several Level IV ecoregions representing variable ecological settings and habitat types, including the Sand Hills, Southern Pine Plains and Hills, Central Florida Ridges and Uplands, Gulf Coast Flatwoods, and Sea Island Flatwoods (EPA 2018b).
Species presence and pseudoabsence data
Unless otherwise noted, we performed all spatial analyses in ArcGIS version 10.4 (ESRI, Redlands, CA) and statistical analyses in R version 3.1 (R Core Team 2016). We compiled a geospatial database of occurrence records for each species (Table 1; Figure 1) from datasets maintained by natural heritage programs and state agencies within all states in our study extent, USFWS, U.S. Forest Service, U.S. Department of Defense, academic researchers, and HerpMapper—an online platform where species records are reported by the public and validated by professional herpetologists (HerpMapper 2018). Records included opportunistic sightings, as well as observations from systematic research and monitoring studies (e.g., Enge et al. 2014).
We applied two data filtering steps to maximize the likelihood that all presence points used in models indicate places with currently suitable habitat conditions that facilitate species occupancy. First, from the full location dataset, we removed any record found before 1981. We chose this cutoff year in consultation with experts because it achieved a balance of including the majority of records for our focal species, most of which were data-limited, while excluding earlier records that typically had lower spatial accuracy and represented locations where the habitat conditions have likely changed since the species was observed. Second, we removed any record found after 1980 that data managers (e.g., state agencies, natural heritage programs) ranked as likely extirpated; this ranking typically signified the location had been developed since the animal was observed. We checked the accuracy of these rankings by overlaying any point deemed likely extirpated with the 2011 National Land Cover Dataset (see below) to confirm the presence of developed or agriculture land. All records emerging from the filtering steps yielded location (decimal degrees) recorded to five or more decimal places (< 2 m). To reduce potential biases that can result from spatially clustered localities (Phillips et al. 2009; Boria et al. 2014), we applied a final filter by randomly removing records for each species so that no records in the final presence dataset occurred within 150 m of each other. We chose this distance by generating correlograms of predictors (see below) at presence points. Correlograms for predictors indicated the minimum distance where spatial autocorrelation was negligible (i.e., where 95% confidence intervals for Moran's I overlapped 0) ranged between 50 and 250 m, so we chose an interpoint distance of 150 m as an average value to indicate spatial independence.
Because robust absence data do not currently exist for these species, we used an HSM approach that uses presence and pseudoabsence points across a species' range to compare the range and variation of available habitat to the subset of habitat conditions found at known species locations (e.g., Engler et al. 2004; Barrett et al. 2014). However, we performed additional steps to use expert judgment to inform where and how pseudoabsence points were generated, which were designed to improve the quality of absence information and overall accuracy of HSM predictions. We conducted an expert elicitation exercise using Google Earth (Google Inc., Googleplex, Mountain View, CA) where species experts identified (outlined) areas in their region of work where they believed each species was 1) present, if they had no presence data but had heard unverified reports or had other reason to think the species was there, or 2) absent, if they had been to a site but had never detected the species or had other reason to think the species was absent. Experts had a minimum of 5 (and usually more than 15) years of experience with research and monitoring on our focal species in their assigned region.
For each species, we generated background points from a boundary that included all areas 1) within the boundary extent of the species' range, 2) outside of a 500-m buffer surrounding presence points, and 3) not within areas where experts believed species were likely present. We generated sampling pools of 100,000 pseudoabsence points for the gopher tortoise and 30,000 points for all other species to adequately sample range-wide environmental conditions. Because expert-defined absence areas likely represented conditions of low habitat suitability, we gave these areas more weight when randomly generating pseudoabsence points so that the percentage of pseudoabsence within these areas was roughly double the percentage of a species' range identified as absence areas. For example, whereas experts identified 9.4% of the range of the gopher tortoise as “absence areas,” we drew 20% of all pseudoabsence points from these areas (Table 1). Additionally, we accounted for potential bias toward road sampling that can impact HSM predictions (Phillips et al. 2009). We used the TIGER US road dataset (U.S. Census Bureau 2016) to identify the percentage of presence points for each reptile species located on or within 30 m from roads. We then used the road layer as a mask to generate an approximately equal percentage of pseudoabsence points on roads outside of the expert-absence areas. For example, 58.7% of Florida pine snake qualifying presence records were associated with a road (Table 1); therefore, from all nonexpert-identified absence areas, we drew 50% of the pseudoabsences from ≤ 30 m of roads. For the two amphibian species, experts did not identify a sufficient number of suspected presence or absence areas in the Google Earth exercise (areas covered < 3% of species' ranges), so we generated pseudoabsence points randomly across each species' range. Additionally, we did not account for road sampling bias for the amphibians because > 99% of records were found ≥ 50 m from roads.
We gathered important environmental and landscape attributes (hereafter, predictors) that were related to species' natural histories based on previous studies and expert input. Specifically, we had 27 species experts complete an online survey where they listed environmental, spatial, and biophysical attributes they associated with ideal habitat and presence of each species at a site and categorized attributes by their degree of influence (see Text S1, Figure S1, Supplemental Material) and discussed responses at an in-person workshop. We used expert responses and previous studies to prioritize spatial datasets of important predictors that were available across the study extent. We obtained spatial data in the form of 30-m rasters for 17 predictors and grouped these into seven categories: 1) geographic ecoregion groups, 2) edaphic (soil) factors, 3) vegetation, 4) disturbance and connectivity, 5) wetlands, 6) climate, and 7) topography (Table 2). We “smoothed” several predictors (see Table 2) by calculating average conditions within neighborhoods from 90 to 2,000 m (i.e., a moving window) to which species presence may be influenced by conditions at various scales (e.g., within a home range), based on expert input and previous studies. We extracted values of all predictors to presence and pseudoabsence points for fitting HSMs. We briefly describe predictors in each category below (but see Text S1, Supplemental Material for additional predictor processing details).
Experts agreed habitats used by each species may vary geographically across their ranges, so we created species-specific geographic units by grouping EPA level IV ecoregions and fit an HSM for each unit. We grouped original EPA level IV ecoregions that had similar ecological characteristics, divided them by major geographic features identified in the literature or by experts as barriers or divisions among distinct genetic groups, and ensured final ecoregion groups had sufficient species presences (≥ 25 locations per group).
We created a soil drainage index predictor using drainage classes from gridded SSURGO (raster) data from the Natural Resources Conservation Service. We also obtained the gopher tortoise soil suitability index previously developed by the Natural Resources Conservation Service (2017) to compare this index with the soil drainage index.
We created four vegetation predictors representing compatible land cover, historically disturbed habitat, tree canopy cover, and a deciduous index. We created a compatible land cover predictor by reclassifying types included in the National Land Cover Dataset and Florida Cooperative Land Cover Dataset (Kawula and Redner 2018) as compatible (1: e.g., evergreen forest, shrub/scrub) or incompatible (0: e.g., developed land, cultivated crops) for each species based on expert input and previous habitat selection, movement, and natural history studies (Hermann et al. 2002; Jones and Dorr 2004; Baskaran et al. 2006; Smith et al. 2006; Roznik et al. 2009; Humphries and Sisson 2012; Miller et al. 2012; Beane et al. 2014; Farmer et al. 2017). We created a historically disturbed habitat predictor using historical land cover rasters representing places that were ever classified as developed or agriculture over the period 1938–2001. Studies have found historical land disturbance can lead to suboptimal vegetation communities for at-risk species in the longleaf system even after restoration (Hedman et al. 2000; Kirkman et al. 2004; Brudvig et al. 2013). We included percentage of tree canopy cover and a deciduous index, which evaluated the difference in winter and summer greenness from the Enhanced Vegetation Index, as predictors to capture vegetation conditions.
Disturbance and connectivity.
We developed a fire frequency predictor representing the percentage of years an area burned during 2001–2016. We created this predictor by combining National Aeronautics and Space Administration MODIS data of annual fire detections (NASA 2016b) with U.S. Department of Agriculture Forest Service and U.S. Department of the Interior (2016) LANDFIRE fuel disturbance data to capture burn frequency within this period and compensate for accuracy limitations of each dataset. We created an edge density predictor layer (ratio of edge length to area of compatible land cover) to represent degree of habitat fragmentation.
We used the USGS National Wetlands Inventory database (USGS 2016d) to calculate the mean wetland cover and count of wetlands within a neighborhood. We filtered the National Wetlands Inventory dataset to only include wetlands compatible for the gopher frog and striped newt: wetlands classified as freshwater emergent, freshwater ponds, scrub/shrub, small (< 40-ha) lakes, or forested wetlands not within river floodplains (contained in the USGS National Hydrography Dataset;USGS 2016a).
We included five climate predictors, using University of Idaho Gridded Surface Meteorological Data, representing average conditions in a 30-y period (1981–2010): mean maximum summer (June 1–August 31) temperature, mean minimum winter (December 1–February 27) temperature, and mean annual, summer, and winter precipitation.
We created a topographic position index (TPI) raster using the USGS Digital Elevation Model (USGS 2017). The TPI represents a location's relative elevation to its local surroundings; positive TPI generally indicates ridges (or in the coastal plain, high sandhills), and negative TPI indicates valleys.
Habitat suitability model processing and validation
We built HSMs using generalized linear models (logistic regression) where the presence (1) or absence (0) of a species was the response variable influenced by a set of predictor variables at a location. We chose logistic regression because it has been extensively used for habitat suitability modeling (Elith and Leathwick 2009), can directly incorporate expected ecological relationships (e.g., linear, unimodal: Engler et al. 2004), and can perform similarly to other commonly used algorithms (e.g., MaxEnt: Brotons et al. 2004; Bradter et al. 2018). We used the predicted probability of presence as a habitat suitability index (HSI) ranging from 0 (unsuitable) to 1 (most suitable habitat).
We followed standard modeling practices for presence-background HSMs for each species (Elith and Leathwick 2009; Franklin 2010): we created a dataset with all species presence points and a random sample of pseudoabsence points (drawn from the previously created sampling pools) at a 1:3 ratio within each ecoregion group (following Barbet-Massin et al. 2012), tested predictors for multicollinearity, performed model comparison via Akaike information criterion (AICc, using the AICcmodavg package; Mazerolle and Mazerolle 2019) to select the combination of predictors that best fit the data (Burnham and Anderson 2002), and conducted model validation to measure overall performance and accuracy of predictions. Because species experts were confident that some species–habitat relationships vary by ecoregion, we partitioned the data and performed model selection for each ecoregion group delineated for each species. We performed model selection in two stages to identify the set of predictors that best fit the data. In the first stage, we performed model selection for single variables derived at multiple scales where we compared neighborhood sizes to identify the best-supported scale (e.g., lowest AICc among models comparing soil drainage derived from 90-, 450-, 900-, and 2,000-m neighborhoods). If two variables were correlated, we also compared models identifying the best-supported variable and scale at this stage. In the second stage, we used the full set of uncorrelated predictors at their best-supported scales, included all possible quadratic and 2-way interaction terms among the predictors, and performed AICc-based backwards step-wise regression where predictor effects in the model were dropped if model fit improved (> 2 reduction in AICc: Burnham and Anderson 2002). In trial runs, we performed model selection procedures 10 times for a species where we randomly selected different sets of pseudoabsence points, but the best-fitting model was consistent across iterations. Therefore, we used a single random sample of pseudoabsence points for each species to fit HSMs and predict habitat suitability.
We used five-fold cross-validation to test the performance of each best-fitting model for each species and ecoregion group, and we recorded several evaluation statistics using the ROCR package (Boyce et al. 2002; Johnson et al. 2006). We calculated the area under the curve, obtained by the receiver-operating characteristic plot method (Fielding and Bell 1997), where values of 0.5 indicate model performance no better than random and values > 0.7 indicate the model can acceptably discriminate between sites where the species is present or absent. We calculated sensitivity (the percentage of correctly classified presences), specificity (the percentage of correctly classified pseudoabsences), the true skill statistic (= sensitivity + specificity − 1; Allouche et al. 2006), and the optimal cutoff value (i.e., the HSI value that maximizes accuracy, as calculated by the percentage of presence and pseudoabsence points classified as suitable or unsuitable habitat, respectively). Lastly, we evaluated the influence of individual predictors on habitat suitability. We used hierarchical partitioning analysis (Mac Nally 2002) through the hier.part package (Walsh and Mac Nally 2015) to measure percentage of variance explained by individual predictors. We generated partial response curves to show relationships between each predictor and habitat suitability by ecoregion. We predicted habitat suitability for a given predictor and ecoregion by varying the predictor's value across its range (e.g., 0–1) while holding all other predictors at their means. We did not extrapolate beyond the range of predictor values in our input data. For example, we varied fire frequency between 0 and 60% for the gopher tortoise (the highest value observed in the data), although fire frequencies up to 100% (i.e., burns every year) were possible.
We conducted further model validation where experts reviewed maps of predicted habitat suitability (see below) to qualitatively assess accuracy of preliminary and final models. We first generated preliminary models using the process described above, printed suitability maps, and had experts review these during an in-person workshop. Expert feedback informed us where models may be under- or overpredicting suitability, prompted us in some cases to recombine ecoregion groups that would better capture spatial variation of species–habitat relationships, and highlighted sites with existing data that had not been previously shared by experts. We incorporated feedback (i.e., additional data, modified ecoregion groups) into final models and again had experts qualitatively assess suitability maps for accuracy. All experts approved final model predictions.
Summarizing habitat suitability predictions
We used the best-fitting model to predict and map habitat suitability across each species' range (in a 30-m resolution) using the raster package (Hijmans 2016). To aid in interpreting patterns of suitable habitat, we converted the continuous HSI to classes of unsuitable, low, moderate, and high suitability in ArcGIS. The threshold value separating low and moderate suitability classes for each species represented the mean optimal cutoff value across ecoregion models that maximized accuracy (i.e., the percentage of presence and pseudoabsence points classified as suitable or unsuitable habitat, respectively). We chose thresholds for low and high suitability classes to examine patterns of habitat suitability using more inclusive or restrictive criteria, and these thresholds generally resulted in a 5% reduction in classification accuracy, relative to the optimal cutoff value.
Our experts expressed an interest—one likely shared by conservation practitioners generally—in the distribution of suitable habitat in patches > 1 km2 (100 ha) that could support viable populations and be managed more efficiently than smaller patches. Therefore, we removed smaller patches and created additional maps that only included suitable habitat in patches > 1 km2, and we summarized the area of habitat by suitability class at state- and range-wide scales. Using the > 1-km2 patch map, we further summarized habitat by its current protection status. Patches of suitable habitat were first split in ArcGIS by boundaries of protected areas included in the USGS Protected Areas Database (GreenInfo Network 2016), Florida Natural Areas Inventory Conservation Lands Database (Florida Natural Areas Inventory 2016), Georgia Department of Natural Resources Conservation Lands Database (Georgia Department of Natural Resources 2016), or North Carolina Natural Heritage Program Managed Areas Database (North Carolina Natural Heritage Program 2016). We classified any portion of a patch of suitable habitat contained in a protected area as protected. These areas include publicly owned and -managed lands as well as private lands registered in state or federal programs where natural resource conservation is one of the management goals. We note that these methods classified habitat on U.S. Department of Defense and other multiuse lands as protected. Although these lands are often actively managed for habitats and wildlife species, there is the potential of land use shifts in the future that could result in the loss or degradation of suitable habitat.
Lastly, conservation efforts may be informed by prioritizing areas suitable for multiple species (e.g., Sutherland and deMaynadier 2012; Barrett et al. 2014). Therefore, we reclassified HSIs for each species into binary (compatible/incompatible) classes based on optimal cutoff values (i.e., compatible: moderate and high suitability classes), overlaid binary rasters for each species and summed values to indicate the number of species (0 to 5) with compatible habitat in each cell, and filtered patches so the final map only included patches of compatible habitat > 1 km2. We summarized the area of habitat compatible for one or more species, as well as the percentage in protected areas, by state and across the study extent.
Influence of environmental predictors on habitat suitability
After data filtering, we used between 251 (striped newt) and 15,410 (gopher tortoise) presence records for HSMs for each species (Table 1). Ecoregion-specific HSMs performed well for each species with adequate to high values for the area under the curve ranging between 0.75 and 0.97 (Table 3). Using the optimal cutoff value for a given species across all of its ecoregion-specific HSMs, models exhibited a high degree of accuracy and correctly classified between 82.3% (Florida pine snake) and 93.3% (striped newt) of presence and pseudoabsence points as suitable/unsuitable habitat, respectively (Table 3).
The number of predictors included in best-fitting models across ecoregions for a species ranged between four (gopher frog) and nine (gopher tortoise: Table 4; Figure 2; Figure S2, Supplemental Material). Besides ecoregion, 12 of 16 predictors were included in at least one species' HSMs (Table 4), and all were significant (P < 0.05) in at least one ecoregion (Table S2, Supplemental Material). The four predictors not included in any model were either correlated with other predictors or were not significant predictors of HSI in any ecoregion. The set of predictors included in best-fitting models and their relationships with habitat suitability varied across ecoregions for all species except the striped newt (Table S2; Figure S2, Supplemental Material). Because preliminary ecoregion-specific models for the striped newt did not fit the data better than one that did not allow for ecoregional differences, we used the full dataset to fit a range-wide model. Predictors included in best-fitting models represented average conditions across varying scales (Table S1, Supplementary Material). Soil drainage, compatible land cover, and fire frequency were the only predictors with percentage contributions greater than 10% for all species (Table 4), and, although the degree of their influence varied by ecoregion, these predictors tended to be positively associated with habitat suitability (Table S2, Figure S2, Supplemental Material). For example, HSI for the gopher tortoise increased with soil drainage, compatible land cover, and as fire frequency increased from 0 to 0.5 (one fire every other year; Figure 2). For amphibian species, total area of wetlands in a neighborhood was included in the best model for gopher frogs, whereas the number of wetlands in a neighborhood was in the best model for striped newts. See Text S1, Supplemental Material for additional results.
Spatial patterns of habitat suitability
We present summarized metrics of predicted suitable habitat range-wide for each species in Table 5 (for summarized metrics by state, see Table S3, Supplemental Material). We present maps of gopher tortoise habitat by predicted suitability class in Figure 3, moderate and high suitability habitat in > 1-km2 patches by protection status for the gopher tortoise in Figure 4, and compatible habitat in > 1-km2 patches for all five species in Figure S3, Supplemental Material. The amount of suitable habitat and percentage in protected areas varied by species and state, and Florida and Georgia contained the most suitable habitat for this suite of species (Table 6). Habitat suitability models classified between 8.6% (gopher frog) and 14.6% (gopher tortoise) of each species' range as moderately or highly suitable habitat, and these percentages were lower for habitat in > 1-km2 patches. Between 28.1% (southern hognose snake) and 47.5% (gopher frog) of suitable habitat in large patches was on protected lands (Table 5). Using higher cutoff values resulted in less area classified as suitable for each species, but a higher percentage of highly suitable habitat was contained in large patches on protected lands (range between 36.2% [southern hognose snake] and 60.1% [gopher frog]).
We present summarized metrics for the area of compatible habitat (moderate or high suitability classes) for multiple species and the percentage of habitat on protected lands in Table 6. The number of species whose range extended into each state limited the amount of compatible habitat for multiple species. The percentage of compatible habitat for multiple species that was in protected lands varied between 11 and 70% among states. Habitat compatible for multiple species rather than a single species was more concentrated in protected areas than elsewhere. Habitat suitability models predicted several large areas to be compatible for multiple species, as well as several areas in Georgia and Florida to be compatible habitat for all five species (Figure 5; Figure S4, Supplemental Material). Most large areas were on protected lands, including U.S. Department of Defense lands, national forests, and state-owned conservation areas, with smaller areas of compatible habitat on unprotected lands (Figure 5). Geospatial datasets (raster and shapefiles) of predicted habitat suitability for each species and multiple species are available from U.S. Geological Survey's ScienceBase-Catalog (Crawford et al. 2020; Data A1, Archived Material).
Our study identified important species–habitat relationships and the range-wide distribution and protection status of suitable habitat for five at-risk species in the longleaf pine ecosystem, which can directly inform several conservation decisions. We used a modeling approach that expert judgment informed over multiple stages of input and validation. Because numerous at-risk species are designated for formal status assessments or are the subjects of conservation planning initiatives, our work adds to previous habitat suitability modeling studies to demonstrate an expert-informed approach that produces decision-relevant information about habitat conditions even when available data, such as true absences, are limited.
Across all five species, the most important predictors of suitable habitat were well-drained soil, compatible land cover, and frequent fires, and habitat suitability tended to increase with each predictor. These relationships are consistent with previous findings that the focal species generally favor longleaf pine sandhills, scrub, flatwoods, and other habitats characterized as xeric and fire-dependent (Auffenberg and Franz 1982; Tuberville et al. 2000; Smith et al. 2006; Roznik et al. 2009; Miller et al. 2012; Beane et al. 2014; Farmer et al. 2017). Specifically, predicted suitability for species reached its peak in most ecoregions when fire frequency was above a threshold between 0.5 and 0.2 (Figure 2; Figure S2, Supplemental Material), corresponding to fires occurring once every 2 to 5 years, respectively. This fire-return interval aligns well with prevailing recommendations regarding frequency of prescribed burning to maintain native diversity in longleaf pine and related systems (Platt et al. 1988; Glitzenstein et al. 1995; Russell et al. 1999; Glitzenstein et al. 2003).
For amphibian species, different wetland predictors (total area of wetlands and number of wetlands in a neighborhood) were included in best models for gopher frogs and striped newts, respectively. These two predictors were highly correlated, so these results reflect which wetland predictor outperformed the other in preliminary models. Both species are found in meta-population complexes with multiple breeding ponds (Greenberg 2001; Johnson 2005; Humphries and Sisson 2012), but additional research is needed to understand the influence of size versus number of wetlands on meta-population dynamics for these species. Lastly, relationships between many predictors and HSI varied across ecoregions. This is not surprising as wide-ranging species are often found in multiple ecological settings (e.g., upland ridges and sandhills, mesic flatwoods) and biological communities depending on their geography (e.g., Castellón et al. 2018), which warrants further study and consideration when conducting localized management.
We could not include all attributes that were ranked by experts or known from previous studies to influence habitat suitability for species (Figure S1, Supplemental Material) because the spatial data were not available at the extent or resolution required for our study. Previous studies have found or suggested positive effects of herbaceous groundcover and negative effects of red imported fire ants Solenopsis invicta and feral hogs Sus scrofa on habitat suitability and survival of focal species across their ranges (Tuberville et al. 2000; Johnson 2005; Smith et al. 2006; Roznik et al. 2009; Miller et al. 2012; McIntyre et al. 2019). However, fine-grain spatial datasets are not currently available for these predictors across broader scales. Moreover, we used 30-y averages of climate predictors in models, but additional climatic variables that better capture the intensity and frequency of events like storms or drought may have stronger influences on habitat suitability and species presence (Smith et al. 2006; Blaustein et al. 2010; Farmer et al. 2017). For the two amphibian species, experts ranked wetland hydroperiod as having the greatest influence on habitat suitability (Figures S1d and S1e, Supplemental Material), which is supported by previous studies (Greenberg 2001; Johnson 2005), but there are currently no practical or reliable data layers for hydroperiod, particularly at the extent of our models. Recent studies have developed models predicting hydroperiod of ephemeral wetlands at local scales using climate, soil, and topographic data (Greenberg et al. 2015; Riley et al. 2017); however, incorporation of hydroperiod models over larger spatial extents necessary for conducting range-wide HSMs will require overcoming challenges of computer processing speed and data storage.
As with all spatial models of habitat suitability and species distributions, the accuracy of datasets used as predictors may vary. Remotely sensed data products and large national datasets may contain inherent errors of omission and commission, and some degree of misclassification has been found in the National Land Cover Dataset (Wickham et al. 2018) and the National Wetland Inventory dataset, especially for small, ephemeral wetlands used by focal amphibian species in this study (e.g., Sullivan et al. 2017). We did not perform field verification or reviews of ancillary datasets/aerial imagery to assess the accuracy of these datasets in this study; however, we did find that 98% of our amphibian locations from wetland sampling were on or ≤ 30 m from a wetland in our filtered National Wetland Inventory dataset. Despite any limitations of predictor datasets, our final HSMs exhibited high degrees of accuracy, produced results consistent with previous studies and expert opinion, and can provide reliable information useful for conservation decisions.
Predicted areas of higher habitat suitability tended to highlight known species strongholds, such as national forests, U.S. Department of Defense lands, and other conservation lands. Habitat suitability models predicted that several large protected lands (including Eglin Air Force Base, Blackwater River State Forest, and Ocala National Forest in Florida; Fort Benning and Fort Stewart in Georgia; and the Savannah River Site in South Carolina) would provide suitable habitat for three or more out of five focal herpetofauna species in our study (Figure 5). The amount and distribution of currently protected suitable habitat likely reflects efforts taken in recent years by federal and state agencies and other regional conservation networks (e.g., America's Longleaf Restoration Initiative 2020) to proactively acquire, restore, and manage key areas for habitats occupied by at-risk species, especially with prescribed fire and thinning of upland forests. However, many patches of suitable habitat fell on unprotected lands, and all states had ≤ 55% of their large patches of moderately suitable habitat on protected lands (Table 6). Our results may inform practitioners about specific parcels within their jurisdictions that 1) are highly suitable and protected, motivating continued management; 2) are highly suitable but unprotected, motivating parcel acquisition or enrollment of private landowners in easements or incentive programs; or 3) are less suitable but may serve as reserves or corridors between highly suitable patches if acquired, managed, or restored as necessary (Guisan et al. 2013). Additionally, predicted suitable habitat patches without known species records could be prioritized for survey efforts to document new populations.
An important limitation of our study was that we did not formally evaluate the performance of models generated by our expert-informed approach with others generated without expert input. Previous studies quantitatively assessing model accuracy have found expert-informed models perform variably well compared to empirically driven models (e.g., Pearce et al. 2001; Brandt et al. 2017; Di Febbraro et al. 2018), but these studies vary in how experts informed model construction. Although outside the scope of this study, further studies could quantitatively measure model accuracy when generating background points at different proportions inside and outside expert-informed absence areas and compare these approaches against those using background points drawn at random without constraints. Still, researchers can achieve more accurate predictions of species distributions and habitat suitability with true absence data, which are often unavailable for at-risk species. Increased coordination among federal and state agencies, along with other researchers and managers, to systematically document and store presence–absence data from surveys may be a means to improve the quality of future HSM predictions for at-risk species.
Despite this limitation, we qualitatively observed or can reasonably suspect several benefits from using this approach that overcame challenges inherent in modeling habitat suitability at large (i.e., range-wide) scales, and our approach could be replicated for at-risk, data-limited species. First, we captured expert judgment, along with published studies, to prioritize sets of predictors to test in models that represent key biophysical needs and threats for species. Although some consultation with experts is common in HSM applications (e.g., Murray et al. 2009; Barrett et al. 2014), we took a formal approach using online surveys and in-person discussions that could be easily replicated for other species. This expert input helped identify areas where future research is needed by revealing factors the experts hypothesized to influence species survival and habitat quality but for which published studies or datasets were lacking.
Second, we used expert input to determine appropriate scales for creating averaged predictors that represented best judgments of species' home ranges or dispersal distances when this information was not available from previous studies (e.g., Humphries and Sisson 2012). Third, we used expert judgment to delineate ecoregion groups that represented ecological or genetic variability of species, even when this information was not available in published studies. Using this expert input allowed us to account for spatial variation in species–habitat relationships, which likely increases accuracy of range-wide habitat suitability predictions (Murphy and Lovett-Doust 2007). Fourth, we addressed the challenge of lacking true absence data by using expert judgment to inform where we should generate pseudoabsence points to better capture areas where species do or do not exist in HSMs. Other studies have found HSM prediction accuracy increases with the quality of absence information provided (Brotons et al. 2004; Gu and Swihart 2004), even when models use pseudoabsences that were informed by expert judgment (Murray et al. 2009; O'Leary et al. 2009; Bradter et al. 2018). Our method, which used an interactive mapping platform (Google Earth) to identify suspected areas of species presence and absence when no formal data were collected, adds another tool to those used in previous studies that may be useful for eliciting spatially explicit expert knowledge to aid future applications of HSMs.
Lastly, many of our experts were researchers and managers who can potentially use information about habitat suitability to make decisions, and involving them throughout all stages of the study allowed us to produce trustworthy, decision-relevant results. For example, experts conveyed through in-person discussions during the premodeling stage an interest in knowing where habitat was suitable for multiple species and contained in large (> 1-km2) patches that could allow for cost-effective management actions, which guided how we summarized HSM predictions. It is also reasonable that involvement of experts before, during, and after HSM development increased transparency of modeling tools, experts' understanding of the overall process, and their trust in the results (Martin et al. 2012; Addison et al. 2013; Guisan et al. 2013). Without their involvement early in the model-building process, decision makers may object to or distrust the use of models entirely for a particular context and disregard model outputs when making decisions (Addison et al. 2013). In our study, experts qualitatively evaluated the accuracy of preliminary results, gave feedback that improved model performance, and approved final results. We used varied methods to frequently communicate between groups of scientists, decision makers, and other experts (i.e., e-mail, online exercises, phone calls, in-person workshops), and although it took sustained effort, this degree of expert and decision-maker involvement was beneficial for generating quality, decision-relevant information that can aid at-risk species conservation. Practitioners applying HSMs in the future may agree the benefits outweigh any time- or effort-related costs of collaboration. Given the number of at-risk species in need of habitat assessments and conservation decisions, our methods and other suggested practices provide sound guidance to those leading modeling efforts involving collaborations with decision makers (e.g., Addison et al. 2013; Guisan et al. 2013). These practices include identifying the problem and objectives that models can address, looking for creative opportunities to capture on-the-ground expert knowledge and fill data gaps, and facilitating effective communication between scientists and stakeholders that builds trust of each other, the approach, and the results.
Please note: The Journal of Fish and Wildlife Management is not responsible for the content or functionality of any supplemental material. Queries should be directed to the corresponding author for the article.
Text S1. Supplemental methods and results for at-risk species habitat suitability models (as of 2018) for each of five at-risk species of herpetofauna (gopher tortoise Gopherus polyphemus, southern hognose snake Heterodon simus, Florida pine snake Pituophis melanoleucus mugitus, gopher frog Lithobates [Rana] capito, and striped newt Notophthalmus perstriatus) in the southeastern United States.
Table S1. Neighborhood size (m) of smoothed predictors included in best-fitting models for species habitat suitability (as of 2018) across species' ranges in the southeastern United States. Values indicate the radius of neighborhoods used to average predictors. Results are given for the gopher tortoise Gopherus polyphemus (GT), southern hognose snake Heterodon simus (SHS), Florida pine snake Pituophis melanoleucus mugitus (FPS), gopher frog Lithobates (Rana) capito (GF), and striped newt Notophthalmus perstriatus (SN).
Table S2. Model estimates of predictor effects and model fit statistics for the best-fitting model for at-risk species habitat suitability (as of 2018) for each ecoregion group across the species' range in the southeastern United States. All effects were significant (P < 0.05) unless otherwise noted with a superscripted “NS.” Species-specific results are given for (a) gopher tortoise Gopherus polyphemus, (b) southern hognose snake Heterodon simus, (c) Florida pine snake Pituophis melanoleucus mugitus, (d) gopher frog Lithobates (Rana) capito, and (e) striped newt Notophthalmus perstriatus.
Table S3. At-risk species summary statistics of occurrence records (1880–2018) and habitat suitability area in km2 (percent of total area within summary unit) across the species' range in the southeastern United States. Metrics include area of suitable habitat in > 1-km2 (100-ha) patches and in protected areas
Figure S1. Influential environmental, landscape, and biophysical attributes for each species' suitable habitat and presence at a site (as of 2018), as identified in surveys of 7 to 17 experts per species (N = 27). Results are given for the (a) gopher tortoise Gopherus polyphemus, (b) southern hognose snake Heterodon simus, (c) Florida pine snake Pituophis melanoleucus mugitus, (d) gopher frog Lithobates (Rana) capito, and (e) striped newt Notophthalmus perstriatus. Attributes are generally ordered from highest (top rows) to lowest (bottom rows) influence on habitat suitability and species presence. Definitions for habitat rankings: highly, attributes must occur at a site for the species to be present; somewhat, attributes occurring on the landscape greatly increase the likelihood of species being present, but species may occasionally use landscapes without these attributes; slightly, attributes occurring on the landscape slightly or variably increase the likelihood of species being present, but species may use landscapes without these attributes.
Figure S2. Relationships from the best-fitting model between habitat suitability (as of 2018) and environmental predictors, by ecoregion group within at-risk species' ranges in the southeastern United States (shown in insets). Species-specific results are given for (a) gopher tortoise Gopherus polyphemus, (b) southern hognose snake Heterodon simus, (c) Florida pine snake Pituophis melanoleucus mugitus, (d) gopher frog Lithobates (Rana) capito, and (e) striped newt Notophthalmus perstriatus.
Figure S3. Spatial distribution of suitable habitat (as of 2018) in > 1-km2 (100-ha) patches for at-risk species within the species' range in the southeastern United States. Species-specific results are given for (a) gopher tortoise Gopherus polyphemus, (b) southern hognose snake Heterodon simus, (c) Florida pine snake Pituophis melanoleucus mugitus, (d) gopher frog Lithobates (Rana) capito, and (e) striped newt Notophthalmus perstriatus.
Figure S4. Spatial distribution of habitat (as of 2018) in > 1-km2 (100-ha) patches classified by the number of focal, at-risk species of herpetofauna (out of the five evaluated in this study: gopher tortoise Gopherus polyphemus, southern hognose snake Heterodon simus, Florida pine snake Pituophis melanoleucus mugitus, gopher frog Lithobates [Rana] capito, and striped newt Notophthalmus perstriatus) for which habitat was predicted to be compatible (moderate or high suitability). The grey background represents the full study extent across species' ranges in the southeastern United States.
Reference S1. Enge KM, Farmer AL, Mays JD, Castellón TD, Hill EP, Moler PE. 2014. Survey of winter-breeding amphibian species. Final report. Florida Fish and Wildlife Conservation Commission, Fish and Wildlife Research Institute, Gainesville, Florida.
Found at DOI: https://doi.org/10.3996/092019-JFWM-075.S2 (8.99 MB PDF); also available at https://f50006a.eos-intl.net/ELIBSQL12_F50006A_Documents/14_Enge.pdf.
Reference S2. Natural Resources Conservation Service. 2017. Gopher tortoise soil suitability. U.S. Department of Agriculture.
Found at DOI: https://doi.org/10.3996/092019-JFWM-075.S3 (1.03 MB PDF); also available at https://www.nrcs.usda.gov/Internet/FSE_DOCUMENTS/stelprdb1101742.pdf.
Reference S3. Outcalt KW, Sheffield RM. 1996. The longleaf pine forest: trends and current conditions. Pages 28 in Resource Bulletin SRS-9. USDA Forest Service, Southern Research Station, Asheville, North Carolina.
Reference S4. Sutherland R, deMaynadier P. 2012. Model criteria and implementation guidance for a Priority Amphibian And Reptile Conservation Area (PARCA) system in the USA.
Found at DOI: https://doi.org/10.3996/092019-JFWM-075.S5 (1.61 MB PDF); also available at http://parcplace.org/wp-content/uploads/2017/08/PARCA_System_Criteria_and_Implementation_Guidance_FINAL.pdf.
Reference S5. [USFWS] U.S. Fish and Wildlife Service. 2016. Long-term listing transformation and five guiding principles of the unified listing team. Director's memorandum 16 March 2016.
Found at DOI: https://doi.org/10.3996/092019-JFWM-075.S6 (98 KB PDF).
Please note: The Journal of Fish and Wildlife Management is not responsible for the content or functionality of any archived material. Queries should be directed to the corresponding author for the article.
To cite this archived material, please cite both the journal article (formatting found in the Abstract section of this article) and the following recommended format for the archived material.
Data A1. Geospatial datasets (raster and shapefiles) of predicted habitat suitability (as of 2018) for each of five at-risk species of herpetofauna (gopher tortoise Gopherus polyphemus, southern hognose snake Heterodon simus, Florida pine snake Pituophis melanoleucus mugitus, gopher frog Lithobates [Rana] capito, and striped newt Notophthalmus perstriatus) in the southeastern United States. Available: https://doi.org/10.5066/P92PZN7G
This work is funded and supported by the U.S. Fish and Wildlife Service through the former Gulf Coastal Plains and Ozarks Landscape Conservation Cooperative (Cooperative Agreement Award F16AC00985), the Natural Resources Conservation Service (Award 68-7482-16-545 and Agreement NR183A750001C001), and the Southeast Climate Adaptation Science Center (Cooperative Agreement Award G16A00402). We thank Mike Harris and D. Todd Jones-Farrand for supporting the development of this work. We thank all individuals who generously contributed data and expert judgment that facilitated this work from the North Carolina Wildlife Resources Commission, North Carolina Natural Heritage Program, North Carolina Museum of Natural Science, Davidson College, South Carolina Department of Natural Resources, University of Georgia Savannah River Ecology Laboratory, Francis Marion University, Georgia Department of Natural Resources, Joseph W. Jones Ecological Research Center, Florida Fish and Wildlife Conservation Commission, Archbold Biological Station, Coastal Plains Institute, Alabama Department of Conservation and Natural Resources, Auburn University, Mississippi Wildlife, Fisheries, and Parks, Louisiana Department of Wildlife and Fisheries, South Atlantic Landscape Conservation Cooperative, Department of Defense, U.S. Army (Fort Stewart), U.S. Air Force (Eglin Air Force Base), Longleaf Alliance, The Nature Conservancy, Partners in Amphibian and Reptile Conservation, and HerpMapper. We thank the editors and three anonymous reviewers for their comments that substantially improved the manuscript. The Georgia Cooperative Fish and Wildlife Research Unit is jointly sponsored by U.S. Geological Survey, the University of Georgia, U.S. Fish and Wildlife Service, Georgia Department of Natural Resources, and the Wildlife Management Institute.
Any use of trade, product, website, or firm names is for descriptive purposes only and does not imply endorsement by the U.S. Government.
Citation: Crawford BA, Maerz JC, Moore CT. 2020. Expert-informed habitat suitability analysis for at-risk species assessment and conservation planning. Journal of Fish and Wildlife Management 11(1):130–150; e1944-687X. https://doi.org/10.3996/092019-JFWM-075
The findings and conclusions in this article are those of the author(s) and do not necessarily represent the views of the U.S. Fish and Wildlife Service.