# Causal networks of phytoplankton diversity and biomass are modulated by environmental context

Chun-Wei Chang, Takeshi Miki, Hao Ye, Sami Souissi, Rita Adrian, Orlane Anneville, Helen Agasild, Syuhei Ban, Yaron Be'Eri-Shlevin, Yin-Ru Chiang, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Chun-Wei Chang, Takeshi Miki, Hao Ye, Sami Souissi, Rita Adrian, et al. Causal networks of phytoplankton diversity and biomass are modulated by environmental context. Nature Communications, 2022, 13 (1), 10.1038/s41467-022-28761-3. insu-03664870

HAL Id: insu-03664870

https://insu.hal.science/insu-03664870

Submitted on 11 May 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

ARTICLE

https://doi.org/10.1038/s41467-022-28761-3

<figure><img src="image_4.png" /></figure>

Check for updates

OPEN

# Causal networks of phytoplankton diversity and biomass are modulated by environmental context

Chun-Wei Chang ID 1,2, Takeshi Miki ID 3,4,5, Hao Ye ID 6, Sami Souissi ID 7, Rita Adrian 8,9, Orlane Anneville 10, Helen Agasild 11, Syuhei Ban ID 12, Yaron Be'eri-Shlevin 13, Yin-Ru Chiang 14, Heidrun Feuchtmayr ID 15, Gideon Gal ID 13, Satoshi Ichise 16, Maiko Kagami ID 17,18, Michio Kumagai ID 16,19, Xin Liu ID 12, Shin-Ichiro S. Matsuzaki 20, Marina M. Manca ID 21, Peeter Nõges ID 11, Roberta Piscia ID 21, Michela Rogora 21, Fuh-Kwo Shiah 2,4, Stephen J. Thackeray ID 15, Claire E. Widdicombe 22, Jiunn-Tzong Wu 14, Tamar Zohary 13 & Chih-hao Hsieh ID 1,2,4,23✉

Untangling causal links and feedbacks among biodiversity, ecosystem functioning, and environmental factors is challenging due to their complex and context-dependent interactions (e.g., a nutrient-dependent relationship between diversity and biomass). Consequently, studies that only consider separable, unidirectional effects can produce divergent conclusions and equivocal ecological implications. To address this complexity, we use empirical dynamic modeling to assemble causal networks for 19 natural aquatic ecosystems (N24°~N58°) and quantified strengths of feedbacks among phytoplankton diversity, phytoplankton biomass, and environmental factors. Through a cross-system comparison, we identify macroecological patterns; in more diverse, oligotrophic ecosystems, biodiversity effects are more important than environmental effects (nutrients and temperature) as drivers of biomass. Furthermore, feedback strengths vary with productivity. In warm, productive systems, strong nitrate-mediated feedbacks usually prevail, whereas there are strong, phosphate-mediated feedbacks in cold, less productive systems. Our findings, based on recovered feedbacks, highlight the importance of a network view in future ecosystem management.
Since ecosystems were first described as delicate feedback systems by Tansley¹, feedback has been a recurring theme among global-scale ecosystem studies²⁻⁵. A feedback is defined as a directed and connected path of causal interactions that ends on the originating node (i.e., a “cycle” in network terminology). As feedbacks have a critical role in dynamical systems over long-term observations⁶,⁷, feedbacks have been suggested as an important consideration for elucidating interactions between biodiversity and ecosystem functioning (BDEF)⁸ and how they regulate natural systems⁹. However, biodiversity and ecosystem functioning are only two components of a larger interconnected network with many causal links and feedbacks¹⁰⁻¹² among a multitude of environmental factors, including nutrient availability¹³ and temperature¹⁴. Ignoring these feedbacks and the role of environmental factors can complicate the interpretation of BDEF relationships. For example, impacts of plant diversity loss on plant biomass cannot be precisely evaluated if feedbacks among plant diversity, biomass, and environment are overlooked⁸.

A holistic view of the causal network associated with biodiversity, integrating posited relationships from previous studies, is shown (Fig. 1). Here, each arrow represents a simple causal interaction (e.g., BD → EF depicts biodiversity effects on ecosystem functioning). In addition to pairwise feedback between biodiversity and ecosystem functioning (BD ↔ EF), more complex triangular feedbacks exist when including endogenous environmental factors (e.g., nutrients) that can affect and be affected by organisms¹⁵. Considering this complexity in natural ecosystems, we need to incorporate feedbacks into the current research framework of biodiversity.

Although causal networks among biodiversity, ecosystem functioning, and the environment have been discussed¹⁰⁻¹²,

<figure><img src="image_5.png" /><figcaption>**Fig. 1** Conceptual causal network of biodiversity and ecosystem functioning (BDEF) relationships. Environmental variables in the causal network (a) can be exogenous (e.g., climate), which influence BD and EF, or endogenous (e.g., nutrients), which influence and can be influenced by BD and EF. Whereas endogenous factors can affect and be affected by organisms⁵, exogenous factors, such as precipitation and temperature, can only affect ecosystems (organisms do not influence precipitation and temperature on the scales considered in a majority of ecological studies, e.g., daily, monthly, or annual scales) and therefore cannot be included in feedbacks²⁶. The causal network can be decomposed into modules (b): (i) individual causal links (e.g., L1-L8), (ii) pairwise feedbacks (e.g., L3-L4), and more complex (iii) triangular feedbacks. Triangular feedbacks connecting biodiversity, ecosystem functioning, and nutrients (gray triangle in (a)) can be classified based on direction: BD → EF → Nutrients (Type I feedback: L4-L5-L6) and EF → BD → Nutrients (Type II feedback: L3-L7-L8).</figcaption></figure>

quantification of these networks is yet to be fully realized. Because manipulating multiple interdependent processes is infeasible¹²,¹⁶, experimental studies often examine individual interactions in isolation¹⁷. As previously noted, however, these studies might ignore the context of other relevant factors¹⁸. Alternatively, empirical reconstructions of causal networks from observational data have often led to equivocal results¹⁹,²⁰. The use of linear statistical methods may be an explanation: correlation and regression-based approaches assume static relationships among variables and are not designed to investigate interdependent feedbacks that produce time-varying interactions observed in natural systems²¹⁻²³. Consequently, there is a lack of network-based approaches that can identify relative contributions of various drivers to a focal ecosystem process and explore conditions under which their contributions might change. For instance, diversity is recognized as the strongest determinant for ecosystem functioning²⁴ in experimental systems (i.e., L4 in Fig. 1). However, in natural systems, debate remains over whether the effect of diversity on ecosystem functioning is stronger than exogenous (L2 in Fig. 1) or endogenous drivers²⁵ (L8 in Fig. 1), both of which affect organisms³, although only endogenous drivers can be affected by organisms and involved in feedbacks²⁶. Similarly, the consensus is lacking about relative contributions among causal determinants for species diversity²⁷⁻²⁹ (L1, L3, and L6 in Fig. 1). Therefore, lack of proper quantification of individual causal links could fail to identify the most critical drivers in ecosystems and to reconstruct complex feedbacks under various environmental contexts.

Better quantification of more complex feedbacks (e.g., pairwise feedbacks and triangular feedbacks in Fig. 1b) is essential to predict responses of ecosystems to external perturbations³⁰. Whereas major pairwise feedbacks indicate interactions in a network that potentially amplify or dampen external perturbations, quantifying triangular feedbacks can provide a better mechanistic understanding of these feedbacks and enable more accurate predictions for how changes in one variable will propagate to other parts of the network, producing more comprehensive, nonadditive impacts on ecosystems than biodiversity effects alone⁸,³¹. Thus, quantifying these complex network modules may have important management implications, in shifting the focus from managing individual state variables toward managing integrated ecosystem processes³².

In this study, we used a combination of nonlinear time series methods, convergent cross-mapping (CCM)³³, and cross-system network analysis to elucidate the role of diversity and ecosystem functioning in natural aquatic ecosystems. CCM is a method rooted in the theory of dynamical systems³⁴ that enables the detection of causation between time series variables³³. It is noteworthy that our current knowledge on interactions among biodiversity, ecosystem functioning, and environmental contexts is mainly derived from terrestrial ecosystems rather than natural aquatic ecosystems¹²,³⁵, despite phytoplankton accounting for >50% of global primary production³⁶. Although effects of phytoplankton species diversity on biomass and resource use efficiency have been examined³⁵,³⁷, the importance of diversity effects and more complex interaction modules remain unclear in natural aquatic systems. Therefore, we employed CCM³³ to assemble causal networks for 19 sites (Supplementary Fig. S1) among 16 globally distributed ecosystems (“Methods” and Supplementary Fig. S1), representing various types of aquatic ecosystems with various morphometrics and trophic states (from oligotrophic to eutrophic systems presented in Supplementary Table S1). Our datasets consisted of long-term monthly measurements (16–41 years) of phytoplankton data, with biodiversity and ecosystem function operationalized as phytoplankton species richness and community biomass (using chlorophyll-a
concentration as a proxy)⁹, respectively. In addition, environmental variables, including concentrations of nitrate (NO₃) and phosphate (PO₄) (endogenous factors), and water temperature (an exogenous factor) were also involved in the reconstruction of causal networks for each system (see Supplementary Fig. S2 for an example, and similarly such a causal network comprised of the same variables as Supplementary Fig. S2 was reconstructed for each of the 19 sites). Based on the reconstructed networks, we aimed to understand biodiversity in aquatic ecosystems by addressing the following questions:

