A Large-Scale MaxEnt Model for the Distribution of the Endangered Pygmy Madtom

Current distribution maps and knowledge of habitat requirements at multiple spatial scales provide critical information for rare species conservation by guiding management actions such as protecting habitat throughout a species' range. However, obtaining comprehensive distributional records and habitat data can be a daunting task for conservation biologists when a species is difficult to sample, little is known about its habitat-use patterns, and it occupies a large geographic area. Lomolino (2004) described this lack of species distributional knowledge on multiple scales as the “Wallecean shortfall,” one of several key knowledge gaps limiting conservation of biodiversity.

Species distribution models (SDMs) are often used to describe potential habitats for rare species in need of conservation management (Liang et al. 2012). These models provide insight into 1) overall distribution for a particular species, 2) historic and future occurrence probabilities for the species, and 3) practical understanding of niche limits of the species (Aguirre-Gutierrez et al. 2013). In short, SDMs spatially connect species presence records at a location with environmental attributes at those sites (Elith et al. 2011).

MaxEnt software, created by Phillips et al. (2006), applies a technique called maximum entropy as a modeling approach for conservation research. MaxEnt is a machine learning method that uses presence-only data to determine habitat suitability in the form of relative probabilities of species occurrence (Liang et al. 2012; Merow et al. 2013). Relative occurrence suitability scores may be derived from MaxEnt results with the assumption that the detection suitability score for the species of interest remains constant throughout the sampling range (Yackulic et al. 2013). The predictive performance of MaxEnt is competitive with other high-performing methods, and the MaxEnt modeling approach is used by many government and nongovernment organizations (Elith et al. 2011). For example, Endries (2011) used MaxEnt to create predictive habitat maps for 247 freshwater aquatic species in North Carolina, by comparing species presence data to a collection of stream and land cover environmental variables (Endries 2011). Liang et al. (2012) provided an example of the use of MaxEnt to model headwater stream habitat for the threatened Blackside Dace Chrosomus cumberlandensis. Their results suggest that natural and anthropogenic disturbances form constraints on the fish's habitat. Liang et al. (2012) also provided maps that display habitat suitability and clustered “hotspots” for the Blackside Dace that are useful to conservation managers. Managers are interested in siluriform catfishes, which represent a diverse array of aquatic animals with a widespread, near-global distribution. As with many groups of fishes, the development and use of SDMs (including MaxEnt) serve a variety of purposes. For example, Roderiguez-Rey et al. (2019) used MaxEnt and four other SDMs to model distribution scenarios for aquatic invasive species in Great Britain, including the Wels Catfish Silurus glanis. In South America, Frederico et al. (2014) used MaxEnt to predict the distribution of six Amazonian fish species, including the cetopsid catfish Helogenes marmoratus. In North America, studies describe the use of SDMs for at least 10 ictalurid catfish species, including several madtom species in the genus Noturus (Endries 2011; Bouska et al. 2015; Huang et al. 2016).

The Pygmy Madtom Noturus stanauli is federally protected as an endangered species under the U.S. Endangered Species Act (ESA 1973, as amended), with occurrence records currently known from a small section of the Upper Clinch River and a larger section of the lower middle Duck River within the Tennessee River Drainage (Figure 1; Etnier and Jenkins 1980; Etnier and Starnes 1993). These rivers are separated by 1,055 river kilometers, leading managers to ask whether there may be other populations within the Tennessee River Basin. The Pygmy Madtom has the ability to bury itself under gravel and rely on interstitial water during the day, which may have made detection difficult in the past (Etnier and Starnes 1993). Since its federal listing, research has been performed on the Pygmy Madtom, including studies on genetic relationships (Bennett et al. 2009), population structure and longevity (Wells and Mattingly 2019), benthic fish community associations (Wells and Mattingly 2020), microhabitat use in the Clinch River (Wells et al. 2020), and detection by using environmental DNA techniques (Paine et al. 2021). Phylogenic and population genetic analyses by Bennett et al. (2009) revealed surprisingly little molecular divergence between the populations in the Clinch and Duck rivers, at least for the loci examined (Wells and Mattingly 2020). Data availability is still limited on Pygmy Madtom, with only 20 known occurrence locations, and no documented applications of SDMs on this species.

