Simulations reveal insidious challenges to artificial community selection and possible strategies for success

Multi-species microbial communities often display “community functions” stemming from interactions of member species. Interactions are often difficult to decipher, making it challenging to design communities with desired functions. Alternatively, similar to artificial selection for individuals in agriculture and industry, one could repeatedly choose communities with the highest community functions to reproduce by randomly partitioning each into multiple “Newborn” communities for the next cycle. However, previous efforts in selecting complex communities have generated mixed outcomes that are difficult to interpret. To understand how to effectively enact community selection, we simulated community selection to improve a community function that requires two species and imposes a fitness cost on one or both species. Our simulations predict that improvement could be easily stalled unless various aspects of selection, including species choice, selection regimen parameters, and stochastic populating of Newborn communities, were carefully considered. When these considerations were addressed in experimentally feasible manners, community selection could overcome natural selection to improve community function, and in some cases, even force species to evolve to coexist. Our conclusions hold under various alternative model assumptions, and are thus applicable to a variety of communities.


Introduction 21
Multi-species microbial communities often display important community functions, defined as biochemical 22 activities not achievable by member species in isolation. For example, a six-species microbial community, 23 but not any member species alone, cleared relapsing Clostridium difficile infections in mice [1]. Com-24 munity functions arise from interactions where an individual alters the physiology of another individual. 25 Thus, to improve community functions, one could identify and modify interactions [2,3]. In reality, this 26 is no trivial task: each species can release tens or more compounds, many of which may influence the 27 partner species in diverse fashions [4,5,6,7]. From this myriad of interactions, one would then need 28 to identify those critical for community function, and modify them by altering species genotypes or the 29 abiotic environment. One could also artificially assemble different combinations of species or genotypes 30 Community selection can enforce species coexistence 330 In most communities, species coexistence may not be guaranteed due to competition for shared resources. 331 Here, we show that properly executed community selection could also improve the functions of such 332 communities, in part by enforcing species coexistence. Consider an H-M community where, unlike the 333 H-M community we have considered so far, H had the evolutionary potential to grow much faster than 334 M. In this case, high community function not only required M to pay a fitness cost of f P , but also 335 required H to grow sufficiently slowly to not out-compete M. 336 We started community selection at ancestral growth parameters, and allowed them and f P to mutate. 337 When community selection was ineffective ("top-dog" with "pipetting"; Figure 6A), H's maximal growth 338 rate evolved to exceed M's maximal growth rate ( Figure 6A, Figure S17 ). This drove M to almost 339 extinction, and community function was very low ( Figure 6A). During effective community selection 340 ("top-dog" with "cell-sorting" or "bet-hedging" with "pipetting", Figure   Robust conclusions under alternative model assumptions 345 We have demonstrated that when selecting for high H-M community function, seemingly innocuous 346 experimental procedures (e.g. choosing the top-functioning Adults and pipetting portions of them to 347 form Newborns) could be problematic. Instead, more precise procedures (e.g. "cell-sorting") or reduced 348 selection strength (e.g. "bed-hedging") might be required. Our conclusions hold when we used a much 349 lower mutation rate (2 × 10 −5 instead of 2 × 10 −3 mutation per cell per generation per phenotype, Figure   350 S18), although lower mutation rate slowed down community function improvement. Our conclusions 351 also hold when we used a different distribution of mutation effects (a non-null mutation increasing or 352 decreasing f P by on average 2% or eliminating null mutants; Figure S19), or incorporating epistasis (i.e. 353 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 354 was low; Figure S20; Figure S21; Methods Section 5). 355 To further test the generality of our conclusions, we simulated community selection on a mutualistic

Discussions
However, if we solely rely on species types, then without a constant influx of new species, community function will likely level off quickly [12,21]. Here, we consider artificial selection of communities 371 with defined member species so that improvement of community function requires new genotypes that 372 contribute more toward the community function of interest at a cost to itself. 373 Community selection can be challenging but is feasible 374 Artificial selection of whole communities to improve a costly community function requires careful consid-  385 We can take obvious steps to mitigate non-heritable variations in community function. For example, 386 we can repeatedly measure community function to increase measurement precision, thereby facilitating 387 selection ( Figure 5). We can also use the bet-hedging strategy so that lower-functioning communities 388 harboring desirable genotypes can have a chance to reproduce (Figure 3). We can use a cell sorter to 389 fix the cell number or biomass of member species in Newborns so that community function suffers less 390 non-heritable variations (Figure 3). 391 The need to suppress non-heritable variations in community function can have practical implications 392 that may initially seem non-intuitive. For example, when shared resource is non-limiting (to avoid 393 stationary phase), we must dilute a chosen Adult community to a fixed target biomass instead of by a 394 fixed-fold. This is because otherwise, selection would fail as we choose larger and larger Newborn size 395 instead of higher and higher f P (Methods, Section 7; Figure S23). 396 The definition of community function is also critical. If we had defined community function as Product  Intra-community selection versus inter-community selection 403 When improving a costly community function, intra-community selection and inter-community selection 404 are both important. Intra-community selection occurs during community maturation, and favors fast 405 growers. Inter-community selection occurs during community reproduction, and favors high community 406 function.

407
For production cost f P , intra-community selection favors low f P , while inter-community selection 408 favors f * P , the f P leading to the highest community function. Thus, when current f P < f * P , inter-For growth parameters (maximal growth rates, affinities for metabolites), depending on their evolutionary upper bounds, the two selection forces may or may not be aligned. For example, using parameters 413 in Table 1, improving growth parameters promoted community function (Figure S8A-C; Figure S24A). 414 This is because with these choices of evolutionary upper bounds, H could not evolve to grow so fast to 415 overwhelm M, Thus, with sufficient Resource and without the danger of species loss (Figure 2 bottom), 416 faster H and M growth resulted in more Byproduct, larger M populations, and consequently higher Prod-417 uct level. If H could evolve to grow faster than M, then increasing growth parameters could decrease 418 community function due to H dominance ( Figure 6A; Figure S17A; Figure S24B), although properly exe-419 cuted community selection can improve community function while promoting species coexistence ( Figure   420 6B, C).

421
Contrasting selection at different levels 422 From a "gene-centric" perspective, selection of individuals bears resemblance to selection of communities,  Figure S1B; Figure S7) Figure S25B). Optimal group function was indeed realized during selection ( Figure S25A). However, 481 optimal f P for group function is much lower than optimal f P for community function (f The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/264689 doi: bioRxiv preprint ratio: one could initiate Newborn community replicates and measure community functions using the most 498 precise method (e.g. cell-sorting during Newborn formation; many repeated measurements of community 499 function). Despite this, some levels of non-heritable variations in community function are inevitable due 500 to, for example, non-genetic phenotypic variations among cells [67] or stochasticity in cell birth and 501 death. If "noises" (variations among community replicates) are small compared to "signals" (variations

512
Our work suggests that if selection for a costly microbial community function should occur in nature, 513 then mechanisms for suppressing non-heritable variations in community function should be in place.

514
Future directions 515 Our work touched upon only the tip of the iceberg of community selection. We expect that certain 516 rules will be insensitive to details of a community. For example, reducing non-heritable variations in 517 community function and judicious bet-hedging can facilitate community selection. Regardless, many 518 fascinating questions remain. Here, we outline a few: 519 1. How might we best "hedge bets" when choosing Adult communities to reproduce? We have 520 chosen top ten communities to contribute an equal number of Newborns, but alternative strategies (e.g.

521
allowing higher-functioning Adults to contribute more Newborns) may work better. This aspect might 522 be explored through applying population genetics theories, which has considered the balance between 523 the strength of selection and variation among the individuals.  The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/264689 doi: bioRxiv preprint concentration. If during Newborn formation, any one species is lost by chance, then community function 541 would be stuck at local sub-optima. In this case, "bet-hedging" with community migration (mixing) 542 could recover the lost species and help community selection reach a higher optimum.   Note that cell-sorting only needs to be performed on several high-functioning communities, and thus 552 would not be cost prohibitive.
Total potential growth rate of M g M depends on the levels of Resource and Byproduct (R andB) 567 according to the Mankad-Bungay model [27] due to its experimental support: Figure S3). 1 − f P fraction of M growth is channeled to 569 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 571 factor, see below and Table 1).

