Co-occurring soil bacteria exhibit a robust competitive hierarchy and lack of non-transitive interactions

Microbial communities are typically incredibly diverse, and this diversity is thought to play a key role in community function. However, explaining how this diversity can be maintained is a major challenge in ecology. Temporal fluctuations and spatial structure in the environment likely play a key role, but it has also been suggested that the structure of interactions within the community may act as a stabilizing force for species diversity. In particular, if competitive interactions are non-transitive as in the classic rock-paper-scissors game, they can contribute to the maintenance of species diversity; on the other hand, if they are predominantly hierarchical, any observed diversity must be maintained via other mechanisms. Here, we investigate the network of pairwise competitive interactions in a model community consisting of 20 strains of naturally co-occurring soil bacteria. We find that the interaction network is strongly hierarchical and lacks significant non-transitive motifs, a result that is robust across multiple environments. Moreover, in agreement with recently proposed community assembly rules, the full 20-strain competition resulted in extinction of all but three of the most highly competitive strains, indicating that higher order interactions do not play a major role in structuring this community. The lack of non-transitivity and higher order interactions in vitro indicates that other factors, such as temporal or spatial heterogeneity, must be at play in enabling these strains to coexist in nature.

area of interest in microbial ecology has been understanding the factors that give rise to 39 the diversity observed within microbial communities 4,5 . A better understanding of the 40 structure of microbial communities is desirable for both managing existing microbial 41 communities 6 and, eventually, engineering them de novo 7 . 42 Many factors can contribute to the generation and maintenance of diversity in ecological 43 communities. Non-transitivity 8 , bistability 9 , weak interactions 10 , facilitation, multiple 44 limiting factors, and spatial or temporal segregation 11 have all been hypothesized to play 45 a role; however, there is little empirical data regarding the relative importance of each of 46 these factors in actual natural communities. By investigating the network of underlying 47 interactions among the members of a given community, we can better understand each 48 factor's relative importance in structuring the community 12 . Since interspecific 49 competition is thought to be a dominant factor in determining whether a given species 50 can persist in a community 13,14 , the network of competitive interactions between species 51 may be particularly informative of the structure of the community within which the 52 interaction takes place. Features of competitive interaction networks that could contribute 53 to community diversity can include non-transitive motifs such as the classic rock-paper-54 scissors triad, network modularity 15 , or overall trends towards weak interactions among 55 species 56 While non-transitivity in particular is often cited as a potential driver of interspecies 57 coexistence 16,17,18 , the degree to which it occurs in natural communities remains largely 58 unknown. Indeed, Levine and colleagues recently asserted that despite the theoretical 59 potential of non-transitive interactions to stabilize community structure, there is scant 60 evidence that they are widespread in natural systems, and that further empirical studies 61 are warranted 19 . Recent experimental work using a field-parameterized model of 62 competition in annual plants 20 and naturally co-occurring Streptomyces bacteria 9 suggest 63 that rock-paper-scissors type interactions may be less common in natural communities 64 than we might assume; however, further studies of competitive interaction networks in 65 diverse ecological communities are warranted, particularly among phylogenetically 66 diverse natural assemblages. 67 Here, we add to this small but growing body of research that suggests that non-transitive 68 interactions may play a less significant role in maintaining species diversity than is 69 commonly assumed. We use a model system composed of heterotrophic bacteria isolated 70 from a single soil grain. By competing in all pairwise combinations in laboratory culture, 71 we find that the overarching feature of the resulting interaction network is a strong 72 competitive hierarchy, a feature that is naturally at odds with a high incidence of non-73 transitivity. Therefore, in the natural environment of these bacteria, other factors must be 74 at play that account for their ability to co-occur. 75

