Environmental fluctuations reshape an unexpected diversity-disturbance relationship in a microbial community

Environmental disturbances have long been theorized to play a significant role in shaping the diversity and composition of ecosystems. However, an inability to specify the characteristics of a disturbance experimentally has produced an inconsistent picture of diversity-disturbance relationships (DDRs). Here, using a high-throughput programmable culture system, we subjected a soil-derived bacterial community to dilution disturbance profiles with different intensities (mean dilution rates), applied either constantly or with fluctuations of different frequencies. We observed an unexpected U-shaped relationship between community diversity and disturbance intensity in the absence of fluctuations. Adding fluctuations increased community diversity and erased the U-shape. All our results are well-captured by a Monod consumer resource model, which also explains how U-shaped DDRs emerge via a novel ‘niche flip’ mechanism. Broadly, our combined experimental and modeling framework demonstrates how distinct features of an environmental disturbance can interact in complex ways to govern ecosystem assembly and offers strategies for reshaping the composition of microbiomes.


Introduction
Biodiversity is a cornerstone of ecosystem stability and function (Tilman et al., 2014). While it is well appreciated that environmental changes influence species diversity in all ecosystems, the exact nature of this critical relationship is unclear. Without a predictive understanding of how ecosystems respond to perturbations, we are poorly prepared for environmental changes of anthropogenic origin, such as rising global temperatures (Hoegh-Guldberg and Bruno, 2010), and unable to design effective and robust interventions in ecosystems, such as microbiomes of medical or agricultural importance (Lemon et al., 2012;Widder et al., 2016). Accordingly, there have been many efforts aimed at understanding the role of environmental disturbances, which are perturbations to the state of an environment. These disturbances are of ecological interest for the impact they have on a community, for example, by bringing about mortality of organisms and a reduction of biomass of a community. Various diversity-disturbance (DDR) relationships have been proposed that draw intuition from observations of natural ecosystems. DDRs describe how community diversity depends on disturbance intensity, which we define as the mean mortality rate imposed by the disturbance over time. A famous example is the Intermediate Disturbance Hypothesis (Connell, 1978;Huston, 1979), in which species diversity peaks at intermediate disturbance intensities (Figure 1a). However, DDRs derived from observational studies of disparate ecosystems and disturbance regimes often have inconsistent results (Mackey and Currie, 2001;Hughes et al., 2007). Earlier assertions that disturbance weakens or interrupts competition (Connell, 1978;Huston, 1979)  In the laboratory, microbial communities can be cultivated and subjected to varying disturbance intensity levels by tuning the dilution rate in chemostats. (c) A bacterial community exhibits a U-shaped diversity dependence on the disturbance intensity. Samples of a soil-derived bacterial community were grown for 6 days in eVOLVER mini-chemostats at four different dilution rates. Top: Optical density over time quantifies biomass for replicate cultures. Middle: Mean relative abundance of bacterial genera from replicate cultures, determined by 16S sequencing. Mean rank abundance is denoted by order, taxonomic similarity is denoted by color. Bottom: Plotting the endpoint number of species (Amplicon Sequence Variants) vs. dilution rate produces a U-shaped curve, rather than a peaked DDR. Shaded window indicates a one standard deviation confidence interval.
theory (Chesson and Huntly, 1997;Fox, 2013) and empirical findings (Violle et al., 2010) that harsher environments instead reinforce dependence on limiting factors. Without a quantitative framework that directly pairs theory and experiment, it has been difficult to determine the source of disagreement between the many conflicting predictions and observations surrounding DDRs.
Importantly, the impact of a disturbance on an ecosystem depends on the disturbance characteristics. For example, although some environmental disturbances do not vary over time, many environmental disturbances occur with lower frequencies and introduce fluctuations that drive an ecosystem in and out of different states. The environmental fluctuations associated with a disturbance may in fact stabilize communities by creating temporal niches, similar to seasonal effects (Levins, 1979;Chesson, 1994). Indeed, coexistence can be promoted in a fluctuation-dependent manner due to storage effects (e.g. dormancy in poor conditions) or if species exhibit relative non-linearities in their competitive responses (e.g. differently shaped growth curves) (Chesson, 1994;Letten et al., 2018). Yet, coexistence might also arise from the overall time-averaged disturbance intensity in a fluctuation-independent manner (Chesson and Huntly, 1997;Fox, 2013). To determine whether the effects of disturbance on diversity are truly fluctuation-dependent (Chesson, 2000), a disturbance should ideally be decomposed into distinct components of mean intensity (e.g. time-averaged disturbance magnitude) and frequency (e.g. temporal profile of fluctuations). Indeed, experimental findings (Benedetti-Cecchi et al., 2006;Hall et al., 2012) and theory (Miller et al., 2011) have suggested that diverse DDRs could arise when considering these factors independently. There is therefore a need for comprehensive, controlled studies which pair theory with experimental methods to produce datasets that can deconvolve the effects of intensity and fluctuation.
Laboratory experiments offer a greater degree of control and throughput compared to field studies, particularly for tractable ecosystem models like microbial communities (De Roy et al., 2014). Microbes are easily quantified with next-generation sequencing (Gohl et al., 2016;Bolyen et al., 2019;Callahan et al., 2016;Janssen et al., 2018), and have been widely used in the laboratory to model community assembly (Friedman et al., 2017;Venturelli et al., 2018;Goldford et al., 2018), cross-feeding relationships (Mee et al., 2014), and succession (Chuang et al., 2019). Laboratory models have also linked changes in diversity in response to fluctuating nutrient levels (Sommer, 1985;Grover, 1988) and disturbances such as sonication (Violle et al., 2010), ultraviolet radiation (Gibbons et al., 2016), osmotic pressures (Letten et al., 2018), or toxic compounds (Santillan et al., 2019). Dilution is perhaps the most common choice for a laboratory disturbance, as it causes species-independent mortality and replenishes the system with fresh nutrients, reminiscent of flow in soil, aquatic, or gut microbiomes. Unlike disturbances with indirect biological impacts (such as pH, temperature, or osmolarity disturbances), there is a direct link between the dilution disturbance event (removal of culture volume) and the biological outcome (mortality of community members). In simple batch culture experiments, where cultures remain undisturbed except for a periodic dilution step, coexistence has been observed at intermediate dilution levels (Gibbons et al., 2016;Abreu et al., 2019), although different DDRs arise under different dilution regimes (Hall et al., 2012), suggesting that the dilution parameter space is vastly under-sampled. For more precise tuning of dilution or other parameters, experimentalists have long relied on continuous culture methods (Sommer, 1985;Grover, 1988;Monod, 1950); unfortunately, these systems have traditionally been intractable to large-scale, multidimensional experiments. Recently, we developed eVOLVER, a flexible and automated continuous culture platform that enables independent control over conditions in a large number of mini-bioreactors (Wong et al., 2018;Heins et al., 2019), thus opening up the possibility to explore microbial community dynamics under controlled, multidimensional environmental disturbances. By programming different dilution profiles with eVOLVER, we set out to independently quantify the effects of disturbance intensity (i.e. average dilution rate) and frequency (i.e. temporal profile of dilution rate) on the composition and diversity of microbial communities.

