Cell size and selection for stress-induced cell fusion in unicellular eukaryotes

In unicellular organisms, sexual reproduction typically begins with the fusion of two cells (plasmogamy) followed by the fusion of their two haploid nuclei (karyogamy) and finally meiosis. Most work on the evolution of sexual reproduction focuses on the benefits of the genetic recombination that takes place during meiosis. However, the selection pressures that may have driven the early evolution of binary cell fusion, which sets the stage for the evolution of karyogamy by bringing nuclei together in the same cell, have seen less attention. In this paper we develop a model for the coevolution of cell size and binary cell fusion rate. The model assumes that larger cells experience a survival advantage from their larger cytoplasmic volume. We find that under favourable environmental conditions, populations can evolve to produce larger cells that undergo obligate binary cell fission. However, under challenging environmental conditions, populations can evolve to subsequently produce smaller cells under binary cell fission that nevertheless retain a survival advantage by fusing with other cells. The model thus parsimoniously recaptures the empirical observation that sexual reproduction is typically triggered by adverse environmental conditions in many unicellular eukaryotes and draws conceptual links to the literature on the evolution of multicellularity. August 16, 2024 1/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint Author summary Sexual reproduction is commonly observed, both in eukaryotic microorganisms and in higher multicellular organisms. Sex has evolved despite numerous apparent costs, including investment in finding a partner and the energetic requirements of sexual reproduction. Binary cell fusion is a process that sets the stage for sexual reproduction by bringing nuclei from different cells into contact. Here, we provide a mathematical explanation of the advantage conferred by binary cell fusion due to increased cell mass. We show that when unicellular organisms have the option to invest in either cell fusion or cell mass, they can evolve to fuse together as rapidly as possible in the face of adverse environments, instead of increasing their mass. These results are consistent with the empirical observation that sexual reproduction is often triggered by environmental stress in unicellular eukaryotes. Our results imply an advantage to cell fusion, which helps to shed light on the early evolution of sexual reproduction itself.

