Assessing Conservation Potential of Streams for Spikedace and Loach Minnow Using Species Distribution Modeling

Decades of aquatic habitat loss, degradation and modification, and introduction of nonnative species present formidable challenges for managers tasked with protecting and restoring native fish communities. These challenges have proven especially difficult to overcome in the arid Gila River basin of Arizona and New Mexico where sparse stream flow has been further altered and reduced because of a legacy of diversions, dams, and groundwater pumping (Rinne et al. 2005). In addition, nonnative fishes have been introduced throughout the basin and have reduced the distribution and abundance of native fish communities (Olden and Poff 2005; Rinne et al. 2005). Consequently, conservation of native fishes endemic to the Gila River basin has relied on strategies that prioritize the least-altered remaining river reaches for conservation, including translocation and reintroductions (Clarkson and Marsh 2010; Hickerson et al. 2021)

Spikedace Meda fulgida and Loach Minnow Rhinichthys cobitis are two cyprinids endemic to the Gila River basin of Arizona and New Mexico that have greatly declined in distribution largely because of human-induced impacts to stream systems (Propst et al. 1986; Propst and Bestgen 1991; Minckley and Marsh 2009). Both species were subsequently federally listed as endangered (USFWS 2012) pursuant to the U.S. Endangered Species Act (ESA 1973, as amended). Both fish are habitat specialists that occupy perennial, low-gradient, intermediate-size streams with relatively intact flow regimes that are largely free of nonnative fishes (Propst et al. 1986; Propst and Bestgen 1991; Rinne 1991). These focal species occupy a relatively broad elevation range (approximately 500–2,100 m for Spikedace and 500–2,500 m for Loach Minnow), and consequently water temperature is not thought to be a primary determinant of their distribution (Propst and Bestgen 1991; Marsh et al. 2003). Both species often occupy the same streams but have relatively specific mesohabitat requirements, with Spikedace preferring runs, eddies, and shearwaters and Loach Minnow preferring riffles with unembedded cobbles (Propst and Bestgen 1991; Rinne 1991).

The ecological processes that constrain the distribution of stream fishes operate within a series of hierarchical scales, ranging from environmental characteristics of entire watersheds at the coarse scale to microhabitat at the fine scale (Poff 1997). The relatively specific habitat requirements at multiple spatial scales along with the historically patchy distribution of Spikedace and Loach Minnow is an additional challenge for managers tasked with restoring populations of these species. With relatively unaltered streams becoming increasingly rare, it has become a priority to understand the distribution and characteristics of streams with the habitat to potentially support these target species. Species distribution modeling is increasingly being used to understand historic and contemporary distributions of stream fishes, the primary ecological drivers of distributions, and to prioritize areas for conservation (Oakes et al. 2005; Pool et al. 2013; Booher and Walters 2021a, 2021b). These models can be used to identify specific stream reaches with environmental conditions similar to those historically occupied by focal species as candidates for further conservation efforts and assessments of habitat at finer spatial scales.

The objective of this study was to identify potential translocation sites for Spikedace and Loach Minnow in the Gila River basin. We used a classification tree approach (Random Forest), informed by historic fish sampling records within the entire historic range for each focal species, and several abiotic variables thought to be important determinants of species distributions at the reach scale, to locate stream reaches with the highest predicted probability of Loach Minnow and Spikedace occurrence. The results of this study may prove valuable for managers looking to direct limited resources to streams with a greater chance for conservation successes.

The Gila River is a tributary to the Colorado River and drains an area of approximately 212,380 km² (Rinne et al. 2005). The hydrology of the Gila River basin is similar to other river basins of the southwestern United States in that there are relatively few perennial reaches, which tend to be patchily distributed on the landscape (Jaeger et al. 2014; Goodrich et al. 2018). Consequently, the flow regimes of many higher-order streams within the basin are classified as either ephemeral or intermittent (Costigan et al. 2016). The lack of reliable stream flow records and potential historical changes in flow regime makes it difficult to assess the occurrence of perennial reaches, and, consequently, streams potentially supporting fish populations across the entire Gila River basin. The study area for each focal species was constrained to their respective historical ranges in the Gila River basin, including all tributaries (except for the Santa Cruz River) upstream of the confluence with the Agua Fria River for Spikedace, and all tributaries upstream of the confluence with the Salt River for Loach Minnow (Figures 1 and 2; Minckley and Marsh 2009).

