Artificial selection of microbial communities can become effective after using evolution-informed strategies

5 Multi-species microbial communities often display “community functions” stemming from interactions of 6 member species. Interactions are often difficult to decipher, making it challenging to design communi7 ties with desired functions. Alternatively, similar to artificial selection for individuals in agriculture and 8 industry, one could repeatedly choose communities with the highest functions to reproduce by randomly 9 partitioning each into multiple “Newborn” communities for the next cycle. However, community selec10 tion is challenging since rapid changes in species and genotype compositions can limit the heritability 11 of community function. To understand how to enact community selection, we used an individual-based 12 model to simulate this process to improve a community function that requires two species and is costly 13 to one species. Improvement was stalled by non-heritable variations in community function, such as the 14 stochastic populating of Newborn communities or measurement errors of community function. Com15 munity function improved when these non-heritable variations were suppressed in experimentally feasible 16 manners. 17


Introduction
In our simulations, cycles of selection were performed on a total of n tot = 100 communities. At the beginning of the first cycle, each Newborn had a total biomass of BM target =100 (60 M and 40 H each of biomass 1). In subsequent cycles, species ratio would converge to the steady state value (Figure2 bottom). Waste (not drawn) was in excess. The amount of Resource in each Newborn (not drawn) was fixed at a value that could support a total biomass of 10 4 . The maturation time T was chosen so that for an average community, Resource was not depleted (in experimental terms, this would avoid complications of the stationary phase). During maturation, Resource R, Byproduct B, Product P , and each cell's biomass were calculated from differential equations (Methods, Section 6). Death occurred stochastically to individual cells. A cell divided into two identical daughter cells once its biomass had reached a threshold of 2. After division, mutations (different shades of oval and rod) occurred stochastically to change a cell's phenotypes (maximal growth rate, affinity for metabolites, and M's f P ). At the end of a cycle (time T ), the top-functioning Adult with the highest Product P (T ) was chosen and diluted into as many Newborns as possible so that on average, each Newborn had a total biomass of approximately BM target . We then proceeded to the next top-functioning Adult until n tot = 100 Newborns were generated for the next selection cycle. Communities with red outlines exemplify one lineage.
6 . CC-BY-NC 4.0 International license not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint (which was this version posted September 21, 2018. . https://doi.org/10.1101/264689 doi: bioRxiv preprint always supplied in excess and thus did not enter our equations. Note that except for the first cycle, the  Figure S5). M cell's actual biomass growth rate was(1 − f P ) fraction of M's potential growth rate (Eq. . When H and M's growth parameters were allowed to mutate, community function was lower than if growth parameters were fixed to the ancestral values (the left bar being lower than the right bar). In these simulations, since H evolved to grow very fast, we adjusted initial Resource R(0) to support 10 5 total biomass (higher than the standard 10 4 total biomass). (B) The evolutionary upper bound of the maximal growth rate of M exceeds that of H (g * Hmax = 0.3 < g * M max = 0.7, Table 1). When H and M's growth parameters were allowed to mutate, community function increased to a higher level compared to when growth parameters were fixed to the ancestral values (the left bar higher than the right bar). In these simulations, R(0) supported standard amount (10 4 ) of total biomass. In both (A) and (B), natural selection improved growth parameters when growth parameters were allowed to mutate. Here the Newborn community has 54 M and 46 H cells of biomass 1, which corresponds to the species composition optimal for P (T ) under these conditions. The corresponding maximal P * (T ) could not be further improved if we allowed all growth parameters and f P to mutate ( Figure S20). Thus, P * (T ) is locally maximal in the sense that small deviation will always reduce P (T ).
by fixing H and M's growth parameters to their upper bounds and only allowing f P to mutate. This 229 simplification is justified for three reasons. First, during these community selection simulations, growth 230 parameters important to community function improved to their upper bounds ( Figure S12C). Second, 231 suppose that a mutation changes a growth parameter already at its upper bound. Since the growth 232 parameter is already at its upper bound, the mutation can only reduce it, and thus, the mutant is disfa-233 vored by natural selection. And since we are studying cases where improving growth parameters improves 234 community function, the mutant will also reduce community function and hence its host community is 235 less likely to reproduce. Overall, the mutant will not persist which means that the growth parameter will 236 remain at its upper bound. Third, our conclusions hold regardless of whether we fix growth parameters 237 or not ( Figure S15). 238 After fixing growth parameters, we can now focus on how f P values of M cells evolve during commu- 239 nity selection. f P is of particular importance because engineered microbes pay a fitness cost to synthesize 240 a product. An M cell with a higher f P will grow slower than a low-producer, and thus be selected against 241 by natural selection during community maturation. However, an M cell with higher f P will result in higher 242 community function, and thus its host community has a higher chance of being chosen to reproduce. As 243 we will demonstrate, proper design of community selection regimen is critical to counter natural selection 244 and successfully improve costly community function. Our conclusions hold even for the more difficult The copyright holder for this preprint (which was this version posted September 21, 2018. . https://doi.org/10.1101/264689 doi: bioRxiv preprint Ineffective community selection is due to non-heritable variations 252 in community function 253 We simulated community selection after fixing all growth parameters to evolutionary upper bounds 254 and allowing only f P to be modified by mutations. As described in "Model structure," we simulated 255 community maturation by tracking the biomass growth of each H and M cells, cell division and death, 256 mutation in f P of each M cell, and metabolite release and consumption. We simulated community 257 reproduction by incorporating stochastic fluctuations associated with volumetric dilution of the Adult 258 community into Newborn communities using a pipette. Both f P and community function P (T ) barely 259 improved over thousands of selection cycles, even though both were far from their theoretical maxima 260 ( Figure 6A and B). Note that community function was above the ancestral value (( Figure 6B, brown 261 star) because growth parameters were fixed to evolutionary upper bounds.

262
To understand why community selection failed to improve community function, we examine the  Newborn f P (0) (Eq 6-10). Note that because the composition of a community varies as it goes from its 271 Newborn stage to its Adult stage, these determinants are all defined at a community's Newborn stage.

272
A determinant is considered heritable if the determinant of a parent community (e.g. the red tube from 273 the top row of Figure 3) and the determinants of its offspring communities (red tubes from the bottom 274 row of Figure 3) are correlated. Among the three determinants, f P (0) can be considered "heritable" as 275 shown in Figure S29: if a parent community has a high f P (0), i.e. a high average f P at its Newborn 276 stage, it will have a high average f P at its Adult stage. Since its offspring communities inherit its M 277 cells, they will have high f P (0). On the other hand, Newborn total biomass BM (0) is not heritable 278 since it is not correlated between a parent community and its offspring communities, as shown in Figure   279 S29. This is because when an Adult community reproduced, the dilution factor was adjusted so that the 280 total biomass of an offspring Newborn community was on average the constant target biomass BM target .

281
The fraction of M biomass of a Newborn φ M (0) is not heritable either, as shown in Figure S29. This  In successful community selection, higher community function should correlate with higher average 289 Newborn f P ( f P (0) ), because the latter is the heritable determinant of the former. However, we 290 observed little correlation between community function P (T ) and f P (0), but strong correlation between 291 community function and its non-heritable determinants (Figure 7). For example, the Newborn that 292 would achieve the highest function (left magenta dot) had a below-median f P (0), but had high total . Magenta dashed lines: f * P optimal for P (T ) and maximal P * (T ) when all five growth parameters are fixed at their upper bounds and φ M (0) is optimal for community function. Black, cyan and gray curves are three independent simulation trials. P (T ) is averaged across the two selected Adults and has the unit of r P , and f P is obtained by averaging within each selected Adult and then averaging across the two selected Adults. Figure 7: When community selection is ineffective, community function correlates weakly with its heritable determinant and strongly with non-heritable determinants. From community selection simulation, we randomly chose a selection cycle where 100 Newborns matured into 100 Adults. We plotted community function P (T ) of each Adult against characteristics of its corresponding Newborn community that determine community function. Each dot represents one community, and the two magenta dots indicate the two "successful" Newborns that achieved the highest community function at adulthood. (A) P (T ) only weakly correlates with f P (0) (f P averaged over all M individuals in a Newborn). (B-C) P (T ) strongly correlates with Newborn total biomass BM (0) and Newborn fraction of M biomass φ M (0).
higher-than-average M biomass would more thoroughly convert Resource (and Byproduct) to Product.

298
Similarly, if a Newborn started with higher-than-average fraction of H biomass (dotted lines in bottom 299 panels of Figure S24), then H would produce higher-than-average Byproduct which meant that M would 300 endure a shorter growth lag, and make more Product.

301
In summary, variation in community function was dominated by variations in non-heritable determi- Reducing non-heritable variations in community function should enable community selection to work.

308
One possibility is to reduce the stochastic fluctuations in non-heritable determinants BM (0) and φ M (0). Section 6). similar P (T ) since "unlucky" communities would have time to "catch up" as "lucky" communities wait in 322 stationary phase. Indeed, with extended T , community function improved without having to fix BM (0)  In summary, non-heritable variations in community function must be sufficiently suppressed for com-335 munity selection to work. During community selection, seemingly innocuous experimental procedures 336 such as pipetting could be problematic, and a more precise procedure such as cell sorting might be 337 required. Our conclusions held when we used a different mutation rate (2 × 10 −5 instead of 2 × 10 −3 338 mutation per cell per generation per phenotype, Figure S26), a different distribution of mutation effects 339 (a non-null mutation increased or decreased f P by on average 2%, Figure S27), or incorporating epistasis 340 (a non-null mutation would likely reduce f P if the current f P was high, and enhance f P if the current f P 341 was low; Figure S28; Figure S9; Methods Section 5). Our conclusions also hold when improved growth   The copyright holder for this preprint (which was this version posted September 21, 2018. . https://doi.org/10.1101/264689 doi: bioRxiv preprint in simulations of M-group selection, f P gradually increased to f * P, M ono = 0.13 ( Figure S22). However, 365 f P optimal for monoculture P (T ) is much lower than f P optimal for community P (T ) ( Figure 5; see 366 Methods Section 8 for an explanation).

367
Artificial selection of whole communities to improve a costly community function requires careful use one culture to inoculate one Newborn, then this chance reduces below 0.01%. However, these two 377 arrangements resulted in almost identical evolutionary dynamics of average f P and P (T ) ( Figure S32).

378
In certain regards, community selection is similar to selection of mono-species groups. Group Figure S1).

385
Community selection and group selection differ in two aspects. First, species interactions in a com-386 munity could drive species composition to a value sub-optimal for community function ([60]). This interfere with community selection ( Figure 6). In the extreme case, a member species may even be lost 393 by chance. Even if a fixed biomass of each species is sorted into Newborns, heredity is much reduced in 394 community selection due to random sampling of genotypes from member species 3 .

395
Although suppressing non-heritable variations in a trait will always increase selection efficacy, here Grow rate g H depends on the level of ResourceR (hat^representing pre-scaled value) as described by Total potential growth rate of M g M depends on the levels of Resource and Byproduct (R andB) 436 according to the Mankad-Bungay model [25] due to its experimental support ( Figure S5): Figure S4). 1 − f P fraction of M growth is channeled to 438 biomass increase. f P fraction of M growth is channeled to making Product: where r P is the amount of Product made at the cost of one M biomass (tilde~representing scaling 440 factor, see below and Table 1).

