Ecoevolutionary dynamics of nested Darwinian populations and the emergence of communitylevel heredity
Abstract
Interactions among microbial cells can generate new chemistries and functions, but exploitation requires establishment of communities that reliably recapitulate communitylevel phenotypes. Using mechanistic mathematical models, we show how simple manipulations to population structure can exogenously impose Darwinianlike properties on communities. Such scaffolding causes communities to participate directly in the process of evolution by natural selection and drives the evolution of celllevel interactions to the point where, despite underlying stochasticity, derived communities give rise to offspring communities that faithfully reestablish parental phenotype. The mechanism is akin to a developmental process (developmental correction) that arises from densitydependent interactions among cells. Knowledge of ecological factors affecting evolution of developmental correction has implications for understanding the evolutionary origin of major egalitarian transitions, symbioses, and for topdown engineering of microbial communities.
Introduction
Thirty years ago, in an article arguing the importance of the ‘superorganism’, Wilson and Sober expressed surprise that biologists had not recognised that communities — in the laboratory — ‘could be treated as entities with heritable variation and selected accordingly’ (Wilson and Sober, 1989). That they might be treated as such, stemmed from recognition that the eukaryotic cell is a tightknit community of two once freeliving microbes (Margulis, 1970), but also from observations in nature of social insect colonies (Wilson, 1985), cellular slime molds (Bonner, 1982; Buss, 1982), and especially of phoretic insect communities (Wilson and Knollenberg, 1987).
Phoretic insect communities comprise a focal organism — often an insect such as a beetle — that moves between patchily distributed ephemeral resources carrying with it a myriad of associated organisms, including mites, nematodes and microbes. Communities associated with each insect differ by virtue of the composite members, with the conceivable possibility that some associations may harm the carrier insect, while others may bring benefit. Given that the role of dispersal is loosely analogous to a communitylevel reproduction event, Wilson and Sober argued that selection at the level of insect communities was likely to trump withincommunity selection leading to communities ‘becoming organised into an elaborate mutualistic network that protects the insect from its natural enemies, gathers food, and so on’.
If this might happen in nature, then why might this not be realised even more potently in the laboratory? Indeed, the logic of Darwinism says it should. Provided there exists heritable variance in fitness at the level of communities, then communities will participate as units (in their own right) in the process of evolution by natural selection (Lewontin, 1970; GodfreySmith, 2009). In nature, the necessary conditions are likely rare (Goodnight and Stevens, 1997), but ecological circumstances can sometimes conspire to ensure that variation among communities is discrete, that communities replicate and that offspring communities show some resemblance to parental communities. Phoretic insect communities are a plausible case in point. In the laboratory, however, the experimenter can readily construct conditions that ensure communities (or any collective of cells) are units of selection (Johnson and Boerlijst, 2002; Day et al., 2011; Xie and Shou, 2018; Xie et al., 2019). A critical requirement is a birthdeath process operating over a time scale longer than the doubling time of individual cells (Hammerschmidt et al., 2014; Rainey et al., 2017; Black et al., 2020).
Empirical support for the prediction that selection really can shape communities was provided by Swenson and colleagues who performed two studies in which artificial selection was imposed on microbial communities from soil (Swenson et al., 2000a; Swenson et al., 2000b). In the first, they selected communities that affected plant growth. In the second, they selected communities for ability to degrade the environmental pollutant 3chloroaniline. In both instances, communities at extreme values of community function were repeatedly propagated. In both studies, a significant response was measured at the level of the community.
Although the finding was a surprise (Goodnight, 2000), it is consistent with expectations that communities of entities — no matter their identity — will participate in the process of evolution by natural selection provided communities are discrete, they replicate, and that offspring communities resemble parental communities (GodfreySmith, 2009). Discreteness is conferred by simply compartmentalising communities via their placement in independent vessels. Replication is achieved by taking a sample of the selected communities with transfer to a new vessel. Heredity, however, is less tangible, especially in the Swenson experiments, where the selected communities were pooled before redistribution into fresh vessels. Nonetheless, intuition says that heredity becomes established through interactions (Wilson and Sober, 1989; Goodnight, 2000). Understanding the mechanistic bases of communitylevel heredity and its emergence motivates our study.
We begin by posing a thought experiment realisable via ever improving capacity to manipulate small volumes of liquid (Baraban et al., 2011; Sackmann et al., 2014; Cottinet et al., 2016). Consider a millifluidic device that controls the composition of emulsions. Consider thousands of microlitresized droplets each harbouring communities comprised of two types of microbes that differ solely in the colour of a fluorescent protein: one type encodes a red fluorescent protein and the other a blue fluorescent protein. Interest is in the evolution of communities that are of the colour purple (an equal ratio of redtoblue cells). Within each droplet, red and blue microbes replicate with growth rate and interaction rates being subject to evolutionary change. In the mean time, the experimenter, via lasers installed on the device, has determined the precise colour of each droplet and a priori decided that half of the droplets with composition furthest from an equal ratio of redtoblue will be eliminated, whilst the fraction whose colour is closest to purple will be allowed to replicate. Replication involves a dilution step during which an aliquot of cells are sampled and nutrients replenished. A further round of growth then ensues along with a further round of dropletlevel selection. The protocol continues thereafter with selection taking place at the level of communities via a birthdeath process. In essence the schema, outlined in Figure 1 and inherent in the work of Swenson and colleagues, involves exogenous imposition of ecological conditions sufficient to cause droplets to function as units of selection. The concept, elaborated in detail elsewhere, is referred to as ‘ecological scaffolding’ (Black et al., 2020).
Under this scaffoldedregime, communities within droplets are endowed with Darwinianlike properties (we use this term to convey the fact that removal of the scaffold leads, at least initially, to complete loss of communitylevel individuality). Collectivelevel variation is discretised by virtue of the bounds provided by the immiscibility of oil and water (communities are thus confined to droplets). Additionally, the device ensures that droplets engage in a birthdeath process: droplets furthest from the collectivelevel trait are extinguished, whereas those closest to the colour purple are diluted and split, thus effecting collectivelevel reproduction. Not determined by the device however is the relationship between parent and offspring droplets. Because the trait of the parent community depends on properties of the cellular constituents, there is — in the absence of interactions between red and blue cells — little chance that purplecoloured communities will reliably give rise to purplecoloured offspring. This is in part due to the stochastic nature of the dilution phase (a droplet with an equal ratio of red to blue is unlikely to give rise to offspring droplets founded with the same equal ratio of types) but also to withindroplet selection favouring fast growing types. Purplecoloured droplets can be maintained, as envisioned by the ‘stochastic corrector’ model (Maynard Smith and Szathmary, 1995; Grey et al., 1995; Johnston and Jones, 2015), provided only those communities with the correct colour are propagated. However, withindroplet selection favours rapidly growing cells resulting in successive reduction of the number of viable droplets.
Here, we show that when cellular interactions are also allowed to evolve, selection imposed at the collective level, leads to evolution of a developmentallike process, which ensures that offspring communities resemble parental communities, irrespective of the initial phenotype at the moment of birth. We illustrate the evolutionary process by means of stochastic simulations for nested populations of cells (particles) and communities (collectives) undergoing a deathbirth process. In order to generalise our findings we derive a deterministic approximation, which we then use to show how selection on community phenotype drives the evolution of ecological interactions that are the basis of communitylevel heredity.
Results
A nested model of collective evolution
As described above and depicted in Figure 1, we consider a nested model in which particles are discretised into a population of collectives. Each collective is comprised of two kinds of selfreplicating particle (red and blue) that together determine collective colour. Colour is important because it is the phenotype upon which collectives succeed or fail. Collectives that are too far from an optimal colour face extinction, whereas those within acceptable bounds persist with the possibility of reproduction. Birthdeath at the level of collectives affects the ecoevolutionary dynamics of particles as particlelevel traits that give rise to unfit collectives are eliminated.
We firstly present numerical simulations where particles undergo a stochastic birthdeath process (Champagnat et al., 2006; Doebeli et al., 2017) and collectives are selected based on colour. The details of the model are described in Appendix 1, but we introduce the main assumptions and underlying principles here. Each particle of type $i\in \{0,1\}$ is characterised by four traits (hereafter particle traits): colour (${c}_{i}$, red or blue), net maximum growth rate ${r}_{i}$, and two competition parameters (${a}_{i}^{\text{intra}}$ and ${a}_{i}^{\text{inter}}$). At any particular instant particles either reproduce or die. Particles of type i reproduce with a constant birth rate ${r}_{i}$ and die as a consequence of competition. The rate of death is densitydependent such that each particle of type increases the death rate of $i$type particles by ${r}_{i}{a}_{j}^{\text{intra}}$ if they share the same colour (${c}_{j}={c}_{i}$), or by ${r}_{i}{a}_{j}^{\text{inter}}$ when colours are different (${c}_{j}\ne {c}_{i}$). All transition rates can be found in Appendix 1, paragraph 'particlelevel ecology'. Competition rates are referred to as ‘interaction’ traits or parameters. We expand on more general types of interaction — from exploitative to mutualistic — in the Discussion and in Appendix 3.
Mutations are introduced at the level of particles. Mutation affects either particle maximum growth rate ($r$) or the intercolour competition parameter (${a}^{\text{inter}}$) by a small random quantity. In the spirit of adaptive dynamics (Geritz et al., 1998), the particle type carrying the new set of traits is referred to as a mutant, and the existing type is designated the resident. Within every collective and at any time, there are four populations composed of resident and mutant types of the two colours. Mutations are assumed to be rare. In order to accelerate numerical simulations, one mutant individual is introduced every time one population of a given colour goes extinct in one of the collectives. The newly added type has the the same colour as the extinct type.
Collectives also undergo a birthdeath process. The number of collectives $D$ is constant and collective generations are discrete and nonoverlapping. Each collective generation begins at time $t=0$ with offspring collectives containing $B$ founding particles. Particles replicate, interact and evolve according to the particle traits. After duration $T$, collectives attain ‘adult’ stage, and a fixed proportion of collectives $\rho $ is marked for extinction. This allows the possibility of selection on collectives based on their properties (the collective phenotype), which is derived from the composing particles. Our focus is collective colour, which is defined as the proportion $\varphi $ of red particles.
Initially, collectives contain red and blue particles in uniformly distributed ratios. Collectives are subject to evolution under two contrasting regimes: one neutral and the other selective. Under the neutral regime, the pool of collectives marked for extinction is sampled at random, whereas under the selective regime, collectives marked for extinction are those whose adult colour departs most from an arbitrarily fixed optimal colour $\widehat{\varphi}$. Extinguished collectives are replaced by offspring from uniformly sampled extant collectives (Figure 1). All other collectives are replaced by their own offspring. Reproduction involves uniformly sampling $B$ particles from the parent collective. Particles from one collective never mix with particles from any other. This establishes an unambiguous parentoffspring relationship (De Monte and Rainey, 2014). The adult colour of offspring collectives depends on the founding frequencies of particles (whose variance is negatively related to bottleneck size $B$), and on ensuing particlelevel population dynamics.
Selection on collectives drives the evolution of particle traits
In the absence of collectivelevel selection (neutral regime), collectives converge to a monochromatic phenotype (Figure 2A). Once collectives are composed of either allred or allblue particles, the contrasting colour cannot be rescued (colour change by mutation or migration is not possible). The distribution of collective colour becomes biased toward fastergrowing particle types, with selection driving a gradual increase in particle growth rate (Figure 2C). The intercolour competition trait (Figure 2E) is irrelevant once collectives become monochromatic (evolution is then governed by drift).
The dynamic is very different once selection is imposed at the level of collectives. By rewarding collectives closest to the colour purple (a fixed $\widehat{\varphi}=0.5$ ratio of red to blue particles), it is possible to prevent fixation of either colour (Figure 2B). Starting, as above, from collectives containing red and blue particles in uniformly distributed ratios, mean collective colour shifts toward red. The time scale is, as in the neutral case, a consequence of the faster initial growth rate of red particles. For a few tens of generations, the population of collectives remains strongly biased towards red. The optimal phenotype is maintained by selection for the least worse collective colour precisely as envisaged by the stochastic corrector model (Maynard Smith and Szathmary, 1995; Appendix 1—figure 5). Subsequently, however, the trend reverses and mean collective colour progressively approaches purple. From generation 1000, variance in the distribution of colour decreases, as a consequence of improvement in the ability of purpleparent collectives to give rise to offspring collectives that at adult age resemble parental types. This is associated with escalating particle growth rate (Figure 2D) and a saturating increase in betweencolour competition (Figure 2F). The latter reflects directional selection that moves the average phenotype in the population of collectives towards the optimal colour $\widehat{\varphi}$ (reached by generation 7000).
By affecting particle traits, selection on colour also modifies dynamics within collectives. Figure 3 shows variation of colour within a single collective growth phase at generation 3 and generation 9000. Prior to selection shaping particle traits, both red and blue particle types follow approximately exponential growth (Figure 3C). The resulting adult collective colour is thus biased towards the fastergrowing red type (Figure 3A). In contrast, at generation 9000 (Figure 3B), both particle types reach a saturating steady state that ensures that adult colour is purple. Initial departures from a 1:1 ratio — caused by the stochasticity of collective reproduction and/or particle growth dynamics — are compensated for during the growth phase (Figure 3D). Compensation is a consequence of the evolution of intercolour competition traits (Figure 2F). Population expansion is in turn dependent upon earlier increases in particle growth rate (Figure 2D). Moreover, selection favours competition trait values for which blue types have no effect on red types: ${a}^{inter}$ of blue types is close to zero by generation 5000 (Figure 2F).
Key features of the evolutionary trajectory discussed so far are representative of replicate realisations of the stochastic individualbased nested model. In the section ‘Variability of the derived particle traits’ of Appendix 1 we show that in repeated simulations, the average adult collective colour always falls within a few percentage of the target colour. Moreover, collectivelevel selection applies more stringently to interaction parameters than to growth rates, with the latter showing a broader distribution of derived values. These conclusions also hold when the target colour $\widehat{\varphi}$ is different from 0.5, with broader variation from one realisation to the other and slower convergence to the target as colour ratios become more extreme (see Appendix 1, paragraph ‘Different target colours’).
Ability of offspring collectives to correct departures from the optimal colour during the course of growth is akin to a developmental, or canalising process: irrespective of the phenotype of the newborn (which will likely be different to that of the adult) the child — as it grows to adulthood — develops a phenotype that closely resembles that of the parent. Evidence of this apparent canalising process can be seen upon removal of collectivelevel selection (Figure 4). Collectives founded by particles with ancestral traits become composed of a single (red or blue) colour in less than 10 generations (Figure 4A). In contrast, derived collectives are comprised of particles whose traits ensure that collectives continue to express phenotypes narrowly distributed around the optimal (purple) phenotype (as long as there is no mutation [Figure 4B]). Even when mutation is allowed to drive within and betweencollective dynamics, stability of phenoytpe holds for more than 200 generations (Appendix 1—figure 6).
From particle ecology to collective phenotype
To understand the mechanistic basis of the canalising process, particle traits must be linked to the evolutionary emergence of collectivelevel inheritance, which we define as the capacity of collectives to reestablish the parental collective colour. Figure 5 shows the relationship between the initial colour of collectives at the moment of birth (the moment immediately following dilution, $t=0$ [the newborn colour]), and collective colour after a single particle growth cycle (the moment immediately preceding dilution, $t=T$ [the adult colour]). Figure 5A shows this relationship at generation 3 while Figure 5B shows this relationship at generation 9000.
At generation 3, the proportion of red particles increases (within a collective generation), irrespective of the initial proportion. This is because red particles grow faster than blue and the primary determinant of particle success is growth rate (interactions are negligible in exponential growth). Thus, the only way that purple collectives can be maintained is if the collective phenotype is sufficiently noisy, to ensure that some collectives happen by chance to be purple, due to, for example, stochastic effects at dilution. Even if offspring collectives do not resemble their parents, purple colour is maintained via strong purifying selection that purges collectives that are either too red or too blue. The mechanism is stochastic correction (Maynard Smith and Szathmary, 1995; Grey et al., 1995; Johnston and Jones, 2015).
This is in marked contrast to the situation at generation 9000. After a single growth cycle, the proportion of red particles increases when the initial proportion is below, and decreases when it is above, the optimal proportion 0.5. Thus, at generation 9000, irrespective of initial conditions, the adult colour of any given collective will be closer to $\widehat{\varphi}=0.5$ than it was on founding. Accordingly, extreme purifying selection is no longer required to maintain the parentoffspring relationship. Indeed, offspring collectives return to the parent phenotype even when the phenotype at birth departs significantly from the parent (adult) phenotype. ‘‘Correction’’ stems from the ecological dynamics of the particles and resembles a developmental process. Hereafter we refer to this correction process as the developmental corrector.
The relationship between newborn and adult colour of collectives shown in Figure 5 can be used to follow the fate of collectives over several cycles of growth and reproduction, provided the stochastic effects associated with the dilution phase are momentarily ignored. The iteration using particle trait values from generation 3 is shown by the dotted line in Figure 5A (the adult colour of a collective is the newborn colour for the next cycle, following a 'staircase' geometric procedure). Because red particles grow faster than blue, it takes just six collective generations for red particles to fix within collectives. Conversely, after particle trait evolution (Figure 5B), the same staircase approach applied to newborn collectives of any colour shows rapid convergence to the colour purple (0.5) irrespective of the starting point. The difference in the relationship between initial and final colour at generation 3 and 9000 is evidence of the emergence of a mechanism for developmental correction.
In order to systematically explore the possible newborntoadult colour map and to understand how it changes through the evolution of particle traits, we use a deterministic approximation (orange line in Figure 5). This approximation is denoted $G$ or growth function (Appendix 2, Definition 2) and stems from an ordinary differential equation model often referred to as the competitive LotkaVolterra system (Appendix 2, Equation 1). This model is the limit for vanishing noise of the stochastic particle ecology, and provides a good approximation of the simulations (Dotted lines in Figure 3C–D). The growth function $G$ captures the outcome of the ecological dynamics (i.e. the fraction of red particles) after founding populations are allowed to grow for a finite time interval $T$. We note similarity between the $G$ function and the recently proposed ‘communityfunction landscape’ (Xie and Shou, 2018). The shape of $G$ depends on the value of particle traits $\theta $ (growth rates ${r}_{0}$ and ${r}_{1}$, and competition parameters ${a}_{00}={a}_{\text{0}}^{\text{intra}}$, ${a}_{10}={a}_{\text{0}}^{\text{inter}}$, ${a}_{01}={a}_{\text{1}}^{\text{inter}}$, ${a}_{11}={a}_{\text{1}}^{\text{intra}}$), but also on the bottleneck size at dilution $B$ and the collective generation duration $T$. The fixed points of $G$ (i.e. $\varphi $ such that $G(\varphi )=\varphi $) are of particular interest: in the deterministic model, these represent colours that are left unchanged during a generation. Such a fixed point is stable if the colours of collectives starting in its neighbourhood all converge to it ($\varphi =1$ in Figure 5A, $\varphi =0.5$ in Figure 5B), and unstable otherwise ($\varphi =0$ in Figure 5A, $\varphi =0$ and $\varphi =1$ in Figure 5B).
Under collectivelevel selection for colour, $T$ and $B$ are constant and particle traits evolve so that $G$ eventually has a stable fixed point, corresponding to the target colour $\widehat{\varphi}$. Progressive change in shape of the $G$ function across collective generations in a simulated lineage (Figure 2B) is illustrated in Figure 6. Note that these changes are continuous: small mutations in particle traits reflect as small changes in the shape of $G$.
The evolutionary trajectory of Figure 2B can now be understood in terms of the progressive evolution of particle traits (see Appendix 2 for a detailed description). At the beginning, particles compete mostly during exponential phase, so that adult colour is biased towards the fastgrowing type. Initial improvement in transmission of colour from parent to offspring arises as exponential growth rates ${r}_{i}$ of the particles align. Correspondingly, the $G$ function approaches linearity. A successive increase in maximal growth rate separates particle and collective time scales, allowing particles to experience densitydependent interactions. Eventually, such interactions evolve towards a regime where the $G$ function is nonlinear and fluctuations are readily compensated, thus developmental correction ensures a reliable colour inheritance.
The $G$ function, which allows characterisation of particle ecology, can now be used as a guide to optimise the ‘life cycle’ of growth and dilution that acts as a scaffold for the evolutionary process. In a typical experiment of communitylevel evolution, collective generation duration $T$ and bottleneck size $B$ are fixed. Some choices of these collectivelevel parameters are however likely to lead to collective phenotypes that are so far from the optimum that collective lineages go extinct. For instance, if in the firstgeneration competitive exclusion occurs rapidly, then distinguishing collectives based on collective colour may be impossible. Intuition suggests that the closer the fixed point of the $G$ function is to the target colour, the more efficient collectivelevel selection will be, and the faster the evolutionary dynamic. It is thus possible to use the distance between the fixed point of $G$ and the target composition $\widehat{\varphi}$ as a proxy for the probability that collective lineages will go extinct before attaining the desired colour. Below, we examine how the position of the fixed point of $G$ changes as a function of collective generation duration $T$ and bottleneck size $B$.
Effect of collective generation duration and bottleneck size
The growth function $G$ is readily computed from the particle traits and collective parameters even though it has in general no analytic expression (but see Appendix 2 for limit cases of exponential and saturating particle growth). There are four possible qualitative shapes of $G$, that differ in the position and stability of the fixed points (illustrated in Appendix 2—figure 31 to 34).
The qualitative dependence of $G$ and its fixed points on collectivelevel parameters varies with the underpinning particle ecology, making it easier for some communities to be starting points for the successful evolution of inheritance. Particle traits can be classified in four broad classes, based of the nature of the corresponding ecological equilibrium. For each of these classes, and when red particles grow faster than blue ${r}_{0}>{r}_{1}$, the fixed points of $G$ are illustrated in Appendix 2—figure 3A to 3D as a function of the collectivelevel parameters $B$ and $T$. Figure 7 refers to the situation where intercolour interaction traits are smaller than intracolour interaction traits. Here, particle populations converge in the long term to a coexistence equilibrium, where collective colour is ${\varphi}^{*}=\frac{{a}_{11}{a}_{01}}{{a}_{11}{a}_{01}+{a}_{00}{a}_{10}}$ (in general, different from the optimum). This equilibrium can be approached within a single collective generation if $T$ and $B$ are large (top right corner). On the other hand, when $T$ and $B$ are small (red region), the only stable fixed point invovles collectives composed solely of fastgrowing particles. This corresponds to cases where individual and collective time scales (quantified by ${r}^{1}$ and $T$, respectively) are insufficiently separated, or newborn size is too small, so that particle demography is essentially exponential and interactions cannot provide sufficient correction. For rapid evolution of collective colour, the most favourable starting position is one where the fixed point is closest to the optimal colour (for $\widehat{\varphi}=0.5$ this occurs for intermediate collective generation duration and bottleneck size for the trait values in Figure 7). Knowledge of the exact values requires, however, some preliminary measure of the ecological dynamics. Even in the absence of such information, the diagram in Figure 7 can be used to optimise experimental design by revealing intrinsic tradeoffs. A decrease in generation time, necessary for practical reasons, may for instance, be compensated by an increase in bottleneck size, without affecting the average collective phenotype.
Even when collectivelevel parameters are optimised so that the attractor of the $G$ function is initially close to the target colour, collectivelevel selection will keep acting on the particle traits, and affect phenotypic variability within the population of collectives. As stability of the fixed point increases, so to does fidelity of phenotype transmission from parent to offspring collectives. Once collectivelevel processes are set as to minimise the probability of collective extinction, the main obstacles to evolving higher inheritance come from constrains acting on particle traits, which may limit the range of attainable $G$ functions. Tradeoffs on particle ecology may prevent the $G$ function to attain an internal fixed point. We discuss two examples on constrained evolution in the following paragraph.
Constrained trajectories
Thus far, we have considered evolution within a fourdimensional parameter space defined by maximum growth rates and intercolour competition parameters. In real systems, however, constrains and tradeoffs may limit the range of achievable variations in particle traits. For instance, even though faster growing particles will always experience positive selection, cell replication rate cannot increase boundlessly. Here, we consider two instances of constrained evolution, where only a subset of particle traits are allowed to mutate.
First, we consider the case where competition parameters are vanishingly small, so that $G$ has no internal fixed point. Under such conditions, particle growth rates evolve to be identical (Figure 8A). In the absence of interactions, this is the only available solution to maintain collectives with an equal number of red and blue type particles. Under these circumstances, $G$ converges to the identity function. In the deterministic approximation, collective composition remains constant in time, but stochastic fluctuations that cause colour to deviate from the optimum are amplified across collective generations. These deviations are nonetheless corrected in the collective population by propagating only those collectives whose colour is closest to the optimum. Such stochastic correction (Maynard Smith and Szathmary, 1995), however, has a high risk of failure if selection is strong and collective population size is small.
Second, we consider the case when mutations only affect the two intertype competition parameters, while growth rates are held constant (to sufficiently high values, so that particles experience densitydependent effects in the growth phase). The evolutionary trajectory can be visualised in the plane of competition parameters $({a}_{01},{a}_{10})$. Figure 8B shows the result of a stochastic simulation superimposed to the value of $G(0.5)$. Independent of the initial values of the interaction parameters, evolution draws the system to the manifold associated with the optimal proportion $\widehat{\varphi}$ (white dashed line). Evolution within this manifold is neutral in the deterministic approximation, but the presence of stochastic fluctuations drives further improvement of the fitness landscape. Correction is indeed more efficient and the distribution of collective phenotypic diversity narrower when the gradient of $G$ in the fixed point is smaller. The condition on particle traits for the latter to vanish only depends on the carrying capacities of the two particle types, and corresponds to the type with smallest carrying capacity having zero interaction rate (see Appendix 2). A similar outcome is observed when, along an evolutionary trajectory, growth rates no longer influence adult colour (Figure 2). Developmental correction thus selects for maximal asymmetry in interactions, whereby one particle type drives the ecological dynamics of the other type, but is itself only affected by its own type (this is fully elaborated in Appendix 2).
Discussion
In nature, communities rarely ever qualify as units of selection in the traditional sense (Lewontin, 1970; GodfreySmith, 2009), because communities in nature rarely manifest heritable variance in fitness. In the laboratory, however, experimenters can exogenously impose (scaffold) Darwinianlike properties on communities such that they have no choice, but to become units of selection (Wilson and Sober, 1989; Xie et al., 2019; Black et al., 2020). This typically involves placement of communities in some kind of container (pot, testtube, flask, droplet, etc.) so they are bounded and variation at the community level is thus discrete. Communities are then allowed time for individual members to replicate and interact. At the end of the 'growth' period, community function is assessed based on predetermined criteria. The experimenter then effects replication of successful communities while discarding those that underperform. Replication typically involves transferring a sample of individuals from a successful community to a new container replete with a fresh supply of nutrients.
Experimental and theoretical studies indicate that artificial selection on microbial communities results in rapid functional improvement (Swenson et al., 2000a; Swenson et al., 2000b; Goodnight, 2000; Wade, 2016; Xie et al., 2019). This is not unexpected given that experimental manipulations ensure that communities engage directly in the process of evolution by (artificial) selection as units in their own right. However, for such effects to manifest there must exist a mechanism of communitylevel inheritance.
Consideration of both the effectiveness of artificial selection and the problem of heredity has led to recognition that the answer likely lies in interactions (Wilson and Sober, 1989; Swenson et al., 2000a; Swenson et al., 2000b; Goodnight, 2000; Rainey et al., 2017). The intuition stems from the fact that in the absence of interactions, communities selected to reproduce because of their beneficial phenotype will likely fail to produce offspring communities with similar functionality. If so, then these communities will be eliminated at the next round. Consider, however, an optimal community in which interactions emerge among individuals that increase the chance that offspring communities resemble the parental type. Such an offspring community will then likely avoid extinction at the next round: selection at the level of communities is thus expected to favour the evolution of interactions because inheritance of phenotype is now the primary determinant of the success (at the community level). Indeed, simulations of multispecies assemblages have shown that evolution of interaction rates not only improves diversitydependent fitness, but also increases collective ‘heritability’, defined as the capacity of randomly seeded offspring communities to reach the same dynamical state as their parents (Ikegami and Hashimoto, 2002; Penn, 2003). Further studies have stressed the role of the extracellular environment and of specific interaction networks, pointing out that microscopic constrains can affect the capacity of communities to participate in evolutionary dynamics at the higher level (Williams and Lenton, 2007; Xie and Shou, 2018; Xie et al., 2019).
Here, inspired by advances in millifluidics, we have developed a minimal mechanistic model containing essential ingredients of multiscale evolution and withincommunity competition. We considered collectives composed of two types of particles (red and blue) that interact by densitydependent competition. By explicitly modelling demographic processes at two levels of organisation, we have obtained mechanistic understanding of how selection on collective character affects evolution of composing particle traits. Betweencollective selection fuels changes in particlelevel traits that feedback to affect collective phenotype. Selection at the level of communities thus drives the evolution of interactions among particles to the point where derived communities, despite stochastic effects associated with sampling at the moment of birth, give rise to offspring communities that reliably recapitulate the parental community phenotype. Such is the basis of communitylevel inheritance. Significantly, it has arisen from the simplest of ingredients and marks an important initial step in the endogenisation of Darwinian properties: properties externally imposed stand to become endogenous features of the evolving system (Black et al., 2020).
The mechanism by which particles interact to establish community phenotype is reminiscent of a development process. We have refered to this as the ‘developmental corrector’. In essence, it is akin to canalisation, a central feature of development in complex living systems (Buss, 1987), and the basis of inheritance (Griesemer, 2002). Developmental correction solves the problem of implementing specific protocols for mitigating nonheritable variations in community function (Xie et al., 2019).
Developmental correction can be viewed as an evolutionary refinement of the stochastic corrector mechanism (Maynard Smith and Szathmary, 1995; Grey et al., 1995; Johnston and Jones, 2015). Both stochastic and developmental correctors solve the problem of producing enough wellformed collectives at each successive generation to prevent communitylevel extinction. The stochastic corrector mechanism relies on a lowfidelity reproduction process coupled to high population sizes. Deviations from successful collective states are corrected by purging collectives that depart significantly from the optimal collective phenotype. However, in the absence of strong collectivelevel selection the optimal community phenotype is rapidly lost. In contrast, the developmental corrector mechanism ensures that the optimal community phenotype is maintained without need for hard selection. Regardless of perturbations introduced by demography or low initial particle number, most collectives reliably reach a successful adult state. In our simulations, we show that community phenotype is maintained even in the absence of communitylevel selection, although ultimately mutational processes affecting particle dynamics result in eventual loss of the developmental corrector mechanism.
An operationally relevant question concerns the conditions (the initial state of the population, the nature of the scaffold and of particlelevel interactions) for selection on a collective character to result in evolution of developmental correction. While we did not detail the probability of collective lineage extinction, it is possible that collectives become monochromatic before evolution has had time to act on particle traits. In such cases, which are more likely if particlelevel traits are far from the region of coexistence, and if timescales of particle and collective generations are not well separated, then collectivelevel evolution will grind to a halt. In all other cases, provided there are no other evolutionary constraints, selection will eventually lead the system toward regions of particle traitspace where the collective phenotype becomes reliably reestablished. The efficiency of this selective process and its transient behaviour depend on collectivelevel parameters that control growth and reproduction.
From our individualbased simulations and ensuing deterministic approximation, it is clear that once densitydependent interactions govern the adult state, then collectivelevel selection for colour is promptly effected. This happens provided the intracollective ecology lasts long enough for nonlinear effects to curb particle growth. When this is not the case, for example when the bottleneck at birth is small, or collectivelevel generation time is too short, evolution of developmental correction will be impeded. The latter favours rapidly growing particles (Abreu et al., 2019) and offers little possibility for the evolution of developmental correction. When the ecological attractor within collectives leads to the extinction of one of the two types, long collectivelevel generation times are incompatible with the maintenance of diversity (van Vliet and Doebeli, 2019). However, in our model, particlelevel evolution changes the nature of the attractor from extinction of one of two types to stable coexistence, and concomitantly particle and collective timescales become separated. Even before developmental correction becomes established, evolution can transiently rely on stochastic correction to ensure the maintenance of particle coexistence.
There are two aspects to heredity: resemblance — the extent to which reproduction and development maintain the average offspring phenotype — and fidelity (or determination) — a measure of phenotypic variance (Jacquard, 1983; Bourrat, 2017). In our model, resemblance is established once densitydependent interactions counter the bias toward fast replicating particles: when the $G$ function has an internal fixed point in $\widehat{\varphi}$, systematic drift of average collective colour is prevented. The increase in resemblance is associated with progressive divergence of particle and collective demographic time scales. As a consequence, the collective phenotype is placed under the control of particle traits rather than demographic stochasticity. On a longer time scale, fidelity improves by subsequent changes in interaction parameters under the constraint that they do not affect average adult colour. The variance of the phenotype around the optimum is reduced by increasing canalisation (flattening of the $G$ function). This is best achieved by a strong asymmetry in the competition traits, whereby one type has a logistic, uncoupled, dynamic, and the second type adjusts its growth to the former’s density. Interestingly, it is always the type with the lower carrying capacity, regardless of its relative growth rate, that acts as the driver.
The relationship between parameters on very long time scales, when the adult colour is essentially dependent on interaction rates, depends critically on the space of possible values particle traits can assume. For instance, our analysis took into account only competitive intractions between colours. Extension of the deterministic approximation to cases when ecological interactions between colours are exploitative or mutualistic indicates that selection for collective coulour can drive changes in the very nature of the interactions so as to make them progress towards less and less reciprocally harmful coexistence (for full elaboration see Appendix 3).
Our goal has been to produce a simple, tractable scenario for studying the effects of artificial selection on collectives, which while theoretical, is firmly connected to plausible biological experiments. The model could be extended in multiple ways in order to analyse the effects of additional factors, including impact of nonoverlapping generations and variation in the timing of reproduction (which would introduce an element of bethedging), of migration and mixing between collectives (which could be akin to gamete production and zygote formation), and inclusion of more than two kinds of particle types. More complex selective regimes can also be envisaged, such as those that reward collectives based on absolute population size of particle types, which would allow less abstract collective functions to be considered. However, regardless of these refinements, we suspect that our core conclusion will stand firm: collectivelevel selection favours particle dynamics that improve collectivelevel heredity. The ability to reliably reestablish successful adult states of pastgenerations from simpler and potentially noisy initial conditions is adaptive at the collective level.
The mechanism of developmental correction is broadly relevant and extends beyond cells and communities to particles of any kind that happen to be nested within higherlevel selfreplicating structures. As such, the mechanism of developmental correction may be relevant to the early stages in each of the major (egalitarian) evolutionary transitions in individuality (Queller, 1997; Maynard Smith and Szathmary, 1995), where maintenance of particle types in optimal proportions was likely an essential requirement. For example, it is hard to see how protocells cells evolved from lower level components (Takeuchi and Hogeweg, 2009; Baum and Vetsigian, 2017), chromosomes from genes (Smith and Szathmáry, 1993), and the eukaryotic cell from independent bacterial entities (Martin and Müller, 1998) without some kind of selfcorrecting mechanism acting at the collective level.
Beyond these fundamental considerations, the mechanism of developmental correction and the ecological factors underpinning its evolution have important implications for topdown engineering of microbial communities for discovery of new chemistries, new functions, and even new organisms. The minimal recipe involves partitioning communities into discrete packages, provision of a period of time for cell growth, selective criteria that lead to purging of suboptimal collectives and reproduction of optimal collectives to establish the next generation of collectives. These manipulations are readily achieved using millifluidic devices that can be engineered to operate in a Turinglike manner allowing artificial selection on communitylevel traits across thousands of independent communities. As mentioned above, critical tuneable parameters beyond number of communities, mode of selection and population size, are duration of collective generation time and bottleneck size at the moment of birth.
The extent to which the conclusions based on our simple abstract model are generally applicable to the evolution of more complex associations, such as symbioses leading to new forms of life, will require future exploration of a broader range of particlelevel ecologies. Possibilities to make community dynamics more realistic by complexifying mathematical descriptions of particlelevel processes are plentiful (Williams and Lenton, 2007; Zomorrodi and Segrè, 2016). Of particular interest for the evolution of efficient developmental correction are cases when community ecology has multiple attractors (Penn and Harvey, 2004), is highly sensitive to initial conditions (Swenson et al., 2000b), or presents finiteeffect mutations sustaining ‘ecoevolutionary tunnelling’ (Kotil and Vetsigian, 2018). Besides enlarging the spectrum of possible withincollective interactions, future relevant extensions might explore the role of physical coupling among particles and of horizontal transmission between collectives (van Vliet and Doebeli, 2019) in enhancing or hampering efficient inheritance of collectivelevel characters.
Methods
Methods are described in Appendix 1, 2 and 3, interspersed with more technical descriptions of the results presented in the main text.
Appendix 1
Stochastic model
This appendix presents an outline for the approximated stochastic simulation of nested population dynamics described in the main text. We first attend to the case where particle populations are monomorphic without mutation and then move to consider the role of particlelevel mutation. Full implementation of the model is available as Source code 1 and at https://gitlab.com/ecoevomath/estaudel (Doulcier, 2020).
Parameters
First, we list the parameters of the numerical model introduced in the main text. The collective selection regime is specified by a set of collectivelevel parameters that are kept constant along the evolutionary trajectory.
Collectives are comprised of two kinds of selfreplicating particle (red and blue) that carry a different set of traits. Traits can mutate (see the mutation section below), and are represented as global variables.
The state variables ($DxM$ matrices) store the adult state of collectives along a trajectory.
Initial conditions
Initial conditions consist in defining the number of red and blue particles in each collective at the beginning of generation 0.
Outline of the main loop
The main loop of the algorithm applies the sequence growthselectionreproduction for each generation.
Particlelevel ecology
The ecological dynamics of particles is expressed by a multitype birthdeath process with a linear birth rate and a linearly densitydependent death rate. Each type of particle $i$ is characterised by four traits: colour (${c}_{i}$, red or blue), maximum growth rate ${r}_{i}$, and two densitydependent interaction parameters. Particles of the same colour interact according to ${a}_{i}^{\text{intra}}$, whereas particles of different colour interact according to ${a}_{i}^{\text{inter}}$. Interaction terms are in the order of 0.1 and scaled by a carrying capacity term $K$. The dynamic is modelled by a continuoustime Markov jump process with rates:
Each particle of type i . . .  With rate. . . 