1. Under what conditions are phytoplankton diversity effects on ecosystem functioning stronger than the effects of environmental drivers?

2. What is the strongest causal determinant for species diversity?

3. What are the most effective pathways through which changes in diversity propagate to other parts of the network, and feedback on themselves?

4. Are there any emerging macroecological patterns explaining how causal links, pairwise feedbacks, and triangular feedbacks vary along large-scale environmental gradients?

To explicitly answer these questions, we performed cross-system comparisons on the reconstructed causal networks to evaluate: (i) the relative importance among causal links affecting phytoplankton biomass; (ii) the relative importance among causal links affecting phytoplankton diversity; (iii) the relative strengths of more complex feedbacks involving biodiversity; and iv) how the strengths of the network modules investigated in (i)–(iii) vary with environmental characteristics.

Overall, our analysis presents quantitative causal networks consisting of causal interactions and feedbacks among phytoplankton diversity, biomass, and environmental drivers and reveals how the network varies along large-scale environmental gradients. Our results indicate that phytoplankton diversity is a more important determinant to phytoplankton biomass than other environmental factors (nutrients and temperature) in more diverse, oligotrophic ecosystems; nutrients play an important role in determining the dynamics of phytoplankton diversity in most systems. In addition, strong nitrate-diversity-biomass feedbacks prevail in warm, productive systems; while strong phosphate-diversity-biomass feedbacks prevail in cold, less productive systems. These findings anticipate the response of aquatic ecosystems to environmental changes from a holistic network view.

## Results and discussion

**Quantification of causal networks.** We first compared the relative strengths of causal links across systems (Supplementary Fig. S3). Phytoplankton species richness was the major controlling factor for phytoplankton biomass (significant in 16 of 19 sites, Fig. 2a) in these diverse aquatic systems, consistent with experimental studies¹⁷. However, the averaged linkage strength for this effect was not significantly different from that of NO₃ (i.e., BD → EF vs. NO₃ → EF; permutation test *P* = 0.501), highlighting that nitrogen availability was equally important in affecting phytoplankton biomass in natural systems.

In the opposite direction, phytoplankton biomass was a significant driver of phytoplankton species richness in most ecosystems (15 of 19 sites, Fig. 2b). However, NO₃ more often had a stronger effect, appearing as the most important driver in 11 of 19 sites compared to phytoplankton biomass (4 of 19 sites) (Fig. 2b). Although the difference in effect strength was not significant (permutation test, *P* = 0.162), these results implicated nitrogen availability as an essential determinant affecting both phytoplankton diversity and biomass. As a sensitivity test, we also

examined the effects of Shannon diversity. The results suggest that the importance of nutrients is robust to the use of other diversity indexes (e.g., Shannon diversity in Supplementary Fig. S4), although the causal effects from phytoplankton biomass became relatively more important compared to biomass effects on species richness (Fig. 2b). Based on these findings, we inferred that processes influencing nutrients (e.g., external loadings and internal cycling³⁸) need to be considered when investigating aquatic biodiversity. Changes in those processes (e.g., climatic³⁹ or anthropogenic⁴⁰ driven nutrient changes) may indeed substantially impact phytoplankton biodiversity, and subsequent ecosystem functioning.

The importance of NO₃ uncovered in our analyses might not be a counter-intuitive result, as many systems analyzed in this study were P-rich. For instance, the average phosphate concentration was 57.5 and 41.7 µgP/L for Lake Mendota (Me) and Lake Monona (Mo) (Supplementary Table S1), respectively. In addition, there were also high total phosphorus (TP) concentrations in shallow lake systems, e.g., average TP was 106.1, 112.5, and 126.4 µgP/L in Lake Inba (Ib), Lake Kasumigaura (Ks), and Müggelsee (Mu), respectively. Phosphorus was not always a limiting factor in eutrophic and mesotrophic systems, e.g., Lake Kasumigaura⁴¹ and Lake Geneva (Gv)⁴². In addition, nitrogen was deficient and limited cyanobacteria bloom in Müggelsee (Mu)⁴³. Nonetheless, we cannot exclude the possibility of colimitation⁴⁴ in N and P and the possibility that P availability also depends on N⁴⁵, which warrants further investigation.

Apart from nutrients and temperature, the causal effects of other important drivers on phytoplankton biomass and diversity were also examined, though not in all 19 systems due to data limitation. The causal effects of physical environmental factors, such as irradiance and water column stability, were presented in Supplementary Fig. S5; the results indicated that the quantified causal strengths on average were not as strong as the effects of diversity and nutrients. Moreover, the effects of consumers (e.g., zooplankton), which have been suggested as important drivers affecting species diversity of phytoplankton communities⁴⁶, were also examined. Based on our analysis of zooplankton, the causal effects of herbivorous crustaceans on phytoplankton biomass and diversity were significant in most of the analyzed systems. However, these effects were on average not as strong as the effects of phytoplankton diversity and nutrients, respectively (Supplementary Fig. S6). Nonetheless, these findings were not generalized to all 19 systems due to a lack of complete datasets as shown in Supplementary Table S3, and thus warrant more detailed investigation in future studies.

In addition to individual causal effects, we investigated feedbacks across systems. Pairwise feedbacks (e.g., BD ↔ EF and NO₃ ↔ EF) were common (Fig. 2c). However, the averaged linkage strength was often stronger in one direction when involving BD (Fig. 3). Specifically, the average strength of BD → EF was stronger than for the opposite direction of EF → BD (permutation test *P* = 0.015); BD → EF was stronger than EF → BD in 14 of the 19 systems (Fig. 3). In addition, biodiversity effects on nutrients (BD → NO₃ and BD → PO₄) were also stronger than their reversed effects (NO₃ → BD and PO₄ → BD) in 12 and 13 systems, respectively. In comparison, the interactions between nutrients and productivity were more symmetrical: nutrient effects on biomass (NO₃ → EF and PO₄ → EF) were stronger than biomass effects on nutrients (EF → NO₃ and EF → PO₄) in only 9 and 8 of 19 systems, respectively. These results supported the previous findings⁸ that biodiversity effects more often operate at short-term scales, which makes effects more observable in our monthly-scale analyses than feedback effects on diversity, which are expected to occur on a more prolonged timescale, e.g., through slowly changing nutrient
<figure><img src="image_6.png" /><figcaption>**Fig. 2** Relative strengths of various modules. Standardized linkage strengths of causal variables affecting (a) phytoplankton biomass and (b) species richness (here, BD) and loop weights for various types of (c) pairwise feedbacks and (d) triangular feedbacks. All statistics were calculated from the 19 independent sites (n = 19) and depicted as joint violins and box plots to present the empirical distribution that labels the maxima and minima at the top and bottom of the violins, respectively, and shows 25, 50, and 75% quantiles in the boxes with whiskers presenting at most 1.5 * interquartile range. The two numbers within the parentheses (S; R₁) above each violin plot report the number of significant results in CCM (S; labeled blue) and the number of systems in which a particular module had the greatest strength (i.e., rank 1; R1; labeled red). Source data are provided as a Source Data file.</figcaption></figure>

cycling³¹ or decomposition⁴⁷. Nevertheless, the timescale dependence of causal interactions in ecosystem networks is a topic that needs further study.

Subsequently, we quantified the strengths of pairwise feedbacks as the geometric mean of the linkage strengths in each direction, following a previous study⁹ (see more details in Methods). Among these feedbacks (Fig. 2c and Supplementary Fig. S7), $BD \leftrightarrow NO_3$ had the highest median and average strength (0.78 and 0.68, respectively) across systems. However, strengths of $BD \leftrightarrow NO_3$ were highly variable among systems (large interquartile range in Fig. 2c), and thus were only significant in 11 of 19 systems, compared to $BD \leftrightarrow EF$ (15 of 19 systems). These findings reinforced the importance of nutrients as key determinants for aquatic biodiversity and implied that nutrient effects are context-dependent. In other words, $BD \leftrightarrow NO_3$ was less common than $BD \leftrightarrow EF$ across systems, despite its stronger average strength. The prevalence of $BD \leftrightarrow EF$ indicated a need for more long-term experiments and process-based/theoretical modeling accounting for *bidirectional* interactions between diversity and biomass¹⁶, because bidirectional interactions and feedbacks may challenge our simple predictions for ecosystem dynamics, based on knowledge of unidirectional interactions³⁰.

