Quantifying the effects of response diversity dynamics on ecosystem stability

Domain: chemistry
Published: 2026-03-17
Editorial outcome: minor_revision

Abstract
Variability in species’ responses to environmental changes (response diversity) is a critical feature in determining stability of ecosystem-level properties. However, quantifying response diversity remains challenging. A key issue is the context-dependent nature of species’ responses: a species reaction to a specific environmental factor can vary depending on other environmental conditions or interactions with other species. Consequently, these responses are dynamic rather than static. We propose a framework to quantify response diversity, incorporating time-variability in species' responses to biotic and abiotic factors and accommodating high-dimensional, nonlinear interactions in complex systems without assuming equilibrium. Using a 40-year set of monthly data from Lake Geneva, we demonstrate that response diversity of the phytoplankton and zooplankton community stabilizes their respective biomass within but not across trophic levels. Importantly, the stabilizing effect of response diversity varies over time, being context-dependent on time-varying environmental conditions. Our method provides an approach for uncovering mechanisms linking biodiversity and ecosystem dynamics.

Introduction
The role of biodiversity in ecosystem stability has been heavily contested, and the strength and direction of biodiversity-stability relationships can vary within and among systems1-9. To better understand diversity-stability relationships, some studies emphasize the importance of adopting a multifaceted perspective of stability2,3,10. Others focus on identifying the mechanism behind the inconsistent relationship between biodiversity and ecosystem stability4,8,9,11. Still others argue that traditional diversity metrics (e.g., species richness or Shannon entropy) are inadequate for uncovering the underlying stabilizing mechanisms12-15. Amid the ongoing debate, ecologists have increasingly embraced the concept of ‘response diversity,’ referring to variation in species’ responses to environmental change within a community, as a more appropriate metric for capturing key aspects of biodiversity most relevant to ecosystem stability12-15.

The history, concept, and impacts of response diversity on ecosystem stability have been thoroughly reviewed12; however, quantification in natural systems continues to pose challenges14. Although previous studies provided valuable insights into quantifying response diversity in experimental settings16, applying these approaches to real-world communities with high species diversity remains challenging. The mainstream approach for inferring response diversity relies on indirectly estimating trait diversity within a multi-dimensional functional trait space17. Defining, quantifying, and validating species' response traits and distinguishing them from effect traits continues to pose difficulties14. Methods have been developed using multivariate autoregressive models to analyze observational time series data of taxa (or functional group) abundances18,19. This approach provides an overall view of abundance variability (a proxy of biomass instability) explained by differentiated species’ responses to environmental variability but cannot quantify the sign and magnitude of individual species’ responses to specific environmental drivers.

Another approach quantifies response diversity by analyzing species’ performance (e.g., growth rate, production, and biomass) across environmental gradients, e.g., temperature15. Specifically, a species’ response to environmental variability is quantified as the first derivative of its performance curve, thus capturing direction and strength of the species’ response. By measuring variations among species-specific derivatives, community-level response diversity is then quantified. This approach can be further extended by considering multiple environmental drivers20.

Regardless of various measures of response diversity, existing statistical approaches and experimental frameworks assume a static value or function, such as a performance curve responding to a changing environmental variable. In reality, species’ responses are dynamic and context-dependent, influenced by various environmental conditions and species’ interactions21-23, as clearly emphasized in previous studies12,13. Thus, a critical knowledge gap lies in quantifying dynamic response diversity in natural systems, where species performance curves are not static, and both biotic and abiotic responses must be integrated.

In this work, we introduce a dynamic framework to quantify how response diversity in plankton communities varies over time. Analyzing long-term data from Lake Geneva, we demonstrate that response diversity stabilizes total phytoplankton and zooplankton biomass, with effects fluctuating seasonally and across years. Importantly, we find that biotic interactions are the primary mechanism mediating this stabilizing effect, revealing response diversity as a dynamic driver of community resilience to environmental changes.

Results
Quantifying temporal variability in response diversity of a community
Here, we introduce a framework to quantify temporal dynamics in response diversity using observational time series data (Fig. 1). Our method leverages nonlinear time series analysis rooted in dynamical system theory, known as Empirical Dynamic Modeling (EDM)24,25 (Methods). This approach accommodates the context-dependent nature of complex ecosystems, enabling quantification of time-varying species' responses to biotic (species’ interactions) and abiotic (environmental changes) factors simultaneously, overcoming limitations of previous methods that assume static species responses.

Our approach requires time series data of species abundances (or other performance metrics) and environmental variables measured at the same time scale (Fig. 1a). We then quantify the causal response of one species to each of the other species (biotic responses) and each of the selected environmental variables, such as physical, chemical, and climate indices, using MDR-Smap22 (Methods). By collecting responses across species, we construct a matrix that summarizes the response profile, encompassing species’ responses to biotic and abiotic factors at every time interval (Fig. 1b). Mathematically, each element in the matrix represents a partial derivative (Jacobian), quantifying the degree of change in a target species’ abundance in response to a small perturbation in a forcing variable, such as the abundance change in another species or perturbation from an environmental factor. This approach uses observational data to quantify how a target species responds to a changing environmental variable or another species at a specific time while accounting for effects of other environmental and biological factors (i.e., the partial derivative conditioned on influences of other forcing variables).

It is noteworthy that our method based on MDR-Smap differs from the classic framework using Multivariate AutoRegression (MAR)26 by relaxing assumptions of static species responses and equilibrium system dynamics. Instead, MDR-Smap enables estimation of time-varying response profiles (Fig. 1b), allowing visualization of temporal variability in these response profiles through graphical representation using dimension-reduction techniques (Fig. 1c).

At each specific time interval, the response profile reveals species in a community that differ in their responses to other species and environmental factors, with these differences providing the basis for describing response diversity. Importantly, temporal variability in response profiles leads to a time series of dynamic response diversity (Fig. 1d). Given sufficient variables (e.g., cross-trophic levels and various environmental variables), the response profile captures information on within-trophic level responses, cross-trophic level responses, and responses to various types of environmental factors (e.g., physical and chemical variables). Thus, response diversity can be categorized into various levels based on research objectives, such as biotic versus abiotic response diversity, top-down versus bottom-up biotic response diversity, and physical versus chemical response diversity (Fig. 1d).

The critical need to disentangle biotic and abiotic contributions to response diversity has long been emphasized15,18, yet operationalizing this distinction remains challenging. Our method enables precise tracking of how a target species i responds to a specific environmental factor j over a time interval from t to t+1 (e.g., J in Fig. 1b) by estimating response profile over each time interval. This allows examination of temporal variability in specific responses (Fig. 1e). Meanwhile, exploring response curves by plotting specific responses against an environmental variable deepens our understanding of context-dependent responses (Fig. 1f). These advances improve the efficacy of our approach for quantifying dynamic, non-static species responses.

Response diversity of Lake Geneva plankton as a case study
To demonstrate our approach, we analyze monthly time series data of phytoplankton and zooplankton together with environmental variables collected from Lake Geneva from 1974 to 2022 (Methods). In this study, we focus on genus-level data to ensure taxonomic consistency in this long-term monitoring project; we consider a total of 33 phytoplankton and 5 zooplankton dominant genera, each occurring in over 20% of time points. We include the following environmental variables: (i) physical variables- water temperature over 0-20m depth (T20m), maximum Brunt-Väisälä buoyancy frequency (MaxBV), mixed-layer depth (MLD); (ii) chemical variables averaged over 0-20m- dissolved inorganic phosphorus (PO4), silicate (Si), and total phosphorus (TP); and (iii) climate indices- Atlantic Multidecadal Oscillation (AMO), Arctic Oscillation (AO), El Niño-Southern Oscillation (ENSO), North Atlantic Oscillation (NAO), and Pacific Decadal Oscillation (PDO). These environmental variables (SI Fig. 1) have been shown to influence the Lake Geneva ecosystem27-29.

The environmental time series data illustrate the long-term warming (Fig. 2a) and eutrophication-reoligotrophication history (Fig. 2b) in Lake Geneva30. Using MDR S-map, we analyzed plankton time series data in combination with environmental variables, generating a response profile for each time interval, i.e., the matrix summarizing how each plankton genus responded to other genera and environmental variables. An example of a response profile at a specific time interval is shown (Fig. 2c). In practice, a time series of response profiles can be estimated (Fig. 1b), with seasonal and inter-annual variability visualized using a nonlinear dimension reduction technique, isomap analysis31 (Fig. 2d). Such dissimilarity in response among genera enables estimation for response diversity, which we quantified using the functional dispersion index32 for calculating response diversity over time (Fig. 2e and SI Fig. 1l-s).