572
ResourceR is consumed proportionally to the growth of M and H; ByproductB is released propor-573 tionally to H growth and consumed proportionally to M growth: In equations above, scaling factors are marked by "~", and will become 1 after scaling. Variables and 581 parameters with hats will be scaled and lose their hats afterwards. Variables and parameters without 582 hats will not be scaled. We scale Resource-related variable (R) and parameters (K M R ,K HR ,ĉ RM , 583 andĉ RH ) against R(0) (Resource supplied to Newborn), Byproduct-related variable (B) and parameters 584 (K M B andĉ BM ) against r B (amount of Byproduct released per H biomass grown), and Product-related 585 variable (P ) against r P (amount of Product made at the cost of one M biomass). For biologists who 586 usually think of quantities with units, the purpose of scaling (and getting rid of units) is to reduce the 587 number of parameters. For example, H biomass growth rate can be re-written as: 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 589 identical forms ( Figure S3). After scaling, the value of R(0) becomes irrelevant (1 with no unit). Similarly, Figure S4).

592
Thus, scaled equations are

16
. CC-BY-NC 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/264689 doi: bioRxiv preprint We have not scaled time here, although time can also be scaled by, for example, the community maturation time. Here, time has the unit of unit time (e.g. hr), and to avoid repetition, we often drop 595 the time unit. After scaling, values of all parameters (including scaling factors) are in Table 1, and 596 variables in our model and simulations are summarized in Table 2.

597
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 The because if a community has φ M (0) = φ M,SS at its Newborn stage, it has the same φ M (T ) = φ M,SS at 606 its Adult stage.