Quantification of the causal network also allowed us to analyze triangular feedbacks. Within the conceptual framework of Fig. 1b, there are four kinds of triangular feedbacks involving biodiversity, ecosystem functioning, and either nitrate or phosphate (Type I: $BD \rightarrow EF \rightarrow NO_3$ and $BD \rightarrow EF \rightarrow PO_4$; Type II: $EF \rightarrow BD \rightarrow NO_3$ and $EF \rightarrow BD \rightarrow PO_4$). There was at least one significant triangular feedback in 14 of 19 sites (Fig. 2d). More specifically, $NO_3$-associated feedbacks (Type I-N and Type II-N) were usually stronger than $PO_4$-associated feedbacks (Type I-P and Type II-P)

(Fig. 2d), although the difference in strength among the four types of feedbacks was not significant (Fig. 2d; Kruskal-Wallis test, *P* = 0.59). The dominance of $NO_3$-associated feedbacks in our study was attributed to many of the sites being marine and eutrophic lakes, which are likely to be N-limited due to an imbalance in external loadings⁴⁸ or strong denitrification⁴⁹. Among both $NO_3$- and $PO_4$-associated feedbacks, there were no significant differences in strength between Type I and Type II feedbacks (Supplementary Fig. S7), suggesting that biodiversity can directly influence biomass (Type I), as well as through a pathway that involves endogenous nutrient variables (Type II) and eventually feeds back on itself.

**Causal networks under environmental contexts.** Our empirical analyses revealed state dependency of the causal links and feedbacks among biodiversity, biomass, and environmental factors in natural systems; that is, their strengths were highly dependent on the state of other variables. Based on a cross-system comparison (Methods), strengths of individual links (e.g., $BD \rightarrow EF$), pairwise feedbacks (e.g., $BD \leftrightarrow EF$), and triangular feedbacks (e.g., $BD \rightarrow EF \rightarrow NO_3 \rightarrow BD$) varied systematically, depending on environmental characteristics (Fig. 4 and Supplementary Fig. S8). Ecosystems with higher species diversity (long-term average species richness) and lower average $PO_4$ concentrations had stronger $BD \rightarrow EF$ links (Fig. 4a; correlation coefficient *r* = 0.600 and −0.513; *P* = 0.007 and 0.025 for species diversity and $PO_4$, respectively). These results were further confirmed by stepwise regression, indicating that the ecosystems characterized by higher diversity, lower average temperature, and oligotrophic conditions had stronger $BD \rightarrow$
<figure><img src="image_7.png" /><figcaption>**Fig. 3 Directional bias in pairwise feedbacks.** The difference in standardized linkage strengths between the two directions was computed for each pairwise feedback and depicted as joint violin and box plots. All statistics were calculated from the 19 independent sites (n = 19) and depicted as joint violins and box plots to present the empirical distribution that labels the maxima and minima at the top and bottom of the violins, respectively, and shows 25, 50, and 75% quantiles in the boxes with whiskers presenting at most 1.5 * interquartile range. The number above the plot indicates the number of systems with a positive difference in linkage strength. For example, BD → EF was stronger than its feedback, EF → BD, in 14 of the systems. In general, the strength of diversity effects (BD → EF, BD → NO₃, BD → PO₄) was usually stronger than feedback effects (EF → BD, NO₃ → BD, PO₄ → BD). Source data are provided as a Source Data file.</figcaption></figure>

EF (best-fit regression model: $BD → EF$ strength = 0.663 + 0.171*BD - 0.139*T - 0.096*PO₄; $F_{3, 15}$ = 9.958 and $P < 0.001$). In contrast, temperature and PO₄ effects on phytoplankton biomass (i.e., $T → EF$ and PO₄ → EF) were negatively associated with long-term average species diversity, but positively associated with average PO₄ (Fig. 4a). Therefore, we inferred that phytoplankton biomass in P-rich systems was more sensitive to warming (due to strong $T → EF$). This synergistic effect of warming and eutrophication on biomass has been reported in other aquatic ecosystems⁵⁰; in this study, this synergistic effect was weaker when species diversity was higher and BD → EF was stronger (Fig. 4a). Perhaps greater diversity and its effects mitigate adverse impacts of global warming⁹, although warming may also weaken biodiversity effects on ecosystem functioning due to strong interspecific competitions under high temperatures⁵¹.

In the opposite direction, stronger EF → BD was associated with lower temperature and PO₄ concentrations, as well as shallower depths (Fig. 4b; best-fit regression model: $EF → BD$ strength = 0.572 - 0.170*PO₄ - 0.101*Depth - 0.096*T; $F_{3, 15}$ = 9.800 and $P < 0.001$). In shallower systems, which are better mixed and less vertically heterogeneous, impacts of species competition on diversity may be more influential⁵². In contrast, the effects of the most important driver on diversity, NO₃ → BD, exhibited an opposite response to temperature (Fig. 4b), suggesting that the dominant determinants of aquatic biodiversity varied along a temperature gradient.

Water temperature and phytoplankton biomass (Chla as a proxy) also critically determine strengths of various pairwise feedbacks. Stronger PO₄-mediated feedbacks (BD ↔ PO₄ and $EF ↔ PO₄$ in Fig. 4c) were usually more associated with cold

($r$ = -0.247 and -0.329, respectively) and less productive systems ($r$ = -0.421 and -0.527, respectively); this contrasted with $BD ↔ NO₃$, which was more associated with warm (Fig. 4c; $r$ = 0.503) and productive environments ($r$ = 0.571). This finding was consistent with the notion that N is more often a limiting element for phytoplankton growth in warm, tropical/subtropical systems than in cold, temperate systems⁵³. However, $EF ↔ PO₄$ and $BD ↔ EF$ had no clear relationship with temperature ($r$ = 0.167 and -0.177, respectively) or biomass ($r$ = 0.284 and 0.163, respectively). Thus, we speculated that climate warming will shift aquatic ecosystems towards stronger coupling in biodiversity-NO₃ feedback than biodiversity-biomass feedback.

Our analyses improved understanding of how more complex regulations varied with environmental characteristics. In our study, strengths of triangular feedbacks, regardless of direction, were positively associated with high diversity environments. Of the four triangular feedbacks, both type I-N feedbacks ($BD → EF → NO₃$) as well as the type II-N feedback ($EF → BD → NO₃$), were positively associated with high-diversity environments (Fig. 4d, $r$ = 0.531 and 0.561, respectively). Therefore, we inferred that the dynamics of biodiversity, biomass, and endogenous variables were tightly coupled in high-diversity systems; that is, dynamics of one component were more responsive to changes in other parts of the feedback. Our findings contrasted with the prevailing view that ecosystem functioning is insensitive to changes in diversity at high levels of diversity (i.e., the redundancy concept¹⁸), prompting a further investigation to clarify the role of biodiversity in regulating ecosystem dynamics⁹.

Interestingly, triangular feedbacks in different directions (i.e., Type I versus II) had distinct responses to biomass levels. The strength of Type I feedbacks had no statistical association with biomass, implying that diversity effects on biomass (i.e., $BD → EF$) can propagate to nutrients (i.e., $EF → nutrients$) and then to diversity itself (i.e., $nutrients → BD$), irrespective of biomass levels. In contrast, Type II feedbacks had associations with biomass; the latter was positively associated with strengths of Type II-NO₃ feedbacks ($r$ = 0.567; $P$ = 0.011), but negatively associated with strengths of Type II-PO₄ feedbacks ($r$ = -0.430; $P$ = 0.066). This highlighted the importance of considering linkage directionality when studying these regulatory feedbacks under various environmental contexts. For example, even when the same components were considered, two causal interactions in the opposite direction (e.g., $BD → EF$ or $EF → BD$) responded to environmental factors differently (e.g., Fig. 4a, b, respectively).