Reproduces (add a particle of type i) Dies (remove a particle of type i)  $r}_{i$ ${r}_{i}(\sum _{j}{\delta}_{{c}_{i}={c}_{j}}{x}_{j}{a}_{j}^{\text{intra}}{K}^{1}+\sum _{j}{\delta}_{{c}_{i}\ne {c}_{j}}{x}_{j}{a}_{j}^{\text{inter}}{K}^{1})$ 
${\delta}_{i=i}=1$ if $i=j$ or 0 if $i\ne j$. Additionally ${\delta}_{i\ne i}=1$ if $i\ne j$ or 0 if $i=j$.
The stochastic trajectory of the system is simulated using a Poissonian approximation used in the basic $\tau \text{leap}$ algorithm (Gillespie, 2001), $dt$ is chosen to be small enough so that population size never becomes negative.
Early tests using an exact stochastic simulation algorithm (DoobGillespie SSA, Gillespie, 1976) did not exhibit qualitative changes in the trajectory, but greatly increased the computation duration.
Collectivelevel selection
The collectivelevel selection phase consists in associating each of the D new collectives from generation $m+1$ with a single parent at generation m. In the main text we contrast two regimes: colourneutral and colourselective. In both cases, a fixed proportion $\rho $ of the collective population at generation $m$ is marked for extinction. In the neutral regime, collectives to be eliminated are selected uniformly at random (Appendix 1—figure 1A), whereas in the selective regime they are those that rank the highest in their distance to the optimal colour $\widehat{\varphi}$ (Appendix 1—figure 1B). Each surviving collective produces offspring. Moreover, the remaining collectives from generation $m+1$ are generated by a parent chosen uniformly at random from the set of surviving collectives.
Other selection procedures can be implemented, such as randomly sampling $\rho D$ collectives with weight based on the colourdeviation to the optimal colour. A nonexhaustive exploration of other selection rules indicates that the qualitative results of the model are robust to changes in the selective regime, as long as collectives with an optimal colour are favoured and the collective population does not go extinct.
Collective reproduction is implemented by seeding an offspring collective with a sample of B particles drawn according to the proportion of the parent collective. We assume that the final particle population sizes are big enough so that each reproduction event can be modelled as an independent binomial sample. Smaller population sizes might require simultaneous multinomial sampling of all offspring.
Mutation of particle traits
The complete model adds the possibility for particle trait ($r,{a}^{inter}$) to mutate. Each collective contains two variants of each colour. Whenever a variant goes extinct, the remaining type is called the ‘resident’, and a mutant type is created as follows. First, traits of the resident are copied in one newborn particle, then one of the mutable traits — either the growth rate ($r$) or the intercolour competition trait (${a}_{\text{inter}}$) — is chosen at random, and finally a random value is added that is taken from a uniform distribution over $[\epsilon ,\epsilon ]$ (in Figure 2, $\epsilon =0.1$). Traits are kept positive by taking the absolute value of the result. This process is in the spirit of adaptive dynamics in which invasion of a single new mutant is repeatedly assessed in a monomorphic population.
We checked that relaxing this assumption (i.e. allowing more than two types of each colour in each collective), or waiting for rare mutations to appear did not change the qualitative results. The pseudocode outlined above was modified in order to track the trait value of both resident and mutant types in each collective, rather than having the trait values as global variables.
Variability of the derived particle traits
The single trajectory of the stochastic model discussed in the main text (Figure 2) is obtained by setting initial conditions for both collective composition (fraction of red and blue particles in each collective at the beginning of the first collective generation) and for the value of particlelevel parameters that are subjected to mutations. We used this specific trajectory to illustrate key features of the evolutionary outcome, that are then explained by the deterministic approximation (Appendix 2). Such approximation, discussed below, allows exhaustive exploration of parameter space, that would otherwise prove impossible because of computation demands. However, it provides only indirect insight on the expected variance of the corresponding particlebased stochastic process. Here, we address variability of the evolutionary outcome of stochastic simulations in a few relevant cases.
Appendix 1—figure 2 illustrates the distribution of average colour in a population of collectives (A) and of the average associated traits (interaction rates [B]; growth rates [C]) after $M=10,000$ collective generations for three different values of the red type’s carrying capacity. For each set of parameters, 20 trajectories where simulated. In none of the 60 simulations do populations became entirely monochromatic.
The evolutionary outcome of the stochastic model is highly reproducible, with only a small percentage variation in the average collective colour achieved in different realizations of the stochastic process. Apart from a small systematic bias toward the type with the higher carrying capacity (i.e. the lower ${a}^{intra}$), changing the carrying capacity does not affect the ability of selection to achieve the target colour (Appendix 1—figure 2A).
The effects of collective selection on particle traits instead change depending on the role of such traits in determining the collective colour. Figure 2 of the main text illustrates the convergence, over one single evolutionary trajectory, of the intercolour interaction parameters to specific values. Correspondingly, the average of these parameters in the population also changes little across multiple realizations of the stochastic evolutionary process (Appendix 1—figure 2B). As expected, the values depend on the parameters that do not evolve, that in the cases illustrated here, are the carrying capacities. However, one can notice that as long as interaction rates are bounded to be positive, one interaction (in the specific case where carrying capacities are equal) will in the long run vanish, indicating that one of the types drives the intracollective dynamics. This corresponds to predictions, as discussed in Appendix 2 and 3.
Maximal growth rates, that were observed to change within one single realization of the evolutionary trajectory, also vary considerably across realizations (Appendix 1—figure 2C). This is a consequence of their quasineutrality with respect to collective selection: once the particle population approaches a steadystate within any collective generation, the speed at which it grows at low density is less important.
Different target colours
In the main text, we have considered that the target of selection was a situation when collectives are composed of half red and half blue particles. Here, we repeated the numerical simulations for different target proportions of the two types and discuss how the choice of $\widehat{\varphi}$ affects the outcome of the evolutionary trajectory. We performed ten independent simulations for 10,000 collective generations, starting from the same initial conditions and parameter values, but varying the target of selection.
Appendix 1—figure 3A shows that different target colours, that is, compositions of the collectives, can be selected with an overall accuracy of a few percent. However, the variability among simulations is larger when the target is strongly skewed toward one or another monochromatic state. This reflects an increased variability in the selected particle traits, which manifests when the derived interaction rates are considered (Appendix 1—figure 3B) and maximal growth rate (Appendix 1—figure 3B). While maintaining the previously observed lower variability in interaction rates than in growth rates, all derived traits change more from run to run when the target is extreme. This is likely the consequence of a lower effectiveness in the action of collectivelevel selection, as a considerable part of the population ends up being monochromatic just because of sampling at birth. As a consequence, convergence towards the target is also considerably slower.
Appendix 2
LotkaVolterra deterministic particle ecology
Model for intracollective dynamics
The stochastic ecological dynamics of red (${N}_{0}$) and blue (${N}_{1}$) particles within a collective simulated by the algorithm described in Appendix 1 is approximated (see Figure 3) by the deterministic competitive LotkaVolterra Ordinary Differential Equation:
Here, $\mathbf{\mathbf{r}}=({r}_{0},{r}_{1})$ is the pair of maximal growth rates for red and blue particles. Since the system is symmetric, we consider only the case where red particles grow faster than blue particles (${r}_{0}>{r}_{1}$). The effect of pairwise competitive interactions between cells are encoded in the matrix $\mathbf{\mathbf{A}}={({a}_{ij})}_{i,j\in {\{0,1\}}^{2}}$. All competitive interactions are considered harmful or neutral ($0\le {a}_{ij}$) (see Appendix 3 for more general types of interactions). Here and in the main text, we consider that the four free parameters can evolve at the same time. In the last section of the Results in the main text we explore cases when their variation is constrained by tradeoffs.
So as to explore the space of all possible interaction intensities, no specific mechanism of interaction is assumed, and thus qualitatively different ecological dynamics of the particle populations. Evolutionary trajectories constrained by particle traits are briefly discussed in the last section of the Results (main text).
The carrying capacity of a monochromatic collective is $\frac{K}{{a}_{00}}$ for red particles and $\frac{K}{{a}_{11}}$ for blue particles. $K$ is a scaling factor for the intensity of pairwise interactions that can be used to rescale the deterministic system to match the stochastic trajectories. Without loss of generality, we thus assume that $K=1$.
A natural set of alternate coordinates for the system in Equation 1 are total population size $N:={N}_{0}+{N}_{1}$ and collective colour, defined as the frequency of red individuals $x:=\frac{{N}_{0}}{N}$. In these coordinates, the deterministic dynamics are the solution to the following ODE.
The functions $g$ and $h$ are polynomials in $x$ and $N$, of coefficients:
Monomial  Coefficient in $h$  Monomial  Coefficient in $g$ 

1  ${r}_{0}{r}_{1}$  1  ${r}_{1}$ 
$N$  ${a}_{11}{r}_{1}{a}_{01}{r}_{0}$  $N$  ${a}_{11}{r}_{1}$ 
$Nx$  ${r}_{1}({a}_{10}{a}_{11})+{r}_{0}({a}_{01}{a}_{00})$  $x$  ${r}_{0}{r}_{1}$ 
$xN$  ${r}_{1}(2{a}_{11}{a}_{10}){r}_{0}{a}_{01}$  
${x}^{2}N$  ${r}_{0}({a}_{01}{a}_{00})+{r}_{1}({a}_{10}{a}_{11})$ 
Appendix 2—figure 1 shows an example of the ODE flow in both coordinate systems. Two unstable trivial equilibria and a stable coexistence equilibrium are located at the intersection of the isoclines. Within one collective generation, the dynamics follow such flows for duration $T$, starting from initial conditions on the line ${N}_{0}=B{N}_{1}$ (that is, $N=B$).
Even though Equation 1 is more directly related to the individualbased stochastic simulation, the dynamics of collective colour are understood more easily using Equation 2. Therefore, in the following we will use the latter formulation.
The dynamics of particle types across collective generations are modelled as a piecewise continuous time change (Appendix 2—figure 2), where ${x}^{m}(t)$ is the fraction of red particles at time $t\in [0,T]$ during collective generation $m$.
We focus in particular on the succession ${N}^{m}(T),{x}^{m}(T)$ of collective ‘‘adult’’ states at the end of each successive generation $m$. In the following we note ${\varphi}^{m}={x}^{m}(T)$, which is the adult colour of the collective at the end of the growth phase. At the beginning of each collective generation, we impose that, regardless of the number of cells in the parent, every collective contains the same number of cells ${N}^{m}(0)=B\in \mathbb{R}\forall m$. In contrast, the newborn collective colour ${x}^{m}(0)$ depends on the colour of the parent. In this deterministic model, we consider that there is no stochastic variation or bias due to sampling at birth, so the collective colour of the newborn is equal to that of its parent: ${x}^{m}(0)={x}^{m1}(T)={\varphi}^{m1}$. In the first generation ($m=0$), the population size is B in every collective and the fraction of red particles is chosen uniformly at random between 0 and 1.
The proportion of red particles at adult stage ${\varphi}^{m}$ is obtained from these initial conditions $(B,{\varphi}^{m1})$ by integrating Equation 2. Since these equations are not explicitly solvable, there is no analytic expression for the result of the transient dynamics. However, having constrained the initial conditions to a single dimension, the adult colour is a singlevalued function of the initial composition of the collective, defined as follows.
Definition (Growth function)
Given the set of particle traits $\theta :=(\mathbf{\mathbf{r}},\mathbf{\mathbf{A}})\in E:={[0,\mathrm{\infty})}^{2}\times {[0,\mathrm{\infty})}^{4}$, the bottleneck size $B\in (0,\mathrm{\infty})$ and the duration of the growth phase $T\in (0,\mathrm{\infty})$, the growth function ${G}_{\theta}$ is defined as the application that maps an initial proportion of red particles $\varphi $ to the proportion of red particles after duration $T$. Thus, ${G}_{\theta}(\varphi ,T,B)=x(T)$, with $x(t)$ is the unique solution to the following Cauchy problem:
To simplify notations in the main text, we set $G(\varphi )={G}_{\theta}(\varphi ,B,T)$.
Relation with the stochastic evolutionary model
As explained in the main text, the growth function ${G}_{\theta}$ defines the recurrence relation between the colour of the newborn offspring and its colour at adulthood. If the reproduction process entails no stochasticity, the latter also defines the composition of collectives at the next generation. As a consequence, the iterative application of ${G}_{\theta}$ approximates the change in time of the adult colour when a population with a given, fixed, set of parameters and traits is transferred across collective generations.
The fixed points of this discretetime system allow understanding and classifcation of behaviours observed along stochastic evolutionary trajectories. When mutation rate of particle traits is sufficiently small and growth rates are not vanishing, collectives approach such fixed points in a few collective generations by iteration of the ${G}_{\theta}$ function (Figure 5). The evolutionary trajectory can thus be seen as a succession of fixed points of ${G}_{\theta}$. The surface of the fixed points of ${G}_{\theta}$ as a function of the particle parameters can be computed numerically, as well as its dependence on the collective generation duration $T$ and bottleneck size $B$. A small number of qualitatively different configurations are possible for the fixed points (illustrated in Appendix 2—figure 3 14). In particular, the case in which ${G}_{\theta}$ possesses an internal, stable fixed point (Appendix 2—figure 3) constitutes the optimal solution to constant selection for collective colour, in that it ensures the highest degree of colour reproducibility, on average, across collective generations.
The parameter values that separate regions with qualitatively different fixed points correspond to transcritical bifurcations, where one of the monochromatic fixed points changes its stability. These lines are analytically computed below, thus allowing generalisation of the conclusions drawn from analysing the representative trajectory of Figure 2. In the following, we detail analysis of how the fixed points of the ${G}_{\theta}$ function depend on particle and collectivelevel parameters.
It is worth stressing here that the deterministic model provides a good quantitative approximation of the system with particlelevel and collectivelevel stochasticity provided that fluctuations at both levels are small. This is the case if populations of particles are large (as for instance in the case of bacterial populations) and if mutations of particle traits are rare and of small magnitude. However, in the numerical simulations we performed, the conclusions drawn from ensuing analysis held qualitatively also in the case of large fluctuations.
Fixed Points of the ${G}_{\theta}$ function and their stability
In this paragraph, we list the key properties of the ${G}_{\theta}$ function, that determine how the asymptotic collective colour depends on particle and collective parameters. Proofs of the propositions are provided in the following paragraph.
Proposition 1 (Fixed points of ${G}_{\theta}$)
Let $\theta \mathrm{\in}E$ be a set of particle traits, $T\mathrm{\in}\mathrm{(}\mathrm{0}\mathrm{,}\mathrm{\infty}\mathrm{)}$ the duration of the growth phase and $B\mathrm{\in}\mathrm{(}\mathrm{0}\mathrm{,}\mathrm{\infty}\mathrm{)}$ the bottleneck size.
Then, ${G}_{\theta}\mathit{}\mathrm{(}\mathrm{0}\mathrm{,}T\mathrm{,}B\mathrm{)}\mathrm{=}\mathrm{0}$, and ${G}_{\theta}\mathit{}\mathrm{(}\mathrm{1}\mathrm{,}T\mathrm{,}B\mathrm{)}\mathrm{=}\mathrm{1}$, hence $\varphi \mathrm{=}\mathrm{0}$ and $\varphi \mathrm{=}\mathrm{1}$ are fixed points of ${G}_{\theta}$ $\mathrm{\forall}\theta $.
Moreover, if the stability with respect to $\varphi $ of 0 and 1 is the same, then ${G}_{\theta}$ has at least one fixed point $\varphi \mathrm{\in}\mathrm{(}\mathrm{0}\mathrm{,}\mathrm{1}\mathrm{)}$.
The stability of the monochromatic fixed points with respect to the fraction $\varphi $ can be numerically assessed. The number and stability of the fixed points as a function of collectivelevel parameters B and T are illustrated in Appendix 2—figure 3. Four different cases are considered, corresponding to the four qualitatively different outcomes of particlelevel ecology (competitive exclusion by one or the other type, bistability, coexistence).
The fixed points can be analytically calculated in certain limit cases, corresponding to parameter values when withincollective particle dynamics are exponential or saturating.
Proposition 2 (Quasiexponential growth)
When T is close to 0 and $B\mathit{}{a}_{\text{\mathit{m}\mathit{a}\mathit{x}}}\mathrm{\ll}\mathrm{1}$ with $a}_{\mathrm{m}\mathrm{a}\mathrm{x}}={\mathrm{m}\mathrm{a}\mathrm{x}}_{i,j}\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}{a}_{ij$ the highest element of the competition matrix $\mathrm{A}$, then ${G}_{\theta}$ has only two fixed points $\varphi \mathrm{=}\mathrm{0}$ and $\varphi \mathrm{=}\mathrm{1}$.
Moreover, the monochromatic fixed point $\varphi \mathrm{=}\mathrm{1}$, corresponding to a population completely composed of the (faster) red type of particle, is stable, whereas the fixed point $\varphi \mathrm{=}\mathrm{0}$, corresponding to a population composed of the slow growing type of particle, is unstable.
Proposition 3 (Saturating growth)
As $T\mathrm{\to}\mathrm{\infty}$ the fixed points ${\varphi}^{\mathrm{*}}$of $G}_{\theta$ and their stability correspond to the equilibria ${x}^{\mathrm{*}}$of Equation 1. These equilibria and their stability range are listed below.
Red alone  Blue alone  Coexistence  

Equilibrium ${x}^{*}$  1  0  $\frac{{a}_{11}{a}_{01}}{Tr(A)CoTr(A)}$ 
Stability range  ${a}_{00}<{a}_{10}$  ${a}_{11}<{a}_{01}$  ${a}_{11}>{a}_{01}$ and ${a}_{00}>{a}_{10}$ 
Here $T\mathit{}r\mathit{}\mathrm{(}\mathrm{A}\mathrm{)}\mathrm{:=}{a}_{\mathrm{00}}\mathrm{+}{a}_{\mathrm{11}}$ is the sum of the diagonal (or trace) of $\mathrm{A}$, and $C\mathit{}o\mathit{}T\mathit{}r\mathit{}\mathrm{(}\mathrm{A}\mathrm{)}\mathrm{:=}{a}_{\mathrm{10}}\mathrm{+}{a}_{\mathrm{01}}$ denotes the sum of the antidiagonal elements of $\mathrm{A}$.
By linearising the system in proximity of the fixed points, it is possible to find exactly the bifurcation parameters where one equilibrium changes stability, thus the limit of the region where there exists an interior fixed point ${\varphi}^{*}$. The bifurcation values in the space of the collective parameters $T$ and $B$ delimit the region in Figure 7 where the ${G}_{\theta}$ function has an internal fixed point.
Even when the fixed point is different from the optimal value $\widehat{\varphi}$, it can nonetheless provide a starting point for evolution to optimise collective colour. Extinction of one of the two colours of particles happens instead very rapidly in the region when the monochromatic fixed points are stable, so that collectives have a higher risk of being extinct before inheritanceincreasing mutations appear.
Proposition 4 (Bifurcations of the monochromatic fixed points of ${G}_{\theta}$)
The stability of 0 changes at $\mathrm{(}{T}_{B}^{\mathrm{*}}\mathrm{,}{B}_{B}^{\mathrm{*}}\mathrm{)}$and the stability of 1 at $\mathrm{(}{T}_{R}^{\mathrm{*}}\mathrm{,}{B}_{R}^{\mathrm{*}}\mathrm{)}$such that:
These results allow understanding of the interplay between time scales of particlelevel ecology and collective reproduction, whose relationship changes along an evolutionary trajectory. Appendix 2—figure 4 illustrates the change of the fixed point of ${G}_{\theta}$ with the collective generation duration T for typical particle traits corresponding to the four qualitative classes of asymptotic equilibria for particle ecology. Of particular relevance for understanding the stochastic trajectory illustrated in Figure 2 is panel D.
When the time scale of exponential particle growth is comparable to T, (such as at the beginning of the evolutionary trajectory displayed in Figure 2) Proposition 2 indicates that the system is expected to converge to an allred solution ($\varphi =1$ is the only stable equilibrium). However, at the same time as this fast dynamic occurs, the growth rates change by mutation. Selection during the exponential phase generally favours fast growing mutants, which means that particle populations achieve highdensity conditions in a shorter time. The system then effectively behaves as if T had increased, thus leading selection to ’see’ interaction traits.
When the time scale of collective reproduction is sufficiently slow with respect to the intracollective dynamics, the system crosses the bifurcation point ${T}_{R}^{*}$ (Proposition 4), so that the function ${G}_{\theta}$ now has an internal fixed point (Figure 6). In the stochastic simulations, this means that more collectives are reproducibly found close to the optimal colour. It takes a relatively short time to adjust the particle traits so that the fixed point is close to the optimum $\widehat{\varphi}$. In this case, the deterministic approximation produces a close to perfect inheritance of the collective colour. However, fluctuations in particle numbers and in composition at birth still result in a large variance of colours among collectives in the stochastic system.
Here, starts the last and slowest phase of the evolutionary trajectory, which results in colour variance reduction through improvement of the ability of particles to correct variations in colour. This is achieved by attaining faster the particle ecological equilibrium, so that fluctuations are more efficiently dampened by demographic dynamics. As a consequence, the conditions described by Proposition 3 will be met. This allows identification of a surface in parameter space, where the fixed point of the ${G}_{\theta}$ function ${\varphi}^{*}$ identifies with the ecological equilibrium ${x}^{*}$, that contains evolutionary equilibria. The ecological equilibrium is displayed in Appendix 2—figure 4 as a function of the crosscolour interaction parameters ${a}_{01}$ and ${a}_{10}$. In the regime where particle and collective time scales are well separated ($r\gg 1/T$), then interaction parameters that correspond to the optimal colour satisfy the following relationship:
This relation identifies the white dotted line in Appendix 2—figure 4 (and in Figure 8B). Once it is attained, mutations cause the deterministic system to move neutrally on this surface. As the stochastic simulation shows, particle parameters keep evolving directionally so as to reduce phenotypic variance. This is achieved by making ${G}_{\theta}$ increasingly flatter in the vicinity of the fixed point, so that the target colour is not only more stable, but it is reached in fewer collective generations.
Successive events of mutation and substitution progressively lead to a growing asymmetry in the ecological relationship between the two types of particles: that with smaller carrying capacity becomes insensitive to the other colour; the latter instead experiences competition, so that its growth is curbed and optimal proportion of colours is eventually realized (see Figure 3D). As illustrated by Appendix 2—figure 6, this conclusion is independent of what is the maximum growth rate, that only affects the advantage of one type at initial stages of growth. Indeed, the position on the intercept between the optimality line and the y axis in Appendix 2—figure 5 and Appendix 2—figure 6 only depends on the difference between intratype competition parameters. A consequence of this is that fastgrowing types systematically display, when they are also those with a larger carrying capacity, a population overshoot.
Proofs
In this section, we present the proof of propositions 2–4 above.
Proof of Proposition 1
In the $(N,x)$ coordinates we have seen that $\frac{df(t)}{dt}=g(x(t),N(t))$ and that $x(t)=0$ and $x(t)=1$ are trivial roots of the polynomial $g$ (Equation 2). Hence 0 and 1 are always fixed points of ${G}_{\theta}$.
Since $g$ and $h$ (from Equation 2) are polynomials of $(N,x)$, they are smooth (of differentiability class ${C}_{\mathrm{\infty}}$) on ${\mathbb{R}}^{2}$. Thus, the global flow corresponding to the Cauchy problem is also smooth on $\mathbb{R}}^{2$ is the partial application of the global flow to the case where ${N}_{0}=B$. Therefore, ${G}_{\theta}$ is continuous on $[0,1]$.
Moreover, let us suppose that 0 and 1 are both unstable. Then ${G}_{\theta}^{\prime}(0,T,B)>1$ and ${G}_{\theta}^{\prime}(1,T,B)>1$, with ${G}_{\theta}^{\prime}$ the derivative of ${G}_{\theta}$ with respect to its first variable. As a consequence there is an $\epsilon \in \mathbb{R}$ such that ${G}_{\theta}(\epsilon ,T,B)>\epsilon $ and ${G}_{\theta}(1\epsilon ,T,B)<1\epsilon $. Since ${G}_{\theta}$ is continuous, there is at least one $c\in [0,1]$ such that ${G}_{\theta}(c,T,B)=c$ by virtue of the intermediate value theorem.
In practice, we never encountered cases when more than one internal fixed point was present. However multistability is expected to occur if the equations describing particle ecology had higherorder nonlinearities.
Proof of Proposition 2
Since the nonlinear terms in Equation 3 are smaller than $B{a}_{\text{max}}$, and this is negligible with respect to 1, particle ecology is approximated by its linearization as long as the population size remains close to the bottleneck value. Around $t=0$, the Cauchy problem can be written as:
In the $(N,x)$ coordinates, the total population size grows exponentially and is decoupled from the colour. On the other hand, $x(t)$ follows the logistic differential equation:
which can be integrated.
For T sufficiently small for the population to be in exponential growth phase, the growth map ${G}_{\theta}$ can be approximated by the solution ${\stackrel{~}{G}}_{\theta}$ of Equation 5:
This function is strictly convex (or concave) on (0, 1) depending on the sign of ${r}_{0}{r}_{1}$:
Since $0<x<1$ and $t>0$, $\frac{{\partial}^{2}{\stackrel{~}{G}}_{\theta}(x,T)}{\partial {x}^{2}}$ is of the same sign as $1{e}^{t({r}_{0}{r}_{1})}$, that is strictly positive if ${r}_{1}>{r}_{0}$, or strictly negative if ${r}_{0}>{r}_{1}$.
Thus, $\stackrel{~}{G}$ is strictly convex on $(0,1)$ if ${r}_{1}>{r}_{0}$, and strictly concave on $(0,1)$ if ${r}_{0}>{r}_{1}$. In the first case, red colour $x=0$ is an unstable equilibrium and blue colour $x=1$ a stable one, and viceversa in the second case. Note that the segment $s=[(0,0),(1,1)]$ is a chord of $\stackrel{~}{G}$. Therefore, the strictly convex (resp. concave) $\stackrel{~}{G}$ do not intersect $s$ except in $(0,0)$ and $(1,1)$.
Proof of Proposition 3
When the collective generation time T is much longer than the demographic time scale, the populations within droplets at the adult stage are well approximated by the equilibrium solution of the LotkaVolterra Equation 1. Solving simultaneously equations $\frac{d{N}_{0}}{dt}=0$ and $\frac{d{N}_{1}}{dt}=0$ (or equivalently $\frac{dN}{dt}=0$, $\frac{dx}{dt}=0$) yields the four equilibria listed below in both coordinate systems. Linear stability analysis allows one to determine the parameter intervals where these are stable, listed below.
Name  Red alone  Blue alone  Coexistence  Extinction 

Position in $[{N}_{0},{N}_{1}]$  $[\frac{1}{{a}_{00}},0]$  $[0,\frac{1}{{a}_{11}}]$  $[\frac{{a}_{11}{a}_{01}}{det(A)},\frac{{a}_{00}{a}_{10}}{det(A)}]$  [0,0] 
Position in $[N,f]$  $[\frac{1}{{a}_{00}},1]$  $[\frac{1}{{a}_{11}},0]$  $[\frac{Tr(A)CoTr(A)}{det(A)},\frac{{a}_{11}{a}_{01}}{Tr(A)CoTr(A)}]$  Undefined 
Condition for stability  ${a}_{00}<{a}_{10}$  ${a}_{11}<{a}_{01}$  ${a}_{11}>{a}_{01}$ and ${a}_{00}>{a}_{10}$  Never 
Proof of Proposition 4
We consider the case when the fixed point 0 changes stability. The stability of the fixed point 1 can be studied analogously.
We aim at identifying the values of the collective parameters $(T,B)$ where a fixed point with $x=0$ changes stability through a transcritical bifurcation. The difficulty lies in the fact that one needs to estimate the dynamics of the second variable $N$ in order to study the stability of the 2D system. Luckily, this can be done in the limit when the collective contains almost exclusively particles of one single colour (in this case, blue).
Total population size is in this case decoupled from colour and Equation 2 can be integrated with initial condition $({N}_{0}{x}_{0})=(B,0)$, yielding the following trajectory:
As long as $x$ is small, the time derivative is approximated by the nonautonomous system:
with $h$ as defined in Equation 2:
Solving this equation allows computation of the adult colour as a function of the parameters $T$ and $B$. At the bifurcation point $({T}^{*},{B}^{*})$, where stability of the 0 fixed point changes, the newborn colour is the same as the adult colour: $x({T}^{\ast},{B}^{\ast})=x(0,{B}^{\ast})$. $({T}^{\ast},{B}^{\ast})$ are then solutions of the integral equation:
Solving for ${B}^{*}$ we get:
With $\alpha ={r}_{0}{r}_{1}\frac{{a}_{00}{a}_{10}}{{r}_{0}{a}_{00}{r}_{1}{a}_{10}}$
Figure 7 shows that this approximation retrieves accurately the numerically computed bifurcation lines.
Appendix 3
Beyond competitive interactions
In the main text, the model is limited to deleterious (competitive) interactions, whereby the intercolour interactions are such that one colour always has a negative densitydependent effect on the growth of the other colour. This choice reflects the expectation that, as in the example of the droplet experiment, competition would be the main driver of interactions among two previously unconnected bacterial species growing on a common resource.
However, it is interesting to consider what happens if interactions of other types are taken into account. This section presents two additional interaction classes: exploitative (e.g. predatorprey or parasitic) interactions (${a}_{01}<0<{a}_{10}$ or ${a}_{10}<0<{a}_{01}$), and mutualism (${a}_{01}<0$ and ${a}_{10}<0$). Note that exploitative interactions correspond to nonobligate predation (or parasitism): the predator (or the parasite) can sustain a nonnull population density even in the absence of the prey (or host). Obligate interactions can be modelled in the limit where the carrying capacity of the exploiters tends to zero. Similarly, unless the carrying capacities of both types tend to zero, mutualism is also not obligate.
The deterministic model (Appendix 2 Equation 1) can be straightforwardly extended to these larger parameter ranges, allowing prediction of the outcome of the evolutionary dynamics. Appendix 2—figure 5 presents the position of the stable fixed point of ${G}_{\theta}$ for negative intercolour interaction parameters when the time scale of collective reproduction is sufficiently slow with respect to the intracollective dynamic $\left(r\gg \frac{1}{T}\right)$. Appendix 3—figure 1 represents the same fixed point for extended values of the interaction parameters. Each quadrant corresponds to a different kind of ecological interaction (competition, exploitation and mutualism). While competition is discussed in the main text, other kinds of interactions are discussed below.
Whenever evolution leads to separation in timescales between particle and collective dynamics (typically by increasing particle growth rates), the evolutionary dynamic within each quadrant will be akin to the competition case: first, successive invasions of mutants bring the system on the manifold of interaction parameters where the equilibrium colour ${x}^{*}$ corresponds to the selected colour $\widehat{\varphi}$ (Appendix 2 Equation 4); second, the system evolves on this manifold toward regions of faster ecological convergence.
Exploitative interactions
Consider first that ${a}_{01}<0<{a}_{10}$, so that red individuals (type 0) are the exploiters (predators or parasites) and blue individuals (type 1) are the exploited (prey or hosts). This corresponds to the bottomright quadrant of Appendix 3—figure 1.
When the exploiter competes with the exploited more than with other exploiters (${a}_{00}<{a}_{10}$), the stable equilibrium for Appendix 2 Equation 1 is monochromatic red (see Appendix 2 Proposition 3). However, when the the competitive interaction is stronger among exploiters (${a}_{00}>{a}_{10}$), the coexistence equilibrium is stable and correspond to a colour of $\frac{{a}_{11}{a}_{01}}{Tr(A)CoTr(A)}$ (see Appendix 2 Proposition 3). In the coexistence region of the bottomleft quadrant, the collective colour can take up values in the range $[\frac{{a}_{11}}{{a}_{11}+{a}_{00}},1]$ bounded on one side by the extinction of the exploited type, and on the other side by the ratio of carrying capacity obtained when there is no competitive interactions. This range does not contain the target colour $\widehat{\varphi}=0.5$ in the example illustrated in Appendix 3—figure 1 because the exploiter has a higher carrying capacity than the exploited $\left(\frac{1}{{a}_{11}}<\frac{1}{{a}_{00}}\right)$. It is thus impossible to achieve the target colour $\widehat{\varphi}=0.5$ in the lower right quadrant. Unless the nature of ecological interactions changes, over evolutionary times the population of collectives will move towards the absence of interactions. Here, collectives of colour close to the optimum can still manifest in the stochastic individualbased model as the effect of sampling at dilution and fluctuations during growth, but the average colour will remain off target.
In the upper left quadrant, where ${a}_{10}<0<{a}_{01}$, blue individuals (type 0) are the exploiters (predators or parasites) and red individuals are the exploited (prey or hosts). In the coexistence region (i.e., when the competition is stronger among exploiters than between exploiters and exploited ${a}_{01}<{a}_{11}$), the collective colour can take up values in the range $[0,\frac{{a}_{11}}{{a}_{11}+{a}_{00}}]$, that contains the target $\widehat{\varphi}=0.5$ in the example illustrated in Appendix 3—figure 1 because the exploiter has a lower carrying capacity than the exploited $\left(\frac{1}{{a}_{11}}<\frac{1}{{a}_{00}}\right)$. Within the boundaries of this quadrant, longterm evolution will tend to reduce the (positive) interaction ${a}_{01}$, thus mitigating exploitation, while at the same time making the (negative) interaction term ${a}_{10}$ smaller, that is increasing the advantage provided by the exploited to the exploiter. Again, if the nature of ecological interactions does not change, the system will evolve towards a case where one particle type (in this case, the exploiter, that also has smaller carrying capacity), gets decoupled from the demography of the other type.
Mutualistic interactions
Consider now that red and blue individuals have nonobligate, mutualistic interactions. This corresponds to both intercolour competition traits ${a}_{01}$ and ${a}_{10}$ being negative, and can be visualised in the bottomleft quadrant of Appendix 3—figure 1.
When interactions are mutualistic, the equilibria of Appendix 2 Equation 1 resulting from the extinction of one type or the other are never stable (${a}_{00}<{a}_{10}$ or ${a}_{11}<{a}_{01}$ are never fulfilled since ${a}_{01}<0<{a}_{00}$ and ${a}_{01}<0<{a}_{11}$ respectively, see Appendix 2 Proposition 3) and the solution is always coexistence.
Depending on the sign of the determinant of the interaction matrix $\mathbf{\mathbf{A}}$, the dynamics within a collective generation can have two qualitatively different behaviours. When mutualistic interactions have small intensity, the intracollective dynamics attains an equilibrium, even though the population size of one of the two particle types is larger than their respective carrying capacities (because of the positive contribution of the other type). However, when $\frac{{a}_{00}{a}_{11}}{{a}_{01}{a}_{10}}<1$ (white line in Appendix 3—figure 1), the populations would undergo unbounded growth if they were not embedded in collectives of finite lifetime.
This unrealistic feature of simplified LotkaVolterra models, that makes them not well fit to describe mutualistic interactions, is not particularly problematic in our case, since the fraction $x$ of red particles converges to its equilibrium value even though the population size $N$ diverges. Appendix 3—figure 1 illustrates how the withincollective particle dynamics changes along an evolutionary trajectory towards increasingly strong mutualistic interactions. If no physical or biological constraint keeps parameters in the region of bounded growth, then adherence of the model predictions from reality would need to be tested by studying their robustness to inclusion of other sources of colimitation, such as, for instance, nutrient exhaustion.
Transition between classes of ecological interactions
So far, we have examined different ecological interactions separately, as they are usually considered to belong to separate categories. It is, however, conceivable that mutations induce a competitive interaction to become exploitative, or an exploitative interaction to become mutualistic (Sørensen et al., 2019), which in our case would be the consequence of the intercolour interaction parameters changing sign. In this case, the prediction of the deterministic model is that, as long as there are no constraints that limit the variation of the interaction parameters, selection at the collective level should lead particles to become more and more benign towards each other, and eventually reach a regime where their exchanges are mutualistic.
Data availability
The source code for all simulations and figures in the manuscript is available as a zip file uploaded with the manuscript and in a public git repository (https://gitlab.com/ecoevomath/estaudel; copy archived at https://github.com/elifesciencespublications/estaudel).
References

Mortality causes universal changes in microbial community compositionNature Communications 10:1–9.https://doi.org/10.1038/s41467019099250

Millifluidic droplet analyser for microbiologyLab on a Chip 11:4057–4062.https://doi.org/10.1039/c1lc20545e

An experimental framework for generating evolvable chemical systems in the laboratoryOrigins of Life and Evolution of Biospheres 47:481–497.https://doi.org/10.1007/s110840169526x

Ecological scaffolding and the evolution of individualityNature Ecology & Evolution 4:426–436.https://doi.org/10.1038/s4155901910869

Evolutionary strategies and developmental constraints in the cellular slime moldsThe American Naturalist 119:530–552.https://doi.org/10.1086/283930

Unifying evolutionary dynamics: from individual stochastic processes to macroscopic modelsTheoretical Population Biology 69:297–321.https://doi.org/10.1016/j.tpb.2005.10.004

Microbial communities as experimental unitsBioScience 61:398–406.https://doi.org/10.1525/bio.2011.61.5.9

Nascent multicellular life and the emergence of individualityJournal of Biosciences 39:237–300.https://doi.org/10.1007/s1203801494205

A general method for numerically simulating the stochastic time evolution of coupled chemical reactionsJournal of Computational Physics 22:403–434.https://doi.org/10.1016/00219991(76)900413

Approximate accelerated stochastic simulation of chemically reacting systemsThe Journal of Chemical Physics 115:1716–1733.https://doi.org/10.1063/1.1378322

BookDarwinian Populations and Natural SelectionOxford: Oxford university Press.https://doi.org/10.1086/656856

Experimental studies of group selection: what do they tell Us about group selection in nature?The American Naturalist 150 Suppl 1:S59–S79.https://doi.org/10.1086/286050

A reexamination of the stochastic corrector modelProceedings of the Royal Society B: Biological Sciences 262:29–35.https://doi.org/10.1098/rspb.1995.0172

What is "epi" about epigenetics?Annals of the New York Academy of Sciences 981:97–110.https://doi.org/10.1111/j.17496632.2002.tb04914.x

Dynamical systems approach to Higherlevel heritabilityJournal of Biological Physics 28:799–804.https://doi.org/10.1023/A:1021215511897

Selection at the level of the community: the importance of spatial structureTrends in Ecology & Evolution 17:83–90.https://doi.org/10.1016/S01695347(01)023850

Closedform stochastic solutions for nonequilibrium dynamics and inheritance of cellular components over many cell divisionsProceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 471:20150050.https://doi.org/10.1098/rspa.2015.0050

Emergence of evolutionarily stable communities through ecoevolutionary tunnellingNature Ecology & Evolution 2:1644–1653.https://doi.org/10.1038/s4155901806557

The units of selectionAnnual Review of Ecology and Systematics 1:1–18.https://doi.org/10.1146/annurev.es.01.110170.000245

BookOrigin of Eukaryotic Cells: Evidence and Research Implications for a Theory of the Origin and Evolution of Microbial, Plant, and Animal Cells on the Precambrian EarthNew Haven, USA: Yale University Press.

BookModelling artificial ecosystem selection: A preliminary investigationIn: Banzhaf W, Ziegler J, Christaller T, Dittrich P, Kim J. T, editors. European Conference on Artificial Life. Springer. pp. 659–666.https://doi.org/10.1007/9783540394327_71

ConferenceThe role of nongenetic change in the heritability, variation, and response to selection of artificially selected ecosystemsArtificial Life IX: Proceedings of the Ninth International Conference on theSimulation and Synthesis of Artificial Life.

Darwin was right: where now for experimental evolution?Current Opinion in Genetics & Development 47:102–109.https://doi.org/10.1016/j.gde.2017.09.003

The origin of chromosomes. I. selection for linkageJournal of Theoretical Biology 164:437–446.https://doi.org/10.1006/jtbi.1993.1165

The role of exploitation in the establishment of mutualistic microbial symbiosesFEMS Microbiology Letters 366:fnz148.https://doi.org/10.1093/femsle/fnz148

Artificial selection of microbial ecosystems for 3chloroaniline biodegradationEnvironmental Microbiology 2:564–571.https://doi.org/10.1046/j.14622920.2000.00140.x

BookAdaptation in MetapopulationsChicago: University of Chicago Press.https://doi.org/10.7208/chicago/9780226129877.001.0001

The sociogenesis of insect coloniesScience 228:1489–1495.https://doi.org/10.1126/science.228.4707.1489

Reviving the superorganismJournal of Theoretical Biology 136:337–356.https://doi.org/10.1016/S00225193(89)801699

Synthetic ecology of microbes: mathematical models and applicationsJournal of Molecular Biology 428:837–861.https://doi.org/10.1016/j.jmb.2015.10.019
Decision letter

Wenying ShouReviewing Editor; Fred Hutchinson Cancer Research Center, United States

Patricia J WittkoppSenior Editor; University of Michigan, United States

Wenying ShouReviewer; Fred Hutchinson Cancer Research Center, United States

Sara MitriReviewer; University of Lausanne, Switzerland

Alvaro SanchezReviewer; Yale University, United States
In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.
Acceptance summary:
Doulcier et al. tackled multilevel community selection: selection at the level of individual cells for fast growth and at the level of whole communities for a desired community property. The authors showed that when selecting communities based on species ratios, the growth rates and the interspecies competition effects can evolve such that a desired species ratio (here 1:1) can be achieved. This work contributes to a growing interest in using multilevel selection to achieve desired community properties.
Decision letter after peer review:
Thank you for submitting your article "Ecoevolutionary dynamics of nested Darwinian populations and the emergence of communitylevel heredity" for consideration by eLife. Your article has been reviewed by three peer reviewers, including Wenying Shou as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen by Patricia Wittkopp as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity:); Sara Mitri (Reviewer #2); Alvaro Sanchez (Reviewer #3).
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
Summary:
Doulcier et al. tackled multilevel community selection: selection at the level of individual cells for fast growth and at the level of whole communities for a desired community property. The authors showed that when selecting communities based on species ratios, the growth rates and the interspecies competition effects can evolve such that a desired species ratio (here 1:1) can be achieved. The evolved mechanism is rather intuitive: species Blue grows faster than species Red, but species Blue is subject to competitive inhibition by species Red. This allows the species ratio to be 1:1. Consistent with intuition, when possibilities for competition coefficients are eliminated, growth rates of the two species evolved to be equal. When basal grow rates were not allowed to evolve, asymmetric competition arose to achieve equal postcompetition growth rates. All three reviewers are positive about this work.
Essential revisions:
1) Allowing ecological interactions to evolve to be positive.
2) Statistics of simulation repeats.
3) What happens if selecting for uneven ratios when noise is present?
4) What does "rapidly" really mean?
To help you revise, the original reviews are attached.
Reviewer 1:
1) Define "interactions".
2) "Developmental correction" is a biologist's way of saying "steadystate species ratio" or "attractor/fixed point". There are multiple ways of achieving steadystate species ratios. The authors had limited species interactions to competition (nonpositive coefficients). If they allow interactions to change signs, then commensal or mutually beneficial interactions can arise to stabilize strain ratios. "Developmental correction thus selects for maximal asymmetry in interactions" may no longer hold for positive interactions. This needs to be checked.
Reviewer #2:
In this manuscript, the authors set out to explore whether it's possible to select for phenotypes at the community level, that arise in spite of competition between community members at the individual level. Using a computational model, they show that this is indeed possible, but relies on communities being properly enclosed and only passing on their individuals to their own offspring communities. They also analyse the strategy allowing for the evolution of community phenotypes. They show that community members evolve individual phenotypes that allow them to robustly converge to the community property under selection, independently of the initial community state. They posit that a similar mechanism might explain major evolutionary transitions, allowing lower level entities to be selected as whole entities.
Overall, we like the manuscript. It makes, in our opinion, a major contribution in understanding how competition between individual community members can be overcome – or even used – to allow for selection of the community phenotype. The idea that selection can favor individuallevel parameters that have a functional form that increases resemblance between parent and offspring and robustness to initial conditions is very interesting and overcomes some of the issues with community selection observed in previous studies (Xie et al). The authors also nicely show that once this strategy has evolved, it is ecologically stable, and they explore how the evolutionary trajectory depends on starting conditions.
We have some general thoughts that we think are worth discussing further in the manuscript.
First, what we find was missing were statistics on what happens on different runs of the simulation. As far as we could tell, all data shown is based on individual "example" runs. The authors claim that these are general, but there is no data to support that. Please add statistics on the repeatability of the results. A related point: is the outcome always qualitatively similar? Specifically that a) the faster grower becomes less competitive with respect to growth rate, and b) that the emerging population regulation consists of the slower grower suppressing the faster grower while c) the faster grower does not affect the slower. It was not clear to us why, and we believe this is important to understanding the conclusion of the paper.
Second, we believe the authors could elaborate more on the specificity of their system and how their findings depend on them. It is not clear to what extent the same result would be observed in a different system that (1) has more than 2 types, (2) where selection acts on total abundance rather than ratios of two colors, (3) where the selective regime is different (the authors say "A nonexhaustive exploration of other selection rules indicates that the qualitative results of the model are robust to changes in the selective regime, as long as collectives with an optimal colour are favoured and the collective population does not go extinct." It would be worth expanding on this in the supplement); and finally, (4) where the phenotype being selected for is not intermediate. Would this work equally well if the phenotype selected is almost blue or almost red? In short, under what general conditions do you expect to see the same?
From the ecological literature, it is known that coexistence between two competing types requires some limiting force for the strongest grower (e.g. predator). Is that simply what you are selecting for? It may be good to make a link to this literature.
We are also unsure about the parallels made to developmental processes and think this would merit some clarification. Does one see this strategy where a slowergrowing type limits the growth of a fastergrowing one? Or are the authors hypothesizing that this might explain how smooth developmental processes occur?
Finally, Figure 2 is lacking in visual clarity (color scheme). The selection method should also be better explained in the main text.
Please note that we did not work through all the math nor run the code.
Reviewer #3:
I very much enjoyed this paper. I do not have any significant concerns.
[Editors' note: further revisions were suggested prior to acceptance, as described below.]
Thank you for resubmitting your article "Ecoevolutionary dynamics of nested Darwinian populations and the emergence of communitylevel heredity" for consideration by eLife. Your article has been rereviewed by two of the three original reviewers, and the evaluation has been overseen by a Reviewing Editor and Patricia Wittkopp as the Senior Editor.
We would like to draw your attention to changes in our revision policy that we have made in response to COVID19 (https://elifesciences.org/articles/57162). Specifically, we are asking editors to accept without delay manuscripts, like yours, that they judge can stand as eLife papers without additional data, even if they feel that they would make the manuscript stronger. Thus the revisions requested below only address clarity and presentation.
The only thing that we believe is still missing is a paragraph in the Discussion on the generality of these findings. We understand that adding more than 2 types, trying different selection regimes or selecting on total abundance is beyond the scope of the current manuscript, but we still believe these questions may remain for a reader of the paper, so it is worth bringing them up in the Discussion as future directions.
https://doi.org/10.7554/eLife.53433.sa1Author response
Essential revisions:
1) Allowing ecological interactions to evolve to be positive.
We have performed further analyses in which interactions evolve to be positive. This is included in a new supplementary document (Appendix 3) and mentioned in the text.
Since the deterministic approximation for the intracollective particle population dynamics is not dependent on the sign of the interaction rates, we were able to use the steadystate of the G function to study how colour can reach its target value by evolution of particle traits. In Appendix 3 we report the existence (or not) of evolutionary solutions for three additional classes of ecological interaction: exploitation of the type that has larger or smaller carrying capacity (interaction parameters of opposite sign), and mutualism (positive interaction parameters).
We firstly considered the four classes of interaction separately, but show that if interaction parameters are free to change sign, maximal asymmetry of competitive interactions will be replaced by exploitative interactions first, then by mutualism of increasing strength, as correctly intuited by reviewer 1.
We believe that in these cases, even more than for competitive interactions, constrains on particlelevel parameters should play a fundamental role in determining the solutions that are actually achievable. For instance, increasing mutualism would lead to the urrealistic situation – typical of LotkaVolterra models – where mutual benefits are so high that the two evolved types would, in the absence of the imposed collective life cycle, grow indefinitely.
We have commented on this point in the Discussion, where we refer to the new Appendix 3, and introduced earlier the concept that different types of ecological interactions can be conceived (Results subsection “A nested model of collective evolution”).
2) Statistics of simulation repeats.
We perfectly agree with the reviewers that the extent to which our reference simulation is representative of other realizations of the same stochastic process needed to be quantified. As mentioned above, the computational demands associated with such simulations limits the possibility of obtaining this information systematically for any parameter value. We thus realized multiple (20) simulations by fixing the nonevolving carrying capacities to values that illustrate different regions in the space of ecological interactions.
We have expanded Appendix 1 so as to include statistics of the variation of the target colour at the end of the simulations, together with the corresponding variation of the underlying intercolour interaction and maximal growth parameters for the two particle types of different colour. As already observed when discussing the persistence of the evolutionary trajectory, at the evolutionary equilibrium growth rate becomes largely irrelevant (as long as it is large enough), and therefore is more variable than intercolour interaction parameters.
We have expanded the first section of the Results in order to include this information and refer to the new paragraph “Variability of the derived particle traits” within Appendix 1.
3) What happens if selecting for uneven ratios when noise is present?
By looking at the deterministic model (Appendix 2—figure 5), it is clear that in the absence of individuallevel stochasticity the target colour influences the relative values of the evolving traits and also the resilience of the fixed point to fluctuations (given by the gradient of the fixed point position). However, we expect that fluctuations associated notably with sampling might have major effects on the evolutionary trajectory.
In order to test the effect of the stochastic individualbased process, we ran multiple (10) simulations for different levels of the target colour and expanded Appendix 1 accordingly (paragraph “Different target colours”). As shown in the new Appendix 1—figure 3, the target colour is always successfully achieved, with variation among independent simulations being slightly larger for more extreme colours. However, the variation in the corresponding particle traits is broader, and increases for more extreme target values, together with the time scale associated to reaching the evolutionary equilibrium. These additional observations have been mentioned in the main text, subsection “Selection on collectives drives the evolution of particle traits”, where we refer to the new paragraph “Different target colours” within Appendix 1.
4) What does "rapidly" really mean?
We realise that this term may appear uninformative when addressing the neutral collective selection regime. This was mostly a consequence of mention of issues pertaining to timescales, prior to the later more complete elaboration. Consequently we have removed this early mention and added later that the time scale of early colour variation under collective selection, driven by differences in growth rate, is comparable to that observed in the neutral regime.
We have also added a term of comparison with the stochastic corrector mechanism, that acts on constant particle traits, therefore in conditions where ecological dynamics are dominated by differences in maximal growth rate (Appendix 1—figure 5). The evolutionary dynamics in this case explains the transient regime on short time scales, when the average colour is still very different from the target.
Taking also into account some of the detailed comments, we have expanded the description of the numerical simulation in subsection “Selection on collectives drives the evolution of particle traits”.
To help you revise, the original reviews are attached.
Reviewer 1:
1) Define "interactions".
We now define interactions in the second paragraph of the Results, where we also mention the possibility of having different kinds of ecological interactions other than competition.
2) "Developmental correction" is a biologist's way of saying "steadystate species ratio" or "attractor/fixed point". There are multiple ways of achieving steadystate species ratios.
The steady state of the growth function corresponds to a trajectory connecting the state of a collective at birth to that at adult stage, and its value does not necessarily correspond to the equilibrium of the particle population dynamics (rather, it captures some property of the transient). We therefore consider that it has many analogies to developmental processes. Moreover, we contrast it with stochastic correction, where correction is not determined by particle ecology that primes, but by fluctuations due to the sampling. We feel therefore that the metaphor is useful to convey our message.
The authors had limited species interactions to competition (nonpositive coefficients). If they allow interactions to change signs, then commensal or mutually beneficial interactions can arise to stabilize strain ratios. "Developmental correction thus selects for maximal asymmetry in interactions" may no longer hold for positive interactions. This needs to be checked.
See Essential revision #1.
Reviewer #2:
[…]
First, what we find was missing were statistics on what happens on different runs of the simulation. As far as we could tell, all data shown is based on individual "example" runs. The authors claim that these are general, but there is no data to support that. Please add statistics on the repeatability of the results. A related point: is the outcome always qualitatively similar? Specifically that a) the faster grower becomes less competitive with respect to growth rate, and b) that the emerging population regulation consists of the slower grower suppressing the faster grower while c) the faster grower does not affect the slower. It was not clear to us why, and we believe this is important to understanding the conclusion of the paper.
See Essential revision #2.
Second, we believe the authors could elaborate more on the specificity of their system and how their findings depend on them. It is not clear to what extent the same result would be observed in a different system that (1) has more than 2 types, (2) where selection acts on total abundance rather than ratios of two colors, (3) where the selective regime is different (the authors say "A nonexhaustive exploration of other selection rules indicates that the qualitative results of the model are robust to changes in the selective regime, as long as collectives with an optimal colour are favoured and the collective population does not go extinct." It would be worth expanding on this in the supplement); and finally, (4) where the phenotype being selected for is not intermediate. Would this work equally well if the phenotype selected is almost blue or almost red? In short, under what general conditions do you expect to see the same?
There are many interesting possible extensions of this model, and we have explored just a few, maintaining the basic assumptions that only two colours are present in the population and that optimal collective function requires both colours to be represented (whereas yield maximization can also be obtained by one type alone).
After having tried different selection schemes and realized that they yielded similar results (as long as the collective population did not go extinct), we chose a selection rule that was simple to implement, both in the individualbased numerical model and in possible applications. A systematic exploration of how different selection rules affect the probability of reaching the target colour within a finite time would however demand a much larger study that goes beyond the scope of this work.
However, as noted in response to Essential revisions #1 and #3, we have investigated the effects of changes in interaction rate and target colour.
From the ecological literature, it is known that coexistence between two competing types requires some limiting force for the strongest grower (e.g. predator). Is that simply what you are selecting for? It may be good to make a link to this literature.
The competitive LotkaVolterra equations are indeed widely used in ecology. They are so widespread that we decided that no reference to specific literature was warranted.
If the “strongest grower” is defined based on its exponential growth rate r, then it is important to note that in the derived population not this parameter, but the interaction traits, are the chief determinants of the adult colour. See our response to Essential revision #2.
We are also unsure about the parallels made to developmental processes and think this would merit some clarification. Does one see this strategy where a slowergrowing type limits the growth of a fastergrowing one? Or are the authors hypothesizing that this might explain how smooth developmental processes occur?
Selection of the optimal colour, acting on adult collectives, changes the continuous dynamics that link, within collectives, the newborn to the adult state. In particular, this dynamic is increasingly “canalized”, as illustrated in Figure 3. We have improved Figure 3 in order to make this concept clearer.
Finally, Figure 2 is lacking in visual clarity (color scheme). The selection method should also be better explained in the main text.
We have improved presentation of the model in the main text (following the detailed comments) and the colour contrast of Figure 2. However, we refer the reader to Appendix 1 for all the details on implementation of the selective process. On the basis of the deterministic system and on nonexhaustive exploration of alternative selection schemes, we believe that those details are not essential to understand the key messages of this study.
[Editors' note: further revisions were suggested prior to acceptance, as described below.]
The only thing that we believe is still missing is a paragraph in the Discussion on the generality of these findings. We understand that adding more than 2 types, trying different selection regimes or selecting on total abundance is beyond the scope of the current manuscript, but we still believe these questions may remain for a reader of the paper, so it is worth bringing them up in the discussion as future directions.
As requested we have added a single new paragraph to the Discussion that reads as follows:
“Our goal has been to produce a simple, tractable scenario for studying the effects of artificial selection on collectives, which while theoretical, is firmly connected to plausible biological experiments. The model could be extended in multiple ways in order to analyse the effects of additional factors, including impact of nonoverlapping generations and variation in the timing of reproduction (which would introduce an element of bethedging), of migration and mixing between collectives (which could be akin to gamete production and zygote formation), and inclusion of more than two kinds of particle types. More complex selective regimes can also be envisaged, such as those that reward collectives based on absolute population size of particle types, which would allow less abstract collective function to be considered. However, regardless of these refinements, we suspect that our core conclusion will stand firm: collectivelevel selection favours particle dynamics that improve collectivelevel heredity. The ability to reliably reestablish successful adult states of pastgenerations from simpler and potentially noisy initial conditions is adaptive at the collective level.”
https://doi.org/10.7554/eLife.53433.sa2Article and author information
Author details
Funding
Agence Nationale de la Recherche (ANR10IDEX00102)
 Guilhem Doulcier
 Paul B Rainey
Agence Nationale de la Recherche (ANR10LABX54)
 Silvia De Monte
Agence Nationale de la Recherche (ANR11IDEX000102)
 Silvia De Monte
MaxPlanckGesellschaft
 Paul B Rainey
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Wenying Shou, Sara Mitri and Alvaro Sanchez for review of the manuscript and valuable comments. We also thank Benjamin Kerr, Will Ratcliff, Eric Libby, Michael Doebeli, Pierrick Bourrat and Philippe Nghe as well as the members of our respective teams for their feedback. GD was funded by the Origines et Conditions d’Apparition de la Vie (OCAV) programme, PSL University (ANR10IDEX001–02), awarded to PBR, SDM and AL; SDM was supported by the French Government under the program Investissements d’Avenir (ANR10LABX54 MEMOLIFE and ANR11IDEX000102PSL), AL thanks the Center for Interdisciplinary Research in Biology (CIRB) for funding.
Senior Editor
 Patricia J Wittkopp, University of Michigan, United States
Reviewing Editor
 Wenying Shou, Fred Hutchinson Cancer Research Center, United States
Reviewers
 Wenying Shou, Fred Hutchinson Cancer Research Center, United States
 Sara Mitri, University of Lausanne, Switzerland
 Alvaro Sanchez, Yale University, United States
Publication history
 Received: November 7, 2019
 Accepted: June 12, 2020
 Accepted Manuscript published: July 7, 2020 (version 1)
 Version of Record published: August 20, 2020 (version 2)
Copyright
© 2020, Doulcier et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 2,765
 Page views

 436
 Downloads

 21
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Ecology
 Evolutionary Biology
Circadian clocks infer time of day by integrating information from cyclic environmental factors called zeitgebers, including light and temperature. Single zeitgebers entrain circadian rhythms, but few studies have addressed how multiple, simultaneous zeitgeber cycles interact to affect clock behavior. Misalignment between zeitgebers (‘sensory conflict’) can disrupt circadian rhythms, or alternatively clocks may privilege information from one zeitgeber over another. Here, we show that temperature cycles modulate circadian locomotor rhythms in Nematostella vectensis, a model system for cnidarian circadian biology. We conduct behavioral experiments across a comprehensive range of light and temperature cycles and find that Nematostella’s circadian behavior is disrupted by chronic misalignment between light and temperature, which involves disruption of the endogenous clock itself rather than a simple masking effect. Sensory conflict also disrupts the rhythmic transcriptome, with numerous genes losing rhythmic expression. However, many metabolic genes remained rhythmic and inphase with temperature, and other genes even gained rhythmicity, implying that some rhythmic metabolic processes persist even when behavior is disrupted. Our results show that a cnidarian clock relies on information from light and temperature, rather than prioritizing one signal over the other. Although we identify limits to the clock’s ability to integrate conflicting sensory information, there is also a surprising robustness of behavioral and transcriptional rhythmicity.

 Ecology
 Evolutionary Biology
Division of labor, where subpopulations perform complementary tasks simultaneously within an assembly, characterizes major evolutionary transitions of cooperation in certain cases. Currently, the mechanism and significance of mediating the interaction between different cell types during the division of labor, remain largely unknown. Here, we investigated the molecular mechanism and ecological function of a policing system for optimizing the division of labor in Bacillus velezensis SQR9. During biofilm formation, cells differentiated into the extracellular matrix (ECM)producers and cheaterlike nonproducers. ECMproducers were also active in the biosynthesis of genomic islandgoverned toxic bacillunoic acids (BAs) and selfresistance; while the nonproducers were sensitive to this antibiotic and could be partially eliminated. Spo0A was identified to be the coregulator for triggering both ECM production and BAs synthesis/immunity. Besides its wellknown regulation of ECM secretion, Spo0A activates acetylCoA carboxylase to produce malonylCoA, which is essential for BAs biosynthesis, thereby stimulating BAs production and selfimmunity. Finally, the policing system not only excluded ECMnonproducing cheaterlike individuals but also improved the production of other public goods such as protease and siderophore, consequently, enhancing the population stability and ecological fitness under stress conditions and in the rhizosphere. This study provides insights into our understanding of the maintenance and evolution of microbial cooperation.