607
In our simulations, because we supplied the H-M community with abundant R to avoid stationary 608 phase, H grows almost at the maximal rate through T and releases B. If f P is not too large (f P < 0.4), 609 which is satisfied in our simulations, M grows at a maximal rate allowed by B and keeps B at a low 610 level. Thus, Eq. 14 is applicable and predicts the steady-state φ M,SS well (see Figure S26). Note that 611 significant deviation occurs when f P > 0.4. This is because when f P is large, M's biomass does not 612 grow fast enough to deplete B so that we cannot approximate B(T ) ≈ 0 anymore. (f P ), we can identify the optimal f * P associated with maximal community function ( Figure 2B).

621
For ancestral H, we set g Hmax = 0.25 (equivalent to 2.8-hr doubling time if we choose hr as the time 622 unit), K HR = 1 and c RH = 10 −4 (both with unit of R(0)) ( 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 "Evolved" column are used for most simulations and figures unless otherwise specified. 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 S25B. In Methods Section 2, we explained our parameter choices (including why we hold some parameters constant during evolution).

Symbols
Definition T Community maturation time, corresponding to the duration of a selection cycle t Time within a selection cycle, 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, The average biomass (length) of an individual M or H cell.
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 D The fold dilution when reproducing an Adult community n tot Total number of communities under selection in each cycle   Table 1.

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

629
Since we have assumed M and H to have the same affinity for R (Table 1), g M max must exceed g Hmax , 630 and M's affinity for Byproduct (1/K M B ) must be sufficiently large. Moreover, metabolite release and 631 consumption need to be balanced to avoid extreme ratios between metabolite releaser and consumer. can yield a maximum of 10 4 total biomass.

651
In some simulations (e.g. Figures 6, S8, S17), growth parameters (maximal growth rates g M max 652 and g Hmax and affinities for nutrients 1/K M R , 1/K M B , and 1/K HR ) and production cost parameter

653
(0 ≤ f P ≤ 1) were allowed to change from ancestral values during community maturation, since these 654 phenotypes have been observed to rapidly evolve within tens to hundreds of generations ([29, 30, 31, 32]).

655
For example, several-fold improvement in nutrient affinity and~20% increase in maximal growth rate 656 have been observed in experimental evolution [32, 30]. We therefore allowed affinities 1/K M R , 1/K HR , 657 and 1/K M B to increase by up to 3-fold, 5-fold, and 5-fold respectively, and allowed g Hmax and g M max changes are likely inconsequential. Below, we present evidence that within our parameter ranges (Table 1), improving growth parameters 690 generally improves community function. When f P is optimal for community function (f * P = 0.41), if we 691 fix four of the five growth parameters to their upper bounds, then as the remaining growth parameter 692 improves, community function increases (magenta lines in top panels of Figure S27). Moreover, mutants 693 with a reduced growth parameter are out-competed by their growth-adapted counterparts (magenta lines 694 in bottom panels of Figure S27).

695
When f P = f * P, M ono = 0.13 (optimal for M-monoculture function in Figure S25; the starting genotype 696 for most community selection trials in this paper), community function and individual fitness generally 697 increase as growth parameters improve (black dashed lines in Figure S27). However, when M's affinity for   (Eq. 16). Therefore M cells with lower affinity for R will grow faster than those with higher affinity 709 ( Figure S28). At the end of each cycle, the opposite is true ( Figure S28). As f P decreases, M diverts 710 more toward biomass growth and the first stage of B limitation lasts longer. Consequently, M can gain 711 a slightly higher overall fitness by lowering the affinity for R ( Figure S28A). when we only allow f P to change (compare Figure S29D with Figure 3A). . We term all these cases as "neutral" mutations.

729
Since a larger fraction of neutral mutations is equivalent to a lower rate of phenotype-altering muta-730 tions, our simulations define "mutation rate" as the rate of non-neutral mutations that either enhance 731 a phenotype ("enhancing mutations") or diminish a phenotype ("diminishing mutations"). Enhancing  Depending on the phenotype, the rate of phenotype-altering mutations is highly variable. Although The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/264689 doi: bioRxiv preprint is dysfunctional [85,86,84]. In our simulations, we assume a high, but biologically feasible, rate of 743 2 × 10 −3 phenotype-altering mutations per cell per generation per phenotype to speed up computation.

744
At this rate, an average community would sample~20 new mutations per phenotype during maturation. 745 We have also simulated with a 100-fold lower mutation rate. As expected, evolutionary dynamics slowed 746 down, but all of our conclusions still held ( Figure S18).

747
Among phenotype-altering mutations, tens of percent create null mutants, as illustrated by experimen-748 tal studies on protein, viruses, and yeast [33,34,35]. Thus, we assumed that 50% of phenotype-altering 749 mutations were null (i.e. resulting in zero maximal growth rate, zero affinity for metabolite, or zero f P ).

750
Among non-null mutations, the relative abundances of enhancing versus diminishing mutations are highly 751 variable in different experiments. It can be impacted by effective population size. For example, with a 752 large effective population size, the survival rate of beneficial mutations is 1000-fold lower due to clonal 753 interference (competition between beneficial mutations) [87]. The relative abundance of enhancing ver-754 sus diminishing mutations also strongly depends on the starting phenotype [33,80,78]. For example 755 with ampicillin as a substrate, the wild-type TEM-1 β-lactamase is a "perfect" enzyme. Consequently, The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/264689 doi: bioRxiv preprint where µ w (j) is short-hand notation for µ w ( s W T = j × 0.04) and so on. Our calculated µ ∆s (∆s) with 783 error bar of δµ ∆s is shown in Figure S6. 784 Our reanalysis demonstrated that distributions of mutation fitness effects µ ∆s (∆s) are largely con-785 served regardless of nutrient conditions and mutation types ( Figure S6B). In all cases, the relative fitness 786 changes caused by beneficial (fitness-enhancing) and deleterious (fitness-diminishing) mutations can be 787 approximated by a bilateral exponential distribution with means s + and s − for the positive and negative 788 halves, respectively. After normalizing the total probability to 1, we have: We fitted the Dunham lab haploid data (since microbes are often haploid) to Eq. 19, using 790 µ ∆s (i)/δµ ∆s (i) as the weight for non-linear least squared regression (green lines in Figure S6B). We 791 obtained s + = 0.050 ± 0.002 and s − = 0.067 ± 0.003.

792
Interestingly, exponential distribution described the fitness effects of deleterious mutations in an RNA 793 virus remarkably well [33].  [87,94]. 798 We have also simulated smaller average mutational effects based on measurements of spontaneous 799 or chemically-induced (instead of deletion) mutations. For example, the fitness effects of nonlethal 800 deleterious mutations in S. cerevisiae were mostly 1%~5% [35], and the mean selection coefficient 801 of beneficial mutations in E. coli was 1%∼2% [90,87]. As an alternative, we also simulated with 802 s + = s − = 0.02, and obtained the same conclusions ( Figure S19).

820
Based on the above experimental and theoretical literature, we modeled epistasis on f P in the 821 following manner. Let the relative mutation effect on f P be ∆f P = (f P,mut − f P ) /f P (note ∆f P ≥ −1).

822
Then, µ(∆f P , f P ), the PDF of ∆f P at the current f P value, is described by a form similar to Eq. 19: The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/264689 doi: bioRxiv preprint Here, s + (f P ) and s − (f P ) are respectively the mean ∆f P for enhancing and diminishing mutations 824 at current f P . We assigned s + (f P ) = s +init /(1 + g × (f P /f P, init − 1)), where f P, init is the f P of the 825 initial background in a community selection simulation ( f P, init = f * P, M ono = 0.13), s +init is the mean From Eq. 7 and 8, within ∆τ , After division, each mutable phenotype of each cell had a probability of P mut to be modified by a mutation (Methods, Section 4). As an example, let's consider mutations in f P . If a mutation occurred, 891 then f P would be multiplied by (1 + ∆f P ), where ∆f P was determined as below.

892
First, a uniform random number u 1 between 0 and 1 was generated. If u 1 ≤ 0.5, ∆f P = −1, 893 which represented 50% chance of a null mutation (f P = 0). If 0.5 < u 1 ≤ 1, ∆f P followed the 894 distribution defined by Eq. 20 with s + (f P ) = 0.05 for f P -enhancing mutations and s − (f P ) = 0.067 for 895 f P -diminishing mutations when epistasis was not considered (Methods, Section 4). In the simulation, 896 ∆f P was generated via inverse transform sampling. Specifically, C(∆f P ), the cumulative distribution 897 function (CDF) of ∆f P , could be found by integrating Eq. 19 from -1 to ∆f P : The two parts of Eq. 25 overlap at C(∆f

899
In order to generate ∆f P satisfying the distribution in Eq. 19, a uniform random number u 2 between 900 0 and 1 was generated and we set C(∆f P ) = u 2 . Inverting Eq. 25 yielded

905
If a mutation increased or decreased the phenotypic parameter beyond its bound (Table 1), the 906 phenotypic parameter was set to the bound value.

907
The above growth/death/division/mutation cycle was repeated from time 0 to T . Note that since  We obtained the same conclusion ( Figure S13, right panels).  Here we describe problems associated with two alternative definitions of community function and one 964 alternative method of community reproduction.
T is long enough for cells to double at least three or four times).

973
If we define community function as P (T )/M (T ) ≈ f P 1−f P , then higher community function requires 974 higher f P 1−f P or higher f P . However, if we select for very high f P , then M can go extinct ( Figure 2).

975
In our community selection scheme, the average total biomass of Newborn communities was set to a Newborn total biomass (instead of higher f p ). This will continue until Resource becomes limiting, and 980 then communities will get into the stationary phase. For groups or communities with a certain T g M dt, we can calculate f P optimal for community function 983 from Eq. ?? by setting We have If T g M dt 1, f P is very small, then the optimal f P for P (T ) is M grows faster in monoculture than in community because B is supplied in excess in monoculture while 988 in community, H-supplied Byproduct is initially limiting. Thus, T g M dt is larger in monoculture than in 989 community. According to Eq. 27, f * P ≈ 1/ T g M dt is smaller for monoculture than for community. i. In our simulations, cycles of selection were performed on a total of n tot = 100 communities with the indicated initial conditions. At the beginning of the first cycle, each Newborn had a total biomass of the target biomass (BM target =100; 60 M and 40 H each of biomass 1). In subsequent cycles, each Newborn's species ratio would be approximately that of its parent Adult. The amount of Resource in each Newborn was fixed at a value that could support a total biomass of 10 4 (unless otherwise stated). ii. The maturation time T was chosen so that for an average community, Resource was not depleted by time T (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). Once a cell's biomass had grown from 1 to 2, it divided into two identical daughter cells. Death occurred stochastically to individual cells (not depicted). After division, mutations (different shades of oval) occurred stochastically to change a cell's phenotypes (e.g. M's f P ). iii. At the end of a cycle, community functions (total Product P (T )) were ranked. iv. During community reproduction, high-functioning Adults were chosen and diluted into Newborns so that on average, each Newborn had a total biomass of approximately the target biomass BM target . A total of n tot = 100 Newborns were generated for the next selection cycle. 31 . CC-BY-NC 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It .  (46 H and 54 M cells). Note that at zero f P , no Product would be made; at f P = 1, M would go extinct. The maximal P * (T ) could not be further improved even if we allowed all growth parameters and f P to mutate ( Figure S10). Thus, P * (T ) is locally maximal in the sense that small deviation will always reduce P (T ). Ancestral f P (grey) is lower than f * P . The central question is: can community selection improve f P to f * P despite natural selection's favoring lower f P ? 32 . CC-BY-NC 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . . Community selection was not effective. Average f P and community function failed to improve to their theoretical optima, and community function poorly correlated with its heritable determinant f P (0). Black and magenta dots: un-chosen and chosen communities from one selection cycle, respectively. (D-F) Adults were chosen using the "top-dog" strategy, and a fixed H biomass and M biomass from the chosen Adults were sorted into Newborns. Community selection was successful. Community function also correlated with its heritable determinant f P (0). Here, Newborn total biomass BM (0) and fraction of M biomass φ M (0) were respectively fixed to BM target = 100 and φ M (T ) of the chosen Adult of the previous cycle. (G-I) When we chose the top ten Adults and let each reproduce 10 Newborns via pipetting, community function improved despite poor correlation between community function and its heritable determinant f P (0). For selection dynamics over many cycles, see Figure S14. (J-L) Evolution dynamics when maturation time was long (T = 20) such that most Resource was consumed by the end of T . Adults were chosen using the "top-dog" strategy, and reproduced via "pipetting". Community selection was successful due to high correlation between community function and its heritable determinant f P (0), assuming that variable duration in stationary phase would not introduce non-heritable variations in community function. Black, cyan and gray curves are three independent simulation trials. P (T ) was the average of P (T ) across all chosen Adults. f P (T ) was obtained by first averaging among M within each chosen Adult and then averaging across all chosen Adults.