Our network-based approach enabled us to describe how responses of triangular feedbacks to environmental changes differed from responses of individual links or pairwise feedbacks. Different regulatory feedbacks were associated with distinctive environmental characteristics (Fig. 4d) and not necessarily similar to that of individual links involved. For instance, $EF → BD$-driven triangular feedbacks (Type II) were associated with average phytoplankton biomass (Fig. 4d), but phytoplankton biomass was associated with neither the $EF → BD$ nor $BD ↔ EF$, individually. Indeed, the cross-system patterns of triangular feedbacks were statistically distinguishable from that of pairwise feedbacks (Supplementary Fig. S9), implying unique responses of complex feedbacks to environmental gradients. Therefore, predicting ecosystem responses to environmental changes is challenging, even if responses of individual links or of pairwise feedbacks can be elucidated, because quantification of individual causal links in isolation might fail to recover more complex network modules. As feedbacks and other network modules⁵⁴ are critical for stabilizing/destabilizing ecosystem dynamics⁵⁵, it is becoming apparent that studying interdependencies among key ecological processes from a holistic network view is needed for predicting ecosystem dynamics under environmental changes.
<figure><img src="image_8.png" /><figcaption>b Drivers of phytoplankton species richness: X→BD</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>c Pairwise feedback X→Y</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>d Triangular feedback</figcaption></figure>

<figure><img src="image_11.png" /></figure>

**Fig. 4** Strengths of causal interactions and feedbacks vary along environmental gradients of aquatic ecosystems. Multivariate ordination illustrating associations between long-term averages of environmental factors (blue) and quantitative network modules, including the strength of links affecting: **a** phytoplankton biomass (X → EF; red) and **b** species richness (X → BD; red), and the loop weight of **c** pairwise feedbacks and **d** triangular feedbacks. Here, water temperature was abbreviated by T. Because biomass was selected throughout the analyses, we used color scales (green) to visualize data points according to their log₁₀ biomass values. A similar figure based on the color scales of water temperature was presented in Supplementary Fig. S8. Black text indicates abbreviated site names (see Supplementary Fig. S1). **a** Stronger BD → EF was observed in environments characterized by low PO₄ but high diversity. In comparison, T → EF and PO₄ → EF were associated with high PO₄ but low diversity, and NO₃ → EF was associated with high temperature. **b** Stronger EF → BD was observed in the environments characterized by intermediate level of phytoplankton biomass. In comparison, stronger T → BD and PO₄ → BD were usually present in systems with a lower phytoplankton biomass. However, NO₃ → BD was associated with high-temperature environments. By including feedbacks of diversity and biomass on endogenous environments, we determined that: **c** stronger diversity-biomass feedbacks were more often observed in phosphate-poor systems. Diversity-nitrate feedbacks (BD ↔ NO₃) were stronger in high-temperature environments; diversity-phosphate (BD ↔ PO₄) feedback were weaker in low biomass environments. **d** Stronger phosphate-associated feedbacks (Type I-P and Type II-P) were more often observed in systems characterized by lower temperature, phytoplankton biomass, and nitrate (NO₃). In contrast, stronger nitrate-associated feedbacks (Type I-N and Type II-N) were more often observed in environments with higher temperature (especially for Type II-N), biomass, and species diversity. Selected environmental variables significantly explained the multivariate ordinations (permutation test *P* = 0.004, 0.029, <0.001, and 0.001 for Panels **a** to **d**, respectively). Source data are provided as a Source Data file.

**Sensitivity analysis using composition-converted biomass measure.** Our conclusions were robust to the use of alternative biomass measure. Specifically, the main findings (Figs. 2 and 4) based on the analysis of phytoplankton biomass inferred by Chla were qualitatively similar with findings based on composition-converted biomass (Supplementary Figs. S10 and 11). Although the relationship between Chla and true phytoplankton biomass varies with environmental conditions (e.g., light), it remains an effective functional index inferring phytoplankton stock with

more emphasis on photosynthesis capacity (i.e., biomass of photosynthetic machines). In contrast, composition-converted biomass, though has similar meaning with overall phytoplankton biomass, contains high uncertainty by assuming species-specific conversion factors, especially when the measurement of individual cells size was lacking. In addition, these conversion factors were determined by various geometrical models, which differ among systems; this is in contrast to the standard chemical approach used in determining Chla, which makes Chla more
suitable to reveal cross-system variations in causal strengths along environmental gradients (Fig. 4 and Supplementary Fig. S11e–h). However, Chla integrates all kinds of photoautotrophs that might not be fully included in counting data (e.g., picoplankton). Thus, investigating the diversity effects of more complete phytoplankton groups requires novel techniques (e.g., metagenomics³⁶). Nonetheless, it remains an open question about how the presented causal links and feedbacks change when considering various types of functional indices and diversity measures.

**Caveats for the reconstruction of causal networks in natural phytoplankton communities.** Several issues warrant further studies in aquatic ecosystem networks involving phytoplankton diversity and biomass. Firstly, the number of marine sites was limited, hindering comparisons of marine versus freshwater systems. Furthermore, our analyses based on CCM cannot access the sign of feedbacks (i.e., positive or negative), although it is known that the sign is important in determining the response of feedbacks to external perturbations (e.g., amplified or dampened). Although methods to estimate the sign of interactions were proposed (e.g., S-map²²,⁵⁷), the robustness of these methods has not been thoroughly examined⁵⁸. Lastly, due to limitations of data availability, our analysis only quantified causal strength across systems at a consensus monthly scale, acknowledging that state-space reconstruction methods (e.g., CCM) are scale-dependent⁵⁹, e.g., one causal driver dominated monthly might not necessarily dominate at other time scales. Therefore, exploring causal feedbacks at other time scales needs further investigations by including more datasets with high temporal resolution and long-duration monitoring.

**Final remarks.** Our findings bridged two popular and contrasting research directions: whereas many studies consider diversity effects on ecosystem functioning⁸,⁶⁰, other studies aim to identify determinants of species diversity⁶¹, which can be traced to Hutchinson's seminal question about species coexistence⁶². Our findings highlighted that these two ecological processes are interdependent and embedded in a complex network. Thus, these ecological processes should not be investigated in isolation¹⁰,¹²,³⁰, but instead be examined in integrated feedbacks, especially for rapid turnover systems (e.g., plankton or microbial communities in aquatic systems) in which feedbacks from biomass to diversity can operate quickly through light shading²⁹ or exploitation on nutrients⁶³. It is noteworthy that our proposed methodological framework can be applied to explore in more detail causal feedbacks or paths if precise mechanistic measures (e.g., nutrient recycling rate rather than nutrient stock) can be monitored over time. For instance, biodiversity was suggested to influence ecosystem functioning via species complementarity or selection effects. Measuring complementarity or selection effects is available for experimental data⁶⁴, but remains a challenging task for observational data. Thus, the incorporation of these detailed mechanistic measures in the causal networks is an important future research topic. More comprehensive surveys are required in future ecological monitoring to improve our understanding of causal mechanisms embedded in causal networks.

Our analyses quantified diversity-associated causal networks among various natural aquatic ecosystems. The selected long-term datasets from various aquatic ecosystems represented a reasonable parallel to long-term biodiversity experiments conducted in terrestrial grassland ecosystems⁶⁵. Revealing the dominance of diversity effects on biomass in these systems (Fig. 2) was enabled by using methods for nonlinear dynamical systems, in lieu of linear statistical analyses¹⁹,²⁵. Indeed, when

linear analyses were applied in our datasets (Supplementary Fig. S12), the importance of diversity effects on biomass (BD → EF) and nitrate effects on diversity (NO₃ → BD) were not clearly identified as that shown in the nonlinear CCM analysis (Fig. 2a, b).

Through cross-system comparison, the strength of causal relationships associated with species diversity varied with nutrient levels and time-averaged mean levels of diversity in aquatic systems (Fig. 4a). Associations with environmental gradients were also present in other network modules (i.e., feedbacks; Fig. 4c, d). These statistical associations revealed the macroecological relationships of how the strength of biodiversity effects and related feedbacks varied with environmental gradients. Moreover, the unveiled macroecological patterns also improved our understanding of how causes and effects of biodiversity are modulated by biotic and abiotic contexts. Although the proposed relationships need to be examined through more long-term experiments, our quantitative and empirical framework for constructing causal networks provided a foundation for better predicting the consequences of biodiversity loss across ecosystems.

## Methods

**Data.** Phytoplankton species composition and environmental data were compiled from 16 aquatic ecosystems with a total of 19 long-term monitoring sites, spanning a large range of freshwater and marine types, from shallow to deep and from oligotrophic to eutrophic, as follows (Supplementary Table S1 and Supplementary Fig. S1): (1) Lake Biwa, Japan, 1978–2010; (2) Feitsui Reservoir, Taiwan, 1986–2017; (3) Lake Geneva, France/Switzerland, 1974–2014; (4) Lake Inba, Japan, 1986–2016 (including three stations in disparate lake basins); (5) Lake Kasumigaura, Japan, 1978–2009 (including two stations in distinct lake basins); (6) Lake Kinneret, Israel, 1996–2012; (7) Lake Maggiore, Italy, 1997–2015; (8) Lake Mendota, USA, 1995–2012; (9) Lake Monona, USA, 1995–2011; (10) Müggelsee, Germany, 1994–2013; (11) Narragansett Bay, USA, 1999–2014; (12) Lake Oneida, USA, 1975–1995; (13) Shin River, Japan, 1986–2016; (14) Lake Võrtsjärv, Estonia, 2001–2016; (15) Station L4, Western English Channel, England, 1992–2009; and (16) Windermere, England, 1993–2010 (South basin).