1 Although the details of the early evolution of sexual reproduction in the last common 2 eukaryotic common ancestor (LECA) are shrouded in mystery, it is argued that the 3 emergence of eukaryotic sex began with the evolution of cell–cell fusion and meiosis [1] 4 in an archaeal ancestor [2,3]. This step can be further broken down into the evolution of 5 binary cell fusion, the one spindle apparatus, homologous pairing and chiasma, and 6 finally reduction, division and syngamy [4]. The vast majority of theoretical studies 7 investigating the evolution of sexual reproduction have focused on later stages of this 8 evolutionary trajectory, namely the conditions that give rise to a selective pressure for 9 genetic recombination [5–8]. However, comparatively few studies have investigated the 10 selective pressures that may have first given rise to binary cell fusion, which may have 11 facilitated the evolution of a host of other eukaryotic traits [9], including the homologous 12 pairing and meiotic recombination, by bringing nuclei together in the same cell. 13 Hypotheses for the evolution of binary cell fusion often rely on hybrid fitness 14 advantage. It has been suggested that selection for cell–cell fusions might have initially 15 been driven by “selfish” transposons and plasmids [10–12], or negative epistatic 16 August 16, 2024 2/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint interactions between mitochondrial mutations [13,14]. However, once a heterokaryotic 17 cell has been formed (binucleate with nuclei from both parental cells), the advantage of 18 hybrid vigor and the masking of deleterious mutations could lead to the maintenance of 19 cell fusion [4]. Such benefits are required to alleviate costs to cell-fusion, which include 20 selfish extra-genomic elements in the cytoplasm [15] and cytoplasmic conflict [16,17]. 21 In these previous studies on the evolution of binary cell fusion, the effect of changing 22 environmental conditions is not considered. However, in many extant unicellular 23 organisms, binary cell fusion (and the karyogamy and genetic recombination that follow) 24 occur in response to challenging environmental conditions [18] such as starvation 25 (Tetrahymena [19,20]) and depleted nitrogen levels ( Chlamydomonas reinhardtii [21] 26 and Saccharomyces pombe [22]). Meanwhile in benign conditions with abundant 27 resources these species reproduce asexually using binary cell fission. The mechanisms 28 that drive selection for genetic recombination under challenging environmental 29 conditions are well-studied [23]; recombination can facilitate adaptation to a novel 30 environment [24,25] and evolving to engage in more sex when fitness is low 31 (fitness-associated sex) can allow an organism to maximise the advantages of sex while 32 minimising the costs [26,27]. However, this focus on the benefits of recombination leaves 33 space to ask whether binary cell fusion itself could be selected for as a stress response, 34 even in the absence of any genetic advantages. 35 In this paper, we do not account for the genetic factors discussed above. Instead, we 36 focus only on how the survival advantage associated with increasing cytoplasmic volume 37 might select for binary cell fusion; this relies on the physiological advantages conferred 38 by cell-cell fusion and is independent of the question of the genetic advantages (and 39 disadvantages) of sexual reproduction. This alternative perspective offers useful new 40 insights that can be compared with empirical observation. 41 That size-based processes could play a role in the early evolution of sexual 42 reproduction has empirical and theoretical support. The “food hypothesis” [28] suggests 43 that metabolic uptake could drive horizontal gene transfer in bacteria and archaea, with 44 DNA molecules providing nutrients for the receiving cell [29,30]. Indeed, horizontal 45 gene transfer has been shown experimentally to be an important source of carbon and 46 nutrients in bacteria [31,32]. Binary cell fusion is possible in bacteria (where it has been 47 shown to come with selective benefits from mixed cytoplasm [33]) and archaea [34,35]. 48 August 16, 2024 3/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint Meanwhile amongst eukaryotes, the benefits of increasing cytoplasmic volume are 49 understood to be strong enough to drive selection for the sexes themselves [36,37]. That 50 early selection for syngamy may have been driven by the survival benefits of larger 51 cytoplasmic volume (and the concomitant increase in solid food reserves) has been 52 argued verbally [38], however has not been explored mathematically. In suggesting a 53 mechanistic hypothesis for the evolution of binary cell fusion, our work has interesting 54 parallels with [39], where an advantage to cell fusion is identified in terms of shortening 55 the cell-cycle. 56 Moving to consider potential physiological benefits of binary cell fusion naturally 57 leads to work on the evolution of multicellularity. While multicellularity and binary cell 58 fusion are clearly biologically distinct, from a modelling perspective they share 59 similarities in that they can involve the “coming together” of cells to produce a larger 60 complex [40]. Multicellularity achieved via aggregation allows organisms to rapidly 61 adapt to novel environments that favour increased size [41,42]. It has also been 62 suggested that the genetic nonuniformity of such aggregates may also make them 63 well-suited to resource limited environments [43,44], echoing the hybrid vigor 64 hypotheses for the evolution of early syngamy [4]. Relatively few theoretical studies 65 have investigated the evolution of facultative aggregation in response to changing 66 environments [42]. However in the context of clonal multicellularity (“staying together”) 67 such environments have been considered more extensively [45,46]. In this clonal context, 68 the evolutionary dynamics act primarily on fragmentation modes [47,48] (e.g. how a 69 “parental” multicellular complex divides to form new progeny). Interestingly the same 70 quality-quantity trade-off arises here [45] as drives selection for the sexes [36] 71 (anisogamy, gametes of differing sizes); larger daughter cells (or gametes) are more able 72 to withstand unfavourable environmental conditions, while smaller cells can be 73 produced in larger quantities. 74 In this paper we adapt the classic Parker-Baker-Smith [36] (PBS [37]) model for the 75 evolution of sexes in order to investigate the evolution of binary cell fusion. This builds 76 on recent work that investigates how the possibility of parthenogenetic reproduction can 77 drive selection for oogamy in eukaryotes [49]. We assume for simplicity that parental 78 cells undergo a number of cell-divisions. The size of daughter cells is a compound 79 evolvable trait determined both by the size of parental cells and the number of cell 80 August 16, 2024 4/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint divisions. Daughter cells are then introduced to a pool in which they can undergo 81 binary cell fusion, with larger fused cells experiencing a survival advantage. Unlike in 82 the classic PBS model, the fusion rate is also a trait subject to evolution. In the 83 following sections we proceed to outline some insights developed from numerical 84 simulations before going on to develop analytical results for the model in a fixed 85 environment. Finally we introduce switching environments and show that under plastic 86 phenotypic responses, facultative binary cell fusion in response to harsh environmental 87 conditions can evolve. 88 Model 89 Insights from Simulations 90 We consider a computational model of a haploid population that reproduces via binary 91 cell fission. This population dynamics proceed as follows: there is an initial growth and 92 binary fission phase where a fixed energy budget E can be used for population growth 93 and binary fission (see Fig 1A). The parameter E thus sets the carrying capacity of the 94 population. This growth phase is followed by an environmentally-induced mortality 95 phase where growth and fission are suspended and survival depends on cell size. 96 Explicitly, we assume that cells reach maturity at size M, so that in the absence of 97 binary fission ( n = 0 rounds of fission) the total number of mature cells at the end of a 98 growth cycle is E/M (see Fig 1A). If, however, the population undergoes n > 0 rounds 99 of binary fission, then each resulting daughter cell has mass m = M/2n and the total 100 number of daughter cells is (2 nE)/M. After this growth and fission phase each 101 daughter cell is subject to an extrinsic mass-dependent mortality, such that larger 102 daughter cells are more likely to survive into the next growth cycle (see Fig 1A). We 103 denote this survival function S(m; β); the parameter β describes the magnitude of the 104 mortality process (i.e. the harshness of the environment). For a given value of m, an 105 increase in β decreases the survival probability. 106 In an initial investigation, we view the mass of daughter cells, m, as a single trait 107 subject to evolution. To increase the number of daughter cells they produce ( E/m) 108 mature cells can grow to smaller sizes (reduced M, which increases the number of 109 August 16, 2024 5/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint Division Survival Growth Fusion Possibility of fusion failure Growth / Veg. Seg. (b) Coevolutionary dynamics of cell mass and cell fusion (a) Evolutionary dynamics of cell mass Fig 1. Schematic for the model dynamics within each growth cycle. Panel (a): Illustrative dynamics for the evolutionary dynamics of cell mass alone. Due to energetic constraints genotypes in the population can either produce fewer, larger mature cells or more numerous, smaller cells (see different shades of green). Daughter cells are produced following cell division. Their survival is dependent on mass, such that smaller cells are more likely to die (see Eq. (2)). Surviving cells seed the next growth cycle. Panel (b): Illustrative dynamics for the coevolutionary dynamics of cell mass and cell fusion rate. The model is similar to that in panel A, but now a fraction of daughter cells are given the opportunity to risk fusing to form binucleated cells; with probability C fusion fails, and both daughter cells are lost. However should a fused cell successfully form, it experiences an enhanced survival probability as a result of its larger cytoplasmic volume. Following growth and vegetative segregation, surviving cells seed the next growth cycle. mature cells in the population) or increase their number of cell-divisions (increased n) . 110 However by decreasing M and increasing n, individuals also produce smaller daughter 111 cells that are more vulnerable to extrinsic mortality. The size of daughter cells is thus 112 subject to a quality-quantity trade-off. For simplicity, we model the mass of daughter 113 cells m as a continuous trait, and explore its evolution using simulations (see 114 Supporting Information 3.3). 115 Fig 2 summarises the outcome of such evolutionary dynamics. Fig 2A shows that the 116 population evolves towards an evolutionarily stable strategy (ESS) in m for a given 117 environment. Should the environment suddenly become harsher (via an increase in β) 118 the population evolves towards a new ESS, in which daughter cells are larger (i.e. 119 daughter cells evolve to become larger to withstand more adverse conditions). 120 We now modify the model to allow for the possibility of binary cell fusion following 121 the cell fission described above. Daughter cells may now fuse to form a binucleated cell 122 (e.g. a dikaryon [50], in which the cytoplasm of the contributing cells are mixed but 123 August 16, 2024 6/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint their nuclei or nucleoids remain distinct [38]) or remain a mononucleated cell (with a 124 single nucleus or nucleoid). The rate of cell fusion is given by α, such that when α = 0 125 all cells remain mononucleated, and cell survival into the next growth cycle is calculated 126 as before. Conversely for α > 0, some proportion of daughter cells will have fused and 127 in the limit α → ∞, all cells will have fused. These fused cells will receive a survival 128 advantage from their increased mass. However they will also pay an additional cost, C, 129 resulting from factors such as cell-fusion failure [51], selfish extra-genomic elements in 130 the cytoplasm [15], cytoplasmic conflict [16,17] and maintenance of a binucleated 131 cell [52]. Together this means that fused cells survive with a total probability 132 (1 − C)S(2m; β). Surviving adults divide to form a new growth cycle of mononucleated 133 haploid daughter cells, with binucleated parental cells producing mononucleated 134 progeny through vegetative segregation [53] (or alternatively through plasmid 135 segregation machinery [54]). Note that although we do not account for the possibility of 136 binucleated cells failing to form mononucleated progeny (i.e. failed segregation), this 137 can be accounted for by their additional survival cost, C (see Fig 1B). 138 We now explore the coevolution of daughter cell mass, m, and fusion rate, α. In 139 Fig 2B, we see that in the benign environment, α remains at zero, and the population 140 evolves towards an ESS in m as in Fig 2A. However now when the population is 141 introduced to a harsher environment, the evolutionary dynamics differ from those in 142 Fig 2B (where α was held artificially at zero). Rather than cells evolving to be larger, 143 we see a different response emerging; selection for binary cell fusion ( α > 0). 144 The result above is in some sense surprising. Despite the presence of additional 145 survival costs associated with binary cell-fusion, selection for non-zero fusion rates 146 (rather than increased daughter cell size) persists in the harsh environment. We explain 147 the emergence of this behaviour mathematically in the Results section. 148 Mathematical Model 149 Our model takes inspiration from the classic PBS model for the evolution of 150 anisogamy [36] (the production of sex cells of differing size). However, whereas such 151 models typically consider the binary cell fusion (fertilization) rate a fixed parameter, we 152 here treat it as a trait subject to evolution. In doing so our work builds on [49], where a 153 August 16, 2024 7/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint 0 500 1000 1500 2000 0.0 0.5 1.0 1.5 2.0 Tim e ( Generat ions,τ) Cell m ass ( m ) A Evolut ionary dynam ics of cell m ass Benign Environm ent Harsh Environm ent 0.0 0.5 1.0 1.5 2.0Cell m ass ( m ) Coevolut ionary dynam ics of cell m ass and cell fusion Benign Environm ent Harsh Environm ent 0 500 1000 1500 2000 0.0 0.2 0.4 0.6 0.8 1.0 Tim e ( Generat ions,τ) Fusion Rat e (α) Benign Environm ent Harsh Environm ent B Fig 2. Selection for cell fusion as an alternative to increased cell size in response to a harsh environment . Stochastic simulations of evolutionary trajectories when the system is subject to a switch from the benign environment (β1 = 0.5, green region) to the harsh environment (β 2 = 2.2, orange region) at growth cycle 500. Panel A illustrates the case where the fusion rate is held at α = 0, representing the scenario where the physiological machinery for fusion has not evolved. Panel B illustrates the case where fusion rate is also subject to evolution. Remaining model and simulation parameters are given in Supporting Information 7 and the initial condition is ( m(0), α(0)) = (1.16, 0). very similar model with a different biological motivation was used to investigate the 154 evolution of anisogamy with parthenogenesis. In order to analyse the dynamics of the 155 model, we use tools from adaptive dynamics [55], assuming that traits are continuous 156 and that mutations have small effect. 157 In addition we will explore the effect of switching environments, another departure 158 from the PBS model. As such it is important to keep track of the hierarchy of 159 August 16, 2024 8/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint timescales at play. The shortest timescale is the timescale of a growth cycle (see 160 Supporting Information 1). The intermediate timescale is that over which the invasion 161 of a rare mutant (taking place over many growth cycles) can occur (see Supporting 162 Information 2.1) and the longest timescale is the evolutionary timescale, representing 163 the cumulative effect of multiple mutations and invasions (see Supporting 164 Information 3.1). Finally we assume that environmental switching can take place on 165 either intermediate or long timescales (see Supporting Information 5 and [49]). 166 Dynamics within each growth cycle 167 A total of (2 nE)/M daughter cells enter a pool in which binary cell fusion can occur. 168 After a finite time window, the resultant cells are subject to a round of mass dependent 169 mortality, such that cells of larger mass are more likely to survive. The surviving cells 170 form the basis of the next growth cycle, as illustrated in Fig 1. 171 F usion Kinetics 172 We assume that all daughter cells may fuse with each other, an assumption consistent 173 with most models of the early evolution of sexual reproduction, which suppose the 174 existence of a “unisexual” early ancestor that mated indiscriminately [56]. Following 175 growth and binary cell fission, the population is comprised of N unfused daughter cells. 176 Fusion between these mononucleated cells occurs at a rate of α, such that the number of 177 unfused cells, N, is given by the solution to 178 dN dt = −αN 2 , N (0) = 2nE M = E m . (1) At the end of the fusion window of duration T there are then N(T ) unfused 179 (mononucleated) cells remaining, and (N (0) − N(T ))/2 fused (binucleated) cells. 180 Survival Probability 181 We assume that both unfused and fused cells are subject to the same extrinsic 182 mass-dependent mortality function, S(m; β), while fused cells pay an additional 183 mass-independent cost C. Many choices for such a function are possible, so long as it is 184 an increasing function of cell size (which we equivalently refer to as cell mass m). 185 August 16, 2024 9/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint However here we assume that S(m; β) is the Vance survival function [57], a common 186 assumption in the literature [58 –60]. We thus have that at the end of the fusion window, 187 the survival probability of unfused and fused cells are given respectively by 188 S(mi; β) = exp  − β mi  , (1 − C)S(mj + mk; β) = exp  − β mj + mk  , (2) where mi is the mass of a particular unfused daughter cell and mj and mk are the 189 masses of two daughter cells that have fused. For a given cell mass, increasing β will 190 decrease the survival probability. We therefore refer to β as the environmental 191 harshness parameter, with high β corresponding to harsh environments in which 192 survival is difficult, and low β corresponding to more benign environments in which 193 even cells of modest mass have a high probability of surviving. 194 Having defined how the survival of a cell depends on its mass, we have the necessary 195 tools to mathematically characterise the fitness of a rare mutant, and whether it can 196 invade the resident population. In the following section, we provide mathematical 197 approximations of the invasion dynamics of such a mutant. 198 Invasion Dynamics 199 Adopting the classical assumptions of adaptive dynamics [55,61] (see also Supporting 200 Information 3.1 and [49]), we mathematically approximate the invasion dynamics of a 201 mutant (which occur over discrete growth cycles). Deriving these invasion dynamics 202 analytically is only possible when we assume that mutations in m and α occur 203 independently. However the evolutionary dynamics we obtain if we consider mutations 204 occurring in both m and α simultaneously remains identical to those obtained by 205 assuming that they occur independently (see Supporting Information 4 and [49]). 206 Denoting by ˆfm the frequency of mutants of size m ± δm in the population where 207 δm is the mutational stepsize in m, which is assumed to be small and tg the number of 208 growth cycles, we find (see Supporting Information 2) 209 d ˆfm dtg = hm(m, α, β, C) ˆfm(1 − ˆfm) , (3) August 16, 2024 10/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint where hm(m, α, β, C) is a constant that depends on the parameters m, α, β, C (see 210 Eq. (S10)). This constant provides the fitness gradient of the mutant. Similarly, 211 denoting by ˆfα the frequency of mutants with fusion rate α ± δα in the population 212 where δα is the mutational stepsize in α, which is again assumed to be small, we find 213 d ˆfα dtg = hα(m, α, β, C) ˆfα(1 − ˆfα) (4) where hα(m, α, β, C) is the fitness gradient of a mutant with fusion rate α + δα. We see 214 that in the case of a single mutant, we have frequency-independent selection for mutants 215 with different masses and fusion rates. We note that in reality, frequency-dependent 216 invasion dynamics can occur when multiple mutants that change both m and α arise in 217 the population (see Supporting Information 4 for mathematical analysis), which can in 218 turn lead to evolutionary branching [62,63]. However, since this branching does not 219 occur in the regimes we are focusing on in this paper [49], we assume for simplicity that 220 mutants encounter a monomorphic resident population (trait substitution) for the 221 remainder of the mathematical analysis. 222 Evolutionary Dynamics 223 We assume that haploid daughter cells are characterised by two genetically 224 non-recombining traits mass m and cell fusion rate α. We assume that mutations occur 225 in m or α independently at a fixed rate µ, where µ is measured in units of 226 (number of growth cycles) −1 (see Supporting Information 3.3 and [49]). A mutation in 227 m represents a change in the mass of the daughter cell produced, and a mutation in α 228 represents a change in the fraction of the population that undertakes either one of the 229 reproductive routes (i.e. binary cell fusion vs strictly binary cell fission). 230 Mutants with a different mass to their ancestor can produce either more or fewer 231 daughter cells than their ancestor (see Eq. (1)), which impacts their survival (see 232 Eq. (2)). When mutants have a different fusion rate to their ancestor, although the 233 number of daughter cells produced does not differ from their ancestor, the number of 234 fused cells at the end of a growth cycle can either increase/decrease, which impacts their 235 survival, since fused cells have greater mass. The survival of fused cells is also 236 influenced by the cost of fusion C (see Eq. (2)). 237 August 16, 2024 11/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint Our mathematical analysis in the remainder of this section assumes that mutants 238 encounter a strictly monomorphic population (i.e. that mutations fixate before the 239

of a new mutant). However in our numerical simulations, we release this 240 restriction and stochastically allow for the coexistence of multiple traits in the 241 population, held under a mutation-selection balance, as described in the subsequent 242 section. 243 Fixed Environment 244 We first consider the evolutionary dynamics in the case where the environment is fixed 245 (i.e. when the parameter β, which measures the harshness of the environment (see 246 Eq. (2)), is constant throughout the evolution. Assuming that δm and δα are small 247 (small mutational step size), we use techniques from adaptive dynamics [49,64] to 248 obtain equations for the evolutionary dynamics of m and α, which are given by 249 dm dτ = Hm(m, α; β, C) = −4m(m − β) + EαT (1 − C)e β 2m (4m − β) 4m2(m + EαT (1 − C)e β 2m ) dα dτ = Hα(m, α; β, C) = −m(1 − (1 − C)e β 2m ) ln(1 + EαT m ) 2α(EαT (1 − C)e β 2m + m) (5) for α ≥ 0. This boundary is imposed to prevent α from becoming negative, which is 250 biologically unrealistic since it corresponds to an increase in daughter cell numbers 251 during the fusion period, as can be seen from Eq. (1). Therefore when α becomes 252 strictly decreasing along this boundary α = 0 boundary (i.e [ dα/dτ]|α=0 < 0 in Eq. (5)), 253 we introduce a discontinuous change in the dynamics, given by 254 dm dτ = Hm(m, 0; β, C) = β − m m2 dα dτ = 0 (6) The derivations of these equations can be found in Supporting Information 3.1. 255 August 16, 2024 12/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint Switching environments with phenotypic plasticity 256 Now, we consider the case where evolution acts on the same traits as before, but the 257 environment is subject to change. We model environmental change as switching 258 between two environments β1 and β2. If β1 > β 2, then β1 is the harsher environment 259 (see Eq. (2)). We also allow for phenotypic plasticity such that the population can 260 evolve different strategies in different environments. The population’s evolutionary state 261 is now described by four traits; the daughter cell mass in environments 1 and 2 ( m1 and 262 m2) and fusion rate in these environments (α 1 and α2). 263 For simplicity we assume that any cost of phenotypic switching or environmental 264 sensing is negligible and that this plastic switching is instantaneous upon detection of 265 the change in environmental conditions. The evolutionary dynamics in each environment 266 are then decoupled. However the evolutionary trajectories in each environment are 267 coupled by the initial trait values for the population in each environment, which we 268 assume are the same (i.e. the population begins in a phenotypically undifferentiated 269 state). With phenotypic plasticity, the evolutionary dynamics are then given by 270 dm1 dτ = P1Hm(m1, α1; β1, C), dm2 dτ = P2Hm(m2, α2; β2, C) dα1 dτ = P1Hα(m1, α1; β1, C), dα2 dτ = P2Hα(m2, α2; β2, C) (7) with initial conditions 271 m1(0) = m2(0) = m0, and α1(0) = α2(0) = α0 . (8) Here, Hm(m, α; β, C) and Hα(m, α; β, C) retain the functional form in Eq. (5). 272 As Eqs. (7) are only coupled through their shared initial conditions, m0 and α0, the 273 choice of these initial conditions is an important consideration. Since we are interested 274 in the initial evolution of binary cell fusion, it is natural to assume that the population 275 evolves from a state of zero fusion, α0 = 0. Deciding on a plausible initial daughter cell 276 mass takes more thought. One parsimonious choice would be that the population is 277 already adapted to either environment 1 or environment 2 and that m0 is given by an 278 evolutionary fixed point in one of these environments (this is the situation illustrated in 279 Fig 2). However if the population has been exposed to both the environments before 280 August 16, 2024 13/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint phenotypic placticity has evolved, it is possible that m0 is instead given by a bet-hedging 281 strategy. We explore that such a strategy would take in the following section. 282 Switching environments without phenotypic plasticity 283 We now consider the case where there is switching between environments but where the 284 population exhibits no phenotypic plasticity. As described in the previous section, we 285 are particularly concerned with the period before the physiological machinery for cell 286 fusion has evolved, and so focus on the case where the cell fusion rate is fixed to zero, 287 α = 0. Evolution then solely acts on the daughter cell mass, m. 288 As in [49], environmental switching is modelled as a discrete stochastic telegraph 289 process, with the time spent in each environment distributed geometrically. The 290 population spends an average of τ1 = 1/λ1→2 in environment 1 and τ1 = 1/λ2→1 in 291 environment 2, where λi→j is the transition rate from environment i to j. 292 The two switching rates most relevant to our model are when the environment 293 switches many times before an invasion can complete, (fast relative to invasion, FRTI) 294 and when each switching event occurs after multiple invasions have completed, (fast 295 relative to evolution, FRTE). More detail of these switching rates are provided in 296 Supporting Information 5. However in [49], we show that the evolutionary dynamics for 297 m in both these regimes can be approximated using the same dynamical equations. 298 Using adaptive dynamics techniques modified to account for such environmental 299 switching [65], we obtain 300 dm dτ = P1Hm(m, 0; β1, C) + P2Hm(m, 0; β2, C) (9) where Hm(m, α; β, C) retains the functional form in Eq. (5) and P1 = τ1/(τ1 + τ2) and 301 P2 = τ2/(τ1 + τ2) are the probabilities of finding the population in the two respective 302 environments. We therefore see that in the absence of phenotypic plasticity, the 303 evolutionary dynamics is the weighted average of the dynamics in the two environments. 304 Obtaining the ESS for Eq. (9) is relatively straightfoward. Substituting for 305 Hm(m, 0; β1, C) and Hm(m, 0; β2, C) using the functional form given in Eq. (5) and 306 August 16, 2024 14/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint setting dm/dτ = 0 in Eq. (9), we obtain the ESS 307 m∗ BH,α=0 = P1β1 + (1 − P1)β2 . (10) This strategy constitutes a bet-hedging strategy in cell mass when the population has 308 yet to evolve phenotypic plasticity nor the capacity for cell-cell fusion. In the limits 309 P1 → 1 and P2 → 1, we can recover the ESS strategies in the two respective 310 environments: 311 m∗ 1,α=0 = β1 , m ∗ 2,α=0 = β2 , (11) which can be verified from a consideration of the equations for d m/dτ = 0 in a fixed 312 environment with α = 0 (see Eq. (6)). We can now proceed to analyse how binary cell 313 fusion can be selected for when the fusion rate α is allowed to increase from zero in the 314 following Results section. 315 Implementation of Numerical Simulations 316 The stochastic simulations of the evolutionary trajectories are also implemented using a 317 Gillespie algorithm [66] where successive mutations and environmental switching events 318 occur randomly with geometrically distributed waiting times. The rates of mutations µ 319 and environmental switching λ are measured in units of (number of growth cycles) −1. 320 In the simulations, multiple traits coexist under a mutation-selection balance (see 321 Supporting Information and [49] and [67] for more detail), which allows us to account 322 for variations in selection strengths in simulations of our evolutionary trajectories. 323

324 In this section we proceed to analyse the evolutionary dynamics derived from the 325 mathematical model and compare our results to numerical simulations of the full 326 stochastic simulations. 327 August 16, 2024 15/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint In a fixed environment the population evolves to either no cell 328 fusion, or to high levels of cell-fusion, dependent on the cost of 329 cell fusion 330 In Fig 3, we see two possible evolutionary outcomes for the fusion rates in a fixed 331 environment; the population can evolve either a high (technically infinite) fusion rate or 332 to zero fusion rates. To which of these fusion rates the population is attracted depends 333 both on the parameters and the initial conditions. 334 When the costs to cell fusion are low ( C ⪅ 0.39), the only evolutionary fixed point is 335 the high fusion rate fixed point (see Fig 3A). In this scenario, obligate fusion is the only 336 evolutionary outcome. 337 For intermediate costs to cell fusion 0 .39 ⪅ C ⪅ 0.86, there are two possible 338 evolutionary outcomes. The outcome depends on the initial conditions (see Fig 3B). If 339 the initial mass on the α = 0 boundary is small, selection acts to increase fusion rate 340 and obligate fusion is the ESS. However, if the initial mass on the boundary is 341 sufficiently large, the state of no cell fusion becomes the evolutionarily stable state. 342 Finally when costs to fusion are extremely high ( C ⪆ 0.86, see Fig 3C), selection for 343 decreased fusion rate acts regardless of the initial value of m on the α = 0 boundary, 344 and a state in which α = 0 (zero fusion rate) is the only evolutionary outcome. Under 345 this high cost regime, dα/dt < 0 along the entire line α = 0 and so fusion rate is never 346 selected to increase given any initial daughter cell mass. 347 A rigorous mathematical analysis that formalises the arguments above are provided 348 in [49], which uses a similar model to investigate the evolution of anisogamy with 349 parthenogenesis. In summary, the set of possible evolutionary attractors ( m∗, α∗), 350 starting from an initial condition ( m(0), α(0)) = (m(0), 0), are given by 351 (m∗, α∗) →           (β/4, ∞) if 1 − 1 √e > C ≥ 0 (β/4, ∞) or ( β, 0) if 1 − 1 e2 > C ≥ 1 − 1√e (β, 0) if C ≥ 1 − 1 e2 (12) where we note 1 − e−1/2 ≈ 0.39 and 1 − e−2 ≈ 0.86. While intermediate costs 352 (0.86 ⪆ C ⪆ 0.39) lead to two potential evolutionary outcomes depending on the initial 353 August 16, 2024 16/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint 0.0 0.5 1.0 1.5 2.0 0.0 0.2 0.4 0.6 0.8 1.0 Daught . cell m ass ( m ) Fusion Rat e (α) Low Fusion Cost ( C≲0.39) 0.0 0.5 1.0 1.5 2.0 0.0 0.2 0.4 0.6 0.8 1.0 Daught . cell m ass(m) Fusion Rat e (α) I nt erm ediat e Fusion Cost ( 0.39≲C≲0.86) 0.0 0.5 1.0 1.5 2.0 0.0 0.2 0.4 0.6 0.8 1.0 Daught . cell mass(m) Fusion Rat e (α) High Fusion Cost ( 0.86≲C)A CB Fig 3. Phase portraits for the co-evolutionary dynamics in a fixed environment (see Eq. (5)). High fusion rates are the only evolutionary outcome when costs to cell fusion are low (panel A), while under intermediate costs (panel B), high fusion rate and zero fusion rate (obligate asex) are both evolutionary outcomes, and under high costs, zero fusion rate becomes the only evolutionary outcome (panel C), as summarised analytically in Eq. (12). The red shaded region shows trajectories leading to points on the α = 0 boundary for which evolution selects for decreasing fusion rate (d α/dτ < 0) and the critical point at which d α/dτ = 0 is marked by the red arrow (see Supporting Information 6.2). The red circles mark a fixed point in the evolutionary dynamics of m (m∗ = β, see Eq. (12)), which may be unstable (open circles) or stable (filled circle) under coevolution with α. The blue circles and arrows illustrate the evolutionary fixed point for high fusion rates ((m ∗, α∗) → (β/4, ∞), see Eq. (12)). Average population trait trajectories, (⟨m⟩(t), ⟨α⟩(t)), from simulation of the full stochastic model are plotted in light gray, and their mean over multiple realisations are dashed. Initial conditions: ( m(0), α(0)) = (1.5, 0.6) and ( m(0), α(0)) = (2, 0.1). Simulation is run for 1.1 × 107 growth cycles in panel A, 1.24 × 107 growth cycles in panel B and 10 7 growth cycles in panel C. Remaining system parameters are given in Supporting Information 7. conditions, it is the second of these, ( m∗, α∗) = (β, 0), that is arguably the most 354 relevant for the evolution of early cell fusion; if evolution had acted on daughter cell size, 355 m, before the physiological machinery necessary for cell fusion had evolved, the initial 356 condition for the co-evolutionary dynamics would be ( m(0), α(0)) = (β, 0), at which the 357 population would be subsequently held by costs to fusion. 358 In Fig 3 we also see that our mathematical analysis is a good predictor of the 359 outcome of stochastic simulations (gray shaded lines). One minor point of departure is 360 that at high fusion rates our simulated trajectories begin to diverge from our analytic 361 prediction. This discrepancy is the result of evolutionary branching in cell mass, which 362 we explore in another paper relating to the emergence of size dimorphism in sex 363 cells [49]. However this branching happens at a later evolutionary stage than the focus 364 of this study, the early emergence of binary cell fusion. 365 We conclude this section by addressing the key biological result that arises from this 366 August 16, 2024 17/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint analysis; cell fusion is uniformly selected for even under moderately high costs (with a 367 fraction of up to C ≈ 0.39 of fused cells failing to survive) and can even be selected for 368 under higher costs (up to C ≈ 0.86) given necessary initial conditions. In the context of 369 the evolution of early binary cell fusion, this provides a surprising nascent advantage to 370 cell fusion. This advantage could even help compensate for other short-term costs 371 arising from the later evolution of sex and recombination. The selective advantage 372 experienced by fusing cells comes from their increased cytoplasmic volume, which leads 373 to increased survival probabilities. 374 In a switching environment with phenotypic plasticity, binary 375 cell fusion can evolve as a facultative stress response to harsh 376 environments 377 Having considered the case of the evolutionary dynamics in a fixed environment, we now 378 move on to consider the evolutionary dynamics of a population exhibiting phenotypic 379 plasticity in a switching environment (see Eq. (7)). We recall that under the 380 assumptions of costless and immediate phenotypic switching, the dynamics of ( m1, α1) 381 and (m2, α2) are decoupled. The evolution of the traits in the respective environments 382 are coupled however through the initial conditions from which they evolve, which must 383 be the same (i.e. a phenotypically undifferentiated state). 384 We consider two parsimonious choices for these initial conditions, both beginning in 385 a state without fusion (α 1(0) = α2(0) = 0). In the first scenario, we assume that the 386 population has evolved to a stable non-fusing mass adapted to a single environment (see 387 Eq. (11)) such that m1(0) = m2(0) = m∗ 1,α=0 or m1(0) = m2(0) = m∗ 2,α=0. This is a 388 situation in which the alternate environment is in some sense novel and one to which 389 the population has not adapted. In the second scenario, we instead assume that the 390 population has evolved to a bet-hedging strategy in mass (optimising the mass of 391 daughter cells across the two environments) such that m1(0) = m2(0) = m∗ BH,α=0 (see 392 Eq. (10)). An illustrative phase portrait is shown in Fig 4. 393 We initially consider first scenario in which a population has initially evolved under 394 environment 1 to reach a stable state ( m1, α1) = (β1, 0) (see red disk and surrounding 395 purple circle, Fig 4A). The population is now exposed to a second, harsher environment 396 August 16, 2024 18/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint (β2 > β 1) and allowed to evolve a phenotypically plastic response to this new 397 environment. Starting from an initial state ( m2(0), α2(0)) = (β1, 0) (see purple circle, 398 Fig 4B), traits m2 and α2 can evolve to ( m2, α2) = (β2/4, ∞). 399 Similarly, we now consider second scenario in which the population has initially 400 evolved under the switching environments to reach a stable state bet-hedging strategy 401 in mass ( m1(0), α1(0)) = (m2(0), α2(0)) = (m∗ α=0, 0) (see orange circles in Fig 4). Upon 402 the evolution of phenotypic plasticity, in the benign environment we see the population 403 traits relax to a stable state ( m1, α1) = (β1, 0) (no cell fusion, see Fig 4A). However in 404 the harsh environment, we see that the bet-hedging strategy in mass becomes unstable, 405 and the population evolves towards traits ( m2, α2) = (β2/4, ∞) in environment 2 (i.e. 406 we observe selection for binary cell fusion, see Fig 4B). 407 0.0 0.5 1.0 1.5 2.0 0.0 0.2 0.4 0.6 0.8 1.0 Daught . cell m ass ( m1) Fusion Rat e (α1) Benign Environm ent 1 A 0.0 0.5 1.0 1.5 2.0 0.0 0.2 0.4 0.6 0.8 1.0 Daught . cell m ass ( m2) Fusion Rat e (α2) Harsh Environm ent 2 B Fig 4. Selection for cell fusion as an alternative to increased cell size in response to a harsh environment . Illustrative phase portrait for co-evolutionary dynamics of ( m1, α1, m2, α2) in a switching environment with phenotypic switching that exhibits facultative binary cell fusion. In both environment 1 (panel A) and environment 2 (panel B) the cost to cell fusion is C = 0.6, purple circles represent the initial condition (m 1(0), α1(0)) = (m2(0), α2(0)) = (β1, 0), and orange circles represent the initial condition (m 1(0), α1(0)) = (m2(0), α2(0)) = (m∗ α=0, 0), with m∗ α=0 taken from Eq. (10). Environmental parameters are β1 = 0.5 and β2 = 2.2 making environment 1 the more benign environment, in which the population typically spends a proportion P1 = 0.7 of its time. In both scenarios described above, we see the emergence of facultative binary 408 cell-fusion as a response to harsh environmental conditions that lower the survival 409 probability of daughter cells. However we note that this is only possible if there is an 410 appreciable increase in environmental harshness, β, between the environments. In Fig 5, 411 we summarise the key results over the β1 − β2 parameter plane. Here we assume that 412 the cost to cell fusion is intermediate (1 − e−2 > C > 1 − e−1/2, i.e. 0 .86 ⪆ C ⪆ 0.39) 413 such that there are regions on the boundary α = 0 at which increased fusion rates are 414 August 16, 2024 19/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint both selected for and against depending on the value of m (see Supporting 415 Information 6.2 and Eq. (12)); this restricts us to the more interesting parameter regime 416 in which different evolutionary outcomes are possible in each environment. 417 In Fig 5, we see that when one of the environments is not appreciably worse than the 418 other, binary cell fusion does not evolve in either environment. However when the 419 difference between the environments grows more substantial, it is possible to evolve cell 420 fusion in the harsher environment from initial condition 421 (m1(0), α1(0)) = (m2(0), α2(0)) = (βi, 0) (where the population has first evolved towards 422 the evolutionary optimum of the more benign environment). Finally when the difference 423 between the environments is extreme, it is also possible to evolve cell fusion in the 424 harsher environment from initial condition ( m1(0), α1(0)) = (m2(0), α2(0)) = (m∗ α=0, 0) 425 (where the population has first evolved towards a bet-hedging strategy in cell mass). 426 0 1 2 3 4 0 1 2 3 4 β1 β2 No Fusion in either environment Fusion evolves in environment given m(0)=m1,α=0 Fusion evolves in environment 2 given m(0)=m1,α=0 and m(0)=mBH,α=0 * * * 2 2 Fig 5. Regions in the β1 − β2 plane where binary cell fusion evolves as a stress-response to environment 2. The region plot is independent of E and T . Here, C = 0.5, P1 = 0.3 and the initial condition is ( m(0), 0). Since C > 1 − 1/√e (see Eq. (12)), fusion can only evolve in at most one of the two environments. In this case it is environment 2 where fusion can evolve since m(0) = m∗ 1,α=0 (see Eq. (11)). A numerical simulation to support this regionplot is shown in Fig S2.

427 The evolution of sexual reproduction and its consequences for the subsequent 428 evolutionary trajectory of populations is of general importance in biology [7,13,68]. In 429 this paper we have illustrated a reversal of the classic two fold cost of sex that appears 430 in organisms with distinct sexes [69]; in unisexual [56], unicellular organisms, binary cell 431 August 16, 2024 20/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint fusion can be selected for, even in the presence of substantial costs, due to a survival 432 benefit that comes from increased mass. These results allow us to quantitatively assess 433 the verbal hypothesis that syngamy evolved allow cells to store more food reserves and 434 thus increase their survival rate [38]. It is particularly interesting that the benefits 435 conferred to cell fusion through increased cytoplasmic mass are sufficient to withstand 436 remarkably high costs; “obligate sexuality” is the only evolutionary outcome with costs 437 equivalent to a loss of ∼ 39% cells that attempt to fuse, and remains a potential 438 outcome with costs of up to ∼ 86% of fused cells dying. 439 Perhaps most interesting is the case of switching environments with phenotypic 440 plasticity. Here we find under a broad set of biologically reasonable conditions (costs to 441 cell fusion equivalent to 39% − 86% additional mortality to fused cells and at least 442 moderate changes in environmental quality) that high fusion rates are selected for in 443 harsh environments and zero fusion rates are maintained in benign environments. This 444 behaviour parsimoniously recapitulates the empirically observed reproductive strategies 445 of numerous facultatively sexual species, including C. reinhardtii [21], S. pombe [22] and 446 D. discoideum [70]. This mechanism, under which cell fusion evolves to increase the 447 survival probability of daughter cells, provides a complementary perspective on the 448 frequent evolution of survival structures (resistant to environmental stress) that form 449 following the formation of a zygote. These include ascospores in fungi [71] and 450 zygote-specific stress-resistant stress wall in C. reinhardtii [72]. Note that such 451 correlations between sexual reproduction and the formation of survival structures are 452 not as easily explained under genetic explanations for the evolution of sexual 453 reproduction, where engaging in both behaviours at once constitutes a simultaneous 454 (and therefore potentially costly) change in genotype and temporal dislocation in 455 environment [73,74]. 456 The results above are particularly interesting in the case of the evolution of early 457 binary cell-fusion as a first step in the evolution of sexual reproduction. While most 458 studies focus on the genetic benefits of cell-fusion [75] (including a functionally-diploid 459 dikaryotic cell [4]), or the genetic benefits of mixed cytoplasm [13,14] (which can also 460 come with costs [15–17,40,76–78]), the mechanism at play here is purely physiological. 461 Yet, as addressed above, it naturally captures the empirical observation of binary 462 cell-fusion as response to challenging environmental conditions, a feature absent in these 463 August 16, 2024 21/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint earlier models. While the mechanism does not explain the evolution of sexual 464 reproduction and genetic recombination itself, it does provide a nascent advantage to 465 binary cell fusion that sets the stage for the evolution of sex by bringing nuclei from 466 different cells into contact for prolonged periods. The mechanism also shows the 467 potential to counter short-term costs associated with the initial formation of a 468 binucleated cell. In this way the mechanism could facilitate the transition from 469 horizontal gene transfer [79,80] to meiotic recombination, which is advantageous when 470 genome sizes increase in length [8]. Conceivably, if genetic recombination is beneficial 471 for myriad genetic reasons in the long-term [8,81], it would seem natural that it would 472 be instigated when the opportunity arises (i.e. when physiological survival mechanisms 473 bring nuclei into close contact). We note that it is obviously possible that the first 474 diploid cells arose by errors in endomitosis [82 –84] (essentially doubling the chromosome 475 number within a single cell) and that meiosis first evolved in this context. Such a 476 sequence of events is still compatible with our very general model, which can alleviate 477 short-term costs of sex such as the energy involved in finding a partner and undergoing 478 fusion [69]. In either scenario sexual reproduction may not be only a direct response to 479 environmental variability [85,86], but also to the correlated formation of a survival 480 structure. 481 More generally, it is interesting to note that the conditions for facultative sexuality 482 (e.g. harsh environmental conditions) broadly coincide with those for facultative 483 multicellularity in both bacteria and eukaryotes, with starvation triggering the 484 formation of fruiting bodies in myxobacteria [87,88] and flocking in yeast [89,90]. 485 Meanwhile in C. reinhardtii, the formation of multicellular palmelloids and aggregates 486 are an alternate stress response to sexual reproduction [91], as are the formation of 487 fruiting bodies in D. discoideum [92]. In this multicellular context, the sexual behaviour 488 of D. discoideum is particularly interesting, as once formed, the zygote attracts 489 hundreds of neighboring cells that are then cannibalised for the provision of a 490 macrocyst [93]. These various survival strategies are unified in our model as a 491 mechanism for the evolution of binary cell fusion. 492 One element absent from our model is the fusion of multiple cells, which is likely to 493 be selected for under the assumptions implicit in our model. There would clearly be an 494 upper-limit on the number of fusions selected for, arising from the likely multiplicative 495 August 16, 2024 22/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint effect of the fusion cost C. However in this context, it is interesting to note that one of 496 the hypotheses for the evolution of self-incompatible mating types is as a signal to 497 prevent the formation of polyploid cells [94]. Such a mechanism could also prevent the 498 formation of trikaryotic cells should the cost of multiple fusions be too great. Thus, the 499 model neatly preempts the second stage in models for the evolution of eukaryotic sex, 500 the regulation of cell–cell fusion [1]. 501 The selection pressure for binary cell fusion is dependent on the formation of a 502 binucleated cell with an increased survival benefit arising from its larger size. Although 503 it is reasonable to assume that a single mononuclear cell has lower total resource 504 requirements than a multicellular complex of the same size [95], we have not considered 505 the detailed energetics of the maintenance of two nuclei [9]. In these respects 506 incorporating dynamic energy budget theory into the model would be an important 507 next step [96] as it would provide a clear distinction between the survival benefits of 508 fused cells and unfused bicellular complexes. Within our modelling framework, these 509 two structures are broadly similar [40]. However as we have shown, increased cell-cell 510 attraction can be selected for even in the presence of large costs that one might expect 511 under binary cell fusion but not associate with the formation of a bicellular complex. 512 We have assumed for simplicity a simple cell division scheme; parental cells undergo 513 n rounds of symmetric division to produce 2 n daughter cells. In the context of 514 multicellularity, switching environments have been shown to promote binary 515 fragmentation [45]. However non-symmetric modes can be selected for [47] reflecting the 516 diverse modes of facultative multicellular life cycles observed in bacteria [97]. It would 517 be interesting to incorporate our results into models of cell division that account 518 explicitly for growth [98,99] to determine how these results for multicellular organisms 519 carry over to the unicellular scenario, and further how they may affect those we have 520 shown here. 521 Finally, we have not explicitly modelled any sources of cytoplasmic or genetic 522 conflict [100], which we have for simplicity included in the fusion cost C. Nevertheless, 523 social conflict does emerge in this model. In a recent paper we have shown how 524 evolutionary branching can arise, with some individuals producing fewer larger cells and 525 others producing more numerous but smaller daughter cells [49]. This branching is 526 driven by the same evolutionary forces that drive selection for anisogamy [101], in which 527 August 16, 2024 23/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint context this can be viewed as sexual conflict [102]. That social conflict should arise in 528 the formation of multicellular aggregates is well understood [40,88,103]. However these 529 models typically assume cells of fixed size [46,104]. Combining the insights derived from 530 the evolution of anisogamy literature with the theory developed in the multicellularity 531 literature represents another promising research direction. 532 As addressed above, trade-offs between cell fusion rate and mass [105], 533 cell-energetics, inbreeding, and the possibility of multiple cell-fusion events offer 534 interesting avenues to extend this analysis. In addition, we have not accounted for the 535 discrete nature of divisions leading to daughter cells, costs to phenotypic switching, 536 non-local trait mutations, or pre-existing mating types. More generally, extending our 537 mathematical approach leveraging adaptive dynamics to switching environments in 538 other facultatively sexual populations might prove particularly fruitful [106,107]. 539 In this paper we have adapted the classic PBS model [36] in two key ways; allowing 540 the fusion rate to evolve and subjecting the population to switching environments. In 541 doing so, we have shown its capacity to parsimoniously capture the evolution of obligate 542 binary cell fusion, obligate binary cell fission, and stressed induced binary cell fusion in 543 unicellular organisms. These results offer particularly interesting implications for the 544 evolution of binary cell-fusion as a precursor to sexual reproduction, as well as 545 suggesting common mechanistic links between the evolution of binary cell fusion and 546 multicellularity. Moreover, our analysis emphasises the importance of exploring the 547 coevolutionary dynamics of a range of evolutionary parameters, and of developing 548 computational and mathematical approaches to elucidate facultative sexual 549 reproduction. 550 Supporting information 551 S1 T ext. Supporting information.pdf 552 Acknowledgments 553 This work has made use of Viking high performance computing service at the University 554 of York. 555 August 16, 2024 24/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint Author Contributions 556 Conceptualization: Jonathan Pitchford, George Constable 557 Methodology: Xiaoyuan Liu, George Constable, Jonathan Pitchford 558 F ormal Analysis: Xiaoyuan Liu, George Constable 559 Software: Xiaoyuan Liu 560 Investigation: Xiaoyuan Liu 561 Supervision: George Constable, Jonathan Pitchford 562 W riting- Original draft preparation: Xiaoyuan Liu, George Constable 563 W riting- review & editing: Jonathan Pitchford, George Constable 564

1. Goodenough U, Heitman J. Origins of eukaryotic sexual reproduction. Cold Spring Harbor perspectives in biology. 2014;6(3):a016154. 2. Eme L, Spang A, Lombard J, Stairs CW, Ettema TJ. Archaea and the origin of eukaryotes. Nature Reviews Microbiology. 2017;15(12):711–723. 3. Eme L, Tamarit D, Caceres EF, Stairs C, de Anda V, Schon ME, et al. Inference and reconstruction of the heimdallarchaeial ancestry of eukaryotes. BioRxiv. 2023; p. 2023–03. 4. Maynard Smith J. The Evolution of Sex. Cambridge UK: Cambridge University Press; 1978. 5. Muller HJ. The relation of recombination to mutational advance. Mutation Research/Fundamental and Molecular Mechanisms of Mutagenesis. 1964;1(1):2–9. 6. Van Valen L. The red queen. The American Naturalist. 1977;111(980):809–810. 7. Hartfield M, Keightley PD. Current Hypotheses for the evolution of sex and recombination. Integr Zool. 2012;7(2):192–209. August 16, 2024 25/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint 8. Colnaghi M, Lane N, Pomiankowski A. Genome expansion in early eukaryotes drove the transition from lateral gene transfer to meiotic sex. Elife. 2020;9:e58873. 9. Lane N. Energetics and genetics across the prokaryote-eukaryote divide. Biology direct. 2011;6:1–31. 10. Rose MR. The contagion mechanism for the origin of sex. Journal of Theoretical Biology. 1983;101(1):137–146. 11. Hickey DA. Selfish DNA: a sexually-transmitted nuclear parasite. Genetics. 1982;101(3-4):519–531. 12. Hickey DA. Molecular symbionts and the evolution of sex. Journal of Heredity. 1993;84(5):410–414. 13. Radzvilavicius AL. Evolutionary dynamics of cytoplasmic segregation and fusion: mitochondrial mixing facilitated the evolution of sex at the origin of eukaryotes. Journal of Theoretical Biology. 2016;404:160–168. 14. Tilquin A, Christie JR, Kokko H. Mitochondrial complementation: a possible neglected factor behind early eukaryotic sex. Journal of Evolutionary Biology. 2018;31(8):1152–1164. 15. Harrison E, MacLean R, Koufopanou V, Burt A. Sex drives intracellular conflict in yeast. Journal of Evolutionary Biology. 2014;27(8):1757–1763. 16. Hurst LD, Hamilton WD. Cytoplasmic fusion and the nature of sexes. Proceedings of the Royal Society of London Series B: Biological Sciences. 1992;247(1320):189–194. 17. Hutson V, Law R. Four steps to two sexes. Proceedings of the Royal Society of London Series B: Biological Sciences. 1993;253(1336):43–51. 18. Speijer D. What can we infer about the origin of sex in early eukaryotes? Philosophical Transactions of the Royal Society B: Biological Sciences. 2016;371(1706):20150530. August 16, 2024 26/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint 19. Phadke S. Sex ratios, mate choice, and diversity in multi-sex systems of ciliates. University of Houston; 2011. 20. Umen JG. Genetics: Swinging ciliates’ seven sexes. Current Biology. 2013;23(11):R475–R477. 21. Bishop C. Life cycle control of Chlamydomonas reinhardtii. Genome Biology. 2003;4(1):1–2. 22. Sieber B, Coronas-Serna JM, Martin SG. A focus on yeast mating: From pheromone signaling to cell-cell fusion. In: Seminars in Cell & Developmental Biology. Elsevier; 2022. 23. Ram Y, Hadany L. Condition-dependent sex: who does it, when and why? Philosophical Transactions of the Royal Society B: Biological Sciences. 2016;371(1706):20150539. 24. Fisher RA. The genetical theory of natural selection. Clarendon; 1930. 25. Waxman D, Peck JR. Sex and adaptation in a changing environment. Genetics. 1999;153(2):1041–1053. 26. Hadany L, Otto SP. The evolution of condition-dependent sex in the face of high costs. Genetics. 2007;176(3):1713–1727. 27. Hadany L, Otto SP. Condition-dependent sex and the rate of adaptation. the american naturalist. 2009;174(S1):S71–S78. 28. Vos M, Didelot X. A comparison of homologous recombination rates in bacteria and archaea. The ISME journal. 2009;3(2):199–208. 29. Redfield R. Genes for breakfast: The have-your-cake and-eat-lt-too of bacterial transformation. Journal of Heredity. 1993;84(5):400–404. 30. Redfield RJ. Do bacteria have sex? Nature Reviews Genetics. 2001;2(8):634–639. 31. Palchevskiy V, Finkel SE. Escherichia coli competence gene homologs are essential for competitive fitness and the use of DNA as a nutrient. Journal of bacteriology. 2006;188(11):3902–3910. August 16, 2024 27/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint 32. Narra HP, Ochman H. Of what use is sex to bacteria? Current Biology. 2006;16(17):R705–R710. 33. Shitut S, Shen MJ, Claushuis B, Derks RJ, Giera M, Rozen D, et al. Generating heterokaryotic cells via bacterial cell-cell fusion. Microbiology Spectrum. 2022;10(4):e01693–22. 34. Kuwabara T, Minaba M, Ogi N, Kamekura M. Thermococcus celericrescens sp. nov., a fast-growing and cell-fusing hyperthermophilic archaeon from a deep-sea hydrothermal vent. International journal of systematic and evolutionary microbiology. 2007;57(3):437–443. 35. Naor A, Gophna U. Cell fusion and hybrids in Archaea: prospects for genome shuffling and accelerated strain development for biotechnology. Bioengineered. 2013;4(3):126–129. 36. Parker GA, Baker RR, Smith V. The origin and evolution of gamete dimorphism and the male-female phenomenon. Journal of Theoretical Biology. 1972;36(3):529–553. 37. Lehtonen J. The Legacy of Parker, Baker and Smith 1972: Gamete Competition, the Evolution of Anisogamy and Model Robustness. Cells. 2021;10(3):573. 38. Cavalier-Smith T. Origin of the cell nucleus, mitosis and sex: roles of intracellular coevolution. Biology Direct. 2010;5:1–78. 39. Mancebo Quintana J, Mancebo Quintana S. A Short-Term Advantage for Syngamy in the Origin of Eukaryotic Sex: Effects of Cell Fusion on Cell Cycle Duration and Other Effects Related to the Duration of the Cell Cycle—Relationship between Cell Growth Curve and the Optimal Size of the Species, and Circadian Cell Cycle in Photosynthetic Unicellular Organisms. International Journal of Evolutionary Biology. 2012;2012. 40. Tarnita CE, Taubes CH, Nowak MA. Evolutionary construction by staying together and coming together. Journal of theoretical biology. 2013;320:10–22. August 16, 2024 28/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint 41. Pentz JT, M´ arquez-Zacar´ ıas P, Bozdag GO, Burnetti A, Yunker PJ, Libby E, et al. Ecological advantages and evolutionary limitations of aggregative multicellular development. Current Biology. 2020;30(21):4155–4164. 42. Staps M, Van Gestel J, Tarnita CE. Emergence of diverse life cycles and life histories at the origin of multicellularity. Nature ecology & evolution. 2019;3(8):1197–1205. 43. Bonner JT. Cellular slime molds. vol. 2127. Princeton University Press; 2015. 44. Hamant O, Bhat R, Nanjundiah V, Newman SA. Does resource availability help determine the evolutionary route to multicellularity? Evolution & Development. 2019;21(3):115–119. 45. Pichugin Y, Park HJ, Traulsen A. Evolution of simple multicellular life cycles in dynamic environments. Journal of the Royal Society Interface. 2019;16(154):20190054. 46. Tang S, Pichugin Y, Hammerschmidt K. An environmentally induced multicellular life cycle of a unicellular cyanobacterium. Current Biology. 2023;33(4):764–769. 47. Pichugin Y, Traulsen A. Evolution of multicellular life cycles under costly fragmentation. PLoS computational biology. 2020;16(11):e1008406. 48. Gao Y, Pichugin Y, Gokhale CS, Traulsen A. Evolution of reproductive strategies in incipient multicellularity. Journal of the Royal Society Interface. 2022;19(188):20210716. 49. Liu X, Pitchford JW, Constable GWA. Adaptive dynamics, switching environments and the origin of the sexes. bioRxiv. 2024; p. 2024–03. 50. Kruzel EK, Hull CM. Establishing an unusual cell type: how to make a dikaryon. Current opinion in microbiology. 2010;13(6):706–711. 51. Hall AE, Rose MD. Cell fusion in yeast is negatively regulated by components of the cell wall integrity pathway. Molecular biology of the cell. 2019;30(4):441–452. August 16, 2024 29/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint 52. Warr J. A mutant of Chlamydomonas reinhardii with abnormal cell division. Microbiology. 1968;52(2):243–251. 53. Rodilla V. Origin and evolution of binucleated cells and binucleated cells with micronuclei in cisplatin-treated CHO cultures. Mutation Research/Genetic Toxicology. 1993;300(3-4):281–291. 54. Ebersbach G, Gerdes K. Plasmid segregation mechanisms. Annu Rev Genet. 2005;39(1):453–479. 55. Br¨ annstr¨ om˚A, Johansson J, Von Festenberg N. The hitchhiker’s guide to adaptive dynamics. Games. 2013;4(3):304–328. 56. Heitman J. Evolution of sexual reproduction: A view from the fungal kingdom supports an evolutionary epoch with sex before sexes. Fungal Biology Reviews. 2015;29(3):108 – 117. doi:https://doi.org/10.1016/j.fbr.2015.08.002. 57. Vance RR. On reproductive strategies in marine benthic invertebrates. The American Naturalist. 1973;107(955):339–352. 58. Levitan DR. Optimal egg size in marine invertebrates: theory and phylogenetic analysis of the critical relationship between egg size and development time in echinoids. The American Naturalist. 2000;156(2):175–192. 59. Bulmer MG, Parker GA. The evolution of anisogamy: a game-theoretic approach. Proc R Soc Lond B. 2002;269:2381–2388. 60. Lehtonen J, Kokko H. Two roads to two sexes: unifying gamete competition and gamete limitation in a single model of anisogamy evolution. Behav Ecol Sociobiol. 2011;65:445–459. doi:10.1007/s00265-010-1116-8. 61. Kisdi ´E, Stefan A, Geritz H. Adaptive dynamics: a framework to model evolution in the ecological theatre. Journal of mathematical biology. 2010;61(1):165. 62. D´ ebarre F, Nuismer S, Doebeli M. Multidimensional (co) evolutionary stability. The American Naturalist. 2014;184(2):158–171. August 16, 2024 30/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint 63. Geritz SA, Metz JA, Rueffler C. Mutual invadability near evolutionarily singular strategies for multivariate traits, with special reference to the strongly convergence stable case. Journal of mathematical biology. 2016;72(4):1081–1099. 64. Constable GWA, Kokko H. Parthenogenesis and the Evolution of Anisogamy. Cells. 2021;10(9):2467. 65. M¨ uller J, Hense BA, Fuchs TM, Utz M, P¨ otzsche C. Bet-hedging in stochastically switching environments. Journal of theoretical biology. 2013;336:144–157. 66. Gillespie DT. Exact stochastic simulation of coupled chemical reactions. The journal of physical chemistry. 1977;81(25):2340–2361. 67. Rafaluk-Mohr C, Ashby B, Dahan DA, King KC. Mutual fitness benefits arise during coevolution in a nematode-defensive microbe model. Evolution Letters. 2018;2(3):246–256. 68. Agrawal AF. Evolution of sex: why do organisms shuffle their genotypes? Current Biology. 2006;16(17):R696–R704. 69. Lehtonen J, Jennions MD, Kokko H. The many costs of sex. Trends in ecology & evolution. 2012;27(3):172–178. 70. Nickerson AW, Raper KB. Macrocysts in the life cycle of the Dictyosteliaceae. II. Germination of the macrocysts. American Journal of Botany. 1973;60(3):247–254. 71. Nieuwenhuis BP, James TY. The frequency of sex in fungi. Philosophical Transactions of the Royal Society B: Biological Sciences. 2016;371(1706):20150540. 72. Waffenschmidt S, Woessner JP, Beer K, Goodenough UW. Isodityrosine cross-linking mediates insolubilization of cell walls in Chlamydomonas. The Plant Cell. 1993;5(7):809–820. 73. Gerber N, Kokko H. Abandoning the ship using sex, dispersal or dormancy: multiple escape routes from challenging conditions. Philosophical Transactions of the Royal Society B: Biological Sciences. 2018;373(1757):20170424. August 16, 2024 31/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint 74. Kokko H. When synchrony makes the best of both worlds even better: How well do we really understand facultative sex? The American Naturalist. 2020;195(2):380–392. 75. Otto SP, Lenormand T. Resolving the paradox of sex and recombination. Nature Reviews Genetics. 2002;3(4):252–261. 76. Hastings. The costs of sex due to deleterious intracellular parasites. Journal of Evolutionary Biology. 1999;12(1):177–183. 77. Togashi T, Cox PA, Bartelt JL. Underwater fertilization dynamics of marine green algae. Mathematical biosciences. 2007;209(1):205–221. 78. Nowak MA. Five rules for the evolution of cooperation. science. 2006;314(5805):1560–1563. 79. Bowler C, Allen AE, Badger JH, Grimwood J, Jabbari K, Kuo A, et al. The Phaeodactylum genome reveals the evolutionary history of diatom genomes. Nature. 2008;456(7219):239–244. 80. Sch¨ onknecht G, Chen WH, Ternes CM, Barbier GG, Shrestha RP, Stanke M, et al. Gene transfer from bacteria and archaea facilitated evolution of an extremophilic eukaryote. Science. 2013;339(6124):1207–1210. 81. Neiman M, Lively CM, Meirmans S. Why sex? A pluralist approach revisited. Trends in ecology & evolution. 2017;32(8):589–600. 82. Hurst LD, Nurse P. A note on the evolution of meiosis. Journal of Theoretical Biology. 1991;150(4):561–563. 83. Cavalier-Smith T. Cell cycles, diplokaryosis and the archezoan origin of sex. Archiv f¨ ur Protistenkunde. 1995;145(3-4):189–207. 84. Wilkins AS, Holliday R. The evolution of meiosis from mitosis. Genetics. 2009;181(1):3–12. 85. Crow J. An advantage of sexual reproduction in a rapidly changing environment. Journal of Heredity. 1992;83(3):169–173. August 16, 2024 32/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint 86. B¨ urger R. Evolution of genetic variability and the advantage of sex and recombination in changing environments. Genetics. 1999;153(2):1055–1069. 87. Kaiser D, Manoil C, Dworkin M. Myxobacteria: cell interactions, genetics, and development. Annual Reviews in Microbiology. 1979;33(1):595–639. 88. Bonner JT. The origins of multicellularity. Integrative Biology: Issues, News, and Reviews: Published in Association with The Society for Integrative and Comparative Biology. 1998;1(1):27–36. 89. Verstrepen K, Derdelinckx G, Verachtert H, Delvaux F. Yeast flocculation: what brewers should know. Applied microbiology and biotechnology. 2003;61(3):197–205. 90. Day TC, M´ arquez-Zacar´ ıas P, Bravo P, Pokhrel AR, MacGillivray KA, Ratcliff WC, et al. Varied solutions to multicellularity: The biophysical and evolutionary consequences of diverse intercellular bonds. Biophysics reviews. 2022;3(2):021305. 91. de Carpentier F, Lemaire SD, Danon A. When unity is strength: the strategies used by Chlamydomonas to survive environmental stresses. Cells. 2019;8(11):1307. 92. Strassmann JE, Queller DC. Evolution of cooperation and control of cheating in a social microbe. Proceedings of the National Academy of Sciences. 2011;108(supplement 2):10855–10862. 93. Douglas T, Queller DC, Strassmann JE. Social amoebae mating types do not invest unequally in sexual offspring. Journal of evolutionary biology. 2017;30(5):926–937. 94. Perrin N. What uses are mating types? The “developmental switch” model. Evolution: International Journal of Organic Evolution. 2012;66(4):947–956. 95. Raven JA. Allometry and stoichiometry of unicellular, colonial and multicellular phytoplankton. New phytologist. 2008;181:295–309. August 16, 2024 33/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint 96. Kooijman SALM. Dynamic energy budget theory for metabolic organisation. Cambridge university press; 2010. 97. Koyama T, Yamada M, Matsuhashi M. Formation of regular packets of Staphylococcus aureus cells. Journal of Bacteriology. 1977;129(3):1518–1523. 98. Iyer-Biswas S, Wright CS, Henry JT, Lo K, Burov S, Lin Y, et al. Scaling laws governing stochastic growth and division of single bacterial cells. Proceedings of the National Academy of Sciences. 2014;111(45):15912–15917. 99. Marantan A, Amir A. Stochastic modeling of cell growth with symmetric or asymmetric division. Physical Review E. 2016;94(1):012405. 100. Christie JR, Schaerf TM, Beekman M. Selection against heteroplasmy explains the evolution of uniparental inheritance of mitochondria. PLoS genetics. 2015;11(4):e1005112. 101. Maire N, Ackermann1 M, Doebeli M. Evolutionary branching and the evolution of anisogamy. Selection. 2002;2(1-2):119–131. 102. Parker GA. The sexual cascade and the rise of pre-ejaculatory (Darwinian) sexual selection, sex roles, and sexual conflict. Cold Spring Harbor perspectives in biology. 2014;6(10):a017509. 103. Tarnita CE, Washburne A, Martinez-Garcia R, Sgro AE, Levin SA. Fitness tradeoffs between spores and nonaggregating cells can explain the coexistence of diverse genotypes in cellular slime molds. Proceedings of the National Academy of Sciences. 2015;112(9):2776–2781. 104. Pichugin Y, Pe˜ na J, Rainey PB, Traulsen A. Fragmentation modes and the evolution of life cycles. PLoS computational biology. 2017;13(11):e1005860. 105. Siljestam M, Martinossi-Allibert I. Anisogamy does not always promote the evolution of mating competition traits in males. The American Naturalist. 2024;203(2):230–253. 106. Krumbeck Y, Constable GWA, Rogers T. Fitness differences suppress the number of mating types in evolving isogamous species. Royal Society Open Science. 2020;7:192126. August 16, 2024 34/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint 107. Czuppon P, Constable GWA. Invasion and extinction dynamics of mating types under facultative sexual reproduction. Genetics. 2019;213(2):567–580. August 16, 2024 35/35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted August 20, 2024. ; https://doi.org/10.1101/2024.08.19.608569doi: bioRxiv preprint