We categorized response diversity (RD) into several levels. For example, we estimated RD-PhyToAll, which focuses on responses of phytoplankton genera to each other, to zooplankton genera, and to environmental variables (SI Fig. 1l); RD-PhyToPhy focusing exclusively on responses of phytoplankton genera to other phytoplankton within the same trophic level (SI Fig. 1m); RD-PhyToZoo focusing on responses of phytoplankton genera to zooplankton genera, considering only cross-trophic level responses (SI Fig. 1n); and RD-PhyToEnv assessing responses of phytoplankton genera to environmental variables apart from biotic interactions (SI Fig. 1o). Likewise, to describe responses of the zooplankton community, we estimated RD-ZooToAll, RD-ZooToPhy, RD-ZooToLZoo, and RD-ZooToEnv, respectively (SI Fig. 1p-s). We note that these categories are not necessarily independent and exhaustive, and their interdependency needs to be carefully considered.

Using Staurastrum (a phytoplankton genus) as an example, we illustrate temporal variation in its abundance responding to increasing total phosphorus (TP) loading in Lake Geneva, as captured by estimated partial derivatives (Fig. 2f). The magnitude of Staurastrum's positive response to elevated TP diminishes as TP levels increased (Fig. 2g). Specifically, at low TP levels, increases in TP level leads to a noticeable rise in Staurastrum abundance, whereas at high TP levels, further increases in TP result in minimal or no change in Staurastrum abundance. Such a saturating response reflects a non-static, context-dependent response of phytoplankton to nutrient availability. As Lake Geneva is currently in a meso-oligotrophic condition, our analysis predicts a noticeable increase in Staurastrum abundance with elevated TP inputs.

Response diversity and its role in stabilizing total biomass amid changing ecosystem conditions in Lake Geneva
With estimates of response diversities, we tested the hypothesis that the community's response diversity contributes to stabilizing the community's total biomass13,16,33 (Fig. 3a) by enhancing resilience to environmental change12. We used Chlorophyll a time series to indicate the total biomass of the phytoplankton community (SI Fig. 1t) and zooplankton sedimentation volume to indicate total biomass of the zooplankton community (SI Fig. 1u). Phytoplankton production measurements were used to represent primary productivity of the system (SI Fig. 1v). We employed the modified dynamic eigen-value (DEV) approach34 to quantify dynamic instability, also known as local Lyapunov instability, of phytoplankton or zooplankton total biomass at each time point. This instability is assessed by computing the dominant eigenvalue of the Jacobian matrix derived from time-delayed coordinate embedding (see Methods for details). A smaller DEV indicates a greater local stability of that variable at a specific time, with a threshold of 1 marking a transition to unstable conditions. In essence, a lower DEV reflects higher system resilience to perturbations. Using this approach, we can quantify the instability of the total biomass of phytoplankton and zooplankton through time, matching the time scale of response diversity. We avoided using temporal variability measures (e.g., coefficient of variation or standard deviation) to quantify instability, because these variability measures are sensitive to observational errors35 and require predefined selection of a moving window, which is hard to decide a priori. Our analyses of Lake Geneva data indicate that dynamic instability of phytoplankton biomass (called PhyInstability) is below 1 throughout the study period (SI Fig. 1w), indicating that the phytoplankton biomass fluctuated between relatively stable and less stable conditions. In contrast, the DEV of zooplankton biomass (called ZooInstability) fluctuated between ~0.8 and ~1.1 (SI Fig. 1x), suggesting shifts between stable and unstable conditions.

Whereas our objective was to examine effects of temporal dynamics in response diversity on biomass instability, we must consider complex interactive effects from other potential contributors, such as trophic interactions and environmental forcings, to accommodate context dependency in nonlinear dynamical systems. When examining PhyInstability, we included the influence of consumers (e.g., zooplankton biomass, ZooInstability, and zooplankton response diversity) and environmental variables. We first identified causal variables for the instability index using Convergent Cross Mapping (CCM) (Table S1). Then, using the MDR S-map, our analytical framework evaluated whether a community's response diversity helped stabilize its total biomass while accounting for these additional factors. Importantly, this approach quantified changes in the instability index caused by increasing response diversity (i.e., the partial derivative; see Methods) for each time interval.

Supporting the theoretical expectation, our results indicate that elevating the response diversity of the phytoplankton community (RD-PhyToAll) consistently contributes to stabilization of total biomass (Fig. 4a and SI Fig. 2). While increasing response diversity consistently enhances biomass stability, the magnitude of this stabilizing effect fluctuates over time with clear seasonal and interannual patterns (SI Fig. 2). For instance, magnitudes of RD-PhyToAll stabilizing effects are stronger in cool seasons (from Oct. to May of the next year) when MLD effects are minimal; in contrast, effects of RD-PhyToAll are weaker in warm seasons (Jun. to Sep.) when MLD effects on instability are stronger (SI Fig. 2b).

Importantly, our approach estimates effect of response diversity on stability at each time interval. Compared to linear analyses that assess the explained variance in stability by response diversity over the entire study period, our approach based on EDM offers a more detailed, time-specific perspective. On average, the magnitude of effect (i.e., effect size regardless of positive or negative) of RD-PhyToAll ranks highest among selected factors (SI Fig. 3a), highlighting the importance of response diversity in determining biomass stability. Interestingly, the response diversity of the zooplankton community (RD-ZooToAll) also helps stabilize total phytoplankton biomass, but this stabilizing effect is limited to the period from Jan to Mar. (SI Fig. 2b), with no consistent influence in other months. However, zooplankton biomass (a potential indicator of consumer effects on phytoplankton) and ZooInstability exhibit no consistent effects (Fig. 4a).

Whereas high response diversity in primary producers generally enhances stability of total biomass, high response diversity in consumers may have complex effects, either stabilizing or destabilizing total biomass of primary producers. Effects of consumer response diversity on producer stability likely depend on species composition, interaction characteristics, and environmental context (e.g., temperature or nutrient status). This is similar to findings in studies of diversity considering cross-trophic levels: high diversity in either predator or prey can lead to complex effects on stability at both prey and predator levels36,37.

Environmental and climate indices, including T20m, MLD, TP, and ENSO, influence phytoplankton biomass stability (Fig. 4a and SI Table 1). However, MLD, TP, and ENSO exhibited inconsistent effects on phytoplankton biomass stability, with positive effects only observed in certain months (SI Fig. 2b). These inconsistencies may arise from complex mechanisms with intertwined environmental conditions and predator communities38,39. In addition, plankton taxonomic composition exhibited strong seasonal succession40,41, further complicating stabilizing mechanisms. Surprisingly, increasing water temperature tends to stabilize phytoplankton biomass, which seems to contrast with large-scale meta-analyses across global aquatic systems indicating long-term warming destabilizing phytoplankton biomass11. However, we caution that our analysis operates at a monthly scale. Short-term temperature increases may favor specific species42 or promote vertical niche segregation among phytoplankton species43, leading to stabilization of total phytoplankton biomass; however, underlying physicochemical and biological processes warrant further investigations.

Similar to findings for the phytoplankton community, increasing the response diversity of the zooplankton community (RD-ZooToAll) also helps stabilize zooplankton total biomass (Fig. 4b and SI Fig. 4). In addition, primary productivity (PP) stabilizes zooplankton biomass, suggesting that enhancing food supply stabilizes zooplankton biomass; however, stabilization effects of PP are stronger in cool seasons (Oct. to May next year, SI Fig. 4b) when total phytoplankton biomass is low, again, highlighting context-dependent effects. TP, AMO, and PhyInstability exhibit inconsistent effects on zooplankton biomass stability, with inter-annual and seasonal trends observed for TP and PhyInstability (SI Fig. 4). These inter-annual trends may be linked to the lake's eutrophication history, highlighting different mechanisms under nutrient-rich versus nutrient-limited conditions. On average, the effect size of RD-ZooToAll ranks second highest among selected factors (SI Fig. 3b), only lower than stabilizing effects of primary productivity.