For all systems, there were five types of variables: (1) phytoplankton species richness (number of species recorded in a sample); (2) chlorophyll-α concentration as a measure of phytoplankton biomass and ecosystem function, a widely used proxy of algal biomass in the BDEF literature¹⁷,⁶⁶; (3) phosphate concentration (PO₄); (4) nitrate concentration (NO₃); and (5) water temperature. Phytoplankton samples were identified to the finest taxonomical level (generally species level if possible) and enumerated under an optical microscope, based on counting methods summarized in Supplementary Table S2. The counting methods used were similar (e.g., Utermöhl⁶⁷ method and relevant approaches). Based on composition data, species richness was derived and defined as the number of species present in the phytoplankton community. In systems with depth-resolved measurements, data were depth-integrated averages in the euphotic zone; otherwise, measurements were from surface layer samples.

We additionally examined other environmental factors that are known to be important to phytoplankton, including water column stability, irradiance, and zooplankton abundance, although these variables were not measured in all systems (Supplementary Table S3). Water column stability was calculated as maximal Brunt-Väisälä frequency from temperature vertical profile data as an index of water column stability; irradiance data was compiled from in situ measurements of weather or buoy stations near the sampling sites. For zooplankton analysis, we compiled composition and density data (individual/L) of crustacean zooplankton in 11 sites (Supplementary Table S3) based on microscopic counting. Specifically, we investigated grazing effects of the following three zooplankton categories: (i) herbivorous cladocerans excluding predatory taxa, such as *Bythotrephes* spp. and *Leptodora* spp.; (ii) herbivorous copepods including all calanoids and naupliar stages of cyclopoids; and (iii) herbivorous crustaceans including both herbivorous cladocerans and herbivorous copepods (i.e., i+ii). It is noteworthy that our zooplankton analysis was based on density instead of biomass data because zooplankton length measurements, which are required for converting individual counts to biomass data, were absent in 5 of the 11 analyzed sites.

**Data treatment.** For consistency, monthly time series were generated by averaging over observations if sampling occurred on a finer timescale. Although such compilation potentially causes some inconsistency in smoothing temporal fluctuations of time series data among systems with various sampling frequencies, it was necessary because our methods based on state-space reconstruction require time series data at equal intervals, dictating the temporal scale of analysis. In our case, the monthly resolution is the only consensus that can be applied to all time series datasets and the monthly average is the most representative measure at this scale.
Nonetheless, causal strengths estimated by CCM analysis were robust to this data averaging according to our comparisons using eight stations where regular and frequent sampling (i.e., sampling frequency higher than monthly) were available (Supplementary Fig. S13). Overall, our data compilation yielded 5554 data points for each variable across the 19 sites (Supplementary Table S1). To ensure stationarity, we removed the long-term linear trend from each time series by using the residuals from a linear regression against time⁹. We accounted for seasonality by scaling against the mean and standard deviation of values occurring in the same month⁹, $D_{m,v}(t_i) = (O(t_i) - \mu_{month i})/\sigma_{month i}$, where $\mu_{month i}$ is the monthly mean, $\sigma_{month i}$ is the monthly standard deviation for each of 12 months, $O(t_i)$ is the original time series, and $D_{m,v}(t_i)$ is the deseasonalized time series, $i = 1, 2, ..., 12$. Finally, each time series was re-scaled to zero long-term mean and unit variance⁶⁸.

**Convergent cross-mapping analysis.** Causal networks were reconstructed among phytoplankton species richness, phytoplankton biomass, and the environment with a method specifically designed for quantifying causality in nonlinear dynamical ecosystems, convergent cross-mapping (CCM)³³. In that regard, CCM is a causality analysis based on Takens' theorem for dynamical systems³⁴,⁶⁹, which infers the causal relationship among variables from their empirical time series. CCM³³ tests for causality between pairs of time series by measuring the extent to which the historical record of an effect variable, X, can reliably estimate states of a causal variable, $Y^{33,69,70}$. Cross-map skill, the quantification of this measure, is defined as the correlation coefficient $\rho$ between estimated states of the causal variable and actual observations⁷¹. CCM is based on information recovery (i.e., effect variables contain encoded information on causal variables), instead of predictive ability (using causal variables to predict future values of effect variables, e.g., Granger's causality). The essential ideas of CCM are summarized in the following brief animations: tinyurl.com/EDM-intro. The aforementioned deseasoning procedure reduces detection of false positives caused by 'dynamical synchronization'⁷²,⁷³ under strong seasonality. Apart from seasonality, dynamical synchronization can also occur when interactions between two variables are very strong³³; nevertheless, very strong interactions are of less concern here because most interactions in real ecosystems are weak to moderate⁷⁴. A modeling study also indicated that CCM was robust against moderate noise from process and observational errors⁷⁵.

Several limitations in applying CCM analysis need to be acknowledged. First, CCM is based on lagged coordinated embedding in which each embedded variable needs to be lagged by a fixed time interval that determines the timescale of CCM analysis. For example, time series analyzed in our study were integrated to the monthly scale as stated in "Data treatment" section. Multi-scale analysis is possible only when data points are measured very regularly across various time scales (e.g., from weekly to monthly). Second, as required in many time series analysis, CCM analysis requires time series data being stationary⁷⁶. Otherwise, CCM likely produces false-positive findings (e.g., caused by strong seasonality); note however, we have removed seasonality in analysis in this work. Third, a time series including too many zero values (or other constant values) is not suitable for CCM analysis (as a general statistical issue in any time series analysis). This is because embedding such a time series potentially produces many zero vectors, which violates the general assumption of EDM that assumes a one-to-one mapping between each embedded vector and the vector on dynamical manifold³⁴ (i.e., zero vector can map to many possibilities on manifold). Thus, the embedded zero vectors need to be excluded or separated from the prediction set⁷³.

CCM analysis accounts for influences of confounding variables implicitly. Specifically, CCM incorporates influences of confounding variables using lagged embeddings, e.g., $(X_{t-1}, X_{t-2}, ...)$, which have accounted for historical effects of other variables in lagged terms, even if those variables were unobserved or difficult to identify. As such, CCM does not require identifying or ruling out influences of confounding variables in order to quantify causations between two variables, and thus can be applied in more general dynamical systems⁶⁸. In addition, CCM is a nonparametric approach, free from assumptions of the specific form of quantitative relationships between causal variables. Although this makes CCM difficult to explore quantitative features, e.g., the minimal number of species required to maintain 80% levels of ecosystem function, it provides high flexibility to infer causations in nonlinear dynamical systems. Such flexibility is important for inferring nonlinear dynamical systems, because quantitative relationships between any two dynamical variables could change, depending on the varying state of other state variables⁷⁷ or environmental contexts. For example, linear associations between two variables will appear then disappear or change sign—so-called mirage correlations³³, making methods based on modeling static, parametric relationships difficult to correctly identify causations⁹.

In this study, the best embedding dimension (E) used in CCM was determined by testing a range from 2 to 20 and selecting the value that optimized the hindcast cross-mapping in which $X(t)$ projected one-step backward to $Y(t-1)$; this avoids overfitting of E when X and Y are unrelated time series⁷³. Note that E can vary for each variable pair: the E selected for X cross-mapping Y can be different from the E selected for Y cross-mapping X or X cross-mapping Z.

The possibility of lagged causal effects⁷⁸ was explored in this study. Specifically, we tested for the causal influence of Y on X by cross-mapping between $X(t + k)$ and

$Y(t)$ using time lags ($k$) of 0, 1, 2, or 3 months (corresponding to the timescale of phytoplankton dynamics) and selecting the mapping with the highest cross-map skill. To determine the convergence in cross-mapping, we followed the procedure in Sugihara et al.³³ and computed the cross-map skill for subsamples of $X(t)$ with varying library lengths (L). Here, the minimal library length, $L_0$, is equal to the embedding dimension, and the maximal library length, $L_{max}$, is equal to the length of the whole time series. To test the convergence of CCM, we applied two statistical criteria. First, we tested whether there was a significant monotonic increasing trend in cross-map skill, $\rho(L)$, using Kendall's $\tau$ test. Next, we tested the significance of improvement in cross-map skill using Fisher's $\Delta\rho Z$ test and compared the cross-map skill for the maximal library size ($\rho(L_{max})$) against the cross-map skill at the minimal library length ($\rho(L_0)$).

