Host Gut Motility and Bacterial Competition Drive Instability in a Model Intestinal Microbiota

The gut microbiota is a complex consortium of microorganisms with the ability to influence important aspects of host health and development. Harnessing this ‘microbial organ’ for biomedical applications requires clarifying the degree to which host and bacterial factors act alone or in combination to govern the stability of specific lineages. To address this we combined bacteriological manipulation and light sheet fluorescence microscopy to monitor the dynamics of a defined two-species microbiota within the vertebrate gut. We observed that the interplay between each population and the gut environment produced distinct spatiotemporal patterns. Consequently, one species dominates while the other experiences dramatic collapses that are well fit by a stochastic mathematical model. Modeling revealed that bacterial competition could only partially explain the observed phenomena, suggesting that a host factor is also important in shaping the community. We hypothesized the host determinant to be gut motility, and tested this mechanism by measuring colonization in hosts with enteric nervous system dysfunction due to mutation in the Hirschsprung disease locus ret. In mutant hosts we found reduced gut motility and, confirming our hypothesis, robust coexistence of both bacterial species. This study provides evidence that host-mediated spatial structuring and stochastic perturbation of communities along with bacterial competition drives population dynamics within the gut. In addition, this work highlights the capacity of the enteric nervous system to affect stability of gut microbiota constituents, demonstrating that the ‘gut-brain axis’ is bidirectional. Ultimately, these findings will help inform disease mitigation strategies focused on engineering the intestinal ecosystem.


Abstract 30
The gut microbiota is a complex consortium of microorganisms with the ability to 31 influence important aspects of host health and development. Harnessing this 'microbial 32 organ' for biomedical applications requires clarifying the degree to which host and 33 bacterial factors act alone or in combination to govern the stability of specific lineages. 34 To address this we combined bacteriological manipulation and light sheet fluorescence 35 microscopy to monitor the dynamics of a defined two-species microbiota within the 36 vertebrate gut. We observed that the interplay between each population and the gut 37 environment produced distinct spatiotemporal patterns. Consequently, one species 38 dominates while the other experiences dramatic collapses that are well fit by a 39 stochastic mathematical model. Modeling revealed that bacterial competition could only 40 partially explain the observed phenomena, suggesting that a host factor is also 41 important in shaping the community. We hypothesized the host determinant to be gut 42 motility, and tested this mechanism by measuring colonization in hosts with enteric 43 nervous system dysfunction due to mutation in the Hirschsprung disease locus ret. In 44 mutant hosts we found reduced gut motility and, confirming our hypothesis, robust 45 coexistence of both bacterial species. This study provides evidence that host-mediated 46 spatial structuring and stochastic perturbation of communities along with bacterial 47 competition drives population dynamics within the gut. In addition, this work highlights 48 the capacity of the enteric nervous system to affect stability of gut microbiota 49 constituents, demonstrating that the 'gut-brain axis' is bidirectional. Ultimately, these 50 findings will help inform disease mitigation strategies focused on engineering the 51 intestinal ecosystem. 52