We utilized hydrography from the National Hydrography Dataset Plus (NHDPlus) High Resolution. Value-added attributes, which include tabular data associated with each NHD reach that enhance analysis, data display, and stream networking, were configured consistent with the NHDPlus High Resolution user guide (Moore et al. 2019). Similarly, we obtained discharge and velocity estimates for each reach from 1971 to 2000 from the enhanced unit runoff method table (EROM) and joined estimates to the hydrography. The enhanced unit runoff method utilizes estimates of mean annual runoff produced by the U.S. Global Change Research Program that are adjusted using excess evapotranspiration estimates, man-made additions and removals of flow, and reference stream gauge data (Moore et al. 2019 p.56–63). We constrained hydrography for the Gila River basin (hydrologic unit code 1504–1507) to the historical range of each species (Data S1 and S2). All reaches not classified as intermittent, perennial, or artificial path were removed from the data set to ensure that ephemeral reaches and man-made conveyances (e.g., canals) were not included in analyses. In addition, we excluded all first-order reaches because neither focal species is known or thought to occur in first-order reaches. We removed a small number of reaches (about 1%) with missing values for any attributes that could not easily be imputed from the data set.

We assembled species occurrence point data sets from the lower Colorado River basin Aquatic Gap Analysis Project data set, which contains fish survey records from many different sources within the lower Colorado River basin (Pool et al. 2010; Data S3 and S4). We used fish collection records from 1890 to 2006 to characterize the historical geographic range of the species. Although stream habitat conditions have certainly changed within this period, our goal is to evaluate habitat availability for Spikedace and Loach Minnow within the context of environmental gradients at present day. We created presence data sets from the LCRB data set by spatially joining all records for each focal species to the nearest NHDPlus reach. To ensure an exact spatial overlap with hydrography we used a snapping tolerance (i.e., search distance) of 50 m to make minor adjustments to the spatial location of point records. Similarly, we compiled absence data sets by spatially joining all records that did not match one of the focal species to the nearest NHDPlus reach using the same procedure. We manually checked snapped locations of points using locality description information and corrected where necessary to ensure that points snapped to the correct reach. Both the presence and absence data sets for each focal species were reduced to a single point per NHDPlus reach to avoid pseudoreplication. Point shapefiles were manipulated in ArcGIS Pro (v. 2.4.2).

We chose a total of seven variables to include in our models. Six of these variables were chosen because of known species-habitat relationships for Spikedace and Loach Minnow (e.g., slope, elevation, estimated discharge, estimated velocity), or the ability to describe stream size (stream order, total upstream drainage area; Table 1; Propst et al. 1986; Propst and Bestgen 1991). We derived variables for each NHD reach from the NHDPlus High Resolution value-added attributes and enhanced unit runoff method tables available with NHDPlus High Resolution hydrography (Table 1). The location and extent of perennial stream reaches is an important constraint on the distribution of stream fishes in the Gila River basin. However, these reaches can be difficult to characterize using NHD attributes alone. In the Gila River basin these perennial stream reaches are often characterized by dense, intact riparian areas (Nguyen et al. 2014). We chose to use normalized difference vegetation index (NDVI) values derived from satellite imagery to identify reaches with dense riparian areas as an additional surrogate for perennial reaches (Goodrich et al. 2018). In addition to providing a measure of vegetation greenness, NDVI values can help distinguish wetted areas (e.g., NDVI values −1 to 0). We calculated mean NDVI values within a 30-m buffer of each NHD reach using Landsat imagery from 1990 to 1992 (Landsat 5). We used imagery from this era as a representation of historical conditions because several Spikedace and Loach Minnow populations that were still extant during this period (e.g., Upper Verde River, Eagle Creek) have potentially become extirpated. We chose a 30-m buffer to match the resolution of the imagery and the relatively narrow riparian areas associated with most streams.

We assessed the probability of Spikedace and Loach Minnow historical occurrence within the Gila River basin using a classification model approach (Random Forest). Random Forest models are widely used for species distribution modeling because they can achieve accurate predictions while accounting for common problems in ecological data analyses (i.e., spatial autocorrelation and nonlinear interactions; Evans et al. 2011). The Random Forest model is a classification tree method that relies on bootstrap replicates (i.e., out of bag [OOB] samples) for internal validation. Specifically, we used a class-balanced Random Forest model (package rfUtilities) because the number of presence points for each species was imbalanced relative to the number of absence points (Evans and Cushman 2009). We performed variable selection using the model improvement ratio (function rf.modelSel, package randomForest), which minimizes mean squared error of the model, limits the number of retained predictors, and maximizes the percentage of variation explained by iteratively assessing the relative importance of each covariate toward improving overall model fit (Murphy et al. 2010). We tuned model hyperparameters using the function tuneRF (package randomForest) to determine the number of candidate variables drawn at each split (i.e., mtry = 1). We ran each model with an increasing number of trees until OOB error converged between 1,500 and 1,700. Models were run with focal species occurrence (present = 1, absent = 0) as the response variable. We fit separate models for Spikedace and Loach Minnow because of known differences in distribution and habitat requirements between the species.

