Methods and Materials

The study ecosystem

We studied a model of anoxic-oxic regime shifts in aquatic ecosystems in response to changes in oxygen diffusivity. Depletion of oxygen has detrimental effects on the survival of many aquatic organisms, leading to changes in community composition and to a decline in biodiversity (Vaquer-Sunyer & Duarte 2008; Hughes et al. 2015). Hypoxic conditions can arise naturally when thermal stratification in summer leads to reduced oxygen diffusivity to deep-water habitats (Wetzel 2001). Anthropogenic increases of nutrients and temperature have enhanced the frequency and duration of low-oxygen episodes (Diaz & Rosenberg 2008; Scavia et al. 2014; Jane et al. 2021). Moreover, severe hypoxia/anoxia can be difficult to reverse because of feedbacks between organisms and their environment (Diaz & Rosenberg 2008; Conley et al. 2009; Bushet al. 2017).
Using a mathematical model, Bush et al. (2017) described anoxic-oxic regime shifts resulting from interactions among microbial organisms. The model ecosystem contains three functional groups of bacteria (Figure 2; cyanobacteria (CB), sulfate-reducing bacteria (SB), and phototrophic sulfur bacteria (PB)), four chemical substrates (oxygen, phosphorous, reduced sulfur, and oxidized sulfur), and four types of flows/interactions (production, consumption, inhibition, and diffusion); there is also abiotic oxidization of the reduced sulfur.
Critical for the presence of alternative stable states in this system is the mutual inhibition of cyanobacteria and sulfur bacteria due to their intolerance of the chemical substrate that the other produces: cyanobacteria inhibit sulfur bacteria by producing oxygen, whereas sulfate-reducing bacteria inhibit cyanobacteria by producing reduced sulfur, e.g. sulfide (Figure 2). Bush et al. (2017) simulated this ecosystem with a set of ordinary differential equations (ODEs) and demonstrated that the state of the system depends on the rate at which oxygen diffuses into (or out of) the system: the oxygen diffusivity. The system is oxic and dominated by cyanobacteria when oxygen diffusivity is high, but anoxic and dominated by the two groups of sulfur bacteria when oxygen diffusivity is low. At intermediate oxygen diffusivity, the system is bistable, i.e. either oxic or anoxic depending on whether cyanobacteria or sulfur bacteria dominated in the past (Bush et al. 2017).
The Bush et al. model (2017) assumes zero diversity within each of the three functional groups. There is, however, empirical evidence of intra- and inter-specific variation in tolerance (Rolfe et al.1978; Knoll & Bauld 1989; Miller & Bebout 2004; Ramel et al.2015). Furthermore, experimental evolution studies demonstrated considerable and rapid increases in tolerance of bacteria (Martín-Clemente et al. 2019; Schoeffler et al. 2019).

The model and our extension of it

Bush et al. (2017) has an accessible and complete description of the model of the ecosystem, including rate equations, parameter values, and initial conditions. Our implementation of the model is documented in the microxanox R package (R Core Team 2021; Krug & Petchey 2022), which has two vignettes (a user guide and a partial reproduction of the analyses in Bush et al. (2017)). The R code to run, analyze, and visualize the numerical experiments of our study is available in a repository of supplementary content (Petchey et al. 2022), which also includes all variables and parameter values.
We extended the Bush et al. model (2017) by modelling multiple strains per functional group. Each of the strains of a functional group has its own parameter set and state variable. For example, a system with nine strains per functional group has 27 state variables for the strains, and four for the chemical substrates. Other than having potentially different parameter values for maximum growth rate and tolerance, the dynamics of the strains in a functional group were described by the same equations. We adjusted the chemical substrate ODEs to correctly account for the number of strains, by summing their substrate consumption and production.

Creating and varying within functional group diversity

