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

Multispecies microbial communities often display “community functions” arising 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 2 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 are carefully considered. These aspects include promoting species coexistence, suppressing noncontributors, choosing additional communities besides the highest functioning ones to reproduce, and reducing stochastic fluctuations in the biomass of each member species in Newborn communities. These considerations can be addressed experimentally. When executed effectively, community selection is predicted to improve costly community function, and may even force species to evolve slow growth to achieve species coexistence. Our conclusions hold under various alternative model assumptions and are therefore applicable to a variety of communities.


Introduction
Multispecies microbial communities often display community functions, defined as biochemical activities not achievable by member species in isolation. Many community functions have important commercial values. For example, a 6-species microbial community-but not any member species alone-cleared relapsing Clostridium difficile infections in mice [1]. Community functions arise from interactions by which an individual alters the physiology of another individual. Thus, to improve community functions, one could identify and modify interactions PLOS Biology | https://doi.org/10.1371/journal.pbio.3000295 June 25, 2019 1 / 52 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 [2,3]. In reality, this is no trivial task: each species can release dozens or more compounds, many of which may influence the partner species in diverse fashions [4,5,6,7]. From this myriad of interactions, one would then need to identify those critical for community function and modify them by altering species genotypes or the abiotic environment. One could also artificially assemble different combinations of species or genotypes to screen for high community function (e.g. [8,9,10]). However, some species may not be culturable in isolation. Moreover, the number of combinations becomes very large even when testing a moderate number of species and genotypes at various ratios, although recent advance has enabled massive parallel screening of synthetic microbial communities in droplets [11].
In an alternative approach, artificial selection of whole communities could be carried out over cycles to improve community trait [12][13][14][15][16][17][115][116][117][118][119][120] (reviewed in [18,19,20]; Fig 1A). A selection cycle starts with a collection of low-density "Newborn" communities with artificially imposed boundaries (e.g., inside culture tubes). These low-density communities are incubated for a period of time ("maturation") to form "Adult" communities. During maturation, community members multiply and interact with each other and possibly mutate, and the community function of interest develops. At the end of maturation, desired Adult communities are chosen to "reproduce" such that each is randomly partitioned into multiple Newborn communities to start the next cycle. Superficially, this process may seem straightforward since "one gets what one selects for." After all, artificial selection on individuals has been successfully implemented to obtain, e.g., proteins of enhanced activities ( [21,22,23]; S1 Fig). However, compared to artificial selection of individuals or monospecies groups, artificial selection of multispecies communities is more challenging (see detailed explanation in S1 Fig). For example, member species critical for community function may get lost during community reproduction.
The few attempts at community selection have generated interesting results. One theoretical study simulated artificial selection on multispecies communities based on the presence or absence of a member species [17]. Communities responded to selection, but only under certain conditions. In another theoretical study, multispecies communities responded to artificial selection for their ability to modify their abiotic environment in user-defined fashions [14]. In both cases, the response to selection quickly leveled off, and could be generated without mutations. Thus, community selection acted entirely on species types instead of new genotypes [14,17]. In experiments, complex microbial communities were selected for various traits [12,13,15,16,115,116]. For example, microbial communities selected to promote early or late flowering in plants were dominated by distinct species types [15]. However in other cases, a community trait may fail to improve despite selection and may improve even without selection [12,13,116].
Because communities used in these selection attempts were complex, much remains unknown. First, was the trait under selection a community function or achievable by a single species? If the latter, then community selection may not be needed, and the simpler task of selecting individuals or monospecies groups could be performed instead (S1 Fig). Second, did selection act solely on species types or also on newly arising genotypes? If selection acted solely on species types [14,15,17,118], then without immigration of new species to generate new variations, community function may quickly plateau [14,17,118]. If selection acted on genotypes, then community function could continue to improve as new genotypes evolve. Finally, why might a community trait sometimes fail to improve despite selection [12,13]?
In this study, we simulated artificial selection on communities with 2 defined species whose phenotypes can be modified by random mutations. Our goal is to improve a "costly" community function. A community function is costly if any community member's fitness is reduced by contributing to that community function (Fig 1B). Costly community functions can arise when microbes are engineered to contribute. Costly community functions are particularly challenging to improve: because contributors to community function grow slower than noncontributors, noncontributors will take over during community maturation. If all Adult communities are dominated by noncontributors, then community selection will fail. To improve a costly community function, intercommunity selection (which occurs once every cycle) must overcome intracommunity selection throughout community maturation (Fig 1B).
By simulating a simplified 2-species community, we could compare the efficacy of different selection regimens with relative ease, unlike experimental comparisons that would require concerted efforts from multiple labs [24]. Additionally, by analyzing evolving communitiestheir functions and member species phenotypes-we could begin to understand evolutionary dynamics during community selection. We also designed our simulations to mimic real lab experiments so that our conclusions could guide future experiments. For example, our simulations incorporated not only chemical mechanisms of species interactions (as advocated by [25,26]) but also experimental procedures (e.g., pipetting cultures during community reproduction). Model parameters, including species phenotypes, mutation rate, and distribution of mutation effects, were based on a wide variety of published experiments. Note that most previous models (e.g., [121]) focused on binary phenotypes (e.g., either contributing or noncontributing) and therefore could not model community function improvement driven by the evolution of quantitative phenotypes. We show that artificial community selection can improve a costly community function, but only after circumventing a multitude of failure traps. In discussions, we will elaborate on (i) challenges and solutions to community selection; (ii) the tension between intracommunity selection versus intercommunity selection; (iii) similarities and dinstinctions among individual selection, group selection, and community selection; (iv) why optimizing monocultures may not lead to optimal community function; (v) implications of our work; and (vi) future directions.

Results
We will first introduce the subject of our community selection simulation: a commensal 2-species community that converts substrates to a valued product. We will then define community function and describe how we simulate artificial community selection. Using simulation results, we will demonstrate critical measures that make community selection effective. Finally, we show that our conclusions are robust under alternative model assumptions, applicable also to mutualistic communities and communities whose member species may not coexist. To avoid confusion, we will use "community selection" or "selection" to describe the entire process of artificial community selection (community formation, growth, selection, and reproduction), and use "choose" or "intercommunity selection" to refer to the step in which the experimentalist decides which communities will reproduce.

A Helper-Manufacturer community that converts substrates into a product
Motivated by previous successes in engineering 2-species microbial communities that convert substrates into useful products [27,28,29], we numerically simulated selection of such communities.
In our community, Manufacturer M can manufacture Product P of value to us (e.g., a biofuel or a drug) at a fitness cost to self, but only if assisted by Helper H (Fig 1C). Specifically, Helper but not Manufacturer can digest an agricultural waste (e.g., cellulose), and as Helper grows biomass, Helper releases Byproduct B at no fitness cost to itself. Manufacturer requires H's Byproduct (e.g., carbon source) to grow. In addition, Manufacturer pays the cost of f P (0 � f P � 1) fraction of its potential growth to make Product P while using the rest (1 − f P ) for its biomass growth. Both species also require a shared Resource R (e.g., nitrogen). Thus, the 2 species together-but not any species alone-can convert substrates (agricultural waste and Resource) into Product.
We define community function as the total amount of Product accumulated as a lowdensity Newborn community grows into an Adult community over maturation time T, i.e., P(T). In the Discussion section, we explain problems associated with an alternative definition of community function (e.g., per capita production; Methods Section 7; S2 Fig). We will initially focus on the scenario in which community function is not costly to Helpers but incurs a fitness cost of f P to M. Later, we will show that our conclusions also hold when community function is costly to both H and M. Below, we will describe how we simulated community selection, followed by how we chose parameters of species phenotypes and parameters of selection regimen.

Simulating community selection
We simulated 4 stages of community selection ( Fig 1D) as follows: (i) forming Newborn communities, (ii) Newborn communities maturing into Adult communities, (iii) choosing highfunctioning Adult communities, and (iv) reproducing the chosen Adult communities by splitting each into multiple Newborn communities of the next cycle. Our simulation was individual-based. That is, it tracked phenotypes and biomass of individual H and M cells in each community as cells grew, divided, mutated, or died. Our simulations also tracked dynamics of chemicals (including Product) in each community and accounted for actual experimental steps such as pipetting cultures during community reproduction. Below is a brief summary of our simulations, with more details in Methods (Section 6).
Each simulation started with n tot (= 100) number of Newborn communities. Each Newborn community always started with a fixed amount of Resource R(0). Agricultural waste was always supplied in excess and thus did not enter our equations. In the first cycle, each Newborn community had a total biomass equal to a target value BM target (= 100; 60 Manufacturers and 40 Helpers each of biomass 1). See Discussion for problems associated with not having a biomass target.
During proportional to its potential growth rate. M grew biomass at (1 − f P ) fraction of potential growth rate (Eq 24), and released Product at a rate proportional to f P fraction of potential growth rate (Eq 25). Once an H or M cell's biomass grew from 1 to 2, it divided into 2 cells of equal biomass with identical phenotypes, thus capturing experimental observations of continuous biomass increase (S5 Fig) and discrete cell division events [31]. At any time, H and M cells died stochastically at a constant death rate. Although mutations can occur during any stage of the cell cycle, we assigned mutations immediately after cell division, such that each mutable phenotype of both cells mutated independently.
Mutable phenotypes included H and M's "growth parameters" (maximal growth rates in excess nutrients; affinities for nutrients), and M's cost f P . These phenotypes have been observed to rapidly change during evolution [32,33,34,35]. Mutated phenotypes could range between 0 and their respective evolutionary upper bounds. Among mutations that alter phenotypes (denoted "mutations"), on average, half abolished the function (e.g., zero growth rate, zero affinity, or zero cost f P ) based on experiments on green fluorescent protein, viruses, and yeast [36,37,38]. Effects of the other 50% of mutations were bilateral-exponentially distributed, enhancing or diminishing a phenotype by a few percent, based on our reanalysis of published yeast data sets (S6 Fig) [39]. We held death rates constant, because death rates were much smaller than growth rates and therefore mutations in death rates would be inconsequential. We also held release and consumption coefficients constant. This is because, e.g., the amount of Byproduct released per H biomass generated is constrained by biochemical stoichiometry.
At the end of community maturation time T, we compared community function P(T) (the total amount of Product accumulated in the community by time T) for each Adult community. We chose high-functioning Adults to reproduce. Each chosen Adult was randomly partitioned into Newborns with target total biomass BM target . For example, if the chosen Adult had a total biomass of 60BM target , then each cell would be assigned a random integer between 1 and 60, and those cells with the same random integer would be allocated to the same Newborn. Experimentally, this is equivalent to volumetric dilution using a pipette ("pipetting"). Thus, for each Newborn, the total biomass and species ratio fluctuated around their expected values in a fashion associated with pipetting (Methods Section 9). From top-functioning Adults, a total of n tot Newborns were obtained to enter the next selection cycle (see below for details).

