Taxonomic identity, biodiversity, and antecedent disturbances shape the dimensional stability of stream invertebrates

The “dimensional stability” approach measures different components of ecological stability to investigate how they are related. Yet, most empirical work has used small‐scale and short‐term experimental manipulations. Here, we apply this framework to a long‐term observational dataset of stream macroinvertebrates sampled between the winter flooding and summer monsoon seasons. We test hypotheses that relate variation among stability metrics across different taxa, the magnitude of antecedent (monsoon) and immediate (winter) floods to stability metrics, and the relative importance of disturbance magnitude and taxonomic richness on community dimensional stability. Cluster analysis revealed four distinct stability types, and we found that the magnitude of floods during the prior monsoon was more important in influencing stability than the winter flood itself. For dimensional stability at the community level, taxonomic richness was more important than disturbance magnitude. This work demonstrates that abiotic and biotic factors determine dimensional stability in a natural ecosystem.

community dimensional stability. Cluster analysis revealed four distinct stability types, and we found that the magnitude of floods during the prior monsoon was more important in influencing stability than the winter flood itself. For dimensional stability at the community level, taxonomic richness was more important than disturbance magnitude. This work demonstrates that abiotic and biotic factors determine dimensional stability in a natural ecosystem.
Climate change is increasing the frequency and magnitude of floods (Blöschl et al. 2019;Marsooli et al. 2019) and droughts (Allen et al. 2019;Zipper et al. 2021), underscoring the importance of understanding ecological stability. The "dimensional stability" approach quantifies multiple stability metrics simultaneously to examine how they are interrelated. Individual metrics measure stability in different ways, which may be positively, negatively, or not related to each other (Donohue et al. 2013(Donohue et al. , 2016Domínguez-García et al. 2019). Understanding how stability metrics relate provides insight into how organisms respond to disturbances, which is fundamental to understanding global change.
Dimensional stability is critical for ecosystem functioning such as biomass production (Hillebrand et al. 2018), so it is imperative to understand the factors control it. We know dimensional stability can vary by disturbance type (Radchuk et al. 2019), species removals can have stabilizing and destabilizing effects on dimensional stability (White et al. 2020), and multiple stressors can produce interactive effects on dimensional stability (Polazzo and Rico 2021). Yet, prior work has leveraged theoretical models or short-term and small-scale experimental studies. Long-term observational studies on dimensional stability are exceedingly rare, but exceedingly important to understand the consequences of global change in natural ecosystems.
Stream invertebrate communities are a model study system for disturbance ecology (Fisher et al. 1982;Resh et al. 1988;Poff and Ward 1989). Streams are dynamic. Episodic floods wash out benthic organisms, and drying events desiccate them. A stream's hydrologic regime is a filter that structures communities, favoring certain species over others (Townsend and Hildrew 1994;Datry et al. 2014). Thus, stream invertebrates are used to test hypotheses that link species identity, disturbance regimes, and ecological stability.
Here, we use a long-term invertebrate dataset to investigate how taxonomic identity and disturbance magnitude influence population dimensional stability, and how disturbance magnitude and taxonomic richness influence community-wide dimensional stability. We calculate the stability of invertebrate populations after winter floods with five different metrics ( Fig. 1): resistance, the capacity to remain unchanged or amount remaining after a disturbance; resilience, the rate of return to predisturbance state; invariability, the ability to remain constant; persistence, the ability to remain near equilibrium; and recovery, the propensity to return to equilibrium after a disturbance. We then test the following hypotheses: Hypothesis 1: The dimensional stability of individual species varies because of tradeoffs between life-history strategies. Some taxa might be considered more stable when measured by one metric while less stable when measured by another. For example, taxa with life-history traits that confer rapid reproduction and growth should have greater resilience, while these same taxa are often the most prone to disturbances and the least resistant (Poff et al. 2018;Sarremejane et al. 2021).
Hypothesis 2: Disturbance magnitude will influence population dimensional stability. In streams with a predictable flooding season, some taxa may be more stable in a year with strong flood, while others may be more stable in years with a weak flood (or absence of a flood), since invertebrate taxa vary in widely in traits that make them prone or resistant to disturbances (Power et al. 2008;Sponseller et al. 2010).
Hypothesis 3: Variation in community-wide dimensional stability is influenced by disturbance magnitude and biodiversity. More extreme disturbances should select for a subset of taxa (Poff et al. 2018), leading to a simpler response and reduced dimensional stability, where stability metrics are tightly correlated. Conversely, a community with more taxa should display more varied responses to that disturbance within the community, leading to a more complex response at the community level, increasing dimensional stability when stability metrics are less correlated.

Study area
Sycamore Creek, AZ, USA, is one of the most extensively studied desert streams (Fisher et al. 1982;Grimm and Fisher 1989;Boulton et al. 1992;Stanley et al. 1997;Sponseller et al. 2010). It drains a mountainous watershed (elevation range 427-2164 m) and its watershed area is 424.8 km 2 at USGS stream gage 09510200 (Fig. S1). Located in the Sonoran Desert, it is characterized by a semi-arid, hot climate. Most precipitation occurs during the winter (December-March) or monsoon (July-September) seasons. Sycamore Creek is flashy, where annual discharge and number of floods is correlated with winter precipitation (Grimm and Fisher 1989;Sponseller et al. 2010).

Macroinvertebrate sampling
Benthic macroinvertebrates were sampled in Sycamore Creek from 1985-1999 to 2010-2019. Data from the first time period have been published in other studies (Grimm and Allen et al. Taxonomic identity, biodiversity, antecedent disturbances Fisher 1989;Boulton et al. 1992;Sponseller et al. 2010). Both sampling periods used the same location (Lat: 33.750524, Long: À111.507666). Samples were collected using a 0.004 or a 0.008 m 2 corer, 10 cm in depth, at five locations throughout a 100-m gravel run and field elutriated through a 250-μm mesh net. Identifications were made at different taxonomic resolutions, so we aggregated data to the lowest common resolution, resulting in 18 taxa (more details are provided in Supporting Information Methods).

Stability metrics
We used taxon abundance data to describe 120 postwinter flood trajectories from 18 taxa, and calculated five different stability metrics: resistance, resilience, invariability, persistence, and recovery (Figs. 1, S2). At the population level, we analyze these metrics separately, and simultaneously in a cluster analysis (see below). To measure stability across all five metrics simultaneously for the entire community, we calculated annual community-wide ellipsoid volumes following Donohue et al. (2013). Here, ellipsoid volumes measure the dimensionality of community stability, where the ellipsoid has five axes (one for each metric), and the volume of the ellipsoid is determined by the range of taxon values observed (each yielding one datapoint in five-dimensional space). Smaller volumes mean that the stability metrics are more closely related, and larger volumes means that they are not. More detailed descriptions are in Supporting Information Methods. In F, metrics are tightly correlated and dimensional stability is low. In G, metrics are uncorrelated, and dimensional stability is high. In F and G we present three metrics to aid in visualization, but here we analyze all five.

Flow metric calculation
We obtained discharge data from USGS stream gage 09510200, located $ 12 km downstream of the study site. Daily values were log transformed and normalized, and we extracted the seasonal signal using discrete fast Fourier transform (Sabo et al. 2008(Sabo et al. , 2017. We calculated mean daily residual values (σ) from the seasonal signal, a measure of flow extremeness ( Fig. S3; Table S1). We used these data to compute the following: flow residual (σ) from last winter pulse, mean flow residuals (σ) from prior winter season, and mean flow residuals (σ) from prior monsoon season.

Statistical analysis
All data wrangling, analyses, and visualizations were performed using R version 4.1.3 (R Development Core Team 2021) and packages are cited in Supporting Information Methods.

Taxonomic variation in individual stability metrics
To investigate how macroinvertebrate taxa varied in stability metrics, we used data from 13 taxa with > 5 years of data.
With mean-centered and standardized metric values, we used mixed effects models to test for differences in stability metrics between taxa while accounting for repeated measures (data from the same taxa across multiple years). Each metric was a single response variable, taxon was a fixed effect, and year was a random effect. To characterize how stability metrics were related to each other, we used K-means clustering analysis, where the unit of analysis was an individual taxon-year combination of metrics. We identified the optimal number of clusters using the NbClust package. To identify if stability cluster assignment varied among taxa, we used Fisher's exact test with Monte Carlo simulation (2000 times).

Disturbance magnitude influences on individual stability metrics
We used mixed effects models to investigate the effect of flow disturbance magnitude on individual stability metrics. For disturbance magnitude, we used mean daily residuals during the winter (1 December-31 March) and monsoon (1 July- taxon-year combinations. Chironomidae, Oligochaeta, and Acari had the greatest affinity for cluster 1; Helicopsychidae, Tipuloidea, and other Ephemeroptera had the greatest affinity for cluster 2; Fallceon, Oligochaeta, and Tricorythodes had the greatest affinity for cluster 3; and Gastropoda, other coleoptera, and Ceratopogonidae had the greatest affinity for cluster 4. 30 September) seasons prior to sampling, as well as the daily residual of the most recent flow pulse. We used the stability metric as the response variable, flow variables as fixed effects, taxon as a random effect. We assumed Gaussian distributed errors, which was confirmed by visually inspecting plots of model residuals. For each stability metric and flow variable combination, we first investigated whether a linear or quadratic model was a better fit to the data in a univariate model as determined by AICc. We chose the better fitting model if there was a significant difference, if not we used a linear model.

Disturbance and taxonomic richness effects on community-wide dimensional stability
We were interested in assessing how disturbance magnitude and taxonomic richness influenced community-wide dimensional stability as measured by ellipsoid volumes. We developed a suite of general linear models to compare in a model selection approach using community-wide ellipsoid volume as the response variable and four predictor variables: last flow-pulse residual, winter mean flow residual, monsoon mean flow residual, and taxonomic richness. For each predictor, we first compared univariate models with linear and quadratic functions, choosing the function with a better fit or linear for parsimony. We assessed the performance of all possible univariate and multivariate models using AICc, and natural log transformed ellipsoid volumes.
K-means cluster analysis identified four clusters of stability metric combinations (Fig. 2). Data visualizations suggest that cluster 1 had high resistance values, cluster 2 had high resilience and high recovery values, cluster 3 had high resilience but low recovery values, and cluster 4 had low persistence but high invariability values. Taxa significantly differed in their affinity for a cluster (Fisher's exact test, p < 0.001), and were often binned into multiple clusters across years (Fig. S5). Table 1. Summary of results from mixed effect models using flow magnitudes of the last pulse event and winter or monsoon season flow magnitude means as fixed effects and taxon as a random effect, and individual stability metrics as response variables. Marginal R 2 presents variance explained by the fixed effects without the random effect, while conditional R 2 values present the variance explained by both fixed and random effects. Significant values are indicated by bold text.

Metric
Flow

Disturbance magnitude influences on stability metrics
We found significant effects of flow magnitude on resistance, resilience, invariability, and recovery; but not persistence (Table 1). Effects of monsoon flow magnitude were significant on all stability metrics except for persistence, while the magnitude of winter flows and the last flow pulse were significant for just one metric each ( Fig. S6; Table 1). Monsoon flow magnitude had a significant convex quadratic effect on resistance, where resistance was greatest at high values (extreme monsoon flood years), lowest at intermediate values (typical monsoon flood years), and greater at low values (monsoon droughts). In addition, monsoon flow magnitude had a negative linear effect on resilience, but positive linear effects on invariability and recovery. Overall, the effects of the most recent flow-pulse magnitude and winter mean flow magnitude were less pronounced, as we observed a significant effect on one stability metric each. In contrast, the magnitude of the last flow pulse had a convex quadratic effect on recovery, where recovery was greatest at high flow magnitudes, lowest at intermediate magnitudes, and rose slightly as magnitudes decreased further. Winter mean flow magnitude had a negative effect on resilience.

Influence of disturbance magnitude and taxonomic richness on community-wide dimensional stability
A model with taxonomic richness as a univariate quadratic predictor of community-wide dimensional stability (measured as ellipsoid volume) outcompeted models that incorporated hydrologic variables (Table S3). Dimensional stability was greatest at intermediate richness values (Fig. 3, R 2 = 0.56).

Discussion
Using a long-term dataset, we found that stream invertebrates showed taxonomic variation in dimensional stability in response to hydrologic disturbances. Analyzing multiple stability metrics simultaneously showed that taxa could be classified into four different response types. Furthermore, we observed that flow magnitude during the prior monsoon season was more important than flow magnitude during the prior winter season in influencing stability metrics. Thus, the legacy of prior monsoon disturbances had stronger effects on invertebrate responses to winter floods than the winter floods themselves. Finally, we found that biodiversity was more important than disturbance magnitude in influencing community-wide dimensional stability.
We found mixed support for a "resilience-resistance" tradeoff (Miller and Chesson 2009). Two clusters were described by high resilience, but only one of these also had low resistance. Interestingly, in our study recovery was only high when resistance was very low, and resilience was high (cluster 2). This finding contrasts with others which found that recovery was high when resistance was also high (White et al. 2020), though others have shown that recovery can be high when resilience is also high (Hillebrand et al. 2018;Hillebrand and Kunze 2020). In our study, high resistance only occurred when persistence was also high (cluster 1), suggesting that in our study system resistance and persistence are related. Regarding other stability metrics, we also observed a trade-off between invariability and recovery. Cluster 4 was characterized by low recovery and high invariability. It was comprised of taxa such as Gastropoda and Ceratopogonoidae that become abundant and invariable in stable conditions but were less likely to fully recover, and less resistant and resilient, in response to disturbances (Fig. S5). In sum, both resistance and resilience are important mechanisms that allow invertebrate populations to persist and recover in response to floods in this flashy desert stream ecosystem.
Patterns of cluster assignments differed among taxa. Some were only classified into clusters that had high resilience, suggesting a more consistent responses to floods, such as the mayfly Fallceon (always in cluster 3) and the caddisfly Helicopsyche (either in cluster 2 or 3). Chironomidae was the taxon most frequently in cluster 1, suggesting it is the most resistant taxon in our study system. In other years it belonged in cluster 2, suggesting that chironomids can be both resistant and resilient. Chironomids are a large family with many diverse taxa, so perhaps it is unsurprising to see both mechanisms present in this taxon. Other taxa showed even more varied responses. For example, the mayfly Tricorythodes was assigned into each of the four clusters throughout the course of the study. Overall our analysis shows that while taxa vary Fig. 3. Plot of taxonomic richness against community dimensional stability, as measured by multivariate ellipsoid volumes of the five stability metrics. A larger ellipsoid volume is indicative of a more complex stability response produced by multiple mechanisms, whereas a smaller ellipsoid volume is indicative of a simpler response produced by fewer mechanisms. The maximum value predicted by the model is 6.60.
in dimensional stability, some are more consistent than others in their response to disturbances. This could be a product of aggregating taxa to coarser levels of taxonomic resolution in our dataset, which could either increase or decrease the observed variation within some taxa. Nevertheless, this finding mirrors that of other work that has found different flood responses by invertebrates with different life-history syndromes (McMullen et al. 2017;Poff et al. 2018).
Variation in flood magnitude is common. Indeed, extremely large floods can have very strong impacts on stream invertebrates, even leading to local extirpations (Poff et al. 2018). Long-term studies of the same system are required, however, to quantify how this variation can influence populations and communities. For example, Power et al. (2008) showed that when floods were large enough to reduce armored caddisfly larvae populations, other invertebrates such as tuft-weaving chironomids also became abundant. But when floods were weak, caddisfly populations became the dominant taxon. In the most recent study from our study system, Sponseller et al. (2010) showed that winter flood impacts on macroinvertebrate community structure were weaker than the overall effects of antecedent flooding and drought conditions. The work we present here extends this finding to directly show that the magnitude of monsoon flows, rather than the magnitude of winter flows, has a greater impact on the stability of invertebrates in Sycamore Creek. Indeed, after controlling for taxonomic identity, we saw impacts of monsoon flood magnitude on the resistance, resilience, invariability, and recovery of invertebrates from winter floods. In the Sonoran Desert, the summers are very hot and dry, so a respite provided by increased monsoon season flows may lead to invertebrate populations that are more stable in response to subsequent winter flooding. Unfortunately, we do not have data collected during this period to further investigate such mechanisms. Nevertheless, monsoon conditions in this system seem to set the stage for how invertebrates respond to winter floods. This is conceptually consistent with Eagle et al. (2021), who also applied the dimensional stability framework to a long-term dataset of stream invertebrates. With annual sampling they were not able to calculate populationspecific responses to specific flood events, focusing instead on community-oriented metrics such as the number of extinctions, colonizations, species turnover, and compositional similarity. Still, our work parallels one of the major findings from Eagle et al. (2021), where a single extreme flood event generated a legacy effect on stability in subsequent years. More broadly, this suggests that antecedent conditions should be considered in future work on disturbance impacts in streams.
There are few studies that have empirically investigated the controls on dimensional stability of ecological communities. When we measured dimensional stability at the community level, neither antecedent nor immediate disturbances influenced dimensional stability, and instead taxonomic richness was the most important factor. This is in direct contrast to our analysis of the five stability metrics calculated at the population level, where antecedent disturbances were an important factor. We should note, however, that the magnitude of the most recent winter flow pulse and taxonomic richness had some correlation in our data (r = 0.66). Our sample sizes are not large enough to investigate the direct effect of disturbances on dimensional stability along with indirect effects of disturbances mediated through taxonomic richness, yet both may exist. Nevertheless, in our analyses different factors were important for stability measured at different scales of biological organization. Polazzo and Rico (2021) studied how multiple stressors (pesticides, herbicides, and nutrient enrichment) influenced stability of phytoplankton, zooplankton, and macroinvertebrate communities in lentic mesocosms. In some ways their results are like ours, as they found that some disturbances influenced dimensional stability, while others did not. Nevertheless, we found a hump-shaped relationship between taxonomic richness and dimensional stability. Though this result is contrary to our hypothesis, it is consistent with other work that has found nonlinear relationships between ecosystem stability and richness. In particular, Pennekamp et al. (2018) manipulated species richness of microbes in microcosms and found a hump-shaped relationship between richness and stability. Although richness increased overall biomass, over the course of the warming experiment the most diverse treatment was also the least resistant to warming. Thus, as diversity increases the likelihood of having species that are sensitive to disturbances increases, which may lead to reduced community stability.
Since the dimensional stability conceptual framework was first presented by Donohue et al. (2013), it has rarely been applied to natural ecosystems. Most empirical work has been conducted in small-scale experimental micro-and mesocosms (Hillebrand et al. 2018;Polazzo and Rico 2021) and field exclosures (White et al. 2020). In our natural study system, we demonstrate that recovery is most strongly influenced by resilience, which contrasts with a meta-analysis of pulse disturbance experiments that showed that resistance was more important (Hillebrand and Kunze 2020). Moreover, we showed that antecedent disturbances can have stronger impacts on stability than the immediate disturbance itself, and that biodiversity has a clear effect on stability at the community level. Given that disturbance regimes are rapidly changing, our work indicates that while the relationship between disturbances and ecosystem stability is complex, long-term data can be employed to disentangle that complexity. Such work is required to predict the ecological consequences of global climate change in natural ecosystems.