We generated diversity within functional groups by creating multiple strains that differed in maximum growth rate and in tolerance to the respective inhibiting substrate (sulfide tolerance of cyanobacteria, and oxygen tolerance of sulfate-reducing bacteria and phototrophic sulfur bacteria). We assumed a trade-off between the two traits, in line with empirical studies that report a negative correlation between maximum growth rate and resistance to environmental stress in microbial organisms (Ferenci 2016). For simplicity, we used a linear trade-off between log10-transformed trait values, with among strain variation distributed evenly on a log-scale (see Supplementary Report Section 3 for more information). In the main text, we report results for nine strains in each functional group. Results were relatively robust to changes in the number of strains, although some effects of strain number were observed at high levels of trait diversity (see Supplementary Report Section 7).
Functional diversity has multiple components, and we manipulated only one of these, namely the range of trait space represented, which has been termed functional richness (Villéger et al. 2008). We compared 20 different amounts of trait diversity ranging from no variation, i.e. all nine strains had the reference trait value (Table 1), to greatest trait variation with the range of trait values according to Table 1. The maximum trait range was selected to ensure that the region of bistability was mostly within the range of oxygen diffusivity values used in Bush et al. (2017) when all functional groups contained the maximum amount of variation. The other 18 trait diversity levels were evenly distributed between the no diversity scenario and the greatest amount of diversity. Twenty diversity levels gave sufficient resolution and required reasonable computation time to simulate. A detailed account is in Supplementary Report Section 3.
We also manipulated the combination of functional groups that contained diversity. Specifically, we examined the effect of adding diversity to only one group (either CB, or SB, or PB), simultaneously to two functional groups (CB-SB, CB-PB, and SB-PB), and simultaneously to all three functional groups (CB-SB-PB). This led to 7 combinations * 20 diversity levels, but with seven of these being the same zero diversity in all functional groups, there were 134 diversity scenarios in total.
We have confirmed that when there is no within functional group strain variation, the dynamics of the model are not affected by the number of strains in each of the functional groups. This is a necessary condition for such a model (Moisset de Espanés et al. 2021), and gives confidence that our methods and their implementation are correct and robust.

Stable state finding

To investigate the effect of trait diversity on tipping point positions, we parameterized the system according to each of the diversity scenarios described above, and then ran two simulations, one with a stepwise increase and one with a stepwise decrease of oxygen diffusivity. This technique was used to find stable states, rather than to simulate the effects of environmental change per se. At each of 300 values of oxygen diffusivity, ranging from 0 to -8 (log10h-1), we simulated the dynamics for 1,000,000 hours, which was sufficient for the position of the tipping points to stabilize, even though some strain dynamics had not yet completely stabilized (Supplementary Report Section 6.2). At the start of the two simulations, total initial abundance of each functional group was 105 cells, equally divided among the strains of each functional group. To avoid very low strain densities which would cause computational problems, we added 1 cell L-1 to every biological state variable every 1,000 hours (an alternative approach produced the same results, see Supplementary Report Section 4.3). At the end of each period of constant oxygen diffusivity, the system state was recorded and used as the initial conditions for the simulation at the next step of oxygen diffusivity. This approach for finding stable states differs from the method used by Bush et al. (2017). See the Supplementary Report (Section 4) for a detailed description of our rationale for using our approach, and also for evidence of qualitatively similar results produced by both approaches.
We calculated the tipping point position of the oxic to anoxic (anoxic-oxic) transition as the minimum (maximum) level of oxygen diffusivity at which the total density of sulfate-reducing bacteria was more than 1000 abundance units different between the increasing and decreasing oxygen diffusivity simulations. We did not calculate tipping point positions when coexistence patterns were atypical of the classical pattern shown in Figure 1a (e.g. when cyanobacteria coexisted with sulfur bacteria at low oxygen diffusivities; see Results and Supplementary Report Section 10.2 for details). We calculated the effect size of trait variation on tipping point position as the difference in oxygen diffusivity between a tipping point without trait variation and with a given level of trait variation (Fig. 1c-d). A positive effect size indicates increased resilience.