441
ResourceR is consumed proportionally to the growth of M and H; ByproductB is released propor-442 tionally to H growth and consumed proportionally to M growth: Here,ĉ RM andĉ RH are the amounts ofR consumed per potential M biomass and H biomass, respectively.
CC-BY-NC 4.0 International license not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint (which was this version posted September 21, 2018. . https://doi.org/10.1101/264689 doi: bioRxiv preprint where R =R/ R(0) and K HR =K HR / R(0). Thus, the unscaled g H (R) and the scaled g H (R) share 458 identical forms. After scaling, the value of R(0) becomes irrelevant (1 with no unit). Similarly, since . 460 Thus, scaled equations are We have not scaled time here, although time can also be scaled by, for example, the community  Table 1, and 464 variables in our model and simulations are summarized in Table 2.

465
For this H-M community, the species ratio at time T at can be estimated in the following manner.

466
From Eq. 10: If we approximate Eq. 6-7 by ignoring the death rates so that dH Eq. 11 becomes If B is the limiting factor for the growth of M so that B is mostly depleted, we can approximate If Newborn inherits parent Adult's φ M , then the steady state φ M,SS is In our simulations, because we supplied the H-M community with abundant R to avoid stationary 475 phase, H grows almost at the maximal rate through T and releases B. If f P is not too large (f P < 0.4), 476 which is satisfied in our simulations, M grows at a maximal rate allowed by B and keeps B at a low 477 level. Thus, Eq. 14 is applicable and predicts the steady-state φ M,SS well (see Figure S7). Note that 478 significant deviation occurs when f P > 0.4. This is because when f P is large, M's biomass does not 479 grow fast enough to deplete B so that we cannot approximate B(T ) ≈ 0 anymore. . This way, we could start with communities of H and M whose growth 486 parameters are already maximized ("growth-adapted"), since mutations that reduce growth parameters 487 will be selected against by both natural selection and community selection. In other words, only f P can 488 mutate, and higher f P will be favored by community selection but disfavored by natural selection. With 489 only one mutable parameter (f P ), we can identify the optimal f P * associated with maximal community 490 function ( Figure 5). Evolutionary modeling is also greatly simplified.

491
For ancestral H, we set g Hmax = 0.25 (equivalent to 2.8-hr doubling time if we choose hr as the time 492 unit), K HR = 1 and c RH = 10 −4 (both with unit of R(0)) ( Table 1) Table 1.

498
To ensure the coexistence of H and M, M must grow faster than H for part of the maturation cycle.

499
Since we have assumed M and H to have the same affinity for R ( initial amount of Resource in Newborn scaling factor, 1 f P fraction of M growth diverted to producing P 0.10 0.13# no change per M biomass grown, scaled against r B P mut mutation probability per cell 2 × 10 −5~2 × 10 −3 division for each mutable phenotype Table 1: Parameters for ancestral and evolved (growth-and mono-adapted) H and M. Parameters in the "Ancestral" column are used for Figure S12, while those in the "Evolved" column are used for all other figures. For maximal growth rates, * represents evolutionary upper bound. For K SpeciesM etabolite , * represents evolutionary lower bound, which corresponds to evolutionary upper bound for Species's affinity for Metabolite (1/K SpeciesM etabolite ). # is from Figure 5B. In Methods Section 2, we explain our parameter choices (including why we hold some parameters constant during evolution).

Symbols
Definition The biomass of M or H in a community at time t The total biomass in a community at time t φ M (t) The fraction of M biomass at time t BM target Pre-set target total biomass of Newborns during community reproduction The integer number of M or H cells in a community at time t ϕ M (t) The fraction of M individuals at time t The biomass (length) of an individual M or H cell at time t, ranged between 1 and 2 P (t) The amount of Product P in a community at time t, scaled by r P R(t) The amount of Resource R in a community at time t, scaled by R(0) The amount of Byproduct B in a community at time t, scaled by r B n dil The fold dilution when reproducing an Adult community n tot Total number of communities under selection T Community maturation time, corresponding to the duration of a selection cycle   Figure 5A; magenta dashed line in Figure 6).

551
Second, growth-adapted H and M are evolutionarily stable in the sense that deviations (reductions) 552 from upper bounds will reduce both individual fitness and community function, and are therefore disfa-553 vored by natural selection and community selection.

554
Below, we present evidence that within our parameter ranges (Table 1), improving growth parameters 555 generally improves community function. When f P is optimal for community function (f * P = 0.41), if we 556 fix four of the five growth parameters to their upper bounds, then as the remaining growth parameter 557 improves, community function increases (magenta lines in top panels of Figure S16). Moreover, mutants 558 with a reduced growth parameter are out-competed by their growth-adapted counterparts (magenta lines 559 in bottom panels of Figure S16).

560
When f P = f * P, M ono = 0.13 (optimal for M-monoculture function in Figure 5B; the starting genotype 561 for most community selection trials in this paper), community function and individual fitness generally 562 increase as growth parameters improve (black dashed lines in Figure S16). However, when M's affinity B M (corresponding to limiting R and abundant B),  cycle (Eq. 16), and therefore M cells with lower affinity for R will grow faster than those with higher 574 affinity ( Figure S17). At the end of each cycle, the opposite is true ( Figure S17). As f P decreases, M 575 has higher growth capacity, and thus the first stage of B limitation lasts longer. Consequently, M can 576 gain higher overall fitness by lowering affinity for R ( Figure S17A). The copyright holder for this preprint (which was this version posted September 21, 2018. . https://doi.org/10.1101/264689 doi: bioRxiv preprint parameters to vary ( Figure S18), then 1/K M R decreases somewhat, yet the dynamics of f P is similar to when we only allow f P to change (compare Figure S18D with Figure 6A). Indeed, allowing both f P and 583 1/K M R to evolve does not change our conclusions as shown in Figure S19.  Our "mutation rate" refers to the rate of mutations that either enhance a phenotype ("enhancing  Figure S8 for interpreting PDF), in the following manner.
where µ w (j) is short-hand notation for µ w ( s W T = j × 0.04) and so on. Our calculated µ ∆s (∆s) with 647 error bar of δµ ∆s is shown in Figure S8. 648 Our reanalysis demonstrates that distributions of mutation fitness effects µ ∆s (∆s) are largely con-649 served regardless of nutrient conditions and mutation types ( Figure S8B). In all cases, the relative fitness 650 changes caused by beneficial (fitness-enhancing) and deleterious (fitness-diminishing) mutations can be 651 approximated by separate exponential distributions with different means s + and s − , respectively. After 652 normalization to have a total probability of 1, we have: We fit the Dunham lab haploid data (since microbes are often haploid) to Eq. 19, using µ ∆s (i)/δµ ∆s (i) 654 as the weight for non-linear least squared regression (green lines in Figure S8B). We obtain s + =   After division, each cell has a probability of P mut to acquire a mutation that changes each of its 767 mutable phenotype (Methods, Section 4). As an example, let's consider mutations in f P . If a mutation 768 occurs, then f P will be multiplied by (1 + ∆f P ), where ∆f P is determined as below.

769
First, a uniform random number u 1 between 0 and 1 is generated. If u 1 ≤ 0.5, ∆f P = −1, which 770 represents 50% chance of a null mutation (f P = 0). If 0.5 < u 1 ≤ 1, ∆f P follows the distribution The copyright holder for this preprint (which was this version posted September 21, 2018. . https://doi.org/10.1101/264689 doi: bioRxiv preprint is generated via inverse transform sampling. Specifically, C(∆f P ), the cumulative distribution function 774 (CDF) of ∆f P , can be found by integrating Eq. 19 from -1 to ∆f P : The two parts of Eq. 25 overlap at C(∆f In order to generate ∆f P satisfying the distribution in Eq. 19, a uniform random number u 2 between 777 0 and 1 is generated and we set C(∆f P ) = u 2 . Inverting Eq. 25 yields When epistasis is considered, s + (f P ) = s +init /(1 + g × (f P /f P, init − 1)) and s − (f P ) = s −init ×

781
(1 + g × (f P /f P, init − 1)) are used in Eq. 26 to calculated ∆f P for each cell with different current f P 782 (Methods Section 5).

783
If a mutation increases or decreases the phenotypic parameter beyond its bound (Table 1), the 784 phenotypic parameter is set to the bound value.

785
The above growth/death/division/mutation cycle is repeated from time 0 to T . Note that since the M (0) (true if T is long enough for cells to double at least three or four times).

847
If we define community function as P (T )/M (T ) ≈ f P 1−f P (total Product normalized against M 848 biomass in Adult community), then higher f P 1−f P or higher f P always leads to higher community function.

849
Higher f P in turn leads to M extinction (Figure 1).

850
If the community function is instead defined as P (T )/M (0), then ( Figure S24B). Thus, we end up selecting communities with small φ M (0) ( Figure S6). This means that 854 Manufactures could get lost during community reproduction, and community selection then fails.

855
If Resource is unlimited, then it will be problematic to reproduce an Adult by diluting it by a fixed-856 fold to Newborns. This is because with unlimited Resource, there is no competition between H and M. For groups or communities with a certain T g M dt, we can calculate f P optimal for community function 862 from Eq. 27 by setting We have where "Cov" means covariance and "Var" means variance, and φ M (T ) is the fraction of M biomass in 873 the Adult community from which Newborns are generated.
The larger B 0 , the less inhibitory effect Byproduct has on H and when B 0 → +∞ Byproduct does not 880 inhibit the growth of H. For simulations in Figure S30, Figure S1: Artificial selection is more challenging for multi-species communities than for individuals or mono-species groups. Artificial selection can be applied to any population of entities [95]. An entity can be an individual (A), a mono-species group (B), or a multi-species community (C). Unlike natural selection which selects for fastest-growing cells, artificial selection generally selects for traits that are costly to individuals. In each selection cycle, a population of "Newborn" entities grow for maturation time T to become "Adults". Adults expressing a higher level of the trait of interest (darker shade) are selected to reproduce. An individual reproduces by making copies of itself, while an Adult group or community can reproduce by randomly splitting into multiple Newborns of the next selection cycle. Successful artificial selection requires that i) entities display trait variations; ii) trait variations can be selected to result in differential entity survival and reproduction; and iii) entity trait is sufficiently heritable from one selection cycle to the next [96]. In all three types of selection, entity variations can be introduced by mutations and recombinations in individuals. However, heredity can be low in community selection. (A) Artificial selection of individuals has been successful [97, 98,99], since a trait is largely heritable so long as mutation and recombination are sufficiently rare. (B, C) In group and community selection, if T is small so that newly-arising genotypes cannot rise to high frequencies within a selection cycle, then Adult trait is mostly determined by Newborn composition (the biomass of each genotype in each member species). Then, variation can be defined as the dissimilarity in Newborn composition within a selection cycle, and heredity as the similarity of Newborn composition from one cycle to the next for Newborns connected through lineage (tubes with red outlines in Figure 3). (B) Artificial selection of mono-species groups has been successful [43,45,13]. Suppose cooperators but not cheaters pay a fitness cost to generate a product (shade). Artificial selection for groups producing high total product favors cooperator-dominated groups, although within a group, cheaters grow faster than cooperators. At a large Newborn population size (top), all Newborns will harbor similar fractions of cheaters, and thus inter-group variation will be small. During maturation, cheater frequency will increase, thereby diminishing heredity. In contrast, when Newborn groups are initiated at a small size such as one individual (bottom), a Newborn group will comprise either a cooperator or a cheater, thereby ensuring variation. Furthermore, even if cheaters were to arise during maturation, a fraction of Newborns of the next cycle will by chance inherit a cooperator, thereby ensuring some level of heredity. Thus, group selection can work when Newborn size is small. (C) Artificial selection of multi-species communities may be hindered by insufficient heredity. During maturation, the relative abundance of genotypes and species can rapidly change due to ecological interactions and evolution, which compromises heredity. During community reproduction, stochastic fluctuations in Newborn composition further reduce heredity. Figure S2: Large Newborn size or long maturation time allows non-producers to accumulate and reduces inter-community variation. From target total biomass of 10 2 or 10 4 wild-type cells, M population expands and mutates during a maturation time T that lasts 100 or 6 generations. Immediately following cell division, wild-type daughter cells mutate to non-producers with a probability of 10 −3 . Wild-type and mutant cells follow exponential growth. The growth rate of wild-type cells is 0.87 that of mutants. The fraction of biomass made up by mutants at each wild-type doubling is shown. Note different scales. At 10 2 total biomass, a small fraction mutation (e.g. 0.005) means that some communities will remain free of non-producers at the end of T . Figure S3: A comparison of different growth models. We model exponential biomass growth in excess metabolites. Thick black line: analytical solution with biomass growth rate (0.7/time unit). Grey dashed line: simulation assuming that biomass increases exponentially at 0.7/time unit and that cell division occurs upon reaching a biomass threshold, an assumption used in our model. Color dotted lines: simulations assuming that cell birth occurs at a probability equal to the birth rate multiplied with the length of simulation time step (∆τ = 0.05 time unit). When a cell birth occurs, biomass increases discretely by 1, resulting in step-wise increase in color dotted lines at early time.

