Supplementary Report for the article Functional diversity can facilitate the collapse of an undesirable ecosystem state

microxanox_version_here <- "0.9.1"
microxanox_version_path <- paste0("microxanox_v", microxanox_version_here)
if (packageVersion("microxanox") < package_version(microxanox_version_here)) {
  stop(paste0("microxanox version needs to be at least ", microxanox_version_here))
}
print(paste0("Version of microxanox: ", microxanox_version_here))

## [1] "Version of microxanox: 0.9.1"

1 Introduction

This document is the Supplementary Report for the report Functional diversity can facilitate the collapse of an undesirable ecosystem state.

2 Growth and inhibition functions

The rate of change in population densities is given by the per capita growth rate - per capita death rate, multiplied by population size (Bush et al 2017). Growth rate is the growth rate in absence of any inhibition multiplied by the amount of inhibition.

Growth rate in the absence of inhibition was modeled with the Monod equation with multiplicative kinetics when two substrates determine growth rate:

\(g_j(X, Y) = g_{max,j} (\frac{X}{K_{j,X} + X})(\frac{Y}{K_{j,Y} + Y})\)

\(g_{max,j}\) is the maximum specific growth rate of strain \(j\). \(X\) is the concentration of substrate 1, and \(Y\) is the concentration of substrate 2. \(K_{j,X}\) and \(K_{j,Y}\) are the half-saturation constants of species \(j\) on substrates \(X\) and \(Y\).

Inhibition of growth by a substrate was modeled with the Haldane equation:

\(h_j(X) = \frac{1}{1+ (X / H_{j,X})}\)

\(H_{j,X}\) can be thought of as a half-inhibition constant: the concentration of \(X\) at which the growth rate of species \(j\) is reduced by 50%. A higher half-inhibition constant therefore means higher tolerance of species \(j\) to substrate \(X\).

For further information on the functions used to model the system, please see the original publication by Bush et al 2017.

3 Creating among strain diversity

Here we provide additional description of the method used to create strain variation within functional groups. The following should be read only after the relevant text in the main report. We generated variation in maximum growth rate \(g_{max}\) and in tolerance to reduced sulfur or oxygen (i.e. in the half-inhibition constants \(H_{SR}\) or \(H_{O}\)) by varying the range of trait values around a reference trait value. To manipulate the amount of variation in a given trait and functional group, we introduced the meta-parameter \(Div\) subscripted by the trait and the functional group. Take for example the cyanobacteria functional group, a reference maximum growth rate parameter \(g_{max,CB}\) and a reference tolerance parameter \(H_{SR,CB}\). The amount of among strain variation in the cyanobacteria is controlled by the two meta-parameters \(Div_{g_{max,CB}}\) and \(Div_{H_{SR,CB}}\). To create a trade-off, the two meta-parameters were always of different sign. The among strain variation was calculated such that the growth rate of strain \(i = 1\) was a factor \(1/({2Div_{g_{max,CB}}})\) of the reference growth rate, and the growth rate of the \(i = n\) strain was a factor \(2Div_{g_{max,CB}}\) of the reference growth rate. The growth rate of the \(i = {2,...n-1}\) other strains was the reference growth rate multiplied by a factor that was equally distributed between the factor of the \(i = 1\) and \(i = n\) strain. This implementation of the variation and trade-off has the assumption of a linear trade-off of log-log transformed among-strain variation. Hence, for example, with the reference growth rate \(g_{max,CB} = 0.05\), among strain variation of \(Div_{g_{max,CB}} = 1\), and nine strains (\(n = 9\)), the growth rates of the nine strains are 0.025, 0.03, 0.035, 0.042, 0.05, 0.059, 0.071, 0.084, 0.1. Put another way, the growth rate parameter of the \(i\)-th of the \(n\) strains in a functional group was calculated as \(g_{max,i} = D_i * g_{max}\) where \(D_i\) is the \(i\)-th element of the set \(D = {2f(x)Div_{g_{max}} |x ∈ N, x ≤ n}\), where \(f(x) = (x − (n + 1)/2)/((n − 1)/2)\).

This procedure ensures that the range of trait values is independent of the number of strains. When there were only two strains, the two strains took the minimum and maximum traits values, as defined by the \(Div\) variable for that trait. When there were three strains, one had the minimum trait value, one the maximum, and one the mean. With nine strains, one had the minimum trait value, one the maximum, one the mean, and three were equally spaced (on log scale) between the minimum and the mean, and the remaining three between the maximum and the mean.

We show three examples of the created variation among strains within functional groups, one for nine strains (Figure 3.1), one for three (Figure 3.2), and one for two (Figure 3.3). In each of these three examples, three graphs show the diversity among strains in the maximum growth rate trait and the tolerance trait, plotted against each other and therefore also showing the trade-off between the two traits.

Notice that in all three cases, the range of trait variation represented by the strains is the same, even though the number of strains differs. This approach allowed us to independently manipulate trait diversity and number of strains. Also note that maximum growth rate of strains is evenly distributed on a log scale, even though it may appear at first site that they are distributed evenly on a linear scale (e.g. Figure 3.2).

All three examples have the following amounts of trait variation: \(Div_{g_{max,CB}}\) = 0.032, \(Div_{H_{SR,CB}}\) = -0.16, \(Div_{g_{max,SB}}\) = 0.095, \(Div_{H_{O,SB}}\) = -1.938, \(Div_{g_{max,PB}}\) = 0.095, \(Div_{H_{O,PB}}\) = -1.938. These are also the maximum amounts of trait variation explored in the simulations.

Figure 3.1: Trait variation with nine strains in each of the functional groups. In all panels colour indicates the functional group, darker shades show more tolerant-lower growth rate strains, lighter shades show less tolerant-higher growth rate strains. Tolerance is given as the concentration of the inhibiting substrate where the growth rate of the strain is reduced by 50%.

Figure 3.2: Trait variation with three strains in each of the functional groups. Note that the three displayed strains are the 1st, 5th, and 9th strains from the nine strain example shown in Figure 3.1. In all panels colour indicates the functional group, darker shades show more tolerant-lower growth rate strains, lighter shades show less tolerant-higher growth rate strains. Tolerance is given as the concentration of the inhibiting substrate where the growth rate of the strain is reduced by 50%.