Our conclusion that increasing community response diversity helps stabilize their total biomass remains valid when examining more detailed categories of response diversity (SI Fig. 5). Nevertheless, response diversity resulting from biotic interaction has a more important role in stabilizing total biomass in Lake Geneva. Specifically, response diversity considering phytoplankton interactions (RD-PhyToPhy, SI Fig. 5a) stabilizes phytoplankton biomass, whereas different categories of zooplankton response diversity exhibit inconsistent effects on phytoplankton biomass stability. Likewise, response diversity based on zooplankton interactions (named RD-ZooToZoo, SI Fig. 5d) stabilizes zooplankton biomass, whereas phytoplankton response diversity in relation to environmental variables (RD-PhyToEnv) has inconsistent effects. Although our results provide strong evidence to support the idea that within-trophic level community response diversity helps stabilize their total biomass, cross-trophic analyses hint at the need for adopting a multi-trophic perspective on response diversity and species interactions to fully assess stability9,37,44.

Note that the response diversity index used in this study is not mathematically decomposable; in other words, RD-PhyToAll is not equal to the sum of RD-PhyToPhy, RD-PhyToZoo, and RD-PhyToEnv. Nevertheless, RD-PhyToAll serves as a comprehensive response diversity index, encompassing effects from both within and between trophic levels, as well as environmental factors. We noticed that response diversity related to environmental variables did not have a significant stabilizing effect on the community's total biomass. For example, RD-PhyToEnv and RD-ZooToEnv did not significantly affect PhyInstability and ZooInstability, respectively (SI Table 1). Rather, response diversity indices based on biotic interactions had a more important role in stabilizing the total biomass in Lake Geneva, emphasizing the importance of considering biotic interactions in estimating response diversity15,18. Although RD-PhyToEnv did not exhibit significant effects, RD-PhyToAll showed a stronger causal relationship with PhyInstability, i.e., higher CCM skill, compared to RD-PhyToPhy and RD-PhyToZoo (SI Table 1), suggesting that the abiotic effects become important when integrated with biotic effects.

Some of the identified effects of response diversity indices and biotic and abiotic variables (i.e., the S-map coefficients) exhibit long-term and seasonal variability (SI Figs. 2 & 4). Note that each S-map coefficient cannot be interpreted in isolation; rather, S-map results represent complex interactions and are context-dependent22,45. For example, although RD-PhyToAll consistently exerted stabilizing effects on phytoplankton biomass at a monthly scale, the strength of these effects varied seasonally and interannually (SI Fig. 2), depending on dynamic ecosystem processes and biotic interactions. A more challenging issue is interdependence of various variables. For example, mixed layer depth is influenced by temperature and wind, which are linked to fluctuations in climate indices, whereas phosphorus, an intrinsic environmental factor, may engage in feedback interacting with plankton dynamics. In fact, an ecosystem’s stability can potentially provide feedback to community and ecosystem processes; for instance, ZooInstability has complex effects on PhyInstability (SI Fig. 2). Here, we do not attempt to provide a comprehensive explanation of how biomass stability is determined by complex interactions of response diversity in combination with biotic and abiotic variables. Future work may consider a dynamic network approach that includes stability as a node, together with response diversity and biotic and abiotic variables as other nodes. Monitoring network dynamics may offer insights into important questions such as “what controls stability and how it operates” in natural systems.

Environmental effects on response diversity in Lake Geneva
Given the limited understanding of abiotic effects on response diversity14, we investigated how environmental conditions influence various categories of response diversity (Fig. 3b, SI Fig. 6, and SI Table S2), following the CCM and then MDR S-map framework. TP consistently emerged as an influential variable across all response diversity categories. Increasing TP positively affects RD-PhyToEnv, RD-PhyToZoo, RD-ZooToEnv, and RD-ZooToPhy, but negatively affects RD-ZooToZoo, and inconsistently affects RD-PhyToPhy. Enhancing phosphorus loading in Lake Geneva elevated various categories of response diversity, except for RD-ZooToZoo, suggesting a bottom-up mechanism promoting response diversity. Moreover, effects of increasing TP on various types of response diversity changed seasonally and interannually (SI Fig. 6). Further investigation is needed to understand why different categories of response diversity have consistent or inconsistent results. For example, some studies reported that eutrophication favors certain phytoplankton species while reducing species diversity and functional diversity46,47.

Phytoplankton biomass effects on RD-PhyToEnv and RD-PhyToZoo shifted from negative to positive during the reoligotrophication of Lake Geneva (SI Fig. 6b & h). In addition, phytoplankton biomass effects on RD-PhyToPhy are generally positive but weaker when TP is high (SI Fig. 6e). In contrast, for zooplankton response diversities, phytoplankton biomass effects are generally positive and seem stronger in the eutrophication period for all three categories of zooplankton response diversity (SI Fig. 6k, n, & h). This contrasting influence of phytoplankton biomass on phytoplankton versus zooplankton response diversity reveals a long-term increasing trend in phytoplankton response diversity while a decreasing trend in zooplankton response diversity. These trends coincide with the eutrophication-reoligotrophication history of the Lake Geneva ecosystem (SI Fig. 1d & f).

Water temperature negatively affected RD-PhyToZoo and the three categories of zooplankton response diversity; a shallow mixing layer depth (MLD) also suppressed RD-PhyToEnv, RD-ZooToEnv, and RD-ZooToZoo (SI Fig. 6). These results suggest that stronger water column stratification (a shallower MLD) potentially decreases plankton diversity, perhaps by impacting nutrient supply from deeper water48. In addition, short-term heat impulses may favor warm-adaptive plankton species and change plankton community structure49. Climate indices had no consistent effects, indicating that climate impacts on response diversity are complex and challenging to interpret.

Response diversity may deteriorate due to human- and climate-induced disturbances, such as warming, fishing, eutrophication, and oligotrophication, calling for targeted mitigation and restoration efforts. Monitoring response diversity, in turn, offers valuable feedback to policymaking and environmental interventions (Fig. 3c). These complex interactions are integral to a coupled social-ecological feedback system, having an important role in maintenance of ecosystem service and sustainability. The advantage of this dynamic framework lies in its ability to track response diversity over time, allowing evaluation of policy impacts on response diversity and ecosystem stability. Certainly, mechanisms underlying how environmental factors influence various categories of response diversity need detailed studies.

Considering positive and negative responses
The dispersion index does not distinguish between negative versus positive responses15, limiting its direct link to the mechanism by which response diversity stabilizes community dynamics. To partially address this limitation, we assessed sign heterogeneity in response (Methods); sign heterogeneity is prevalent in species’ responses (SI Figs. 7 and 8), supporting the presence of portfolio effect as a stabilizing mechanism.

In addition, we examined the response divergence index15,20 that accommodates signs of species’ responses. Instability of total phytoplankton biomass was, on average, buffered by increasing response divergence of phytoplankton (Divg-PhytoToAll in SI Fig. 9). These results are qualitatively consistent with those based on the response dispersion index, though divergence results had greater variability and occasional subtle destabilizing effects (comparing SI Fig. 9 versus Fig. 4a). Response divergence of zooplankton (Div-ZooToAll) also exhibited, on average, stabilizing effects on total phytoplankton biomass (SI Fig. 9), although both stabilizing and destabilizing effects were observed. Therefore, at cross-trophic levels, high response diversity in either predator or prey can have complex effects on stability of the phytoplankton community. In contrast, instability of total zooplankton biomass was not causally affected by response divergence (Divg-PhytoToAll and Div-ZooToAll; Table S1). Collectively, stabilizing effects of divergence are less clear than effects of dispersion. However, a key methodological limitation of the divergence index must be considered: it incorporates only extreme values (maximum and minimum) of species responses and ignores variance within the community, which is a critical component for the portfolio effect to operate and stabilize ecosystem properties. Consequently, although valuable for capturing sign heterogeneity, divergence provides an incomplete picture for demonstrating stabilizing effects of response diversity. Therefore, developing a comprehensive index that integrates both the sign and variance of responses remains a crucial yet challenging frontier for future research.

