A Hierarchical Approach to Fish Conservation in Semiarid Landscapes: A Need to Understand Multiscale Environmental Relationships

Abstract

A multiscale perspective is essential for conservation planning of riverine fishes. Coarse-scale habitat (e.g., basis) can influence both finer-scale habitat characteristics (e.g., reaches and microhabitat) and associated species distributions. Finer-scale management and habitat rehabilitation efforts can fail without the consideration of coarser-scale constraints. We provide a conceptual hierarchical framework for multiscale fish conservation strategies in the semiarid Great Plains. The Great Plains stream network is highly fragmented due to dam construction, water withdrawals, and increased drought severity. Our framework uses relationships with basin-scale connectivity and streamflow and reach-scale physicochemical characteristics in the context of aiding species reintroduction and stream habitat improvements.

1. Introduction

The importance of multiscale habitat use by aquatic organisms is well recognized and central to the development of meaningful fisheries conservation actions. The distribution of fishes relies on natural features such as the appropriate climate and geology that comprise the physicochemical characteristics that are typically tolerated by species. These ultimate determinants constrain intermediate and proximate determinants on aquatic organisms (e.g., benthic algae, [1]). For example, the pH of a river is dictated, in part, by the underlying lithology of the region [2], and fishes have specific pH tolerances that regulate a variety of life-history attributes (e.g., egg hatching in salmonids [3]). A combination of other physicochemical factors at finer spatial and temporal scales contribute to a heterogeneous riverscape that shape species assemblages [4, 5], where fish use a set of variables that are assumed to maximize fitness [6] or describe behavioral responses (e.g., changes in cover use, [7]). The habitat needs of fishes are often used as the foundation of conservation and recovery plans [8]. For example, priority use areas can be identified and restoration actions planned. As human pressures on fish populations increase, establishing multiscale relationships is more important than ever for guiding conservation actions.

1.1 Landscape change and anthropogenic pressures on fish populations

Human pressures increase the threats on freshwater ecosystems and taxa. The modification of landscapes from historical land cover to agriculture and urban uses has resulted in significant physicochemical changes and water demands on rivers. Urban rivers, for example, are often associated with flashy hydrographs, increased contaminants, and degraded channel morphology [9] including channel incision and erosion [10]. Agriculture land use is also associated with higher sediment and nutrient loads [11, 12], more homogenous substrates and water depths [13], and bank instability [14, 15]. Pressure on water resources needed for human uses has resulted in rivers being dammed, leveed, and pumped thereby disrupting both flow and sediment regimes [16, 17]. The magnitude, duration, timing, frequency, and rate of change of stream discharge (hereafter flow regime) is considered a primary driver of ecosystem processes [18, 19] and biotic integrity [20, 21]. Flow regimes globally have been altered due to human activities [19, 22, 23]. Current and future changes to our climate and associated weather patterns will only exacerbate threats facing freshwater ecosystems.

North American freshwater fishes are experiencing the highest extinction rates among vertebrates [24]. Flow regime alteration and fragmentation due to dam construction (hereafter damming) are often cited as primary causes [22, 25, 26]. Dams alter fish-assemblage structure [27, 28, 29] and prey availability [30, 31, 32] by creating streamflow conditions and instream habitats favorable to lentic (i.e., lake and reservoir) species. The changes to flow magnitude, in particular, due to damming, have negatively affected many fishes leading to declines in diversity [22, 24, 26]. Altered flow regimes and fragmentation caused by damming are particularly detrimental to lotic (i.e., river and stream) fishes due to their mobility and requirement of multiple habitat types to complete a life cycle [4]. Dams disrupt spawning cues and block migration routes, prevent access to spawning and nursery habitats, and alter nutrient cycles [26, 27, 32]. Damming has largely ceased in North America, but the long-term effects of modified ecosystems result in numerous stream-fish conservation challenges [22, 33, 34]. Existing levels of flow regime alteration and fragmentation are also exacerbated due to climate change and additional anthropogenic pressures, particularly in arid and semiarid regions.

1.2 Climate change and multiscale fish conservation in arid and semiarid landscapes

Flow-regime alteration and stream fragmentation are expected to increase due to climate change and growing human water demands. The combination of extended periods of drought and increased human water demands magnify threats to the long-term persistence of many stream fishes [35, 36, 37]. Native fish species in arid and semiarid ecoregions are particularly vulnerable because they have both the highest level of damming [33, 38] and naturally harsh environmental conditions [39, 40, 41]. Fishes native to arid and semiarid streams have evolved to tolerate intermittent drying, flashy changes is flow, and extreme physicochemical conditions [38, 41, 42]. However, interactions of dams, loss of base flows due to water withdrawals, and increased drought have intensified already harsh environments [37, 38, 43]. Identifying multiscale stream-fish relationships with fragmentation and flow regimes is essential to developing successful conservation strategies in arid and semiarid ecoregions. Understanding coarse-scale constraints on species distributions is particularly important to guiding finer-scale conservation and recovery efforts. For example, the stream reach is a natural scale for fish conservation and management activities, but effective strategies are dependent on basin-scale relationships [5, 44].

1.3 Great Plains small-bodied minnows