Results
Measuring microbial community diversity in constantly applied disturbances reveals a U-shaped relationship First, we sought to measure microbial diversity at various mean disturbance intensities in the absence of fluctuations by culturing a community under different mean dilution rates at the same near-continuous frequency. We cultivated replicate samples of a soil-derived microbiome in separate eVOLVER bioreactor arrays in dilute Nutrient Broth for six days (comprising 20-90 generations), during which continuously diluted cultures approached a stable composition (Figure 3-figure supplements 1 and 2). In a chemostat, the flow of media into the vessel is matched by flow of spent media and cells out of the vessel, so disturbance intensity is directly related to dilution rate ( Figure 1b). We thus varied the disturbance intensity by varying the dilution rate across the arrays (see Materials and methods, Figure 3-figure supplements 1 and 2). We sampled cultures daily and used 16S sequencing to quantify composition and diversity over time. As expected, we observed decreasing biomass of the cultures at increasing dilution rates (Figure 1c and Figure 3-figure supplement 3). Surprisingly, after quantifying the composition of each culture, we observed a U-shaped diversitydisturbance relationship (Figure 1c), with the number of surviving species at intermediate dilution rate at roughly half of the number at either low or high dilution rate. We were particularly intrigued because a U-shaped relationship between diversity and disturbance intensity does not agree with historical wisdom (Huston, 1979), although U-shaped DDRs have been reported for other disturbance characteristics such as frequency (Miller et al., 2011). Although U-shaped DDRs are uncommon in empirical observations (Mackey and Currie, 2001;Hughes et al., 2007), the conditions under which we observed it were quite straightforward: constantly applied disturbance. Thus, to better understand our observation, we sought a modeling framework in which a U-shaped DDR could emerge from constantly applied disturbance, while still capturing other reported DDR shapes.
A consumer-resource model captures experimentally observed DDRs and uncovers a niche-flip mechanism for conditional coexistence To develop a theoretical framework for understanding and predicting DDR behavior, we simulated microbial co-cultures in bioreactor conditions. Before examining higher order systems, we started by examining the simplest case that could give rise to a U-shaped DDR: a two-species competition where coexistence breaks down at intermediate disturbance levels. This simple system illustrates a sufficient condition that leads to the U-shaped DDR, and extensive analysis on a broader set of conditions is provided in Appendix 1 and 2. To link changes in disturbance intensity to changes in competitive outcomes, we turned to consumer resource models (Tilman, 1982;MacArthur, 1970). In consumer resource models, species growth rates are a function of resource concentrations. Resource depletion rates are in turn a function of species abundances, so species interact by competing for the same limiting resources. The range of resource concentrations that can support growth of a species can be graphically analyzed on a Tilman diagram (Tilman, 1982) by defining a Zero Net Growth Isocline (ZNGI) (Figure 2a). Briefly, the ZNGIs intersect with each axis in the Tilman diagram at c* values, which are the concentration of a particular resource necessary for a species to survive at the specified mortality rate. Similarly, the ZNGIs define the boundary of supplied resource combinations that will support growth of a species. As a population consumes resources, the resource concentrations move from the resource supply point toward the ZNGI, where growth rate is equal to the mortality rate (i.e. disturbance intensity) and the population is at equilibrium. The shape and location of these ZNGIs also predict the outcome of competition, as resource consumption by the population can cause equilibrium resource levels to cross the ZNGI of one species, leading to exclusion by the other (Figure 2a). If the ZNGIs intersect with each other, coexistence is possible. The borders of the coexistence region are invasion boundaries where either species can invade the other; in this region, consumption brings equilibrium resource concentrations towards the intersection of the ZNGIs (Figure 2a).
Disturbance affects the ZNGIs in multiple ways. First, since higher mortality rates require higher growth rates to maintain equilibrium, more resources are required to maintain equilibrium. Accordingly, the ZNGIs shift toward higher resource levels at higher disturbance intensities. Second, for saturating growth models, the growth rate of each species may change relative to each other as  Figure 2. A U-shaped diversity disturbance relationship can emerge from a consumer resource model if species undergo niche-flip as disturbance increases. (a) The survival of a species in a consumer resource model depends on the supplied resource levels and the mortality rate (i.e. disturbance intensity). A Zero Net Growth Isocline (ZNGI) may be defined for each species at a given disturbance intensity, delineating the range of resource supply levels where growth can meet or exceed the mortality rate. The ZNGIs intersect the axis at c* values, the minimum concentration needed to support Figure 2 continued on next page resource levels increase. For example, for two species competing for resource 1, species A may outgrow species B at low resource concentrations but not at high concentrations ( Figure 2a). Competitive outcomes at low disturbance intensities will depend on the relative growth rates at low resource concentrations, while competitive outcomes at high disturbance intensities will depend on the relative growth rates at high concentrations ( Figure 2b). The c* values and ZNGIs will change accordingly for each species. We investigated whether these disturbance-linked differences in competitive outcomes could explain the emergence of U-shaped DDRs. In a competition between two species, a U-shaped DDR can be generated if this coexistence region disappears at intermediate disturbance intensities (Figure 2b). We propose that this is possible if the ZNGIs and invasion boundaries flip as disturbance intensity increases, such that at some intermediate intensity the invasion boundaries align and the coexistence region disappears (Figure 2b). We term this behavior 'niche-flip'. Under niche-flip, the winner of competition at a given supply point changes as disturbance intensity varies concomitantly with resource availability.
To look for conditions under which niche flip emerges, we conducted simulations using Monod growth kinetics (Monod, 1949), a commonly used model for microbial growth with Type II functional response. Here, growth scales with the concentration of resource, but saturates to a maximal growth rate r according to a half-saturation constant K (Figure 2c and Appendix 1). Accordingly, a species with high maximum growth rate r may be outcompeted at low resource levels by a species with a low saturation constant K, such that the outcome of competition varies depending on nutrient levels (and thus dilution rate d) as described above (Figure 2-figure supplement 1). Theoretical analysis of the Monod model with two species and two resources shows that ZNGIs and invasion boundaries undergo niche-flip as dilution rate increases ( Figure 2b). We next asked whether these interactions between pairs of species produced the desired DDRs in a more complex system. We simulated sets of 10 species and seven resources with randomly chosen r and K values, with per-capita growth rates composed of a sum across nutrient-specific growth rates. Excitingly, we found that the Monod consumer resource model recapitulates the U-shaped diversity dependence on disturbance intensity that we observed in our chemostat experiments ( Figure 2c). Finally, we emphasize that the illustrated tradeoff is not the only way to obtain niche-flip and U-shaped DDRs (explored further in Appendix 2). To summarize briefly, in the absence of tradeoffs on any one resource, at sufficiently high dilution rates, niche flip and a U-shaped DDR can arise from combinations of two or more resources such that the winner of competition is different at low vs. high equilibrium resource concentrations. In systems of larger numbers of species and resources, niche-flip exclusion events between pairs of species can co-occur, yielding U-shaped DDRs (Appendix 2).