33
. CC-BY-NC 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/264689 doi: bioRxiv preprint Figure 5: Community selection impeded by community function measurement noise can be rescued by "bet-hedging" acting in synergy with "cell-sorting". Adult communities were chosen to reproduce based on "measured community function P (T )" -the sum of actual P (T ) and an "uncertainty term" randomly drawn from a normal distribution with zero mean and standard deviations of 5% or 10% of the ancestral P (T ). Dynamics of average f P and average community function of the chosen Adult communities (f P (T ) and P (T )) are plotted. When uncertainty in community function measurement is low (5%), cell-sorting largely rescues ineffective community selection (A-D, left panels). When uncertainty in community function measurement is high (10%), both cell-sorting and bet-hedging are required (A-D, right panels). Black, cyan and gray curves are three independent simulation trials. P (T ) was averaged across the chosen Adults. f P (T ) was obtained by first averaging among M within each chosen Adult and then averaging across the chosen Adults.

35
. CC-BY-NC 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/264689 doi: bioRxiv preprint Figure 6: Effective community selection can encourage species coexistence. Here, the evolutionary upper bound for g Hmax (g * Hmax = 0.8) was larger than that for g M max (g * M max = 0.7), opposite to that in Figures 2-5. (A) When the "top-dog" strategy and "pipetting" were used to choose and reproduce Adult communities, M was almost outcompeted by H as H evolved to grow faster than M (third panel). Although M would ordinarily go extinct, community selection managed to maintain M at a very low level (bottom). This imbalanced species ratio resulted in very low community function (top). When community selection was effective, either using the "top-dog" strategy with "cell-sorting" (B), or the "bet-hedging" strategy with "pipetting" (C), community selection successfully improved community function and f P . Strikingly in both B and C, H's growth parameter did not increase to its evolutionary upper bound ( Figure S17B and C), allowing a balanced species ratio (bottom) and high community function (top). Resource supplied to Newborn communities here supports 10 5 total biomass to accommodate faster growth rates (and hence community function is larger than in other figures). Black, cyan and gray curves are three independent simulation trials. P (T ) and φ M (T ) (fraction of M biomass in Adult communities) were averaged across the chosen Adults. f P (T ) was obtained by first averaging among M within each chosen Adult and then averaging across all chosen Adults.

