Characteristics of the natural flow regime paradigm explain occurrence of imperiled Great Plains fishes

The natural flow regime (NFR) paradigm posits that naturally occurring temporal fluctuations in streamflow are necessary for maintaining natural ecological communities. Conservation actions guided by the paradigm have contributed to stream organism conservation on a global scale, yet NFR applications in highly altered Great Plains streams are lacking. We analyzed historical (1980–2017) fish occurrence and flow data from sixteen gage locations across three Great Plains river basins with the goal of relating flow indices calculated for one year prior to each fish collection with the occurrence of fishes belonging to the highly imperiled pelagicbroadcast spawning (PBS) reproductive guild. We fit random forest models using flow indices, gage identification, and watershed area as predictor variables and PBS fish occurrence as the response variable for all suspected or confirmed PBS fishes in each basin. Results revealed that flow indices from all NFR characteristics were useful for predicting PBS fish occurrence, but the identity of the most important flow indices varied by species and river basin. We also found that gage location was an important predictor variable, indicating flow– ecology relationships are spatially contingent and suggesting the influence of other factors that vary spatially (e.g., river fragmentation and water quality). Partial dependence plots identified thresholds in flow indices such as the timing of maximum flows and fall rate that may be associated with the presence or absence of PBS fishes, and these plots can be used to identify flow index target values to benefit conservation efforts. Time series predictions for likelihood of occurrence at a single gage in the Brazos River basin revealed a positive correlation between likelihood of occurrence and non-drought years for three of four species present, supporting recent findings that drought suppresses PBS fish populations. This study provides insight into the flow requirements of some of the most imperiled stream fishes in North America and contributes to environmental flow science on a global scale by establishing a robustmethodology for delineating flow–ecology relationships.


