Group-selection via aggregative propagule-formation enables cooperative multicellularity in an individual based, spatial model

The emergence of multicellularity is one of the major transitions in evolution that happened multiple times independently. During aggregative multicellularity, genetically potentially unrelated lineages cooperate to form transient multicellular groups. Unlike clonal multicellularity, aggregative multicellular organisms do not rely on kin selection instead other mechanisms maintain cooperation against cheater phenotypes that benefit from cooperators but do not contribute to groups. Spatiality with limited diffusion can facilitate group selection, as interactions among individuals are restricted to local neighbourhoods only. Selection for larger size (e.g. avoiding predation) may facilitate the emergence of aggregation, though it is unknown, whether and how much role such selection played during the evolution of aggregative multicellularity. We have investigated the effect of spatiality and the necessity of predation on the stability of aggregative multicellularity via individual-based modelling on the ecological timescale. We have examined whether aggregation facilitates the survival of cooperators in a temporally heterogeneous environment against cheaters, where only a subset of the population is allowed to periodically colonize a new, resource-rich habitat. Cooperators constitutively produce adhesive molecules to promote aggregation and propagule-formation while cheaters spare this expense to grow faster but cannot aggregate on their own, hence depending on cooperators for long-term survival. We have compared different population-level reproduction modes with and without individual selection (predation) to evaluate the different hypotheses. In a temporally homogeneous environment without propagule-based colonization, cheaters always win. Predation can benefit cooperators, but it is not enough to maintain the necessary cooperator amount in successive dispersals, either randomly or by fragmentation. Aggregation-based propagation however can ensure the adequate ratio of cooperators-to-cheaters in the propagule and is sufficient to do so even without predation. Spatiality combined with temporal heterogeneity helps cooperators via group selection, thus facilitating aggregative multicellularity. External stress selecting for larger size (e.g. predation) may facilitate aggregation, however, according to our results, it is neither necessary nor sufficient for aggregative multicellularity to be maintained when there is effective group-selection.


Statistical analysis of division count
We have determined the number of cell division events (  ) between successive colonization events.We then compared the   for the various colonization mechanisms (A-F) and with and without predation (Figure A).In case of no colonization (A), we measured the number of divisions between successive resource replenishment events (i.e.whenever the initial amount of resources in the lattice goes below 40%, it is replenished).For the different mechanism, see Fig 2 of the main text.We have conducted the statistical analysis with R, using the tidyverse, ggpubr, gridExtra, dunn.test and outliers packages.respectively).Each box represents 140 independent simulations.The  axis is on a logarithmic scale.The inset shows the zoomed-in region of interest.Points denote outliers.There is significant difference, though it has practically no effect on the absolute number of cell divisions between the different propagation modes.There is a difference between the cases when predation is or is not in effect, but the difference is still negligible overall.Under predation, cells divide faster on average, because due to the increased mortality, there is always more empty space to divide to.The parameters are as in Table 1 in the main text.
First, we determined whether the data series are normally distributed or not.We have used three methods: asymptotic one-sample Kolmogorov-Smirnov test, Shapiro-Wilk normality test and have visually analyzed the Q-Q plots.In order to perform the Shapiro-Wilk normality test in R, we had to reduce our sample size randomly to 5000 data points.Neither the dataset with, nor the one without predation is normally distributed.Due to the non-normal distribution of our data, we have used Kruskal-Wallis rank sum test to compare whether there is a significant difference between the cases with and without predation to see if predation influences the number of divisions or not.We have compared the datasets of with and without predation, ignoring the background variable (colonization type).We have found significant difference in almost all cases, except the (F, D) pair ( 2 = 927.37, = 1,  < 2.2 • 10 −16 ), corresponding with Figure A.
To compare the different colonization modes (A-F), we have used Wilcoxon-Mann-Whitney test with continuity correction (Table A) and Dunn test (Table B).In these tests, the effect (or lack of) predation was the background variable that we have ignored.We have also compared colonization modes using Dunn's test.The results indicate that in most cases, the difference is significant ( 2 = 14826.0887, = 5,  < 2.2 ⋅ 10 −10 ).For detailed comparison of groups, see Table B.Although the data series differ significantly, there is no difference in order of magnitude between the different colonization mechanisms.In addition, pairs of random refugium and random fragmentation (C, D), random refugium and aggregation-based propagule formation (C, F) and aggregation-based dispersion and propagule formation (E, F) are the most similar (demonstrated by the Wilcoxon-Mann-Whitney and the Dunn tests; in case of Dunn, there was no significant difference between C and D).Random dispersion (B) and random refugium (C) are also similar, although there is some difference.Since there is no significant difference in absolute magnitude, these differences can be disregarded having no effect on our main results.