Our study therefore used a MaxEnt modeling approach to predict the distribution of the Pygmy Madtom throughout the larger Tennessee River, which includes the Clinch and Duck rivers of Tennessee by using presence-only data made available from previously conducted research. The use of MaxEnt allowed for an SDM approach, even with limited sampling data, to generate suitability scores that may guide further research efforts. The specific objective of our study was to develop SDMs to predict and prioritize suitable habitat for the Pygmy Madtom in rivers that fall within the Tennessee River Basin by using presence-only occurrence data within a MaxEnt framework. This information will be beneficial to guide future research efforts of state and federal biologists and other managers who are working to monitor, protect, and prevent the extinction of the Pygmy Madtom. Furthermore, our findings are relevant for other researchers seeking to better understand large-scale environmental variables of potential importance to other madtom species.

The Pygmy Madtom is known to only occur in Tennessee within two river systems: the Clinch River and the Duck River (Figure 1). The Clinch River is located within the Ridge and Valley, Level III Ecoregion of eastern Tennessee and southwestern Virginia. This river contains one of the most diverse fish and mussel assemblages in North America (Starnes and Starnes 1985). The Duck River is located entirely within the Interior Plateau, Level III Ecoregion of Central (Middle) Tennessee. The Duck River also contains a diverse aquatic community of 66 mussel species, 32 of which occur downstream of Centerville (Schilling and Williams 2002).

The original range of the Pygmy Madtom is unknown, but a compilation of historical findings of the species was possible from reports and sampling efforts from the Tennessee Department of Environment and Conservation and the Tennessee Valley Authority (Etnier and Jenkins 1980; Bennet et al. 2009; Wells and Mattingly 2019). This compilation resulted in 20 site locations of known presence of the Pygmy Madtom for use as the input for sample locations in the MaxEnt model.

The Pygmy Madtom inhabits medium-sized rivers, where it has been collected from shallow shoals with moderate-to-strong current over a substrate of pea-sized gravel or fine sand (Etnier and Starnes 1991; U.S. Fish and Wildlife Service [USFWS] 1994; Wells and Mattingly 2019). We chose environmental variables that may allow for more insight into the range of the variables driving Pygmy Madtom habitat suitability, from expert-based understanding of the available data (Etnier and Starnes 1993; Pearce et al. 2001; Johnson and Gillingham 2004; Bennet et al. 2009; Wells and Mattingly 2019). We selected and obtained environmental variable data from the National Hydrography Dataset Plus Version 2 (NHD) (USGS 2019) and the 2011 National Land Cover Database (NLCD) dataset. We used the Region Tennessee (06) within the NHD. We used ComID (common identifier of the NHD area feature) and FeatureID (ComID of the associated NHD flowline) fields to join data attribute files to catchment layers. Environmental variables included canopy cover, drainage area, flow, precipitation, slope, stream order, temperature, and velocity measured at the segment scale of Frissell et al. (1986; Table 1). Flow rate was the mean annual flow (cubic feet per second [cfs]) at the bottom (downstream edge) of the flowline using the Vogel method (Vogel et al. 1999). Velocity was the mean annual velocity (feet per second) at the bottom of flowline by using the Jobson method (Jobson 1996). We determined stream order from modified Strahler classifications (McKay et al. 2012). Slope was the slope of the flowline based on smoothed elevation; we initially calculated slopes as meters per kilometer but then converted them to unitless meters per meter (McKay et al. 2012). Temperature was the mean annual air temperature in degrees Celsius × 100, and precipitation was the mean annual precipitation within the segment, upstream of the bottom of the flowline in millimeters × 100. We defined drainage area as the catchment drainage area in square kilometers (AreaSqKm), that is, the land area around the segment draining directly into that segment (not the cumulative drainage area upstream in the watershed; McKay et al. 2012). Canopy cover was the mean canopy cover in a 30-m buffer surrounding the stream segment. We obtained canopy cover data from the 2011 U.S. Forest Service dataset, and we generated the buffer by using the buffer tool in ArcGIS Pro® 2.9.1 (ESRI, Redlands, CA) with a 30-m buffer set to dissolve based on the FeatureID of the flowline. We calculated mean canopy cover for each segment within ArcGIS Pro 2.9.1 by using the focal statistics tool.