**Quantification of causal interaction and loop weight.** As in previous studies⁷⁵,⁷⁶, the strength of causal interaction was quantified based on the cross-mapping skill at the maximal library length, $\rho(L_{max})$. That is, stronger causal effects result in convergence to high cross-map skill⁹,³³,⁷⁶. For instance, a strong causal effect of species richness on phytoplankton biomass revealed by CCM indicated that dynamics of phytoplankton biomass (magnitude or variability) responded strongly to changing species richness. To address systematic differences in cross-map skill among study sites, which may arise due to differences in noise or time series length, causal strength⁹ was standardized. The standardized linkage strength (SLS) was calculated by dividing linkage strength (LS) by the maximum within each system: $SLS = LS/max(LS)$; thus, SLS varied between 0 and 1 and indicated the relative importance with respect to the strongest causal link within the system. It is noteworthy that causal networks were constructed and standardized separately for each system; i.e., it was not assumed that each system had equivalent dynamics and belonged to the same attractor.

Pairwise and triangular feedbacks were quantified using Neutel's loop weight⁵⁵, the geometric mean of SLS for all links within a given feedback. In a pairwise feedback ($X \leftrightarrow Y$), the loop weight is the geometric mean of the SLS in both directions (i.e., $X \rightarrow Y$ and $Y \rightarrow X$). We classified two types of triangular feedbacks, based on the directionality of the involved interactions: "Type I", richness→biomass→nutrients→richness and "Type II", biomass→richness→nutrients→biomass. The Type I feedbacks occur in the direction that includes biodiversity effects on ecosystem function (BD → EF), whereas Type II feedbacks are in the opposite direction and include EF → BD. With two types of nutrients stocks (phosphorus P or nitrogen N), there are a total of four triangular feedbacks (I-N, I-P, II-N, and II-P). To determine the uncertainty of our estimates in causal strength and loop weight, we calculated their standard errors using resampling method that reconstructed sampling distributions from 500 random samples of embedded data points with replacement.

**Linking strengths of links and feedbacks with ecosystem characteristics.**

Multivariate redundancy analysis (RDA)⁷⁹ was used to illustrate how the strength of individual links and feedbacks varied in association with ecosystem characteristics: depth and area of the study site, as well as long-term averages of species richness, temperature, phosphate, and nitrate. Specifically, we conducted RDA analyses for: (i) causal effects on phytoplankton biomass (Fig. 4a); (ii) causal effects on species richness (Fig. 4b); and (iii) pairwise and triangular feedbacks (Fig. 4c, d), respectively. For each multivariate RDA ordination, we constructed the biplot using the first two RDA scores to demonstrate how the prevalence of various links or feedbacks were statistically associated with various environmental characteristics. The use of RDA instead of CCA was justified based on our analysis on coenoclines (Supplementary Fig. S14), with more linear coenoclines and a short gradient length (<3)⁸⁰. RDA significance was evaluated using a permutation test⁷⁹. All permutation tests performed in this study were based on the null distribution generated from 10,000 random permutations. To further support the RDA results, correlation analysis was performed between the strength of key network modules and environmental characteristics and tested significance using a permutation test. These analyses were not intended to examine causation between long-term environmental characteristics and strength of network modules, but rather to provide a picture describing under what environmental conditions a module of interest (e.g., BD → EF) prevailed.

**Computation.** All analyses were done with R (ver. 4.0.3). The CCM analyses and the multivariate RDA analysis were implemented using the rEDM⁸¹ and vegan⁸² packages, respectively.

**Reporting summary.** Further information on research design is available in the Nature Research Reporting Summary linked to this article.

**Data availability**

Raw time series datasets from all research sites are available, on request, through the paths listed in Supplementary Table S2 due to various data use policy. Source data for all figures are provided with this paper and available from GitHub online repository, https://github.com/biozoo/Chang_etal_2022_SI_CausalFeedback⁸³. Source data are provided with this paper.
## Code availability

Documentation of all analytical procedures provided as R codes are available from Github, https://github.com/biozoo/Chang_etal_2022_SI_CausalFeedback⁸³.

Received: 29 December 2020; Accepted: 11 February 2022;
Published online: 03 March 2022

## References

1. Tansley, A. G. The use and abuse of vegetational concepts and terms. *Ecology* **16**, 284–307 (1935).

2. Lovelock, J. E. & Margulis, L. Atmospheric homeostasis by and for the biosphere: the Gaia hypothesis. *Tellus* **26**, 2–10 (1974).

3. Berendse, F. Litter decomposability—a neglected component of plant fitness. *J. Ecol.* **82**, 187–190 (1994).

4. Watson, A. J. & Lovelock, J. E. Biological homeostasis of the global environment: the parable of Daisyworld. *Tellus B: Chem. Phys. Meteorol.* **35**, 284–289 (1983).

5. Miki, T., Ushio, M., Fukui, S. & Kondoh, M. Functional diversity of microbial decomposers facilitates plant coexistence in a plant–microbe–soil feedback model. *Proc. Natl Acad. Sci. USA* **107**, 14251–14256 (2010).

6. Charlson, R. J., Lovelock, J. E., Andreae, M. O. & Warren, S. G. Oceanic phytoplankton, atmospheric sulphur, cloud albedo and climate. *Nature* **326**, 655–661 (1987).

7. Mangan, S. A. et al. Negative plant–soil feedback predicts tree-species relative abundance in a tropical forest. *Nature* **466**, 752–755 (2010).

8. Tilman, D., Isbell, F. & Cowles, J. M. Biodiversity and ecosystem functioning. *Annu. Rev. Ecol. Evol. Syst.* **45**, 471–493 (2014).

9. Chang, C.-W. et al. Long-term warming destabilizes aquatic ecosystems through weakening biodiversity-mediated causal networks. *Glob. Change Biol.* **26**, 6413–6423 (2020).

10. Chapin, F. S. III et al. Consequences of changing biodiversity. *Nature* **405**, 234–242 (2000).

11. Dee, L. E. et al. Operationalizing network theory for ecosystem service assessments. *Trends Ecol. Evolution* **32**, 118–130 (2017).

12. Snelgrove, P. V. R., Thrush, S. F., Wall, D. H. & Norkko, A. Real world biodiversity–ecosystem functioning: a seafloor perspective. *Trends Ecol. Evolution* **29**, 398–405 (2014).

13. Martin, J. H. & Fitzwater, S. E. Iron deficiency limits phytoplankton growth in the north-east Pacific subarctic. *Nature* **331**, 341–343 (1988).

14. Behrenfeld, M. J. et al. Climate-driven trends in contemporary ocean productivity. *Nature* **444**, 752 (2006).

15. Sunda, W. Feedback interactions between trace metal nutrients and phytoplankton in the ocean. *Front. Microbiol.* **3**, 204 (2012).

16. Cardinale, B. J., Weis, J. J., Forbes, A. E., Tilmon, K. J. & Ives, A. R. Biodiversity as both a cause and consequence of resource availability: a study of reciprocal causality in a predator–prey system. *J. Anim. Ecol.* **75**, 497–505 (2006).

17. Cardinale, B. J. Biodiversity improves water quality through niche partitioning. *Nature* **472**, 86–89 (2011).

18. Fetzer, I. et al. The extent of functional redundancy changes as species’ roles shift in different environments. *Proc. Natl Acad. Sci. USA* **112**, 14888–14893 (2015).

19. Adler, P. B. et al. Productivity is a poor predictor of plant species richness. *Science* **333**, 1750–1753 (2011).

20. Waide, R. B. et al. The relationship between productivity and species richness. *Annu. Rev. Ecol. Syst.* **30**, 257–300 (1999).

21. Moran, P. A. P. The statistical analysis of the Canadian lynx cycle. *Aust. J. Zool.* **1**, 291–298 (1953).

22. Deyle, E. R., May, R. M., Munch, S. B. & Sugihara, G. Tracking and forecasting ecosystem interactions in real time. *Proc. R. Soc. B* **283**, 20152258 (2016).

23. Ushio, M. et al. Fluctuating interaction network and time-varying stability of anatural fish community. *Nature* **554**, 360–363 (2018).

24. Hooper, D. U. et al. A global synthesis reveals biodiversity loss as a major driver of ecosystem change. *Nature* **486**, 105–108 (2012).

25. Mittelbach, G. G. et al. What is the observed relationship between species richness and productivity? *Ecology* **82**, 2381–2396 (2001).

26. Enoki, T. et al. Progress in the 21st century: a roadmap for the Ecological Society of Japan. *Ecol. Res.* **29**, 357–368 (2014).

27. Harley, C. D. G. Climate change, keystone predation, and biodiversity loss. *Science* **334**, 1124–1127 (2011).