Statistical analysis of colonization times
Besides the number of cell divisions, we have also determined the elapsed time (number of updates) between two successive colonization events (colonization time, Δ), for the different colonization mechanisms (A-F), and with/without predation (Figure C).We performed Kruskal-Wallis rank sum test to compare the Δ of the different colonization mechanisms.We performed Wilcoxon-Mann-Whitney test to compare the effect of predation for each colonization mechanism.
In most cases, colonization times significantly differ from each other depending on transfer mechanism (A-F).However, since there is no significant difference in absolute magnitude, these differences can be disregarded having no effect on our main results.In case of no colonization (A), we have measured the elapsed time between two successive resourcereplenishment events (occurring when the amount of food in the lattice decreases below 40% of the initial food quantity).Similarly as before, we have tested whether the data is normally distributed or not.We have used three methods: asymptotic one-sample Kolmogorov-Smirnov test, Shapiro-Wilk normality test and have visually analyzed the Q-Q plots.In order to perform the Shapiro-Wilk normality test in R, we had to reduce our sample size randomly to 5000 data points.Neither the dataset with, nor the one without predation is normally distributed.Due to the non-normal distribution of our data, we have used Kruskal-Wallis rank sum test to compare whether there is a significant difference between the cases with and without predation to see if predation influences the number of divisions or not.We have compared the datasets with and without predation, ignoring the background variable (colonization type).Results indicate that the mechanisms significantly differ from each other ( 2 = 5431.37, = 1,  < 2.2 • 10 −16 ).
To compare the different colonization modes (A-F), we have used Wilcoxon-Mann-Whitney test with continuity correction (Table C).

Table C. Results (p values) of the Wilcoxon-Mann-Whitney test, comparing the colonization times for the different colonization mechanisms (A-F). Cells with grey shading indicate insignificant values. Otherwise, the test indicates that most colonization mechanisms produce different colonization times
.
We have also compared colonization modes using Dunn's test ( 2 = 22917.4949, = 5,  = 0), see Table D.While data series may differ significantly, there is no difference in order of magnitude between the different colonization mechanisms.Regarding the elapsed time between colonization events, random refugium differs mostly from the rest.Since there is no significant difference in absolute magnitude, these differences can be disregarded having no effect on our main results.   1 in the main text, with the following difference:  = 150 000 000.Only in case of random fragmentation and aggregation-based propagule formation manifests a considerable increase in the survival probability of cooperators, indicating that in every other case, the threshold for cooperator survival is either the minimum value (at least a single cooperator initially) or that there should be no defectors at all (100% cooperator).Note, that the last column shows <100%, as the initial population contains a few defectors.The parameters are as in Table 1 in the main text, with the following difference:  = 150 000 000.In case of random refugium (C), the benefit of association is negligible at a low association probability (C5), and only dissociated cells divide.Since defectors are better competitors when free, they divide more and ultimately win over cooperators.At medium association and high dissociation rates (C6), associated cooperators divide more, often leaving the aggregate, and the resulting dissociated cooperators will further divide, winning over defectors (C7).However, if the dissociation rate is smaller than the association rate, cooperators cannot leave the aggregate, hindering effective division in dense aggregates (C7, C8).In case of aggregation-based propagule formation (F), at low dissociation rates (F9) defectors easily associate at the edge of aggregates, preventing neighbouring cooperators to effectively divide.If the dissociation probability is high enough, cells more easily leave aggregates and thus can divide into empty space more often (F10).This is more pronounced, as there is predation and defectors at the edge die more frequently.