Streams of the semiarid North American Great Plains ecoregion (hereafter Great Plains) are unique, dynamic ecosystems and home to endemic fish species. The Great Plains is one of the most impounded areas of the world [33, 38, 45]. Damming largely ceased in the 1980s, but resulting fragmentation is extensive and flow regimes remain severely altered [33]. Substrate in Great Plains streams is predominately sand and silt, with constantly changing streambed formations [44]. Channel characteristics differ from gravel-bed streams with stable riffle-run-pool formations except during high flows. Environmental conditions are extreme with periods of long drought followed by large flooding and highly variable and wide-ranging water temperatures and salinities [42, 46]. Natural stream drying has been exacerbated by harsher drought periods due to climate change and excessive groundwater pumping [37, 42, 43]. Small-bodied fishes are common in the Great Plains, with brackish (i.e., salt-tolerant) species dominating areas of higher salinity [46, 47]. True minnows (Leuciscidae), particularly smaller-bodied species, also occur in varying abundances throughout the ecoregion. True minnows are a large, diverse family of fishes (∼700 species) within the order Cypriniformes [48, 49]. In particular, true minnows display a wide range of life-history traits. This includes a group of species, some endemic to the Great Plains, that have a unique reproductive strategy in stream ecosystems.

Pelagic-broadcast spawning [50] is a common reproductive strategy globally in marine and coastal species, but rare in inland freshwater systems. In inland streams of North America, pelagic-broadcast spawning is restricted to mooneyes (Hiodontidae) and three genera of true minnows (Hybognathus, Notropis, and Macrhybopsis, hereafter pelagophils) that occur in the Great Plains [51]. Pelagophils typically spawn in higher-order streams and release transparent, non-adhesive ova that are semi-buoyant [47]. The downstream displacement of eggs and larvae relies on drift [52]. Thus, both minimal fragmentation and higher flow magnitude are essential to successful recruitment [53, 54, 55]. There are similarities between pelagophils and marine pelagic-broadcast spawners in juvenile dispersal strategies to microhabitats with high nutrient concentrations and reduced predation pressure [51, 56]. High-flow events increase nutrient loads and create disconnected temporary slackwater habitats that serve as nurseries for juvenile pelagophils. The pelagophil life cycle is completed through extensive upstream movement by juveniles and adults [53, 57, 58]. Great Plains pelagophils have been strongly negatively affected by damming, water withdrawals, and climate change due to disrupted stream networks, altered flow patterns with reduced magnitude, and loss of floodplain habitats [42, 52, 59]. Numerous studies have reported declines in pelagophil relative abundance and range reductions in the Great Plains (e.g., [43, 54, 60, 61, 62]).

2. Modeling at different spatial scales

We modeled occurrence and relative abundance of true minnows (hereafter minnows) in the upper Red River basin. Species occurrence (i.e., distributions) and abundance (i.e., population size) are fundamental ecological state variables used both in research studies and for conservation and management problems [63, 64, 65]. Both fish distributions and population sizes may be constrained by coarser-scale characteristics [1, 5, 44]. One analytical approach is to model variation in occurrence or abundance as a function of environmental variables at multiple scales (e.g., [66, 67, 68]). However, state variables are quantified differently (e.g., a binary response for occurrence and integers for abundance); thus, it is typically not possible to model multiple states in a single analysis. Different state variables also more naturally align with certain spatial scales. For example, abundance is often not ecologically meaningful at very coarse scales (e.g., basins) and measuring and managing population size at these scales is not realistic. Therefore, occurrence is typically examined at the basin scale, where the stream reach (i.e., a series of representative habitat complexes nested in tributary complexes) is the natural scale for studies and management of abundance [5, 44]. An alternative multiscale analytical approach (our approach here) is to model variation in occurrence and abundance separately at relevant spatial scales.

2.1 Study area

Our study area was in the upper Red River basin of the Great Plains. The area comprised portions of the Central Great Plains and Southwestern Tablelands sub-ecoregions of Oklahoma and Texas (Figure 1). The eastern boundary with the Cross Timbers sub-ecoregion corresponded with a transition from sand-bed to gravel-bed streams with lower salinity, increased vegetation, and fish species hybridization zones [69, 70, 71]. The western boundary corresponded with the higher-elevation, more-arid High Plains sub-ecoregion. The Central Great Plains is characterized by mixed-grass prairie vegetation, cropland, and landforms that include sand dunes, low mountains, and salt flats [70]. The Southwestern Tablelands has a more rugged terrain, with dissected plain, hill, and canyon landforms, sparse short grass prairie vegetation, and less cropland. Annual precipitation in the study area, though highly variable, increases from east to west eastward (mean rainfall 56–97 cm). In addition to minnows, the brackish plains killifish Fundulus zebrinus and Red River pupfish Cyprinodon fluviatilis are also common in the study area.

Figure 1.

Study area in the upper Red River basin. The shading in the upper panel denotes level-three sub-ecoregion boundaries (from west to east: High Plains, Southwestern Tablelands, Central Great Plains, and Cross Timbers). Inner borders on polygons show the delineation of hydrologic response units included in the occurrence modeling (Section 2.2). White circles are stream reaches surveyed for prairie chub in 2019 (Section 2.3). The polylines show the Red River mainstem (thicker line) and select major tributaries: A is the North Fork, B is the Pease River, C is the Wichita River, and D is Cache Creek.

2.2 Small-bodied minnow occurrence among hydraulic response units