Choosing species: Enhancing species coexistence
In order to improve community function through community selection, species should ideally coexist throughout selection cycles. That is, all species should grow at a similar average growth rate within each cycle. Furthermore, species ratio should not be extreme because otherwise, the low-abundance species could be lost by chance during Newborn formation. Species coexistence at a moderate ratio has been experimentally realized in engineered communities [27,28,40,41].
To achieve species coexistence at a moderate ratio in the H-M community, three considerations need to be made. First, M's cost f P for making Product must not be too large, otherwise M would always grow slower than H and thus eventually go extinct (Fig 2A, bottom). Second, H and M's growth parameters (maximal growth rates in excess nutrients; affinities for nutrients) must be balanced. This is because, upon Newborn formation, H can immediately start to grow on agricultural waste and Resource, whereas M cannot grow until H's Byproduct has accumulated to a sufficiently high level. Thus, to achieve coexistence, M must grow faster than H at some point during community maturation. Third, to achieve a moderate steady-state species ratio, metabolite release and consumption need to be balanced [40]. Otherwise, the ratio between releaser and consumer can be extreme.
Based on these considerations and published measurements on Saccharomyces cerevisiae and Escherichia coli, we chose H and M's ancestral growth parameters and their evolutionary upper bounds, as well as release, consumption, and death parameters (Table 1, Methods Section 2 and Section 3). Our choice of parameters ensured that throughout evolution, different species ratios would converge toward a moderate steady-state value during community maturation (Fig 2A, top). Note that if species were not chosen properly, selection might fail due to insufficient species coexistence (Fig 6A), although we will demonstrate that under effective community selection, requirements on species coexistence could be relaxed (Fig 6B-6D).

Choosing selection regimen parameters: Avoiding known failure modes
After choosing member species with appropriate phenotypes, we need to consider the parameters of our selection regimen ( Fig 1D). These parameters include the total number of communities under selection (n tot ), the number of Adult communities chosen to reproduce (n chosen , the "bottleneck size" when choosing a fraction of Adults to reproduce), Newborn target total biomass (BM target , the "bottleneck size" when splitting an Adult into Newborns), the amount of Resource added to each Newborn (R(0)), the amount of mutagenesis that controls the rate of phenotype-altering mutations (μ), and maturation time (T). Compared to the well-studied problem of group selection in which the unit of selection is a monospecies group [42][43][44][45][46][47][48][49][50][51][52][53][54][55][56], community selection is more challenging (Discussion; S1 Fig). However, the 2 types of selections do share some common aspects (Discussion; S1 Fig). Thus, we can apply group selection theory, together with other practical considerations, to better design community selection regimen.
Let's begin by considering the total number of communities under selection. The highest community function achieved among a large number of communities will likely exceed that achieved among a small number of communities. However, experimentally achieving a very large number of communities can be challenging. We chose a total of 100 communities (n tot = 100). n chosen , the number of Adults chosen by the experimentalist to reproduce, reflects intercommunity selection strength. Since the top-functioning Adult is presumably the most desirable, we reproduced it into as many Newborns as possible and then reproduced the second best,  (Table 1, "Preadapted"). (A) Stable species coexistence at moderate to low cost. Bottom: When f P , the fraction of potential growth Manufacturer diverts for making Product, is high (e.g., f P = 0.7), M will eventually go extinct (i.e., fraction of M < 1 � BM target ). Top: At moderate and low cost f P (e.g., f P = 0.4 and f P = 0.1), H and M can stably coexist. That is, different initial species ratios will converge to a steady-state value. At the end of the first cycle (time T = 17), Byproduct and Resource were reset to the initial conditions at time zero (0 and 10 4 , respectively), and total biomass was reduced to the target value BM target = 100, while the fraction of M biomass ϕ M remained the same as that of the parent community. See main text for how values of maturation time and Resource were chosen. (B) Optimal community function occurs at an intermediate cost f P . Community functions at various combinations of f P and fraction of M biomass (out of BM target = 100 total biomass) were computed by integrating Eqs 6-10. Maximal community function P(T) is achieved at an intermediate cost f � P ¼ 0:41 (magenta dashed line) when Newborn species composition is also optimal (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 (S10 Fig). Thus, P � (T) is locally maximal in the sense that small deviation will always reduce P(T). Ancestral f P (gray) is lower than f � P . The central question is this: can intercommunity selection improve ancestral f P to f � P despite intracommunity selection favoring lower f P ? The Matlab code can be found in S1 Code. H, Helper; M, Manufacturer. https://doi.org/10.1371/journal.pbio.3000295.g002 etc., until we obtained n tot Newborn communities for the next cycle (the "top-dog" strategy). In most simulations, the top-functioning community contributed approximately 60-70 Newborns. We then reproduced the next highest-functioning Adult in the same way and randomly chose enough (approximately 30-40) Newborns so that a total of n tot = 100 Newborns were generated for the next selection cycle. Later, we will compare the top-dog strategy with other strategies employing weaker intercommunity selection strengths.
Whether we mix chosen Adults before splitting them into Newborns could impact the efficacy of community selection. In our simulations, chosen Adults were not mixed to limit the spread of noncontributors. In other words, we created a spatially-structured environment, which is known to discourage noncontributors [57][58][59][60].
If the mutation rate is very low, then community function cannot rapidly improve. If the mutation rate is very high, then noncontributors will be generated at a high rate, and as the fast-growing noncontributors take over during community maturation, community function will likely collapse. Here, we chose μ, the rate of phenotype-altering mutations, to be biologically realistic (0.002 per cell per generation per phenotype, which is lower than the highest values observed experimentally; Methods Section 4).
If Newborn target total biomass BM target is very large, or if the number of doublings within maturation time T is very large, then noncontributors will take over in all communities during maturation (S7 Fig, compare B-D with A), as predicted by group selection theory. On the other hand, if BM target and the number of generations within T are both very small, then mutations will be rare within each cycle, and many cycles will be required to improve community function. Finally, if BM target is very small, then a member species might get lost by chance

Community selection may not be effective under conditions reflecting common lab practices
In initial simulations, we only allowed M's cost f p to be modified by mutations, and we fixed H and M's growth parameters (maximal growth rates in excess metabolites; affinities for metabolites) to their evolutionary upper bounds. Such a simplification is justified with our particular parameter choices (Table 1)  Later, we will show an opposite case in which growth parameters were kept below their evolutionary upper bounds by effective community selection (Fig 6). A great advantage can be gained by allowing only cost f p to mutate: we can now calculate the theoretical maximal community function P � (T) and its associated optimal cost f � P ð¼ 0:41Þ and optimal species ratio at a fixed total Newborn biomass ( Fig 2B).
We started simulations with M's cost f p lower than the optimal value f � P . Could intercommunity selection increase f p to f � P , despite intracommunity natural selection favoring lower f p ? As expected, in control simulations in which Adult communities were randomly chosen to reproduce, community function was driven to zero by intracommunity natural selection as fast-growing nonproducing M took over (S9 Fig).
When we chose Adults using the top-dog strategy (starting from the top-functioning Adult followed by the runners-up) and split them into Newborns as if via pipetting, M's cost f p and community function P(T) did not decline to zero but they barely improved, and both were far below their theoretical optima (Fig 3A and 3B).

Common lab practices can generate sufficiently large nonheritable variations in community function that interfere with selection
Why did community selection fail to increase M's cost f P and community function? One possibility is that community function was not sufficiently heritable from one cycle to the next (S1 Fig). We therefore investigated the heredity of community function by examining the heredity of community function determinants.
Community function P(T) was largely determined by phenotypes of cells in the Newborn community. This is because maturation time was sufficiently short (approximately 6 doublings), and therefore newly arising genotypes could not rise to high frequencies within one cycle to significantly affect community function. Because all phenotypes except for f P were fixed, community function had 3 independent determinants: Newborn's total biomass BM(0), Newborn's fraction of M biomass ϕ M (0), and Newborn's average M cost � f P ð0Þ (f P averaged across all M cells in the Newborn). Note that the first 2 values can stochastically fluctuate as expected from pipetting; e.g., to pipette 100 cells, one can end up pipetting between 70 and 130 cells even with a precise pipette.
A community function determinant is considered heritable if it is correlated between Newborns of one cycle (Fig 4A, bottom row) and their respective progeny Newborns in the next cycle ( Fig 4A, color-matched top row). Among the 3 determinants, � f P ð0Þ was heritable ( Fig  4B): if a Newborn community had a high average f P , so would the mature Adult community and Newborn communities reproduced from it. On the other hand, Newborn total biomass BM(0) was not heritable (Fig 4C). This is because when an Adult community reproduced via pipette dilution, the dilution factor was adjusted so that the total biomass of a progeny Newborn community was, on average, the target biomass BM target . Newborn's fraction of M biomass ϕ M (0), which fluctuated around ϕ M (T) of its parent Adult, was not heritable either ( Fig  4D). This is because, regardless of the species composition of Newborns, Adults would have similar steady-state species composition (Fig 2A top panel), and so would their offspring Newborns.
In successful community selection, variations in community function should be mainly caused by variations in its heritable determinants. However, community function P(T) weakly correlated with its heritable determinant � f P ð0Þ but strongly correlated with its nonheritable determinants (Fig 4E-4G). For example, the Newborn that would achieve the highest function had a below-median � f P ð0Þ (left magenta dot in Fig 4E) but had high total biomass BM(0) and low fraction of M biomass ϕ M (0) (Fig 4F and 4G). In other words, variation in community function is largely nonheritable, because it largely arises from variation in nonheritable determinants.
The reason for strong correlations between community function and its nonheritable determinants became clear by examining community dynamics. Recall that to avoid stationary phase, we had chosen maturation time so that Resource would be in excess by the end of maturation. Thus, a "lucky" Newborn community starting with a higher-than-average total biomass would convert more Resource to Product (dotted lines in top panels of S11 Fig). Similarly, if a Newborn started with higher-than-average fraction of Helper H biomass, then H would  (Fig 3B), community function P(T) correlates weakly with heritable determinant but strongly with nonheritable determinants. Each dot represents one community. Magenta dots: "successful" Newborns that achieved the highest community function at adulthood and therefore were chosen to reproduce in the top-dog strategy. The Matlab code for B-D can be found in S2 Code, and the data for E-G can be found in S1 Data. Adult, Adult community; M, Manufacturer; Newborn, Newborn community. https://doi.org/10.1371/journal.pbio.3000295.g004 Artificial selection for microbial community functions: Challenges and solutions produce higher-than-average Byproduct, which meant that M would endure a shorter growth lag, grow more, and make more Product (dotted lines in bottom panels of S11 Fig).
To summarize, when community function significantly correlated with its nonheritable determinants (Fig 4F and 4G), community selection failed to improve community function ( Fig 3B).