We assessed model accuracy, which is the proportion of records predicted correctly by the model (i.e., sum of presence reaches classified as 1 and absence reaches classified as 0, divided by the total number of reaches in training data set). We also calculated area under the receiver operating characteristic curve (AUC), which provides a measure of the proportion of true positives relative to false negatives (Fawcett 2006). Guidance in Evans et al. (2011) suggests that values near 0.5 indicate no classification ability, whereas values of 1.0 indicate perfect discrimination. We evaluated cross-classification error using the kappa statistic (Cohen 1968; Evans and Cushman 2009). In addition, we evaluated model fit with OOB error estimates derived from withheld data for each bootstrap replicate (Evans and Cushman 2009). Because of the limited number of presence points, we did not split fish occurrence data into training and validation data sets for additional cross-validation procedures. We evaluated model significance (P; α ≤ 0.05) following the randomization test procedure outlined in Murphy et al. (2010) with 999 permutations. For each focal species, we applied the final model to modified hydrography for the historic range to obtain a predicted probability of occurrence at each NHD reach. We carried out all analyses in Program R (R Core Team 2021, version 4.0.4).

There were a total of 115 presence and 1,116 absence points (imbalance ratio = 1:9.70) for Spikedace (Figure 1). The Loach Minnow data set contained 136 presence and 1,045 absence points (imbalance ratio = 1:7.68; Figure 2). Overall predicted probability of Spikedace occurrence ranged from 0 to 0.96, with predicted probability at historically occupied reaches ranging from 0.46 to 0.96 (Figure 3, Table 2). Classification performance metrics suggest that the Spikedace model predicted species occurrence moderately well (OOB error = 0.24, accuracy = 76.32, κ = 0.48, AUC = 0.75). Variable importance scores suggest that total upstream watershed area was the most important variable, followed by estimated discharge, estimated velocity, slope, and stream order (Table 3). Predicted probability of Spikedace occurrence generally followed expectations for the species, with higher probability of occurrence associated with reaches with low slopes, higher stream order (i.e., farther distance from the stream source), and relatively high estimated velocity (Figure 4). Predicted probability of occurrence was influenced more subtly by upstream watershed area and estimated discharge.

Overall predicted probability of Loach Minnow occurrence ranged from 0 to 0.98, with predicted probability at historically occupied reaches ranging from 0.45 to 0.98 (Figure 5, Table 4). Classification performance metrics suggest that the Loach Minnow model predicted species occurrence with less accuracy than the Spikedace model (OOB error = 0.25, accuracy = 74.59, κ = 0.44, AUC = 0.72). Variable importance scores suggest that estimated discharge was the most important variable, followed by total upstream watershed area, elevation, and stream order (Table 3). Probability of Loach Minnow occurrence also followed expectations for the species, with positive associations for reaches with relatively small upstream watershed area, moderate elevation and stream order, and relatively low discharge (Figure 6). Both the Spikedace and Loach Minnow models were significant (P = 0.001) and fit the data well.

Our application of species distribution models to assess conservation potential of streams for endangered fishes in the Gila River basin aligns with Spikedace and Loach Minnow habitat associations suggested from previous studies (Propst et al. 1986; Propst and Bestgen 1991), while providing new information to potentially guide their management. The predictor variables included in the model were useful for elucidating coarse habitat associations and identifying potentially suitable reaches for translocations but were limited in their ability to accurately predict species occurrence back to the landscape. Considering that only a handful of attributes available with NHD High Resolution data set were utilized, the final models performed surprisingly well at classifying stream reaches (i.e., AUC ≥ 0.72). Our approach may be useful for managers attempting to improve the conservation status of other native fish species because of the ready availability of NHD data and long-term fish sampling data sets available for other species.