Figure A .
Figure A. Average number of cell divisions in the population between two successive colonization events for the different colonization modes (A-F), with and without predation (size-dependent selection; red and teal boxes, respectively).Each box represents 140 independent simulations.The  axis is on a logarithmic scale.The inset shows the zoomed-in region of interest.Points denote outliers.There is significant difference, though it has practically no effect on the absolute number of cell divisions between the different propagation modes.There is a difference between the cases when predation is or is not in effect, but the difference is still negligible overall.Under predation, cells divide faster on average, because due to the increased mortality, there is always more empty space to divide to.The parameters are as in Table1in the main text.

Figure B .
Figure B. Visual analysis of the distribution of the number of cell divisions between successive colonization events plotted against the theoretical normal distribution.The points denote our sample.The straight line with the grey envelope shows where our data points would be if they were normally distributed.A: data without predation, B: data when predation was in effect.

Figure
Figure C. Average time elapsed between successive colonization events for the different colonization mechanisms (A-F), with and without predation (size-dependent selection; red and teal boxes, respectively).Each box represents 140 independent simulations.The parameters are as in Table1in the main text.

Figure D .
Figure D. Visual analysis of the distribution of the number of cell divisions between successive colonization events plotted against the theoretical normal distribution.The points denote our sample.The straight line with the grey envelope shows where our data points would be if they were normally distributed.A: data without predation, B: data when predation was in effect.

Figure E .
Figure E. The effect of size-dependent selection strength on cooperator survival.Individual plots show the change in the relative frequency of cooperators over time according to the strength of size dependent selection (predation; columns) and propagation mechanism (rows).Each plot shows the time-evolution of 140 independent simulations (coloured lines).The dashed line represents the mean of the relative frequency of cooperators over all simulations, the dotted lines enveloping the grey area represent the mean +/-standard deviation over all 140 simulations.The bar charts give a visual cue of the distribution of independent simulations, showing the terminal ratio of cooperators of individual trajectories categorized into five discrete groups (see legend). Figure F shows only the bars, with scales, for better comparison.The selection strength increases nonlinearly with decreasing  to the right (as per Eq. 1 in main text, see Figure K).As the strength of predation decreases, the survival probability of cooperators decreases too, as defectors are less and less stressed by predators.This effect is most pronounced in the case of aggregationbased propagule formation, where increasing predation quickly gives an advantage to cooperators.The parameters are as in Table1in the main text, with the following difference:  = 150 000 000.

Figure F .
Figure F. Relative frequency of cooperators in the last generation of the simulations of Figure E, depending on propagation mechanism (A-F) and predation strength  (columns, increasing to the right).The  axis denotes the relative frequency of cooperators in the last generation.Top row: Terminal frequencies are categorized into five bins, as per the legend.Bottom row: Box plots show the distribution of terminal cooperator frequencies.Propagation mechanisms are: A: no colonization, B: random dispersion, C: random refuge, D: random fragmentation, E: aggregation-based dispersion, F: aggregation-based propagule formation.The parameters are in Table1in the main text, with the following difference:  = 150 000 000.

Figure G .
Figure G.The effect of initial cooperator ratio on cooperator survival.Individual plots show the change in the relative frequency of cooperators over time, under predation stress, according to the initial ratio of cooperators (columns; note the nonlinear scaling of values) and propagation mechanism (rows A-F, see Figure 2 in main text for details).Predation is in effect ( = 5).Each plot shows 140 independent simulations (coloured lines).The dashed line represents the mean, dotted lines enveloping the grey area represent the mean +/-standard deviation over all 140 simulations.The bar charts give a visual cue of the distribution of independent simulations, showing the terminal ratio of cooperators of individual trajectories categorized into five discrete groups.Figure H shows only the bar charts with scales, for better comparison.The initial cooperator ratio increases to the right.As the initial ratio of cooperators increases, their survival probability increases too in all cases of propagation.Only in case of random fragmentation and aggregation-based propagule formation manifests a considerable increase in thesurvival probability of cooperators, indicating that in every other case, the threshold for cooperator survival is either the minimum value (at least a single cooperator initially) or that there should be no defectors at all (100% cooperator).Note, that the last column shows <100%, as the initial population contains a few defectors.The parameters are as in Table1in the main text, with the following difference:  = 150 000 000.