Reducing nonheritable variations in an experimentally feasible manner promotes artificial community selection
Reducing nonheritable variations in community function should facilitate community selection. One possibility would be to reduce stochastic fluctuations in nonheritable determinants. For example, a cell sorter could allocate to each Newborn community a fixed biomass or cell number from each species, if different species have different fluorescence [61]. Indeed, in simulations, when we fixed Newborn's total biomass and species fraction ("cell sorting"; Methods, Section 6), community function became strongly correlated with its heritable determinant ( Fig  3F). In this case, average cost � f P and community function P(T) both increased under selection (Fig 3D and 3E) to near the optimal. Improvement was not seen if either Newborn total biomass or species fraction was allowed to fluctuate stochastically (S12A Fig). Community function also improved if fixed numbers of H and M cells (instead of biomass) were allocated into each Newborn, even though each cell's biomass fluctuated between 1 and 2 (S13C Fig; Methods, Section 6).
Nonheritable variations in community function could also be curtailed by reducing the dependence of community function on nonheritable determinants. For example, we could extend the maturation time T to nearly deplete Resource. In this selection regimen, Newborns would still experience stochastic fluctuations in total biomass and species composition. However, "unlucky" communities would have time to "catch up" as "lucky" communities wait in stationary phase. Indeed, with this extended maturation time T, community function became strongly correlated with its heritable determinant � f P ð0Þ; and community function improved without having to fix Newborn total biomass or species composition (Fig 3J-3L). However, at long maturation time, nonheritable variations in community function could still arise from stochastic fluctuations in the duration of stationary phase (which could affect cell survival or recovery time in the next selection cycle). Thus, for most simulations, we used short maturation time unless stated otherwise.