28. Hsu, S. B., Hubbell, S. & Waltman, P. A mathematical theory for single-nutrient competition in continuous cultures of micro-organisms. *SIAM J. Appl. Math.* **32**, 366–383 (1977).

29. Huisman, J., Jonker, R. R., Zonneveld, C. & Weissing, F. J. Competition for light between phytoplankton species: experimental tests of mechanistic theory. *Ecology* **80**, 211–222 (1999).

30. Hughes, R. A., Byrnes, J. E., Kimbro, D. L. & Stachowicz, J. J. Reciprocal relationships and potential feedbacks between biodiversity and disturbance. *Ecol. Lett.* **10**, 849–864 (2007).

31. Mueller, K. E., Tilman, D., Fornara, D. A. & Hobbie, S. E. Root depth distribution and the diversity–productivity relationship in a long-term grassland experiment. *Ecology* **94**, 787–793 (2013).

32. Pikitch, E. K. et al. Ecosystem-based fishery management. *Science* **305**, 346–347 (2004).

33. Sugihara, G. et al. Detecting causality in complex ecosystems. *Science* **338**, 496–500 (2012).

34. Takens, F. Dynamical systems and turbulence. *Lect. Notes Math.* **898**, 366–381 (1981).

35. Ye, L. et al. Functional diversity promotes phytoplankton resource use efficiency. *J. Ecol.* **107**, 2353–2363 (2019).

36. Behrenfeld, M. J. et al. Biotic primary production during an ENSO transition. *Science* **291**, 2594–2597 (2001).

37. Duffy, J. E., Godwin, C. M. & Cardinale, B. J. Biodiversity effects in the wild are common and as strong as key drivers of productivity. *Nature* **549**, 261–264 (2017).

38. Anneville, O., Ginot, V. & Angeli, N. Restoration of Lake Geneva: expected versus observed responses of phytoplankton to decreases in phosphorus. *Lakes Reserv. Res. Manag* **7**, 67–80 (2002).

39. North, R. P., North, R. L., Livingstone, D. M., Köster, O. & Kipfer, R. Long-term changes in hypoxia and soluble reactive phosphorus in the hypolimnion of a large temperate lake: consequences of a climate regime shift. *Glob. Change Biol.* **20**, 811–823 (2014).

40. Smith, V. H. & Schindler, D. W. Eutrophication science: where do we go from here? *Trends Ecol. Evolution* **24**, 201–207 (2009).

41. Matsuzaki, S.-iS., Suzuki, K., Kadoya, T., Nakagawa, M. & Takamura, N. Bottom-up linkages between primary production, zooplankton, and fish in a shallow, hypereutrophic lake. *Ecology* **99**, 2025–2036 (2018).

42. Anneville, O. et al. Temporal mapping of phytoplankton assemblages in Lake Geneva: annual and interannual changes in their patterns of succession. *Limnol. Oceanogr.* **47**, 1355–1366 (2002).

43. Shatwell, T. & Köhler, J. Decreased nitrogen loading controls summer cyanobacterial blooms without promoting nitrogen-fixing taxa: Long-term response of a shallow lake. *Limnol. Oceanogr.* **64**, S166–S178 (2019).

44. Paerl, H. W. et al. It takes two to tango: when and where dual nutrient (N & P) reductions are needed to protect lakes and downstream ecosystems. *Environ. Sci. Technol.* **50**, 10805–10813 (2016).

45. Schindler, D. W. The dilemma of controlling cultural eutrophication of lakes. *L. R. Soc. B* **279**, 4322–4333 (2012).

46. Lewandowska, A. M., Hillebrand, H., Lengfellner, K. & Sommer, U. Temperature effects on phytoplankton diversity—The zooplankton link. *J. Sea Res.* **85**, 359–364 (2010).

47. Gessner, M. O. et al. Diversity meets decomposition. *Trends Ecol. Evol.* **25**, 372–380 (2010).

48. Kosten, S. et al. Lake and watershed characteristics rather than climate influence nutrient limitation in shallow lakes. *Ecol. Appl.* **19**, 1791–1804 (2009).

49. Seitzinger, S. P. & Nixon, S. W. Eutrophication and the rate of denitrification and N₂O production in coastal marine sediments. *Limnol. Oceanogr.* **30**, 1332–1339 (1985).

50. Hsieh, C.-h et al. Phytoplankton community reorganization driven by eutrophication and warming in Lake Biwa. *Aquat. Sci.* **72**, 467–483 (2010).

51. Parain, E. C., Rohr, R. P., Gray, S. M. & Bersier, L.-F. Increased temperature disrupts the biodiversity–ecosystem functioning relationship. *Am. Naturalist* **193**, 227–239 (2019).

52. Holt, R. D. Spatial heterogeneity, indirect interactions, and the coexistence of prey species. *Am. Naturalist* **124**, 377–406 (1984).

53. Lewis, W. M. Tropical limnology. *Annu. Rev. Ecol. Syst.* **18**, 159–184 (1987).

54. Allesina, S. & Levine, J. M. A competitive network theory of species diversity. *Proc. Natl Acad. Sci. USA* **108**, 5638–5642 (2011).

55. Neutel, A.-M., Heesterbeek, J. A. P. & de Ruiter, P. C. Stability in real food webs: weak links in long loops. *Science* **296**, 1120–1123 (2002).

56. Cuvelier, M. L. et al. Targeted metagenomics and ecology of globally important uncultured eukaryotic phytoplankton. *Proc. Natl Acad. Sci.* **107**, 14679–14684 (2010).

57. Chang, C.-W. et al. Reconstructing large interaction networks from empirical time series data. *Ecol. Lett.* **00**, 1–12 (2021).

58. Yu, Z., Gan, Z., Huang, H., Zhu, Y. & Meng, F. Regularized S-Map reveals varying bacterial interactions. *Appl. Environ. Microbiol.* **86**, e01615-20 (2020).

59. Gibson, J. F., Doyne Farmer, J., Casdagli, M. & Eubank, S. An analytic approach to practical state space reconstruction. *Phys. D: Nonlinear Phenom.* **57**, 1–30 (1992).

60. Cardinale, B. J., Duffy, J. E., Gonzalez, A., Hooper, D. U. & Perrings, C. Biodiversity loss and its impact on humanity. *Nature* **486**, 59–67 (2012).
61. Rosenzweig, M. L. *Species Diversity in Space and Time* (Cambridge University Press, 1995).

62. Hutchinson, G. E. Homage to Santa Rosalia or why are there so many kinds of animals? *Am. Naturalist* **93**, 145–159 (1959).

63. Tilman, D. *Resource Competition and Community Structure*, Vol. 17 (Princeton University Press, 1982).

64. Loreau, M. & Hector, A. Partitioning selection and complementarity in biodiversity experiments. *Nature* **412**, 72–76 (2001).

65. Jochum, M. et al. The results of biodiversity-ecosystem functioning experiments are realistic. *Nat. Ecol. Evolution* **4**, 1485–1494 (2020).

66. Lewandowska, A. M. et al. The influence of balanced and imbalanced resource supply on biodiversity-functioning relationship across ecosystems. *Philos. Transact. Royal Soc. B: Biol. Sci.* **371**, 20150283 (2016).

67. Utermöhl, H. Zur Vervollkommnung der quantitativen Phytoplankton-Methodik: Mit 1 Tabelle und 15 abbildungen im Text und auf 1 Tafel. *Int. Ver. f.ür. theoretische und Angew. Limnologie: Mitteilungen* **9**, 1–38 (1958).

68. Ye, H. et al. Equation-free mechanistic ecosystem forecasting using empirical dynamic modeling. *Proc. Natl Acad. Sci. USA* **112**, E1569–E1576 (2015).

69. Sauer, T., Yorke, J. & Casdagli, M. Embedology. *J. Stat. Phys.* **65**, 579–616 (1991).

70. Cummins, B., Gedeon, T. & Spendlove, K. On the efficacy of state space reconstruction methods in determining causality. *SIAM J. Appl. Dyn. Syst.* **14**, 335–381 (2015).

71. Clark, A. T. et al. Spatial convergent cross mapping to detect causal relationships from short time series. *Ecology* **96**, 1174–1181 (2015).

72. Sugihara, G., Deyle, E. R. & Ye, H. Reply to Baskerville and Cobey: misconceptions about causation with synchrony and seasonal drivers. *Proc. Natl Acad. Sci. USA* **114**, E2272–E2274 (2017).

73. Deyle, E. R., Maher, M. C., Hernandez, R. D., Basu, S. & Sugihara, G. Global environmental drivers of influenza. *Proc. Natl Acad. Sci. USA* **113**, 13081–13086 (2016).

74. McCann, K., Hastings, A. & Huxel, G. R. Weak trophic interactions and the balance of nature. *Nature* **395**, 794 (1998).