We used geographic information system (GIS) layers as inputs for the MaxEnt model to analyze the relative importance of the variables. We conducted all data processing in ArcGIS Pro 2.9.1. We used MaxEnt 3.4.4 to run the analysis. MaxEnt software required all attribute variables be in American Code for Information Interchange (ASCII) format with matching cell size, geographic projection, and extent. Therefore, we used a standard grid-cell size of 30 × 30 m, with a projection of NAD 1983 (2011) StatePlane Tennessee FIPS 4100 (meters) and a raster file of the state of Tennessee as the processing extent and snap raster. We joined the NHDPlusV2 variables to the catchment shapefile. For canopy cover, we processed the NLCD 2011 U.S. Forest Service tree canopy analytical layer by projecting to NAD 1983 (2011) StatePlane Tennessee FIPS 4100 (meters) and setting the processing extent and applying a snap raster. We created a 30-m buffer around the stream flowline to represent the riparian zone. We calculated the mean canopy cover within the 30-m buffer by using the zonal statistics table tool. We then joined the table to the catchment shapefile. This approach allowed for a stream-segment level of analysis (Frissell et al. 1986) with uniform grid-cell values throughout each stream segment catchment for each environmental variable. Finally, we used the polygon-to-raster tool to produce a raster for each variable in the model, and then we converted each raster into an ASCII file (.ASC) by using the raster-to-ASCII tool. We measured correlation between environmental variables within ArcGIS Pro 2.9.1 by using the Band Collection Statistics tool and the correlation matrix output that contained the sample correlation coefficients. To calculate the critical value for the test of hypothesis, we used an α significance level of 0.05 with 18 df (n − 2), resulting in a critical value of 0.378. The null hypothesis is that there is not a significant correlation between two variables. The alternate hypothesis is that there is significant correlation between two variables. Any correlation coefficient greater than the critical value of 0.378 results in a rejection of the null hypothesis and acceptance the alternate hypothesis that the two variables are correlated. We found that stream order and velocity correlated with each other and with flow; therefore, we removed stream order and velocity from the model due to high correlation. We then ran the MaxEnt model with the remaining variables (canopy cover, drainage area, flow, precipitation, slope, and temperature) for the entire geographic range of the Tennessee River Basin.

The MaxEnt modeling approach entailed the following steps and software settings. We determined relative estimates of probabilities of occurrence by allowing the software to use a complementary log-log (cloglog) transformation (Phillips et al. 2017). To calculate percent contribution of each variable to the model, we added the increase in regularized gain for each of a maximum of 500 iterations of the training algorithm to the contribution of the corresponding variable, or subtracted it if the absolute value of lambda (allow regularization) was negative (Phillips et al. 2006). The permutation importance for each variable was determined through that variable's grid-cell value on training presence and background data being randomly permutated. We then reevaluated the model on permutated data, and this resulted in the reduction of training area under curve (AUC). AUC is typically used to determine how well a model is able to distinguish between presence and absence locations; however, in the case of MaxEnt modeling with presence-only data, the AUC is a measure of the probability that the model will rank a presence location higher than a background location (Merow et al. 2013). We normalized the AUC values to percentages. We used jackknife statistical testing procedures to help understand the relative importance of individual variables (Elith et al. 2011; Phillips 2017).

Finally, we constructed maps based on model forecasts for suitable habitat. The known geographic range of the Pygmy Madtom lies within the Upper Clinch River and the Upper and Lower Duck River. We identified river segments outside of the known geographic range as other potential areas within the Tennessee River Basin for Pygmy Madtom occurrence. We identified and marked catchments with a suitability score greater than or equal to a 99% probability by using ArcGIS Pro 2.9.1 for the purpose of prioritizing segments for future field sampling to detect the species in additional locations within its known geographic range.

Pygmy Madtoms were present in a range of environmental conditions in the Clinch and Duck rivers (Table 2). The four strongest variables that contributed to the model were canopy cover, drainage area, flow rate, and slope. The mean canopy cover within the 30-m riparian buffer was 34.81%. The mean drainage area in was 2.22 km². The mean flow rate was 2,634.4 cfs (697.57 m³/s), and the mean slope was 0.00036 from occurrence data. The training AUC was 0.996 when run with canopy cover, drainage area, flow, precipitation, slope, and temperature. The four strongest contributors to the multiple-variable MaxEnt model were flow rate, canopy cover, drainage area, and slope (Table 3). By contrast, temperature and precipitation were relatively minor contributors, with a percent contribution of 3% or less. For flow rate modeled by itself, estimated suitability score increased to a maximum near 4,000 cfs (113.27 m³/s) and then tapered off at higher flows (Figure 2). For canopy cover modeled by itself, the estimated suitability score increased rapidly to a peak at approximately 22% canopy coverage and then declined gradually as coverage increased. For drainage area modeled by itself, the estimated suitability score declined as drainage area increased, as was the same for slope. We noted very similar responses in the marginal response curves for canopy cover, drainage area, and slope when modeling these variables in a multiple-variable setting. When modeled with all variables, estimated suitability score did not decrease in flows greater than 4,000 cfs, but instead continued to gradually increase.