Methodological considerations
Our approach can complement trait-based approaches for quantifying response diversity. An intriguing research direction is to explore whether and to what extent the response profile derived from this approach can be linked to the trait space of a community. However, it is noteworthy that the response profile estimated from any dataset will critically depend on taxa and environmental variables included in the analysis. As emphasized in previous studies, careful consideration should be given to selection of environmental variables to ensure the calculated diversity reflects responses to realistic threats for the focal system12. In the Lake Geneva case study, plankton genera were used as the analytical unit. However, depending on ecological context and research objectives, users must determine the appropriate taxonomic resolution, with the option to use ecological guilds as the analytical unit50. Furthermore, we considered only dominant taxa, defined as those occurring in >20% of the monthly time points in this study. This choice reflects a methodological limitation in applying EDM that cannot accommodate taxa with time series including many zeros. Thus, effects (responses) of rare taxa on total biomass stability were omitted. Given that the biomass contribution of each rare taxon to total biomass is small, their influences on biomass stability are also likely to be minor. Although our example focused on the abundance of taxa, this analytical framework can be extended to analyze responses of other trait variables (e.g., body size, growth rate, etc.), provided that time series data for these traits are available.

In this case study, functional dispersion was used to measure response diversity. Our results, based on this widely used index, provide strong evidence that response diversity stabilizes total biomass. We have also evaluated the only existing response diversity index considering sign in response in the published literature- Divergence- but found inconclusive/ less consistent results. Although several approaches exist for calculating the dissimilarity of responses among taxa12,15,20,51, we did not attempt a comprehensive comparison in this study. Determining the more suitable approach for computing the response diversity index summarizing the responses (J in Fig. 1b) remains an open question.

Finally, we caution that all limitations of EDM apply to our approach. For example, sufficient time series data are needed, and rare taxa cannot be included in the analysis. More detailed statistical and mathematical discussions of the limitations of EDM are in relevant reviews24,25.

Future research directions
Several important future research questions warrant discussion. When the environment changes, the response diversity (or, to be more precise, the response profile) changes, not only because species composition changes but also because the species’ responses to environmental variables change, likely owing to variation in other environmental variables as well as species interactions. Being able to quantify contributions of changes in species composition versus changes in species’ responses for determining response diversity of a community would benefit policy-making and intervention of environmental conditions. In addition, research is needed to identify critical species contributing more to the response diversity, as well as critical environmental variables driving changes in response diversity. As such, the approach may help estimate if the system is under important stress and if the system would be able to remain resilient. As response diversity plays an important role in system stability, response diversity may serve as an early warning signal that informs potentially drastic changes in community dynamics. Long-term time series monitoring is warranted.

This study quantifies the response diversity (dispersion) at the beta level that considers dissimilarity among taxa in responses. However, this approach does not consider the capacity of responses of individual taxa to respond to biotic and abiotic variables, even though such response capacity could buffer destabilizing effects. We suggest that buffering effects at alpha level for a taxon could be quantified by carefully investigating response curves for a target species (e.g., Fig. 1f).

Concluding remarks
Building upon the theory of nonlinear dynamical systems, we developed an analytical framework to quantify taxon-specific, time-varying responses by integrating time series data from community and environmental variables. Based on the dissimilarity of responses among taxa, our approach quantifies time-varying response diversity, explicitly considering both biotic and abiotic responses. Empirical analyses of the Lake Geneva plankton ecosystem demonstrate that response diversity helped stabilize dynamics of total biomass, supporting a long-standing theoretical prediction that is hitherto difficult to verify in natural systems. Future work should assess the generality of the stabilizing effects of response diversity in other ecosystems.

Importantly, response diversity is not a static characteristic of the community but rather a dynamic property that changes over time in response to shifting ecosystem conditions or disturbance. Although our approach is developed to address the biodiversity-stability issue in ecology, this framework can be extended to other research fields and systems (e.g., social systems) in which response diversity is vital for maintaining system resilience. Since management strategies aimed at enhancing response diversity are crucial for sustainability33,52, our proposed framework provides a way to quantify context-dependent, dynamic response diversity, which offers insights into underlying mechanisms (Fig. 3 and section “Future research directions”) and supports policy-making for managing interconnected ecological and socio-economic systems.

Methods
Lake Geneva Data
We analyzed the average monthly time series data of phytoplankton and zooplankton communities and environmental variables from 1974 to 2022 for Lake Geneva. Sampling, counting, and measuring procedures are described in Rimet et al.53. Phytoplankton and zooplankton individuals were identified at the species level when possible. However, to ensure taxonomic consistency in this long-term monitoring project, we focused on genus-level data by aggregating species-level data in this study. We focused on the dominant genera for analysis, as EDM cannot be applied to time series with numerous zero values24. A total of 33 phytoplankton genera and 5 zooplankton genera were included, selected based on their occurrence in >20% of monthly time points. The temporally restricted species excluded from the analysis were also rare, with very limited relative abundance at each time point, accounting for only 0.038 ± 0.125% (mean ± SD) and 0.168 ± 0.371% of the total phytoplankton and zooplankton abundance, respectively. We included the following environmental variables: (i) physical variables- water temperature over 0-20 m depth, maximum buoyancy frequency (the highest value of buoyancy frequency observed in the temperature profile), mixed-layer depth (the depth where the water temperature differs by < 0.5°C from the surface temperature); (ii) chemical variables- phosphate, silicate, and total phosphorus; (iii) climate indices: Atlantic Multidecadal Oscillation, Arctic Oscillation, Niño-Southern Oscillation, North Atlantic Oscillation, and Pacific Decadal Oscillation. These environmental variables (SI Fig. 1) influence the Lake Geneva ecosystem27-29. This study did not include nitrogen nutrients because Lake Geneva is P-limited rather than N-limited, as nitrate concentrations have increased and stabilized whereas phosphorus has been decreasing54.

Empirical dynamic modeling (EDM)
Evidence abounds, indicating that ecosystem and population dynamics are nonlinear55-58. Moreover, a species’ responses to other species and environmental variables are nonlinear and context-dependent21,22,59,60. Thus, we employed the framework of Empirical Dynamic Modeling (EDM), specifically designed to analyze nonlinear dynamical systems24.

According to the theory of dynamical systems, a system can be described as a set of states (e.g., attractor) whose change over time is determined by a set of rules (e.g., equations). These rules describe how the state of the attractor changes as a function of interacting variables. Motion on the attractor can be projected onto the coordinate axis of associated variables, forming the time series; conversely, the time series of the interacting variables can be projected back to the multidimensional state space to recover the attractor61. Knowing the state and rules governing the underlying system is equivalent to knowing the time series of all interacting variables. When studying empirical systems, the set of rules is unknown. Meanwhile, very likely, only one or a few variables are measured, and their interactions are unknown a priori. However, it is possible to reconstruct a topologically invariant shadow version of the attractor using time-delayed scalar measurements of a single dynamical variable: X = (x, x ,…, x )T where x are observations at time point i, E is the embedding dimension, and tau is the time lag. This idea, called state space reconstruction (SSR), is based on Sauer, Yorke, and Casdagli’s extension62 of Takens’ theorem63 for dynamical systems. It is also feasible to reconstruct the attractor using combinations of multiple variables and their delays61,64. Interestingly, due to topological invariance, the reconstructed attractor preserves the essential mathematical features of the original dynamical system, such as the dominant eigenvalue65.

Convergent Cross-Mapping (CCM)
Convergent Cross-Mapping (CCM)66, a method designed for quantifying causality in nonlinear dynamical ecosystems, was used to identify causal variables (e.g., response diversity) for the target variable (e.g., stability index). CCM infers the causal relationship between pairs of variables from their empirical time series by measuring the extent to which the historical record of an effect variable, X, can reliably recover the states of a causal variable, Y66,67. That is, CCM is based on information recovery (i.e., effect variables contain encoded information on causal variables). This analysis implicitly accounts for influences of confounding variables. Specifically, CCM accounts for influences of confounding variables through lagged embeddings, e.g., (x , x , …), which represent historical effects of other variables in lagged terms, i.e., x = f(x , y , z , …), even if those variables are unobservable or not sampled.

The significance of CCM causality was obtained by comparing predictability of cross-mapping from the effect variable Y to the causal variable X against the predictability of a null CCM model, built upon surrogate causal variable X data. Such surrogate data were generated by randomly shuffling the residual from the fitted seasonal smooth spline. We considered the seasonality in null model generation because our analyses were operated at a monthly scale, where seasonal variation was strong and might obscure a true causal relationship68.

S-map
S-map involves generating an E-dimensional embedding (univariate or multivariate). The state space at time t is given by X = (x, x ,…, x )T. For each target time point t , the S-map algorithm computes a local linear model C that predicts the future value Y using the vector X from the reconstructed state space. That is,

x = C + sum E C x (1)