36
. CC-BY-NC 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/264689 doi: bioRxiv preprint 37 . CC-BY-NC 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/264689 doi: bioRxiv preprint 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 [103]. 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 chosen 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 [104]. 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 [18,19,20,105,106], 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 same-colored outlines in Figure 4A). (B) Artificial selection of mono-species groups has been successful [43, 45, 15]. 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 [62]. 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.

38
. CC-BY-NC 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/264689 doi: bioRxiv preprint  The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/264689 doi: bioRxiv preprint Figure S4: 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  The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/264689 doi: bioRxiv preprint Figure S6: Probability density functions of changes in relative fitness due to mutations (µ ∆s (∆s)). We derived µ ∆s (∆s) from the Dunham lab data [36] 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. 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 area of the hatched region drawn in the leftmost panel. Figure S7: Large Newborn group size or long maturation time allows non-contributors to accumulate and reduces inter-group variation. For simplicity, we modeled the growth of Newborn groups of M cells. From a Newborn biomass of 10 2 or 10 4 wild-type M cells, M population multiplied for 6 or 100 generations. Immediately following cell division, wild-type daughter cells mutated to noncontributors with a probability of 10 −3 . Wild-type and mutant cells followed exponential growth.The growth rate of wild-type cells was 0.87 times that of mutants. The fraction of biomass made up by mutants at each wild-type doubling is shown. Note the different scales. 41 Figure S8: Improved individual growth can promote community function. Here, we allowed mutations to alter M's f P and H and M's growth parameters. Communities are chosen using the "top-dog" strategy. (A-C) Community reproduction via pipetting (i.e. Newborn biomass and species composition can fluctuate). Community function P (T ) increased upon community selection (A). Since f P remained unchanged (B), this increase in P (T ) must be due to improved growth parameters (C). (D-F) Community reproduction via biomass sorting (i.e. fixed Newborn total biomass and species composition). In both cases, the five growth parameters increased to their respective evolutionary upper bounds (green dashed lines). Magenta dashed lines: f P optimal for community function and maximal community function P (T ) when all five growth parameters are fixed at their evolutionary upper bounds and φ M (0) is also optimal for P (T ). Black, cyan, and gray curves show three independent simulations. P (T ) is averaged across chosen Adults.ḡ M max ,ḡ Hmax , and f P are obtained by averaging within each chosen Adult and then averaging across chosen Adults. K SpeciesM etabolite are averaged within each chosen Adult, then averaged across chosen Adults, and finally inverted to represent average affinity. Note different x axis scales. The maximal growth rates (g M max and g Hmax ) have the unit of 1/time. Affinity for Resource (1/K M R , 1/K HR ) has the unit of 1/ R(0), where R(0) is the initial amount of Resource in Newborn. Affinity for Byproduct (1/K M B ) has the unit of 1/ r B , where r B is the amount of Byproduct released per H biomass produced. Product P has the unit of r P , the amount of Product released at the cost of one M biomass. More details on parameters can be found in Table 1. Figure S9: Community function declines to zero in the absence of inter-community selection for higher community function. When random Adults were chosen to reproduce ("pipetting"), natural selection favored zero f P (A). Consequently, P (T ) decreased to zero (B). Black, cyan and gray curves are three independent simulation trials. P (T ) was averaged across the two randomly chosen Adults. f P (T ) was obtained by first averaging among M within each randomly chosen Adult and then averaging across the two randomly chosen Adults. Figure S10: P * (T ) is a local optimum because it cannot be further improved. We started each Newborn community with total biomass BM (0) = 100, all five growth parameters at their evolutionary upper bounds, and f * P = 0.41 and φ * M (0) = 0.54 to achieve P * (T ). We then allowed all five growth parameters and f P to mutate while applying community selection. To ensure effective community selection ( Figure 3D-F), BM (0) was fixed to 100, and φ M (0) was fixed to φ M (T ) of the chosen Adult community from the previous cycle during community reproduction. We found that all five growth parameters remained at their respective evolutionary upper bounds. At the end of the first cycle (Cycle = 1 in insets), even though f P did not change, P (T ) had 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.72, compare with φ M,SS represented by the green dashed line in Figure 2A bottom panel). Later, over hundreds of cycles, f P gradually increased while P (T ) was still below maximal. This is because species composition gravitated toward steady state φ M,SS which deviated from the optimal φ * M (0) ([64]). Other legends are the same as Figure S8.
43 . CC-BY-NC 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/264689 doi: bioRxiv preprint  The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/264689 doi: bioRxiv preprint Figure S14: Bet-hedging can improve selection efficacy when community function experiences non-heriable variations, but not as effectively as reducing non-heritable variations in community function. Simulation setup is identical to that in Figure 3G-H, except that selection here lasts more cycles. Compared to Figure 3A-C, bet-hedging was more effective. However, compared to Figure 3D-F, bet-hedging did not improve community function to the same extent, even over 10 4 cycles. Black, cyan and gray curves are three independent simulation trials. P (T ) was averaged across the chosen Adults. f P (T ) was obtained by first averaging among M within each chosen Adult and then averaging across all chosen Adults. Figure S15: Bet-hedging strategies promoted community selection under a wide range of selection strengths. In a bet-hedging strategy, top k Adults each contributed 100/k Newborns into the next cycle. Here, Adults were reproduced (split) into Newborns via pipetting. When all Adults contributed one Newborn each, selection strength was zero and thus natural selection quickly reduced average f P and community function to zero (rightmost column). Black, cyan and gray curves are three independent simulation trials. P (T ) was averaged across the two chosen Adults. f P (T ) was obtained by first averaging among M within each chosen Adult and then averaging across all chosen Adults.

