Steering ecological-evolutionary dynamics during artificial selection of microbial communities

Microbial communities often perform important functions that arise from interactions among member species. Community functions can be improved via artificial selection: Many communities are repeatedly grown, mutations arise, and communities with the highest desired function are chosen to reproduce where each is partitioned into multiple offspring communities for the next cycle. Since selection efficacy is often unimpressive in published experiments and since multiple experimental parameters need to be tuned, we sought to use computer simulations to learn how to design effective selection strategies. We simulated community selection to improve a community function that requires two species and imposes a fitness cost on one of the species. This simplified case allowed us to distill community function down to two fundamental and orthogonal components: a heritable determinant and a nonheritable determinant. We then visualize a “community function landscape” relating community function to these two determinants, and demonstrate that the evolutionary trajectory on the landscape is restricted along a path designated by ecological interactions. This path can prevent the attainment of maximal community function, and trap communities in landscape locations where community function has low heritability. Exploiting these observations, we devise a species spiking approach to shift the path to improve community function heritability and consequently selection efficacy. We show that our approach is applicable to communities with complex and unknown function landscapes.


Introduction
Multi-species microbial communities often display community functions -biochemical activities not achiev-24 able by any member species alone. For example, a community of Desulfovibrio vulgaris and Methanococcus 25 maripaludis, but not either species alone, converts lactate to methane in the absence of sulfate [1]. Com-26 munity function arises from "interactions" where one community member influences the physiology of other 27 community members. Interactions are typically complex and difficult to characterize, making it challeng-28 ing to design communities with desired functions [2,3]. In a different approach, one could mutagenize 29 individual community members, assemble them at various ratios, and screen the resultant communities for 30 high community function. However, this requires community members being culturable, and the number  In each cycle, Newborns (blue shading) mature into Adults (yellow shading) during which intracommunity selection (olive) favors fast growers. The highest-functioning Adults are then chosen to reproduce, and this inter-community selection (red) favors high function. Note that both parent (top half) and offspring (bottom half) communities have Newborn and Adult stages. Thus, the four rows from top to bottom represent parent Newborn, parent Adult, offspring Newborn, and offspring Adult, respectively. (B) Heritability of community function. Community function is heritable if parent function and average offspring function are positively correlated (top). Community function is highly heritable if variations in function are largely attributable to variations in the heritable determinants (bottom). Heritability can be quantified by the slope of the least squares regression between the parent variable and the average offspring variable. Note that while community function is quantified at the Adult stage, community function determinants can be quantified at the Newborn stage if newly-arising mutations have negligible effects on community function within a cycle.

Results
can quantify selection efficacy as the progress in the heritable determinant f P (0).  ) to Adult (open triangle) over one maturation cycle. All arrows are "attracted" to the species composition attractor (orange curve). When M's cost f P is high, M goes extinct (broken orange line), and at lower f P values, M coexists with H (solid orange line). (C) Community function landscape. Community function P (T ) is plotted in 3D as a function of Newborn community's fraction of M biomass (φ M (0)) and average cost among M (f P (0)). To ensure that community function only has these two determinants, we made the following simplifications: i) Newborn's total biomass and the initial abiotic environment are both fixed; ii) stochastic death rate is relatively low to ensure largely deterministic population dynamics; and iii) mutation rate (0.002/cell/generation) and community maturation time T (~6 population doublings) are relatively small so that newly-arising genotypes remain rare within a cycle and not impact community function ( Figure  S7 of [15]). Indeed, community functions observed in simulations are well predicted by those calculated from the two determinants φ M (0) and f P (0) ( Figure S1). Species composition attractor and community function landscape were calculated using Eqs. 1-5 with the Newborn total biomass=100 and maturation time T = 17. (D) Newborn average cost f P (0) is largely heritable (positive slope), but Newborn's fraction of M biomass φ M (0) is not heritable (zero slope). Each open circle is obtained by averaging over~60 offspring Newborns descending from the same parent community. Error bars indicate one standard deviation. Note that for f P (0), linear regression did not pass through the origin. Offspring costs tend to be less than parent costs because during community maturation, the average cost declines as cheaters take over.

Constrained by attractor, community selection fails to reach maximal function
We now consider the case where ancestral M pays a cost smaller than what is optimal for maximal community 146 function. While inter-community selection will favor cooperator M paying a higher cost, intra-community 147 selection will always favor cheater M paying a lower cost ( Figure 1A). If we randomly select Adults to 148 reproduce, or allow each Adult to reproduce one offspring, then community function rapidly declines to zero 149 as cheaters take over [15]. If we apply community selection, can we improve community function to the 150 maximal value? 151 We find that although community selection can improve community function, community function fails to  The average cost paid by M f P stayed above the optimal f * P ( Figure 3B), which is surprising because high 159 cost would be disfavored by both intra-community selection (which favors fast growth and low cost) and 160 inter-community selection (which favors f * P ).