INTRODUCTION 53
Trillions of microbial cells representing hundreds of species make up the human 54 intestinal microbiota. This multispecies symbiont supports activities as diverse as host 55 development, nutrient acquisition, immune system education, neural function, and 56 defense against pathogens [1][2][3][4][5]. Changes in microbiota diversity and functional 57 composition have been linked with a variety of human disorders, including obesity, 58 colon cancer, opportunistic infection, and inflammatory bowel disease [6,7]. A major 59 goal of host-microbe systems biology is to clarify the ecological factors that determine 60 microbiota integrity by meshing experimental techniques and quantitative modeling. 61 Insights derived from such efforts will inspire the design of novel therapeutic strategies 62 for microbiota-associated diseases. 63 An unresolved question is whether the host and microbiota function 64 independently or together to govern the dynamics and stability of individual bacterial 65 lineages. Addressing this requires identifying the interactions that arise within the 66 spatially complex and heterogeneous environment of the vertebrate gut. However, 67 progress toward this goal has been hindered due to the technical limitations associated 68 with directly observing intestinal communities. Typical interrogation of vertebrate 69 intestinal microbiota involves phylogenetic profiling of fecal material using high-70 throughput sequencing of 16S ribosomal RNA (rRNA) genes. This technique is blind to 71 the spatial structure of microbial communities, which is known in general to strongly 72 influence interactions [8-10] and has recently been predicted to be important for 73 microbiota stability [11]. Sequencing-based studies also have low sensitivity to temporal 74 structure owing to both experimental and analytic constraints. Experimentally, 75 nervous system (ENS) via a mutation in ret, a gene locus associated with human 122 Hirschsprung disease (OMIM 164761), which stabilized the Aeromonas population and 123 neutralized competition with Vibrio. This work reveals a synergy between bacterial 124 competition and host-mediated spatial structuring of microbiota in determining 125 population dynamics and stability-a feature that is likely mirrored in more complex 126 host-microbe systems such as the human gut. 127 periods. Statistical significance of Aeromonas abundances after Vibrio challenge 153 compared to respective mono-association reference populations (i.e. '5-6' vs. '4-6'; '6-7' 154 vs. '4-7'; '5-7' vs. '4-7') was determined by an unpaired t-test. (C) Time course analysis 155 of Aeromonas and Vibrio abundances determined by dissection and plating at three-156 hour intervals over a 12-hour period starting at 6 dpf. Additionally plotted are an 157 Aeromonas mono-association reference population and 24 hour Aeromonas and Vibrio 158 populations previously plotted in 1B ('4-6' and '6-7', respectively). Statistical significance 159 of Aeromonas abundances to the mono-association reference population (ref.) was 160 determined by an unpaired t-test. CFU=colony-forming units; ***=p<0.0001; ns=not 161 significant; N>19/condition. Gray and black dashed lines in panels B and C denote limits 162 of quantification and detection, respectively. 163

