Increasing drought favors nonnative fishes in a dryland river: evidence from a multispecies demographic model

Understanding how novel biological assemblages are structured in relation to dynamic environmental regimes remains a central challenge in ecology. Demographic approaches to modeling species assemblages show promise because they seek to represent fundamental relationships between population dynamics and environmental conditions. In dryland rivers, rapidly changing climate conditions have shifted drought and flooding regimes with implications for fish communities. Our goals were to (1) develop a mechanistic multispecies demographic model that links native and nonnative species with river flow regimes, and (2) evaluate demographic responses in population and community structure to changing flow regimes. Each fish species was represented by a stage-structured matrix, and species were coupled together into a multispecies framework through density-dependent relationships in reproduction. Then, community dynamics were simulated through time using annual flow events classified from gaged streamflow data. We parameterized the model with vital rates and flow–response relationships for a community of native and nonnative fishes using literature-derived values. We applied the simulation model to the Verde River (Arizona, USA), a major tributary within the Colorado River Basin, for the past half century (1964–2017). Model validation revealed a match between model projections and relative abundance trends observed in a long-term fish monitoring dataset (1994–2008). At the beginning of the validation period (1994), model and survey observations showed that native species comprised approximately 80% of total abundance. Model projections beyond the survey data (2008–2017) predicted a shift from a native dominant to a nonnative dominant assemblage, coinciding with increasing drought frequency. Trade-offs between native and nonnative species dominance emerged from differences in mortality in response to the changing sequence of major flow events including spring floods, summer high flows, and droughts. In conclusion, the demographic approach presented here provides a flexible modeling framework that is readily applied to other stream systems and species by adjusting or transferring, when appropriate, species vital rates and flow-event thresholds.