161
To understand this sub-optimal selection outcome, we overlaid the community function landscape with the

166
This explains the sub-optimal selection outcome: Like a hiker who is restricted to a trail that does not 167 traverse the mountain top, community function can only climb to the highest value along the attractor 168 (magenta circle), which is lower than the global maximum (black star). Consequently, the corresponding 169 f P (0) and φ M (0) are higher than their respective optima. In other words, selection outcome is constrained 170 along the species composition attractor by ecological forces.

171
The local geometry of community function landscape predicts the efficacy of 172 inter-community selection 173 The local geometry of community function landscape predicts how much inter-community selection can  184 We can now examine the efficacy of H-M community selection within community function landscape ( Figure   185 5Bi and ii). We allowed Newborn species composition to stochastically fluctuate, as if pipetting samples 186 from an Adult to seed Newborns while keeping the total Newborn turbidity fixed. This experimentally-facile 187 community reproduction method resulted in a "cloud" of Newborns around the species composition attractor  Average heritable determinant before and after inter-community selection are marked by a red square and a red star, respectively. Thus, progress made by inter-community selection can be measured as the distance from red square to red star.
heritable determinant (short red arrow in Figure 5Bii; statistics in the red box plot of Figure 5Biv). This small progress was just enough to counter the decline due to intra-community selection (olive box in Figure   195 5iv), resulting in a net of zero progress (black box in Figure 5iv). Consequently, despite 1000 selection cycles, 196 average community function and average heritable determinant among chosen communities barely improved 197 ( Figure 6A).
Since the contour of landscape around the attractor was not conducive to effective selection, we improved 199 selection by shifting the attractor to a more favorable location. Specifically, when we replaced 30% of total 200 Newborn biomass with ancestral non-evolving H ("30%-H spiking", teal box in Figure 5A i), the attractor 201 was forced into a region where community function contour lines were largely perpendicular to the axis of 202 heritable determinant (Figure 5Aii). Community function became more heritable ( Figure 5C iii). Inter-203 community selection made a larger improvement in the heritable determinant (longer red arrow in Figure   204 5Biii compared to Bii; red box in Figure 5Bv), and total progress was positive (black box in Figure 5Bv).

205
Consequently, community selection was effective ( Figure 6B). Resource. If we reproduce the Adult via "pipetting", then Newborn's total biomass and species composition 226 will fluctuate stochastically, adding two nonheritable determinants. If we additionally consider community 227 function measurement noise (a normal random variable with mean 0 and standard deviation comparable 228 to ancestral community function), we will have yet another nonheritable determinant. Overall, community 229 function will have 9 determinants, 6 heritable and 3 nonheritable. 230 We performed community selection in the above complex scenario (schematic in Figure S4). We started

237
Every 100 cycles, we evaluated spiking strategy, and updated the strategy when appropriate. Specifically, 238 we quantified community function heritability for 5 candidate spiking strategies (no spiking, 30%-H spiking,

239
60%-H spiking, 30%-M spiking and 60%-M spiking) by regressing parent function with median offspring 240 function (similar to Figure 5C). The current spiking strategy was then updated if an alternative strategy 241 conferred significantly higher community function heritability ( Figure S4B, Methods).

242
With the no-spiking strategy, community selection moderately improved community function and her-243 itable determinants ( Figure 7A and B; Figure S5A). Selection outcome was significantly improved when Under the no-spiking strategy, community function heritability is low. iii: Under the 30%-H spiking strategy, community function heritability is high. Error bars indicate one standard deviation. Empirical determination of this relationship between offspring and parent function can be used in lieu of landscape visualization to guide spiking strategy. Figure 6: Evolution dynamics of community function P (T ) and its heritable determinant f P (0) averaged over the chosen communities. Species composition but not total biomass of Newborns was allowed to stochastically fluctuate during community reproduction, as if pipetting an inoculum of a fixed turbidity from an Adult to seed its Newborns. In all cases, community selection is more effective under the 30%-H spiking strategy (where 30% of Newborn biomass was replaced by the non-evolving H) compared to the no-spiking strategy. In "top-2" strategy, top 2 communities each reproduced 50 Newborns. In "top-10" strategy, top 10 communities each reproduced 10 Newborns. In "with measurement noise", measured community function was the sum of the real community function and a normal random variable with a mean of 0 and a standard deviation of 100. Dashed lines correspond to values optimal for community function.  D; Figure S5B). In contrast, adopting the spiking strategy with the lowest heritability led to much worse 246 selection outcome ( Figure S5C), while randomly choosing spiking strategy led to variable results. It is also 247 noteworthy that during dynamic spiking, the adopted spiking strategy was not static ( Figure 7E). A static 248 60%-H or 30%-H spiking strategy offered negligible improvement over no spiking ( Figure S6). Therefore, it 249 is important to evaluate heritability periodically and update accordingly.

