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).

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.

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).

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 < 0.05; Fig. 4B); if the time gap was greater than ~17,000 time steps, or the prior occupied space was larger than 129 grid boxes, the first emerged AFCP would win the competition (Fig. 4C, D). Together with the lineage analysis, these results confirmed that the order of emergence by the different genotypes would largely decide the formation of different interdependent patterns in the final communities, i.e., the earlier an AFCP emerged, the greater the chance for it to dominate the equilibrium communities.

Next, we investigated the role of spatial positioning of new genotypes. AFCPs that emerge first would have an opportunity to occupy more space which would favor its winning, but what if different AFCPs evolved simultaneously? Previous results indicated that if there was just a very short time gap between the emergence of the first and second TFLG, the advantage of prior space occupancy was not significant (Fig. 4C). In our simulations, the occasional fortuitous proximity of a focal OFLG cell to its functional complementary partner cell may also benefit the formation of the interdependent spatial aggregate. Therefore, higher association degree at the beginning of self-organization might facilitate space occupation, and help the AFCP gain the ‘region size advantage’ when the segregated spatial pattern is initially established. To test this hypothesis, we defined ‘partner association degree’ (PAD) and ‘partner association index’ (PAI) in a manner similar to a previous study [29], in order to quantify the different relative association degrees of AFCPs. PAD is the average number of the weighted functional complementary partner cell in the immediate neighborhood of a focal TFLG cell (Supplementary Fig. 3A). PAI is the ratio of PADs of the two AFCPs, which quantify the relative PAD between them. For example, PAI_001:010 equals PAD₀₀₁/PAD₀₁₀, which quantify the relative PAD between pair [0, 0, 1] & [1, 1, 0], and pair [0, 1, 0] & [1, 0, 1], while PAI_001:010 > 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 < 0.05). This analysis demonstrated that the random spatial proximity between the positions of emergence for two members of an AFCP was important in determining the competition outcome.

To further prove this idea, we conducted additional simulations. Four genotypes, [0, 0, 1], [1, 1, 0], [0, 1, 0], and [1, 0, 1], were inoculated with the same initial number. We randomly distributed the genotypes in the 100 × 100 grids to obtain diverse initial spatial distribution (see Methods; Supplementary Fig. 3B). We then selected a series of distributions with a gradient of relative association degrees of AFCPs (PAI_001:010), and ran simulations initialized with those values. The spatial self-organization also occurred in these cases, accompanied by the increase of PAD₀₀₁ and PAD₀₁₀ (Fig. 5B; Supplementary Fig. 4). For a given initial distribution, we conducted 100 replicate simulations to calculate the frequency of winning. A strong positive relationship between initial association degree of an AFCP and its winning frequency was found (Fig. 5B, C, P < 0.01), which further confirmed the important role of relative positioning on deciding the ecological outcome.

One might ask, in absence of emergence order and relative positioning, whether other random events impacted the outcomes of spatial competition? To address this issue, we simulated communities starting from a symmetric initial distribution of cells, and two types of AFCPs, [0, 0, 1] & [1, 1, 0] and [0, 1, 0] & [1, 0, 1], were equally inoculated (Supplementary Fig. 5A). Eliminating initial spatial asymmetry and distinct space occupancy, the segregated interdependent spatial pattern still emerged, and one of the AFCPs dominated the final community in each run (Supplementary Fig. 5B). However, the frequencies of winning for the two AFCPs were nearly identical in the 100 repeated simulations (Supplementary Fig. 5C), suggesting other random events did not change the relative winning probability of the two AFCPs.

Interdependent pattern of community is associated with increased stability

Do the communities dominated by functional complementary LOF genotypes exhibit distinct community properties, e.g., different community stability? Our mathematical framework allowed us to make comparisons among communities at different time points. To investigate whether community stability changes with the succession, we imposed nutrient disturbances on the communities at initial and final time points. We stopped the nutrient supply for 12 h, and after restoring it, explored community response to the disturbance. Although all communities declined after nutrient depletion, the resistance of the evolved communities to disturbances was better than the original community composed of only the autonomous population, suggesting the evolved communities containing interdependent patterns were more stable (Fig. 6A). The improved community stability was attributed to the more effective resource allocation of the LOF genotypes. As shown in Fig. 6B, although the average available substrate for each cell was similar for both scenarios during the disturbance (solid line in the shaded region), the amount of resources required for producing public goods significantly decreased when an interdependent pattern formed (dashed line), indicating more resources can be assigned to growth or maintenance when the community is facing environment perturbations. This suggests that after the interdependent pattern evolved, the resources that were originally wastefully allocated to produce redundant public goods, were saved to fight against the harsh environmental change. This result can explain why the interdependent pattern is selectively favored at the community level.

Source: Ecology - nature.com