Relatively low classification accuracy of the Random Forest models suggests that the patchy historical distribution of both Spikedace and Loach Minnow proved relatively difficult to accurately predict at the reach scale. Classification accuracy may have been slightly lower for Loach Minnow because of the broader range of environmental conditions and stream sizes represented in reaches occupied by this species. Despite the relatively low classification accuracy, predicted probability of occurrence from these models may still be useful for guiding further on-the-ground conservation efforts. Given that the final models utilized relatively few NHD variables, other factors that we were not able to account for likely influence species distribution at the reach scale, and finer spatial scales may have a more important role in regulating distribution. These results are not surprising, as the distribution of these fishes is likely heavily constrained by habitat conditions at the microhabitat scale and by biotic interactions with both native and nonnative fishes (Propst et al. 1986; Propst and Bestgen 1991). Other mechanisms potentially contributing to low classification accuracy may be the relatively narrow temporal range of variables included in the model, bias in the historical data set, and contemporary changes to stream systems (e.g., reduced discharge), which influence important habitat characteristics like riparian cover.

The predictor variables we selected limit our ability to learn detailed information on the ecology of these fishes, which was not the primary goal of this study; however, we were able to obtain information about general habitat associations for both species at the reach scale. Probability of occurrence of both species was most influenced by variables that describe watershed position or stream size. Loach Minnow tend to occupy mid-watershed reaches, as they are less likely to occur in headwater tributaries (i.e., stream order < 3), in main-stem rivers with large watersheds and high annual discharge, or at elevational extremes (both low and high). Spikedace are similar to Loach Minnow in their tendency to occupy reaches in mid-size streams, but are less sensitive to watershed size or discharge and show stronger associations with low-gradient reaches that still have relatively high velocity. Both species tend to follow the “Goldilocks principle” where smaller, headwater streams farther upstream in watersheds and larger, main-stem river farther downstream have a low probability of occurrence, but reaches that have intermediate watershed position and size with specific flow characteristics in the middle are “just right.”

Although classification accuracy of these models was relatively low, predicted probability of occurrence for both species was relatively high for reaches in historically occupied streams. For example, there was at least one reach with a predicted probability of occurrence > 0.46 in every stream historically occupied by each focal species. In general, predicted probability of occurrence in historically occupied streams was lower for streams with few, spatially clustered records (e.g., Cibecue Creek, Salt River, Tonto Creek [Spikedace], Bonita Creek, Pace Creek, Beaver Creek [Loach Minnow]), and highest in streams with relatively abundant and more evenly distributed populations (e.g., Aravaipa Creek). Notably, predicted probability of occurrence was also relatively high in several streams without historical records of focal species. Unoccupied reaches with the highest predicted probability of occurrence may exhibit ecological similarity to historical and extant habitat for Spikedace and Loach Minnow. These areas may have high conservation value due to their potential to support translocated populations of focal species. A recent successful translocation project supports this idea. Spikedace were historically absent from the Blue River, yet the model predicted a maximum probability of occurrence of 0.84. An effort to translocate Spikedace to the Blue River starting in 2012 resulted in the establishment of a robust population of Spikedace in reaches with the highest predicted probability of occurrence (Hickerson et al. 2021). These results suggest that unoccupied reaches with the greatest predicted probability of occurrence may be important to prioritize for further conservation efforts.

The relationship between predicted probability of occurrence and the ability to support a translocated population of either target species should be informed by on-the-ground assessments of stream habitat and fish community composition in addition to model results. For example, two streams with similar maximum predicted probability of Spikedace occurrence (Big Chino Wash [P = 0.86] and the San Simon River [P = 0.82]) to the Blue River are typically dry for most of the year and do not currently support fish populations. These cases highlight limitations of geospatial data sets like the NHD, the importance of ground-truthing model predictions, and the potential role of factors and processes regulating fish distributions at finer spatial scales (Poff 1997; Durance et al. 2006). We emphasize that modeling at the reach scale should be thought of as an initial step in the process of prioritizing stream reaches for conservation. On-the-ground assessments of macrohabitat and fish community composition are necessary to further prioritize reaches most likely to support translocated populations of target species using existing approaches (Hickerson and Walters 2019; Booher and Walters 2021b).

Although the lower Colorado River basin data set is one of the best available sources of historical distribution information for the Gila River basin, the data set is susceptible to bias resulting from nonstandardized sampling, spatial resolution, and errors of omission, which has potential implications for model results. Historical sampling effort typically was not standardized or evenly distributed across the landscape, and some populations may have been extirpated before detection (Olden and Poff 2005). In addition, these small-bodied species may have been overlooked or not recorded in early surveys because of their small size and may have been present but not detected at some reaches (Marsh et al. 2003). These potential sources of uncertainty about the true historic distribution of these species could affect model classification accuracy in situations where cases classified as false positives may actually be true positives that remained undetected in the historical record (Strecker et al. 2011). Understanding the limitations of the data set and the consequences for classification accuracy metrics is important to consider when reviewing model predictions.