46
. CC-BY-NC 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/264689 doi: bioRxiv preprint Figure S16: The "top-dog" strategy is superior to the "bet-hedging" strategy when nonheritable variation in community function is low. 20 replicas of selection simulations were performed using either the "top-dog" strategy (black curves) or the "bet-hedging" strategy (top ten Adults chosen to reproduce; red curves). Community reproduction was through cell-sorting. Community functions improved slightly faster and to a slightly higher level using the "top-dog" strategy. Thus, when nonheritable variations in community function were suppressed, the "top-dog" strategy was superior to the "bet-hedging" strategy.

47
. CC-BY-NC 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/264689 doi: bioRxiv preprint Figure S17: Community function can be improved even if it is costly to both species. Identical to Figure 6, the evolutionary upper bound for g Hmax (g * Hmax = 0.8) was larger than that of g M max (g * M max = 0.7), opposite to that in Figure 3. (A) Chosen Adult communities were reproduced through pipetting such that both BM (0) and φ M (0) could stochastically fluctuate. Eventually, g Hmax and g M max evolved to their respective upper bounds, and thus g Hmax > g M max (compare i and iv). This would ordinarily lead to extinction of M. However, community selection managed to maintain M at a very low level ( Figure 6A bottom panel). (B, C) Adult communities were chosen using the top-dog strategy and reproduced through cell biomass sorting (B) or chosen using the bet-hedging strategy and reproduced through pipetting (C). Community selection worked in the sense that both f P and P (T ) improved over cycles ( Figure 6). Strikingly, the maximal growth rate of H g Hmax did not increase to its upper bound g * Hmax = 0.8, and H's affinity for Resource even decreased from the ancestral level in some cases. Here, Resource supplied to Newborn communities could support 10 5 total biomass to accommodate faster growth rate. Other legends are the same as Figure S8.
48 . CC-BY-NC 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/264689 doi: bioRxiv preprint  Figure 3, except that the null mutations were eliminated. Specifically, the distribution of mutation effects is specified by Eq. 19 where s + = 0.05 and s − = 0.067. f P as well as P (T ) were more stable compared to those in Figure 3. Black, cyan and gray curves are three independent simulation trials. P (T ) was averaged across the chosen Adults. f P (T ) was obtained by first averaging among M within each chosen Adult and then averaging across all chosen Adults.