Other classes of microbial competition model do not produce U-shaped DDRs
We asked whether niche flip or U-shaped DDRs could emerge from other classes of microbial competition model. We simulated 10-species, 7-resource competitions with linear consumer resource models where the growth of a species increases linearly with resource concentration (Appendix 1). We also simulated 10-species competitions with the Lotka-Volterra model, where the maximal Resources are consumed until reaching equilibrium along the ZNGI of the species requiring the fewest resources to survive at the specified mortality rate, the winner of the competition. Different supply points can be found that support neither species (gray), a single species (green or purple) or both species (black). Invasion boundaries indicate regions where one species can increase in density in the presence of the other, outlining the region of coexistence (maroon). (b) In a consumer resource simulation with 2-species / 2-resources and Monod growth, the outcome of competition at the given supply point (black) depends on disturbance intensity. Both ZNGIs and invasion boundaries flip as disturbance intensity increases. At intermediate disturbance, invasion boundaries align and the coexistence region collapses, reducing diversity relative to low/high disturbance intensities at the black supply point. (c) Left: Monod consumer resource model for growth of species i with additive non-linear growth on each resource. Right: Shannon diversity of randomly generated 10-species communities, after six days of simulated growth on seven resources at varying dilution rates. For each model, mean diversity was computed for 100 randomly initialized communities, across each mean dilution rate, 50 of which are shown as individual traces. The online version of this article includes the following figure supplement(s) for figure 2: growth rates (r) of different species are reduced depending on the abundance of other species via interaction coefficients (a) (Appendix 1). Growth rates and interaction coefficients were selected to match the growth rates and diversity of the consumer resource models described above. U-shaped DDRs were not observed for either linear consumer resource or Lotka-Volterra models; instead, diversity generally did not vary with disturbance intensity, even under alternative parameter regimes (Appendix 1-figures 1 and 2). Thus, we concluded that the Monod consumer resource model best captured the U-shaped behavior of our microbial communities under constant disturbance among simple models of competition.
Experimentally measured growth parameters are consistent with Monod growth To confirm that Monod growth is a reasonable choice for modeling our experimental system, we measured how nutrient concentrations affected the growth of bacterial isolates from our community experiments. In 96-well plates, we first measured growth rates at varying concentrations of diluted Nutrient Broth, the same media type used in our eVOLVER experiments. Growth rates varied as a function of nutrient composition and indeed saturated at relevant nutrient concentrations, exhibiting tradeoffs between r and K such that no species had both a fast growth rate and low saturation coefficient ( Figure 3-figure supplement 4). Second, as Nutrient Broth is undefined, we repeated the above experiment in carbon-limited M9 minimal media supplemented with varying concentrations of amino acids as a sole carbon source. Although we were unable to fit Monod parameters due to limited growth, we observed that the isolate with highest growth rate varied for different amino acids (Figure 3-figure supplement 4). These two experimental observations provide evidence for Monod growth kinetics under these conditions, and suggest that species could undergo niche-flip to generate U-shaped DDRs under constant disturbance.

The consumer-resource model predicts fluctuation-dependent changes in DDRs
Both mean disturbance intensity and fluctuations (e.g. low-frequency discrete disturbance events) are hypothesized to play a role in the assembly of communities, but how these two disturbance components interact to reshape DDRs is unclear. Using our modeling framework, we sought to independently vary these two components, simulating a two-dimensional dilution profile. Specifically, we introduced fluctuations into the model by permitting d to vary with time, compressing disturbance into discrete time windows ( Figure 3a); this was done while keeping the time-averaged d equal, thereby allowing us to vary mean disturbance intensity and frequency independently. The Monod consumer resource simulations predict significantly higher diversity in fluctuating conditions comprised of one or more dilution events per day, with the lowest-frequency (i.e. largest-fluctuation) regime predicted to maintain the most diversity ( Figure 3e and Figure 3-figure supplement 5). This is consistent with intuition that fluctuations introduce temporal structure into environments, which may create new niches that promote diversity. Furthermore, the DDR is reshaped entirelyfrom U-shaped to largely uniform -indicating that community composition in the Monod model is conclusively fluctuation-dependent. Notably, neither Lotka-Volterra nor linear consumer resource models predict differences in the DDR between fluctuation frequencies (Appendix 1-figures 1 and 2). The overlap of DDRs of different frequencies indicates that in these models the relevant metric is the time-averaged overall intensity, rather than the frequency, of disturbance (Chesson and Huntly, 1997;Fox, 2013).