INTRODUCTION
Flow regimes of rivers and streams are altered on a global scale, especially in the USA (Carlisle et al. 2019). Flow regime alteration often imperils freshwater fishes and other riverine species adapted to naturally occurring dynamics in flow (Bunn and Arthington 2002, Lytle and Poff 2004, Lehner et al. 2011, Grill et al. 2015. The notion that flow is the master variable (Power et al. 1995, Boltz et al. 2019) maintaining natural ecological communities, geomorphic functions, water quality, and habitat connectivity is at the center of the natural flow regime (NFR) paradigm (Poff et al. 1997, Annear et al. 2004). The NFR paradigm predicts that reestablishment of natural flow dynamics in altered systems should support recovery of imperiled populations and native species assemblages, and this notion is now a central concept to guiding environmental flow management and river restoration worldwide (Anderson et al. 2006, Jowett and Biggs 2006, Suen and Eheart 2006, Assani et al. 2011, Opperman et al. 2019. Implementation of the NFR paradigm as a conservation tool first requires insight into departures from natural flow conditions. The NFR paradigm identifies five hydrologic characteristics that are useful for describing streamflow, including magnitude, frequency, duration, timing, and rate of change (Poff et al. 1997;Fig. 1). These, in turn, are quantified based on specific hydrologic indices that can be used to measure alteration to NFRs. For example, the Indicators of Hydrologic Alteration (IHA) tool quantifies natural and altered flow regimes using 31 indices (described in detail by Richter et al. 1996). Additional indices have been incorporated in regional hydrologic evaluation tools, such as the Hydrology-based Environmental Flow Regime (HEFR) tool developed for rivers in Texas (Opdyke et al. 2014) as well as the Great Plains Environmental Flows Information Toolkit (EFIT; Valente et al. 2019). Metrics based on the NFR paradigm are used to quantify flow alterations and guide flow regime restoration on a global scale (Mathews and Richter 2007), yet our understanding of the ecological consequences of flow alteration and restoration requires continued research (Davies et al. 2014). Fig. 1. Conceptual diagram illustrating flow characteristics from the natural flow regime (NFR) paradigm (duration, frequency, magnitude, rate of change, and timing) for an example unmodified historical year (gray area) and an example altered contemporary year (bold black line). Portions of each hydrograph are colored by flow characteristic to illustrate changes between historical (dashed) and contemporary (solid) periods. The data illustrated are from the Canadian River at Canadian, TX (U.S. Geological Survey gage 07228000; C04 in Fig. 2) for the years 1945 (historical) and 1985 (contemporary) to illustrate changes associated with the completion of Stanford Dam and creation of Lake Meredith reservoir in 1965. See Bonner and Wilde (2000) and Costigan and Daniels (2012) for additional details regarding long-term changes in fish assemblages and flow at this location.
Flow-ecology relationships based on the NFR paradigm can be used to guide management of flows and direct water transactions that aim to restore the ecological integrity of riverine ecosystems (Horne et al. 2017, Arthington et al. 2018, Wheeler et al. 2018a. For example, in the flow assessment known as the Ecological Limits of Hydrologic Alteration (ELOHA), users define and quantify flow-ecology relationships when designing flow restoration guidelines ). However, defining flow-ecology relationships can be difficult given the large number of available ecological (Poff and Zimmerman 2010) and flow Poff 2003, Arthington et al. 2006) indices, many of which are inter-correlated (Gao et al. 2009, Kennard et al. 2010. Another challenge comes with selection of flow indices prior to analyses (e.g., characteristics related only to frequency, Perkin and Bonner 2011) instead of relying on data-driven approaches capable of uncovering relationships involving multiple NFR characteristics. Finally, flow-ecology relationships are not always linear, but instead may be nonlinear or exhibit threshold responses that confound attempts to integrate flow-ecology relationships into models for conservation and management (Rosenfeld 2017). However, new statistical tools such as tree-based machine learning algorithms are well suited for analyzing large numbers of predictor variables and can incorporate nonlinear or threshold relationships (Liaw and Wiener 2002, Elith et al. 2008, Barker et al. 2014, Hain et al. 2018. Such adaptable models provide potential for establishing flow-ecology relationships for regions where flow and ecological data exist but previous challenges have slowed progress toward implementation of the NFR paradigm. The Great Plains region of North America contains some of the most altered rivers and fish assemblages on Earth (Lehner et al. 2011, Costigan andDaniels 2012). Reservoir construction and operation, channelization and levee installation, groundwater overexploitation, and widespread land use transformation for agriculture have collectively altered most Great Plains Rivers (Gido et al. 2010, Perkin et al. 2015a. Consequently, native fishes have declined throughout the region (Hoagstrom et al. 2011), particularly species belonging to the pelagicbroadcast spawning (PBS) guild characterized by release of semi-buoyant ova that require downstream drift during development within pelagic zones of larger rivers (Perkin and Gido 2011, Hoagstrom and Turner 2015, Worthington et al. 2018. The flow requirements for some PBS fishes are well studied, and previous works suggest flow magnitude during summer spawning months is critical for population persistence (Bonner and Wilde 2000, Perkin et al. 2015a). However, the importance of other months of the year or the broader suite of hydrologic indices has not been rigorously quantified to the extent that flow-ecology relationships exist. A recent study on the Brazos River, Texas, USA (Rodger et al. 2016), assessed the relationship between flow magnitude and recruitment for a single species, Shoal Chub Macrhybopsis hyostoma, and results from that study were used to develop flow recommendations (Texas Instream Flow Program 2018). However, such recommendations do not exist for the majority of river segments or fish species in the Great Plains, aside from a few notable exceptions , Durham and Wilde 2009, Wilde and Urbanczyk 2013. Consequently, an improved understanding of flow-ecology relationships is necessary for the conservation and restoration of fish assemblages in rivers of the southern Great Plains (Worthington et al. 2018).
The goal of this study was to define flowecology relationships for PBS fishes in three southern Great Plains river basins. We used a data-driven methodology, rather than relying solely on professional judgment (Poff and Zimmerman 2010), to begin the process of flow index selection. We hypothesized that flow indices describing flow magnitude, particularly those pertinent to flows during the PBS reproductive period of May to September, would be most influential in predicting the occurrence of PBS fishes compared with flow indices describing the other characteristics of the NFR: frequency, duration, timing, or rate of change. We tested this hypothesis using species-specific random forest models fit to a suite of flow indices and species occurrence data to compare the relative importance of NFR characteristics in describing species occurrence. v www.esajournals.org

Study area
The southern Great Plains is a semi-arid region of North America with tallgrass prairie to the east and shortgrass prairie further west (Kunkel et al. 2013). In an unaltered state, the streams and rivers of the southern Great Plains were dynamic and flashy, characterized by strong floods and seasonal stream drying (Dodds et al. 2004; Fig. 1). Native fish assemblages of the region are well adapted to these dynamic flows, including opportunistic life history strategies such as rapid maturation as well as the ability to recolonize habitats that undergo drying and rewetting (Lytle andPoff 2004, Mims et al. 2010). During the past five decades however, flow regimes have been disrupted by diversion of water for agricultural and municipal uses by rapidly growing human populations (Gido et al. 2010, Perkin et al. 2015a Sotola et al. 2019) river basins, with distributions and abundances declining for each species in each basin, in part due to an inability to recolonize former habitat (Perkin et al. 2015b, see review in Worthington et al. 2018).

Streamflow data
Daily streamflow data were downloaded for 16 United States Geological Survey (USGS) gages (USGS 2019; Fig. 2) and used to quantify hydrological indices for the period 1979-2017. For the 68 annual IHA indices (Richter et al. 1996), we defined a water year as October 1 to September 30. We calculated 35 additional indices specific to the region following the HEFR framework described by Opdyke et al. (2014). Of these 103 flow indices, those missing more than 25% of the values across all gages were excluded, resulting in 63 flow indices (i.e., 34 from IHA, 29 from HEFR; Appendix S1: Table S1). We addressed multicollinearity within the reduced set of flow indices by calculating Pearson correlation coefficients using the function cor from the R package stats and removed those with a coefficient greater than |0.70| (Appendix S1: Fig. S1). This resulted in retention of eighteen flow indices (i.e., 9 from IHA, 9 from HEFR) as inputs for our models in addition to stream gage (as a factorial variable) and drainage area (Table 1). Flow indices measured from one year prior to fish collections (see next section) were retained for further analysis because these flows are most linked with the short-lived species in the collections (Worthington et al. 2018). Flows from the same year as fish collections were not used because index calculations included data measured after fish collections took place (e.g., a fish collection taking place in May would not be impacted by flows measured later in the same water year). We further transformed flow magnitude indices by dividing by drainage area (km 2 ) to standardize flow data among gages and then z-scored transformed all indices prior to analysis (Kennard et al. 2010).

Fish occurrence data
Fish occurrence (i.e., present and detected) data, predominantly from seining collections, included 105 Great Plains fish species from 222 collections between 1980 and 2018 (Appendix S2: Table S1). Each fish collection was assigned to a USGS gage when a gage was within a distance of 25 river km; otherwise, collections were not included in analysis (Fig. 2). We retained streamflow gages and their associated fish assemblage data for analysis if at least three fish collections were linked to a gage, resulting in four gages in the Brazos, five in the Canadian, and seven in the Red (Table 2). All fish species were classified by their reproductive strategy as either a pelagic, benthic, or live-bearing species based on the Fish-Traits database (Frimpong and Angermeier 2010) or additional information for species not included in that database (Balon 1975, Perkin et al. 2015a, Worthington et al. 2018. Of these, we focused on nine species classified by Worthington et al. (2018) as known or suspected pelagic-broadcast spawning (PBS) fishes and analyzed species independently by basin using presence (class = 1) and absence (class = 0) data from each location across multiple samples. The guild of PBS fishes was chosen because of their responsiveness to flow as an environmental cue v www.esajournals.org for reproduction wherein the species release semi-bouyant ova to be passively carried downstream (reviewed in Balon 1975, Hoagstrom and Turner 2015, Worthington et al. 2018). These filtering steps resulted in 11 PBS species-specific datasets across the three basins (Table 3).

Statistical analyses
We used random forest (RF) models to assess flow-ecology relationships (Liaw and Wiener 2002). We fit models to each species within each basin independently using annual indices, gage ID as a factor, and drainage area. We included the gage ID factor to represent unexplained variance among sites within a basin, and drainage area to represent differences among sites along the river continuum because flow-ecology relationships are known to vary across space (Bruckerhoff et al. 2019). These variables thus allowed for assessing gage-specific occurrences of PBS fishes. Because PBS species are sometimes rare, we utilized a synthetic minority oversampling technique (SMOTE) to address class imbalance (Chawla et al. 2002). We used the SMOTE function from the DMwR package in R to generate equivalent numbers of absences and presences by doubling the lesser of the two categories while under-sampling a subset of the other class. For example, if a basin dataset had 10 presence and 40 absence records for a given species, we would double the lesser (i.e., double the 10 presences) and subsample the greater to match (subset from  Table 2. v www.esajournals.org 40 to 20 absence records) so that there were equal absence and presence data points (i.e., 20 each class). All RF models were fit to SMOTEadjusted datasets using the randomForest function from the randomForest package with 500 trees, 4 variables tried at each split, and crossvalidated using the rf.crossValidation function from the rfUtilities package in R. This latter function allowed for fivefold cross-validation such that model performance could be assessed across training and testing subsets of the data. We assessed model performance by considering accuracy, kappa, specificity, and sensitivity across five subsets of the data. Model accuracy is the average between specificity (i.e., the proportion of absences that are correctly predicted) and sensitivity (i.e., the proportion of presences that are correctly predicted). Kappa is the ratio of the difference between observed and expected accuracy over 1 minus the expected accuracy and ranges 0-1. Kappa values ranging 0.41-0.60 are considered indicative of moderate classification performance, 0.61-0.80 are indicative of good classification performance, and values 0.81-1.0 are indicative of very good classification performance (Landis and Koch 1977). We tested the significance of each RF model using the rf.significance function from the rfUtilities package in R.
Variable importance and fish responses for each model were assessed using the relative influence index and partial dependence plots. The relative influence of a variable for each of these species-specific models is measured by its mean decrease in Gini (MDG). Gini impurity is a measure of the likelihood of a particular observation being incorrectly classified based on the distribution of classes in a data subset. When Gini impurity reaches zero, it suggests that all the v www.esajournals.org observations in a group have been correctly classified together. In other words, MDG is a proxy for variable importance such that the greater the MDG, the more important the variable. We illustrated relative variable importance from each species-specific model by dividing each MDG value by its maximum and plotting a heat map of all relative importance values across all species Note: Basin ID labels correspond to locations shown in Fig. 2   and river basins. This approach allowed for visually assessing trends in variable importance across multiple taxa. For each NFR flow characteristic class, we plotted partial dependence plots (PDPs) to illustrate variable-specific relationships with likelihood of occurrence (i.e., continuous range 0-1 from the RF model fit). The benefit of PDPs was that individual indices could be explored while holding all other indices at their mean, that is, the partial influence of each index on likelihood of occurrence could be explored. We used the rf.partial.prob function from the rfUtilities package in R to create PDPs that illustrate smoothed LOESS regression lines fit to splines as described in detail by Baruch-Mordo et al. (2013). We combined PDPs for all species within a basin to illustrate consistencies and inconsistences in species-specific likelihood-ofoccurrence profiles for each selected flow metric. We used untransformed flow metric data to construct PDPs to improve interpretation. We tested our hypothesis that magnitude indices were the most important flow characteristic for predicting PBS fish distributions using outputs from RF models. We calculated relative importance from MDG values for each flow index included in each model and combined all species within each basin during testing. For each basin, we then compared the base mean (i.e., mean across all flow indices) with distributions within each flow metric using the nonparametric Wilcoxon signed rank test with pvalues adjusted according to the method of Holm (1979) and implemented with the stat_com-pare_means function from the ggpubr package. Finally, we used the predict.randomforest function to estimate likelihood-of-occurrence values for all years (sampled and unsampled) at the Salt Fork of the Brazos River near Aspermont, Texas, USA (USGS gage 08082000). We plotted the time series of likelihood of occurrence during 1980-2017 and placed points on the y-axis at 0 when a species was not collected during sampling and 1 when a species was collected during sampling to illustrate temporal changes in likelihood of occurrence and observed species occurrence. We then regressed likelihood-of-occurrence values (dependent) against Palmer drought severity index for March-September as described by Perkin et al. (2015b) to assess the relationship between likelihood of occurrence and water availability determined by drought. We then tested the significance of the regression line slope for each species independently and calculated the coefficient of determination. All analyses were conducted in R version 4.0.3 (R Development Core Team 2020).

RESULTS
The number of fish collections retained for modeling included 130 collections from four gages in the Brazos River Basin, 30 collections from five gages in the Canadian River Basin, and 59 collections from seven gages in the Red River Basin ( Table 2). The number of collections per gage ranged 14-56 for the Brazos River Basin, 3-9 for the Canadian River Basin, and 3-27 for the Red River Basin. Models were fit for four species in the Brazos, three species in the Canadian, and four species in the Red (Table 3). Within the Red River Basin, we combined Shoal Chub and Prairie Chub to form the Macrhybopsis species complex because the species were not yet split during the time of collections (Eisenhour 1997(Eisenhour , 1999. Presence records for individual species ranged from 45-83 for the Brazos River basin, 10-19 for the Canadian River basin, and 23-48 for the Red River basin. Following SMOTE adjustment, balanced occurrence ranged 90-114 for the Brazos River basin, 20-28 for the Canadian River basin, and 28-48 for the Red River basin. All of the eleven RF models were significant (Table 3).
Model accuracy values ranged 0.74-0.88, 0.78-1.00, and 0.72-1.00 for the Brazos, Canadian, and Red River basins, respectively (Table 3). Kappa values were within the range of good classification performance (i.e., 0.61-0.80) for all models except Red River Shiner (Kappa = 0.44; moderate classification performance) and Sharpnose Shiner (Kappa = 0.81; very good classification performance). The most important indices based on relative MDG differed among species and basins without any uniformly important indices across all species and basins (Fig. 3). There was no support for our hypothesis that magnitude indices alone were the most important indices for predicting PBS fish occurrence (Fig. 4). In the Brazos River, gage ID and drainage area had significantly higher relative importance compared with the base mean, whereas one pulse per fall (i.e., v www.esajournals.org one large flow pulse-taking place during September to November), one pulse per summer (i.e., one large flow pulse-taking place during June to August), large flood frequency, and zero flow days during winter had significantly lower relative importance. In the Canadian River, gage ID and date of maximum had significantly higher relative importance, while zero flow days in fall, spring, and winter had significantly lower relative importance values compared with the base mean. In the Red River basin, gage ID, one pulse per winter, large flood frequency, and zero flow days during winter had significantly higher relative importance, while 30-day maximum and zero flow days during spring and summer had significantly lower relative importance compared with the base mean.
Partial dependent plots for select flow indices within the NFR characteristic classes revealed some consistencies and idiosyncrasies among Fig. 3. Heat map illustrating relative importance of NFR paradigm flow characteristics and geographic variables (rows) from random forest models fit to fishes (columns) where darker blue represents higher relative importance and orange represents lower relative importance. v www.esajournals.org 9 September 2021 v Volume 12(9) v Article e03669 The horizontal dashed line in each plot represents the mean relative importance value, and asterisks denote significantly higher relative importance (boxes above dashed line) or significantly lower relative importance (boxes below dashed line) compared with the mean for the basin ( Ã P < 0.05, ÃÃ P < 0.01). species within and among basins. In the Brazos River, species-specific responses to flow magnitude during spring (Fig. 5a) and high flow frequency (Fig. 5b) were consistent, while responses to 30-day maximum flows (Fig. 5c), date of maximum (Fig. 5d), and fall rate (Fig. 5e) were most consistent among Plains Minnow and Sharpnose Shiner. In the Canadian River, species-specific likelihood-of-occurrence values were inconsistent for spring magnitude (Fig. 5f), 30-d maximum (Fig. 5h), and fall rate (Fig. 5j), but much more consistent for high flow frequency (Fig. 5g) and date of maximum (Fig. 5i). In the Red River, Red River Shiner was more responsive to spring magnitude (Fig. 5k) and high flow frequency (Fig. 5l), Chub Shiner was more responsive to 30-d maximum (Fig. 5m), Red River Shiner and Plains Minnow were more responsive to date of maximum (Fig. 5n), and all species responded consistently to slower fall rates (Fig. 5o).
Predicted likelihood-of-occurrence values for the Salt Fork Brazos River near Aspermont, Texas, USA, were temporally dynamic. High likelihood-of-occurrence years typically correlated with observed occurrences during sampling, and low likelihood-of-occurrence years typically correlated with observed absences during sampling. This was most pronounced during the 2011-2012 drought period when a consistent effort was afforded but species did not occur (Fig. 6). Correlations between Palmer drought severity index and likelihood-of-occurrence values revealed significant positive correlations Fig. 5. Partial dependence plots for select flow indices from each NFP flow characteristic groups based on random forest models fit to fish species from the Brazos River (top row), Canadian River (middle row), and Red River (bottom row) basins. Panels show partial (i.e., all other variables held at their mean) flow indices (x-axis) versus likelihood of occurrence (y-axis).

DISCUSSION
Our findings suggest that occurrence of PBS fishes in the southern Great Plains is related to multiple NFR characteristics beyond simply magnitude. Despite magnitude often being identified as the most important flow characteristic (reviewed by Poff and Zimmerman 2010) as well as the known importance of spring flow magnitude for PBS reproduction (reviewed by Hoagstrom and Turner 2015), we found that the contributions of magnitude in explaining PBS fish occurrence were generally comparable to other NFR characteristics such as rate of change and timing. We also found that stream gage ID and drainage area were consistent predictors of PBS fish occurrence, confirming that the spatial context of a river basin must be considered in flow-ecology relationships (Arthington et al. 2018, Bruckerhoff et al. 2019). The consistently high importance of stream gage ID indicates the potential influence of other environmental factors not included in this analysis, such as river fragmentation, water quality, and other landscape variables (Perkin et al. 2015a). Although we found evidence of emergent patterns across taxa and basins, we also documented more nuanced flow-ecology relationships that preclude implementing any one strategy across all PBS fishes (Hoagstrom et al. 2011, McKenna et al. 2018. Although theory and our empirical work suggest restoring all characteristics of the natural flow regime will benefit the maximum number of species (Poff et al. 1997), our work provides guidance for which characteristics of the natural flow regime might be prioritized for instances in which not all flow characteristics can be restored at once (e.g., Propst and Gido 2004) or where nonstationarity among hydrologic regimes prevents the restoration of some flow characteristics (Poff 2018). For example, prioritizing flow management during periods of drought may serve as a useful target for environmental water transaction strategies that seek to restore or preserve flow regimes during a time when Great Plains fishes are vulnerable to water stress (Perkin et al. 2019).
Our initial hypothesis that flow magnitude would be the single most influential NFR class was not supported. Our hypothesis was based on magnitude being the focus of previous Great Plains flow-ecology research (e.g., Wilde 2009) and because declining streamflow magnitude is a known driver of Great Plains stream fish declines (Hoagstrom et al. 2011). In reality, flow regimes are more complex than magnitude alone, even though reduced magnitude and subsequent drying are the most visually obvious changes to Great Plains streams. In our case study of the Salt Fork Brazos River, likelihood of occurrence for fishes declined during drought, a period typically associated with lower discharge magnitudes (Nielsen-Gammon 2012). However, our models predicted reduced likelihood of occurrence for PBS fishes using a variety of flow metrics beyond magnitude, namely frequency, timing, and rate of change. More broadly, we found that fall rate was consistently among the important indices across species, and partial dependence plots consistently illustrated greater likelihood of occurrence at intermediate but relatively slow fall rates. In their review of hydrologic alterations to Great Plains Rivers, Costigan and Daniels (2012) noted that faster fall rates were among the most dramatic changes to natural flow regimes across the region. Increased fall rates could be particularly detrimental to PBS fishes that require obligatory drift periods during early development after synchronous spawning during flow pulses . We also found that the timing of flow events was consistently important. For example, in the Red River all species exhibited intermediate optima in likelihood of occurrence as a function of date of maximum, albeit some species (Plains Minnow and Red River Shiner) were more responsive. Flow magnitude optima related to recruitment of single species were previously explored , Durham and Wilde 2009, Rodger et al. 2016, and our approach extends the breadth of species and flow indices for which preliminary flowecology relationships exist. The presence of intermediate optima suggests that ecological flow v www.esajournals.org requirements do not follow a pattern of more water is always better, but rather are bound by upper and lower limits for multiple flow indices (Rodger et al. 2016).
Despite differences in methodologies, our results are consistent with the limited number of existing flow-ecology analyses conducted in the Great Plains. Wheeler et al. (2018a) described two complementary classes of methodologies used to establish flow-ecology relationships, including rate-based and state-based approaches. Ratebased approaches generally assess changes in demographic properties through time, such as dispersal, recruitment, or survival; whereas, state-based approaches assess snapshot conditions that represent the state of a system, such as species abundance, richness, or occurrence. These two approaches represent endpoints on a continuum, between which hybrid strategies involving both methods exist (Wheeler et al. 2018a). Our results suggest a repeated-states approach in which the state of species occurrence (presence or absences) is assessed repeatedly through time is capable of producing flow-ecology relationships comparable to purely ratebased approaches. For example, rate-based studies suggest Peppered Chub require summer flows averaging 11.9 m 3 /s or more to maintain positive population growth in the Canadian River , Smalleye Shiner require an average of 6.4 m 3 /s during the reproductive season to maintain positive population growth in the upper Brazos River (Durham and Wilde 2009), and Shoal Chub recruitment is greatest when summer pulses ranging between 79 and 139 m 3 /s in the lower Brazos River occur (Rodger et al. 2016). We found likelihood of occurrence for Peppered Chub in the Canadian River based on 30-d maximum flows increased linearly before reaching a plateau, supporting previously documented declines in likelihood of occurrence as flows are depleted Durham 2008, Perkin et al. 2019). Furthermore, we found that likelihood of occurrence for Smalleye Shiner in the Brazos River was significantly negatively correlated with drought intensity, supporting previous assessments that minimal flows are required for persistence  and that populations are sensitive to extirpation when minimum summer flows are not maintained (Mayes et al. 2019). Such consistencies among these patterns despite the contrasting methodologies employed suggest flow-ecology relationships derived from ratebased, state-based, and repeated-states assessments are indeed complementary (Wheeler et al. 2018a).
Caveats and limitations to our study approach should be considered. First, our study relied on only contemporary flow indices spanning the period after the most significant changes to Great Plains streamflow regimes occurred (post-1980;Perkin et al. 2015a, Yarnell et al. 2020. We acknowledge this means basin-wide historical flow alterations are not addressed in our work (Eng et al. 2017), but we point out restoring flows to their historical levels may no longer be possible in some streams (Poff 2018). By modeling contemporary fish data collected solely in a postalteration context, concerns regarding decoupled flow-ecology relationships were alleviated as fishes were responding to the existing altered flow conditions (Patrick andYuan 2017, Hain et al. 2018). Second, limited data points at certain gages may change the shape of flow-ecology relationships for those locations and it may be necessary to evaluate the ecological relevance of some indices. For example, Arkansas River Shiner and Plains Minnow models from the Canadian River suggested the highest likelihood of occurrence at minimum 30-day maximum flows. This may be more of an artifact of the flow characteristics of the years for which biological data are available, especially given documented extirpations of these same species during extreme low summer flows elsewhere in their range (Perkin et al. 2015b). A third caveat is the general criticism of states approaches that rely on temporally restricted (i.e., snapshot) data and therefore may result in patterns that arise from spatiotemporal correlation and not due to a mechanistic flow-ecology relationship or process (Wheeler et al. 2018b). We removed flow indices with high multicollinearity prior to analysis and used RF models in which a random subset of predictor variables was selected during tree construction; therefore, collinearity was less of a concern (in a statistical sense) compared with other parametric models such as multiple regression (Strobl et al. 2009). The fourth consideration of our study design was the use of species-specific models over a single model with species terms v www.esajournals.org that allow for assessing interactions (Fujiwara et al. 2019). Our use of single-species models was largely out of necessity, because PBS fishes can be rare at times and occurrence data are therefore zero-inflated. Such imbalance between presence and absence classes typically existed, and we used species-specific SMOTE-adjusted data to address imbalance. We also point out that biological interactions are generally less intense, compared with abiotic factors, in dynamic ecosystems such as Great Plains streams (Dodds et al. 2004).
The findings of this study strengthen the imperative need to put ecology into flow management through consideration of the NFR paradigm. Although the field of conservation biology was borne out of the implied application of ecological principles to conservation issues, recent research suggests such linkages are rarely made (Hintzen et al. 2019). Previous flow management literature has often focused on magnitude to the detriment of riverine organisms that require dynamic flow regimes (White andHowe 2004, Bunch et al. 2011). Our findings suggest that holistic restoration of flow regimes beyond just magnitude characteristics is needed within the reaches of our study and the species we modeled. Other studies have suggested that restoration of natural flow regimes results in ecological recovery (Propst andGido 2004, Kiernan et al. 2012) and our findings lend credence to this approach for PBS fishes within the southern Great Plains. Given the reproductive strategy of PBS fishes and their streamflow requirements, such species are particularly vulnerable to modified flows (Perkin et al. 2015a), but also presumably more responsive to natural flow restoration. Although our work was contextualized in the southern Great Plains, the methods utilized here provide additional support for applications of the NFR to benefit conservation of native fishes elsewhere, so long as both contemporary fish and flow data are present (Poff et al. 1997, Propst and Gido 2004, Kiernan et al. 2012). Our study findings provide guidance for more effective and meaningful future flow preservation and restoration endeavors by providing the necessary flow-ecology linkages needed to more broadly implement NFR-based conservation tools such as the ELOHA ). More broadly, flow-ecology relationships that ultimately lead to the implementation of streamflow management targets provide a mechanism for integrating ecological considerations into water planning and permitting processes.