The linear model is fitted to other vectors in the state space. However, in contrast to linear regression models, linear approximation in S-map is done only locally by giving greater weighting to points that are close to the target point, X . The model C is the singular value decomposition (SVD) solution to the equation,

B = A∙C (2)

where B is an n-dimensional vector (n is number of observations) of weighted future values of Y for each historical point, t, given by

B = w(||X -X || )Y (3)

and A is then n × E dimensional matrix given by

A = w(||X -X || )X (4)

The weighting function w is defined by

w(d) = exp(-theta*d/d_bar) (5)

where ||.|| denotes the Euclidean distance and d_bar is the average distance between Y t* and all other vectors on the attractor. The extent of this weighting is tuned by the nonlinear parameter theta >= 0. Note that the model C is separately calculated for each time point, t*. Thus, C potentially differs for each time point of the time series. Coefficients of the local linear model provided by C are an approximation of local Jacobians45,59 that can be used to proxy interaction strength between embedded variables.

MDR-Smap
The original S-map algorithm requires that the system's optimal embedding dimension (E) equals the number of nodes used for inference (e.g., number of species). This condition does not always hold, particularly in diverse communities. To address this limitation, we employed the MDR S-map algorithm22, briefly summarized as follows: The MDR S-map consists of two steps. The aim of the first step is to obtain the multiview distance that operates at the optimal embedding dimension (E) that may not be equal to the dimensionality of a large interaction network (m=number of nodes). Under the curse of dimensionality, neighborhood relationships among data points cannot be precisely inferred from the distance, d (.) among high-dimensional data points in the manifold M that summarizes the dynamics of entire interaction networks. Therefore, (1-1) neighborhood relationships among high-dimensional data points X(t) are recovered from numerous low-dimensional state space reconstructions (SSR; M , M …, M ) (i.e., multiview SSR). Within these SSRs, (1-2) the distances (L2 norm, i.e., Euclidean distance) between data points are computed under the optimal embedding dimension E, enabling a more accurate evaluation of neighborhood relationships. By aggregating these distances, (1-3) multiview distances (dE) are obtained, and (1-4) the corresponding data weights (WE) are determined for subsequent S-map analysis. In the second step, the high-dimensional interaction strength (B) is estimated by locally weighted least square optimization used in the S-map algorithm (argmin (|| . || )2) with regularization (lambda[.]) applied. Specifically, (2-1) the weights, wE=exp(-theta*dE/mean(dE)), derived in the first step, are incorporated into the S-map optimization. (2-2) Under the regularization constraints (with lambda and alpha as penalty factors), (2-3) the high-dimensional local linear coefficients approximating interaction strengths are solved at each time step. Overall, the estimating equation of the MDR S-map is as follows:

B_hat(i) = argmin (||sqrt(W E i)(Y(i) - X B (i) )||^2 + lambda[alpha||B (i) ||_1 + (1-alpha)||B (i) ||_2^2]) (6)

Here, lambda is the penalized factor, and alpha is the adjusted parameter that balances regularization using either the L1 (||.||_1) or L2 (||.||_2) norm of the parameter vector B (i).

Estimating the response profile (Jacobian matrix) by extending MDR-Smap to include environmental variables
MDR S-map was originally developed to quantify time-varying species interaction22. Here, we extend the algorithm to include environmental variables as forcing factors. Computationally, this extension is straightforward- environmental variables are incorporated into X in the estimating equation of the MDR S-map, with the procedure for solving the S-map coefficients remaining unchanged. This allows us to quantify a species' response to both other species (i.e., biotic responses) and selected environmental variables (e.g., physical, chemical, and climate indices). Collectively, responses across all species in a community form a matrix that summarizes the response profile encompassing both biotic and abiotic responses at a given time; importantly, this response profile is computed for every time interval (Fig. 1b).

Temporal changes in response profile
We visualize the temporal changes of the response profiles using isomap analysis31 (a nonlinear dimension reduction technique); the biplot of isomap illustrates seasonal and interannual variation of response profiles (Fig. 2d). Here, we used isomap, as this approach accommodates nonlinear features of the data. Other ordination methods are possible.

Response diversity
One response profile (Fig. 1b) summarizes responses of each taxon to all the other taxa and environmental variables at a given time interval. The dissimilarity in response among taxa allows estimation for response diversity. Prior to computing the dissimilarity, we muted the J item to zero, following the previous study22, because J does not indicate intraspecific interaction (direct negative self-regulation mechanism) but represents the local, aggregated, intraspecific effect summarizing all potential pathways through which one species impacts itself via all species69. Here, we used the functional dispersion index32 for calculating response diversity through time (Fig. 2e and SI Fig. 1l-s). We adopted the unweighted version of the functional dispersion index, to avoid emphasizing any specific variable70. First, we define the multidimensional centroid C = sum J / n, where J indicates the element in the response matrix (Jacobian matrix), and n represents the number of taxa. Then, we calculated Z as the distance from each taxon to the centroid. Finally, response diversity for each time point is sum Z / n. Other indices to quantify dissimilarity at the community level are possible.

Response curve
For a given taxon, we can visualize temporal variation of changing abundance with respect to the increment of an environmental variable or another taxon (i.e., estimated partial derivative) (e.g., Fig. 2e). Then, plotting the partial derivative against an environmental variable or another taxon signals the context-dependence of that specific response (e.g., Fig. 2f).

Instability index of total biomass based on dynamical eigen-value (DEV) analysis
We employed the dynamic eigen-value (DEV) approach34 to quantify the dynamic instability of total phytoplankton or zooplankton biomass. The original S-map algorithm proposed in Deyle et al.45 aimed to estimate local Jacobians of the current state of the system using multivariate time series of interacting variables (as explained in the section on S-map above). However, the same algorithm can be employed for embedded time series from a single variable following the concept of lagged coordinate embedding (i.e., univariate embedding)34. Here, we aimed to estimate the coefficients of X_t+tau = J X_t + v with respect to matrix J and offset vector v. Since the reconstructed state space is based on a delay embedding of a univariate time series, most elements of the Jacobian, J, are zero, except the lower-off diagonal, which are unity, and the first row, which contains the non-trivial local linear approximations obtained from S-map. The same argument holds for the offset vector v. Here, the linearized dynamics read (7).

Using the S-map algorithm, matrix J of parameters j , …, j can be calculated as the associated S-map coefficients45,59. Because J is evaluated at each time point, it represents a time-varying local Jacobian J(t). From these local Jacobians, we calculate the dominant eigenvalue at each time step, enabling us to evaluate the system’s local Lyapunov stability. However, unlike Grziwotz et. al.34, our estimation of DEV was performed without using moving-window56,59 to ensure consistency with the MDR S-map framework.

Quantifying effects of response diversity on instability amid other biotic and abiotic variables
Although our objective is to examine the dynamical effect of response diversity on biomass instability, we cannot evaluate effects of response diversity in isolation. We must consider complex interactive effects from other potential factors contributing to biomass instability, in order to accommodate context dependency in nonlinear dynamical systems. Using MDR S-map, we examined whether a community's response diversity helps stabilize that community's total biomass while considering other variables. That is, in addition to response diversity, we considered effects of consumers (including zooplankton biomass, ZooInstability, and zooplankton response diversity) and physiochemical variables and climate indices. Some variables had a high correlation; for example, total phosphorus is highly correlated with silicate and phosphate, and maximum buoyance frequency is highly correlated with water temperature. To avoid redundancy in ecological interpretation, silicate, phosphate, and maximum buoyance frequency were excluded from this MDR S-map analysis. Prior to MDR-Smap analyses, CCM was used to identify significant causal variables for the instability index (SI Table 1). Specifically, the instability index is the target variable, whereas response diversity and other biotic and abiotic factors are forcing variables. Results of MDR-Smap provide the time-varying strength (a.k.a., Smap coefficients) of each selected variable (e.g., SI Fig. 2a). Violin plots were used to summarize strengths (e.g., Fig. 4). In addition, seasonal variations of those effects were visualized by making a violin plot for each calendar month (e.g., SI Fig. 2b).