We removed the two least-contributing variables, temperature and precipitation, from the model, and the MaxEnt model was run again with the four strongest contributing variables—flow rate, canopy cover, drainage area, and slope in that order (Table 3). We calculated the final suitability score of occurrence for the Pygmy Madtom by this final model run with canopy cover, drainage area, flow, and slope. Flow and canopy cover remained the strongest contributing variables, with flow having the greatest percent contribution of 68.5%. The jackknife test also revealed that the variable with the highest gain when used in isolation was flow (Figure 3). The environmental variable that decreases the training gain most when the variable omitted was also flow, showing that flow appears to have the most information that is not present in the other model variables. The suitability score of occurrence for the Pygmy Madtom in the Upper Clinch River is shown in Figure 4, and the suitability score of occurrence for the Pygmy Madtom in the Duck River is shown in Figure 5, with the catchments having the highest suitability score of 99% or above denoted by a red pin. Warmer colors indicate a higher suitability score. One hundred and seventy-six of the catchments within the Tennessee River Basin NHDPlusV2 dataset had a suitability score of 99% or greater. Of the 176 catchments, 85 fall within the Upper Duck and Lower Duck river watersheds. Eight of the 176 catchments fall within the Upper Clinch and Lower Clinch river watersheds. Eighty-three catchments with a greater than 99% suitability score fall within watersheds outside of the known range of the Pygmy Madtom. Other rivers with the highest density of suitable catchments for Pygmy Madtom occurrence include stretches of the Nolichucky River, the French Broad River, the Hiwassee River, and the Elk River.

We used MaxEnt model outputs (probabilities of Pygmy Madtom occurrence) to produce suitability maps based on selected environmental variables for the Tennessee River Basin. We gained insight regarding areas most suitable for Pygmy Madtom occurrence within its known distribution within the Clinch and Duck rivers as well as areas outside of the known range where Pygmy Madtom occurrence may be possible. MaxEnt heat maps resulted in colors of intensity with corresponding values based on suitability score of occurrence.

Because of sampling difficulty and the size of the rivers that the Pygmy Madtom inhabits, an SDM was an appealing approach to predict the distributional probabilities for this species. Using a modeling approach such as MaxEnt was very helpful given the impossibility of thoroughly sampling these systems (Liang et al. 2012). Our study produced successful “management-friendly” products (maps and tables) for state and federal agencies similar to the products generated by Liang et al. (2012) in their study on the Blackside Dace in the Upper Cumberland River in Kentucky. Our results provide an up-to-date distribution of predicted habitat locations for the Pygmy Madtom in the Clinch and Duck rivers. River stretches identified by MaxEnt will be useful for future research efforts by agencies and universities. The occurrence maps highlight where new sampling efforts could start throughout the Tennessee River Basin and throughout the known range of the Clinch and Duck rivers. Future MaxEnt efforts need to extend into Virginia because of the possibility of Pygmy Madtom occurrence in the Clinch River extending into that state. Discovering populations beyond the historic sampling site at Frost Ford would be encouraging for the Pygmy Madtom in the Clinch River. The Nolichucky, French Broad, Hiwassee, and Elk rivers also have stretches of high suitability scores for the Pygmy Madtom. Historical presence data are limited for the Pygmy Madtom and limit the capabilities of the model to predict the Pygmy Madtom distribution. However, these suitability predictions will facilitate more efficient survey efforts to determine whether there may be additional populations of the Pygmy Madtom throughout the Tennessee River Basin, and the model results may be refined through future surveys.

