Selfishness driving reductive evolution shapes interdependent patterns in spatially structured microbial communities
The logic of the model
Our spatially-resolved model was simulated in discrete grid boxes of a 100 × 100 array, which included four basic assumptions: (1) Initial individuals were assumed to secrete three public goods but may randomly mutate to lose any of those functions with a certain probability; (2) Secreting a public good created a corresponding metabolic burden, therefore in losing a function the individual would gain a benefit; (3) All public goods were essential for growth. The net growth rates of individuals were dependent on the local concentrations of public goods; (4) Substrate and public goods diffused between two grid boxes at rates proportional to the concentration gradient.
For the 1st assumption, we included three functions because it is the minimal unit and tersest design to simulate complex communities, allows for the emergence of three categories of interaction patterns, and a single cooperative LOF genotype might evolve from differential evolutionary paths (Fig. 1A, B). The genotypes were described by bit strings containing 1 and 0 which indicated the genotype could produce the corresponding public good or not, respectively. Eight genotypes could emerge during the simulations, which were the initial autonomous producer [1, 1, 1], three one-function loss genotypes (OFLGs, i.e., [1, 1, 0], [1, 0, 1], and [0, 1, 1]), three two-function loss genotypes (TFLGs, i.e., [1, 0, 0], [0, 1, 0], and [0, 0, 1]), and a nonproducing cheater [0, 0, 0] (Fig. 1A).
Fig. 1: Logic of the individual-based model.
A Possible genotypes and evolutionary relationships among them emerging from reductive evolution when starting with an autonomous genotype that performs three essential public functions. Note that in this three-function model, some genotypes, i.e., Two-function loss genotypes and cheaters, might evolve from different mother genotypes. B Interaction patterns that could possibly be established in the spatially structured communities. C Schematic of the individual-based simulations. A 100 × 100 array initialization with all autonomous phenotypic individuals (left) was conducted with a long-term stepwise iteration to investigate if diverse interaction patterns could form (right). At each time step, calculations were done from the level of individual grids (top) to whole lattice (bottom). Within each grid box, Monod equation modified by basic assumptions of the Black Queen Hypothesis was used to calculate the microbial growth, while minimum and maximum thresholds of biomass were defined to decide the division and death of individuals (top middle). Microbial individuals were allowed to randomly mutate to lose functions (top middle). Classical discretization of the diffusion equation gave local rules for updating the concentrations of public goods and nutrients in each box (middle). State changes at the individual level lead to the evolutionary dynamics of the communities, which may give rise to the formation of diverse interaction patterns (bottom).
Full size image
The 2nd and 3rd assumptions were developed from the basic mathematical assumption of the BQH [19], and defined individual growth by integrating the benefit and cost of function loss (Fig. 1C). To conceptualize the cost of performing a function, we supposed a parameter (α) which is the fraction of biomass used to produce a public good per unit time of an individual. In addition, we defined a second parameter (β) as the ratio of the amount of public goods required during each step to account for the produced public goods. Therefore the redundant fraction of public goods production was 1−βj, and lower βj reflected a higher amount of redundant public goods that could be gained from the producers by the LOF genotypes, resulting in decreased risk in association with function loss (see Supporting Information S1 for more details). During the model simulation, spatiotemporal dynamic variables, i.e., positioning of genotypes and the time points at which genotypes evolved, would be collected. We initiated the simulations by randomly distributing 100 ancestor cells [1, 1, 1] into the grid boxes and iterated for at least 1,500,000 time steps. During each time step, individuals grew, decayed, reproduced, and mutated according to the previously mentioned assumptions (Fig. 1C). We paid attention to whether stable communities with various interdependent patterns could be formed after a specified number of iterations, as well as recorded the spatiotemporal dynamics of the communities.
Diverse interdependent patterns emerged with high level of function cost and varied level of functional redundancy
For model simulations, the function cost (parameter α) and functional redundancy (parameter β) were assigned to 0.0001, 0.0005, 0.001, and 0.4, 0.6, 0.8, respectively. A total of 2891 independent simulations with 9 parameter sets displayed different community structures (Fig. 2A). When the function cost was assigned to a low level, i.e., 0.0001, the autonomous ancestor dominated the community. When function costs were assigned to higher levels, 0.0005 and 0.001, new genotypes evolved and later interacted to form three distinct types of interdependent patterns even within the same α and β combination, i.e., asymmetric functional complementation (AFC), complete functional division pattern, and one-way dependency, with the relative amounts of 1677/2891, 143/2891, and 48/2891, respectively. In addition, higher functional redundancies favored the loss of more functions, increasing complexity of the community structures.
Fig. 2: Reductive evolution shapes diverse interdependent patterns in microbial communities.
A The final (steady state) community structures across gradients of function cost (α) and functional redundancy (1-β). Results were summarized from at least 300 interdependent runs for each parameter set. Community structures were assessed after simulation for 170,000 iterations, where 98.9% (2891/2923) of runs reached steady state. According to the structures, replicates were clustered into several scenarios for each parameter set, which are shown separately in the area plots. Note that the values of β is the proportion of public goods that is required for growth, and thus 1-β reflects the level of function redundancy. B, C Six representative community dynamics on the spatial lattices were selected from one interdependent simulation with the given conditions (mut = 10−5, α = 0.001, β = 0.8), showing the evolution of three types of asymmetric functional complementary pairs (AFCPs) (B), three different paths for the evolution of pairs [0, 0, 1] & [1, 1, 0] (C). Left images indicate the distribution of different genotypes at different points in evolutionary time. Curve plot in the middle describes the community dynamics of the corresponding simulation. Schematics at right briefly summarize the spatiotemporal dynamics of each simulation: the arrays in (B) indicate one type of AFCP directly dominated the communities without competition from others; the boxes in (C) indicates the composition of ancestor or AFCP in the related time points, while the windows inside indicate the spatial coexistence of multiple AFCPs and the size of the windows represents the relative fraction of different AFCPs.
Full size image
Among the three possible kinds of interactions, the AFC pattern was the most widespread, which was the combination of a two-function-loss genotype (TFLG) and its complementary one-function-loss genotype (OFLG). For example, [0, 0, 1], which produced a single essential public good, depended on its functional complement one-function-loss partner [1, 1, 0], for the other two public goods. Specifically, three types of the asymmetric functional complementary pairs (AFCPs), that is, [0, 0, 1] coupled with [1, 1, 0], [0, 1, 0] coupled with [1, 0, 1], and [1, 0, 0] coupled with [0, 1, 1], colonized most of the grid with a similar frequency of emergence. Interestingly, under the condition of high level of cost, the emergence of AFC patterns was accompanied by some nonproducing cheaters, whose relative abundance rose with the increase in functional redundancy (Fig. 2A top row). The addition of cheaters significantly reduced the total biomass of the communities, suggesting that high functional redundancy favors the evolution of cheaters which may decrease the community productivity. In addition, function loss happened more easily with high function cost. As the function cost parameter α increased from 0.0005 to 0.001, relative abundance of TFLGs increased approximately from 55 to 70% (Fig. 2A).
Besides the AFC patterns, two additional types of interdependent patterns evolved at a relatively lower frequency. The complete functional division pattern, that is, coexistence of [0, 0, 1], [0, 1, 0], and [1, 0, 0], only evolved when both factors were at high levels (α = 0.001, β = 0.4) with a frequency of approximately 45% (143 of 319 simulations, Fig. 2A, top right), which described a scenario with high benefit and low cost of function loss, favoring the loss of more functions and consequently more likely to maintain the evolution of TFLPs. Another form of interactions that emerged was one-way dependency, where one partner performs all functions and other none (i.e., coexistence of [1, 1, 1] and [0, 0, 0]). This form emerged at a low frequency (48 out of all 2891 simulations shown in Fig. 2A), but evolved with a higher probability under the condition of a mid-level function cost and low level of functional redundancy (α = 0.0005, β = 0.6, Fig. 2A, middle left), where the extinction of [1, 1, 1] was ~2.5 times slower than in other scenarios (Supplementary Fig. 1), leading to a higher potential for the spatial proximity between [1, 1, 1] and [0, 0, 0] during evolution.
Taken together, these phenomena demonstrated that the mutualistic exchange of complementary functions happened only when function cost was high. The emergence of different interdependent interaction patterns was related to the function cost and function redundancy, especially for the complete functional division and one-way dependency pattern, which only emerged within a limited parameter range. However, even for a given combination of α and β, it still remained possible for the evolution of distinct interaction patterns, suggesting that stochastic processes may play a role.
Same interdependent patterns might evolve via different modes
Because the evolution of three kinds of AFCPs were the most common scenarios in our simulations, we then focused on the role of stochastic processes, i.e., the key random events, in deciding the winning complementary pair among the three similar but different AFCPs. As a first step, we traced the variation in the spatiotemporal dynamics, trying to cluster the numerous evolutionary dynamics into limited modes and divide the complex evolutionary courses into several stages. These simplifications would facilitate the search for key random events.
Therefore, we analyzed the dynamics of 296 simulations with a typical parameter set (α = 0.001, β = 0.8), because under this condition, only the three types of AFCPs evolved, with a similar frequency of emergence (Fig. 2A, Top left), in order to avoid interference from the other interaction patterns. As described above, any of the three types of AFCPs could potentially take over the final community under this condition (Fig. 2B; Supplementary video 1–3). Using the emergence of AFCP [0, 0, 1] & [1, 1, 0] as an example, three categories of dynamic modes could give rise to its final domination. (1) After pair [0, 0, 1] & [1, 1, 0] emerged and formed a spatial aggregation, it rapidly expanded and took over the entire grid (Fig. 2C, first line; Supplementary video 3). (2) In addition to the pair [0, 0, 1] & [1, 1, 0], spatial aggregations of another AFCP also emerged (e.g., pair [0, 1, 0] & [1, 0, 1] in Fig. 2C, second line and Supplementary video 4). In this scenario, a special spatial pattern was established in a short period after the evolution of both AFCPs e, where pairs of two complementary members exhibited strong spatial mixing, while the two different AFCPs were totally segregated. Community succession was then governed by spatial competition between the two AFCPs. If pair [0, 0, 1] & [1, 1, 0] won the competition, it would dominate the final community. (3) Spatial aggregations of all three AFCPs emerged, and then pair [0, 0, 1] & [1, 1, 0] dominated the community after outcompeting the other two AFCPs (Fig. 2C, third line; Supplementary video 5). The clustering of these three possible modes of AFC patterns was also shown by the temporal dynamics of the α-diversity across different parameter sets (Supplementary Fig. 2), where the evolution modes of the AFC patterns were clearly clustered into three possible categories, suggesting that this clustering is independent of the determined factors α and β.
In sum, the succession of interdependent patterns could be divided into two stages: (1) the emergence of spatial aggregations composed of two interdependent members with strong connections; (2) spatial competition among different aggregations drive the community to evolve to the final state, composed of only one type of interdependent interactions. Of course, if only one type of AFCP emerged, the spatial competition stage would be unnecessary during succession.
Evolutionary random events play important roles in deciding the dominant AFCP in equilibrium communities
The presence of two evolutionary stages lead us to hypothesize that the random events affecting ecological outcomes should arise from two aspects. First, in the initial evolutionary stage, the emergence of interdependent spatial aggregations should be related to the order in which new genotypes emerge. Second, the outcome of the spatial competition should be also influenced by the initial positioning of the new genotypes.
The fact that each TFLP had two possible evolutionary paths (e.g., [1, 0, 0] could inherit its function from [1, 1, 0] or [1, 0, 1]), suggested that the effects of the random order of emergence for different genotypes were highly correlated with the evolutionary lineage. Therefore, to investigate the effects of this, we analyzed the evolutionary lineage of emergence, colonization, and loss of every genotype within the 296 simulations with the typical parameter set (α = 0.001, β = 0.8). In total, there were 24 evolutionary branches leading to the evolution of the three forms of AFC patterns (8 for each, Fig. 3). Among all these branches, we summarized two key random events (Fig. 3, red and blue boxes).
Fig. 3: The evolutionary trajectories of 296 independent simulations with the typical parameter set (mut = 10−5, α = 0.001, β = 0.8).
We analyzed the evolutionary trajectories of every interdependent run and clustered them into 24 types of branches (top, see Methods). The area plot shows the final community structures and the frequencies of each branch (bottom). Blue dashed box shows the evolutionary diversification into four scenarios after the first key event occurs, while the red dashed box indicates the 24 different evolutionary trajectories that diverged after the second key event occurs. Solid boxes with colored circles represent the genotypic composition of communities at different evolutionary time points. Red arrows indicate the branches where one type of asymmetric functional complementary pair (AFCP) directly dominated the communities without competition with other AFCPs, while the blue arrows indicate the branches where one type of AFCP took over the entire space after competitions with other AFCPs. Dashed boxes at the figure labels (right) indicate different AFCPs.
Full size image
The first event occurred after two types of OFLGs emerged. After this evolutionary time point, all three public functions were included in OFLGs. With the benefit of the function loss, these two OFLGs would expand and gradually outcompete the autonomous genotype [1, 1, 1]. Thus, the first key event was whether all three OFLGs could emerge before the autonomous genotype entirely disappeared (Fig. 3, blue box). If not, the third type of AFCP would never evolve; if so, all three types of AFCPs would still have a chance to dominate the final community. In the 296 simulations, the frequencies of these two scenarios were nearly same, that is, 147 simulations were clustered to the former, while 149 simulations were clustered to the latter. The 147 simulations, where the third type of AFCP never evolved, could be then divided into three categories with similar frequencies, where two of the three OFLGs occupied the whole space and excluded the ancestral population.
The second key evolutionary event was the emergence of TFLGs (Fig. 3, red box). After the two or three types of OFLGs successfully colonized, whose functional complementary TFLGs first to emerge in the next evolutionary time would lead to the prior formation of the spatial aggregation of the AFCP. It is obvious that if no other AFCP aggregations formed later, this AFCP would dominate the final community (Fig. 3, red arrow indicated branches). Alternatively, if other AFCP aggregations formed during the expansion process, the spatial competition between different AFCPs would decide the dominant AFCP in the equilibrium communities (Fig. 3, blue arrow indicated branches). In our analysis, the chance of only one AFCP evolving reached 64.7% (198 of the 296 simulations). If only two OFLGs evolved after the first event, the frequency of only one AFCP evolving reached 79.6% (121 of the 152 simulations). In contrast, if three OFLGs evolved after the first event, there could be a relative higher possibility of two or three AFCPs evolving (47.4%), meaning that spatial competition could then be an important process.
What decided the winner of the competition? We observed that after the segregated interdependent spatial pattern was newly established, the relative region sizes occupied by different AFCPs were the key to determining the winner (Fig. 2C, the second and third lines; Supplementary video 4 and 5). We analyzed the time gaps between the emergence of the two AFCPs in the second categories of succession modes and the size of the regions they occupied (Fig. 4A). The result indicated a significantly positive correlation between the length of the time gaps and the region size the prior AFCP occupied (t-test, p 1 indicates pair [0, 0, 1] & [1, 1, 0] is more spatially associated than pair [0, 1, 0] & [1, 0, 1]. Applying these definitions, the simulation results where the advantage of prior space occupancy was not significant (left side of blue line in Fig. 4C, 33 replicates) were selected for analysis, and we found a significantly positive correlation between the relative PAD at the beginning of spatial self-organization and the ‘region size advantage’ (Fig. 5A; p 1 means the prior emerged AFCPs are more spatially associated than the second AFCPs. Red dots indicate the first to emerge AFCP won the competition in the corresponding replicate, while the green dots indicate the second to emerge AFCP won the competition. B Two typical examples of simulations initialized with premixing the two types of AFCPs, [0, 0, 1] & [1, 1, 0] and [0, 1, 0] & [1, 0, 1], which represent scenarios when initial PAI001:010 1, respectively. C The significant positive correlation between the winning frequency of pair [0, 0, 1] & [1, 1, 0] and the initial value of PAI001:010. When initial PAI001:010 > 1, final communities were more likely to be dominated by pair [0, 0, 1] & [1, 1, 0], oppositely, pair [0, 1, 0] & [1, 0, 1] were more favorable when PAI001:010 More