Environmental effects on different levels of response diversity
For each of the six categories of response diversity (SI Fig. 1l-s), we examined time-varying ecosystem effects using MDR S-map. We considered biotic effects, such as phytoplankton and zooplankton biomass, and abiotic effects, such as physiochemical variables and climate indices. The analytical procedure is the same as explained in the section above. Prior to MDR S-map analyses, CCM was used to identify significant causal variables for the target category of response diversity (SI Table 2). Then, the target category of response diversity is the target variable, whereas biotic and abiotic factors are forcing variables in the MDR S-map analysis. Violin plots were used to summarize strengths (SI Fig. 6). Then, seasonal variations of those effects were visualized by making a violin plot for each calendar month (SI Fig. 6).

Index of sign heterogeneity in response at the community level
We quantified the sign heterogeneity of community responses to a given variable j using the following metric: Sign Heterogeneity = 1 - | sum sgn(J_ij,non-zero) | / n_j,non-zero (8).

Here, J_ij,non-zero denotes the non-zero response of species i to variable j within the Jacobian matrix, and n_j,non-zero is the number of all non-zero responses to variable j. The sign function sgn(⋅) returns 1, 0, or -1 for negative, zero, or positive inputs, respectively. The index ranges from 0 to 1: a value of 1 indicates an equal number of positive and negative responses (maximum sign heterogeneity), whereas 0 indicates all non-zero responses share the same sign (minimum sign heterogeneity). Cases where all responses are zero are excluded from the analysis.

Response divergence
Sign heterogeneity in response plays an important role in the insurance effect of response biodiversity15. Here, we adopted the “response divergence” proposed in Ross et al15 to accommodate the effects of positive and negative responses. The index is calculated as follows, Divergence_j = (max[J_ij] - min[J_ij]) - ||max[J_ij]| - |min[J_ij]||) / (max[J_ij] - min[J_ij]) (9).

where J_ij is the Jacobian element estimated from the MDR-Smap (Fig. 1b). If the range of J_ij does not span zero, then all species respond in the same direction to the environment (no sign heterogeneity in response); in such cases, Divergence is 0. In contrast, when J_ij are perfectly symmetric around zero, Divergence is 1. In summary, “response divergence” ranges between 0 and 1, with 1 indicating perfect sign symmetry around zero of responses to the environment. The multivariate version of response divergence can be computed simply by averaging divergence measures across all forcing variables20. We note that this index of response divergence considers only the maximum and minimum of the J_ij but ignores the variance of J_ij that is critical for the portfolio effect to operate for stabilizing total biomass.