We modeled occurrence probability of nine minnows species among hydraulic response units (HRUs) nested in the upper Red River basin (Figure 1). HRUs are sub-basins that represent 10-digit hydrologic units with refined boundaries for flow modeling based on local characteristics [72, 73]. The focal species included four pelagophils (emerald shiner Notropis atherinoides, plains minnow Hybognathus placitus, prairie chub Macrhybopsis australis, and Red River shiner N. bairdi, [52]), bullhead minnow Pimephales vigilax, fathead minnow P. promelas, red shiner Cyprinella lutrensis, sand shiner N. stramineus, and suckermouth minnow Phenacobius mirabilis. Most species occur elsewhere in North America east of the Rocky Mountains, with emerald shiner, fathead minnow, and red shiner widely distributed (www.iucnredlist.org). Red River shiner is endemic to most of the Red River basin and introduced to other basins of Arkansas, Oklahoma, and Texas [61]. Prairie chub, a species of conservation concern (see Section 2.3), is restricted to the upper Red River basin [69, 71]. We did not include minnows that only occur near the ecotone with the Cross Timbers.

Our study period was 2002–2015. The temporal range encompassed a relatively dry climatic period (2002–2014) and one year of heavy rainfall (Figure 3 in [74]). We assumed static species occurrence states (i.e., no colonization or extirpations of HRUs) across the study period and at least a one-year lag time for any changes in minnow occurrence states following the end of the dry period.

2.2.1 Fish surveys

We compiled stream-fish surveys from state conservation and management agencies, data repositories, and online databases (Appendix 1, Table A1). For online databases, we used the terms “fish” and “fishes” to search all Oklahoma and Texas counties within the study area from 2002–2015. Data were processed to remove duplicate surveys. Each unique survey was spatially referenced to an HRU based on the latitude and longitude. We compiled capture histories for each species (i.e., one for detection and zero for nondetection) at each HRU. Repeat surveys at HRUs were treated as spatial replicates with replacement [76]. We also compiled the date, sampling gear (reported in 92% of surveys), and collector (e.g., agency or individual, Appendix 1, Table A2). We assumed seining for surveys that did not report the gear type because it is the most-common stream-fish sampling method [77] and comprised the majority of reported surveys (64%). The additional gear types were backpack electrofishing (18% of surveys) and boat or barge electrofishing (9% of surveys).

2.2.2 Flow regime and fragmentation metrics

We characterized the flow regime of each HRU using mean daily discharge estimates at the outlet. The discharge estimates were obtained from a precipitation-runoff modeling system [72] adapted from the National Hydrologic Model [78] for the Red River basin [79, 80]. We calculated flow regime metrics [81, 82] using EflowStats (version 5.0.1, median option, [83]). Due to the size of the dataset and inherent correlations among metrics, we limited the set to five with each metric representing one flow regime component (Table 1). The selected flow metrics were based on expected ecological relationships with minnow occurrence, particularly pelagophils (see Section 1.3), and maintaining the absolute value of Pearson’s pairwise correlation coefficient r < 0.5.