Figure 3.3: Trait variation with two strains in each of the functional groups. Note that the two displayed strains are the 1st and 9th strains from the nine strain example shown in Figure 3.1. In all panels colour indicates the functional group, darker shades show more tolerant-lower growth rate strains, lighter shades show less tolerant-higher growth rate strains. Tolerance is given as the concentration of the inhibiting substrate where the growth rate of the strain is reduced by 50%.

A consequence of this variation and of the trade-offs is shown in Figure 3.4, which shows the relationship between realised growth rate and the concentration in the environment of the respective inhibiting substance, e.g. reduced sulfur for the cyanobacteria. In these graphs only the least and most tolerant strains are shown. The lines cross due to the trade-off. At high concentrations of the inhibiting substance the more tolerant strain has higher realised growth rate, while at lower concentrations it is the less tolerant strain that has the higher realised growth rate.

Figure 3.4: Response of realised growth rate to variation in concentration of the inhibiting substrate for two strains, namely the most and the least tolerant strain. As always colour indicates the functional group, darker shades show more tolerant-lower growth rate strains, lighter shades show less tolerant-higher growth rate strains.

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.

4 Method for finding stable states

4.1 General method

We now discuss two approaches for finding how the stable state of the system varies along the gradient of oxygen diffusivity when there exists the possibility of multiple stable states. In the main article we describe and report results of the method that we term the Temporal method.

Here we also describe and give results of what we term the Replication method. In this, we made separate simulations for each of many values of oxygen diffusivity, ran the simulations for long enough to ensure transients had passed, and recorded the state. We ran multiple simulations per value of oxygen diffusivity, and started them with different initial conditions (either low or high total cyanobacterial abundance), so as to check for presence of alternate stable states. This method was used by Bush et al 2017. We term it the replication method since it is achieved by having many independent replicates of the ecosystem, each with a different oxygen diffusivity, and also a separate replicate for each of the different chosen initial conditions.

As is shown in the figures below, the temporal and replication methods result in very similar regions of bistability.

4.1.1 Replication method, no diversity

Figure 4.1: Stable states found with the replication method. The left column of panels refers to the stable state when cyanobacteria density is at first high (favouring the oxic state). The right column refers to when cyanobacteria density is at first low (favouring the anoxic state). Points (which sometimes are so close as to appear as thick lines) show the stable state at the end of each period of constant oxygen diffusivity. Thin lines are for visualisation only, and join the points. Grey arrows show the oxygen diffusivity at which the system flips, and the direction of the flip. The length of the grey arrows is without meaning.

4.1.2 Temporal method, no diversity

Figure 4.2: Stable states found with the temporal method. The left column of panels refers to the stable state when oxygen diffusivity is at first high (favouring the oxic state). The right column refers to when oxygen diffusivity is at first low (favouring the anoxic state). Points (which sometimes are so close as to appear as thick lines) show the stable state at the end of each period of constant oxygen diffusivity. Thin lines are for visualisation only, and join the points. Grey arrows show the oxygen diffusivity at which the system flips, and the direction of the flip. The length of the grey arrows is without meaning.

4.1.3 Both methods overlaid, no diversity

Figure 4.3: Stable states of the two methods overlaid, the temporal method with a dashed line, and the replication method with a solid line.

4.2 Effect of diversity, replication method

4.2.1 Overall effects of diversity

The replication method shows qualitatively the same effects of functional diversity (Figure 4.4) as those revealed by the temporal method. Namely

• Introducing functional diversity in only the cyanobacteria increases the resilience of the oxic state and slightly decreases the resilience of the anoxic state.
• Introducing functional diversity in only the sulfate-reducing bacteria increases the resilience of the anoxic state and reduces the resilience of the oxic state.
• Introducing functional diversity in only the phototrophic sulfur bacteria reduces the resilience of the anoxic state and has no effect on the resilience of the oxic state.
• Introducing functional diversity in more than one functional group reduces, reverses, or does not change the effect of diversity in individual functional groups. For example, introducing functional diversity in all three functional groups reduces the resilience of the oxic state when trait variation is low but increases the resilience of the oxic state when trait variation is high.

Figure 4.4: Replication method results. Left column of panels: the effect of trait diversity (y-axis) on the positions of the tipping points (anoxic to oxic in blue, oxic to anoxic in red), which determine the extent of the region of bistability (green shaded region) in the oxygen diffusivity gradient. The left-most region is anoxic only, the middle is oxic or anoxic depending on historical conditions, and the right-most is oxic only. The thick grey vertical lines at -8 and 0 on the x-axis show the extent of the oxygen diffusivity gradient investigated in the simulations. In panels e, and i, the orange region shows when patterns of coexistence were atypical of the classical pattern shown in Figure 1a of the main text (e.g. when cyanobacteria coexisted with sulfur bacteria at low oxygen diffusivities). Right column of panels: The effect size of trait variation on resilience, with positive values on the y-axis indicating that diversity has increased resilience.

4.2.2 An example of stable states

Figure 4.5 shows an example of stable states along the oxygen diffusivity gradient produced by the replication method for finding stable states. The patterns are similar to those produced by the temporal method for finding the stable states (see Figure 4 in the main report), though they differ in one respect: close to the transition back to the state a group dominates, some of the more tolerant strains occur in the stable states found by the temporal method (Figure 4 in the main report), while they do not occur in stable states found by the replication method (Figure 4.5 here). We understand that the appearance of these strains is transient, and the difference in their appearance results from the differences in starting conditions of the two methods. Specifically, in the temporal method, the initial condition at a new level of oxygen diffusivity is the stable state at the previous level of oxygen diffusivity, while in the replication method the starting conditions are identical at each level of oxygen diffusivity. These differences in starting conditions apparently result in longer duration of transients in the temporal method than the replication method.

An interesting feature of these results is that although the position of the oxic to anoxic tipping point is affected by diversity in the sulfate-reducing bacteria (Figure 4.4 and Figure 3 in the main report) there is no evidence of more tolerant strains of sulfate-reducing bacteria being present in the final system state (Figure 4.5). This phenomenon is explained in Section 9.2.

Figure 4.5: Stable states found with the replication method with intermediate levels of diversity in all functional groups. Figure features are otherwise the same as those in Figure 4.1).

4.3 Effect of abundance change event