Figure H .
Figure H.Relative frequency of cooperators in the last generation of the simulations of Figure G, depending on propagation mechanism (A-F, see Fig. 2 in main text for details) and the initial ratio of cooperators (columns; note the nonlinear scaling of values).The  axis denotes the relative frequency of cooperators in the last generation.Top row: Terminal frequencies are categorized into five bins, as per the legend.Bottom row: Box plots show the distribution of terminal cooperator frequencies.The parameters are in Table1in the main text, with the following difference:  = 150 000 000.

Figure I .
Figure I. Relative frequency of cooperators at selected association and dissociation rates.Predation was in effect ( = 5).The six columns (1-6) correspond to specific association-dissociation rate pairs (  , , respectively).Rows A-F correspond to colonization mechanisms (see Fig. 2 in main text for details).Inset maps at the right reproduce the respective plots ofFig.4 of the main text indicating the positions of the numbered (  , ) pairs.

Figure J .
Figure J.The total number of cell divisions in case of selected propagation mechanisms and association/dissociation pairs (  , ), categorized by the cell type (bluish cooperator, reddish defector).Only those propagation mechanisms are shown, where there was an interesting pattern to explain in the respective heatmap of Fig. 4 of the main text (right).The coordinates for various association/dissociation probability pairs ( axis) are indicated by indexed points overlaid on the respective heatmaps.Predation was in effect ( = 5).The  axis represents the number of cell divisions.In the left column, all divisions are counted; in the right column, only those divisions are counted where the daughter cell occupied an empty space.Note, that in most cases, the boxes corresponding to dissociated cooperators and associated defectors are too thin to see their color.In case of no colonization (A), the advantage of cooperators is the result of predation stress strongly affecting defector mortality.As the association probability increases, the survival (and division) probability of cooperators also increases.Cooperators at smaller association probabilities mostly divide when dissociated (A1, A2, A3).Increasing the association rate above a threshold causes associated cooperators to divide more than dissociated ones, because most cooperators are now associated (A4).In case of random refugium (C), the benefit of association is negligible at a low association probability (C5), and only dissociated cells divide.Since defectors are better competitors when free, they divide more and ultimately win over cooperators.At medium association and high dissociation rates (C6), associated cooperators divide more, often leaving the aggregate, and the resulting dissociated cooperators will further divide, winning over defectors (C7).However, if the dissociation rate is smaller than the association rate, cooperators cannot leave the aggregate, hindering effective division in dense aggregates (C7, C8).In case of aggregation-based propagule formation (F), at low dissociation rates (F9) defectors easily associate at the edge of aggregates, preventing neighbouring cooperators to effectively divide.If the dissociation probability is high enough, cells more easily leave aggregates and thus can divide into empty space more often (F10).This is more pronounced, as there is predation and defectors at the edge die more frequently.

Figure K .
Figure K.The probability of death from size dependent selection (predation) depending on the number of associated neighbours ( axis) and the strength of predation  (differently colored bars), according to Eq. 1 of the main text.Note the logarithmic  scale.The smaller  is, the higher the selection strength will be.

Table D .
Results (p values) of Dunn's test for comparing the colonization times for the different colonization mechanisms (A-F).Cells with grey shading indicate non-significant values.Dunn's test:  2 = 22917.4949, = 5,  = 0.

Table E .
The life cycles emerged in the computer simulations of [1], according to model parameters.