INTRODUCTION
Novel biological assemblages, comprised of species combinations differing from the past, are increasingly widespread across the world (Hobbs et al. 2009). Mounting evidence suggests that novel assemblages do not arise from a random reshuffling of species, but from heterogeneous rates of species losses and gains over time and space (Zavaleta et al. 2009, Dirzo et al. 2014, Moore and Olden 2017. Species invasions and climate change are primary contributors to novel assemblages, causing changes in species composition and non-systematic reductions in species richness (Dornelas et al. 2014). Consequently, modeling changes in community structure in response to non-stationary climate regimes, particularly in light of on-going species invasions, is a primary research challenge.
Demography-based approaches to modeling communities show promise because they represent fundamental relationships between population dynamics and environmental conditions (Keith et al. 2008, van de Pol et al. 2010, Yen et al. 2013. Demographic models utilize knowledge regarding the autecology of a species, such as rates of birth, growth, fecundity, and mortality of individuals (i.e., vital rates) to project population dynamics through time. Causal mechanisms can be readily incorporated by allowing vital rates to change as a function of environmental stochasticity or density dependence (Caswell 2001). Unlike modeling approaches that largely rely on combining single-species predictions to infer community change between two points in time, demographic approaches can model population responses to specific sequences of environmental events. This makes demographic approaches useful for understanding community responses to environmental change (Lytle and Merritt 2004, Yen et al. 2013, Wheeler et al. 2018. Representative communities can be modeled by linking individual species together via density dependence in space requirements, food, or some other limiting resource. The demographic community models are interaction neutral in the sense that pairwise species interactions arise from the model structure itself, rather than being specified a priori as parameters , Tonkin et al. 2018). This community-wide approach has demonstrated an ability to recover realistic patterns of community dynamics in freshwater ecosystems and shows promise for revealing how species interact under novel environmental conditions.
A rapidly changing climate that includes more frequent and severe droughts and flooding is poised to reshape fish communities of temporary and perennial rivers in dryland regions (Datry et al. 2014, Kominoski et al. 2018. Climate models project that decreased snow accumulation and higher evapotranspiration rates in spring and summer months will lead to more frequent and severe droughts in the southwestern United States, especially when combined with growing human water demands (Christensen et al. 2004, Seager et al. 2013, Udall and Overpeck 2017. In fact, climate-driven changes to streamflow in the Colorado River Basin have already been observed (Solander et al. 2017), where increased low-flow anomalies and decreasing habitat connectivity threaten native fish persistence (Jaeger et al. 2014, Ruh ı et al. 2015 and may favor nonnative fishes in the future (Ruh ı et al. 2016). Lifehistory traits have proven useful to understand past (Gido et al. 2013) and predict future responses of native and nonnative fishes to environmental change (Whitney et al. 2017). Thus, we expect that demographically based community models that account for species-specific relationships with hydrology will help to predict past and future changes in fish assemblages.
Dryland rivers of the southwestern United States are a flashpoint for the conservation challenges associated with changing river hydrology and a proliferation of introduced species. Widespread dam construction, flow diversions, and surface and groundwater abstraction for growing human populations have significantly altered environmental regimes in the region, creating conditions that threaten native species persistence and promote nonnative fishes (Minckley and Deacon 1968, Olden and Poff 2005, Strecker et al. 2011. As a result of nonnative species introduction and their establishment and proliferation from reservoirs, the number of nonnative species equals or exceeds native species in most watersheds throughout the southwestern United States (Pool et al. 2010, Walsworth andBudy 2015).
Here, we developed a mechanistic multispecies demographic model to evaluate how native and nonnative fish populations change in response to changing flow regimes. We modeled a fish community by examining species-specific vital rates that varied as a function of flow regimes describing patterns of drought and flooding. Density-dependent relationships coupled populations together into a multispecies framework, and community dynamics were simulated through a sequence of flow events defined by the streamflow record. We validated the model by comparing model projections to empirical data from a long-term fish monitoring program. We then examined how contemporary climate-driven change in the sequence of drought and flood events affects the composition of native and nonnative species in the community. Results from this study provide important insight into how changing flow regimes and invasive species may threaten the future of endemic native fishes in the Colorado River Basin, and broadly highlight the utility of multispecies demographic modeling approaches in ecology.

Study system and species
The Verde River, a tributary within the Colorado River Basin, drains over 17,000 km 2 of central Arizona (Fig. 1). The perennial mainstem river runs approximately 270 km through private, state, tribal, and United States Forest Service lands, originating in Big Chino Wash (1325 m a.s.l.) and flowing to its confluence with the Salt River north of Phoenix, Arizona (402 m a.s.l.). We focused our study on the unregulated upper Verde River mainstem, where development is primarily limited to livestock grazing and reductions in baseflows are a result of groundwater withdrawals (Garner et al. 2013).
Highly valued for its natural beauty and management priorities as a Wild and Scenic River, the Verde River is a focal point for the conservation of endemic native fishes (Averitt et al. 1994, Turner andList 2007). At least 12 fish species were historically native to the system, but the fish assemblage is changing rapidly, and only five native species have been observed since 1997 (Rinne 2005). By contrast, numerous nonnative fishes are present in the Verde system, including several species of centrarchids, ictalurid catfishes, and minnows (Rinne 2005). The Verde River has been the focus of detailed monitoring efforts starting in the 1990s, where fish community composition in relation to flow and habitat requirements have been examined annually for the period 1994-2008 (Stefferud and Rinne 1995, Rinne and Miller 2006, Neary et al. 2012). These surveys included seven sites in the upper Verde River, encompassing a spatial extent of approximately 60 river kilometers and representing all valley types and habitats occurring within the Verde River ( Fig. 1; Neary et al. 2012). River discharge representative of our study area has been measured in the mainstem upper Verde River continuously since 1963 (USGS gage 09503700).
We examined the seven most common fish species in the upper Verde River, collectively representing, on average, 87% of stream reach biomass (Gibson et al. 2015). Native fishes included desert sucker (Catostomus clarki), Sonora sucker (Catosotmus insignis), and roundtail chub (Gila robusta); species that are endemic to the Colorado River Basin. Nonnative fishes included yellow bullhead (Ameirus natalis), green sunfish (Lepomis cyanellus), smallmouth bass (Micropterus dolomieu), and red shiner (Cyprinella lutrensis); species with known ecological impacts (Ruppert et al. 1993, Dudley and Matter 2000, Propst et al. 2015. These seven species represent a range of body sizes and major life-history trade-offs between size and age at maturity (growth), juvenile survivorship (survival), and fecundity (reproduction; Olden et al. 2006a). This range of functional traits is represented in dryland streams throughout the western USA.

Modeling framework
All seven species populations were modeled simultaneously to represent the fish community of a 1-km river reach. The foundation was a stage-structured matrix population model for each species, modified to incorporate environmental variability and density-dependent relationships (Caswell 2001). The general model structure was adapted from a multispecies matrix population model originally designed to model riparian vegetation population dynamics as a function of river hydrology (Lytle and Merritt 2004, Tonkin et al. 2018), but with a number of important modifications described below. Model implementation followed four major steps: parameterization, simulation, validation, and perturbation analysis. Parameter values were based on the flow regime, species biology, and biomass estimates. Model simulations projected fish community composition at annual time steps for each water year in the flow record . Model validation compared population and community model projections against the entire record of long-term fish surveys for the seven upper Verde sites from 1994 to 2008. Finally, perturbation analysis evaluated the effect of uncertainty in parameter estimates on the model output. These steps are described in detail below.
Parameterization: streamflow, fish, and flowresponse relationships.-Hydrology in the upper Verde River mainstem is characterized by relatively steady, spring-fed baseflow, with highflow events that vary in magnitude and timing among years in response to winter and summer storm runoff ( Fig. 2A; Goetz and Schwarz 2018). Each time-step in the model represented one year. Using the historical flow record, each year was classified into a flow-event year type according to the timing, magnitude, and duration of flows for a water year (1 October-30 September) from USGS gage number 09503700 near Paulden, Arizona, USA. Flood events were defined using discharge thresholds based on recurrence intervals, and drought events were defined by the duration of baseflow (Fig. 2B). Spring high-flood events were years in which the maximum of late winter/early spring discharge (1 January-30 April) exceeded 19.8 m 3 /s, which has a 4-yr return interval during the spring time window (following Brouder 2001). Spring medium-flood events, corresponded roughly to bankfull flows, had a maximum discharge that exceeded 6.2 m 3 /s, with a 2.5-yr recurrence interval in the spring Ingersoll 1998, Neary et al. 2012). Summer (and monsoon season) high-flow events (1 May-30 September) exceeded a maximum discharge of 5.7 m 3 /s, representing a 4-yr recurrence interval in the summer. Drought events were categorized by the absence of floods, when lowflow conditions (i.e., 25th percentile of flows following Bêche et al. 2009) persisted for a continuous duration of 40 or more days (i.e., exceeding the 75th percentile duration of lowflow events). Nonevents occurred by default if years were otherwise not defined by flood or drought. A year type with both spring flood and summer high-flow events was possible, but all other flow events were mutually exclusive. This resulted in six possible years: spring high flood, spring medium flood, summer high flow, spring flood and summer high flow, drought, and nonevent (Fig. 2B).
The life-history adaptations of fishes to the flow regime are directly related to their vital rates (Lytle and Poff 2004). Because of the relationship between growth, survival, reproduction, and the flow regime, the use of vital rates inherent to each species supports the transferability of our model to other riverine systems (Mims andOlden 2012, Chen and). We conducted an extensive literature search of peerreviewed articles, graduate theses, and professional reports for each species, or closely related congeners, to determine parameter estimates for the vital rates used in the model (Table 1; Appendix S1).
Our integrated approach to demographic modeling requires knowledge of how these vital rates vary according to key components of the hydrologic regime. In another literature search, we reviewed studies about fish responses to flow ❖ www.esajournals.org components, specifically low-flow (drought) and high-flow (flooding) events in the Colorado River Basin as much as possible. Limiting our review to the Colorado River Basin minimized variability in vital rates that would be introduced by nonnative fish responses observed in other physiographic regions (Chen and Olden 2018). The breadth of information revealed in this literature review included observed relationships between abundance and discharge, effect sizes on changes in abundance in response to highand low-flow events, and timing of reproductive behavior (Table 2). This analysis allowed us to assign different vital rates according to types of flow-event years.
The riverine flow regime affects two key life stages that ultimately shape fish population structure and dynamics: juvenile survival and recruitment, and adult survival to reproduction (Schlosser 1985, Humphries et al. 1999). Thus, we allowed the vital rates for each species to vary according to literature-informed relationships between major flow events and fish abundance. To implement this, we used species-specific flow modifiers to adjust baseline vital rates for each year type in the flow record (Appendix S1: Table S1). For example, juvenile survival and recruitment of roundtail chub and smallmouth bass are influenced by the magnitude and timing of high-flow events. Spring flooding increases recruitment of juvenile roundtail chub, whereas elevated summer flows increase mortality of juvenile smallmouth bass (Brouder 2001, Smith et al. 2005. Adult survival is most affected by extended droughts. Low-flow events reduce survival and abundances of species Stefferud 1998, Ruh ı et al. 2015). Droughts typically create conditions where fish suffer because of limited resources and degraded water quality conditions as stream reaches are reduced to shallow isolated pools Minckley 1974, Lake 2003). The magnitude of these flow modifiers is set to reduce or increase mortality by a factor related to trends and effect sizes that were obtained from studies conducted in the region, and occasionally from other watersheds when data for nonnative species were otherwise unavailable (Table 2; Appendix S1: Table S1).
Model structure and simulation.-Each fish species was represented by a three-stage demographic matrix containing species-specific vital rates (Fig. 3, Table 1). The three life stages in the life-cycle model represented important ontogenetic shifts for each species. These stages corresponded to: year-1, juvenile recruitment into the population; year-2, sub-adults at first maturity; and year-3 or older adults of fully mature and reproductive individuals (Fig. 3). The parameters within the life-cycle model were adjusted by species to reflect real differences in population stage structure and traits such as lifespan and age at maturity (details follow). Individuals in each stage of the model occupied biomass calculated from length-weight relationships (Appendix S1: Table S2). Lengths at each stage corresponded to literature-reported values of young-of-year and/or immature fish (stage 1), average length at age of maturity (stage 2), and the average length of mature adult individuals from samples in the Upper Verde River (Appendix S1; unpublished data from Gibson et al. 2015).
Transition probabilities in the matrix differed according to the major flow events, which allowed recruitment and survival to vary Notes: Parameters included the Gonadal-Somatic Index (GSI ij ), the conversion factor from individuals to biomass for eggs (D ej ), stage 1 juveniles (D 1j ), stage 2 sub-adults (D 2j ), stage 3 (D 3j ), egg and larval mortality (M ej ), and mortality (M ij ) for each stage i and species j, transition probability of adult mortality in stage 3 (a 3j ) for each species j, and starting adult population biomass (B 3j ) for each species j selected from a negative binomial distribution with mean (k) and dispersion parameter (j). References for vital rate parameter estimates can be found in Appendix S1.
according to the hydrologic conditions in a particular year. The seven single-species demographic matrices were coupled via densitydependent relationships in reproduction, limited by the total biomass carrying capacity of a representative river reach. In this way, the reproductive output of each species declined as the total aggregate fish biomass of the entire community approached a reach-wide carrying capacity (Eq. 2). This generalized density dependence is analogous to exploitative competition experienced by organisms competing for a single limiting resource (Hardin 1960). A similar approach has been implemented in multispecies models using spatial density dependence .
The model had several assumptions regarding carrying capacity and vital rates. First, we assumed that the carrying capacity of the reach was limited by the amount of total fish biomass that could be sustained. Carrying capacity, K, was set to the average total biomass found in a 1km reach from surveys in nine replicate 100-m sampling sites located in the Upper Verde River and conducted in 2012 (Gibson et al. 2015). We chose to use average total biomass because some reaches will naturally be more or less productive and suitable for fishes than others. Second, population growth was limited by a densitydependent function in the reproductive term (fecundity, F) so no recruitment occurred if total biomass in the reach was greater than or equal to K (Eq. 2). Third, baseline mortality was the same for all life stages within each species, except for the egg and larval phase. The combination of egg and larval mortality was calculated as part of fecundity so that the starting population of adults produced sufficient offspring to equal their replacement after stage 1 and stage 2 baseline mortality was taken into account. Last, all species had a 1:1 sex ratio, and an individual's life cycle could be completed in a 1-km reach for all species. This is a reasonable assumption given that fishes in the Verde River typically occupy  (1987) Notes: Information presented includes the directionality of response to generalized flow attributes reported in the literature (+/À), list of relevant fish species, the dependent response variable to the flow attribute, and supporting references. Flow modifier values assigned according to these relationships are provided in Appendix S1: Table S1. Fig. 3. Generic life-cycle graph for fish species, with one-year projection intervals. Each circle represents a life stage (S i ): juveniles (S1), sub-adults (S2), and mature individuals (S3). Each arrow represents a transition between life stages with probabilities of growth and survival (G i ), surviving and remaining in the same stage (P i ), or reproduction (F i ). All fish reproduce by S3. Dotted lines represent arrows for species who begin reproducing in their first or second year of life, at S1 or S2, respectively. small home ranges and most fish move <1 km (Jaeger et al. 2014, Comte and. Species' biomass, rather than species abundances, was the currency for the model framework (although abundance and biomass values were interchangeable using the stage-specific weight of individuals). For each species j at each stage i for a given year type k, the biomass change from time t to time t + 1 was given by: where B ij was the total biomass (g) of the species at the corresponding life stage. Stage-specific fecundity (F ij ) was linearly density dependent on the total fish biomass in the community: where i = 1, 2, . . ., m was an index representing stages for each of j = 1, 2, . . ., n species in the community. Reach-wide carrying capacity K was the maximum aggregate biomass attainable for all species in all stages combined. Species began spawning at the stage that corresponded to the age of first maturity (Appendix S1: Fig. S3). Therefore, in each species' matrix, the number of fecundity (F) terms corresponded to the number of years the species can reproduce and was zero otherwise. For example, red shiner may begin reproducing in its first year of life, so the matrix had three fecundity terms, whereas yellow bullhead begin reproducing in their third year of life, so they had one fecundity term at stage 3. The proportion of egg biomass (f j ) produced by all females in the population after accounting for mortality during the first year depended on the Gonadal-Somatic Index (GSI), the proportion of gonad mass to total body mass, the proportion of female individuals (0.5), and the vital rate or combined egg and larval mortality (M ej ). The egg biomass of each species was converted to stage 1 biomass to account for fish development and growth. We converted egg biomass using egg density (D ej ) in units of number of eggs per gram and the average weights of stage 1 individuals (D À1 1j ; Table 1). Egg density of each species was the average value calculated from literature-reported values (Appendix S1) of the number of mature eggs in ripe females, divided by the total ovary mass. Ovary mass was reported or calculated using her length-weight relationship and GSI. If the relationship between total length and number of mature eggs was not published, egg density was estimated from the means of reported values.
The growth rate and survival probability (G) transitions from stage 1 to stage 2 and stage 2 to stage 3 were (4) Survival was one minus mortality (M i,j ) multiplied by a flow modifier Y ijk for stage i at species j for flow event k (Table 2; Appendix S1: Table S1).
The probability of surviving and remaining in the adult stage was related to the literaturereported lifespan of each species: where the probability that adults stay in stage 3 (a 3j ) is the reciprocal maximum age of species j after accounting for the first two life stages of the model. At the end of each model run, biomass output was converted to abundance for subsequent data analysis and interpretation. We ran 1000 iterations of the model, simulating the community from 1964 through 2017. Each iteration began with different initial population biomass to account for spatial variation. Initial biomass for each species j was sampled from a negative binomial distribution with mean k and dispersion parameter j converted from the mean ❖ www.esajournals.org and variance in relative abundance of the longterm Verde dataset for all seven sites and 15 yr of data (n = 105). As with a Poisson distribution, the negative binomial is appropriate for counts of organisms that occur randomly over time or space, but the negative binomial allows the variance to exceed the mean. Model validation.-Model validation was performed by comparing population and community projections to data from long-term fish surveys from seven sites in the Upper Verde River between 1994 and 2008 (Stefferud and Rinne 1995, Rinne et al. 1998, Neary et al. 2012. Data from 2002 were omitted from the analysis because only two of seven long-term sites were surveyed. Relative annual species abundances, averaged across the seven sites, were compared to model projected relative abundances according to the sum of non-juvenile (stages 2 and 3) individuals because juveniles (stage 1) are underrepresented by electrofishing survey methods due to their more cryptic behavior and small body size (Bonar et al. 2009). Relative abundance was chosen over absolute abundance because survey effort was not reported and was not always consistent between years, despite the use of a standardized collecting protocol. The strength of association between observed and modeled relative abundances was assessed using Spearman rank correlations. Spearman rank correlations are reported for each species across years and each year across species between 1994 and 2008, again omitting 2002. We also reported root mean square error (RMSE) and coverage (C), the percentage of 95% confidence intervals that overlap true values, to evaluate model performance.
Perturbation analysis.-The influence of model uncertainty was evaluated with respect to mortality rate, GSI, flow event thresholds, and biomass carrying capacity. These values were chosen because they are the most likely to affect transitions between life stages (i.e., mortality rate), fecundity/reproduction (i.e., GSI), species responses to environmental changes (i.e., flowevent thresholds), and outcomes in species dominance (i.e., biomass carrying capacity). We conducted a direct perturbation analysis (sensu Regan et al. 2003, Bond et al. 2014 to estimate the effects of parameter uncertainty by adjusting each value by AE10% of the starting value. Flow events were evaluated by adjusting thresholds of spring and summer high-flow events and duration of low-flow drought events. Uncertainty was quantified by the proportional change in species' relative abundances in response to parameter adjustments. Perturbation analysis was favored over analytical sensitivity analysis (via partial differentiation of a vital rate with respect to population growth rate) because it can be used to evaluate uncertainty in non-matrix elements, such as carrying capacity, and it does not require calculation of the long-term population growth rate (

RESULTS
Model population projections reflected compositional trends observed in the long-term fish surveys conducted from 1994 to 2008 (Fig. 4). During this validation period, trends in observed relative abundances of green sunfish, smallmouth bass, yellow bullhead, and roundtail chub were aligned with the simulated model populations (Fig. 4C, D, F, G). Model performance was best for green sunfish, with highest values for correlations (r), coverage (C), and lowest prediction errors (RMSE; Fig. 4F). The remaining species demonstrated varied correlations, coverage, and prediction errors. Predictions for relative abundance for roundtail chub and smallmouth bass were significantly correlated with observed data and had low to moderate prediction errors (Fig. 4C, G). Both smallmouth bass and yellow bullhead had high coverage (Fig. 4D, G). Desert sucker predictions correlated moderately well with observations and had moderate prediction errors and coverage (Fig. 4A). By contrast, the model tended to overestimate Sonora sucker and underestimate red shiner relative abundances (Fig. 4B, E). Relative abundances of red shiner, a species with small maximum body size that often occurs patchily in large schools of individuals, varied greatly among sites, evidenced by the large error bars (Fig. 4E). Sonora sucker had the lowest correlation with survey data and lowest coverage because the model did not capture an observed shift toward lower relative abundance in the middle of the survey period.
Community structure shifted from native dominant to nonnative dominant over the course of the model simulation period. Native species comprised approximately 80% of total abundance at the beginning of the validation period (1994) for both model and survey data (Fig. 5).
Then, after a series of drought and nonevent flow years, there was a marked transition toward nonnative dominance. At the end of the observed survey data, after a spring high-flood event in 2005, native species rebounded to represent 50-60% of the assemblage abundance (Fig. 5). Fig. 4. Modeled (lines) and observed (dots) relative abundances through time for native (A-C) and nonnative (D-G) fish species. Model and observed data include 95% confidence intervals, gray bands, and error bars, respectively. In the upper right corner of each panel, root mean square error (RMSE) is a measure of the difference, and Spearman rank correlations (r) is a measure of the association strength, and coverage is a measure of the percent overlap, between observed and model values.
Model simulations followed the general trend of decreasing, then increasing native abundances, but the magnitude of the model fluctuations was dampened compared to survey observations (Fig. 5). Despite some differences in magnitude, overall there was a strong correlation between rank order abundances of species between model projections and survey observations (mean r AE SE: 0.67 AE 0.018) with significant annual correlations for half (7 out of 15 yr) of the validation period (Appendix S1: Table S3). The model performed poorly between 1999 and 2005, coinciding with a period in the flow record from 1996 to 2004 that lacked spring flood events, which are important for juvenile (stage 1) survivorship of native Sonora sucker, desert sucker, and roundtail chub (Table 1; Appendix S1: Table S1). Beyond the period of the fish survey (post-2008), the modeled community projection, in response to observed streamflows, continued to gradually shift from native to nonnative dominant assemblage to the end of the flow record in 2017. Put in the perspective of flow events, the last ten years of the model had a drought frequency of 30%, compared to 9% drought frequency of the full flow record (Fig 2B).
Perturbation analysis on species vital rates had the largest effects on overall community composition, a larger effect than was seen by perturbing flow thresholds or biomass carrying capacity ( Table 3). Perturbation of desert sucker and green sunfish mortality and red shiner reproduction (GSI) had the largest influence on community composition overall, but the most affected species was always the one whose parameter was being perturbed. For example, when desert sucker's mortality rate experienced a 10% decrease, their relative abundance increased by 58% while other species compensated for this increase with a decline of 14-23% relative abundance. A decrease in red shiner's GSI parameter resulted in a 44% decrease in their relative abundance with a compensatory increase of 6-30% for the other species. This species, with its short lifespan, small body size, and high reproductive rate, was the most influential among changes in GSI, but also the most influenced by increased carrying capacity. When carrying capacity for a reach was increased, red shiner relative abundance increased by 17%, due to its high reproductive rate. Other members of the community had smaller changes in relative abundance, with declines of 4-9%.
Perturbation analysis in flow event thresholds demonstrated that community composition was most influenced by increasing the threshold of spring medium-flood events and decreasing the threshold of summer high-flow events. By increasing the spring threshold, two fewer spring flow events and two more nonevent years occurred during the hydrologic record. This influenced community composition by increasing yellow bullhead and green sunfish relative abundance by 39% each, accompanied by smaller increases in relative abundance for the other nonnative species, and decreases in relative abundance for the native species (Table 3). By contrast, decreasing the threshold for summer high-flow events resulted in three fewer nonevent years and a decrease in yellow bullhead and green sunfish relative abundance by 35% and 52%, respectively, but with minor or no changes to the relative abundance of the other species (Table 3).

DISCUSSION
Demographic models provide new opportunities to better understand species responses to shifting environmental conditions such as Fig. 5. Ratio of native species to nonnative species abundance. Presented are model predicted mean (black solid line) and 95% confidence intervals (gray band) and observed mean (points) and 95% confidence intervals (whiskers). Correlations between observed and modeled data for each year are available in Appendix S1: Table S3. climate-induced changes to hydrologic regimes and human-caused flow alteration (Shenton et al. 2012, Bond et al. 2018. Approaches that enable temporally specific predictions of species responses to specific flow sequences are likely to be most useful in practice (Wheeler et al. 2018). Here, we used demographic modeling to simulate how hydrologic drought and flood events interact with species' vital rates to shape native and nonnative fish composition in a dryland river. Although biotic interactions, such as competition and predation, were not specified a priori in the model framework, pairwise species interactions arose from the model structure itself due to the assumption of aggregate density dependence. The approach captured community trends using trade-offs in flow-related mortality of different life stages, as evidenced by predicting ranked abundances of species over a decadeand-a-half time period. Mismatches between modeled and observed numbers were most noticeable for species with vital rates that are difficult to obtain, or arose from potential life histories that lead to observation error associated with field survey methodology. Mechanistically representing trends in community response to environmental drivers using independently published vital rates opens up the possibility for hypothesis testing and scenario analysis for exploring management options for multiple species at once. In our simulated model of the Verde River, overall patterns in community structure, over a 54-yr period, demonstrated that more frequent drought events and fewer spring flood events created conditions where native fishes fared poorly compared to nonnative fishes. This finding is supported by empirical research in dryland rivers of the same region (Propst et al. 2008, Gido and Propst 2012, Gido et al. 2013). In the model, both nonevent years and droughts supported population growth by nonnative species. Nonevents represent years of steady baseflows that favored the survival of all juvenile nonnative fishes and adult life stages of yellow bullhead, red shiner, and green sunfish. Similarly, during drought years, native suckers (Sonora sucker and desert sucker) and roundtail chub experienced high mortality rates in both juvenile and adult life stages. By contrast, nonnative species had lower mortality rates for juvenile life stages, and depending on the species, had higher or lower mortality rates for adult life stages during drought. Higher mortality rates of large native species (i.e., the suckers) during drought years facilitated a compensatory response by nonnative species. That is, newly available portions of carrying capacity in the form of biomass were taken up by small-and medium-bodied fishes with high reproductive and low mortality rates during drought years (i.e., red shiner and green sunfish).
Past research demonstrates that droughts modify the spatiotemporal connectivity of riverine habitats, ultimately driving patterns in the composition and trophic structure of fishes (Matthews and Marsh- Matthews 2003, Rolls et al. 2012. As streamflows decline, fishes move to seek deep refuge pools (Labbe and Fausch 2000, Magoulick and Kobza 2003, Marshall et al. 2016. As pool habitats contract, the density of organisms initially increases with several consequences. First, smaller volumes of water concentrate prey, providing a food subsidy and increasing survival probability for young-of-year fishes for species spawning in warm months, or for extended breeding seasons (Schlosser 1985, Craven et al. 2010). In the Verde River, all nonnative species included in this study spawn during the warm spring-summer months. Of these, red shiner has the longest breeding season being a serial spawner, thereby performing well in drought years. Second, crowding intensifies predation and competition for resources among fishes (Magoulick andKobza 2003, Matthews andMarsh-Matthews 2006). In rivers of the southwestern United States, nonnative fishes such as smallmouth bass, green sunfish, and yellow bullhead outcompete or consume native fishes, especially juvenile life stages or smallbodied species, leading to local extirpations and community change when streamflow variability declines (Eby et al. 2003, Stefferud et al. 2011. Eventually, the combination of low flows and negative species interactions may lead to local native species extirpation. For instance, in the Verde River, the small-bodied native spikedace (Meda fulgida) has already been lost from the system but maintains viable populations where it has been repatriated in streams without nonnative fishes (Neary et al. 2012). As a corollary, high flows tend to favor native species, sometimes at the cost of nonnative species. In particular, spring floods may displace nonnative species and delay their reproduction (Propst and Gido 2004). This suggests that low-flow conditions favor species invasions (Bêche et al. 2009, Diez et al. 2012 and can have prolonged consequences even after the cessation of drought (Humphries and Baldwin 2003).
In the model results, low frequency or absence of spring flood events was detrimental to native fish populations. As observed in the perturbation analysis, minor increases in thresholds for medium spring flood events increased the number of nonevent years in the hydrologic period of record. This resulted in lower relative abundances of native species, including both sucker species and roundtail chub. The absence of spring flood events for nine years during the model simulation led to lower population growth rates for native species. The lack of spring flood events, which favored juvenile recruitment, resulted in native species experiencing higher mortality rates compared to nonnative species during these times. Therefore, given that the average lifespan of native species in this study range from 5 to 10 yr, it is likely that low recruitment potential for the nine-year period was detrimental to the persistence of these fish populations.
We found that native fish populations trended downward as nonnative fishes became dominant during a period of more droughts (three occurrences) and nonevents (two occurrences) starting in 2008. This decline was only punctuated by minor increases in native fish abundances following 2010 and 2015 spring floods. Because our model simulates a closed community with no rescue effects, once native fish population growth rates experienced a precipitous decline for several consecutive years after droughts, species within the model community were more likely to go extinct than in observational studies Notes: Relative changes resulted from AE10% change in parameters for species mortality (M ij ), Gonadal-Somatic Index (GSI ij ), flow-event thresholds, and maximum biomass carrying capacity (K). For each species, the table displays the directional change (i.e., + or À) that had the greatest effect on species responses. Ranked responses represent the largest absolute change in relative abundance across species in descending order. Bold text highlights the species with the greatest response to each affected parameter. Abbreviations for flow event year (SP_HF, spring high flood; SU_HF, summer high flood; SP_MF, spring medium flood; DR_length, length of the drought period). of riverine fish communities. Unlike our model community, fish communities observed under prolonged and supra-seasonal droughts can be resilient in spatially connected landscapes. For example, communities may eventually recover to pre-drought composition if deep refuge pools persisted during drought or other locations acted as sources for recolonization following the cessation of drought (Davey and Kelly 2007, Matthews et al. 2013, Rolls et al. 2016). However, persistent low-flow conditions have caused longterm changes in community structure where species adapted to low-flow conditions expanded their range and increased in abundance with complementary reductions in species favored by high flows (Lawson and Johnston 2016).
Demographic models that link species vital rates with hydrology have the advantage of addressing differential responses of native and nonnative fishes to streamflow variability because they leverage mechanistic associations between populations and specific flow events. This modeling approach implicitly includes species interactions via the density-dependent relationship between reproduction and biomass carrying capacity but does not explicitly include trophic interactions such as predation. For this reason, as with most modeling approaches, the strength of our conclusions depends on how much variability in population and community structure is explained by the unknown true contribution of flow events to demographic responses compared to other factors. Separating the influence of flow conditions independent of other extrinsic environmental factors, including species interactions, is inherently difficult (Chen and Olden 2018). However, a number of observational studies have demonstrated that the importance of biotic interactions in shaping fish communities is often overridden by environmental forcing (Grossman et al. 1998, Ruh ı et al. 2015, Giam and Olden 2016, Bino et al. 2017). Our goal was to develop and present a transferable approach to modeling multiple species within a community, an approach that allows for the incorporation of non-stationary environmental change and demographic variability. Continued research will help establish the relative strength of different drivers to demographic and community responses and help improve model performance and interpretation of model results.
A persistent challenge of community modeling is to simultaneously represent all individual species (Olden et al. 2006b). Individual species correlations between observed and simulated model data were moderate, but it is important to note that our model predictions were derived from independently reported parameters. In other words, the demographic multispecies model we have presented is unlike a statistical model where empirical data are used to fit model parameters. Rather, model performance depended solely on independent information, including that from literature-based vital rates, initial population sizes taken from a probability distribution, and empirically informed flow modifiers. Discrepancies between observed and modeled results were at least partially due to uncertainty in parameter estimates. Despite long interest in flow-ecology relationships, the strength of our understanding of these associations for all species in a community remains limited (Davies et al. 2014, Rosenfeld 2017. Vital rate and flow-dependent transition parameters can be adjusted as new empirical knowledge is gained. Although no one method is perfect, natural observations and in-stream or mesocosm experiments remain critical for acquiring vital rate information as a function of streamflow and will help further improve model predictions of community structure (Shenton et al. 2012, Poff 2018, Wheeler et al. 2018).
Discrepancies between model predictions and observed values may also occur because of observation error. For example, we found that model predictions underestimated the abundance of red shiner and overestimated the abundance of Sonora sucker. Together, these discrepancies led to a period where native dominance in the community was overestimated. Models for red shiner demonstrated the poorest performance, which may be due to highly variable sampling efficiencies due to the schooling behavior of this species. Sampling bias may also affect the observation of species with large maximum body sizes such as Sonora sucker and smallmouth bass. Older age groups of these species may be observed less than expected by model predictions because they occupy habitats too deep to sample with standard seining and backpack electrofishing techniques. Although fish surveys document low-recruitment periods during dry years with smaller sample sizes of Sonora sucker in wadeable habitats, model projections of adult Sonora sucker were relatively stable. This apparent stability may reflect mechanisms related to the temporal storage effect, whereby species are able to store up gains during favorable periods to persist during non-favorable periods, enabling multiple species to coexist in variable environments Warner 1981, Warner andChesson 1985). Here, high-recruitment periods (i.e., spring flood years) allow Sonora sucker populations to remain stable because long-lived adults survive through periods more favorable to other species.
In our model, the sequence of flow events appeared to support long-term multispecies coexistence, by favoring different species in years with different environmental conditions, at least over the 54-yr time frame examined in this analysis. In explorations over longer timescales (centuries or longer), the model is expected to eventually predict the complete dominance of a single species-the species with traits that resulted in the highest stochastic population growth rate for that particular hydrologic regime -due to the fact that we are modeling finite populations in a finite reach under stochastic conditions. Body size, age at maturity, and fecundity represent strong trade-offs in the life-history strategies of fishes (Olden et al. 2006a). Therefore, one might expect that the largest or earliest maturing species with the highest fecundity or reproductive rate, respectively, would always outcompete all other species in the community (Cushing 1992). However, because modeled vital rates are tightly coupled to flow conditions, patterns reflected observed species coexistence dynamics. In the model, the sequence and frequency of particular flow events played an important role in determining which species persisted over longer timescales.
Demographic models can be applied to any component of river ecosystems to explore species responses to changing hydrology, or other drivers of persistence and mortality. Where vital rates are known, they have been applied to individual species or functional guilds of invertebrates (McMullen et al. 2017), riparian plants , Tonkin et al. 2018, and fish (Yen et al. 2013). Demographic models are also flexible because different parameters may be applied to forecast the effects of environmental change in other ecosystems (McMullen et al. 2017). For example, one could apply our model to another stream system by adjusting or transferring, when appropriate, species vital rates, and flow event thresholds specific to the hydrograph and ecology of the system of interest. The flexibility of these models, and ability to predict non-stationary temporal dynamics, makes them useful for exploring scenarios of environmental change and the outcomes of various flow management interventions into the future. In-stream flow management plans to benefit species and ecosystems are often challenged by other competing demands on river flows, and population models that forecast the outcomes of various flow futures are a valuable and much-needed tool for decision makers.

ACKNOWLEDGMENTS
We are grateful to Jerome Stefferud and John Rinne for providing us with the long-term fish survey data in the Upper Verde River, Arizona. Tim Walsworth and Keith Gido provided supplemental data that helped verify information on fish growth and body size. Thanks to two anonymous reviewers for valuable comments and Christian Torgersen, Laura Prugh, and Gordon Holtgrieve for insightful suggestions. Funding was provided by the U.S. Department of Defense (SERDP RC-2511).