Based on MaxEnt results and preliminary microhabitat results by Wells et al. (2020), the Pygmy Madtom was more likely to occupy nearshore habitats in November such as shallow runs associated or near channel bars found in reaches of medium-sized rivers. Changes in mean annual flow rates and mean annual velocities appear to alter the suitability score of occurrence of the Pygmy Madtom, as seen in the response curves of our results. Alterations of river reaches into reservoirs throughout the Tennessee Valley occurred during the 20th century. We noted alterations in flow rates and velocities because of these impoundments, and the results were eradication of fish species and fragmentation of fish populations, which might explain the Pygmy Madtom's current disjunct distribution to the mainstem of the Clinch and Duck rivers. Localized Pygmy Madtom populations may have been more widely distributed historically. Populations in the Tennessee River mainstem and its tributaries may have been extirpated with the rise of anthropogenic modifications, rather than geologic events (Bennett et al. 2009). Bennett et al. (2009) disagree with the hypothesis that the Upper Tennessee River separated from the Duck River during the Cretaceous and Tertiary millions of year ago, because of the current low genetic divergence between Pygmy Madtom populations in the Clinch and Duck rivers. Unfortunately, freshwater ecosystems worldwide suffer many of the same broad threats: anthropogenic activities that cause habitat degradation, fragmentation, and loss; flow modifications; invasive species; overexploitation; and pollution (Jelks et al. 2008). Our study is the first to use MaxEnt to predict suitable habitat reaches and distribution for the Pygmy Madtom. Response curves for each of the environmental variables provide insight into the factors that may be driving Pygmy Madtom occurrence (Figure 2). Suitability scores of Pygmy Madtom increased as mean canopy cover in the 30-m riparian zone buffer increased, up to the threshold of approximately 22%, after which there was a decrease in the suitability score. These results reflect the observed pattern that Pygmy Madtoms have been found along channel bars, most abundantly in areas of intact riparian buffers to minimize sedimentation and loss of habitat heterogeneity (Wells et al. 2020). Pygmy Madtoms were less likely to occur if flow rate increased above a mean annual flow rate of 4,000 cfs. Pygmy Madtoms were also less likely to be found as catchment area or slope increased. Increases in catchment area or slope are related to the river continuum concept (Vannote et al. 1980), which conceptually describes how streams exhibit a continuous gradient of physical characteristics and available energy that correspond to biological communities (Endries 2011). The Pygmy Madtom has only been reported in Tennessee, but due to high suitability scores upstream into Virginia and throughout other river drainages in Tennessee, sampling efforts could be conducted in additional river drainages with high suitability scores. Important results from the MaxEnt maps are new locations to examine within the known range throughout the Clinch River upstream of Norris Reservoir and in the Duck River from Hurricane Mills to Columbia in Maury County. In addition, rivers throughout the Tennessee River Basin that have high probability for suitable Pygmy Madtom habitat have been identified for future sampling efforts. These rivers could have other populations of Pygmy Madtoms that may have been connected before habitat fragmentation.

Please note: The Journal of Fish and Wildlife Management is not responsible for the content of functionality of any supplemental material. Queries should be directed to the corresponding author.

Data S1. Environmental variables (flow rate, canopy cover, drainage area, temperature, precipitation, slope, stream order, and velocity) used in our MaxEnt modeling approach to determine the distribution of the Pygmy Madtom Noturus stanauli in the Clinch (n = 5) and Duck (n = 11) rivers, Tennessee. Flow rate is the mean annual flow in cubic feet per second at the bottom of the flowline by using the Vogel method, and canopy cover is the mean percent canopy cover in the 30-m riparian buffer for each stream segment. Drainage area is the area of the catchment in square kilometers. Temperature is the mean annual temperature in degrees Celsius multiplied by 100. Precipitation is the mean annual precipitation in millimeters multiplied by 100. The slope is unitless meters per meter. Stream order is the modified Strahler method for classification of stream order. Velocity is mean annual velocity in feet per second at the bottom of the flowline by using the Jobson method. Canopy cover data are from the 2011 National Land Cover Database dataset. All other values are from the National Hydrography Dataset Plus Version 2 (USGS 2019).

Available: https://doi.org/10.3996/JFWM-21-057.S1 (26 KB DOCX)

This project was supported by USFWS; Tennessee Wildlife Resources Agency; the School of Environmental Studies, Tennessee Technological University (TTU); Center for the Management, Utilization and Protection of Water Resources (TTU); and the Department of Biology (TTU). We thank the technicians, undergraduates, and graduate students from TTU who assisted with this project. We would also like to thank the reviewers and the Associate Editor for their very helpful comments and suggestions that helped to improve the paper.

Any use of trade, product, website, or firm names in this publication is for descriptive purposes only and does not imply endorsement by the U.S. Government.