Biotic interactions, particularly with nonnative fishes, are often highly important determinants of the distribution and abundance of stream fishes (Olden and Poff 2005; Cucherousset and Olden 2011). Spikedace and Loach Minnow are no exception to this pattern, as their distribution and abundance are constrained by interactions with other fishes, particularly nonnative fishes (Propst et al. 1986, 2014). Nonnative fishes gradually replaced Spikedace and Loach Minnow over time in some locations, which could result in few to no records at these locations and has potential consequences for model prediction (Olden and Poff 2005). Unfortunately, nonnative fish distribution data are not available for the Gila River basin at the resolution necessary to inform our modeling approach. Consequently, we did not include fish community information in our models despite their known importance as a constraint on distribution of the target species.

In addition to the introduction of nonnative fishes, other changes to stream environments have influenced the distribution of both species during the historical record. Habitat fragmentation and alteration accompanied by the creation of dams and diversions has hastened the spread of nonnative fishes and altered hydrologic processes over time (Olden and Poff 2005; Pool et al. 2010). Contemporaneous groundwater depletion, flow reduction, and stream diversions have eliminated or reduced habitat that may have once been suitable for fish (Ruhí et al. 2016). Consequently, stream reaches that may have historically been capable of supporting both focal species may now lack the perennial flow necessary to support fish populations. Our study highlights the importance of maintaining and updating long-term data sets that can be leveraged to inform conservation strategies for current stream conditions using modern analytical methods.

Our results can be used as a first step for locating and prioritizing reaches most likely to support translocated populations of Spikedace and Loach Minnow within their respective historical ranges. More broadly, our study represents a useful exercise to evaluate conservation potential of stream reaches for endangered fishes across a large river basin. Within the lower Colorado River basin fisheries managers could use these same data sources (Aquatic Gap Analysis Project and NHD) to develop similar exercises when considering diverse conservation actions such as habitat restoration, translocation, fish inventories, and population monitoring. Further development of our approach could incorporate additional hydrologic and climatic data that can be used to assess the relative suitability of candidate translocation sites under variable climate scenarios (Booher and Walters 2021b). Because Spikedace and Loach Minnow share many reaches with relatively high predicted probability of occurrence, opportunities likely exist to improve the status of both species in the same location, which could allow managers to use limited resources for conservation more efficiently. The presence of nonnative fishes within high predicted-probability-of-occurrence reaches remains an important consideration for managers planning conservation efforts, including fish passage barrier construction, nonnative fish removals, and translocation of native fishes (Hickerson et al. 2021). We recommend using existing frameworks for on-the-ground assessments of stream habitat and fish community composition as a necessary next step for further prioritizing stream reaches for conservation opportunities (Galloway et al. 2016; Hickerson and Walters 2019). Many southwestern native fishes are facing similar challenges to Spikedace and Loach Minnow and are also in need of conservation efforts. Our approach may be an applicable first step for other species of conservation concern with available historic records in need of population replication.

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.

Data S1. National Hydrography Dataset Plus High Resolution flow lines and value added attributes for Spikedace Meda fulgida historical range in the Gila River basin as of 2022. Also included is estimated probability of Spikedace occurrence for each National Hydrography Dataset reach.

Available: https://doi.org/10.3996/JFWM-21-093.S1 (54.408 MB ZIP)

Data S2. National Hydrography Dataset Plus High Resolution flow lines and value added attributes for Loach Minnow Rhinichthys cobitis historical range in the Gila River basin as of 2022. Also included is estimated probability of Loach Minnow occurrence for each National Hydrography Dataset reach.

Available: https://doi.org/10.3996/JFWM-21-093.S2 (18.361 MB ZIP)

Data S3. Species occurrence (presence, absence) point data sets for Spikedace Meda fulgida from 1890 to 2006 were assembled from the lower Colorado River basin Aquatic Gap Analysis Project data set and spatially joined to the nearest National Hydrography Dataset Plus High Resolution flow lines.

Available: https://doi.org/10.3996/JFWM-21-093.S3 (126 KB ZIP)

Data S4. Species occurrence (presence, absence) point data sets for Loach Minnow Rhinichthys cobitis from 1890 to 2006 were assembled from the lower Colorado River basin Aquatic Gap Analysis Project data set and spatially joined to the nearest National Hydrography Dataset Plus High Resolution flow lines.

Available: https://doi.org/10.3996/JFWM-21-093.S4 (120 KB ZIP)

We thank Dr Shannon Albeke for providing excellent courses in spatial data analysis that helped lay the foundation for this present work. We also thank the reviewers and Associate Editor for helpful comments to improve this manuscript.

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.