76
To probe the network of pairwise interactions in a community of diverse microbes, we 77 isolated a collection 20 strains of naturally co-occurring heterotrophic bacteria from a 78 single grain of soil. This strain collection is phylogenetically diverse and spans 16 species 79 across seven genera and five families ( Fig. 1a and Methods). Similar to ref 21 , we co-80  inoculated all pairwise combinations of the 20 strains at varying initial fractions and  81  propagated them through at least five growth-dilution cycles. During each growth cycle,  82  cells were cultured for 24 hours and then diluted by a factor of 100 into fresh media. The  83  final outcome of competition was determined by plating the cultures on solid agar and  84 counting colonies, which are morphologically distinct ( Fig. 1c and Supplementary Fig.  85 1). Plating results were confirmed via next-generation sequencing for a random subset of 86 the pairs (Supplementary Fig. 6). 87 Pairwise competitions resulted in one of three qualitatively different outcomes: exclusion, 88 coexistence, or bistability (Fig. 2a-c  The strains exhibit a strong competitive hierarchy. Very few strains were able to exclude 101 a strain with a higher competitive score; out of 187 pairwise competitions measured, only 102 five resulted in the lower-ranked strain excluding the higher-ranked one (Fig. 2d). The 103 degree of hierarchy in this interaction network is highly significant when compared to 104 networks with randomized outcomes (p < 10 -19 ; Fig. 2e). To assess whether the 105 hierarchical pattern was specific to a particular environment, we repeated the 106 competitions with subsets of the full 20-strain collection in different growth media and 107 with different dilution rates ( Supplementary Fig. 2). We found that the resultant 108 interaction networks in these different environments were also highly hierarchical, 109 despite changes in which strains were most competitive ( Fig. 2f and Supplementary Fig.  110 3). Thus, we conclude that hierarchy in pairwise competition is a robust feature of this 111 model community. 112 Next, we asked what characteristics of a strain might best predict its performance in 113 competition. We hypothesized that strains that grow well in monoculture will have 114 competitive advantages over strains that grow more poorly. Indeed, we found that 115 exponential growth rate (r) was positively correlated with competitive score (Spearman's 116 rho = 0.77; p < 10 -4 ; Fig. 3a) and that the typical outcome was for the strain with the 117 higher r to exclude the strain with the lower r, which occurred for 67% of pairs (Fig 3c). 118 Carrying capacity (K) in monoculture was less predictive of competitive superiority, but 119 was still significantly correlated (Spearman's rho = 0.55, p < 0.05; Fig. 3b). In general, 120 the likelihood of outcomes other than the stronger grower outcompeting the weaker 121 grower decreases for large differences in r and K ( Supplementary Fig. 4). While 122 differences in these two parameters can be indicators of the likelihood of a given 123 competitive outcome, there are many exceptions, and, indeed, some of the stronger 124 competitors do not necessarily have correspondingly strong single-species growth 125 parameters. Thus, while the each species' intrinsic growth ability correlates with 126 competitive ability, the significant number of exceptions indicates that growth ability 127 alone does not fully explain the hierarchical competitive structure that we observe. 128 An important corollary of the high degree of hierarchy we observed in the interaction 129 network is that non-transitive motifs are vanishingly rare. Non-transitive motifs are 130 instances in which a clear competitive hierarchy among members of a sub-group does not 131 exist, the classic example being a rock-paper-scissors (RPS) triad. Of the 987 triads in our 132 collection for which complete pairwise outcome data are available, only three (0.3%) 133 display the RPS topology. This number is significantly less than is found in randomized 134 networks, where on average 14% of triads were RPS (p < 10 15 ; Fig. 4). Furthermore, the 135 three triads that we classify as RPS each feature strains that display unusually high 136 variability from experiment to experiment, possibly due to rapid evolution, and further 137 efforts to characterize these triads failed to reproduce the non-transitive network 138 topology. As dictated by its hierarchical structure, our network is also highly enriched for 139 perfectly hierarchical feedforward loops, which were observed in over 50% of triads ( Fig.  140 4). Due to the paucity and irreproducibility of observable non-transitive relationships 141 among our strains in vitro, we conclude that such relationships are unlikely to be a 142 significant contributor to their coexistence in a natural environment. 143 Given the hierarchical structure of the pairwise interaction network, we wondered about 144 the potential of higher-order interactions and indirect effects among our strains to give 145 rise to a diverse community. To address this, we inoculated three replicate cultures with 146 equal proportions of all 20 strains and propagated them through five growth-dilution 147 cycles (Fig. 5b). The resulting assemblages were highly replicable, and consisted of three 148 strains representing some of the strongest competitors in pairwise experiments (Fig. 5a,c), 149 all of which were found to coexist with each other in pairwise competition. Notably, this 150 combination of survivors was consistent with the simple community assembly rule we 151 recently developed 21 : namely, that a strain is expected to survive in multispecies 152 competition if and only if it is not excluded by any other surviving species. Since 153 pairwise outcomes alone are sufficient to predict the outcome of multispecies competition 154 in this environment, we conclude that higher-order interactions are unlikely to play a 155 major role in structuring this community. 156

Discussion 157
Many factors can contribute to the generation and maintenance of diversity in ecological 158 communities. Non-transitivity, facilitation, bistability, weak interactions, multiple 159 limiting factors, and spatial or temporal segregation have all been hypothesized to play a 160 role 22 ; however, there is little empirical data regarding the relative importance of each of 161 these factors in actual natural communities. Here, we explored one such community. In 162 this work, we explored the network of pairwise interactions for a community of naturally 163 co-occurring bacteria. Our results indicate that diversity in this community is likely 164 maintained primarily due to factors including and spatial or temporal segregation or 165 multiple limiting factors, rather than frequent bistability, non-transitivity, or higher order 166 interactions, all of which have been hypothesized to play a role in generating and 167 maintaining diversity. Nonetheless, we still do not completely understand the processes 168 that give rise to the diversity we observe in nature. 169 Given that soil is a heterogeneous mixture with a multitude of microhabitats, microbial 170 co-occurrence in soil may be facilitated by niche separation and spatial de-mixing. This 171 would allow the coexistence of strains that display strong inhibitory interactions in well-172 mixed environments. Microbes in soil also experience a strongly fluctuating environment, 173 which can lead to coexistence of multiple strains over time via the soil spore bank. 174 Members of the genus Bacillus are particularly well known for their spore-forming 175 ability, which may allow them to persist in a non-vegetative, and therefore non-176 competitive state, until conditions favor their growth 23 . Finally, our experimental 177 approach clearly requires that the strains to be competed be culturable in the laboratory, 178 so it is impossible for us to exclude the possibility that other strains present within the 179 soil might behave very differently. 180 Simulations of our experimental system using the generalized Lotka-Volterra model 181 (gLV) predicted that, if the underlying ecological interactions among species are assigned 182 at random, the pairwise interaction network should become less hierarchical at lower 183 death rates, corresponding to a lower daily dilution rate in our experimental setup 184 ( Supplementary Fig. 5). In order to test this hypothesis, we competed a subset of pairs 185 while experimentally reducing the dilution rate from 1:100 to 1:10 (Fig. 2f). The 186 hierarchical network structure was robust to this manipulation, and remained highly 187 correlated with growth rates in monoculture. While it is possible that reducing the death 188 rate further could weaken the hierarchy by reducing the importance of a growth rate 189 advantage in determining survival, the most straightforward interpretation of our data is 190 that the hierarchy is not simply due to differences in growth rates. 191 This experimental system also gives us the opportunity to test the importance of higher 192 order interactions in shaping communities. Higher order interactions are said to take 193 place when the presence of an additional species changes the interaction between two 194 existing species 24 , and have the potential to contribute to the maintenance of species 195 diversity 25 . In bacterial systems, this can be driven by complex networks of selective 196 antibiotic production and sensitivity 26 . Despite the potential for higher order interactions 197 in our model community, our simple assembly rule 21 , which disregards higher order 198 interactions entirely, accurately predicted the survivors in all-versus-all competition in 199 vitro, suggesting that higher order interactions are not a major driver of community 200 structure in this instance. 201 The observation of high levels of diversity in communities of competing organisms is a 202 long-standing paradox in community ecology 27  on the type strain with the highest seqmatch score relative to the query strain. Three 228 strains (B. toyonensis 1, 2, and 3) had identical 16S rRNA sequences, and were therefore 229 differentiated using a 404-bp fragment of the pyrE gene amplified using the primers 5′-230 TCGCATCGCATTTATTAGAA-3′ and 5′-CCTGCTTCAAGCTCGTATG-3′ following 231 protocols described in ref 29 . A list of the strains used, their GenBank accession numbers, 232 competitive scores, and inferred growth parameters is given in Supplementary Table 1.  233 For phylogenetic analysis, sequences were aligned using MUSCLE 30 and a tree was 234 constructed using PhyML 3.0 31,32 . 235

Estimation of single-species growth parameters 236
The carrying capacity of each individual strain was estimated to be its optical density at 237 600 nm (OD 600 ) in 0.2X nutrient broth after five repeated growth-dilution cycles, starting 238 from an initial OD 600 of 3×10 -3 . Growth curves at OD 600 were measured in flat-bottomed 239 96-well microtiter plates (BD Biosciences) with lids sealed with Parafilm in a Tecan 240 Infinite M200 Pro plate reader over 48 hr at 25 C with maximum shaking. An 241 approximation of the exponential growth rate of each individual strain was extracted from 242 the growth curves using the time each strain took to reach a threshold optical density. The 243 time-to-threshold method was chosen over other estimates of growth rate due to wide 244 variations in growth patterns across the strains, which led to difficulties in fitting 245 parameters to other population growth models. 246

Competition experiments 247
Prior to competition experiments, cells were streaked out on nutrient agar plates, grown 248 for 48 hr at room temperature, and then stored at 4 C for up to two weeks. Single 249 colonies were picked from these plates and grown for 24 hr at room temperature in 0.2X 250 nutrient broth. 251 cycle, the cultures were thoroughly mixed and then diluted by a factor of 100 into fresh 258 medium. OD 600 was measured at the end of each cycle, and final species fractions were 259 estimated after five (or, in the case of initially low plating density, seven) cycles. 260 To measure the final species fractions, the co-cultures were diluted by a factor of 10 4 -10 6 261 (depending on OD 600 ) in PBS. Seventy-five L of the diluent was plated onto 10 cm 262 Petri dishes containing 25 mL of nutrient agar and incubated at room temperature for 48 263 hr. All but a small fraction of the strain pairs have distinct colony morphologies, so 264 species fractions were estimated by counting colonies of each type (median: 51 colonies 265 per plate). Next-generation sequencing of a subset of the co-cultures affirmed the overall 266 accuracy of the plating technique ( Supplementary Fig. 6). 267

Determining the outcome of competition 268
The result of competition was classified as one of three outcomes: exclusion of a single 269 strain, coexistence of both strains, or bistability. A strain was said to exclude its 270 competitor if it was the sole strain observed from both starting frequencies after 5 cycles, 271 or if it excluded its competitor when starting from an initial frequency of 0.95 and 272 achieved a frequency of 0.85 or greater when starting from an initial frequency of 0.05. 273 Pairs were considered bistable if the strain that started out at a frequency of 0.95 excluded 274 the competitor. All other outcomes were classified as coexistence. 275 Calculating competitive score and network hierarchy score 276 The competitive score s i of each strain i was defined as its mean fraction f ij after co-277 culture with each of the n -1 competitor strains: 278 The hierarchy score (h) for an n-member network is calculated as: 279 The network hierarchy score for the observed set of competitive outcomes was then 280 compared against the distribution of scores for 10,000 simulated networks in which each 281 pair was randomly assigned an outcome of exclusion, coexistence, or bistability with 282 probability proportional to the incidence of each outcome in the empirical dataset. The 283 resulting distribution of hierarchy scores was approximated using the normal distribution 284 to determine p-values. 285 Identifying network motifs 286 The frequencies of distinct topologies among the The data that support the findings of this study are available from the corresponding 295 authors upon reasonable request. An implementation of the routine for estimating the 296 distribution of hierarchy scores and motifs in randomized networks is also available upon 297 reasonable request. parentheses. Error bars represent +/-1 s.d. Differences in observed versus randomized 329 scores were were highly significant in all environments (p < 10 -7 ). 330 331