As mentioned in the main report, to avoid very low strain densities, we added 1 unit of density to every biological state variable every 1,000 hours. We here show that a slightly different method produced the same results. Specifically, we tested the effect of making an abundance floor of 1 unit that was implemented every 1,000 hours. I.e. every 1,000 hours, any abundances that were lower than 1 unit were set to 1 unit. As shown in Figure 4.6, the results with this floor abundance are similar as when 1 was added every 1,000 hours.

Figure 4.6: Results with an abundance floor. Left column of panels: the effect of trait diversity (y-axis) on the positions of the tipping points (anoxic to oxic in blue, oxic to anoxic in red), which determine the extent of the region of bistability (green shaded region) in the oxygen diffusivity gradient. The left-most region is anoxic only, the middle is oxic or anoxic depending on historical conditions, and the right-most is oxic only. The thick grey vertical lines at -8 and 0 on the x-axis show the extent of the oxygen diffusivity gradient investigated in the simulations. In panels e, and i, the orange region shows when patterns of coexistence were atypical of the classical pattern shown in Figure 1a of the main text (e.g. when cyanobacteria coexisted with sulfur bacteria at low oxygen diffusivities). Right column of panels: The effect size of trait variation on resilience, with positive values on the y-axis indicating that diversity has increased resilience, i.e. delayed the shift to the other state.

4.4 Further thoughts

Here we give some further thoughts regarding the comparison of the temporal and replication method for stable state finding.

A particularly notable difference between the two approaches is in the initial condition. In the temporal approach, the initial condition at a new level of oxygen diffusivity is the stable state at the previous level of oxygen diffusivity, while in the approach of Bush et al 2017 the initial conditions are fixed and arbitrary (though this approach is still useful for illustrating the bistability of the system).

The temporal approach is particularly useful when investigating the effects of strain richness. In preliminary simulations using the Bush et al 2017 approach, we found the counterintuitive result that higher strain richness led to lower resilience. We then recalled that we chose to divide the initial total abundance of the cyanobacterial (CB) functional group equally among the strains (termed a substitutive design in biodiversity manipulation experiments). This meant that the most tolerant CB strain in a three-strain simulation had higher initial abundance than the same most tolerant CB strain in a nine-strain simulation. The greater abundance of the most CB tolerant strain made it possible for a three-strain system to develop to be oxic, whereas the nine-strain system developed to be anoxic, despite identical range of trait variation. We could have chosen to do otherwise, for example start with the least tolerant strains of cyanobacteria as more abundant when the system started with high cyanobacterial abundance, and the most tolerant strain as more abundant when the system starts with low cyanobacterial abundance. Using the temporal approach removed the need to specify initial abundances (except the very first one), and therefore results were not so determined by the assumptions used to determine the initial conditions.

5 Finding the separatrix

In the original publication (Bush et al 2017) Figure 3a shows a separatrix in the cyanobacteria abundance - oxygen diffusivity dimensions of the system. If the system has initial cyanobacterial abundance greater than the value at the separatrix then the cyanobacterial abundance increases, and the system goes to the oxic stable state. Conversely, if the system has initial cyanobacterial abundance lower than the value at the separatrix then the cyanobacterial abundance decreases, and the system goes to the anoxic stable state.

It is critical to note, however, that the position of the separatrix depends on both the parameters of the system, and the initial conditions of variables other than cyanobacterial abundance. If the system were started with higher initial concentration, then the position of the separatrix would change (as would the values of oxygen diffusivity at which the system tipped from anoxic to oxic, and from oxic to anoxic).

This issue of dependence of the position of the separatrix on the value of initial condition of multiple state variables has particular significance when there are multiple strains per functional group. Then one cannot just vary initial total cyanobacterial abundance (as was done to create the separatrix on Figure 3a in the original publication) but one must also define how the initial abundance of each of the strains varies. One might choose to make each strain have equal initial abundance, and to vary the total (as we did in our implementation of the replication method in the main article).

We could, however, and arguably should, give the most tolerant strain of cyanobacteria higher initial relative abundance when total initial cyanobacterial abundance was low, and the least tolerant strain higher initial relative abundance when total cyanobacterial abundance was high. And we would need to specify the total and relative abundance of all cyanobacteria (and sulfur bactria) strains. Thus, with a system as complex as the 9 strain system, which has 31 state variables the separatrix is very high dimensional, and evaluating (and plotting) it in one dimension requires many assumptions. We chose to not attempt to find a meaningful set of assumptions in this study, though not because we think it would be unprofitable to do so, but rather because we did not want to risk distracting attention from findings we consider to be significant enough in their own right.

6 Effect of simulation duration

6.1 Diversity effects

Recall that the temporal method for finding stable states involves running the simulation for a certain amount of time, termed the “wait time”, at a constant level of oxygen diffusivity. At the end of that wait time the state of the system is recorded, and then the oxygen diffusivity is slightly increased or decreased, and the system then simulated again for the duration of the wait time.

The following figures 6.1-6.4 show the effects of diversity with long to shorter wait times.