50
. CC-BY-NC 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/264689 doi: bioRxiv preprint , mutational effects on f P depend on the current value of f P . 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.

51
. CC-BY-NC 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . Figure S21: Evolutionary dynamics of chosen Adults when epistasis is considered. When we incorporated different epistasis strengths (epistasis factor of 0.3 and 0.8), we obtained essentially the same conclusions as when epistasis was not considered (Figure 3). Other legend details can be found in The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/264689 doi: bioRxiv preprint Figure S23: Artificial selection on community function in excess Resource failed under fixed-fold pipetting dilution scheme. Excess Resource was supplied to each Newborn (R(0)/K M R = 10 6 ), and chosen Adults were reproduced via a fixed-fold (100-fold) pipetting dilution into Newborns. Because of pipetting, Newborns with larger total biomass will tend to be selected ( Figure 4F). Community selection quickly failed as Newborn total biomass increased exponentially (bottom) while non-producing M cells with f P = 0 quickly took over (top; Figure S7B). Black, cyan and gray curves are three independent simulation trials. P (T ) was averaged across chosen Adults. f P (T ) was obtained by first averaging among M within each chosen Adult and then averaging across all chosen Adults.

54
. CC-BY-NC 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/264689 doi: bioRxiv preprint Figure S24: Improving growth can improve or impair community function, depending on evolutionary upper bounds of growth parameters. Plotted here are plateaued community function after 1500 cycles when simulation did or did not allow mutations in growth parameters or f P . (A) When evolutionary upper bound for g Hmax (g * Hmax = 0.3) was lower than that of g M max (g * M max = 0.7), improving growth parameters improved community function. Compared to community function where no mutations were allowed (i), community function improved when growth parameters and f P were allowed to mutate (ii). Preventing mutations in growth parameters diminished community function improvement (iii). In this case, improved growth of M and H resulted in higher community function. Evolutionary dynamics are shown in Figure S8C. (B) When evolutionary upper bound for g Hmax (g * Hmax = 0.8) was larger than that of g M max (g * M max = 0.7), improving growth parameters could decrease community function. Compared to community function where no mutations were allowed (i), community function decreased when growth parameters and f P were allowed to mutate (ii). Preventing mutations in growth parameters diminished reduction in community function (iii). In this case, improved growth of M and H resulted in lower community function. Evolutionary dynamics are shown in Figure S17A. In B, Resource supplied to Newborn communities could support 10 5 total biomass to accommodate faster growth rate. In both A and B, community reproduction occurred through volumetric dilution via pipetting, and the top-dog strategy was used. Error bars are calculated form three independent selections.