75. BozorgMagham, A. E., Motesharrei, S., Penny, S. G. & Kalnay, E. Causality analysis: identifying the leading element in a coupled dynamical system. *PLoS ONE* **10**, e0131226 (2015).

76. Chang, C.-W., Ushio, M. & Hsieh, C.-H. Empirical dynamic modeling for beginners. *Ecol. Res.* **32**, 785–796 (2017).

77. Clark, T. J. & Luis, A. D. Nonlinear population dynamics are ubiquitous in animals. *Nat. Ecol. Evolution* **4**, 75–81 (2020).

78. van Nes, E. H. et al. Causal feedbacks in climate change. *Nat. Clim. Change* **5**, 445–448 (2015).

79. Legendre, P. & Legendre, L. F. *Numerical Ecology*. Vol. 24 (Elsevier, 2012).

80. Smilauer, P. & Lepš, J. *Multivariate Analysis of Ecological Data Using CANOCO* 5 (Cambridge University Press, 2014).

81. Ye, H. et al. rEDM: Applications of empirical dynamic modeling (EDM) from time series. R package version 1.2.3. https://CRAN.R-project.org/package=rEDM (2020).

82. Oksanen, J. et al. vegan: Community ecology package. R package version 2.5-7. https://CRAN.R-project.org/package=vegan (2020).

83. Chang, C.-W. et al. Source data and R code from: causal networks of phytoplankton diversity and biomass are modulated by environmental context. Github online repository https://github.com/biozoo/Chang_etal_2021 SI CausalFeedback (2022).

## Acknowledgements

This work was supported by the National Center for Theoretical Sciences, National Taiwan University, Academia Sinica, Foundation for the Advancement of Outstanding Scholarship, and the Ministry of Science and Technology, Taiwan (to C.H.H.). Data for Oneida Lake were collected with support from Cornell University's Brown Endowment, New York State Department of Environmental Conservation, and United States Department of Agriculture, National Institute of Food and Agriculture, Hatch Project 0226747. Research on Lake Maggiore is within the framework of the LTER Italian and European networks, site 'IT-08 Southern Alpine lakes' and funded by the International Commission for the Protection of Swiss-Italian Waters (CIPAIS). Data for Lake Geneva were contributed by the Observatory of alpine LAKes (OLA), © SOERE OLA-IS, AnaEE-France, INRAE of Thonon-les-Bains, CIPEL. Data for Lake Vörtsjärv were provided by

Estonian Environment Agency and by the Centre for Limnology at Estonian University of Life Sciences funded by the Estonian Research Council grants PRG 1167 and PRG709, and by the European Union's Horizon 2020 research and innovation programme under grant agreement No. 95196. Data for Lake Kinneret were collected by the Kinneret Limnological Laboratory, Israel Oceanographic & Limnological Research, and funded by the Israel Water Authority, with phytoplankton counts made by T. Fishbein, chlorophyll determined by Y. Yacobi, physical data collected by Y. Lechinsky, chemical analyses conducted by Mekorot Water Company, Watershed Unit, under oversight by A. Nishri, and database management services provided by M. Shlichter. Data from Lake Inba were provided by Chiba Prefecture. Data from Lake Kasumigaura were provided by the National Institute for Environmental Studies (NIES). Data from Lake Biwa were provided by Shiga Prefecture. Data for Müggelsee were provided by the Leibniz-Institute of Freshwater Ecology and Inland Fisheries within their long-term research programme. Data collection at Windermere was supported by Natural Environment Research Council award number NE/R016429/1 as part of the UK-SCAPE program delivering National Capability. Data for Station L4, Western Channel Observatory were collected by Plymouth Marine Laboratory as part of the UK's Natural Environment Research Council's National Capability CLASS Programme grant number NE/R015953/1, and is a contribution to Theme1.3 - Biological Dynamics. This is a contribution of GEISHA project, which was jointly supported by the French Foundation for Research on Biodiversity (FRB) through its synthesis center CESAB (http://www.cesab.org/) and the John Wesley Powell Center for Analysis and Synthesis (https://powellcenter.usgs.gov/). Comments from B. Kraemer, M. Kondoh, M. Ushio, and P. J. Ke on an earlier draft and English editing by J. Kastelic improved the manuscript.

## Author contributions

C.W.C., C.H.H., and T.M. conceived the research idea. C.W.C. analyzed the data with help from C.H.H., H.Y., and S.S. R.A., O.A., S.J.T., Y.B., Y.R.C., H.F., S.I., Ma.K., Mi.K., S.I.M., P.N., M.R., F.K.S., C.E.W., J.T.W., T.Z., H.A., G.G., S.B., X.L., M.M.M., and R.P. collected the data. C.W.C., C.H.H., T.M., and H.Y. wrote the draft of the manuscript with critical comments from all co-authors.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-022-28761-3.

**Correspondence** and requests for materials should be addressed to Chih-hao Hsieh.

**Peer review information** *Nature Communications* thanks Donald L. DeAngelis, Nathalie Niquil and the other anonymous reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_12.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022

<footnote n="1">¹National Center for Theoretical Sciences, Taipei 10617, Taiwan. ²Research Center for Environmental Changes, Academia Sinica, Taipei 11529, Taiwan. ³Faculty of Advanced Science and Technology, Ryukoku University, Otsu, Shiga 520-2194, Japan. ⁴Institute of Oceanography, National Taiwan University, Taipei 10617, Taiwan. ⁵Center for Biodiversity Science, Ryukoku University, Otsu, Shiga 520-2194, Japan. ⁶Health Science</footnote>
Center Libraries, University of Florida, Gainesville, FL 32611, USA. ⁷Univ. Lille, CNRS, Univ, Littoral Côte D'Opale, IRD, UMR 8187, LOG-
Laboratoire D'Océanologie et de Géosciences, Station Marine de Wimereux, F- 59000 Lille, France. ⁸Leibniz Institute of Freshwater Ecology and
Inland Fisheries, IGB, 12587 Berlin, Germany. ⁹Freie Universität Berlin, Department of Biology, Chemistry and Pharmacy, 14195 Berlin, Germany.
¹⁰National Research Institute for Agriculture, Food and Environment (INRAE), CARRTEL, Université Savoie Mont Blanc, 74200 Thonon les
Bains, France. ¹¹Centre for Limnology, Institute of Agricultural and Environmental Sciences, Estonian University of Life Sciences, Kreutzwaldi 5D,
51014 Tartu, Estonia. ¹²Department of Ecosystem Studies, School of Environmental Science, The University of Shiga Prefecture, Hikone 522-8533
Shiga, Japan. ¹³Kinneret Limnological Laboratory, Israel Oceanographic & Limnological Research, P.O. Box 447, 14950 Migdal, Israel.

¹⁰National Research Institute for Agriculture, Food and Environment (INRAE), CARRTEL, Université Savoie Mont Blanc, 74200 Thonon les Bains, France. ¹¹Centre for Limnology, Institute of Agricultural and Environmental Sciences, Estonian University of Life Sciences, Kreutzwaldi 5D, 51014 Tartu, Estonia. ¹²Department of Ecosystem Studies, School of Environmental Science, The University of Shiga Prefecture, Hikone 522-8533 Shiga, Japan. ¹³Kinneret Limnological Laboratory, Israel Oceanographic & Limnological Research, P.O. Box 447, 14950 Migdal, Israel.

¹⁴Biodiversity Research Center, Academia Sinica, Taipei 11529, Taiwan. ¹⁵UK Centre for Ecology & Hydrology, Lancaster Environment Centre, Library Avenue, Bailrigg, Lancaster, Lancashire LA1 4AP, UK. ¹⁶Lake Biwa Environmental Research Institute, Otsu 520-0022, Japan. ¹⁷Faculty of Environment and Information Sciences, Yokohama National University, Yokohama 240-8502 Kanagawa, Japan. ¹⁸Department of Environmental Science, Faculty of Science, Toho University, Funabashi, Chiba 274-8510, Japan. ¹⁹Research Center for Lake Biwa & Environmental Innovation, Ritsumeikan University, Kusatsu 525-0058 Shiga, Japan. ²⁰Biodiversity Division, National Institute for Environmental Studies, 16-2 Onogawa, Tsukuba, Ibaraki 305-8506, Japan. ²¹CNR Water Research Institute (IRSA), L.go Tonolli 50, 28922 Verbania, Pallanza, Italy. ²²Plymouth Marine Laboratory, Prospect Place, West Hoe, Plymouth PL1 3DH, UK. ²³Institute of Ecology and Evolutionary Biology, Department of Life Science, National Taiwan University, Taipei 10617, Taiwan. ²⁴email: chsieh@ntu.edu.tw