Computation
MDR-Smap and CCM were executed using rEDM (ver. 1.2.3) and glmnet (ver. 4.1) packages in R. Response diversity was calculated using Matlab. Buoyancy frequency was calculated using rLakeAnalyzer (ver. 1.11.4.1). Data processing and visualization were conducted using R or Matlab.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This manuscript introduces a novel framework for quantifying dynamical response diversity in ecological systems, integrating nonlinear time series analysis (MDR-Smap) and applying it to a 40-year dataset from Lake Geneva. The authors propose that time-varying, context-dependent species responses, aggregated as functional dispersion (FDis), enhance biomass stability within trophic levels. The approach is methodologically ambitious and addresses a critical gap in biodiversity-stability theory. One core assumption and analytical choice warrant further clarification or revision to strengthen the study’s mechanistic grounding.\n1. The core response diversity metric used—functional dispersion (FD)—does not account for the direction (sign) of species' responses, thereby limiting its mechanistic alignment with ecological stability theory, which emphasizes the importance of compensatory (i.e., opposing) dynamics for stabilizing ecosystem functioning. While the manuscript briefly acknowledges this limitation (p.13), citing Ross et al. (2023), a more substantive discussion and methodological consideration would enhance the rigor of the stability argument. Specifically:\n(a) Could FD be recalculated with sign-sensitive weighting or decomposed into components reflecting directional divergence?\n(b) Are there observable patterns of sign heterogeneity within trophic levels, and do these correspond with observed stability trends?\n(c) Could supplementary figures be added to visualize the temporal distribution of positive vs. negative species-level responses?\nAddressing these questions would strengthen the causal interpretation linking response diversity to stability.\n2. The manuscript concludes that response diversity stabilizes biomass within but not across trophic levels. However, given that FD implemented is sign-agnostic, this claim may overstate the evidence. If consumer species respond in a synchronized or unidirectional manner, FD may still yield high values without reflecting compensatory potential—possibly obscuring destabilizing dynamics. A more careful interpretation acknowledging this limitation would better reflect the scope and constraints of the current metric.", "reviewer_2": "I found this a fascinating paper with the potential to have substantive impact on the field of community ecology. The most profound result is that diversity in the dynamics of responses to environmental change act in a way which stabilises food webs. This is an important result and has the potential to open up a new area of endeavour in community ecology, including achieving a better understanding of the mechanisms underpinning BD-EF relationships.\nThe analyses used are novel in this context but well established in a general sense, and provide a more permissive framework that accommodates more biology than would multivariate auto regression. The application of the DEV approach is a useful advance for describing resilience.\nThe focus is on ‘dominant’ genera (L134) in Lake Geneva. You need more detail on what was considered dominant and at least 2-3 sentences on the implications of omitting rarer taxa. This should be added to the summary of effects of taxonomic resolution at or around L359.\nThe paper is very dense and would benefit from the use of an experienced editor to try and simplify sentence structure and condense the text.\nThe paper is not very well written, and I have attempted to address the numerous grammatical errors. I propose a further detailed copyedit for grammar and clarity.\nTitle\nI find ‘dynamical response diversity’ to be a clunky phrase. Could it be ‘Quantifying the effects of response diversity dynamics on ecosystem stability’\nL27 The opening sentence is oddly phrased. Simplify as below:\n“Variability in species’ responses to environmental changes (response diversity) is a critical feature in determining stability of ecosystem-level properties.”\nL30 “context-dependence” should be “context-dependent”\nL30 “species responses” should be “species’ responses”\nL31 insert “other” between on and environmental\nL33 Change “time-varying” to “time-variability in”\nL35 Change “Using a 40-year monthly dataset from Lake Geneva” to Using a 40-year set of monthly data from Lake Geneva”\nL42 Change “While the role of biodiversity in ecosystem stability has been widely demonstrated, the strength and even sign of biodiversity-stability relationships can vary within and among systems” to “The role of biodiversity in ecosystem stability has been heavily contested and the the strength and direction of biodiversity-stability relationships can vary within and among systems”\nL44 Change “further the understanding of” to “better understand”\nL45 Change “of considering the multifaceted concept of stability” to “of adopting a multifaceted perspective of stability”.\nL50 Change “changes’ to “change”\nL54 Change “; however, quantifying it” to “; although quantification in”\nL57 Delete “Consequently,”\nL57 Change “the main-stream” to “The mainstream”\nL59 Delete “Nevertheless,” capitalise “Defining…”\nL60 Delete “Alternatively,” capitalise “Methods…:\nL63 Change “variations” to “variability”\nL64 Change “variations” to “variability”\nL66 Change “species” to “species’”\nL70 Change “variations” to “variability”\nL71 Change “environment” to “environmental”\nL74 Change “the experimental framework” to “experimental frameworks”\nL75 Delete “, however,”\nL78 Change “species” to “species’”\nL81 Change “time-varying” to “temporal variability in”\nL82 Change “dynamical” to “temporal dynamics in”\nL88 Change “species” to “species’”\nL106 Change “variations” to “variability”\nL111 Change “quantifying” to “describing”\nL111 Change “time-varying estimations of” to “temporal variability in”\nL112 Change “lead” to “leads”\nL116 Delete “the”\nL120 Delete “Furthermore,” capitalise “Our”\nL123 “Examining temporal variations in a specific response provides critical insights into the\nL124 variability of that response over time (Fig. 1e).” This is a truism. Change to “This allows examination of temporal variability in specific responses (Fig. 1e)”.\nL126 Change “response” to “responses”\nL126 Change “advancements demonstrate” to “advances improve”\nL148 Delete “their”\nL160 Change “emphasize” to “describe”\nL170 “We note that the response estimated by MDR S-map implicitly accounts for the effect from other factors and builds on which the saturating response reflects a non-static, context-dependent response of phytoplankton to nutrient availability” I don’t understand this sentence. Reword.\nL180 Change “changes” to “change”\nL182 Delete “In addition,” capitalise “Phytoplankton”.\nL185 Change “quantify the dynamic” to “quantify dynamic”\nL202 Change “dynamical effect of response” to “effects of temporal dynamics in response”\nParagraph beginning L202. This should all be in past tense.\nL215 Delete “their”\nL220 New paragraph at “Importantly…”\nL231 New paragraph at “Our results…”\nL234 Delete “the”\nL264 Delete “their”\nL292 Change “variations” to “variability”\nL378 “Forward-looking” is an odd section title. Change to “Future research directions”"}, "author_response": "Comments from reviewers are copied below with our point-by-point responses. We use black text for the reviewer’s comments and blue text for our responses. Changes in the text of the revised manuscript are underlined.\nReviewer #1 (Remarks to the Author):\nThis manuscript introduces a novel framework for quantifying dynamical response diversity in ecological systems, integrating nonlinear time series analysis (MDR-Smap) and applying it to a 40-year dataset from Lake Geneva. The authors propose that time-varying, context-dependent species responses, aggregated as functional dispersion (FDis), enhance biomass stability within trophic levels. The approach is methodologically ambitious and addresses a critical gap in biodiversity-stability theory. One core assumption and analytical choice warrant further clarification or revision to strengthen the study’s mechanistic grounding.\nResponse:\nThank you very much for your positive comments. Following the advice of Reviewer #1, we have provided additional clarification and conducted further analyses to strengthen the mechanistic grounding of the study.\n1. The core response diversity metric used—functional dispersion (FD)—does not account for the direction (sign) of species' responses, thereby limiting its mechanistic alignment with ecological stability theory, which emphasizes the importance of compensatory (i.e., opposing) dynamics for stabilizing ecosystem functioning. While the manuscript briefly acknowledges this limitation (p.13), citing Ross et al. (2023), a more substantive discussion and methodological consideration would enhance the rigor of the stability argument. Specifically:\n(a) Could FD be recalculated with sign-sensitive weighting or decomposed into components reflecting directional divergence?\nResponse:\nWe agree that accounting for the sign of species' responses is a critical aspect of response diversity. Following the advice, we have incorporated the “Divergence” index as proposed in Ross et al. (2023), implemented with the updated multivariate methodology (Polazzo et al 2024) into our analysis in the revision.\nPage 12:\n“Considering positive and negative responses\nThe dispersion index does not distinguish between negative versus positive responses15, limiting its direct link to the mechanism by which response diversity stabilizes community dynamics. To partially address this limitation, we assessed sign heterogeneity in response (Methods); sign heterogeneity is prevalent in species’ responses (SI Figs. 7 and 8), supporting the presence of portfolio effect as a stabilizing mechanism.\nIn addition, we examined the response divergence index15,20 that accommodates signs of species’ responses. Instability of total phytoplankton biomass was, on average, buffered by increasing response divergence of phytoplankton (Divg-PhytoToAll in SI Fig. 9). These results are qualitatively consistent with those based on the response dispersion index, though divergence results had greater variability and occasional subtle destabilizing effects (comparing SI Fig. 9 versus Fig. 4A). Response divergence of zooplankton (Div-ZooToAll) also exhibited, on average, stabilizing effects on total phytoplankton biomass (SI Figure 9), although both stabilizing and destabilizing effects were observed. Therefore, at cross-trophic levels, high response diversity in either predator or prey can have complex effects on stability of the phytoplankton community. In contrast, instability of total zooplankton biomass was not causally affected by response divergence (Divg-PhytoToAll and Div-ZooToAll; Table S1). Collectively, stabilizing effects of divergence are less clear than effects of dispersion. However, a key methodological limitation of the divergence index must be considered: it incorporates only extreme values (maximum and minimum) of species responses and ignores variance within the community, which is a critical component for the portfolio effect to operate and stabilize ecosystem properties. Consequently, although valuable for capturing sign heterogeneity, divergence provides an incomplete picture for demonstrating stabilizing effects of response diversity. Therefore, developing a comprehensive index that integrates both the sign and variance of responses remains a crucial yet challenging frontier for future research.”\nThe explanation of the “Response divergence”, Page 23:\n“Sign heterogeneity in response plays an important role in the insurance effect of response biodiversity15. Here, we adopted the “response divergence” proposed in Ross et al15 to accommodate the effects of positive and negative responses. The index is calculated as follows,\n𝐷𝑖𝑣𝑒𝑟𝑔𝑒𝑛𝑐𝑒 = (max[𝐽𝑖,𝑗]−min[𝐽𝑖,𝑗])−||max[𝐽𝑖,𝑗]|−|min[𝐽𝑖,𝑗]||) / (max[𝐽𝑖,𝑗]−min[𝐽𝑖,𝑗])\nwhere 𝐽 is the Jacobian element estimated from the MDR-Smap (Fig. 1b). If the range of 𝐽𝑖,𝑗 does not span zero, then all species respond in the same direction to the environment (no sign heterogeneity in response); in such cases, Divergence is 0. In contrast, when 𝐽𝑖,𝑗 are perfectly symmetric around zero, Divergence is 1. In summary, “response divergence” ranges between 0 and 1, with 1 indicating perfect sign symmetry around zero of responses to the environment. The multivariate version of response divergence can be computed simply by averaging divergence measures across all forcing variables20. We note that this index of response divergence considers only the maximum and minimum of the 𝐽𝑖,𝑗 but ignores the variance of 𝐽𝑖,𝑗 that is critical for the portfolio effect to operate for stabilizing total biomass.”\nIn summary, we recognize that a reliable sign-sensitive index for response diversity is not currently available. The development of an index that reflects directional divergence remains a considerable methodological challenge and an important direction for future research, beyond the aims of the present work.\nSI Figure 9. Violin plots showing effects of response divergence and environmental variables on total phytoplankton biomass instability. (a) Instability of total phytoplankton biomass was, on average, buffered by increasing response divergence of phytoplankton, zooplankton, and lake water temperature, although considerable variation was observed. Effects of other influential variables showed no consistent sign over time. Total biomass instability was quantified using dynamic eigen-values (Methods). Influential variables were identified with Convergent Cross Mapping to assess dynamic causality (Table S1), and their effects on instability were estimated by MDR S-map across time intervals. Because instability is the target variable, negative effects estimated by MDR S-map indicate buffering against instability. Violin color and the position of the black point represent the median effect, while the grey dots beneath each violin represent the monthly estimates. Divg-PhyToAll represents the response divergence of the phytoplankton community to all other plankton genera and environmental variables, and likewise for Divg-ZooToAll. For the definition of other variable names, see SI Figure 1.\n(b) Are there observable patterns of sign heterogeneity within trophic levels, and do these correspond with observed stability trends?\nResponse:\nWe have developed a new index to quantify sign heterogeneity in response. Our analysis reveals that sign heterogeneity is indeed prevalent in response to variables within and across trophic levels as well as to environmental variables (new SI Figures 7 & 8, included below). However, we did not find a simple, consistent correspondence between these patterns of sign heterogeneity and observed stability trends beyond those driven by seasonality. This result is theoretically expected: while sign heterogeneity is a necessary condition for compensatory dynamics and stability, it is not a sufficient condition. Other factors, such as the strength and variance of responses, are critical for stabilizing ecosystem properties.\nWe included this supportive analysis in the Discussion, Page 12, middle:\n“The dispersion index does not distinguish between negative versus positive responses15, limiting its direct link to the mechanism by which response diversity stabilizes community dynamics. To partially address this limitation, we assessed sign heterogeneity in response (Methods); sign heterogeneity is prevalent in species’ responses (SI Figs. 7 and 8), supporting the presence of portfolio effect as a stabilizing mechanism.\nThe explanation of the index of sign heterogeneity in response, Page 23:\n“Index of sign heterogeneity in response at the community level\nWe quantified the sign heterogeneity of community responses to a given variable 𝑗 using the following metric:\nSign Heterogeneity = 1− | ∑ sgn(𝐽𝑖𝑗,non−zero) | / 𝑛𝑗,non−zero\nHere, 𝐽𝑖𝑗,non−zero denotes the non-zero response of species i to variable j within the Jacobian matrix, and n is the number of all non-zero responses to variable j. The sign function sgn(⋅) returns 1, 0, or -1 for negative, zero, or positive inputs, respectively. The index ranges from 0 to 1: a value of 1 indicates an equal number of positive and negative responses (maximum sign heterogeneity), whereas 0 indicates all non-zero responses share the same sign (minimum sign heterogeneity). Cases where all responses are zero are excluded from the analysis.”\nSI Figure 7. Sign heterogeneity in phytoplankton responses. Panel a, d, and g illustrates the monthly time series of community-level sign heterogeneity in phytoplankton’s response to each zooplankton genus, each phytoplankton genus, and each environmental variable, respectively. Values consistently above zero confirm the existence of sign heterogeneity in response. Violin plots (b, e, and h) summarize the distribution of the sign heterogeneity in response to each of the forcing variables. See Index of sign heterogeneity in response at the community level for the explanation of the computation. The overall long-term distribution of 𝐽𝑖,𝑗 (c, f, and j) reveals that both positive and negative responses are prevalent. The sharp peak around zero indicates that weak responses are the most common. Zero responses are excluded from the analyses.\nSI Figure 8. Sign heterogeneity in zooplankton’s responses. Panel a, d, and g illustrate the monthly time series of community-level sign heterogeneity in zooplankton’s response to each zooplankton genus, each phytoplankton genus, and each of the environmental variables, respectively. Values consistently above zero confirm the existence of sign heterogeneity in response. Violin plots (b, e, and h) summarize the distribution of the sign heterogeneity in response to each of the forcing variables. See Index of sign heterogeneity in response at the community level for the explanation of the computation. The overall long-term distribution of 𝐽𝑖,𝑗 (c, f, and j) reveals that both positive and negative responses are prevalent. The sharp peak around zero indicates that weak responses are the most common. Zero responses are excluded from analyses.\n(c) Could supplementary figures be added to visualize the temporal distribution of positive vs. negative species-level responses?\nResponse:\nYes, we have included SI Figures 7 and 8 to visualize the sign heterogeneity in response. See the response above.\nAddressing these questions would strengthen the causal interpretation linking response diversity to stability.\nResponse:\nThank you for the very constructive suggestions. We have included new results and discussion for these issues in the revision. See our responses above.\n2. The manuscript concludes that response diversity stabilizes biomass within but not across trophic levels. However, given that FD implemented is sign-agnostic, this claim may overstate the evidence. If consumer species respond in a synchronized or unidirectional manner, FD may still yield high values without reflecting compensatory potential—possibly obscuring destabilizing dynamics. A more careful interpretation acknowledging this limitation would better reflect the scope and constraints of the current metric.\nResponse:\nAs shown in SI Figures 7 and 8, sign heterogeneity in response is prevalent in the cross-trophic level response. Thus, the concern of the reviewer#1 regarding the potential synchronized responses is not an issue in this case study of the Lake Geneva plankton ecosystem. To clarify this, we have added the following discussion:\nPage 12, bottom:\n“Response divergence of zooplankton (Div-ZooToAll) also exhibited, on average, stabilizing effects on total phytoplankton biomass (SI Figure 9), although both stabilizing and destabilizing effects were observed. Therefore, at cross-trophic levels, high response diversity in either predator or prey can have complex effects on stability of the phytoplankton community. In contrast, instability of total zooplankton biomass was not causally affected by response divergence (Divg-PhytoToAll and Div-ZooToAll; Table S1).”\nReviewer #2 (Remarks to the Author):\nI found this a fascinating paper with the potential to have substantive impact on the field of community ecology. The most profound result is that diversity in the dynamics of responses to environmental change act in a way which stabilises food webs. This is an important result and has the potential to open up a new area of endeavour in community ecology, including achieving a better understanding of the mechanisms underpinning BD-EF relationships.\nThe analyses used are novel in this context but well established in a general sense, and provide a more permissive framework that accommodates more biology than would multivariate auto regression. The application of the DEV approach is a useful advance for describing resilience.\nResponse:\nThank you very much for your positive comments.\nThe focus is on ‘dominant’ genera (L134) in Lake Geneva. You need more detail on what was considered dominant and at least 2-3 sentences on the implications of omitting rarer taxa. This should be added to the summary of effects of taxonomic resolution at or around L359.\nResponse:\nIn the original version, the definition of ‘dominant’ genera was provided in Methods. In the revision, we have additionally clarified this definition in the main text and discussed the caveats of omitting rarer taxa, as requested.\nPage 13, middle:\n“Furthermore, we considered only dominant taxa, defined as those occurring in >20% of the monthly time points in this study. This choice reflects a methodological limitation in applying EDM that cannot accommodate taxa with time series including many zeros. Thus, effects (responses) of rare taxa on total biomass stability were omitted. Given that the biomass contribution of each rare taxon to total biomass is small, their influences on biomass stability are also likely to be minor.”\nThe paper is very dense and would benefit from the use of an experienced editor to try and simplify sentence structure and condense the text.\nThe paper is not very well written, and I have attempted to address the numerous grammatical errors. I propose a further detailed copyedit for grammar and clarity.\nTitle\nI find ‘dynamical response diversity’ to be a clunky phrase. Could it be ‘Quantifying the effects of response diversity dynamics on ecosystem stability’\nL27 The opening sentence is oddly phrased. Simplify as below:\n“Variability in species’ responses to environmental changes (response diversity) is a critical feature in determining stability of ecosystem-level properties.”\nL30 “context-dependence” should be “context-dependent”\nL30 “species responses” should be “species’ responses”\nL31 insert “other” between on and environmental\nL33 Change “time-varying” to “time-variability in”\nL35 Change “Using a 40-year monthly dataset from Lake Geneva” to Using a 40-year set of monthly data from Lake Geneva”\nL42 Change “While the role of biodiversity in ecosystem stability has been widely demonstrated, the strength and even sign of biodiversity-stability relationships can vary within and among systems” to “The role of biodiversity in ecosystem stability has been heavily contested and the the strength and direction of biodiversity-stability relationships can vary within and among systems”\nL44 Change “further the understanding of” to “better understand”\nL45 Change “of considering the multifaceted concept of stability” to “of adopting a multifaceted perspective of stability”.\nL50 Change “changes’ to “change”\nL54 Change “; however, quantifying it” to “; although quantification in”\nL57 Delete “Consequently,”\nL57 Change “the main-stream” to “The mainstream”\nL59 Delete “Nevertheless,” capitalise “Defining…”\nL60 Delete “Alternatively,” capitalise “Methods…:\nL63 Change “variations” to “variability”\nL64 Change “variations” to “variability”\nL66 Change “species” to “species’”\nL70 Change “variations” to “variability”\nL71 Change “environment” to “environmental”\nL74 Change “the experimental framework” to “experimental frameworks”\nL75 Delete “, however,”\nL78 Change “species” to “species’”\nL81 Change “time-varying” to “temporal variability in”\nL82 Change “dynamical” to “temporal dynamics in”\nL88 Change “species” to “species’”\nL106 Change “variations” to “variability”\nL111 Change “quantifying” to “describing”\nL111 Change “time-varying estimations of” to “temporal variability in”\nL112 Change “lead” to “leads”\nL116 Delete “the”\nL120 Delete “Furthermore,” capitalise “Our”\nL123 “Examining temporal variations in a specific response provides critical insights into the\nL124 variability of that response over time (Fig. 1e).” This is a truism. Change to “This allows examination of temporal variability in specific responses (Fig. 1e)”.\nL126 Change “response” to “responses”\nL126 Change “advancements demonstrate” to “advances improve”\nL148 Delete “their”\nL160 Change “emphasize” to “describe”\nL170 “We note that the response estimated by MDR S-map implicitly accounts for the effect from other factors and builds on which the saturating response reflects a non-static, context-dependent response of phytoplankton to nutrient availability” I don’t understand this sentence. Reword.\nL180 Change “changes” to “change”\nL182 Delete “In addition,” capitalise “Phytoplankton”.\nL185 Change “quantify the dynamic” to “quantify dynamic”\nL202 Change “dynamical effect of response” to “effects of temporal dynamics in response”\nParagraph beginning L202. This should all be in past tense.\nL215 Delete “their”\nL220 New paragraph at “Importantly…”\nL231 New paragraph at “Our results…”\nL234 Delete “the”\nL264 Delete “their”\nL292 Change “variations” to “variability”\nL378 “Forward-looking” is an odd section title. Change to “Future research directions”\nResponse:\nAll the suggestions have been incorporated in the revision, and the English was edited by Prof. John Kastelic, a professional editor. Thank you very much."}, {"round": 2, "reviewer_comments": {"reviewer_1": "All of reviewers’ concerns on the earlier version of the manuscript were addressed in a satisfying manner. I do not have anything else to add. Congratulations for this nice work."}, "author_response": ""}]

Decision letter