55
. CC-BY-NC 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/264689 doi: bioRxiv preprint Then, maximal group function is achieved at f P = f * P, M ono = 0.13 (dashed line), lower than the optimal f P for the community function f * P = 0.41 ( Figure 2B). Here, the growth parameters of M are all fixed at their evolutionary upper bounds and P (T ) has the unit of r P . Figure S26: Comparison between the steady-state φ M,SS calculated from Eqs. 6-10 (black curve) and from Eq. 14 (red line).

56
. CC-BY-NC 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/264689 doi: bioRxiv preprint Figure S27: Improved maximal growth rates and nutrient affinities generally, but do not always, improve individual fitness and community function. In all figures, solid and dashed lines respectively represent calculations with f P = f * P = 0.41 (optimal for community function; Figure 2B) and f P = f * P, M ono = 0.13 (optimal for M monoculture production when Byproduct is in excess; Figure S25). Except for the indicated growth parameter, all other growth parameters were set to their respective upper bounds. Dynamics within one selection cycle is plotted. (A-D) Community function increases as the indicated growth parameter increases. For example in (A), all growth parameters except for g M max were set to their upper bounds. For each g M max , the steady-state φ M,SS was calculated using equations in Methods Section 1. This steady-state φ M,SS was then used to calculate P (T ). (F-I) The ratio between mutant population (whose indicated growth parameter was 10% lower than the upper bound) and growth-adapted population over maturation time T = 17. The decreasing ratio indicates that the mutant has a lower fitness compared to the growth-adapted cells. For example in (F), a Newborn community had 70 M and 30 H. 90% of M were growth-adapted and had upper bound g M max = 0.7 ("upper bound"). 10% of M had g M max = 0.63, 10% less than the upper bound ("mutant"). The ratio between "mutant" and "upper bound" cells declined over maturation time, indicating that mutant M cells had a lower fitness. (E, J) When f P = 0.13 (black dashed line) but not when f P = 0.41 (magenta line), increasing M's affinity for Resource (1/K M R ) slightly decreases individual fitness, and barely affects community function.

57
. CC-BY-NC 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It .   Figure S25). When we simulated community selection while allowing all growth parameters and f P to vary, M's affinity for R 1/K M R decreased slightly because at low f P = 0.13, M with a lower affinity for R (lower 1/K M R ) has a slightly improved individual fitness ( Figure S28). 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 were selected against by intra-community selection and by inter-community selection ( Figure S27). Other legends are the same as Figure S8.

58
. CC-BY-NC 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/264689 doi: bioRxiv preprint Figure S30: Different methods of pre-growth had limited impact on selection dynamics. (Top Panels) Histograms of the number of Newborn communities free of non-contributor M mutants when Newborn communities of the first cycle were inoculated from a single M monoculture (Left panel) or from independently-grown M monocultures (Right panel). To generate the histograms, the pre-growth and inoculation process was repeated 100 times. (Middle and Bottom panels) Improvement in f P (T ) and P (T ) was only slightly slower when Newborn communities from the first cycle were inoculated by the same M monoculture (Left panel) than by distinct monocultures (Right panel). Here we assumed that each M 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 (f P = 0) at a fixed probability of 10 −3 . If a non-null M cell has f P = 0.13, then it will grow at a rate 87% of that of a null cell. After~23 doublings, the M monocultures have on average~3% null mutants. 60 randomly-chosen M cells from the same monoculture or from distinct monocultures, together with 40 H cells, were used to inoculate each of the 100 Newborns for the first selection cycle. The top-dog strategy was used to choose Adults which were then reproduced via cell-sorting. Black, cyan and gray curves are three independent simulation trials. P (T ) was averaged across the two chosen Adults. f P (T ) was obtained by first averaging among M within each chosen Adult and then averaging across the two chosen Adults.

59
. CC-BY-NC 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/264689 doi: bioRxiv preprint 60 . CC-BY-NC 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/264689 doi: bioRxiv preprint 61 . CC-BY-NC 4.0 International license is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the author/funder. It . https://doi.org/10.1101/264689 doi: bioRxiv preprint