Reducing inter-community selection strength can promote selection when nonheritable variations hinder selection
Since the highest community function may not correspond to the highest f P (Fig 3C), we examined whether a "top-tier" strategy might outperform the top-dog strategy. In the top-tier strategy, we chose, for example, the top 10% Adults, allowing each to reproduce 10 Newborns. When we partitioned Adults into Newborns as if using pipetting, although community function failed to improve under the top-dog strategy, it improved somewhat under the toptier strategy (Fig 3, compare G-I with A-C). In fact, the top-tier strategy improved community function under a wide range of selection strengths (e.g., top 5% to top 50%; S15 Fig). When we chose top 2% (2 Adults, each contributing 50 Newborns), community function did not improve (S15 Fig). When we chose all 100 Adults (each contributing 1 Newborn), community function declined to zero as expected (S15 Fig), because there was no intercommunity selection.
The superiority of top-tier over top-dog strategy rests on giving "unlucky" Adults a chance to reproduce. We reached this conclusion by noting that if we minimized nonheritable variations in community function (by fixing total biomass and species fraction in Newborns), then top-dog is superior to top-tier strategy (S16 Fig As expected, community function measurement noise-another source of nonheritable variation-interferes with community selection (compare Fig 3A and 3B with Fig 5A; Fig 3D  and 3E with Fig 5C). In this case, community selection can be improved by leveraging both the top-tier strategy and cell sorting (Fig 5, compare B and C with D). This is presumably because (i) cell sorting reduces nonheritable variations in community function and (ii) the top-tier strategy gives "unlucky" communities suffering unfavorable community function measurement errors a chance to reproduce.

Effective community selection can enforce species coexistence
Properly executed community selection could even improve the functions of communities whose member species may not coexist. Consider an H-M community in which, unlike the H-M community we have considered so far, H had the evolutionary potential to grow much faster than M. In this case, high community function not only required M to pay a fitness cost of f P but also required H to pay a fitness cost by growing sufficiently slowly to not outcompete M.
We started community selection at ancestral growth parameters and allowed them and f P to mutate. When community selection was ineffective (top-dog with pipetting; Fig 6A), H's maximal growth rate evolved to its upper bound and exceeded M's maximal growth rate ( . This is because if H's maximal growth rate in a community had evolved to be too high, then H would drive M to low abundance, and the resulting low community function would be disfavored by intercommunity selection. Because H was constrained to grow slower than M, H and M can coexist at a moderate ratio, and community function improved (Fig 6B-6D, row 1). Note that, unlike Fig 5, here top-tier and cell sorting did not show synergism. This is because when nonheritable variation in community function is minimized by cell sorting, the top-tier strategy does not seem to be helpful (S16 Fig).

Robust conclusions under alternative model assumptions
We have demonstrated that when selecting for high H-M community function, seemingly innocuous experimental procedures (e.g., choosing the top-functioning Adults and pipetting portions of them to form Newborns) could be problematic. Instead, more precise procedures (e.g., cell sorting) or moderate intercommunity selection strength (e.g., the top-tier strategy) might be required. Our conclusions held when we used a much lower mutation rate (2 × 10 −5 instead of 2 × 10 −3 mutation per cell per generation per phenotype, S18 Fig), although lower mutation rate slowed down community function improvement. Our conclusions also held when we used different distributions of mutation effects (S19 Fig) or incorporated epistasis (i.e., a non-null mutation would likely reduce f P if the current f P was high and increase f P if the current f P was low; S20 and S21 Figs; Methods Section 5).
To further test the generality of our conclusions, we simulated community selection on a mutualistic H-M community. Specifically, we assumed that Byproduct was inhibitory to H. Thus, H benefited M by providing Byproduct, and M benefited H by removing the inhibitory Byproduct, similar to the syntrophic community of Desulfovibrio vulgaris and Methanococcus maripaludis [62]. We obtained similar conclusions in this mutualistic H-M community (S22 Fig). We have also shown that similar conclusions hold for communities in which member species may not coexist (Fig 6).
In summary, our conclusions seem general under a variety of model assumptions and apply to a variety of communities. Ineffective selection due to community function measurement noise can be rescued by the top-tier strategy 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 a "noise 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 community function measurement noise is low (5%), cell sorting largely rescues ineffective community selection (A-D, left panels). When community function measurement noise is high (10%), both cell sorting and top-tier strategy are required (A-D, right panels). Black, cyan, and gray curves represent 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. The simulation codes can be found in S3 Code, and the data can be found in S2 Data. Adult, Adult community; M, Manufacturer. https://doi.org/10.1371/journal.pbio.3000295.g005

Discussion
A desired community function can be attained by identifying appropriate combinations of species types. For example, by using cellulose as the main carbon source in a process called "enrichment," Kato and colleagues obtained a community consisting of a few species that together degrade cellulose [63]. However, if we solely rely on species types, then without a Artificial selection for microbial community functions: Challenges and solutions constant influx of new species, community function will likely level off quickly [14,17,118]. Here, we consider artificial selection of communities with defined member species so that improvement of community function requires new genotypes that contribute more toward community function at a higher cost to itself.

Community selection can be challenging but is feasible
Artificial selection of whole communities to improve a costly community function requires careful considerations. These considerations include species choice (Figs 2 and 6), mutation rate (compare Fig 3 and S18 Fig) (Fig 3, S15 and S16 Figs), how we reproduce Adults (e.g., pipetting versus cell sorting; Fig 3), and the uncertainty in community function measurements (Fig 5).
Many of these considerations face dilemmas. For example, a large Newborn size (BM target ) might lead to reproducible takeover by noncontributors (S7 Fig), but a small Newborn size would mean that large nonheritable variations in community function can readily arise and interfere with selection unless special measures are taken (Fig 3).
We can take obvious steps to mitigate nonheritable variations in community function. For example, we can repeatedly measure community function to increase measurement precision, thereby facilitating selection ( Fig 5). We can also use the top-tier strategy so that unlucky communities harboring desirable genotypes can have a chance to reproduce (Figs 3 and 5). Note that a top-tier strategy had often been implemented during artificial individual selection [114] and artificial community selection [119] to maintain variations among entities. Although some community functions (such as steady-state species ratio or steady-state growth rate of mutualistic communities) are not sensitive to fluctuations in Newborn biomass compositions [40,64]; for those that are, we can use a cell sorter to fix Newborn species biomass compositions to reduce nonheritable variations in community function (Figs 3 and 5).
The need to suppress nonheritable variations in community function can have practical implications that may initially seem nonintuitive. For example, when shared Resource is nonlimiting (to avoid stationary phase), we must dilute a chosen Adult community to a fixed target biomass instead of by a fixed fold. This is because otherwise, selection would fail as we choose larger and larger Newborn size instead of higher and higher f P (Methods, Section 7; S23 Fig).
The definition of community function is also critical. If we had defined community function as Product per M biomass in the Adult community P(T)/M(T) (which is approximately proportional to f P 1À f P : see Methods Section 7), then we would be selecting for higher and higher

Intracommunity selection versus intercommunity selection
Intracommunity selection and intercommunity selection are both important. Intracommunity selection occurs during community maturation and favors fast growers. Intercommunity selection occurs prior to community reproduction and favors high community function.
For M's cost f P , intracommunity selection favors low f P , while intercommunity selection favors f � P , the f P value leading to the highest community function. Thus, when current f P < f � P , intercommunity selection runs against intracommunity selection. When current f P > f � P , intracommunity and intercommunity selections are aligned.
For growth parameters (maximal growth rates, affinities for metabolites), depending on their evolutionary upper bounds, intercommunity selection may or may not be aligned with intracommunity selection. For example, using parameters in Table 1 (Fig 3), improving growth parameters promoted community function (S8A- S8C Fig; S24A Fig). This is because with these choices of evolutionary upper bounds, H could not evolve to grow so fast to overwhelm M. Thus, with sufficient Resource and with species coexistence (Fig 2A top), faster H and M growth resulted in more Byproduct, larger M populations, and consequently higher Product level. If H could evolve to grow faster than M, then increasing growth parameters could decrease community function due to H dominance ( Fig 6A; S17A Fig; S24B Fig). Note that even in this case, properly executed community selection can promote species coexistence and improve community function (Fig 6B-6D).

Contrasting selection at different levels
Selection of individuals bears some resemblance to selection of communities. Whereas community function relies on interactions between different species, an individual's fitness relies on interactions between different genes. To ensure sufficient heredity between an individual and its offspring, elaborate cellular mechanisms have evolved. They include cell cycle checkpoints to ensure accurate DNA replication and segregation [65], small RNA-mediated silencing of transposons [66], and clustered regularly interspaced short palindromic repeats and the CRISPR-associated protein (CRISPR-Cas) degradation of foreign viral DNA [67]. In community selection, heredity-enhancing mechanisms such as stable species ratio (Fig 2A top panel) could already be in place due to ecological interactions or arise due to evolution (e.g., mutations that affect ecological interactions). If endosymbiosis should evolve in response to community selection (i.e., one microbe stably living inside another microbe much like chloroplasts living inside plant cells), then community selection would transition to individual selection.
Group selection is connected to individual selection and community selection. Group selection, and in a related sense kin selection [42][43][44][45][46][47][48][49][50][51][52][53][54][55]68], have been extensively examined to explain, e.g., the evolution of traits that lower individual fitness but increase the success of a group (e.g., sterile ants helping the survival of an ant colony). Note that the term "group selection" has often been used to describe individual selection in spatially structured populations without group births or deaths, although such usage has been criticized [69]. Artificial group selection can sometimes be viewed as artificial individual selection. For example, when Newborn groups start with a single founder producing a product of interest, then artificial group selection for high group production is equivalent to artificial individual selection for the founder's ability to produce over time as it grows into a population. On the other hand, if group function relies on synergistic interactions between populations with distinct phenotypes, then group selection can be thought of as a special case of community selection. The difference is that, unlike community function, group function can arise as a single founder multiplies and differentiates into distinct populations (e.g., filamentous cyanobacteria growing and differentiating into nitrogen-fixing cells and photosynthetic cells [70]).
Group selection and community selection display additional similarities and distinctions. First, group selection and community selection are similar in that Newborn size must not be too large [71,72], and maturation time must not be too long. Otherwise, all entities (groups or communities) will accumulate noncontributors in a similar fashion, and this lack of variation among entities impedes selection (Price equation [56]; S1B Fig; S7 Fig). Second, species interactions in a community could drive species composition to a value suboptimal for community function [73]. A similar constraint could also occur during artificial group selection if the founder genotype gives rise to interacting subpopulations. Otherwise, the problem of suboptimal composition does not exist for group selection. Finally, in group selection, when a Newborn group starts with a small number of individuals (e.g., one individual), a fraction of Newborn groups of the next cycle will be highly similar to the original Newborn group (S1B Fig, bottom panel). This heredity facilitates group selection. In contrast, when a Newborn community starts with a small number of total individuals, large stochastic fluctuations in Newborn composition can interfere with community selection (Figs 3 and 4). In the extreme case, a member species may even be lost by chance. Even if a fixed biomass of each species is allocated into Newborns during community selection, heredity is much reduced due to random sampling of genotypes from multiple species. For example, if Newborn communities start with one contributor from each of the 2 species and if the highest-functioning Adult community has accumulated 50% noncontributors in each species, then only 50%×50% = 25% Newborn communities of the next cycle will be similar to the original Newborns. In contrast, if Newborn groups are initiated with a single contributor and if the highest-functioning Adult group has accumulated 50% noncontributors, then 50% Newborn groups of the next cycle will be similar to the original Newborn.

Community function may not be maximized through pre-optimizing member species in monocultures
If we know how each member species contributes to community function, might we pre-optimize member species in monocultures before assembling them into high-functioning communities? This turns out to be challenging due to the difficulty of recapitulating community dynamics in monocultures. For example, artificial group selection on M failed to increase M's cost f P to f � P optimal for community function. Specifically, we started simulations with n tot of 100 Newborn M groups, each inoculated with one M cell (to facilitate group selection, S1B However, the associated optimal f P was much lower than that for community function (f � P ¼ 0:41; see Methods Section 8 for an explanation). Thus, optimizing monoculture activity does not necessarily lead to optimized community function.

Implications of our work
Our work sheds new light on previous work. In the work of Swenson and colleagues [12], authors tested 2 selection regimens with Newborn sizes differing by 100-fold. The authors hypothesized that smaller Newborns would have a high level of variation that should facilitate selection. However, the hypothesis was not corroborated by experiments. As a possible explanation, the authors invoked the "butterfly effect" (the sensitivity of chaotic systems to initial conditions). Our results suggest that even for nonchaotic systems like the H-M community, selection could fail due to interference from nonheritable variations. This is because in Newborns with small sizes, fluctuations in community composition can be large, which compromises heredity of a community trait.
A general implication of our work is that before launching a selection experiment, one should carefully design the selection regimen. For example, one may want to check the sensitivity of community function to fluctuations in Newborn biomass composition. The first method of checking involves estimating the "signal to noise" ratio: one could use the most precise method to initiate Newborn community replicates and measure community functions (e.g., cell sorting during Newborn formation; many repeated measurements of community function). Despite this, some levels of nonheritable variations in community function are inevitable due to, e.g., nongenetic phenotypic variations among cells [74] or stochasticity in cell birth and death. If "noises" (variations among replicate communities) are small compared to "signals" (variations among communities of different genotypes that affect the community function of interest), then one can test and possibly adopt less precise procedures (e.g., cell culture pipetting during Newborn formation; fewer repeated measurements of community function). The second checking method involves empirically estimating the heritability of community function, especially if significant variations in community function naturally arise within the first few cycles (due to, e.g., preexisting mutations). In this case, one could experimentally evaluate whether community functions of the previous cycle are correlated with community functions of the current cycle across independent lineages (similar to Fig 4). Given the ubiquitous nature of nonheritable variations in community function, the top-tier strategy could be useful (Figs 5 and 6).
Our work also contributes to how we might think about community selection in the natural environment. Microbes can co-evolve with each other and with their host [75,76,77]. Some have proposed that complex microbial communities such as the gut microbiota could serve as a unit of selection [19]. Our work suggests that if selection for a costly microbial community function should occur in nature, then certain mechanisms may need to be in place. These mechanisms include (i) suppressing nonheritable variations in community function, and (ii) exerting an appropriate strength of intercommunity selection.

Future directions
Our work touched upon only the tip of the iceberg of community selection. We expect that certain rules will be insensitive to details of a community. For example, community selection can be facilitated by reducing nonheritable variations in community function, or by mitigating the effects of nonheritable variations via top-tier strategies. Still, much more awaits further investigation. Here, we outline a few: 1. Explore the best strategies for choosing and reproducing Adult communities. We have chosen top 10% communities to contribute an equal number of Newborns, but alternative strategies (e.g., allowing higher-functioning Adults to contribute more Newborns) may work better.
2. Examine the impact of migration (community mixing) on community selection. Here, we did not allow migration. Excessive migration could deter community selection by allowing fast-growing noncontributors to spread. However, by combining the best genotypes of multiple member species, migration could speed up community selection, much like the effects of sexual recombination on the evolution of finite populations [78].
3. Investigate how interaction structure might affect selection efficacy. We have shown that our conclusions hold for 2-species communities engaging in commensalism or mutualism.
We have also shown that our conclusions hold regardless of whether the 2 species can evolve to coexist or not. The next step would be to test other types of interactions and complex interaction networks. For example, when species mutually inhibit each other, multistability could arise. In this case, species dominance [79] and thus community function could be highly sensitive to stochastic fluctuations in Newborn species composition. How might multistability affect community selection [117,118]?
A complex community might also have multiple local optima of community function. This can arise when, e.g., the maximal community function requires multiple species with partially redundant activities. Consider the community function of waste degradation. Suppose that one species is good at degrading waste at high concentration while the other is good at degrading waste at low concentration. Then, maximal waste degradation requires both species. If any one species is lost by chance during Newborn formation, then community function would be stuck at local suboptima. In this case, community migration (mixing) could recover the lost species and help community selection reach a global optimum.
4. Develop a general theory to understand the rate of community function improvement. Community function improvement depends on many experimental parameters, as we have demonstrated here. Ultimately, the rate of improvement will depend on variation and heredity of community function, which are affected by intracommunity and intercommunity selection. How might experimental parameters affect variation and heredity of community function under selection? This aspect might be explored through expanding population genetics theories which have so far focused on individual selection [114].
5. Experimentally test model predictions. Effective community selection will require a fast and precise assay for community function. If community function is sensitive to species biomass composition in Newborns, then species should ideally be distinguishable by flow cytometry (e.g., different fluorescence or different scattering patterns) so that fixed biomass of each species can be sorted into Newborns. Note that cell sorting only needs to be performed on several high-functioning communities and thus would not be cost prohibitive if the total number of Newborn communities is moderate.
6. Discover drivers of community function. Once high-functioning communities are obtained through selection, one could compare metagenomes of evolved communities with those of ancestral communities. This could illuminate genes and species interactions that are important for community function.

Equations of H-M community
H, the biomass of H, changes as a function of growth and death, Growth rate g H depends on the level of ResourceR (hat^representing prescaled absolute value) as described by the Monod growth model whereK HR is theR at which g Hmax /2 is achieved. δ H is the death rate of H. Note that because agricultural waste is in excess, its level does not enter the equation. M, the biomass of M, changes as a function of growth and death, Total potential growth rate of M g M depends on the levels of Resource and Byproduct (R and B) according to the Mankad-Bungay model [30] which has received experimental support (S4 Fig) Fig). Out of total potential growth rate of M, 1 − f P fraction is channeled to biomass increase, while f P fraction is channeled to making Product dP dt wherer P is the amount of Product made at the cost of one M biomass (tilde "~" represents scaling factor; see below and Table 1). ResourceR is consumed proportionally to the growth of M and H; ByproductB is released proportionally 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.ĉ BM is the amount ofB consumed per potential M biomass.r B is the amount ofB released per H biomass grown. Our model assumes that Byproduct or Product is generated proportionally to H or M biomass grown, which is reasonable given the stoichiometry of metabolic reactions and experimental support [80]. The volume of community is set to 1, and thus cell or metabolite quantities (which are considered here) are numerically identical to cell or metabolite concentrations. In the equations above, scaling factors are marked by "~" and will become 1 after scaling. Variables and parameters with hats will be scaled and lose their hats afterwards. Variables and parameters without hats will not be scaled. We scale Resource-related variable (R) and parameters (K MR ;K HR ;ĉ RM , andĉ RH ) againstRð0Þ (Resource supplied to Newborn), Byproductrelated variable (B) and parameters (K MB andĉ BM ) againstr B (amount of Byproduct released per H biomass grown), and Product-related variable (P) againstr P (amount of Product made at the cost of one M biomass). For biologists who usually think of quantities with units, the purpose of scaling (and getting rid of units) is to reduce the number of parameters. For example, H biomass growth rate can be rewritten 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 identical forms (S3 Fig). After scaling, the value ofRð0Þ becomes irrelevant (1 with no unit).
Thus, scaled equations are as follows: We have not scaled time here, although time can also be scaled by, e.g., the community maturation time. Here, time has the unit of unit time (e.g., hour), and to avoid repetition, we often drop the time unit. Scaling factors and values of species phenotypes after scaling are in Table 1. Additional symbols, including state variables and selection scheme parameters, are summarized in Table 2.
Below, we calcuate the steady-state species ratio in the H-M community. From Eq 10: If we approximate Eqs 6-7 by ignoring the death rates so that dH dt � g H R ð ÞH and In our simulations, with our parameter choices, if f P is not too large (f P < 0.
and the M:H ratio at time T is

MðTÞ
HðTÞ Thus when Byproduct B(T) � 0, ϕ M,SS , the steady-state fraction of M biomass, is then For f P < 0.4, Eq 14 is applicable and predicts the steady-state ϕ M,SS well (see S26 Fig). Note that significant deviation occurs when f P > 0.4. This is because when f P is large, M's biomass does not grow fast enough to deplete B so that we cannot approximate B(T) � 0 anymore.

Parameter choices
Our parameter choices are based on experimental measurements from a variety of organisms. Additionally, we chose growth parameters (maximal growth rates and affinities for metabolites) of ancestral and evolved H and M so that (i) the 2 species can coexist at a moderate ratio for a range of f P over selection cycles and (ii) improving all growth parameters up to their evolutionary upper bounds generally improves community function (Methods Section 3). This way, we could simplify our simulation by fixing growth parameters at their respective evolutionary upper bounds. With only 1 mutable parameter (f P ), we can identify the optimal f � P associated with maximal community function (Fig 2B).
To ensure the coexistence of H and M, M must grow faster than H for part of the maturation cycle since M has to wait for H's Byproduct at the beginning of a cycle. Because we have assumed M and H to have similar affinities for R (Table 1), g Mmax must exceed g Hmax , and M's affinity for Byproduct (1/K MB ) must be sufficiently large. Moreover, metabolite release and consumption need to be balanced to avoid extreme species ratios. Thus, for ancestral M, we chose g Mmax = 0.58 (equivalent to a doubling time of 1.2 hours). We set c BM ¼ 1 3 (units ofr B ), meaning that Byproduct released during one H biomass growth is sufficient to generate 3 potential M biomass, which is biologically achievable [40,81]. When we chose K MB ¼ 5 3 � 10 2 (units ofr B ), H and M can coexist for a range of f P (Fig 2A). This value is biologically realistic. For example, suppose that H releases hypoxanthine as Byproduct. A hypoxanthine-requiring S. cerevisiae strain evolved under hypoxanthine limitation could achieve a Monod constant for hypoxanthine on the order of 0.1 μM [64]. If the volume of the community is 10 μL, then K MB ¼ 5 3 � 10 2 (units ofr B ) corresponds to an absolute release rater B ¼ 0:1 mM � 10 mL= 5 3 � 10 2 À � ¼ 6 fmole per releaser biomass born. At 8-hour doubling time, this translates to 6 fmole/(1 cell×8 h)� 0.75 fmole/cell/h, within the ballpark of experimental observation (approximately 0.3 fmole/cell/h [64]). As a comparison, a lysine-overproducing yeast strain reaches a release rate of 0.8 fmole/cell/h [64], and a leucine-overproducing strain reaches a release rate of 4.2 fmole/cell/ h [81]. Death rates δ H and δ M were chosen to be 0.5% of H and M's respective upper bound of maximal growth rate, which are within the ballpark of experimental observations (e.g., the death rate of a lys-strain in lysine-limited chemostat is 0.4% of maximal growth rate [64]).
We assume that H and M consume the same amount of Resource R per new cell (c RH = c RM ) because the biomass of various microbes share similar elemental (e.g., carbon or nitrogen) compositions [82]. Specifically, c RH = c RM = 10 −4 (units ofRð0Þ), meaning that the Resource supplied to each Newborn community can yield a maximum of 10 4 total biomass.
In some simulations (e.g., Fig 6, S8 and S17 Figs), growth parameters (maximal growth rates g Mmax and g Hmax and affinities for nutrients 1/K MR , 1/K MB , and 1/K HR ) and production cost parameter (0 � f P � 1) were allowed to change from ancestral values during community maturation because these phenotypes have been observed to rapidly evolve within dozens to hundreds of generations [32,33,34,35]. For example, several-fold improvement in nutrient affinity and 20% increase in maximal growth rate have been observed in experimental evolution [33,35]. We therefore allowed affinities 1/K MR , 1/K HR , and 1/K MB to increase by up to 3-fold, 5-fold, and 5-fold, respectively, and allowed g Hmax and g Mmax to increase by up to 20%. These evolutionary upper bounds also ensured that evolved H and M could coexist for f P < 0.5 and that Resource was, on average, not depleted by T to avoid stationary phase.
We also simulated community selection in which improved growth parameters could reduce community function (Fig 6, S17 and S24B Figs). In this simulation, g Hmax was allowed to increase by up to 220%, and each Newborn community was supplied with Resource that can support up to 10 5 cells (10 units ofRð0Þ).
Although empirical studies sometimes find trade-off between maximal growth rate and nutrient affinity (e.g., [33]), for simplicity we assumed here that the two traits are independent of each other. We held metabolite consumption (c RM , c BM , c RH ) constant because conversion of essential elements such as carbon and nitrogen into biomass is unlikely to evolve quickly and dramatically, especially when these elements are not in large excess [82]. Similarly, we held the scaling factorsr P (Product released at the cost of one M biomass) andr B (Byproduct released per H biomass grown) constant, assuming that they do not change rapidly during evolution. Indeed, metabolite release rate and metabolite consumption amount per biomass (biomass measured as cell size) did not change significantly for a commonly-arising mutant in a yeast evolution experiment [83]. We held death rates (δ M , δ H ) constant because they are much smaller than growth rates in general and thus any changes are likely inconsequential.

Choosing growth parameter ranges so that we can fix growth parameters to upper bounds
Improving growth parameters (maximal growth rate and affinity for metabolites) does not always lead to improved community function (Fig 6, S17 and S24B Figs). However, for most simulations, we have chosen H and M growth parameters so that improving them from their ancestral values up to evolutionary upper bounds generally improves community function (see below). H and M with growth parameters at upper bounds are called "preadapted". When Newborn communities are assembled from preadapted H and M, 2 advantages are apparent.
First, after fixing growth parameters of H and M to their upper bounds, we can identify the locally maximal community function. Specifically, for a Newborn with total biomass BM(0) = 100 and fixed Resource R, we can calculate community function P(T) under various cost f P and fraction M biomass in Newborn ϕ M (0), assuming that all M cells have the same f P . Since both numbers range between 0 and 1, we calculate P(T, f P = 0.01 × i, ϕ M (0) = 0.01 × j) for integers i and j between 1 and 99. There is a single maximum for P(T) when i = 41 and j = 54. In other words, if M invests f � P ¼ 0:41 of its potential growth to make Product and if the fraction of M biomass in Newborn � � M ð0Þ ¼ 0:54, then maximal community function P � (T) is achieved (Fig 2B; magenta dashed line in Fig 3).
Second, preadapted H and M are evolutionarily stable in the sense that deviations (reductions) from upper bounds will reduce both individual fitness and community function (S27 Fig), and are therefore disfavored by intracommunity selection and intercommunity selection.
Below, we present evidence that within our parameter ranges (Table 1), improving growth parameters generally improves community function. When f P is optimal for community function (f � and thus @g M @K MR < 0. This is the familiar case in which growth rate increases as the Monod constant decreases (i.e., affinity increases). However, if and thus @g M @K MR > 0. In this case, growth rate decreases as the Monod constant decreases (i.e., affinity increases). In other words, decreased affinity for the abundant nutrient improves growth rate. Transporter competition for membrane space [84] could lead to this result, since reduced affinity for abundant nutrient may increase affinity for rare nutrient.
At the beginning of each cycle, R is abundant and B is limiting (Eq 16). Therefore M cells with lower affinity for R will grow faster than those with higher affinity (

Mutation rate and the distribution of mutation effects
Literature values of mutation rate and the distribution of mutation effects are highly variable. Below, we briefly review the literature and discuss rationales of our choices.
Among mutations, a fraction is neutral in that they do not affect the phenotype of interest. For example, the vast majority of synonymous mutations are neutral [85]. Furthermore, mutations with small effects may appear neutral, which can depend on the effective population size and selection condition. For example, at low population size due to genetic drift (i.e., changes in allele frequencies due to chance), a beneficial or deleterious mutation may not be selected for or selected against and is thus neutral with respect to selection [86,87]. As another example, the same mutation in an antibiotic-degrading gene can be neutral under low antibiotic concentrations but deleterious under high antibiotic concentrations [88]. We term all mutations with zero effects on phenotypes as "neutral" mutations.
Since a larger fraction of neutral mutations is equivalent to a lower rate of phenotype-altering mutations, our simulations define "mutation rate" as the rate of non-neutral mutations that either enhance a phenotype ("enhancing mutations") or diminish a phenotype ("diminishing mutations"). Enhancing mutations of maximal growth rates (g Hmax and g Mmax ) and of nutrient affinities (1/K HR , 1/K MR , 1/K MB ) enhance the fitness of an individual ("beneficial mutations"). In contrast, enhancing mutations in f P diminish the fitness of an individual ("deleterious mutations").
Depending on the phenotype, the rate of phenotype-altering mutations is highly variable. Mutations that cause qualitative phenotypic changes (e.g., drug resistance) occur at a rate of 10 −8~1 0 −6 per genome per generation in bacteria and yeast [89,90]. In contrast, mutations affecting quantitative traits such as growth rate occur much more frequently. For example in yeast, mutations that increase growth rate by � 2% occur at a rate of * 10 −4 per genome per generation (calculated from Fig 3 of Ref. [91]), and mutations that reduce growth rate occur at a rate of 10 −4 * 10 −3 per genome per generation [38,92]. Moreover, mutation rate can be elevated by as much as 100-fold in hypermutators where DNA repair is dysfunctional [92,93,94]. In our simulations, we assume a high, but biologically feasible, rate of 2 × 10 −3 phenotypealtering mutations per cell per generation per phenotype to speed up computation. At this rate, an average community would sample approximately 20 new mutations per phenotype during maturation. We have also simulated with a 100-fold lower mutation rate. As expected, evolutionary dynamics slowed down, but all of our conclusions still held (S18 Fig). Among phenotype-altering mutations, around 10% to 40% create null (loss-of-function) mutants, as illustrated by experimental studies on protein, viruses, and yeast [36,37,38]. Thus, we assumed that 50% of phenotype-altering mutations were null (i.e., resulting in zero maximal growth rate, zero affinity for metabolite, or zero f P ). Among non-null mutations, the relative abundances of enhancing versus diminishing mutations are highly variable in different experiments. It can be impacted by effective population size. For example, with a large effective population size, the survival rate of beneficial mutations is 1,000-fold lower due to clonal interference (competition between beneficial mutations) [95]. The relative abundance of enhancing versus diminishing mutations also strongly depends on the starting phenotype [36,86,88]. For example, with ampicillin as a substrate, the wild-type TEM-1 β-lactamase is a "perfect" enzyme. Consequently, mutations were either neutral or diminishing, and few enhanced enzyme activity [88]. In contrast, with a novel substrate such as cefotaxime, the enzyme had undetectable activity, and diminishing mutations were not detected, whereas 2% of tested mutations were enhancing [88]. When modeling H-M communities, we assumed that the ancestral H and M had intermediate phenotypes that can be enhanced or diminished.
We based our distribution of mutation effects on experimental studies in which a large number of enhancing and diminishing mutants have been quantified in an unbiased fashion. An example is a study from the Dunham lab in which the fitness effects of thousands of S. cerevisiae mutations were quantified under various nutrient limitations [39]. Specifically for each nutrient limitation, the authors first measured Δs WT ¼ ðw WT À � w WT Þ= � w WT ¼ w WT = � w WT À 1, the deviation in relative fitness of thousands of bar-coded wild-type control strains from the wild-type mean fitness � w WT (i.e., selection coefficients). Due to experimental noise, Δs WT is distributed with zero mean and nonzero variance. Then, the authors measured thousands of Δs MT , each corresponding to the relative fitness change of a bar-coded mutant strain with respect to the mean of wild-type fitness (i.e., Δs MT ¼ ðw MT À � w WT Þ= � w WT ). From these 2 distributions, we used convolution to derive μ ΔS , the probability density function (PDF) of relative fitness change caused by mutations Δs = Δs MT − Δs WT (S6 Fig), in the following manner.
First, we calculated μ m (Δs MT ), the discrete PDF of the relative fitness change of mutant strains, with bin width 0.04. In other words, μ m (Δs MT ) = counts in the bin of [Δs MT − 0.02, Δs MT + 0.02] � total counts � 0.04 where Δs MT ranges from −0.6 and 0.6, which is sufficient to cover the range of experimental outcome. The Poissonian uncertainty of μ m is dm m ðΔs MT Þ ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi counts per bin p � total counts � 0.04. Repeating this process for the wild-type collection, we obtained the PDF of the relative fitness change of wild-type strains μ w (Δs WT ). Next, from μ w (Δs WT ) and μ m (Δs MT ), we derived μ Δs (Δs), the PDF of Δs with bin width 0.04: where μ w (j) is short-hand notation for μ w (Δs WT = j × 0.04) and so on. Our calculated μ Δs (Δs) with error bar of δμ Δs is shown in S6 Fig. Our re-analysis demonstrated that distributions of mutation fitness effects μ Δs (Δs) are largely conserved regardless of nutrient conditions and mutation types (S6 Fig). In all cases, the relative fitness changes caused by beneficial (fitness-enhancing) and deleterious (fitnessdiminishing) mutations can be approximated by a bilateral exponential distribution with means s + and s − for the positive and negative 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 μ Δs (i)/δμ Δs (i) as the weight for nonlinear least squared regression (green lines in S6 Fig). We obtained s + = 0.050 ± 0.002 and s − = 0.067 ± 0.003.
Exponential distribution described the fitness effects of deleterious mutations in an RNA virus remarkably well [36]. Based on extreme value theory, the fitness effects of beneficial mutations were predicted to follow an exponential distribution [96,97], which has gained experimental support from bacterium and virus [98,99,100] (although see [91,101] for counterexamples). Evolutionary models based on exponential distributions of fitness effects have shown good agreement with experimental data [95,102].
We have also simulated smaller average mutation effects based on measurements of spontaneous or chemically induced (instead of deletion) mutations. For example, the fitness effects of nonlethal deleterious mutations in S. cerevisiae were mostly between 1% to 5% [38], and the mean selection coefficient of beneficial mutations in E. coli was approximately 1% to 2% [95,98]. Thus, as an alternative, we simulated with s + = s − = 0.02. We obtained the same conclusions (S19A Fig).

Modeling epistasis on f P
Epistasis, in which the effect of a new mutation depends on prior mutations ("genetic background"), is known to affect evolutionary dynamics. Epistatic effects have been quantified in various ways. Experiments on viruses, bacteria, yeast, and proteins have demonstrated that if 2 mutations were both deleterious or both random, viable double mutants experienced epistatic effects that distributed nearly symmetrically around a value close to zero [103,104,105,106]. In other words, a significant fraction of mutation pairs show no epistasis, and a small fraction show positive or negative epistasis (i.e., a double mutant displays a stronger or weaker phenotype than expected from additive effects of the 2 mutations). Epistasis between 2 beneficial mutations can vary from being predominantly negative [104] to being symmetrically distributed around zero (Fig 2 of [107]). Furthermore, a beneficial mutation tends to confer a lower beneficial effect if the background already has high fitness ("diminishing returns") [107,108,109].
A mathematical model by Wiser and colleagues incorporates diminishing-return epistasis [102]. In this model, beneficial mutations of advantage s in the ancestral background are exponentially distributed with PDF αe −αs , where 1/α > 0 is the mean advantage. After a mutation with advantage s has occurred, the mean advantage of the next mutation would be reduced to 1/[α(1 + gs)], where g > 0 is the "diminishing returns parameter." Wiser and colleagues estimate g � 6. This model quantitatively fits the fitness dynamics of evolving E. coli populations.
Based on the above experimental and theoretical literature, we modeled epistasis on f P in the following manner. Let the relative mutation effect on current f P be Δf P = (f P,mut − f P )/f P (note Δf P � −1). 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: Here, s + (f P ) and s − (f P ) are, respectively, the mean Δf P for enhancing and diminishing mutations 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 initial background in a community selection simulation (f P;init ¼ f � P;Mono ¼ 0:13), s +init is the mean enhancing Δf P occurring in the initial background, and 0 < g < 1 is the epistatic factor. Similarly, s − (f P ) = s −init × (1 + g × (f P /f P,init − 1)) is the mean |Δf P | for diminishing mutations at current f P . In the initial background, because f P = f P,init , we have s + (f P ) = s +init and s − (f P ) = s −init (s +init = 0.050 and s −init = 0.067 in S6 Fig). Consistent with the diminishing returns principle, for subsequent mutations that alter f P , if current f P > f P,init , then a new enhancing mutation became less likely and its mean effect smaller, while a new diminishing mutation became more likely and its mean effect bigger (ensured by g > 0; S20 Fig right panel). Similarly, if current f P < f P,init , then a new enhancing mutation became more likely and its mean effect bigger, while a diminishing mutation became less likely and its mean effect smaller (ensured by 0 < g < 1; S20 Fig left panel). In summary, our model captured not only diminishing returns epistasis, but also our understanding of mutational effects on protein stability [86].

Simulation code of community selection
As described in the main text, our simulations tracked the biomass and phenotypes of individual cells as well as the amounts of Resource, Byproduct, and Product in each community throughout community selection. Cell biomass growth, cell division, and changes in chemical concentrations were calculated deterministically. Stochastic processes, including cell death, mutation, and the partitioning of cells of a chosen Adult community into Newborn communities were simulated using the Monte Carlo method.
Specifically, each simulation was initialized with a total of n tot = 100 Newborn communities with identical configuration: • Each community had 100 cells of biomass 1. Thus, total biomass BM(0) = 100. At the beginning of each selection cycle, a random number was used to seed the random number generator for each Newborn community. This number was saved so that the maturation of each Newborn community can be replayed. In most simulations, the initial amount of Resource was 1 unit ofRð0Þ unless otherwise specified, the initial Byproduct was B(0) = 0, and the initial Product was P(0) = 0.
The maturation time T was divided into time steps of Δτ = 0.05. Resource R(t) and Byproduct B(t) during each time interval [τ, τ + Δτ] were calculated by solving the following equations (similar to Eqs 9-10) using the initial condition R(τ) and B(τ) via the ode23s solver in Matlab (MathWorks; www.mathworks.com): Thus, Similarly, let the length of an M cell be L M (τ). The continuous growth of M could be described as Thus for an M cell, its length L M (τ + Δτ) could be described as From Eqs 7-8, within Δτ since death is modeled separately, we have where M(τ + Δτ) = SL M (τ + Δτ) represented the sum of the biomass (or lengths) of all M cells at τ + Δτ.
At the end of each Δτ, each H and M cell had a probability of δ H Δτ and δ M Δτ to die, respectively. This was simulated by assigning a random number between 0 and 1 for each cell. Cells assigned with a random number less than δ H Δτ or δ M Δτ then got eliminated. For surviving cells, if a cell's length � 2, this cell would divide into 2 cells with half the original length.
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, then f P would be multiplied by (1 + Δf P ), where Δf P was determined as below.
If a mutation increased or decreased the phenotypic parameter beyond its bound (Table 1), the phenotypic parameter was set to the bound value.
The above growth, death, division, and mutation cycle was repeated from time 0 to T. Note that because the size of each M and H cell is larger than or equal to 1, the integer numbers of In the top-dog strategy, we started with the Adult with the highest community function. We first calculated the fold by which this Adult would be diluted as n D ¼ b MðTÞþHðTÞ BM target c. Here, BM target = 100 was the preset target for Newborn total biomass, and bxc is the floor (round down) function that generates the largest integer that is smaller than x. Note that rounding down causes negligible errors. For example, with 7,000 total biomass in an Adult, fold dilution n D ¼ 7000 100 ¼ 70 and each Newborn would have an average 100 total biomass. In comparison, with 6999 total biomass in an Adult, n D = 69 and each Newborn would have an average 101.4 total biomass. This difference of 1% is much smaller than stochastic fluctuations caused by pippetting. If n D was less than n tot , the total number of communities within a cycle, then all n D Newborn communities were kept, and the Adult with the next highest function was partitioned to obtain an additional batch of Newborns. From the last batch of Newborns, we would randomly choose enough to obtain n tot Newborns. The next cycle then began.
During reproduction of the top-tier strategy, n chosen Adults with the highest functions would each contribute n tot /n chosen Newborns for the next cycle. We have always picked n chosen so that n tot /n chosen is an integer.
Before community reproduction, the current random number generator state was saved so that the random partitioning of Adult communities could be replayed.
To mimic partitioning Adult communities via pipetting into Newborn communities with an average total biomass of BM target , we first calculated the fold of dilution n D as described above.

If the Adult community had I H (T) H cells and I M (T) M cells, I H (T) + I M (T) random
integers between 1 and n D were uniformly generated so that each M and H cell was assigned a random integer between 1 and n D . All cells assigned with the same random integer were then assigned to the same Newborn, generating n D Newborn communities. This partition regimen can be experimentally implemented by pipetting 1/n D volume of an Adult community into a new well.
To fix BM(0) to BM target and ϕ M (0) to ϕ M (T) of the parent Adult (cell sorting), the code randomly assigned M cells from the chosen Adult until the total biomass of M came closest to BM target ϕ M (T) without exceeding it. H cells were assigned similarly. Because each M and H cells had a length between 1 and 2, the biomass of M could vary between BM target ϕ M (T)−2 and BM target ϕ M (T), and the biomass of H could vary between BM target (1 − ϕ M (T))−2 and BM target (1 − ϕ M (T)). Variations in BM(0) and ϕ M (0) were sufficiently small so that community selection worked (Fig 3D and 3E). We also simulated sorting cells such that H and M cell numbers (instead of biomass) were fixed in Newborns. Specifically, bBM target φ M (T)/1.5c M cells and bBM target (1 − φ M (T))/1.5c H cells were sorted into each Newborn community, where we assumed that the average biomass of a cell was 1.5, and φ M (T) = I M (T)/(I M (T) + I H (T)) was calculated from cell numbers in the parent Adult community. We obtained the same conclusion (S13C Fig). To fix Newborn total biomass BM(0) to the target total biomass BM target while allowing ϕ M (0) to fluctuate (S12 Fig first and third columns), H and M cells were randomly assigned to a Newborn community until BM(0) came closest to BM target without exceeding it (otherwise, P(T) might exceed the theoretical maximum). For example, suppose that a certain number of M and H cells had been sorted into a Newborn so that the total biomass was 98.6. If the next cell, either M or H, had a biomass of 1.3, this cell would go into the community so that the total biomass would be 98.6 + 1.3 = 99.9. However, if a cell of mass 1.6 happened to be picked, this cell would not go into this community so that this Newborn had a total biomass of 98.6 and the cell of mass 1.6 would go to the next Newborn. Thus, each Newborn might not have exactly the biomass of BM target , but rather between BM target −2 and BM target . Experimentally, total biomass can be determined from the optical density, or from the total fluorescence if cells are fluorescently labeled [61]. To fix the total cell number (instead of total biomass) in a Newborn, the code randomly assigned a total of bBM target /1.5c cells into each Newborn, assuming an average cell biomass of 1.5. We obtained the same conclusion, as shown in S13A Fig. To fix ϕ M (0) to ϕ M (T) of the chosen Adult community from the previous cycle while allowing BM(0) to fluctuate (S12 Fig second and fourth columns)

Problems associated with an alternative definition of community function and an alternative method of reproducing Adults
Here, we describe problems associated with an alternative definition of community function and an alternative method of community reproduction.
An alternative definition of community function is Product per M biomass in an Adult community: P(T)/M(T). To illustrate problems with this definition, let's calculate P(T)/M(T) assuming that cell death is negligible. From Eqs 7 and 8, where biomass growth rate g M is a function of B and R. Thus, T is long enough for cells to double at least 3 or 4 times). In this case, If we define community function as P T ð Þ=M T ð Þ � f P 1À f P , then higher community function requires higher f P 1À f P or higher f P . However, if we select for very high f P , then M can go extinct (Fig 2 and S2 Fig).
In our community selection scheme, the average total biomass of Newborn communities was set to a constant BM target . Alternatively, each Adult community can be partitioned into a constant number of Newborn communities (i.e., fixed-fold dilution). If Resource is not limiting, then there is no competition between H and M, and P(T) increases as M(0) and H(0) increase. Therefore, selection for higher P(T) results in selection for higher Newborn total biomass (instead of higher f P , S23 Fig). This will continue until Resource becomes limiting, and then communities will get into the stationary phase.

f � P is smaller for M group than for H-M community
For groups or communities with a certain R T g M dt, we can calculate f P optimal for community function from Eq 28 by setting If R T g M dt � 1, f P is very small, then the optimal f P for P(T) is R T g M dt is larger in monoculture than in community because M grows faster in monoculture than in community. This is because B is supplied in excess in monoculture, whereas in community, H-supplied Byproduct is initially limiting. Thus, according to Eq 29, f � P � 1= R T g M dt is smaller for monoculture than for community.

Stochastic fluctuations during community reproduction
where "Cov" means covariance (equal to zero) and "Var" means variance.

Mutualistic H-M community
In the mutualistic H-M community, Byproduct inhibits the growth of H. According to Luli and colleagues [110], the growth rate of E. coli decreases exponentially as the exogenously added acetate concentration increases. Thus, we only need to modify the growth of H by a factor of exp(−B/B 0 ) where B is the concentration of Byproduct and B 0 is the concentration of Byproduct at which H's growth rate is reduced by e −1~0 .37: The larger B 0 , the less inhibitory effect Byproduct has on H, and when B 0 ! +1, Byproduct does not inhibit the growth of H. For simulations in S22 Fig, we Supporting information S1 Fig. Artificial selection is more challenging for multispecies communities than for individuals or monospecies groups. Artificial selection can be applied to any population of entities [111]. An entity can be an individual (A), a monospecies group (B), or a multispecies 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 by the experimentalist 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 [112]. In all 3 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 [21,22,23,113], since a trait is largely heritable so long as mutation and recombination are sufficiently rare. (B, C) In group selection and community selection, if maturation time 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). In this case, variation can be defined as the dissimilarity in Newborn compositions within a selection cycle, while heredity can be defined as the similarity of compositions between Newborns connected through lineage across consecutive selection cycles (tubes with same-colored outlines in Fig 4A). (B) Artificial selection of monospecies groups has been successful [18,46,48]. Suppose cooperators but not cheaters pay a fitness cost to generate a product (dark shade). Artificial selection for groups producing high total product favors cooperator-dominated groups, although within a group, cheaters grow faster than cooperators. If a Newborn group has a large number of cells (top), all Newborns will harbor similar fractions of cheaters, and thus intergroup variation will be small [71]. During maturation, cheater frequency will increase, thereby diminishing heredity. Low variation and low heredity interfere with group selection. 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 of the group, 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 multispecies 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 species and genotype composition further reduce heredity. Suppose that cell growth rate depends on each substrate 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 2 limiting substrates in a multiplicative fashion. (B) In the model proposed by Mankad and Bungay, 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-like fashion. However, when  Δs)). We derived μ Δs (Δs) from the Dunham lab data [39] 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 nonlinear least squared fitting of data to Eq 19 using all 3 sets of data. Note that data with larger uncertainty are given less weight and thus deviate more from the fitting. For an exponentially distributed PDF p(x) = exp(−x/r)/r where x, r > 0, and the average of x is r. When plotted on a semi-log scale, we get a straight line with slope −1/r, which 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. The Matlab codes can be found in S7 Code. PDF, probability density function.  (Table 1). Communities were 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 (panel B), this increase in P(T) must be due to improved growth parameters (panel C). (D-F) Community reproduction via biomass sorting (i.e., fixed Newborn total biomass and species composition). Community function improved to a much higher level (panel D). In both strategies, the 5 growth parameters increased to their respective evolutionary upper bounds (green dashed lines). Magenta dashed lines: optimal f P for community function and maximal community function P(T) when all 5 growth parameters are fixed at their evolutionary upper bounds and ϕ M (0) is also optimal for P(T). Black, cyan, and gray curves show independent simulations. � PðTÞ is averaged across chosen Adults. � g Mmax ; � g Hmax , and � f P are obtained by averaging within each chosen Adult and then averaging across chosen Adults. K SpeciesMetabolite are averaged within each chosen Adult, then averaged across chosen Adults, and finally inverted to represent average affinity. Note different horizontal axis scales. The maximal growth rates (g Mmax and g Hmax ) have the unit of 1/time. Affinity for Resource (1/K MR , 1/K HR ) has the unit of 1=Rð0Þ, whereRð0Þ is the initial amount of Resource in Newborn. Affinity for Byproduct (1/K MB ) has the unit of 1=r B , wherer B is the amount of Byproduct released per H biomass produced. Product P has the unit ofr P , the amount of Product released at the cost of 1 M biomass. More details on parameters and variables can be found in Tables 1 and 2. The simulation codes can be found in S9 Code, and the data can be found in S5 Data. (TIF) S9 Fig. Community function declines to zero in the absence of intercommunity selection for higher community function. An Adult was randomly chosen to reproduce as many Newborns as possible, and a second Adult was randomly chosen to reproduce more Newborns until 100 Newborns were obtained. Natural selection favored zero f P (panel A). Consequently, P(T) decreased to zero (panel B). Here, reproduction was done by pipetting. Black, cyan, and gray curves are independent simulation trials. � PðTÞ was averaged across the randomly chosen Adults. � f P ðTÞ was obtained by first averaging among M within each randomly chosen Adult and then averaging across the chosen Adults. The simulation codes can be found in S10 Code, and the data can be found in S6 Data. (TIF) S10 Fig. P � (T) is a local optimum because it cannot be further improved by small changes. We started each Newborn community with total biomass BM(0) = 100, all 5 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 5 growth parameters and f P to mutate while applying community selection. To ensure effective community selection (Fig 3D-3F), the strategy of top-dog with cell sorting was implemented. We found that all 5 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.64, compare with ϕ M,SS represented by the blue dashed line in Fig 2A top 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Þ. Other legends are the same as S8 Fig. The simulation codes can be found in S11 Code, and the data can be found in S7 Data.  Fig 3G-3I, respectively). Black, cyan, and gray curves are 3 independent simulation trials. � PðTÞ was averaged 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. The simulation codes can be found in S2 Code, and the data can be found in S8 Data. (TIF) S13 Fig. Fixing H Fig 3G-3H, except that selection here lasted more cycles. Compared to Fig 3A-3C (top-dog, pipetting), the top 10% strategy was more effective. However, compared to Fig 3D-3F (top-dog, cell sorting), the top 10% strategy was less effective, even over 10 4 cycles. Black, cyan, and gray curves are 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. The simulation codes can be found in S2 Code, and the data can be found in S1 Data. (TIF) S15 Fig. Top-tier strategies promoted community selection under a wide range of selection strengths. In a top-tier strategy ("top 2%" to "top 50%"), top n chosen (= 2*50) Adults each contributed 100/n chosen Newborns into the next cycle. Here, Adults were reproduced (split) into Newborns as if via pipetting. Note that "top 2%" yielded qualitatively similar results as "topdog". Note also that when all Adults contributed one Newborn each ("all 100%"), intercommunity selection strength was zero, and thus natural selection quickly reduced average cost f P and community function to zero. Black, cyan, and gray curves are 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. The simulation codes can be found in S2 Code, and the data can be found in S10 Data. (TIF) S16 Fig. The top-dog strategy is superior to the top 10% strategy when nonheritable variation in community function is low. Twenty replicas of selection simulations were performed using either the top-dog strategy (black curves) or the top-tier strategy (top 10 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 top-tier strategy. The simulation codes can be found in S2 Code, and the data can be found in S11 Data. (TIF) Fig 6, the evolutionary upper bound for g Hmax (g � Hmax ¼ 0:8) was larger than that of g Mmax (g � Mmax ¼ 0:7), opposite to that in Fig 3. (A) When using the top-dog strategy with pipetting, g Hmax and g Mmax evolved to their respective upper bounds, and thus g Hmax > g Hmax (compare i and iv). This would ordinarily lead to extinction of M. However, community selection managed to maintain M at a very low level (Fig 6A bottom panel). (B-D) When using the topdog strategy with cell sorting (panel B), the top 10% strategy with pipetting (panel C), or the top 10% strategy with cell sorting (panel D), community selection worked in the sense that both � f P and P(T) improved over cycles (Fig 6B-6D). 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 S8 Fig PðTÞ were more stable compared to when null mutations were present (Fig 3). Black, cyan, and gray curves are 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. For panel A, the simulation codes can be found in S14 Code, and the data can be found in S13 Data. For panel B, the simulation codes can be found in S15 Code, and the data can be found in S14 Data. , 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 (= 1/slope) becomes larger. If current f P is high (right), the opposite is true. (TIF) S21 Fig. 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 (Fig 3). Black, cyan, and gray curves are 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. For panel A, the simulation codes can be found in S16 Code, and the data can be found in S15 Data. For panel B, the simulation codes can be found in S16 Code, and the data can be found in S16 Data. 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. For panel A, the Matlab codes can be found in S17 Code. For panel B, the simulation codes can be found in S18 Code, and the data can be found in S17 Data. (TIF)

S23 Fig. Artificial selection in excess
Resource failed under fixed-fold pipetting dilution scheme. Excess Resource was supplied to each Newborn (R(0)/K MR = 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 (Fig 4F). Community function initially increased as Newborn total biomass increased exponentially (middle and bottom panels), while nonproducing M cells with f P = 0 quickly took over (top panel; S7B Fig). Black, cyan, and gray curves are 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. The simulation codes can be found in S19 Code, and the data can be found in S18 Data. (TIF) S24 Fig. Improving growth parameters can improve or impair community function, depending on evolutionary upper bounds of growth parameters. Plotted here are plateaued community function after 1,500 cycles when simulation did or did not allow mutations in growth parameters or f P . The top-dog strategy and pipetting were used. (A) When evolutionary upper bound for g Hmax (g � Hmax ¼ 0:3) was lower than that of g Mmax (g � Mmax ¼ 0:7), improving growth parameters improved community function. Compared to community function where no mutations were allowed (i), community function improved when both 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 S8A-S8C Fig. (B) When evolutionary upper bound for g Hmax (g � Hmax ¼ 0:8) was larger than that of g Hmax (g � Mmax ¼ 0:7), improving growth parameters could decrease community function. Compared to community function in which no mutations were allowed (i), community function decreased when both 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 S17A Fig and  Fig 6A. In panel B, Resource supplied to Newborn communities could support 10 5 total biomass to accommodate faster growth rate. Error bars are calculated form 3 independent selections. The simulation codes can be found in S2 and S4 Codes. The plot can be generated by S20 Code from S19 Data. . Then, maximal group function is achieved at f P ¼ f � P;Mono ¼ 0:13 (middle panel), lower than the optimal f P for the community function f � P ¼ 0:41 ( Fig 2B). Here, the growth parameters of M are all fixed at their evolutionary upper bounds, and P(T) has the unit ofr P . For panel A, the simulation codes can be found in S21 Code, and the data can be found in S20 Data. For panel B, the Matlab codes can be found in S22 Code.  Table 1, 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; Fig 2B) and f P ¼ f � P;Mono ¼ 0:13 (the starting point for most of our simulations; optimal for M monoculture production when Byproduct is in excess-see S25B Fig). Except for the growth parameter indicated on the horizontal axis, all other growth parameters were set to their respective upper bounds. (A-D) Community function increases as the indicated growth parameter increases. For example, in (A), all growth parameters except for g Mmax were set to their upper bounds. For each g Mmax , 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 preadapted population (with all growth parameters at upper bounds) over maturation time T = 17. The decreasing ratio indicates that the mutant has a lower fitness compared to the growth-adapted cells. doublings and therefore multiplied into approximately 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 approximately 23 doublings, the M monocultures have on average about 3% null mutants. Sixty 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. (Top panels) Histograms of the number of Newborn communities of the first cycle that are free of noncontributor M mutants when inoculated from a single M monoculture (Left panel) or from independently grown M monocultures (right panel). To generate the histograms, the pregrowth and inoculation process was repeated 1,000 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). The top-dog strategy was used to choose Adults that were then reproduced via cell sorting. Black, cyan, and gray curves are independent simulation trials. � PðTÞ was averaged across the 2 chosen Adults. � f P ðTÞ was obtained by first averaging among M within each chosen Adult and then averaging across the chosen Adults. The simulation codes for evolution dynamics can be found in S26 Code, the simulation codes for the histograms can be found in S27 Code, and the data can be found in S22 Data. (TIF) S1 Code. Matlab codes that generated Fig 2. (ZIP) Resources: Li Xie, Alex E. Yuan.