164
We first enumerated total bacterial abundance by gut dissection and standard 165 plating techniques. In mono-associations beginning at 4 dpf and extending 24, 48, or 72 166 hours Aeromonas populations consistently reach 10 4 colony-forming units (CFU) per 167 host ( Fig 1B). In contrast, challenge of established Aeromonas populations with Vibrio 168 over various 24 or 48 hour temporal windows leads to dramatically lower Aeromonas 169 abundance as well as increased host-to-host variation and frequent extinction events 170 ( Fig 1B). Under these conditions Vibrio exhibits only modest deviations in abundance 171 compared to mono-association (S1A and S1B Fig indicate that subsequent waves of colonizing bacteria alone do not account for the 175 observed competitive interaction and that it is inter-specific in nature. We also verified 176 that competition between Aeromonas and Vibrio is dependent on the host environment, 177 as there was not an appreciable difference between their abundances during in vitro 178 mono-and co-culture experiments (S1E Fig). To determine if the abundance of Vibrio 179 correlates with a reduction of Aeromonas populations in vivo, we performed a time 180 course experiment in which zebrafish were sacrificed and assayed every 3 hours for 12 181 hours after inoculation with Vibrio. We found that Vibrio rapidly infiltrates Aeromonas-182 colonized intestines and steadily increases in number over the 12-hour assay period 183 ( Fig 1C). Surprisingly, we did not detect a concomitant decline in Aeromonas, implying 184 that Aeromonas populations do not merely respond proportionally to the abundance of 185 Vibrio ( Fig 1C).   We detected clear instances of Aeromonas population collapses under these conditions 254 (Fig 3D and 3E). However, in contrast to Vibrio-associated collapses, Aeromonas was 255 found to consistently recover from collapses during mono-association. Additionally, 256 Aeromonas collapse events associated with mono-association were smaller and more 257 uniform. Defining collapses as events in which the population decreases by at least a 258 factor of ten within one hour, and assigning their magnitude f as the ratio of the 259 population after collapse to that before, we find that for Aeromonas challenged by 260 Vibrio, log 10 (f) = -1.9 ± 1.0 (mean ± std. dev.) (Fig 3F and 3G). The ratio of the Vibrio 261 population before and after Aeromonas collapse events within the same fish is 262 approximately 1 (Fig 3F), corroborating observations from imaging and plating data that 263 Vibrio is resistant to the perturbations that affect Aeromonas. We found that during 264 Aeromonas mono-associations the magnitude of collapse events was about half of that 265 observed in the presence of Vibrio, log 10 (f) = -1.6 ± 0.4 ( Fig 3G). We also detected a rates (f and p c ). We also note that during challenge experiments, the gross spatial 315 distribution of Vibrio is similar to its distribution during mono-association, while there is 316 considerable broadening in the spatial distribution of Aeromonas when challenged (S5 317 Finally, a conceptually minimal model of interaction is that with Vibrio present, the 318 resources available to Aeromonas post-collapse are less than with Vibrio absent, 319 thereby placing a limit on its potential for recovery. We assess this possibility 320 quantitatively below by estimating carrying capacities, and we also examine the 321 synergistic consequences of changes to the carrying capacity and collapse properties. 322 323 Synergy between stochastic collapse events and competition with Vibrio 324 underlies Aeromonas population dynamics. 325 Thus far our data suggest that Aeromonas is susceptible to stochastic 326 disturbances mediated by host intestinal motility, and that its recovery from these 327 disturbances is altered in the presence of Vibrio. If this is the case, we should be able to 328 build a quantitative model that reflects these data, explains the high variance observed 329 in plating assays (Fig 1), and offers insights into the differential outcomes between 330 mono-association and challenge experiments. The model we construct is illustrated 331 schematically in Fig 5. Consider a bacterial species exhibiting logistic growth, with 332 growth rate r and carrying capacity K (Fig 5A and 5B); in other words, the population N 333 grows with time t according to: 334 Superimposed on this are rare collapses, during which the population drops to f times 335 its pre-collapse value, where f is between 0 and 1, and after which it resumes logistic 336 growth ( Fig 5C). The collapses are stochastic and modeled as Poisson processes; i.e. 337 they occur at random with some probability per unit time p c (Fig 5D). This model arises 338 in many ecological contexts, and some of its mathematical properties have been 339 explored in various studies [28]. Of course, this model incorporates stochastic 340 population collapses by construction, and so does not predict them from first principles. 341 However, the parameter values that emerge from fitting such a model to the data can, 342 as shown below, yield quantitative insights into the mechanisms underlying the 343 observed dynamics that are not evident from mere visual inspection of the raw data. Simulating ensembles of populations that exhibit the above dynamics, we 354 examine the mean and, importantly, the standard deviation of the population at discrete 355 terminal time points, as these are statistics that allow direct comparison to results from 356 plating assays. As shown in detail in the Supporting Text (S1 Text), the apparent 357 dependence of the model on the parameters, r, K, p c , and f, collapses to two effective 358 parameters. The growth rate, r, is both independently known and irrelevant for the 359 conditions considered, and the dynamics depend on the combination z = -p c log 10 (f) 360 rather than on p c and f independently. Values of the two remaining relevant parameters, 361 K and z, which characterize the carrying capacity and the collapse dynamics, 362 respectively, determine the model predictions for the mean and variance of populations. 363 A grid search through the (K, z) space for the values that minimize the distance 364 between the predicted and observed Aeromonas population statistics gives the best-fit 365 model parameters. Additional details and discussion are provided in the Supporting Text 366 (S1 Text). It is important to note that because our imaging data revealed that 367 Aeromonas is often in a state of experiencing or recovering from collapse events, the 368 observed population is likely never close to K, and thus we cannot simply use the mean 369 of the bacterial abundance to estimate K. Rather, we must use a model to infer the 370 carrying capacity that would yield the observed populations. 371 Using Aeromonas abundance data obtained by gut dissection and plating 24 372 hours post Vibrio challenge (Fig 1B, '6-7'), we find best-fit parameters log 10 (K) = 3.2 ± 373 0.5 and z = 0.13 ± 0.05 hr. -1 , the latter providing a constraint on p c and f together. We 374 can independently estimate p c and f from imaging-derived data (Fig 3). As noted 375 previously, for Aeromonas challenged by Vibrio, we find p c = 0.10 ± 0.03 hr. -1 and 376 log 10 (f) = -1.9 ± 0.3 (mean ± std. error), yielding z = 0.19 ± 0.06 hr. -1 , which is consistent 377 with the plating-derived value. The agreement between the separately determined 378 measures of z is remarkable, as it indicates that the statistical properties inferred from 379 an ensemble of populations at a discrete time point are consistent with the properties 380 inferred from the temporal dynamics within individual hosts. As expected, log 10 (K) is 381 greater than the observed mean Aeromonas abundance at 24 hours, since the model-382 derived K represents an upper bound for the population in the absence of any 383 stochastic collapses or Vibrio competition. As another test of consistency, we note that 384 simulating our stochastic model for 48 hours post-challenge using the best-fit 385 parameters determined from plating experiments 24 hours post-challenge predicts a 386 mean and standard deviation of log 10 (population+1) of 1.3 ± 0.3 and 1.5 ± 0.2, 387 respectively, in agreement with the observed plating values of 1.7 ± 0.3 and 1.6 ± 0.3 388 ( Fig 1B). All of these assessments support the conclusion that the observed population 389 dynamics are governed by a mechanism of stochastic collapse. 390 We can also apply this model to Aeromonas mono-association data. Here, the 391 variance of the plating-derived populations is small (Fig 1B), likely due to comparatively 392 rare and/or weak collapses as discussed earlier. For reasons described in detail in the 393 Supporting Text, this hinders robust determination of z, though K remains well fit. We 394 find that z = 0.01 ± 0.01 hr. -1 and log 10 (K) = 4.2 ± 0.1. From live imaging data, p c = 0.04 395 ± 0.01 hr. -1 and log 10 (f) = -1.6 ± 0.2, from which z = 0.06 ± 0.02 hr. -1 (S1 Text). Our 396 identification of thresholds is, by construction, only sensitive to collapses of a factor of 397 10 or more in magnitude (i.e. log 10 (f) ≤ -1), so our estimate of f, and therefore z, is 398 biased toward larger values. 399 The above analysis yields insights into the nature of the competition between 400 Aeromonas and Vibrio that are not obvious from simple visual inspection of the data. 401 The carrying capacity (K) experienced by Aeromonas, as estimated by our model, is 402 only one order of magnitude lower in the presence of Vibrio (log 10 (K) = 3.2 ± 0.5) than 403 when Vibrio is absent (log 10 (K) = 4.2 ± 0.1). However, the observed abundance of 404 Aeromonas is suppressed by more than two orders of magnitude: 405 mean(log 10 (population+1)) = 1.7 ± 0.3 and 4.1 ± 0.1 when challenged by Vibrio and 406 mono-associated, respectively ( Fig 1B). These results suggest that the combined effect 407 of stochastic collapses, which are likely driven by the host environment, and a reduced 408 carrying capacity, as a result of competition with Vibrio, has a far greater influence on 409 population dynamics than either mechanism would provide alone. Aeromonas populations would be more stable despite the presence of Vibrio. To test 418 this hypothesis, we carried out succession assays in mutant zebrafish hosts essentially 419 lacking a functional enteric nervous system (ENS) because of disruption of the gene 420 encoding the Ret tyrosine kinase, which is critical for ENS development [29]. Using DIC 421 microscopy to assess intestinal motility, we found that ret mutant larvae (ret -/-) still 422 exhibit rhythmic contractions, but with different characteristics than wild-type (ret +/+ ) and 423 heterozygous siblings (ret +/-) (S8 and S9 Movie). Because we observed that ret +/+ and 424 ret +/animals are phenotypically similarity with regard to gut motility and that the 425 ret1 hu2846 mutant allele is recessive we further designate ret +/+ and ret +/as 'wild type'. 426 Computational analysis of time-series DIC images allows quantification of the 427 displacement of intestinal tissue during contractile waves (Methods). The average peak 428 amplitude of longitudinal contractions is greater in wild-type than in ret mutant larvae, 429 and in both genotypes declines with age (Fig 6A). At 6 dpf a considerable fraction of ret 430 mutant larvae show low amplitudes, similar to the quiescent state observed in both 431 genotypes at 7 dpf (Fig 6A). Though the amplitude of intestinal contractions might not 432 be directly related to the magnitude or rate of Aeromonas collapse events, it is 433 reasonable to expect some monotonic relationship between the two, as they both reflect 434 intestinal activity. Therefore, we would expect to observe stabilization of  challenged Aeromonas populations in ret mutant hosts only during challenge periods 436 starting at 6 dpf when the difference in intestinal motility between the genotypes is 437 greatest. Indeed, Vibrio challenge of established Aeromonas populations between 5 and 438 6 dpf yielded the same decrease in Aeromonas abundance in both ret mutant hosts and 439 wild types (Fig 6B). In contrast, Aeromonas populations were significantly stabilized 440 during Vibrio challenge from 6 to 7 dpf in ret mutant hosts and in fact were statistically 441 indistinguishable from a reference Aeromonas mono-association (Fig 6B). These results 442 provide strong evidence that ENS-driven intestinal motility contributes to the shaping of 443 this model two-member community by facilitating their apparent competitive interaction. For these reasons we set out to investigate bacterial population dynamics within 478 the vertebrate intestine using a combination of absolute abundance measurements, 479 time-series imaging, and quantitative modeling. Though our system is minimal, 480 consisting of two bacterial species and a larval zebrafish host, it has revealed factors we 481 expect to be of broad relevance to other animal-associated microbiota. Most notably, in 482

this model system the emergence of the apparent competition between Aeromonas and 483
Vibrio is driven in large part by the physical activity of the host, namely the motility of the 484 intestine. In mutant zebrafish hosts that have reduced intestinal motility due to mutation 485 of the gene ret, which impairs ENS development and function, competition between 486 these bacterial species is offset (Fig 6). Motility and mass transport are, of course, key Hirschsprung disease, which is a gastrointestinal motility disorder commonly associated 496 with mutation of ret, have been found to harbor dysbiotic microbial communities [40,41]. 497 The differential susceptibility of our two model bacterial species to intestinal 498 motility can be explained by their distinct community architectures. Highly motile Vibrio 499 are relatively unaffected by intestinal contractions, which is in contrast to the large, non-500 motile aggregates of Aeromonas (Fig 4). Earlier observations of a related A. veronii 501 strain showed higher growth rates for aggregated bacteria compared to planktonic [21], 502 suggesting a tradeoff between enhanced growth and resistance to population level 503 perturbations. In general, we suspect that the spatial structure of microbial communities 504 within the intestine will be an important determinant of their dynamics and a key 505 consideration for the generation of successful predictive models. 506 We are able to construct a quantitative model of the observed Aeromonas 507 dynamics that consists of growth punctuated by stochastic collapses. Data derived from 508 gut dissection, in which many fish are sampled at a single time point, can be fit to the 509 model to determine its two relevant parameters, the bacterial carrying capacity, K, and a 510 factor that characterizes the collapses, z. In itself, this is trivial. However, we can also 511 determine z from independent, and quite different, data, namely image-based time-512 series of individual fish. The two measures agree, which provides strong support for the 513 proposed stochastic-collapse-driven model of inter-species competition. Furthermore, 514 the fit of the model to the data reveals that the impact of Vibrio on Aeromonas 515 populations is twofold: reducing the overall carrying capacity and increasing sensitivity 516 to physical perturbations-the combined effect of the two being much greater than 517 either alone. More generally, our analysis provides evidence that quantitative, data-518 based models of interactions among species within the gut are possible, and that 519 stochastic, rather than purely deterministic, dynamics can play a major role in shaping 520  From an ecological perspective, it is unsurprising that the physical environment 526 and stochastic perturbations influence species abundance; these concepts are 527 mainstays of our understanding of macroscopic multi-species communities [43]. A rich 528 literature describes various stochastic population models and the characteristics, such 529 as extinction probabilities, that emerge from them [44-47]. As shown here, it is likely that 530 such models will in general be useful for providing a conceptual and predictive 531 framework for understanding inter-species bacterial competition. Again mirroring well-532 established ecological concepts, we can frame our understanding of Vibrio and 533 Aeromonas dynamics in the intestine as a study of these species' differential resistance 534 and resilience to environmental perturbations. Aeromonas during mono-association is 535 not resistant to disturbances related to intestinal motility, but it is resilient, able to grow 536 to high abundances despite sporadic collapses. Vibrio, in contrast, is highly resistant to 537 perturbations; it shows smooth growth unfazed by the environmental perturbations that 538 affect Aeromonas (Fig 3). In the presence of Vibrio, both the resistance and resilience of 539 Aeromonas are compromised, as the magnitude of collapses is greater and the carrying 540 capacity to which to recover is diminished. 541 While ecological concepts can help us characterize microbial dynamics, data on 542 microbial systems can, conversely, enhance our understanding of ecological theory. the challenges of performing high-precision field studies. We expect, therefore, that data 550 of the sort presented here, which yield collapse statistics as well as fits to stochastic 551 models, will have utility in contexts far removed from microbiota research. 552 Aeromonas population collapses are well described by stochastic dynamics, but 553 the underlying mechanism by which Vibrio compromises resistance and resilience of 554 The combination of gnotobiotic manipulation and imaging-based analyses can be 566 further elaborated in larval zebrafish, both by increasing the diversity of monitored 567 microbial species and by examining interactions with particular aspects of the host such 568 as its immune system [25]. As illustrated here, we expect that such studies will yield 569 additional insights into the factors that drive the dynamics of complex, natural host-570 associated microbiota. 571

Ethics statement 573
All experiments with zebrafish were done in accordance with protocols approved by the 574

Imaging experiments 651
Sample mounting is done as previously described [21]. Larval zebrafish were removed 652 from culture flasks and anaesthetized using 120 µg/ml tricaine methanesulfonate 653 (Western Chemical, Ferndale, WA). Individual specimens were then briefly immersed in 654 0.5% agar (maximum temperature: 42° C) and drawn into a glass capillary, which was 655 then mounted onto a sample holder. The agar-embedded specimens were partially 656 extruded from the capillary so that the excitation and emission optical paths did not pass 657 through glass interfaces. The specimen holder can hold up to six samples, all of which 658 are immersed EM maintained at 28°C, with tricaine present as an anaesthetic. All long-659 term imaging experiments were done overnight, beginning in the late afternoon.

Measuring intestinal motility 679
Larval intestinal motility was assessed from images captured using differential 680 interference contrast (DIC) microscopy, performed as previously described [27]. The 681 displacement field from frame to frame in time-series was determined using particle 682 image velocimetry (PIV) algorithms [60], which calculate the motions necessary for 683 regions in one frame to be mapped onto regions in another. We focused our analysis on 684 the frequency and amplitude of these motions, restricting our analysis to components of 685 displacement along the intestinal axis. Fourier spectra of the displacements, averaged 686 over location in the intestine, yielded in all cases a clear peak whose frequency and 687 magnitude are indicative of the characteristic frequency and amplitude of intestinal 688 motility, respectively. This method is described in greater detail in a forthcoming paper. 689

ACKNOWLEDGEMENTS 690
We thank Michael Taormina    We describe a simple model of growth and collapse behavior and examine its 909 predictions for population sizes. We also fit the model to experimental data on bacterial 910 abundance. 911 912 1 The model 913 914 Consider a species with population N at time t that exhibits logistic growth, with 915 growth rate r and carrying capacity K: 916 We superimpose on these dynamics events in which the population collapses to a value 918 f times its pre-collapse value, where f is between 0 and 1, after which it resumes logistic 919 growth. We model the timing of the collapses as a Poisson process: collapses are 920 uncorrelated and stochastic, occurring with a probability per unit time p c . Formally, one 921 can write this as a stochastic differential equation: 922 where dM is a Poisson process of unit step. (In other words, dM =1 with probability p c dt, 924 and dM = 0 with probability 1 -p c dt.) N dM refers to N immediately before the collapse. 925 An illustration of the roles of the parameters r, K, f, and p c is provided in Figure 5. As 926 noted in the main text, this model is not new; it has been invoked and studied in many 927 ecological contexts [S1]. However, the particular treatment presented here is, to the 928 best of our knowledge, novel, especially with respect to determining relevant 929 parameters for fits to experimental data. We determine statistical properties of the 930 model using numerical simulations. For infinite carrying capacity, these properties can 931 be calculated analytically, but for the biologically relevant case of finite carrying 932 capacity, exact solutions do not at present exist. 933 934 2 Simulations 935 936 The model described above is simple to simulate by numerical integration, which 937 yields the population x t at time t. Two typical x t are shown in Figure ST1, with 938 parameters as noted in the caption. The model has four parameters, r, K, f, and p c , and a boundary condition set by x 0 (the 949 initial population). The value of x 0 is irrelevant for the experimental conditions 950 considered: the populations start from a small value and grow rapidly. In our simulations 951 x 0 is taken to be 10. 952 The growth rate, r, is known from measurements. Moreover, the model dynamics 953 are fairly insensitive to r, since the experimental timescales of ~10 hours are 954 considerably larger than the timescale set by the growth rate (1/r ~ 1 hour). 955 The key determinants of the population statistics, therefore, are the collapse 956 properties (p c and f) and the carrying capacity, K. The carrying capacity may exhibit 957 considerable variation between fish. Typically, the final populations of Aeromonas in 958 mono-associations are found to be approximately log-normally distributed ( Figure ST2), 959 as is commonly the case for species abundances, and so in simulations we draw K from 960 log-normal distributions. In other words, log 10 (K) for a given simulation is drawn from a 961 Gaussian distribution with some mean value and standard deviation σ K , where σ K is 962 typically 0.5, discussed further below. We note that in the absence of collapse (e.g. p c = 963 0 or f=1) this model is completely deterministic, and the variance in final bacterial 964 populations between fish is solely due to the variance in K. 965 For particular parameter values, we simulate many instances of the above 966 dynamics (typically 1,000 to 10,000) and examine the statistical properties of the final 967 population, x t , assessed at t = 24 hours. For the values used in Figure ST1 above, for 968 example, the mean and standard deviation of the final x t are (6.4 ± 11.5) × 10 3 . The 969 distributions span orders of magnitude, including zero, so it is useful to consider the 970 mean and standard deviation of log 10 (x t +1), similar to a geometric mean. For these 971 parameters, this gives a mean and standard deviation of log 10 (x t +1) of 2.8 ± 1.5. We will 972 define y as 973 y = log 10 (x t + 1) , We can vary the model parameters to determine the relationship between the 986 mean and the variance of the final population, which will allow direct comparison 987 between our model and measurements of bacterial abundance (e.g. Figure 1). The 988 dependence of the mean and standard deviation (std.) of y on p c and f is plotted in 989 Figure ST3. We can intuitively understand its behavior: for small p c or f near 1, the 990 properties of x t are largely set by the mean and variance of the carrying capacity. 991 However, as p c increases (or f decreases), the mean of x t decreases, because larger 992 collapses are more likely to occur, and the standard deviation of x t increases, because 993 the stochastic collapses play a more significant role in the dynamics. For still larger p c 994 (or smaller f), the final population becomes more uniformly small, because the 995 population is dominated by very frequent collapses and cannot grow appreciably. 996 Treating f as a random, rather than a fixed, parameter has little effect on the 997 behavior of the model. Drawing f from a beta distribution, chosen because it is 998 continuous, spans [0, 1], and has two parameters that can be mapped onto a mean and 999 variance, gives the curve shown in Figure ST4. The mean collapse magnitude is chosen 1000 over the same range as f in Figure ST3, and for each mean f, several f values are 1001 drawn from a beta distribution with standard deviations relative to the mean spanning [0, 1002 0.8]. All the resulting population characteristics are plotted in Figure ST4; the resulting 1003 curve is nearly identical to that of Figure ST3. 1004 1005 Figure ST3. The mean and standard deviation of simulated populations at t = 24 hrs., with r = 0.8 hr. -1 and K drawn from a lognormal distribution with mean 10 4 and a standard deviation of half a decade. Blue crosses: p c is fixed at = 0.1 hr. -1 , and f varies between 10 -4 and 10 -0.3 . Red circles: f is fixed at 10 -2 and p c varies between 10 -2 and 10 -0.7 hr. -1 . Each point is calculated from 10,000 simulated runs. 1006 Figure ST4. The mean and standard deviation of simulated populations at t = 24 hrs., with r = 1 hr. -1 and K drawn from a lognormal distribution with mean 10 4 and a standard deviation of half a decade. The collapse probability p c is fixed at = 0.1 hr. -1 , and f is drawn from a beta distribution with mean between 10 -4 and 10 -0.3 , and relative standard deviation between 0 and 80%. Each point is calculated from 1,000 simulated runs. Figure ST7. Contours of χ 2 , the squared distance between simulated (mean(y), std(y)) and the measured value from monoassociation experiments (4.1, 0.6), for a range of z and K, with K drawn from log-normal distributions of width 0.1 decades.

S9 Movie. Example of intestinal motility in a ret mutant larval zebrafish. Differential 1184
interference contrast (DIC) microscopy video of intestinal motility in a conventionally 1185 raised 6 dpf ret mutant larval zebrafish. Scale bar: 50 µm. 1186