33
. CC-BY-NC 4.0 International license not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint (which was this version posted September 21, 2018. . https://doi.org/10.1101/264689 doi: bioRxiv preprint  Figure S5: A comparison of dual-substrate models. Suppose that cell growth rate depends on each of the two substrates S 1 and S 2 in a Monod-like, saturable fashion. When S 2 is in excess, the S 1 at which half maximal growth rate is achieved is K 1 . When S 1 is in excess, the S 2 at which half maximal growth rate is achieved is K 2 . (A) In the "Double Monod" model, growth rate depends on the two limiting substrates in a multiplicative fashion. In the model proposed by Mankad and Bungay (B), growth rate takes a different form. In both models, when one substrate is in excess, growth rate depends on the other substrate in a Monod-fashion. However, when S 1 K 1 = S 2 K 2 = 1, the growth rate is predicted to be   Figure S7: Comparison between the steady-state φ M,SS calculated from Eqs. 6-10 (black curve) and from Eq. 14 (red line).

35
. CC-BY-NC 4.0 International license not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint (which was this version posted September 21, 2018. . https://doi.org/10.1101/264689 doi: bioRxiv preprint Figure S8: Probability density functions of changes in relative fitness due to mutations ( µ ∆s (∆s) ). We derived µ ∆s (∆s) from the Dunham lab data [34] where bar-coded mutant strains were competed under sulfate-limitation (red), carbon-limitation (blue), or phosphate-limitation (black). Error bars represent uncertainty δµ ∆s (the lower error bar is omitted if the lower estimate is negative). In the leftmost panel, green lines show non-linear least squared fitting of data to Eq. 19 using all three sets of data. Note that data with larger uncertainty are given less weight, and thus deviate more from the fitting lines. For an exponentially-distributed probability density function p(x) = exp(−x/r)/r where x, r > 0, the average of x is r. When plotted on a semi-log scale, we get a straight line with slope 1/r, and inverting this gets us the average effect r. From the green line on the right side, we obtain the average effect of enhancing mutations s + = 0.050 ± 0.002, and from the green line on the left side, we obtain the average effect of diminishing mutations s − = 0.067 ± 0.003. The probability of a mutation altering a phenotype by ±α is the shaded area.

36
. CC-BY-NC 4.0 International license not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint (which was this version posted September 21, 2018. . If current f P is low (left), enhancing mutations are more likely to occur (the area to the right of f P = 0 becomes bigger) and their mean mutational effect becomes larger (mean=1/slope becomes larger due to smaller slope), while diminishing mutations are less likely to occur and their mean mutational effect is smaller. If current f P is high (right), the opposite is true.    Figure S14: Reducing non-heritable variations improves community function even when improved growth parameters impair community function. Similar to Figure S13, the upper bound for g Hmax (g * Hmax = 0.8) is larger than that of g M max (g * M max = 0.7). When both BM (0) and φ M (0) were allowed to fluctuate stochastically, community function declined to very low levels due to low abundance of M ( Figure S13F). Note that M did not go extinct because communities without any M would not be chosen to reproduce. When both BM (0) and φ M (0) were fixed, both f P and P (T ) improved over cycles. Here, Resource supplied to Newborn communities could support 10 5 total biomass to accommodate faster growth rate.

41
. CC-BY-NC 4.0 International license not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint (which was this version posted September 21, 2018. . https://doi.org/10.1101/264689 doi: bioRxiv preprint Figure S15: Community selection succeeds when controlling the right experimental variables even if growth parameters are allowed to be modified by mutations. Dynamics of (A) f P (T ) and (B) P (T ) of selected communities when the maturation time T = 17, g * Hmax = 0.3 and g * M max = 0.7. All other simulation parameters are in Table 1. Compared to simulations whose results are presented in    Figure 5B), and perform community selection while allowing all growth parameters and f P to vary. M's affinity for R 1/K M R decreases slightly because at low f P = 0.13, M with a lower affinity for R (lower 1/K M R ) slightly improves individual fitness while slightly decreasing community function ( Figure  S17). Other growth parameters (ḡ M max ,ḡ Hmax , 1/K M B and 1/K HR ) remain mostly constant during community selection because mutants with lower-than-maximal values are selected against by natural selection and by community selection ( Figure S16). Other legend details can be found in Figure S10.

44
. CC-BY-NC 4.0 International license not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint (which was this version posted September 21, 2018. . https://doi.org/10.1101/264689 doi: bioRxiv preprint Figure S19: Evolution dynamics of selected Adult communities when both f P and K M R are allowed to mutate. The dynamics are similar to when only f P is allowed to vary ( Figure 6). Other legend details can be found in Figure S10.

45
. CC-BY-NC 4.0 International license not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint (which was this version posted September 21, 2018. . https://doi.org/10.1101/264689 doi: bioRxiv preprint Figure S20: Local optimality of community function P * (T ). We start each Newborn community with total biomass BM (0)=100, all five growth parameters at their upper bounds, and f * P = 0.41 and φ * M (0) = 0.54 to achieve P * (T ). We then allow all five growth parameters and f P to mutate while applying community selection. To ensure effective community selection (Figure 6), BM (0) is fixed to 100, and φ M (0) is fixed to φ M (T ) of the previous cycle during community reproduction. We find that all five growth parameters remain at their respective evolutionary upper bounds. At the end of the first cycle (Cycle = 1 in insets), even though f P has not changed, P (T ) has already declined from the original magenta dashed line. This is because species interactions have driven φ M (0) from the optimal φ * M (0) (=0.54) to near the steady state value (φ M =0.73, compare with φ M,SS represented by the green dashed line in Figure 1C bottom panel). Later, over hundreds of cycles, f P gradually increases, which increases P (T ). However, P (T ) is still below maximal. This is because species composition gravitates toward steady state φ M,SS which deviates from the optimal φ * M (0) ([60]). Other legend details can be found in Figure S10.

47
. CC-BY-NC 4.0 International license not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint (which was this version posted September 21, 2018. . https://doi.org/10.1101/264689 doi: bioRxiv preprint Figure S24: Variations in community function can arise from non-heritable variations in Newborn compositions. An average Newborn community (solid lines) has a total biomass of 100 with 75% M. (A) A "lucky" Newborn community (dotted lines), by stochastic fluctuations, has a total biomass of 130 with 75% M. Even though the two communities share identical f P = 0.1, the Newborn with 130 total biomass has its M growing to a larger size (left), depleting more Resource (middle), and making more Product (right) by the end of short T (=17). (B) A "lucky" Newborn community (dotted lines), by stochastic fluctuations, has 100 total biomass with 65% M. Even though the two communities share identical f P = 0.1, the Newborn with lower φ M (0) (dotted) has its M enjoying a shorter growth lag and growing to a larger size (left), depleting more Resource (middle), and making more Product (right) by the end of short T (=17). In both cases, the difference between lucky (dotted) and average (solid) communities is diminished at longer T (T = 20) compared to shorter T (T = 17, dash dot line).

48
. CC-BY-NC 4.0 International license not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint (which was this version posted September 21, 2018. . https://doi.org/10.1101/264689 doi: bioRxiv preprint The copyright holder for this preprint (which was this version posted September 21, 2018. . https://doi.org/10.1101/264689 doi: bioRxiv preprint Figure S28: Evolutionary dynamics of selected Adults when epistasis is considered. When we incorporate different epistasis strengths (epistasis factor of 0.3 and 0.8), we obtain essentially the same conclusions as when epistasis is not considered (Figure 6). Other legend details can be found in Figure  6. Figure S29: Correlation of the three determinants of the community function between parent communities and offspring communities. The scatter plots show the correlation between the offspring communities' determinants and their parent community's determinants. For example, the abscissa of each point in (A) indicates f P (0) of a parent community; the ordinate and error bar of each point in (A) indicate the mean and standard deviation of f P (0) among the offspring communities formed out of the parent community. 100 communities from the 100th cycle of one of the simulations shown in Figure  6A and B are analyzed to generate this plot.

52
. CC-BY-NC 4.0 International license not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint (which was this version posted September 21, 2018. . https://doi.org/10.1101/264689 doi: bioRxiv preprint The copyright holder for this preprint (which was this version posted September 21, 2018. . https://doi.org/10.1101/264689 doi: bioRxiv preprint Figure S31: Selection dynamics in the presence of measurement uncertainty in P (T ). Evolution of f P (T ) and P (T ) when Adult communities are chosen to reproduce based on "measured P (T )" -the sum of actual P (T ) and an "uncertainty term" randomly drawn from a normal distribution with zero mean. The amplitude of the noise is characterized by the standard deviation of the normal distribution. In the left, center, and right panels, the noise terms were drawn from normal distributions with standard deviations of 5%, 7.5%, and 10% of the ancestral P (T ), respectively. The middle and lower panels show the average actual P (T ) and the average measured P (T ), respectively.

54
. CC-BY-NC 4.0 International license not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint (which was this version posted September 21, 2018. . https://doi.org/10.1101/264689 doi: bioRxiv preprint Figure S32: Different ways of inoculating the Newborn communities of the first cycle had limited impact on selection dynamics. (Top Panel) The number of communities initially free of non-producer M mutants depends on whether each Newborn community from the first cycle is inoculated from a distinct M monoculture. Each Newborn community for the first selection cycle was then inoculated with 60 M cells, either from the same M monoculture (Left panel), or from distinct M monocultures (Right panel). This pre-growth process is repeated 100 times, and the frequency of total numbers of Newborn communities out of 100 without non-producers is plotted. Selection dynamics are almost the same when the Newborn communities from the first cycle are inoculated by (Left panel) the same M monoculture or by (Right panel) distinct monocultures. Here we assumed that each monoculture grew from a single non-null M cell. This M cell went through~23 doublings and therefore multiplied into~10 7 cells. Every time a non-null M cell divides the mother and daughter cells can independently mutate and become a null M cell with f P = 0 at a fixed probability of 10 −3 . Assuming that all non-null M cells have identical f P = 0.13, non-null M cells grow at a rate 87% of that of a null cell. As a result, after~23 doublings, the M monocultures have on average~3% null mutants.

55
. CC-BY-NC 4.0 International license not certified by peer review) is the author/funder. It is made available under a The copyright holder for this preprint (which was this version posted September 21, 2018.