Metric	Description
Low flow duration variability (DL17, %)	CV of annual Q below 25^th PCTL
†Flood pulse count (FH3, # of d/yr)	Median of annual # of days Q is above 3 * median daily Q
†Median daily flow (MA2, m^3/s)	Median daily Q for entire flow record
†Variability of reversals (RA9, %)	CV of flow reversals
Annual maxima variability (TH2, %)	CV of Julian day of annual maximum Q

Table 1.

Flow-regime metrics used in the species occurrence model.

^†

natural-log transformed prior to analysis due to a right-skewed distribution.

The alphanumeric codes for metrics correspond to detailed descriptions in Appendix 7 of Kennen et al. [84].

Q: stream discharge.

PCTL: percentile.

We used the upstream network density of major dams (UNDR, #/100 rkm) at the outlet of each HRU obtained from an online database [85] to characterize fragmentation. We used UNDR because it was highly correlated with numerous other fragmentation metrics and upstream dams have been shown to be strongly associated with pelagophil distributions in the upper Red River basin [56]. We natural log transformed UNDR prior to modeling due to a right-skewed distribution. Correlation levels were reasonable between UNDR and flow regime metrics (|r|<0.30).

2.2.3 Occupancy modeling methods

We modeled minnow occurrence relationships, while accounting for variable detection probability, using the hierarchical framework described by [86]. The latent occurrence state for minnow i at HRU j was treated as partially observed, with z_ij = 1 if the species was truly present and z_ij =0 if the species was truly absent. Each z_ij followed a Bernoulli distribution with occurrence probability Ψ:

$${\mathit{z}}_{\mathit{ij}}\sim \text{Bernoulli}\left({\mathrm{\Psi }}_{\mathit{ij}}\right)$$E1

The detection of species i at HRU j for survey k was conditional on both the true occurrence state and detection probability p (the probability of detecting a species in a single survey if present), where y_ijk followed a Bernoulli distribution:

$${\mathit{y}}_{\mathit{ijk}}\sim \text{Bernoulli}\left({{\mathit{z}}_{\mathit{ij}}}^{\ast }{\mathit{p}}_{\mathit{ijk}}\right)$$E2

We modeled variation in Ψ and p as a function of covariates [86]. Detection covariates comprised HRU surface area (km², hereafter area) and drainage area. Spatially replicated surveys can result in a violation of the closed-system assumption for occupancy modeling because a species may occur at a site, but not be locally present at the time of the survey [76]. Thus, we used area to account for variation in p associated with patchier species distributions in larger HRUs. Drainage area characterized the stream order of the mainstem for each HRU to account for variation in p associated with species abundance. We natural-log transformed both detection covariates due to right-skewed distributions. Detection relationships with covariates were allowed to vary by species as deflections around the group mean hyperparameter governed by a probability distribution [64, 87, 88]. More common minnows may have inherently higher detection probability. Thus, we modeled the correlation (ρ) between species occurrence probability intercepts α_i and species detection probability intercepts υ_i. The intercepts were jointly distributed as [α_i , υ_i | Σ] ∼ N(0, Σ), where Σ is a 2 x 2 matrix comprising variance components ${\mathbf{\sigma }}_{\mathit{\alpha }}^{\mathbf{2}}$ and ${\mathbf{\sigma }}_{\mathit{v}}^{\mathbf{2}}$ and covariance σ_ασ [88, 89]. We also allowed each species detection intercept to vary by both sampling gear type g (1, seining, 2,backpack electrofishing, 3, boat and barge electrofishing) and collector c (1–6, Appendix 1, Table A2) using a grouping factor [90, 91]. The detection component of the occupancy model can be written as:

$$\text{logit}\left({\mathit{p}}_{\mathit{\text{ijkgc}}}\right)={\mathrm{\upsilon }}_{\mathit{i}}+{\mathrm{\beta }}_{1\mathit{i}}{\mathit{X}}_{1\mathit{jk}}+{\mathrm{\beta }}_{2\mathit{i}}{\mathit{X}}_{2\mathit{jk}}+{\mathrm{\gamma }}_{\mathit{ig}}+{\mathrm{\tau }}_{\mathit{ic}},\text{for}\mathit{i}=9,\text{for}\mathit{j}=97,\text{for}\mathit{k}=1,2\dots \mathit{K},$$E3

1. β_ni ∼ Gaussian (μ_βn, ${\mathrm{\sigma }}_{\mathrm{\beta }n}^2$), for n = 2

2. γ_ig ∼ Gaussian (0, ${\mathrm{\sigma }}_{\mathrm{\gamma }}^2$), for g = 3

3. τ_ic ∼ Gaussian (0, ${\mathrm{\sigma }}_{\mathrm{\tau }}^2$), for c =6

where υ is the detection probability intercept, β₁ and β₂ are slopes for associated detection covariates area X₁ and drainage area X₂, γ is the gear type grouping factor, τ is the collector grouping factor, and μ is species group mean. Occurrence covariates comprised UNDR and five flow regime metrics (see Section 2.2.3). We allowed each occurrence covariate to vary by species using the same model structure described for the detection component of the model. The occurrence component of the occupancy model can be written as:

$$\text{logit}\left({\mathrm{\Psi }}_{\mathit{ij}}\right)={\mathrm{\alpha }}_{\mathit{i}}+{\mathrm{\beta }}_{\mathit{ni}}{\mathit{X}}_{\mathit{nj},}\text{for}\mathit{i}=9,\text{for}\mathit{j}=97,\text{for}\mathit{n}=6,{\mathrm{\beta }}_{\mathit{ni}}\sim \mathit{\text{Gaussian}}\left({\mathrm{\mu }}_{\mathrm{\beta }\mathit{n}},{\mathbf{\sigma }}_{\mathbf{\beta }\mathit{n}}^{\mathbf{2}}\right),$$E4

where α is the detection probability intercept and β₁–β₆ are slopes for associated occurrence covariates X₁–X₆. All detection and occurrence covariates were standardized to a mean of zero and standard deviation (SD) of one. Model coefficients are reported as the mode (most likely value) with a 90% highest density interval (HDI, [92]). Model specifications, diagnostics, and fit tests are provided as supplemental information (Appendix 2).

2.2.4 Occupancy modeling results

There were varying detection relationships among minnows. Detection probability across all minnows (i.e., the group mean) at mean levels of covariates was 0.47 (90% HDI: 0.31, 0.64). Detection probability was lower in larger HRUs for all minnows, with the strength of the relationship similar to the group mean (Appendix 1, Table A3). The strength of the detection relationship with drainage area was higher than the group mean for emerald shiner, plains minnow, prairie chub, and Red River shiner and lower for the remaining five minnows. There was more unexplained variation in detection probability attributed to collector (SD = 0.84) than gear type (SD = 0.25). Detection probability and occurrence probability intercepts were moderately positively correlated (ρ = 0.63).

The direction of occurrence relationships with flow-regime characteristics and fragmentation varied both among all minnows and within pelagophils (emerald shiner, plains minnow, prairie chub, and Red River shiner). Occurrence probability across all minnows at mean levels of covariates was 0.76 (90% HDI: 0.60, 0.89).

Occurrence probability for all pelagophils increased sharply with increasing daily streamflow magnitude (MA2, Table 2 and Figure 2a). There was a weak positive occurrence relationship with MA2 for bullhead minnow and suckermouth minnow, a weak negative relationship for sand shiner, and no relationship for fathead minnow and red shiner. Occurrence probability decreased with increasing variability of reversals (RA9) for all minnows, including the cosmopolitan red shiner (Table 2 and Figure 2b). Emerald shiner occurrence probability also increased sharply with increasing variability in annual maxima (TH2, Table 2 and Figure 2c). The positive occurrence relationship with TH2 was weaker for bullhead minnow, fathead minnow, red shiner, sand shiner, and suckermouth minnow. Prairie chub and Red River shiner had a weak negative occurrence relationship with TH2, while plains minnow had no relationship. The group mean and all minnow occurrence relationships were similar with low flow duration variability (DL17) and flood pulse count (FH3). Occurrence probability increased with both increasing DL17 and FH3 (Table 2). Plains minnow, prairie chub, and Red River shiner occurrence probability decreased sharply with increasing upstream dam density (UNDR, Table 2, Figure 2d). However, there was no occurrence relationship with UNDR for emerald shiner. Sand shiner occurrence was also negatively associated with UNDR, but the strength of the relationship was weaker than that for the three pelagophils. Suckermouth minnow and bullhead minnow occurrence probability increased sharply with increasing UNDR (Figure 2d). Fathead minnow occurrence was also positively associated with UNDR, but the strength of the relationship was weaker. There was no occurrence relationship with UNDR for red shiner.

Species	Intercept	UNDR	DL17	FH3	MA2	RA9	TH2
Mean	1.1; (0.4, 1.8)	−1.1; (−3.1, 0.9)	0.5; (0.2, 0.9)	0.4; (0.1, 0.9)	0.8; (–0.1, 1.9)	–1.0; (–1.6, –0.4)	0.3; (–0.4, 0.8)
BUM	2.0; (1.3, 2.8)	1.3; (0.4, 2.2)	0.5; (0.1, 1.0)	0.5; (0.1, 1.0)	0.2; (–0.4, 0.7)	–1.4; (–2.4, –0.6)	0.3; (–0.2, 0.9)
EMS	0.3; (–0.8, 1.4)	–0.1; (–1.9, 1.9)	0.6; (0.1, 1.2)	0.4; (–0.2, 1.0)	1.6; (0.5, 2.7)	–1.5; (–2.9, –0.5)	1.2; (0.3, 2.8)
FAM	0.5; (–0.1, 1.1)	0.4; (–0.2, 0.9)	0.6; (0.2, 1.1)	0.3; (−0.2, 0.8)	0.0; (−0.5, 0.5)	−1.0; (−1.8, −0.3)	0.3; (−0.2, 0.8)
PLM	1.7; (0.8, 3.0)	−4.6; (−7.7, −2.0)	0.6; (−0.1, 1.5)	0.5; (−0.1, 1.1)	1.4; (0.3, 2.6)	−0.6; (−1.4, 0.3)	0.0; (−1.1, 0.8)
PRC	0.8; (−0.4, 2.2)	−2.5; (−5.3, −0.3)	0.5; (−0.2, 1.2)	0.4; (−0.2, 1.1)	2.2; (0.3, 4.2)	−0.8; (−1.6, 0.2)	−0.2; (−1.7, 0.8)
RRS	0.5; (−0.9, 1.8)	−4.7; (−9.3, −1.0)	0.5; (−0.2, 1.2)	0.4; (−0.3, 1.0)	2.4; (0.5, 4.7)	−0.4; (−1.3, 0.7)	−0.2; (−1.9, 0.7)
RES	2.6; (1.9, 3.6)	−0.1; (−0.9, 0.6)	0.4; (−0.1, 0.9)	0.5; (0.1, 1.0)	−0.0; (−0.6, 0.6)	−1.0; (−1.6, −0.4)	0.5; (−0.1, 1.0)
SAS	0.4; (−0.3, 1.1)	−0.7; (−1.4, −0.1)	0.5; (0.1, 1.0)	0.3; (−0.2, 0.8)	−0.6; (−1.4, 0.2)	−0.8; (−1.7, −0.1)	0.3; (−0.3, 0.8)
SUM	1.5; (0.7, 2.4)	1.2; (0.2, 2.2)	0.6; (0.2, 1.2)	0.5; (0.2, 1.2)	0.4; (−0.4, 1.2)	−1.0; (−2.0, −0.2)	0.7; (0.2, 1.4)

Table 2.

Minnow occurrence model coefficients reported on the logit scale as the mode with associated 90% highest density interval (HDI) from the posterior distribution. The intercept is interpreted as estimated occurrence probability at mean levels of covariates. Each covariate coefficient is interpreted with others held constant.

BUM, bullhead minnow; EMS, emerald shiner; FAM, fathead minnow; PLM, plains minnow; PRC, prairie chub; RRS, Red River shiner; RES, red shiner; SAS, sand shiner; SUM, Suckermouth minnow; UNDR, upstream network dam density; DL17, low flow duration variability; FH3, flood pulse count; MA2, median daily flow; RA9, variability of reversals; and TH2, annual maxima variability.

Figure 2.

Line graphs showing relationships between occurrence probability and daily streamflow magnitude (MA2, panel a), variability in reversals of flow (RA9, panel b), variability in timing of annual maximum streamflow (TH2, panel c), and upstream dam density (UNDR, panel d) for four small-bodied minnows in hydraulic response units of the upper Red River basin. The yellow lines represent emerald shiner, the bluelines represent prairie chub, the orange lines represent red shiner, and the black lines represent suckermouth minnow. Mean daily discharge was measured as m³/s.

2.2.5 Projected minnow distributions

We projected distributions adjusted for detection probability across the study area for emerald shiner, prairie chub, red shiner, and suckermouth minnow. These minnows represented varying occurrence relationships with flow regimes and fragmentation among focal species (see Section 2.2.4). A species was considered present at all HRUs with a detection. We calculated occurrence probabilities for each species at HRUs where either there were no surveys, or they were not detected using occurrence model coefficients and covariate values (Table 2). We emphasize that a high occurrence probability more appropriately represents suitable conditions for occurrence, not an assurance the species is present, and a species might be present at an HRU with a low occurrence probability. Opportunity (i.e., biogeography) and other spatial and environmental factors (i.e., biogeography) not considered here also play a role in aquatic species distributions [1, 5, 93] and, like any modeled relationship, there is inherent error. This is particularly true for HRUs closer to western ecotone with the desert-like High Plains (Figure 1). Nevertheless, the projected distributions reflect underlying ecological relationships based on where species were either detected during sampling or likely present but not detected. Further, these modeled relationships can provide an initial step to guide species reintroduction or translocation efforts, which do not depend on present occurrence state (also see Section 3.1).

All four minnows had similar distributions based on naïve occurrence (i.e., where the species was detected) and high occurrence probability along the downstream portion of the Wichita River; however, pelagophils were less likely to occur in HRUs elsewhere (Figure 3). As expected, red shiner had a high occurrence probability throughout the study area (Figure 3a). The lower red shiner occurrence probability in the southwest portion of the study area reflects the negative relationship with higher variability in flow reversals shared with all upper Red River minnows. Suckermouth minnow had a high occurrence probability along the upstream portion Red River mainstem and in the northern portion of the study area (Figure 3b). The higher suckermouth minnow occurrence probability corresponds to HRUs with more dams upstream. Suckermouth minnow occurrence probability was low in the southern portion of the study area with fewer upstream dams. Conversely, prairie chub had a low occurrence probability in the northern portion of the study area and along the upstream mainstem (Figure 3c). Prairie chub occurrence probability was higher than both emerald shiner and suckermouth minnow in the southern portion of the study area, particularly along the upper Wichita River. Although emerald shiner is more widespread overall than prairie chub, its projected distribution was narrower in the upper Red River (Figure 3d). The higher emerald shiner occurrence probability in the northern portion of the study area and along the upper mainstem is reflective of not sharing the negative upstream dam relationship with other pelagophils.

Figure 3.

Predicted species distribution maps for hydraulic response units (HRUs) in the Cross Timbers and Southwestern Tablelands level-three sub-ecoregions in the upper Red River basin (Section 2.1, Figure 1). Occurrence probabilities were estimated using modeled relationships with flow regime metrics and upstream dam density (Section 2.3). Panel a is red shiner, panel b is suckermouth minnow, panel c is prairie chub, and panel d is emerald shiner. Species were assumed to not occur at HRUs filled with black, which are completely under reservoirs.

2.3 Prairie chub relative abundance at the stream reach scale

Historically, the endemic prairie chub was abundant in the upper Red River mainstem and its higher-order tributaries [52]. Suspected population declines and poorly understood ecology has resulted in conservation concerns for prairie chub in multiple states [94, 95]. At the federal level, prairie chub is currently threatened and included as a potential endangered species on the 2021–2025 National Domestic Listing Workplan [96].

We show relationships between reach-scale adult prairie chub counts and the predictor variables longitude, stream discharge, and salinity. Prairie chub populations were surveyed at 44 stream reaches of the upper Red River basin in early autumn 2019 (Figure 1). A reach was defined as a longitudinal distance of twenty times mean wetted width constrained by a minimum of 100 m and a maximum of 500 m. Adult prairie chub counts were obtained using multi-pass removal sampling with a seine. Sampling occurred west of the overlap and hybridization zone with the morphologically similar shoal chub Macrhybopsis hyostoma [71] to minimize misidentified prairie chub (see [97] for a detailed description of sampling methods). Capture probability was fairly constant, which allowed for relative abundance comparisons among reaches. We present relationships as descriptive scatterplots of adult prairie chub counts versus each predictor variable; thus, relationships are independent and not additive. Eastern longitude and stream discharge were more highly correlated (r = 0.40) than salinity and longitude (r = 0.08) and salinity and discharge (r = −0.02).

Adult prairie chub counts of zero were most common (n = 32 reaches), with a count of only one at two reaches (Figure 4). All counts >1 were in the eastern half of the study area, with the highest counts furthest east (Figure 4a). There was a general positive linear trend in relative abundance with increasing reach discharge for counts >1 (Figure 4b). Counts >1 also increased with increasing salinity (Figure 4c). However, the highest counts were associated with intermediate salinity (∼2–6 PPT), which suggested a quadratic relationship.

Figure 4.

Scatterplots depicting relationships between reach-scale adult prairie chub counts and longitude (panel a), stream discharge (panel b), and salinity (panel c). A constant of one was added to each count prior to natural-log transforming. Reach discharge was measured as m³/s and salinity was measured as PPT.

3. Conservation implications

Our study highlights three important fish conservation aspects in a river basin: (1) the consideration of multiple spatial scales for directing conservation, (2) the tradeoffs of assemblage level (i.e., multiple species) conservation, and (3) the implications of ignoring detection error. Although beyond the scope here, the minnow occurrence relationships can also be used for predictive simulations under different flow regime and fragmentation scenarios. For example, changes in species distributions could be predicted under different levels of dam removal or long-term changes in flow magnitude with increased drought (or both).

3.1 Multiscale fish conservation strategies

The occurrence and relative abundance relationships for prairie chub can be used to identify target areas in the stream network with a higher chance of habitat restoration or reintroduction success. Spatial position has been shown to be strongly associated with the structure of fish populations and assemblages [5, 98, 99]. Prairie chub’s distribution in the upper Red River basin is severely constrained by upstream dams (Figure 2d). There is essentially no probability of prairie chub occurring below heavily dammed HRUs, presumably due to connectivity requirements for pelagophil reproduction (see Section 1.3). Thus, finer-scale conservation actions (e.g., reintroduction or instream habitat enhancements) in these HRUs would be futile and waste available resources. The most favorable HRUs in the upper Red River stream network for prairie chub occurrence are along the mainstem or higher-order tributaries (i.e., higher long-term flow magnitude), with low upstream dam density and more constancy in rate of change. In particular, HRUs along the Wichita River without prairie chub detections had a high occurrence probability (Figures 1 and 3c). If feasible, increasing flow in HRUs with lower dam density could increase the range of the favorable area. Spatial position was also associated with reach-scale prairie chub relative abundance. All high adult counts were associated with reaches in the eastern portion of the study area and higher discharge (Figure 4a and c). Longitude and discharge were somewhat confounded, and we did not consider the effect of each with the other held constant (see Section 2.3). The number of low-flow days increases further east in the upper Red River basin [75]. However, stream discharge was variable at both the HRU and reach scale. Thus, management actions targeting prairie chub would likely be most effective at reaches in high occurrence probability HRUs in the eastern portion of the study area with higher average flow magnitude. Higher prairie chub relative abundance was also associated with intermediate salinity levels (Figure 4c). Salinity has been shown to be strongly associated with fish assemblage structure in Great Plains streams [100], and a quadratic relationship with population size makes sense ecologically for a freshwater species adapted to semiarid streams. However, salinity is highly variable across both space and time in the upper Red River basin [101]. There is also a salinity gradient at the ecotone with the Cross Timbers that constrains prairie chub’s eastern distribution and forms a hybrid zone with shoal chub Macrhybopsis hyostoma [71]. Thus, to effectively implement a salinity target for prairie chub conservation, improved salinity monitoring would be needed. The findings for our study period are likely reflective of multiscale prairie chub ecological relationships in both wet and dry periods. Occurrence probability among HRUs has been shown to be similar in both wet and dry periods [56], and reach-scale adult counts were collected during a relatively wet period [97].

3.2 Assemblage-level fish conservation

The mixed occurrence relationships among minnows with flow regime and fragmentation have implications for upper Red River basin conservation strategies. It is important that managers consider conservation actions that benefit target species without detrimental effects to other native fishes in the assemblage. Reducing variability in annual maxima timing and removing upstream dams might be beneficial for prairie chub. However, emerald shiner and suckermouth minnow have high occurrence probabilities in the northern portion of the study area where these changes might take place. Unless the mechanism(s) driving the distributions of emerald shiner and suckermouth minnow is better understood, it is possible for conservation actions designed to improve conditions for prairie chub to incidentally harm other species. Fragmentation might prevent prairie chub from successful upstream movement that is important for completing the pelagophil life cycle [58, 102]. However, emerald shiner and suckermouth minnow are more widespread and might be able to adapt to a wider variety of conditions including fragmented river systems (e.g., phenotypic plasticity, [103, 104]). Although prairie chub and emerald shiner are both pelagophils, there is evidence that emerald shiner is less sensitive to flow disturbances than some other pelagophil species [105] and portions of some pelagophil populations are residents that do not make upstream spawning movements [106, 107]. Emerald shiner is also adapted to lentic environments [108] and may benefit competitively in habitats near reservoirs. Because suckermouth minnow is not a pelagophil species, it might not require long unimpeded lengths of river for spawning. It is also possible that another unmeasured or confounding habitat metric is the driver of the emerald shiner and suckermouth minnow distributions. It is prudent to balance conservation efforts to benefit target species while maintaining habitat for other natives. For example, strategic dam removal in the northern portion of the upper Red River basin could benefit prairie chub while preserving habitat favorable to emerald shiner and suckermouth minnow. Also, consideration of flow-regime patterns that benefit numerous minnows (e.g., increased flow constancy) or may not affect other species (e.g., increasing flow magnitude) provide a balanced assemblage-level conservation approach.

3.3 Imperfect and variable detection

Species occurrence is never perfectly observed (i.e., detection is imperfect), and detection probability varies differently among species and sampling methods across sampling conditions [86, 109, 110]. Thus, ignoring detection error results in only naïve occurrence and apparent species distributions (i.e., true distributions are always larger than observed). In addition to underestimating true occurrence, a high detection probability can be misinterpreted as high habitat suitability and lead to misinformed conservation strategies [110, 111, 112]. Species-specific fish detection probability varies in relation to numerous environmental characteristics (e.g., water depth, water temperature, conductivity, water clarity, and flow [113, 114]). We show ignoring detection probability at the HRU scale resulted in similar apparent distributions for all minnows in the upper Red River basin (Figure 3). There is, of course, uncertainty in the predicted distributions (also see Section 2.2.50). Nevertheless, adjusting for detection error resulted in a clearer picture of true distributions and underlying ecological relationships. Detection probability has also been shown to vary among species and across sampling conditions at finer scales in Great Plains streams [115, 116]. Given the highly variable nature of the Great Plains stream environment, accounting for detection error in fish species distribution studies at all spatial scales is particularly important for sound river basin conservation.

Acknowledgments

This research is a contribution of the Oklahoma Cooperative Fish and Wildlife Research Unit (U.S. Geological Survey, Oklahoma Department of Wildlife Conservation, Oklahoma State University, and Wildlife Management Institute cooperating). Funding was provided by the Oklahoma Department of Wildlife Conservation (F18AF00575) and the U.S. Geological Survey, Water Smart Program. This study was performed under the auspices of Oklahoma State University protocol # AG-19-11. We thank the numerous scientists for their sampling and data collection efforts that made this study possible. We also thank Glena Osban for technical assistance. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government. The Red River fish survey dataset compiled by the authors is available at https://doi.org/10.22488/okstate.22.000006.

Data source
Fishes of Texas (www.fishesoftexas.org/home/)
iDigBio (www.idigbio.org/)
MARIS (www.sciencebase.gov/catalog/item/51c45ef1e4b03c77dce65a84)
Oklahoma Conservation Commission (www.ok.gov/conservation/)1
Oklahoma Museum of Natural History (www.samnoblemuseum.ou.edu/)
Oklahoma Water Resources Board (www.owrb.ok.gov/)2
Perkin et al. [75]3
Texas State University4
VertNet (www.vertnet.org/index.html)

Table A1.

Data sources for stream-fish assemblage surveys compiled from 1983 to 2015.

¹

Cheryl Cheadle ([e-mail]) and Jason Ramming ([e-mail]).

²

Chris Adams ([e-mail]).

³

DOI:10.1890/14-0121.1.

⁴

David Ruppel ([e-mail]) and Tim Bonner ([e-mail]).

Footnotes denote contact(s) or source for datasets not available online.

ID	Collector	Proportion of surveys
1	Oklahoma Conservation Commission	0.55
2	Oklahoma Water Resources Board	0.15
3	Oklahoma Department of Environmental Quality	0.08
4	Perkin et al.	0.06
5	Texas Tech University	0.08
6	Miscellaneous	0.08

Table A2.

Collector descriptions used for the grouping factor in the detection model and the proportion of surveys.

Surveys were pooled by predominant collectors (i.e., ≥ 10 surveys, also see Table A1).

Species	Intercept	HRU area	Drainage area
Mean	−0.1 (−0.7, 0.5)	−0.3 (−0.5, −0.1)	0.1 (-0.4, 0.5)
BUM	0.9 (0.2, 1.6)	−0.3 (−0.6, −0.1)	−0.2 (−0.5, 0.1)
EMS	−0.8 (−1.5, −0.1)	−0.4 (−0.6, −0.1)	0.6 (0.2, 1.0)
FAM	−0.1 (−0.8, 0.6)	−0.2 (−0.5, 0.1)	−0.3 (−0.7, 0.1)
PLM	−0.0 (−0.7, 0.7)	−0.3 (−0.5, −0.1)	0.4 (0.1, 0.7)
PRC	−1.2 (−2.0, −0.4)	−0.4 (−0.7, −0.2)	1.0 (0.6, 1.4)
RRS	−0.8 (−1.6, −0.1)	−0.3 −0.6, −0.2)	0.5 (0.2, 0.8)
RES	1.9 (1.0, 2.7)	−0.3 (−0.6, 0.1)	−0.1 (−0.4, 0.3)
SAS	−0.5 (−1.3, 0.2)	−0.1 (−0.4, 0.3)	−0.4 (−0.9, −0.1)
SUM	−0.2 (−0.9, 0.4)	−0.2 (−0.4, 0.1)	−0.7 (−1.0, −0.3)

Table A3.

Minnow detection model coefficients reported on the logit scale as the mode with an associated 90% highest density interval from the posterior distribution. The intercept is interpreted as estimated detection probability at mean levels of covariates. Each Covariate coefficient is interpreted with other covariates held constant.

HRU, hydraulic response unit; BUM, bullhead minnow; EMS, emerald shiner; FAM, fathead minnow; PLM, plains minnow; PRC, prairie chub; RRS, Red River shiner; RES, red shiner; SAS, sand shiner; SUM, Suckermouth minnow; UNDR, upstream network dam density; DL17, low flow duration variability; FH3, flood pulse count; MA2, median daily flow; RA9, variability of reversals; and TH2, annual maxima variability.

We fit the occupancy model with the program JAGS [117] called from the statistical software R [118] using the package jagsUI [119]. We used vague truncated normal priors for species coefficients and vague gamma priors for associated standard deviations [64]. Posterior distributions for coefficients were estimated with Markov chain Monte Carlo methods using four chains of 20,000 iterations each run in parallel after a 5,000-iteration burn-in phase (thinning = 10). We considered adequate convergence a potential scale reduction factor ($\widehat{R}$) <1.05 [120] and “grassy” trace plots for all parameters [64]. We calculated the 90% highest density intervals using the R package HDInterval [121].

We examined model fit using a posterior predictive check [64] based on the goodness-of-fit test described by [122]. We simulated expected species encounter histories under model parameters to compare discrepancies with observed encounter histories and calculated a Bayesian p-value (0.47). A Bayesian p-value near 0.5 suggests adequate fit and extreme values (i.e., <0.05 or >0.95) indicate a lack of fit [64, 123].