Experimentally introducing disturbance fluctuations increases community diversity and reshapes DDRs
We returned to experiments to see whether the U-shaped DDR observed in constant dilution chemostats is reshaped by fluctuations, as predicted by the Monod model. To comprehensively test for fluctuation-dependency, we implemented dilution profiles with 1, 4, or 16 fluctuations per day (alongside the constant dilution conditions) at varying mean dilution rates in an eVOLVER experiment comprised of 64 simultaneous cultures ('DDR64 Experiment') ( Figure 3b). We cultivated replicate samples of the soil-derived microbiome in eVOLVER for 6 days, taking samples every 24 hr to quantify composition (   (1) We observed U-shaped diversity curves in regimes of constant disturbance and small frequent disturbances, in both experiment and simulations.
(2) Larger fluctuations preserved higher levels of diversity, and (3) larger fluctuations reshaped the diversity-disturbance curves towards a more uniform relationship. These general results were consistent across different diversity metrics: Shannon diversity of 16S ASVs (Amplicon Sequence Variants), richness of 16S ASVs, and richness of colony morphotypes (Figure 3-figure supplement 11).
To test the reproducibility of our results, we ran a subsequent and independent 48-vial experiment from frozen inoculum at dilution rates ranging from 0.3 h À1 to 1.5 h À1 (Figure 3-figure supplements 1 and 3). In overlapping parts of the parameter space, we were able to reproduce the results of the DDR64 experiment ( Figure 3-figure supplements 7, 10 and 11). We further used this replicate experiment to extend the parameter space and examine effects at extreme levels of Source data 1. Population size and diversity metrics for DDR64 and DDR washout. Source data 2. Genus level composition for DDR64 and DDR washout.            disturbance. However, despite imposing dilution rates that exceed the maximum growth rate of our cultures (Figure 3-figure supplement 4), we did not observe complete washout (Figure 3-figure  supplement 3). In principle, cells could escape dilution by forming biofilm or due to incomplete mixing at extreme dilution rates, so we were cautious to draw strong conclusions from this part of the parameter space.
Although other measurables varied across the disturbance parameter space, they do not explain the differences in diversity as clearly as the Monod consumer resource model does (Figure 3-figure  supplements 6, 12 and 13). For example, we measured differences in biofilm accumulation across the parameter space using optical density measurements in eVOLVER, but these did not correlate with diversity. Similarly, changes in pH, oxygenation, phage, or antimicrobial levels could conceivably be linked to diversity changes and were not measured in this experiment. Notably, modifying the Monod consumer resource model to include asymmetric mortality rates did not qualitatively change the shape of DDRs in simulations (Appendix 1-figure 4), suggesting that our results could be robust to these types of effects. Overall, we found it striking that the model captures the features of our results so well while being relatively simple and non-parameterized.
The consumer-resource model can produce all major classes of DDRs and is robust to alternative formulations We set out to develop a conceptual framework that not only captures our results, but also provides mechanistic intuition about how the characteristics of disturbance, such as frequency, can produce and reshape DDRs more generally. In this case, we demonstrated how a U-shaped DDR emerging  under constantly-applied disturbance was reshaped to a more neutral relationship by the addition of low-frequency fluctuations. This framework further identified a niche-flip mechanism as a dominant factor in community assembly under constantly applied disturbance. We next wondered whether this simple model could produce other classes of DDRs observed in nature (Mackey and Currie, 2001), as a path toward identifying the organizational mechanisms in other ecosystems. To explore a wider range of behaviors, we extended the range of simulations to include extreme disturbance intensities that eventually lead to community extinction, and prevented artificial extinction of all species simultaneously by introducing noise into the normalization of growth rates (creating small differences in relative fitness between species). As depicted in a contour plot (Figure 4), we again found that fluctuations reshaped the DDRs to produce complex diversity landscapes. These simulations revealed that a diverse repertoire of DDRs can emerge from the combination of these two disturbance characteristics (intensity and frequency). By exploring subsets of this parameter space, we observed every major class of DDR: positive, negative, neutral, peaked, and U-shaped (Figure 4). Different behaviors can emerge from similar communities maintained under different portions of the disturbance parameter range, or when similar disturbances are applied to communities with different growth parameters. By simply examining these two disturbance characteristics in a systematic way, we developed a unifying framework that can reconcile disparate DDR observations without invoking more complex phenomena.
To check how well these behaviors were preserved under variations of the Monod model, we performed simulations in which we (1) modestly varied parameter values, (2) significantly changed parameter distribution assumptions, or (3) changed the functional form of the model altogether.
Notably, as detailed below, the U-shaped DDR and frequency-dependence were maintained for nearly all model variations, indicating that disturbance intensity and fluctuations remained important predictive factors of community diversity.
First, we made quantitative changes to parameters that did not qualitatively change the relationship between species. In one example, we explored larger relative fitness differences between species by removing growth rate normalization between species, but the U shape and fluctuationdependent DDR changes were preserved (Appendix 1-figure 3e). In another, we chose parameters in ranges that diverged from experimental conditions and found the model results to be robust to increased simulation length (Appendix 1-figure 3d) and 10-fold increased resource supply (Appendix 1-figure 3b), but not 10-fold decreased resource supply (Appendix 1-figure 3c). At increased number of species and resources, a U-shaped DDR is still observed, although its relative impact on total diversity is reduced (Appendix 1-figure 3i).
Second, we modified parameters in ways that changed the underlying assumptions about community structure. We found no qualitative changes to model results when including 'specialist' species or eliminating r/K tradeoffs with sorted r and K values (Appendix 1-figure 3f and g). However, when we held K constant for all species on each resource, both the U shape and frequency dependence were eliminated (Appendix 1-figure 3h). For a mathematical discussion of the importance of these parameter choices, see Appendix 2.
Finally, we sought to determine whether changing the functional form of the summed Monod model would affect the U shape. First, we employed a non-additive formulation of Monod growth on mixed substrates (Sears and Chesson, 2007); similar U-shaped DDRs and frequency dependence were observed (Appendix 1-figure 4). Next, we incorporated the possibilities of migration or species-specific death rates in different models. These modifications respectively increased or decreased the maximum diversity of the system, but the U-shaped DDR and frequency-dependence were maintained (Appendix 1-figure 4). Again, we found it striking that over a wide range of model variations, disturbance intensity and fluctuations remained important factors shaping the diversity of ecosystems in a predictable manner.

Discussion
In this work, advances in automated continuous culture technology and next-generation sequencing enabled laboratory microbial ecology studies to systematically dissect the role of environmental disturbance (intensity and fluctuation) with fine resolution at scale. We found replicable patterns in composition and diversity of a soil-derived microbial community across different disturbance regimes. Notably, we observed an unexpected U-shaped DDR under constant disturbance, and found that adding fluctuations (in the form of low-frequency discrete disturbance events) increased community diversity and reshaped the DDR. All of these results are well captured by the Monod consumer resource model, which subsequently led us to describe and propose a novel niche-flip mechanism for structuring these ecosystems. Taken together, these combined experimental and modeling results (1) provide new insight into how microbial community assembly depends on environmental conditions and (2) demonstrate a role for environmental fluctuations in promoting diversity. Examining the behavior of the Monod model over a broad range of disturbance intensities and frequencies, we see that diverse classes of DDRs (including increasing, decreasing, and peaked) could emerge when only subsets of the parameter space are sampled (Figure 4). Understanding how disturbance intensity and frequency interact is essential for reconciling disparate observations of DDRs under a single quantitative unifying framework (Mackey and Currie, 2001;Hughes et al., 2007;Hall et al., 2012;Miller et al., 2011).
While we focused on two characteristics, a number of different temporal features of a disturbance could be considered (Mackey and Currie, 2001;Hall et al., 2012;Miller et al., 2011). In our present work, by focusing on mean disturbance intensity rather than maximum disturbance intensity, we have sampled conditions with a negative correlation between frequency and maximum intensity (e.g. comparing frequent small dilution events to infrequent large dilution events). Furthermore, our definition of constant disturbance may in fact represent the undisturbed state of some ecosystems. For natural ecosystems with constant flow, such as lentic aquatic systems, it may be more accurate to consider a drought disturbance with a sudden reduction in flow rate (with distinct biological outcomes). Performing an experiment in which the maximum flow rate, minimum flow rate, duration at maximum and duration at minimum were all independently varied could help to understand a wider range of disturbance profiles.
With further study and evaluation of underlying assumptions, the findings of this work may be extended to other systems. Although our experimental results were reproducible (Figure 3-figure supplements 7, 10 and 11), it remains to be seen whether other species and disturbance types (including asymmetric disturbances like toxins or heat shock) behave similarly in experiments. Our simulation findings indicate that the niche-flip mechanism is quite generalizable across varying assumptions and formulations, and thus would predict similar roles for other disturbance types and conditions. Furthermore, the growth saturation feature of the Monod model that enables niche-flip is also characteristic of Type II and III functional responses used in consumer resource models across different scales of ecology. It remains to be seen whether niche-flip mechanisms could arise from non-resource-based models. It is plausible that similar growth tradeoffs arising in response to other disturbance-correlated features could lead to loss of coexistence at intermediate disturbance intensities. Therefore, niche flip could be a more general principle extending beyond relative growth nonlinearities explored in this work to systems driven by dynamic abiotic stresses and/or storage effects (Sears and Chesson, 2007;Angert et al., 2009;Usinowicz et al., 2017;Hallinen et al., 2020). We also note however, that these types of disturbances do not share the direct link between environmental change and biological outcome that is characteristic of dilution disturbance, so the impact may be less clear.
Broadly, our results highlight that the structure of ecosystems and their response to perturbation is contextual. We demonstrated that increasing the disturbance intensity can increase, decrease, or have no effect on the diversity of a system. Critically however, we found that under a unifying framework that considers both disturbance intensity and frequency, these relationships become predictable rather than idiosyncratic (Hall et al., 2012;Miller et al., 2011). Intriguingly, these complex ecosystem assembly rules can emerge from temporal structure alone, without invoking other organizing principles, such as spatial structure (Gude et al., 2020) or network structure (e.g. cross-feeding and antagonism) (Mee et al., 2014;Kelsic et al., 2015). While we considered both constant and oscillating temporal patterns of disturbance, additional complexity may arise from temporal patterns that are random or those which change which types of resources available over time. Since complexity can arise from temporal, spatial, and network effects, the relative importance of any one organizing factor is likely ecosystem dependent. If predictable response to perturbation depends on context, then designing predictable interventions to ecosystems (in medicine, agriculture, and conservation) will require the ability to measure and understand the environmental context. With the staggering amount of compositional data being generated with high-throughput sequencing of microbiomes (Thompson et al., 2017), inference of environmental context and design of robust ecological interventions may not be far off.

Preparing inoculum
Two g of dirt from the Communications Lawn of Boston University (collected on 09/15/2018) was vortexed in 10 mL PBS + 200 mg/mL cycloheximide, then incubated in the dark at room temperature for 48 hr. For pre-enrichment, 16 eVOLVER vials were prepared with 25 mL of 0.1X Nutrient Broth (NB) media (0.3 g/L yeast extract + 0.5 g/L peptone (Fisherbrand)) with 200 mg/mL cycloheximide, inoculated with 350 uL of PBS immersion, and grown for 18 hr in eVOLVER at 25˚C. All 16 pre-enrichment cultures were mixed together to form the experiment inoculum. Aliquots in 15% glycerol were stored frozen at À80˚C.
Running eVOLVER Experiments eVOLVER lines were sterilized using 10% bleach and ethanol (Heins et al., 2019;Wong et al., 2018), then autoclaved vials were loaded with 23 mL of 0.1X NB. Each vial was inoculated with 1 mL of inoculum, and grown for 5 h at 25˚C with stirring prior to the first dilution disturbance. eVOLVER was operated in chemostat mode with 0.5 mL bolus size, with dilutions either evenly distributed over time (constant disturbance) or concentrated in fluctuation periods lasting 15 minutes. For these cultures, the flow rate during a fluctuation d f depended on the number of fluctuations per day f and mean dilution rate d according to the following equation: At the end of each experiment, vials were flushed with media, and 10 optical density measurements were taken in eVOLVER to measure the biofilm levels.
Bottles and lines were routinely checked for contamination. This occurred to only a single vial of the experiment, which was excluded from statistical analysis. For the follow-up washout experiment, the glycerol stock inoculum was thawed at room temperature, 1 mL was inoculated into each vial, then the cultures were allowed to grow for 5.7 hr prior to initiating disturbances. For the washout experiment, a software bug caused a few incorrectly executed dilution events; these vials were excluded from statistical analysis. Code required to execute these experiments will be available on Github (https://github.com/khalillab/DDR-eVOLVER-model).

Sampling cultures
At each timepoint, a 2 mL culture aliquot was removed from each vial with an extended length pipette tip.
For plating, 20 mL of the sample was used for a 10-fold serial dilution series, and 100 uL of diluted cultures at three concentrations were plated on 18 mL Nutrient Broth Agar plates (3 g/L yeast extract, 5 g/L peptone, 15 g/L agar [Fisherbrand]), which were grown at room temperature for 48-60 hr, then imaged on an on an Epsom Perfection 550 scanner. Image analysis was performed with the aid of Cellprofiler 3.1.8 (Lamprecht et al., 2007) and Cellprofiler Analyst 2.2.1 Classifier (Bray et al., 2015) tools.
For DNA extraction, the remainder of the sample was pelleted and frozen at À20˚C. 60-72 hr after freezing, pellets were lysed at 37˚C for 1 hr in 200 uL of lysozyme buffer (25 mM Tris HCl pH 8.0, 2.5 mM EDTA, 1% Triton X-100 with 20 mg/mL lysozyme (Fisher), prepared fresh daily). Lysates were processed using DNEasy Blood and Tissue Kit according to manufacturer specifications, eluted into 10 mM Tris buffer, and normalized to 5 ng/mL DNA based on measurements in a Qubit fluorometer.

Library preparation and sequencing
Briefly, we performed amplicon sequencing of the 16S v4 region based on established protocols (Gohl et al., 2016). Primers prCM543 (TCGTCGGCAGCGTCAGATGTGTATAAGAGAC-AGGTG YCAGCMGCCGCGGTAA) and prCM544 (GTCTCGTGGGCTCGGAGATGTG-TATAAGAGACAGG-GACTACNVGGGTWTCTAAT), adapted from EMP515F (Apprill et al., 2015) and EMP806R (Parada et al., 2016) were used to isolate a 290 bp 16S v4 region, using Kapa Hifi ReadyMix polymerase and the following cycling conditions: (i) denaturation: 95˚C for 5 min; (ii) amplification (25 cycles): 98˚C for 20 s, 55˚C for 15 s, 72˚C for 1 m; (iii) elongation: 72˚C for 5 min. For the negative control and biofilm samples, the number of cycles was increased to 35 to amplify from low biomass. Illumina NexteraXT primers (or equivalents) were used to form a final library 427 bp in length, with the following conditions: (i) denaturation: 95˚C for 5 min; (ii) amplification (eight cycles): 98˚C for 20 s, 55˚C for 15 s, 72˚C for 1 m; (iii) elongation: 72˚C for 10 min. DNA was purified with AMPure XP beads or SequalPrep plates, then samples were multiplexed in groups of 192 alongside control samples at a higher fraction, and spiked with PhiX or whole genome DNA libraries to a final concentration of 50% to increase sequence diversity. Library pools were sequenced at the Harvard Biopolymers Facility across five 250 bp paired end MiSeq v2 runs.

Sequencing analysis
Samples were demultiplexed using the Illumina BaseSpace demultiplexer analysis tool. All subsequent bioinformatic analysis was performed in QIIME2 v2020.2 (Bolyen et al., 2019). Demultiplexed samples were dereplicated using DADA2 sample inference to tabulate Amplicon Sequencing Variants (ASVs) (Callahan et al., 2016). Next, for qualitative description of composition, taxonomy (to the genus level) was assigned to each feature by alignment to the SILVA 132 database (Quast et al., 2013) using the taxa-barplot plugin. For quantitative analysis, samples with technical issues (e.g. contamination, low biomass, poor sequence quality, etc.) were removed and the remaining 698 samples were rarefied to 6840 reads. The fragment-insertion plugin was used to generate a rooted phylogenetic tree using the SEPP algorithm (Janssen et al., 2018). The diversity plugin was used to calculate Shannon diversity, ASV richness, and weighted UniFrac distance (Lozupone et al., 2011), which was used to perform Principle Coordinate Analysis (PCoA) and PERMANOVA analysis.

Isolating strains
Colonies from imaging plates for the inoculum and endpoint timepoints of the DDR64 and follow-up washout experiments were re-struck on NB agar plates, grown overnight in 0.1X NB media, and stored frozen at À80˚C in 15% glycerol. Primers prCM215 (CCATTGTAGCACGTG-TGTAGCC) and prCM216 (ACTCCTACGGGAGGCAGC) were used to amplify the v3-v7 16S region for Sanger sequencing and identification.

Growth characterization of isolates
Several strains with different taxonomic background were selected from the collection of isolates for further study. For measurement of r/K values, nine strains were struck onto NB plates, and single colonies were grown to stationary phase overnight at 30˚C, then diluted to an OD600 of 0.001 in triplicate 200 uL cultures in dilute NB media (ranging from 0.1X to <0.0016X, plus a water-only negative control) in 96-well microplates grown in a Tecan spectrophotometer for 36 hr at 30˚C, shaking, with lid on. Growth rates were fit to log transformed OD600 data, Lineweaver-Burke plots were constructed across media concentration, and Monod curves were fit to each species. Three strains that could grow in M9 defined minimal media were selected for characterization on different limiting resources. Strains were struck onto NB plates, and single colonies were grown to stationary phase overnight at 30˚C in carbon-limited M9 supplemented with aspartate, glutamate, or proline, then diluted and grown in microplates as above, with amino acid concentrations ranging from 10 mM to <0.16 mM, plus a no-carbon control. Due to limited growth on these media, optical densities at subsaturating resource concentrations were below the limit of detection, so only maximum growth rate at 10 mM amino acid was reported.

Data and materials availability
All sequencing data is deposited in the Sequence Read Archive (SRA) accessible with a BioProject accession code PRJNA719465. Agar plate images are being deposited on Figshare accessible at https://doi.org/10.6084/m9.figshare.15117558. Computer code used to run eVOLVER experiments and for theoretical modeling is available at https://github.com/khalillab/DDR-eVOLVER-model (copy archived at swh:1:rev:7019f031598169723a5828f3909eba1e199794d2; Khalil, 2021). All other datasets required to produce the results in the current study are included as supplemental data.
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Data availability
All sequencing data is deposited in the Sequence Read Archive (SRA) accessible with a BioProject accession code PRJNA719465. Agar plate images are deposited on Figshare accessible at https:// doi.org/10.6084/m9.figshare.15117558. Computer code used to run eVOLVER experiments and for theoretical modeling is available at https://github.com/khalillab/DDR-eVOLVER-model (copy archived at https://archive.softwareheritage.org/swh:1:rev:7019f031598169723a5828f3909eb-a1e199794d2). All other datasets required to produce the results in the current study are included as supplemental data. Source data files have been provided.
The following dataset was generated: N i represents the abundance of species i modified by its carrying capacity (we rescale the carrying capacities of all species to one), r i represents its maximal growth rate, and a ik represents inhibition of species i by species k. The model is parameterized such that a ii ¼ 1.
Simulations mimicked experimental conditions as much as possible. Beginning with equal abundances of all 10 species, we integrated equations for a six-day period using the function ode45 in MATLAB. Species abundances were diluted during equally spaced 15 min intervals by integrating a version of Equation 1 modified to include dilution: d, the death or dilution rate, was calculated by distributing the mean dilution rate D (per hour) over equally spaced intervals of 15 minutes (per 1/4 hour) at a frequency f (per 24 hours): Frequency f of dilution ranged from 1 to 64 per day, and mean dilution rate D ranged from 0.1 to 0.8 per hour. Maximal growth rates r i were randomly sampled from a normal distribution with mean 1 and standard deviation 0.1. This range was selected to roughly match measured growth rates of isolates on 0.1X Nutrient Broth (Figure 3-figure supplement 4). Competition coefficients a ik were randomly sampled from a lognormal distribution with parameters ¼ À0:7 and s ¼ 0:2, which are the mean and standard deviation of the associated normal distribution. The mean of the lognormal distribution is exp þ s 2 2 ¼ 0:51 and the standard deviation is exp 2 ð Þ þ s 2 exp s 2 ð Þ À 1 ð Þ ¼ 0:1. The competition coefficients were selected to match the diversity of the resource-explicit simulation results at the zero-dilution condition (see below). Other ranges of r i and a ik did not alter the DDR shape (Appendix 1-figure 1).
We simulated 100 competitions with randomly drawn parameters, across all dilution-frequency combinations. All 10 species began at equal abundances, and we used their final relative abundances p i to calculate the Shannon diversity index of each outcome: In order to be counted, the abundance of a species p i had to exceed a threshold of 0.0001. Finally, we took the average of the 100 values of for each dilution-frequency combination to obtain average diversity.

Linear consumer resource simulations
We simulated 10-species, seven-resource competitions with a linear consumer resource model: N i represents concentration of species i, r ij its growth rate per unit resource on resource j, and c j the concentration of resource j. Units of species and resources per volume are the same, because we assume that one unit of resource is fully converted into one unit of biomass; this assumption is equivalent to assuming the biomass yield is equal to one for all species, and therefore we do not include a yield parameter.
Similar to the Lotka-Volterra simulations, we integrated the dilution-modified equations over the same time period and range of frequencies and dilutions: d, as defined by Equation 3, is not only the dilution rate of cells, but also the influx of fresh resources at source concentration c jo : Simulations began with equal abundances of all species and resources (c jo ¼ 1), except for one resource, which had a slightly different supply concentration (c jo ¼ 1.2), in order to move the resource supply point away from a unique central position. Growth rates per unit resource r ij were randomly sampled from a normal distribution with mean 1 and standard deviation 0.1, and then normalized by dividing by the sum of growth rates for each species across all resources, P 7 j r ij . Other parameter ranges did not alter the DDR shape (Appendix 1-figure 2).

Monod consumer resource simulations
We simulated 10-species, 7-resource competitions with a Monod resource model: The Monod constant, K ij , is the concentration of resource j at which species i reaches its halfmaximal growth rate on that resource. Similar to the linear resource concentration model (Equation 6 and 7), the model can be modified to include dilution: We performed simulations using the same procedure described above. The Monod constants K ij were randomly sampled from a uniform distribution: [0.001, 0.01]. The width of the selected range is consistent with the range of Monod constants measured on Nutrient Broth (Figure 3-figure supplement 4), and similar DDR shapes are generated with alternative parameter choices (Appendix 1- figure 3). Maximal growth rates r ij were sampled and normalized as described above. To check whether the U-shaped DDR was preserved under variations of the Monod model, we performed other simulations.
First, we introduced larger relative fitness differences between species would affect our results or yield different DDRs. We added Gaussian noise to maximal growth rates r ij with standard deviation of 0.02, which created small differences in relative fitness between species. We observed similar results, with a gradual decrease in the diversity as dilution rates exceed individual species maximal growth rates (Figure 4). For larger relative differences, we removed normalization between species, and instead normalized maximal growth rates per resource by dividing by the number of resources, rather than the sum of growth rates across all resources (Appendix 1-figure 3e).
Diverging from the conditions of the experiment, we searched for alternative parameter ranges that could qualitatively alter the DDR. We found the model to be robust to increased simulation length (Appendix 1-figure 3d) and increased resource supply (Appendix 1-figure 3b). However, when reducing each of the resource supply concentrations, c jo , by a divisive factor of 10, the U-shape was almost completely eliminated except for a small effect at the highest dilution frequency (Appendix 1- figure 3c). This modification is equivalent to increasing Monod constants, K, by a multiplicative factor of 10, placing them closer to the range of resource supply concentration, which does not reflect our measurements (Figure 3-figure supplement 4) and is less biologically relevant. This extreme parameter range essentially changes the model choice, moving the growth dynamics closer to those of the linear consumer resource model.
Next, we made modifications to parameter choices under the summed Monod model that change the underlying assumptions to see if any would eliminate the U shape. First, we included the possibility of "specialist" species by widening the maximal growth rate sampling range. We randomly sampled r from a uniform distribution (0.1,1), before normalizing by dividing by the sum of growth rates for each species across all resources, P 7 j r ij , as before. The U shape was preserved (Appendix 1-figure 3f). Next, we eliminated r/K tradeoffs on all resources by sampling r and K as before, and then sorting them such that the species with the highest value of r also had the lowest value of K, and the next-highest and next-lowest, and so on, for each resource. The U shape was preserved (Appendix 1- figure 3g). Finally, we sampled K only once per resource, keeping it constant for all species on each resource. This modification eliminated the U shape (Appendix 1- figure  3h and Appendix 2-figure 1). For a mathematical discussion of the importance of these parameter choices, see Appendix 2.
Finally, we sought to determine whether changing the functional form of the summed Monod model would affect the U shape. First, we employed a non-additive formulation of Monod growth on mixed substrates (Sears and Chesson, 2007). In this formulation, the per-capita growth rate on multiple resources is a saturating function of Monod growth on individual resources, rather than a simple sum (as in Equation 10 above): Here, f ij ¼ rijcj Kijþcj , which is the Monod growth of species i on resource j. The other parameter, l c , is the horizontal intercept in a plot of a species' catabolic enzyme expression as a function of its growth rate. For simplicity, we assumed l c to be equal for all species. Because the species consume resources non-additively in this model, the resources are also depleted non-additively: We can write the dilution-modified model as follows (substituting l c g i for per-capita growth of species i in the absence of dilution): Simulating this model with values of l c ranging between 0.3 and 30 yielded DDRs similar to those of the summed Monod model, where the U shape was preserved (Appendix 1-figure 4a).
Next, we incorporated the possibility of asymmetric death rates by introducing a species-specific death rate in addition to a uniform dilution rate: The species-dependent death rate d i , randomly sampled from a uniform distribution of [-0.1, 0.1], introduced asymmetry to mortality of each species in the community. During simulations, species with high mortality were severely penalized compared to species with low mortality, and the species pool was effectively reduced to only species with low mortality. This modification led to a lower maximum diversity of the system, but U-shaped DDR was maintained (Appendix 1-figure 4b).
Finally, we incorporated the possibility of migration by introducing a species-specific migration from a common source: The migration rate m i , set as 0.001 for all species, is applied continuously regardless of the fluctuation in disturbance, modeling a constant population influx from a highly diverse common source. This uniform migration boosted the diversity overall, yet the U-shaped DDR was preserved (Appendix 1-figure 4c). Appendix 1-figure 4. Comparison between Monod consumer resource models with variant assumptions can retain U-shape. Each plot displays Shannon diversity after 6 days of simulated growth, mean across 100 communities for varying mean dilution rate and dilution frequency (see legend). Variations on the basic Monod model are discussed in Materials and methods. (a) Growth rates across multiple resources were normalized using an alternative mixed-substrate utilization model, retaining the DDR shapes. (b) Rather than having equal mortality rates across species, each species was assigned an independent mortality rate d i in addition to the universal dilution rate d.
Diversity is reduced as species with higher mortality rates are penalized, but this does not qualitatively change the DDRs observed. (c) To account for migration, each species was assigned a migration rate from a universal source community, independent of the dilution rate. Though species are constantly reintroduced into the community, if the migration rates are low compared to growth and dilution rates, then there are no major qualitative changes to the observed DDRs.
Furthermore, we can extend this condition to overall parameters measured in a complex media, rather than on each component resource. Summed Monod growth on all resources can be approximated by a single Monod growth curve with: To satisfy the condition for avoiding a U-shaped DDR, the community should show a positive correlation between r i and r i =K i : In conclusion, in the summed Monod model, a U-shaped DDR is guaranteed when r and r=K in complex media do not have a positive correlation. To be certain, the strength of the U-shape should depend on the prevalence of tradeoffs, and by extension, a weak correlation would result in a less dramatic U-shape than would a strong negative correlation. Additionally, observing a U-shape resulting from a small number of rare tradeoffs would require measuring the community at the dilution rates where the few instances of niche flip occur. The lack of positive correlation exhibited by our isolates (Figure 3-figure supplement 4), however, indicates that niche flip is a relevant mechanism that warrants further investigation.
Thus far, we derived the required conditions for niche flip in a two-species, two-resource community: r 21 >r 11 ; K 21 >K 11 r 21 r 11 ; r 22 <r 12 ; K 22 <K 12 r 22 r 12 These inequalities represent trade-offs between r and K that cause the Monod growth curves of the two species in both resources to cross, but such that the species that grows faster at low resource concentration flips between the two resources ( Figure 2-figure supplement 1).
As discussed in the main text, these r/K tradeoffs make coexistence possible at low and high disturbances, while at some level of intermediate disturbance the niche flip eliminates the possibility for coexistence. This mechanism could lead to a U-shaped diversity-disturbance relationship among a community of many species consuming many resources, assuming that trade-offs between r and K were prevalent. In the following pages, we show that this requirement can be overcome, as r/K trade-offs across different resources can also lead to niche flip. To illustrate this concept, we must consider an even higher disturbance regime than explored previously.
Let us first check whether eliminating r/K tradeoffs on each resource would eliminate the U-shaped DDR. In a simulation of a 10-species, 7-resource community, we sorted r and K such that the Monod growth curves on each of the seven resources do no intersect, in order to avoid tradeoffs in each of the seven resources (Appendix 2- figure 1a). Surprisingly, this sorting did not eliminate the overall U-shaped DDR (Appendix 2-figure 1b). To explain how a lack of r/K tradeoffs on each resource could nevertheless lead to niche flip and the U shape, we return to the example of a two-species, two-resource community. Let us consider an even higher disturbance regime than explored previously.
Condition 3: Very high dilution regime: d>r ij In this regime, the dilution rate is greater than any species' maximum growth rate on any single resource. In the summed Monod model, species can survive when the sum of growth rates across resources can overcome the dilution rate.
In the two-species, two-resource example, neither species can survive on a single resource, but both could survive in the presence of both resources at a sufficient concentration. The niche of each species is now determined not by the intersection of the ZNGIs and the resource axis, but by the asymptotic distance between the ZNGI and the resource axis. A species that requires less additional resource 2 in order to survive would win in media containing mostly resource 1. For this reason, growth on a second resource now plays a role in determining the winner in the (predominantly) first resource. Under the summed Monod model, the necessary concentration of resource 2 can be solved explicitly: Coexistence Appendix 2-figure 1. Generating species with sorted r and K parameters prevents intersection of growth curves but not niche flip. (a) Growth curves of 10 species on a single resource are shown. For these species in a single resource environment, the outcome of competition does not depend on resource concentration or dilution rate. (b) Communities with sorted K per resource still exhibit U-shaped DDR. Shannon diversity of communities vs. mean dilution rate. Color indicates frequency of fluctuations. Calculations performed after 6 days of simulated growth, and averaged across 100 communities with randomly generated and sorted r/K parameters. (c) Niche flip occurs in a community with sorted K parameters at high dilution rate. At low dilution rate, ZNGIs intersect the axes and the resource consumption vectors outline a coexistence region. At high dilution rate, species cannot survive on a single resource, so ZNGIs do not intersect the axis. The relative positions of ZNGIs and consumption vectors (i.e. invasion boundaries) have flipped at high dilution rate, consistent with niche-flip. Therefore, since r/K trade-offs on separate resources can also lead to niche flip, we must keep K constant for all species grown on any particular resource in order to eliminate the possibility of niche flip. In a simulation of a ten-species, seven-resource community, with one K per resource assigned to all ten species, the diversity-disturbance relationship is not U-shaped (Appendix 1-figure 3).
In conclusion, we observe that when the dilution rate is smaller than the maximal growth rate on each resource, r/K tradeoffs on each resource are required for niche-flip and a U-shaped diversitydisturbance relationship. If we consider an even higher dilution rate, r/K tradeoffs on each resource are no longer required for niche-flip because r/K tradeoffs across different resources can also lead to niche-flip. This opens up additional parameter space that would lead to a U-shape diversity-disturbance relationship.
In a multi-resource multi-species community, this implies that even a positive correlation between r and r/K in complex media does not necessarily guarantee the absence of a U-shaped DDR. Our results show that niche flip and a resulting U-shaped DDR under constant disturbance are difficult to avoid in the Monod model, and that their absence requires a very constrained parameter space. Note however, that in this section we have considered only chemostat continuous dynamics. As shown experimentally and in simulations, introducing fluctuations into the system entirely reshapes the DDRs (Figure 3e and f). Finally, while the simulations are robust to a variety of parameter choices and alternative formulations (Appendix 1- figure 3 and 4), populating the model with growth parameters from different species-media combinations may yield a diverse range of DDR shapes. Nevertheless, we believe the interesting features of the niche-flip mechanism could warrant its study in other ecosystems.