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.