250
When we varied parameters of dynamic spiking, selection efficacy did not vary significantly. For example, 251 when the number of colonies used to make the spiking mix was changed from 5 to 1, 2, or 10 clones,  Intriguingly, community trait evolved to be stable even though stability was not directly selected for. This 263 is consistent with the notion that only heritable community trait can respond to selection.

264
A challenging aspect of optimizing community selection is that variation, selection strength, and heri- selection outcomes will be suboptimal ( Figure 3B). Furthermore, if the species composition attractor sits in 278 a region where community function heritability is low, then selection efficacy will be low.

279
These concepts have prompted us to devise species spiking strategies that might force the attractor to 280 a more favorable region with higher community function heritability. When we can visualize the landscape, 281 we can design spiking strategy to force the attractor to a favorable region (i.e., where contour lines of 282 community function are perpendicular to the axis of heritable determinant, Figures 4 and 5). When we can 283 not visualize the landscape, we can choose spiking strategy based on community function heritability ( Figure   284 7), a quantity that can be estimated in experiments ( Figure 5C). Spiking strategy needs to be dynamically 285 adjusted since evolving communities move in the landscape ( Figure S6). spiking, 30%-H spiking, and 30%-M spiking), selection efficacy can be higher than when we also included the 293 60%-H spiking and 60%-M spiking strategies ( Figure S10). This is because 60%-H spiking was so extreme that 294 species composition did not return to the attractor within one cycle ( Figure S11)  have the same f P = f P (0)), that death and birth processes are deterministic, and that there is no mutation. 320 P (T ) can then be numerically integrated from the following set of scaled differential equations for any given 321 pair of φ M (0) and f P (0) [15]: and at which g M max /2 is achieved when R is in excess.

335
To construct the landscape in Figure   where ∆s is a random variable with a distribution Newborn total biomass will be very close to the target BM target = 100, and species ratio will be very close per M biomass grown P mut mutation probability per cell 2 × 10 −3 division for each mutable phenotype

Community reproduction in the complex scenario when both H and M cells evolve 407
In the more complex scenario, both H and M evolve. We thus need to spike with evolved H and M clones.

408
Additionally, Newborns are spiked with H or M clones from their own lineage as demonstrated in Figure   409 S4A. Below we describe the the simulation code for the experimental procedure ( Figure S4A) we simulated.  . Bottom: Adults were mixed, and the mixture gave birth to 100 Newborns Figure S3: Spiking H at a wide range of percentages improves selection outcome. Adults ranking top 10% in measured community functions (= true community function + a normal random variable with zero mean and standard deviation of 100 -~10% of the ancestral community function) were chosen to reproduce. Each of the 10 chosen Adults reproduced 10 Newborns. The total biomass of a Newborn fluctuated around the target value (BM target = 100), of which the indicated percentage was "pipetted" from an H monoculture while the rest was "pipetted" from the parent Adult. Consequently, Newborn total biomass and species composition fluctuated stochastically according to sampling error. Top panels show the dynamics of heritable determinant f P (0), the cost paid by M averaged first within and then among chosen Adults. Bottom panels show P (T ), the average community function of chosen Adults. Black, cyan and gray curves represent three independent replicates. offspring Newborns under the respective spiking strategy, and when these mature, the median of 6 functions defines the "Offspring function" (brown axis label of the scatter plot). Cycle C: At the end of cycle C, we can scatter plot the function of each parent against the average function of its 6 offspring. The heritability is estimated from the slope of the least squares regression line. In parallel, the current spiking strategy continues (purple box) until Cycle C, where spiking strategy is updated. In this example, the updated current strategy is spiking strategy 3.      Figure S11: Spiking a species at too high a fraction can make a non-heritable determinant seem heritable. Compare 60%-H spiking (A) versus 30%-H spiking (B). (i) Species composition dynamics of two communities (blue and green). Starting at Newborns immediately after spiking (crosses; parent Newborns), species compositions moves toward the attractor (dotted line). In 30%-but not 60%-H spiking, species composition reaches the attractor by the end of maturation. After reproducing with respective spiking percentages, average offspring Newborn compositions are plotted with circles. In 30%-H spiking, average offspring Newborns have the same species composition (green circle overlapping with blue circle). In contrast, in 60%-H spiking, average offspring Newborns have different species composition (green circle above blue circle). (ii) In 60-but not 30%-H spiking, the determinant φ M (0) displays heritability. (iii) Consequently, although 60%-H spiking strategy confers similar community function heritability as 30%-H spiking strategy (similar slopes of the red dashed lines), community function improved faster under 30%-H spiking strategy ( Figure S10B) than under 60%-H spiking strategy ( Figure S10A). This is because in 60%-H spiking, a portion of community function heritability originates from the misleading heritability created by a nonheritable determinant. a microbial mutualism. Proceedings of the National Academy of Sciences, 107(5):2124-2129, 2010.