Figure 6.1: Effects of diversity with wait time of 1e6 hours (same figure as in main report). (All else as in (Figure 4.4).

Figure 6.2: Effects of diversity with wait time of 1e5 hours. (All else as in (Figure 4.4).

Figure 6.3: Effects of diversity with wait time of 1e4 hours. (All else as in (Figure 4.4).

Figure 6.4: Effects of diversity with wait time of 1e3 hours. (All else as in Figure 4.4).

6.2 1e6 vs 2e6 hours wait time

The results below show the stable states as a function of oxygen diffusivity for a system with an intermediate level of diversity in all functional groups, with either 1’000’000 hours at each oxygen diffusivity level (Figure 6.5), or 2’000’000 hours (Figure 6.6). This shows that 1’000’000 hours duration of the simulations are sufficient to reach stable state, except in regions of strain replacement, including at some state transitions. E.g. in the region of strain replacement in the sulfate-reducing bacteria there is difference in strain abundances at 1’000’000 and 2’000’000 hours. However, the position of the tipping points (which were calculated from the total SB abundance) and the substrate concentrations did not differ between 1’000’000 and 2’000’000 hours.

Figure 6.5: Stable states found with the temporal method with duration at each oxygen diffusivity value of 1e6 hours. The left column of panels refers to the stable state when oxygen diffusivity is at first high (favouring the oxic state). The right column refers to when oxygen diffusivity is at first low (favouring the anoxic state). Points (which sometimes are so close as to appear as thick lines) show the stable state at the end of each period of constant oxygen diffusivity. Thin lines are for visualisation only, and join the points. Grey arrows show the oxygen diffusivity at which the system flips, and the direction of the flip. The length of the grey arrows is without meaning.

Figure 6.6: Stable states found with the temporal method with duration at each oxygen diffusivity value of 2e6 hours. All other aspects are as in Figure 6.5.

7 Effect of number of strains

In the main text, we present results of simulations with nine strains per functional group. Here we also explore the effect of number of strains (2, 3, 6, 9) on the position of the two tipping points (Figure 7.1). At low to moderate amounts of trait variation, the number of strains did not affect the position of the tipping points. However, at high amounts of trait variation, resilience was lower when the functional groups had only two or three strains rather than six or nine. This effect was particularly pronounced when there was high variation in only SB or in SB-CB: here, simulations with two or three strains led to dramatically lower resilience of the anoxic state than simulations with six or nine strains. We believe that this effect is due to the greater distance in trait space between strains in simulations with two or three strains than in simulations with more strains (also see Section 9).

Figure 7.1: Effects of number of strains and of diversity among strains. In panels PB and CB-PB the region with dashed lines shows when patterns of coexistence were atypical of the classical pattern shown in Figure 1a of the main text (e.g. when cyanobacteria coexisted with sulfur bacteria at low oxygen diffusivities).

8 CB diversity effects

Increased CB diversity causes a slight decrease in the resilience of the anoxic state. The following two figures show the state responses to oxygen diffusivity with low (Figure 8.1) and high (Figure 8.2) levels of CB diversity. The difference in the position of the anoxic to oxic tipping point is small. It is, in fact the smallest it could be, as it is the smallest step size in oxygen diffusivity used in the simulations.

Figure 8.1: Stable state of the system with low level of diversity in only CB, and no diversity in the other two groups. All other aspects are as in Figure 6.5.

Figure 8.2: Stable state of the system with high level of diversity in only CB, and no diversity in the other two groups. All other aspects are as in Figure 6.5.

9 SB diversity effects

9.1 SB diversity effect on the anoxic to oxic transition for short wait times (1e4)

We showed above that wait times of 1e5 (Figure 6.2) and 1e4 (Figure 6.3) produce qualitatively and often quantitatively very similar results to those for wait times of 1e6 (Figure 6.1). There is, however, a notable difference when there is variation in only SB, or in CB-SB with wait times of 1e4 (Figure 6.3). Here, moderate increases in diversity have the same effect as with wait times of 1e6. However, larger increases in diversity suddenly cause a large reduction in resilience of the anoxic-oxic transition. Figure 9.1 gives an example of the stable states when diversity is moderate – increases in oxygen diffusivity cause strain replacement in the SB to the more tolerant strains, which then delay the shift to the oxic state. Suprisingly, when there is just a little more diversity in the SB, then the anoxic to oxic transition happens earlier, i.e. there is a large decrease in resilience (Figure 9.2).

Figure 9.1: Stable state of the system with moderate level of diversity in only SB, and no diversity in the other two groups. All other aspects are as in Figure 6.5.

Figure 9.2: Stable state of the system with higher than moderate level of diversity in only SB, and no diversity in the other two groups. All other aspects are as in Figure 9.1.

We believe that this phenomenon is caused by an effect that we explain with reference to an analogy of using stepping stones to cross a river. Increasing diversity in our simulations is akin to increasing the width of the river, but keeping the number of stepping stones (strains) constant. So long as the river is not too wide, we can still cross it, but once the width of the river is such that the stepping stones are too far apart, we can’t cross the river. And this happens suddenly, the width increases just a tiny bit more, and we suddenly can’t cross the river when we could before.

In the simulated microbial ecosystem, the strains are the stepping stones, and it is harder (takes longer) for strains to replace each other when there is greater distance (in trait space) between them. When they can no longer replace each other in the time available, the system behaves as if only the least tolerant, highest growth rate strain (lightest shade of colour) is present. In the case illustrated above, this leads to much lower resilience of the anoxic system state.

This phenomenon is dependent on having static trait values for the strains. If strain traits were dynamic (e.g. due to evolution), then we expect this result to not occur, and rather for the effects of diversity with two or three strains to be the same as those with six or nine.

9.2 SB diversity effect on the oxic to anoxic transition

We show in the main text that diversity in the sulfate-reducing bacteria results in earlier recovery of the anoxic state (Figure 3 c). However, there is only limited evidence of strain replacement on the trajectory of decreasing oxygen diffusivity: the most tolerant strain does increase in density prior to the tipping point but never reaches dominance – it is the strain with lowest tolerance that dominates immediately after the shift to the anoxic state. We here analyze temporal dynamics rather than the final state and show that the more tolerant strains do play a role, but only a transient one. The following dynamics are for oxic starting conditions and with a value of oxygen diffusivity \(log_{10}\)(oxygen diffusivity) = -4.8. When there is no diversity within functional groups (i.e. two strains that are identical) the system remains oxic (Figure 9.3). When there is diversity in the two groups of sulfur bacteria the system switches to the anoxic state (Figure 9.4). (Note that the average trait value [on a log scale] of the two different strains is the same as the trait value of both strains in Figure 9.3.) The more tolerant strain of sulfate-reducing bacteria increases in density in the early phase of the simulation, reduces the concentration of oxygen by inhibiting cyanobacteria, and thus generates an environment in which the less tolerant strain is able to grow. Then, the more tolerant strain is outcompeted by the less tolerant strain as the concentration of the inhibiting substrate (oxygen) declines.

Figure 9.3: Dynamics of the system with no diversity in any functional group, and log10(oxygen diffusivity) = -4.8. The system remains oxic when starting in the oxic state. Although two strains are shown in the legend, they have identical parameters, and so their dynamics are identical and only one is visible in the dynamics.

Figure 9.4: Dynamics of the system with diversity in the two groups of sulfur bacteria and otherwise exactly the same conditions as in Figure 9.3. The system flips from oxic to anoxic, whereas it did not when there was no diversity.

10 PB diversity effects

10.1 Effect of low to moderate PB diversity

Increases in the diversity within the PB functional group leads to lower resilience of the anoxic state (Figure 6.1, Panel PB, blue line). For example, a small amount of variation in only the PB functional group, and no variation in other functional groups, results in transition from anoxic to oxic at approximately \(log_{10}\)(oxygen diffusivity) = -2.76 (Figure 10.1). Whereas a moderate amount of variation in only the PB functional group, and no variation in other functional groups, results in transition from anoxic to oxic at approximately \(log_{10}\)(oxygen diffusivity) = -3.34 (Figure 10.2).

Figure 10.1: Right-hand panels show the transition from anoxic to oxic state when oxygen diffusivity decreases. The results here are for a low amount of variation in only the PB functional group. The transition from anoxic to oxic occurs at approximately \(log_{10}\)(oxygen diffusivity) = -2.76. All other aspects are the same as in Figure 6.5.

Figure 10.2: Right-hand panels show the transition from anoxic to oxic state when oxygen diffusivity decreases. The results here are for a moderate amount of variation in only the PB functional group. The transition from anoxic to oxic occurs at approximately \(log_{10}\)(oxygen diffusivity) = -3.34. All other aspects are the same as in Figure 6.5.

This result is different to the effect of diversity within the CB and SB functional groups, which both cause increases in resilience of the state that they dominate (Figure 6.1). Why does diversity in the PB functional group reduce the resilience of the state it dominates?

The counter-intuitive effect of diversity in the phototrophic sulfur bacteria (PB) possibly results from their positive effects on conditions and processes that favour the oxic state. The phototrophic sulfur bacteria oxidise reduced sulfur to oxidised sulfur. Since reduced sulfur inhibits the cyanobacteria (CB), PB can have a positive effect on CB via their consumption of reduced sulfur. Hence, high growth rates of PB will have a positive effect on CB, and therefore favour the oxic state. We would expect to see that higher diversity of PB is associated with lower concentrations of reduced sulfur and higher concentrations of oxidised sulfur before the transition.

Another pathway by which a higher growth rate of PB can favour the oxic state is by indirectly increasing the concentration of oxygen, thus inhibiting themselves and sulfate-reducing bacteria. This effect occurs because PB consume reduced sulfur, such that less reduced sulfur is available for abiotic oxidation. Consequently, a larger amount of oxygen remains in the environment and inhibits sulfur bacteria. We would expect to see that higher diversity of PB is associated with higher concentrations of oxygen before the transition.

Among these two pathways, effects of phototrophic sulfur bacteria via abiotic oxidation of sulfide are probably of subordinate importance in natural environments because abiotic oxidation rates of sulfide are very small compared to biotic oxidation rates (Luther et al 2011).

There is, however, at least one pathway by which PB can favour the anoxic state: PB produce oxidised sulfur which is consumed by SB, and SB produce reduced sulfur, which inhibits CB.

Looking at the strain replacement that occurs across the oxygen diffusivity gradient, we see that as oxygen diffusivity increases, the more oxygen tolerant strains of PB replace the less tolerant ones because they have higher realised growth rates under the prevailing environmental conditions. Greater PB diversity leads to the availability of strains with higher oxygen tolerance and thus even higher realised growth rates, which via the pathways described above, can favour the oxic state.

Substrate concentrations close to the anoxic-oxic transition in simulations with low and high PB diversity, respectively, corroborate the expected effects (Figure 10.3). As the transition is approached, higher PB diversity leads to higher concentration of oxygen, higher concentration of oxidised sulfur (SO), and lower concentration of reduced sulfur (SR), than a simulation with lower PB diversity. All these differences in substrate concentrations are in the direction of favouring the oxic state, resulting in the earlier transition to the oxic state.

Figure 10.3: Concentrations of three substrates (O: oxygen; SO: oxidised sulfur; SR: reduced sulfur) in a low PB diversity simulation (x-axis) versus substrate concentration in a moderate PB diversity simulation (y-axis). (The same levels of PB diversity as in Figures 10.1 and 10.2, respectively). The figures show only a small range of oxygen diffusivity on the trajectory of increasing oxygen diffusivity. The black line is the 1:1 line.

10.2 Effect of medium to high PB diversity

When there was medium to high variation in only the phototrophic sulfur bacteria (PB), the strain dynamics deviated from the standard pattern: at low oxygen diffusivity, all three functional groups coexisted. Specifically, on the trajectory of decreasing oxygen diffusivity, cyanobacteria reappeared at low oxygen diffusivity after their collapse (medium variation in PB, Figure 10.4) or never fully collapsed (high variation in PB, Figure 10.5). On the trajectory of increasing oxygen diffusivity, cyanobacteria coexisted with the two groups of sulfur bacteria at low oxygen diffusivity, collapsed at medium levels of oxygen diffusivity, and then increased to sole dominance at the tipping point to the oxic state (Figure 10.4, Figure 10.5).

The same pattern occurred for medium to high simultaneous variation in phototrophic sulfur bacteria and cyanobacteria, albeit only on the trajectory of increasing oxygen diffusivity (Figure 10.6). On the trajectory of decreasing oxygen diffusivity, cyanobacteria dominated over the entire gradient of oxygen diffusivity, such that effects of PB diversity did not play out on this trajectory.

The reason for coexistence of all three groups at low oxygen diffusivity is likely that the PB strain that dominated at low oxygen diffusivity (the least tolerant strain) lowered the sulfide concentrations to such low levels that cyanobacteria were able to grow. The higher the trait variation in PB, the higher the maximum growth rate of the least tolerant PB strain, and the higher its consumption of sulfide. Consequently, moderate variation in PB lowered the sulfide concentration such that the tipping point from anoxic to oxic occurred earlier (as discussed in the previous section), and medium to high variation in PB lowered the sulfide concentration even further such that cyanobacteria were able to grow and coexist with the two groups of sulfur bacteria at low oxygen diffusivity.

At intermediate levels of oxygen diffusivities, the cyanobacteria were excluded (Figure 10.4). We speculate that this might be associated with the increase in sulfide prior to the collapse of the phototrophic sulfur bacteria. This increase in sulfide occurred simultaneously with turnover to PB strains with lower maximum growth rate and with a decline in PB density; these two processes likely resulted in less consumption of sulfide.

wait_time <- "1e+06"
ss_9s <- readRDS(here::here("data",
                            microxanox_version_path, 
                            "temporal_method",
                            "event_definition_2",
                            paste0("ss_data_9strains_waittime", wait_time, "_event_definition_2.RDS")))



max_CB_var_gmax <- max(ss_9s$CB_var_gmax_s)
max_SB_var_gmax <- max(ss_9s$SB_var_gmax_s)
max_PB_var_gmax <- max(ss_9s$PB_var_gmax_s)

ss_9s <- ss_9s %>%
  mutate(standvar_CB_gmax = CB_var_gmax_s / max_CB_var_gmax,
         standvar_SB_gmax = SB_var_gmax_s / max_SB_var_gmax,
         standvar_PB_gmax = PB_var_gmax_s / max_PB_var_gmax)

Figure 10.4: Coexistence of all three groups at low oxygen diffusivities when only PB contains strain variation, and the amount of strain variation is medium to high. This example is when the standardised amount of trait variation is 0.74 (0 being no variation, and 1 being the maximum amount in any simulation.)

Figure 10.5: High trait variation in only the PB group results in coexistence of all groups at low oxygen diffusivity levels (as in Figure 10.4) and also no collapse in CB abundance in the decreasing oxygen diffusivity trajectory. This example is when the standardised amount of trait variation is 0.95 (0 being no variation, and 1 being the maximum amount in any simulation.)

Figure 10.6: Stable states when there is high trait variation in both CB and PB groups. This example is when the standardised amount of trait variation is 0.95 (0 being no variation, and 1 being the maximum amount in any simulation.)

11 Underlying mechanisms

In the main text we show that trait variation (functional diversity) can affect resilience. In this supplementary section, we discuss the mechanism(s) underlying the observed diversity effects. We first explain the two major classes of mechanisms known to drive diversity effects, and then use additional simulations to elucidate the mechanisms operating in our study.

11.1 Complementarity and sampling effect

Research about the relationship between biodiversity and ecosystem functioning has identified two major mechanisms by which diversity affects ecosystem functioning: the selection/sampling effect and the complementarity effect (Loreau & Hector 2001). The sampling effect occurs when more diverse communities have higher functioning because they have a higher probability of containing a high performing species, and when only one (or few) of the species in a community determine the community and ecosystem level properties, such as biomass production or stability. In this case the observed effect of diversity is driven by only one (or few) species (Cardinale et al. 2011). In contrast, the complementarity effect occurs when higher ecosystem functioning in more diverse communities is due to niche differences or positive interactions among species, i.e. the diversity effect is driven by combinations of the species in a community (Cardinale et al. 2011). Both effects are considered valid and important mechanisms by which diversity affects ecosystem functioning (Cardinale et al. 2006, 2011, 2012; Fargione et al. 2007).

Similarly, the complementarity effect and the sampling effect are important mechanisms by which diversity affects stability (Loreau & Mazancourt 2013; Loreau et al. 2021) and resilience (Steiner et al. 2006; Hughes & Stachowicz 2011). For example, diversity can stabilize aggregate ecosystem properties (e.g. community biomass) when species differ in their responses to environmental factors and therefore respond asynchronously to environmental fluctuations (temporal complementarity effect; Loreau & Mazancourt 2013). However, the sampling effect might be of particular importance when considering ecosystem responses to directional environmental change or pulse disturbances rather than to environmental fluctuations (Steiner et al. 2006): more diverse systems might absorb a greater amount of change by having a higher chance of containing species with high environmental tolerance. In this scenario, the positive diversity effect on resilience is driven by one (or few) species. From a conservation perspective, however, one would nevertheless strive to maintain diversity and not only the most tolerant species because of the uncertainties associated with environmental change: we do not know with sufficient certainty which stressors will change in a given ecosystem and which traits will therefore be relevant in the future.

11.2 Investigating the mechanisms operating

We investigated if the diversity effects observed in our study were in line with the sampling effect, that is, driven entirely by the presence of the most tolerant strain or if other strains played a role as well. To this end, we ran additional simulations where functional groups contained only the strain with greatest tolerance (which due to the tradeoff also had the lowest maximum growth rate). That is, rather than increasing diversity we only increased tolerance of a single strain. Specifically, for each of the 20 levels of trait variation, we made simulations where each functional group had only one strain, and its trait values corresponded with those of the most tolerant strain for the given level of trait variation. We also manipulated in which of the functional groups tolerance increased. For example, in a simulation to be compared with variation in only the SB group (and no variation in the CB or PB groups), the SB group contained the strain with highest tolerance, whereas CB and PB each contained the strain with mean tolerance. We then compared the tipping point positions resulting from this tolerance manipulation with the tipping point positions resulting from the diversity manipulation presented in the main text. In addition, we compared the strain dynamics across the oxygen diffusivity gradient of the tolerance manipulation and the diversity manipulation. If the observed diversity effects are caused by the sampling effect, we expect there to be no difference between the tolerance manipulation and the diversity manipulation.

We found no difference between the diversity and the tolerance manipulation when trait values were low (Figure 11.1). When trait values were medium to high, however, we did find differences between the diversity and tolerance manipulation for most combinations except for diversity/tolerance in only CB and in SB-PB (Figure 11.1). Specifically, the diversity and tolerance manipulation differed either in tipping point positions (CB-SB and CB-SB-PB) or in whether coexistence patterns were atypical (SB, PB, CB-PB, CB-SB). We describe and discuss these differences in detail in the sections below. Briefly, in these scenarios less tolerant strains did play a role because of their high maximum growth rates. For example, absence of SB strains with high growth rate led to reduced production of sulfide and therefore to coexistence with cyanobacteria at low oxygen diffusivity or to increased resilience of the oxic state. In contrast, absence of PB strains with high growth rate led to reduced consumption of sulfide and therefore precluded coexistence with cyanobacteria at low oxygen diffusivity (as observed when there was high variation in PB).

Taken together, at low amounts of trait variation the diversity effects were driven entirely by the most tolerant strain, in line with a sampling effect. At medium to high trait variation, however, in most combinations both tolerant and less tolerant strains contributed to the diversity effects on tipping point positions or coexistence patterns. The reason was that in the two groups of sulfur bacteria the trade-off between tolerance and growth rate had ecosystem-level effects, i.e. effects on the sulfide concentration and therefore on the cyanobacteria.

Figure 11.1: Effects of strain variation on position of tipping points compared between simulations with all (nine) strains present and simulations where only the most tolerant strain is present. When tipping points are shown with dashed lines then patterns of coexistence were atypical of the classical pattern shown in Figure 1a of the main text (e.g. when cyanobacteria coexisted with sulfur bacteria at low oxygen diffusivities).

11.3 Increasing tolerance vs. increasing diversity in only the SB

When SB contained multiple strains varying in growth rate (and tolerance) there were two stable states at intermediate levels of oxygen diffusivity, only one stable state at low oxygen diffusivity with only sulfur bacteria, and only one stable state at high oxygen diffusivity with high cyanobacterial density (Figure 11.2), as was previously reported. When, however, only one SB strain with very high tolerance and low maximum growth rate was present, the stable states change: at low oxygen diffusivity all three functional groups coexisted (Figure 11.3) – this never occurred when there were multiple SB strains varying in tolerance.

Apparently, at low oxygen diffusivity levels, the absence of SB strains with high growth rate precludes the formation of a stable state that does not include cyanobacteria. The most tolerant SB strain has a low maximum growth rate, which leads to less production of sulfide. Presence of only the most tolerant SB strain therefore results in lower sulfide concentrations (Figure 11.2 vs Figure 11.3), allowing growth of cyanobacteria at low oxygen diffusivity.

However, cyanobacteria are excluded at intermediate oxygen diffusivity. We speculate that this might be associated with the increase in sulfide at intermediate oxygen diffusivity, which might be due to the decline and collapse of PB and therefore reduced consumption of sulfide.

Figure 11.2: Stable states found when there is variation among nine strains of SB and the amount of variation is high. All other aspects are the same as in Figure 6.5.

Figure 11.3: Stable states found when SB contained only the most tolerant strain, with trait values corresponding with the diversity level of the previous figure. All other aspects are the same as in Figure 6.5.

In the ecosystem shown in Figure 11.3, where the SB contained one strain with very high tolerance, there was large variation in state variables along the oxygen diffusivity gradient at values less than about -4.5. To examine what might be causing this, we simulated the ecosystem at a constant oxygen diffusivity of \(log_{10}\) oxygen diffusivity of -6 (Figure 11.4). We see that the dynamics are qualitatively different from other cases in which there is a stable state reached in such a simulation (Figure 9.4). It appears that the high tolerance SB strain is able to initiate a transition to an anoxic state, but is not able to create a stable anoxic state. This dynamical situation was not present in the ecosystem in which high tolerance and low tolerance strains were present.

Figure 11.4: Temporal dynamics of an ecosystem with a single SB strain with high tolerance, the same one as in Figure 11.3 simulated with constant \(log_{10}\) oxygen diffusivity of -6.

11.4 Increasing tolerance vs. increasing diversity in CB-SB

In the CB-SB case, the presence of only the most tolerant strain led to coexistence of all groups at low levels of oxygen diffusivity when trait values were moderately high (Figure 11.5) or to absence of alternate stable states and sole dominance of cyanobacteria when trait values were high (Figure 11.6). As discussed in the previous section, the low maximum growth rate of the most tolerant SB strain led to less production of sulfide such that SB could not exclude CB.

Figure 11.5: Stable states found when CB and SB each contain only the single most tolerant strain with moderate tolerance. All other aspects are the same as in Figure 6.5.

Figure 11.6: Stable states found when CB and SB each contain only the single most tolerant strain with high tolerance. All other aspects are the same as in Figure 6.5.

11.5 Increasing tolerance vs. increasing diversity in PB and in PB-CB

Whereas medium to high diversity in only PB led to coexistence of all three functional groups at low oxygen diffusivity levels (Figure 11.7), presence of only the most tolerant PB strain led to the standard-pattern of either dominance of sulfur bacteria or of cyanobacteria (Figure 11.8). The most tolerant PB strain has a low maximum growth rate, leading to less consumption of sulfide and therefore to higher sulfide concentrations compared to a system that contains both slow- and fast-growing strains. Due to the higher sulfide concentrations in the system with only the most tolerant PB strain, cyanobacteria were not able to grow at low oxygen diffusivity levels.

For the same reasons we found a similar pattern on the trajectory of increasing oxygen diffusivity for medium to high simultaneous trait variation in PB and CB but not for the corresponding tolerance manipulation (data not shown).

wait_time <- "1e+06"
ss_9s <- readRDS(here::here("data",
                            microxanox_version_path, 
                            "temporal_method",
                            "event_definition_2",
                            paste0("ss_data_9strains_waittime", wait_time, "_event_definition_2.RDS")))



max_CB_var_gmax <- max(ss_9s$CB_var_gmax_s)
max_SB_var_gmax <- max(ss_9s$SB_var_gmax_s)
max_PB_var_gmax <- max(ss_9s$PB_var_gmax_s)

ss_9s <- ss_9s %>%
  mutate(standvar_CB_gmax = CB_var_gmax_s / max_CB_var_gmax,
         standvar_SB_gmax = SB_var_gmax_s / max_SB_var_gmax,
         standvar_PB_gmax = PB_var_gmax_s / max_PB_var_gmax)

Figure 11.7: Coexistence of all three groups at low oxygen diffusivities when only PB contains strain variation, and the amount of strain variation is medium to high. This example is when the standardised amount of trait variation is 0.89 (0 being no variation, and 1 being the maximum amount in any simulation.)

wait_time <- "1e+06"
ss_tol <- readRDS(here::here("data",
                            microxanox_version_path, 
                            "temporal_method",
                            "event_definition_2",
                            paste0("ss_data_-1strains_waittime", wait_time,
                                   "_event_definition_2.RDS")))



max_CB_var_gmax <- max(ss_tol$CB_var_gmax_s)
max_SB_var_gmax <- max(ss_tol$SB_var_gmax_s)
max_PB_var_gmax <- max(ss_tol$PB_var_gmax_s)

ss_tol <- ss_tol %>%
  mutate(standvar_CB_gmax = CB_var_gmax_s / max_CB_var_gmax,
         standvar_SB_gmax = SB_var_gmax_s / max_SB_var_gmax,
         standvar_PB_gmax = PB_var_gmax_s / max_PB_var_gmax)

Figure 11.8: Dynamics when PB contain only the most tolerant strain, with this being the same tolerance level as the most tolerant strain in Figure 11.7 immediately above.

11.6 Increasing tolerance vs. increasing diversity in CB-SB-PB

At high levels of diversity / high levels of tolerance, resilience of the oxic state is greater in the tolerance than diversity manipulation. We speculate that this is because the most tolerant SB strain has reduced ability to suppress the cyanobacteria, as we have also observed in the scenarios of tolerance in only SB and in SB-CB (described above).

Figure 11.9: Stable states found with nine strains of CB-SB-PB and high diversity. All other aspects are the same as in Figure 6.5.

Figure 11.10: Stable states found with the single most tolerant CB strain, SB strain, and PB strain, with trait values corresponding with the diversity level in the previous figure. All other aspects are the same as in Figure 6.5.

12 (Non-)additivity

We investigated if the effects of trait variation of individual functional groups were additive. We compared the tipping point positions as determined by our simulations with the expected position of the respective tipping point if effects of diversity in different functional groups were additive. The expected tipping point was calculated by adding the effect sizes (i.e. the difference in tipping point position between the no-diversity scenario and a given amount of trait variation) on a linear scale.

Figure 12.1: Assessing the additivity or non-additivity of the effects of trait variation. Solid black lines depict the position of the two tipping points at different levels of trait variation, as simulated with the model. The colored, dashed lines show the expected position of tipping points if effects of diversity in different functional groups were additive. The expected tipping point was calculated by adding the effect sizes of trait variation of different functional groups on a linear scale. For example, in (d) the effect size of variation in only CB was added to the effect size of variation in only SB. In (g), the effect size of variation in only CB was added to the effect size of simultaneous variation in both SB and PB. (CB: cyanobacteria, SB: sulfate-reducing bacteria, PB: phototrophic sulfur bacteria). The dotted vertical lines at -8 and 0 on the x-axis show the extent of the oxygen diffusivity gradient investigated in the simulations.

13 Log and non-log oxygen diffusivity

For comparability with the Bush et al publication (2017), we use a log-scale of the oxygen diffusivity gradient when displaying the effect of trait variation on tipping point positions. In our work, we only make comparisons where the outcome is the same when a log- or an arithmetic-scale of oxygen diffusivity was used. For example, we compare the direction of the diversity effects of the three functional groups on the resilience of the ecosystem state that they dominate (CB: oxic state, i.e. oxic-anoxic tipping point; SB and PB: anoxic state, i.e. anoxic-oxic tipping point). The direction of the effects does not depend on the log/arithmetic scale decision.

In contrast, we avoid comparisons where the outcome would differ when a log- or an arithmetic-scale of oxygen diffusivity was used. An example of a comparison that we avoid is of the magnitude of the diversity effect between the oxic-anoxic transition and the anoxic-oxic transition. We avoid this comparison of magnitudes because with a log scale the diversity effect is larger for the oxic-anoxic transition (Figure 13.1n), while on an arithmetic scale the diversity effect is larger for the anoxic-oxic transition (Figure 13.1o).

Figure 13.1: The effect of trait diversity (y-axis) viewed with a log axis (left panels) or arithmetic axis (right panels) for the oxygen diffusivity gradient (x-axis). The positions of the tipping points (anoxic to oxic in blue, oxic to anoxic in red) and the extent of the region of bistability (green or orange shaded region) are shown. The left-most region is anoxic only, the middle is oxic or anoxic depending on historical conditions, and the right-most is oxic only. In the left panels the thick grey vertical lines at -8 and 0 on the x-axis show the extent of the oxygen diffusivity gradient investigated in the simulations. In panels e, f, i, and j the orange region shows when patterns of coexistence were atypical of the classical pattern shown in Figure 1a of the main text (e.g. when cyanobacteria coexisted with sulfur bacteria at low oxygen diffusivities). In the right-hand panels (arithmetic x-axis scale) the blue anoxic to oxic tipping point line often has steps. This is because the simulated oxygen diffusivity levels were evenly distributed on the log-oxygen diffusivity scale, and so are relatively far apart at high oxygen diffusivities on the arithmetic scale.

14 References

Bush, T., Diao, M., Allen, R.J., Sinnige, R., Muyzer, G. & Huisman, J. (2017). Oxic-anoxic regime shifts mediated by feedbacks between biogeochemical processes and microbial community dynamics. Nat. Commun., 8, 789.

Cardinale, B.J., Duffy, J.E., Gonzalez, A., Hooper, D.U., Perrings, C., Venail, P., et al. (2012). Biodiversity loss and its impact on humanity. Nature, 486, 59–67.

Cardinale, B.J., Matulich, K.L., Hooper, D.U., Byrnes, J.E., Duffy, E., Gamfeldt, L., et al. (2011). The functional role of producer diversity in ecosystems. Am. J. Bot., 98, 572–592.

Cardinale, B.J., Srivastava, D.S., Emmett Duffy, J., Wright, J.P., Downing, A.L., Sankaran, M., et al. (2006). Effects of biodiversity on the functioning of trophic groups and ecosystems. Nature, 443, 989–992.

Fargione, J., Tilman, D., Dybzinski, R., Lambers, J.H.R., Clark, C., Harpole, W.S., et al. (2007). From selection to complementarity: shifts in the causes of biodiversity–productivity relationships in a long-term biodiversity experiment. Proc. R. Soc. B Biol. Sci., 274, 871–876.

Hughes, A.R. & Stachowicz, J.J. (2011). Seagrass genotypic diversity increases disturbance response via complementarity and dominance. J. Ecol., 99, 445–453.

Loreau, M., Barbier, M., Filotas, E., Gravel, D., Isbell, F., Miller, S.J., et al. (2021). Biodiversity as insurance: from concept to measurement and application. Biol. Rev., 96, 2333–2354.

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

Loreau, M. & Mazancourt, C. de. (2013). Biodiversity and ecosystem stability: a synthesis of underlying mechanisms. Ecol. Lett., 16, 106–115.

Luther, G.W., Findlay, A., MacDonald, D., Owings, S., Hanson, T., Beinart, R., et al. (2011). Thermodynamics and kinetics of sulfide oxidation by oxygen: a look at inorganically controlled reactions and biologically mediated processes in the environment. Front. Microbiol., 2, 62.

Moisset de Espanés, P., Ramos-Jiliberto, R. & Soto, J.A. (2021). Avoiding artifacts when varying the number of species in ecological models. Ecol. Lett., 24, 1976–1987.

Steiner, C.F., Long, Z.T., Krumins, J.A. & Morin, P.J. (2006). Population and Community Resilience in Multitrophic Communities. Ecology, 87, 996–1007.