Cell size and selection for stress-induced binary cell fusion

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


Introduction
Although the details of the early evolution of sexual reproduction in the last common eukaryotic common ancestor (LECA) are shrouded in mystery, it is argued that the emergence of eukaryotic sex began with the evolution of cell-cell fusion and meiosis [1] in an archaeal ancestor [2,3].This step can be further broken down into the evolution of binary cell fusion, the one spindle apparatus, homologous pairing and chiasma, and finally reduction, division and syngamy [4].The vast majority of theoretical studies investigating the evolution of sexual reproduction have focused on later stages of this evolutionary trajectory, namely the conditions that give rise to a selective pressure for genetic recombination [5][6][7][8].However, comparatively few studies have investigated the selective pressures that may have first given rise to binary cell fusion, which may have facilitated the evolution of a host of other eukaryotic traits [9], including the homologous pairing and meiotic recombination, by bringing nuclei together in the same cell.
Hypotheses for the evolution of binary cell fusion often rely on hybrid fitness advantage.It has been suggested that selection for cell-cell fusions might have initially been driven by "selfish" transposons and plasmids [10][11][12], or negative epistatic August 16, 2024 2/35 interactions between mitochondrial mutations [13,14].However, once a heterokaryotic cell has been formed (binucleate with nuclei from both parental cells), the advantage of hybrid vigor and the masking of deleterious mutations could lead to the maintenance of cell fusion [4].Such benefits are required to alleviate costs to cell-fusion, which include selfish extra-genomic elements in the cytoplasm [15] and cytoplasmic conflict [16,17].
In these previous studies on the evolution of binary cell fusion, the effect of changing environmental conditions is not considered.However, in many extant unicellular organisms, binary cell fusion (and the karyogamy and genetic recombination that follow) occur in response to challenging environmental conditions [18] such as starvation (Tetrahymena [19,20]) and depleted nitrogen levels (Chlamydomonas reinhardtii [21] and Saccharomyces pombe [22]).Meanwhile in benign conditions with abundant resources these species reproduce asexually using binary cell fission.The mechanisms that drive selection for genetic recombination under challenging environmental conditions are well-studied [23]; recombination can facilitate adaptation to a novel environment [24,25] and evolving to engage in more sex when fitness is low (fitness-associated sex) can allow an organism to maximise the advantages of sex while minimising the costs [26,27].However, this focus on the benefits of recombination leaves space to ask whether binary cell fusion itself could be selected for as a stress response, even in the absence of any genetic advantages.
In this paper, we do not account for the genetic factors discussed above.Instead, we focus only on how the survival advantage associated with increasing cytoplasmic volume might select for binary cell fusion; this relies on the physiological advantages conferred by cell-cell fusion and is independent of the question of the genetic advantages (and disadvantages) of sexual reproduction.This alternative perspective offers useful new insights that can be compared with empirical observation.
That size-based processes could play a role in the early evolution of sexual reproduction has empirical and theoretical support.The "food hypothesis" [28] suggests that metabolic uptake could drive horizontal gene transfer in bacteria and archaea, with DNA molecules providing nutrients for the receiving cell [29,30].Indeed, horizontal gene transfer has been shown experimentally to be an important source of carbon and nutrients in bacteria [31,32].Binary cell fusion is possible in bacteria (where it has been shown to come with selective benefits from mixed cytoplasm [33]) and archaea [34,35].
August 16, 2024 3/35 Meanwhile amongst eukaryotes, the benefits of increasing cytoplasmic volume are understood to be strong enough to drive selection for the sexes themselves [36,37].That early selection for syngamy may have been driven by the survival benefits of larger cytoplasmic volume (and the concomitant increase in solid food reserves) has been argued verbally [38], however has not been explored mathematically.In suggesting a mechanistic hypothesis for the evolution of binary cell fusion, our work has interesting parallels with [39], where an advantage to cell fusion is identified in terms of shortening the cell-cycle.
Moving to consider potential physiological benefits of binary cell fusion naturally leads to work on the evolution of multicellularity.While multicellularity and binary cell fusion are clearly biologically distinct, from a modelling perspective they share similarities in that they can involve the "coming together" of cells to produce a larger complex [40].Multicellularity achieved via aggregation allows organisms to rapidly adapt to novel environments that favour increased size [41,42].It has also been suggested that the genetic nonuniformity of such aggregates may also make them well-suited to resource limited environments [43,44], echoing the hybrid vigor hypotheses for the evolution of early syngamy [4].Relatively few theoretical studies have investigated the evolution of facultative aggregation in response to changing environments [42].However in the context of clonal multicellularity ("staying together") such environments have been considered more extensively [45,46].In this clonal context, the evolutionary dynamics act primarily on fragmentation modes [47,48] (e.g.how a "parental" multicellular complex divides to form new progeny).Interestingly the same quality-quantity trade-off arises here [45] as drives selection for the sexes [36] (anisogamy, gametes of differing sizes); larger daughter cells (or gametes) are more able to withstand unfavourable environmental conditions, while smaller cells can be produced in larger quantities.
In this paper we adapt the classic Parker-Baker-Smith [36] (PBS [37]) model for the evolution of sexes in order to investigate the evolution of binary cell fusion.This builds on recent work that investigates how the possibility of parthenogenetic reproduction can drive selection for oogamy in eukaryotes [49].We assume for simplicity that parental cells undergo a number of cell-divisions.The size of daughter cells is a compound evolvable trait determined both by the size of parental cells and the number of cell

Model Insights from Simulations
We consider a computational model of a haploid population that reproduces via binary cell fission.This population dynamics proceed as follows: there is an initial growth and binary fission phase where a fixed energy budget E can be used for population growth and binary fission (see Fig 1A).The parameter E thus sets the carrying capacity of the population.This growth phase is followed by an environmentally-induced mortality phase where growth and fission are suspended and survival depends on cell size.
Explicitly, we assume that cells reach maturity at size M , so that in the absence of 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) .
However by decreasing M and increasing n, individuals also produce smaller daughter cells that are more vulnerable to extrinsic mortality.The size of daughter cells is thus subject to a quality-quantity trade-off.For simplicity, we model the mass of daughter cells m as a continuous trait, and explore its evolution using simulations (see Supporting Information 3.3).daughter cells evolve to become larger to withstand more adverse conditions).
We now modify the model to allow for the possibility of binary cell fusion following the cell fission described above.Daughter cells may now fuse to form a binucleated cell (e.g. a dikaryon [50], in which the cytoplasm of the contributing cells are mixed but August 16, 2024 6/35 their nuclei or nucleoids remain distinct [38]) or remain a mononucleated cell (with a single nucleus or nucleoid).The rate of cell fusion is given by α, such that when α = 0 all cells remain mononucleated, and cell survival into the next growth cycle is calculated as before.Conversely for α > 0, some proportion of daughter cells will have fused and in the limit α → ∞, all cells will have fused.These fused cells will receive a survival advantage from their increased mass.However they will also pay an additional cost, C, resulting from factors such as cell-fusion failure [51], selfish extra-genomic elements in the cytoplasm [15], cytoplasmic conflict [16,17] and maintenance of a binucleated cell [52].Together this means that fused cells survive with a total probability (1 − C)S(2m; β).Surviving adults divide to form a new growth cycle of mononucleated haploid daughter cells, with binucleated parental cells producing mononucleated progeny through vegetative segregation [53] (or alternatively through plasmid segregation machinery [54]).Note that although we do not account for the possibility of binucleated cells failing to form mononucleated progeny (i.e.failed segregation), this can be accounted for by their additional survival cost, C (see Fig 1B).
We now explore the coevolution of daughter cell mass, m, and fusion rate, α.In The result above is in some sense surprising.Despite the presence of additional survival costs associated with binary cell-fusion, selection for non-zero fusion rates (rather than increased daughter cell size) persists in the harsh environment.We explain the emergence of this behaviour mathematically in the Results section.

Mathematical Model
Our model takes inspiration from the classic PBS model for the evolution of anisogamy [36] (the production of sex cells of differing size).However, whereas such models typically consider the binary cell fusion (fertilization) rate a fixed parameter, we here treat it as a trait subject to evolution.In doing so our work builds on [49], where a August 16, 2024 7/35 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 evolution of anisogamy with parthenogenesis.In order to analyse the dynamics of the model, we use tools from adaptive dynamics [55], assuming that traits are continuous and that mutations have small effect.
In addition we will explore the effect of switching environments, another departure from the PBS model.

Dynamics within each growth cycle
A total of (2 n E)/M daughter cells enter a pool in which binary cell fusion can occur.
After a finite time window, the resultant cells are subject to a round of mass dependent mortality, such that cells of larger mass are more likely to survive.The surviving cells form the basis of the next growth cycle, as illustrated in Fig 1.

Fusion Kinetics
We assume that all daughter cells may fuse with each other, an assumption consistent with most models of the early evolution of sexual reproduction, which suppose the existence of a "unisexual" early ancestor that mated indiscriminately [56].Following growth and binary cell fission, the population is comprised of N unfused daughter cells.
Fusion between these mononucleated cells occurs at a rate of α, such that the number of unfused cells, N , is given by the solution to At the end of the fusion window of duration T there are then N (T ) unfused (mononucleated) cells remaining, and (N (0) − N (T ))/2 fused (binucleated) cells.

Survival Probability
We assume that both unfused and fused cells are subject to the same extrinsic mass-dependent mortality function, S(m; β), while fused cells pay an additional mass-independent cost C.Many choices for such a function are possible, so long as it is an increasing function of cell size (which we equivalently refer to as cell mass m).
August 16, 2024 9/35 However here we assume that S(m; β) is the Vance survival function [57], a common assumption in the literature [58][59][60].We thus have that at the end of the fusion window, the survival probability of unfused and fused cells are given respectively by where m i is the mass of a particular unfused daughter cell and m j and m k are the masses of two daughter cells that have fused.For a given cell mass, increasing β will decrease the survival probability.We therefore refer to β as the environmental harshness parameter, with high β corresponding to harsh environments in which survival is difficult, and low β corresponding to more benign environments in which even cells of modest mass have a high probability of surviving.
Having defined how the survival of a cell depends on its mass, we have the necessary tools to mathematically characterise the fitness of a rare mutant, and whether it can invade the resident population.In the following section, we provide mathematical approximations of the invasion dynamics of such a mutant.

Invasion Dynamics
Adopting the classical assumptions of adaptive dynamics [55,61] (see also Supporting Information 3.1 and [49]), we mathematically approximate the invasion dynamics of a mutant (which occur over discrete growth cycles).Deriving these invasion dynamics analytically is only possible when we assume that mutations in m and α occur independently.However the evolutionary dynamics we obtain if we consider mutations occurring in both m and α simultaneously remains identical to those obtained by assuming that they occur independently (see Supporting Information 4 and [49]).
Denoting by fm the frequency of mutants of size m ± δm in the population where δm is the mutational stepsize in m, which is assumed to be small and t g the number of growth cycles, we find (see Supporting Information 2) August 16, 2024 10/35 where h m (m, α, β, C) is a constant that depends on the parameters m, α, β, C (see Eq. (S10)).This constant provides the fitness gradient of the mutant.Similarly, denoting by fα the frequency of mutants with fusion rate α ± δα in the population where δα is the mutational stepsize in α, which is again assumed to be small, we find where h α (m, α, β, C) is the fitness gradient of a mutant with fusion rate α + δα.We see that in the case of a single mutant, we have frequency-independent selection for mutants with different masses and fusion rates.We note that in reality, frequency-dependent invasion dynamics can occur when multiple mutants that change both m and α arise in the population (see Supporting Information 4 for mathematical analysis), which can in turn lead to evolutionary branching [62,63].However, since this branching does not occur in the regimes we are focusing on in this paper [49], we assume for simplicity that mutants encounter a monomorphic resident population (trait substitution) for the remainder of the mathematical analysis.

Evolutionary Dynamics
We assume that haploid daughter cells are characterised by two genetically non-recombining traits mass m and cell fusion rate α.We assume that mutations occur in m or α independently at a fixed rate µ, where µ is measured in units of (number of growth cycles) −1 (see Supporting Information 3.3 and [49]).A mutation in m represents a change in the mass of the daughter cell produced, and a mutation in α represents a change in the fraction of the population that undertakes either one of the reproductive routes (i.e.binary cell fusion vs strictly binary cell fission).
Mutants with a different mass to their ancestor can produce either more or fewer daughter cells than their ancestor (see Eq. ( 1)), which impacts their survival (see Eq. ( 2)).When mutants have a different fusion rate to their ancestor, although the number of daughter cells produced does not differ from their ancestor, the number of fused cells at the end of a growth cycle can either increase/decrease, which impacts their survival, since fused cells have greater mass.The survival of fused cells is also influenced by the cost of fusion C (see Eq. ( 2)).
August 16, 2024 11/35 Our mathematical analysis in the remainder of this section assumes that mutants encounter a strictly monomorphic population (i.e. that mutations fixate before the introduction of a new mutant).However in our numerical simulations, we release this restriction and stochastically allow for the coexistence of multiple traits in the population, held under a mutation-selection balance, as described in the subsequent section.

Fixed Environment
We first consider the evolutionary dynamics in the case where the environment is fixed (i.e. when the parameter β, which measures the harshness of the environment (see Eq. ( 2)), is constant throughout the evolution.Assuming that δm and δα are small (small mutational step size), we use techniques from adaptive dynamics [49,64] to obtain equations for the evolutionary dynamics of m and α, which are given by for α ≥ 0. This boundary is imposed to prevent α from becoming negative, which is biologically unrealistic since it corresponds to an increase in daughter cell numbers during the fusion period, as can be seen from Eq. (1).Therefore when α becomes strictly decreasing along this boundary α = 0 boundary (i.e [dα/dτ ]| α=0 < 0 in Eq. ( 5)), we introduce a discontinuous change in the dynamics, given by The derivations of these equations can be found in Supporting Information 3.1.
August 16, 2024 12/35 Switching environments with phenotypic plasticity Now, we consider the case where evolution acts on the same traits as before, but the environment is subject to change.We model environmental change as switching between two environments β 1 and β 2 .If β 1 > β 2 , then β 1 is the harsher environment (see Eq. ( 2)).We also allow for phenotypic plasticity such that the population can evolve different strategies in different environments.The population's evolutionary state is now described by four traits; the daughter cell mass in environments 1 and 2 (m 1 and m 2 ) and fusion rate in these environments (α 1 and α 2 ).
For simplicity we assume that any cost of phenotypic switching or environmental sensing is negligible and that this plastic switching is instantaneous upon detection of the change in environmental conditions.The evolutionary dynamics in each environment are then decoupled.However the evolutionary trajectories in each environment are coupled by the initial trait values for the population in each environment, which we assume are the same (i.e. the population begins in a phenotypically undifferentiated state).With phenotypic plasticity, the evolutionary dynamics are then given by C) with initial conditions Here, H m (m, α; β, C) and H α (m, α; β, C) retain the functional form in Eq. ( 5).
As Eqs. (7) are only coupled through their shared initial conditions, m 0 and α 0 , the choice of these initial conditions is an important consideration.Since we are interested in the initial evolution of binary cell fusion, it is natural to assume that the population evolves from a state of zero fusion, α 0 = 0. Deciding on a plausible initial daughter cell mass takes more thought.One parsimonious choice would be that the population is already adapted to either environment 1 or environment 2 and that m 0 is given by an evolutionary fixed point in one of these environments (this is the situation illustrated in phenotypic placticity has evolved, it is possible that m 0 is instead given by a bet-hedging strategy.We explore that such a strategy would take in the following section.

Switching environments without phenotypic plasticity
We now consider the case where there is switching between environments but where the population exhibits no phenotypic plasticity.As described in the previous section, we are particularly concerned with the period before the physiological machinery for cell fusion has evolved, and so focus on the case where the cell fusion rate is fixed to zero, α = 0. Evolution then solely acts on the daughter cell mass, m.
As in [49], environmental switching is modelled as a discrete stochastic telegraph process, with the time spent in each environment distributed geometrically.The population spends an average of τ 1 = 1/λ 1→2 in environment 1 and τ 1 = 1/λ 2→1 in environment 2, where λ i→j is the transition rate from environment i to j.
The two switching rates most relevant to our model are when the environment switches many times before an invasion can complete, (fast relative to invasion, FRTI) and when each switching event occurs after multiple invasions have completed, (fast relative to evolution, FRTE).More detail of these switching rates are provided in Supporting Information 5. However in [49], we show that the evolutionary dynamics for m in both these regimes can be approximated using the same dynamical equations.
Using adaptive dynamics techniques modified to account for such environmental switching [65], we obtain where H m (m, α; β, C) retains the functional form in Eq. ( 5) and P 1 = τ 1 /(τ 1 + τ 2 ) and are the probabilities of finding the population in the two respective environments.We therefore see that in the absence of phenotypic plasticity, the evolutionary dynamics is the weighted average of the dynamics in the two environments.
Obtaining the ESS for Eq. ( 9) is relatively straightfoward.Substituting for H m (m, 0; β 1 , C) and H m (m, 0; β 2 , C) using the functional form given in Eq. ( 5) and August 16, 2024 14/35 setting dm/dτ = 0 in Eq. ( 9), we obtain the ESS This strategy constitutes a bet-hedging strategy in cell mass when the population has yet to evolve phenotypic plasticity nor the capacity for cell-cell fusion.In the limits P 1 → 1 and P 2 → 1, we can recover the ESS strategies in the two respective environments: which can be verified from a consideration of the equations for dm/dτ = 0 in a fixed environment with α = 0 (see Eq. ( 6)).We can now proceed to analyse how binary cell fusion can be selected for when the fusion rate α is allowed to increase from zero in the following Results section.

Implementation of Numerical Simulations
The stochastic simulations of the evolutionary trajectories are also implemented using a Gillespie algorithm [66] where successive mutations and environmental switching events occur randomly with geometrically distributed waiting times.The rates of mutations µ and environmental switching λ are measured in units of (number of growth cycles) −1 .
In the simulations, multiple traits coexist under a mutation-selection balance (see Supporting Information and [49] and [67] for more detail), which allows us to account for variations in selection strengths in simulations of our evolutionary trajectories.

Results
In this section we proceed to analyse the evolutionary dynamics derived from the mathematical model and compare our results to numerical simulations of the full stochastic simulations.
August 16, 2024 15/35 In a fixed environment the population evolves to either no cell fusion, or to high levels of cell-fusion, dependent on the cost of cell fusion In Fig 3, we see two possible evolutionary outcomes for the fusion rates in a fixed environment; the population can evolve either a high (technically infinite) fusion rate or to zero fusion rates.To which of these fusion rates the population is attracted depends both on the parameters and the initial conditions.
When the costs to cell fusion are low (C ⪅ 0.39), the only evolutionary fixed point is the high fusion rate fixed point (see Fig 3A).In this scenario, obligate fusion is the only evolutionary outcome.A rigorous mathematical analysis that formalises the arguments above are provided in [49], which uses a similar model to investigate the evolution of anisogamy with parthenogenesis.In summary, the set of possible evolutionary attractors (m * , α * ), starting from an initial condition (m(0), α(0)) = (m (0) , 0), are given by where we note 1 − e −1/2 ≈ 0. 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 × 10 7 growth cycles in panel A, 1.24 × 10 7 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 relevant for the evolution of early cell fusion; if evolution had acted on daughter cell size, m, before the physiological machinery necessary for cell fusion had evolved, the initial condition for the co-evolutionary dynamics would be (m(0), α(0)) = (β, 0), at which the population would be subsequently held by costs to fusion.
In Fig 3 we also see that our mathematical analysis is a good predictor of the outcome of stochastic simulations (gray shaded lines).One minor point of departure is that at high fusion rates our simulated trajectories begin to diverge from our analytic prediction.This discrepancy is the result of evolutionary branching in cell mass, which we explore in another paper relating to the emergence of size dimorphism in sex cells [49].However this branching happens at a later evolutionary stage than the focus of this study, the early emergence of binary cell fusion.
We conclude this section by addressing the key biological result that arises from this August 16, 2024 17/35 analysis; cell fusion is uniformly selected for even under moderately high costs (with a fraction of up to C ≈ 0.39 of fused cells failing to survive) and can even be selected for under higher costs (up to C ≈ 0.86) given necessary initial conditions.In the context of the evolution of early binary cell fusion, this provides a surprising nascent advantage to cell fusion.This advantage could even help compensate for other short-term costs arising from the later evolution of sex and recombination.The selective advantage experienced by fusing cells comes from their increased cytoplasmic volume, which leads to increased survival probabilities.
In a switching environment with phenotypic plasticity, binary cell fusion can evolve as a facultative stress response to harsh environments Having considered the case of the evolutionary dynamics in a fixed environment, we now move on to consider the evolutionary dynamics of a population exhibiting phenotypic plasticity in a switching environment (see Eq. ( 7)).We recall that under the assumptions of costless and immediate phenotypic switching, the dynamics of (m 1 , α 1 ) and (m 2 , α 2 ) are decoupled.The evolution of the traits in the respective environments are coupled however through the initial conditions from which they evolve, which must be the same (i.e. a phenotypically undifferentiated state).
We consider two parsimonious choices for these initial conditions, both beginning in a state without fusion (α 1 (0) = α 2 (0) = 0).In the first scenario, we assume that the population has evolved to a stable non-fusing mass adapted to a single environment (see Eq. ( 11) . This is a situation in which the alternate environment is in some sense novel and one to which the population has not adapted.In the second scenario, we instead assume that the population has evolved to a bet-hedging strategy in mass (optimising the mass of daughter cells across the two environments) such that m 1 (0) = m 2 (0) = m * BH,α=0 (see Eq. (10)).An illustrative phase portrait is shown in Fig 4.
In both scenarios described above, we see the emergence of facultative binary cell-fusion as a response to harsh environmental conditions that lower the survival probability of daughter cells.However we note that this is only possible if there is an appreciable increase in environmental harshness, β, between the environments.In both selected for and against depending on the value of m (see Supporting Information 6.2 and Eq. ( 12)); this restricts us to the more interesting parameter regime in which different evolutionary outcomes are possible in each environment.
In Fig 5, we see that when one of the environments is not appreciably worse than the other, binary cell fusion does not evolve in either environment.However when the difference between the environments grows more substantial, it is possible to evolve cell fusion in the harsher environment from initial condition (m 1 (0), α 1 (0)) = (m 2 (0), α 2 (0)) = (β i , 0) (where the population has first evolved towards the evolutionary optimum of the more benign environment).Finally when the difference between the environments is extreme, it is also possible to evolve cell fusion in the harsher environment from initial condition (m (where the population has first evolved towards a bet-hedging strategy in cell mass).Here, C = 0.5, P 1 = 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.

Discussion
The evolution of sexual reproduction and its consequences for the subsequent evolutionary trajectory of populations is of general importance in biology [7,13,68].In this paper we have illustrated a reversal of the classic two fold cost of sex that appears in organisms with distinct sexes [69]; in unisexual [56], unicellular organisms, binary cell August 16, 2024 20/35 fusion can be selected for, even in the presence of substantial costs, due to a survival benefit that comes from increased mass.These results allow us to quantitatively assess the verbal hypothesis that syngamy evolved allow cells to store more food reserves and thus increase their survival rate [38].It is particularly interesting that the benefits conferred to cell fusion through increased cytoplasmic mass are sufficient to withstand remarkably high costs; "obligate sexuality" is the only evolutionary outcome with costs equivalent to a loss of ∼ 39% cells that attempt to fuse, and remains a potential outcome with costs of up to ∼ 86% of fused cells dying.
Perhaps most interesting is the case of switching environments with phenotypic plasticity.Here we find under a broad set of biologically reasonable conditions (costs to cell fusion equivalent to 39% − 86% additional mortality to fused cells and at least moderate changes in environmental quality) that high fusion rates are selected for in harsh environments and zero fusion rates are maintained in benign environments.This behaviour parsimoniously recapitulates the empirically observed reproductive strategies of numerous facultatively sexual species, including C. reinhardtii [21], S. pombe [22] and D. discoideum [70].This mechanism, under which cell fusion evolves to increase the survival probability of daughter cells, provides a complementary perspective on the frequent evolution of survival structures (resistant to environmental stress) that form following the formation of a zygote.These include ascospores in fungi [71] and zygote-specific stress-resistant stress wall in C. reinhardtii [72].Note that such correlations between sexual reproduction and the formation of survival structures are not as easily explained under genetic explanations for the evolution of sexual reproduction, where engaging in both behaviours at once constitutes a simultaneous (and therefore potentially costly) change in genotype and temporal dislocation in environment [73,74].
The results above are particularly interesting in the case of the evolution of early binary cell-fusion as a first step in the evolution of sexual reproduction.While most studies focus on the genetic benefits of cell-fusion [75] (including a functionally-diploid dikaryotic cell [4]), or the genetic benefits of mixed cytoplasm [13,14] (which can also come with costs [15-17, 40, 76-78]), the mechanism at play here is purely physiological.
Yet, as addressed above, it naturally captures the empirical observation of binary cell-fusion as response to challenging environmental conditions, a feature absent in these August 16, 2024 21/35 earlier models.While the mechanism does not explain the evolution of sexual reproduction and genetic recombination itself, it does provide a nascent advantage to binary cell fusion that sets the stage for the evolution of sex by bringing nuclei from different cells into contact for prolonged periods.The mechanism also shows the potential to counter short-term costs associated with the initial formation of a binucleated cell.In this way the mechanism could facilitate the transition from horizontal gene transfer [79,80] to meiotic recombination, which is advantageous when genome sizes increase in length [8].Conceivably, if genetic recombination is beneficial for myriad genetic reasons in the long-term [8,81], it would seem natural that it would be instigated when the opportunity arises (i.e. when physiological survival mechanisms bring nuclei into close contact).We note that it is obviously possible that the first diploid cells arose by errors in endomitosis [82][83][84] (essentially doubling the chromosome number within a single cell) and that meiosis first evolved in this context.Such a sequence of events is still compatible with our very general model, which can alleviate short-term costs of sex such as the energy involved in finding a partner and undergoing fusion [69].In either scenario sexual reproduction may not be only a direct response to environmental variability [85,86], but also to the correlated formation of a survival structure.
More generally, it is interesting to note that the conditions for facultative sexuality (e.g.harsh environmental conditions) broadly coincide with those for facultative multicellularity in both bacteria and eukaryotes, with starvation triggering the formation of fruiting bodies in myxobacteria [87,88] and flocking in yeast [89,90].
Meanwhile in C. reinhardtii, the formation of multicellular palmelloids and aggregates are an alternate stress response to sexual reproduction [91], as are the formation of fruiting bodies in D. discoideum [92].In this multicellular context, the sexual behaviour of D. discoideum is particularly interesting, as once formed, the zygote attracts hundreds of neighboring cells that are then cannibalised for the provision of a macrocyst [93].These various survival strategies are unified in our model as a mechanism for the evolution of binary cell fusion.
One element absent from our model is the fusion of multiple cells, which is likely to be selected for under the assumptions implicit in our model.There would clearly be an upper-limit on the number of fusions selected for, arising from the likely multiplicative August 16, 2024 22/35 effect of the fusion cost C. However in this context, it is interesting to note that one of the hypotheses for the evolution of self-incompatible mating types is as a signal to prevent the formation of polyploid cells [94].Such a mechanism could also prevent the formation of trikaryotic cells should the cost of multiple fusions be too great.Thus, the model neatly preempts the second stage in models for the evolution of eukaryotic sex, the regulation of cell-cell fusion [1].
The selection pressure for binary cell fusion is dependent on the formation of a binucleated cell with an increased survival benefit arising from its larger size.Although it is reasonable to assume that a single mononuclear cell has lower total resource requirements than a multicellular complex of the same size [95], we have not considered the detailed energetics of the maintenance of two nuclei [9].In these respects incorporating dynamic energy budget theory into the model would be an important next step [96] as it would provide a clear distinction between the survival benefits of fused cells and unfused bicellular complexes.Within our modelling framework, these two structures are broadly similar [40].However as we have shown, increased cell-cell attraction can be selected for even in the presence of large costs that one might expect under binary cell fusion but not associate with the formation of a bicellular complex.
We have assumed for simplicity a simple cell division scheme; parental cells undergo n rounds of symmetric division to produce 2 n daughter cells.In the context of multicellularity, switching environments have been shown to promote binary fragmentation [45].However non-symmetric modes can be selected for [47] reflecting the diverse modes of facultative multicellular life cycles observed in bacteria [97].It would be interesting to incorporate our results into models of cell division that account explicitly for growth [98,99] to determine how these results for multicellular organisms carry over to the unicellular scenario, and further how they may affect those we have shown here.
Finally, we have not explicitly modelled any sources of cytoplasmic or genetic conflict [100], which we have for simplicity included in the fusion cost C.Nevertheless, social conflict does emerge in this model.In a recent paper we have shown how evolutionary branching can arise, with some individuals producing fewer larger cells and others producing more numerous but smaller daughter cells [49].This branching is driven by the same evolutionary forces that drive selection for anisogamy [101], in which August 16, 2024 23/35 context this can be viewed as sexual conflict [102].That social conflict should arise in the formation of multicellular aggregates is well understood [40,88,103].However these models typically assume cells of fixed size [46,104].Combining the insights derived from the evolution of anisogamy literature with the theory developed in the multicellularity literature represents another promising research direction.
As addressed above, trade-offs between cell fusion rate and mass [105], cell-energetics, inbreeding, and the possibility of multiple cell-fusion events offer interesting avenues to extend this analysis.In addition, we have not accounted for the discrete nature of divisions leading to daughter cells, costs to phenotypic switching, non-local trait mutations, or pre-existing mating types.More generally, extending our mathematical approach leveraging adaptive dynamics to switching environments in other facultatively sexual populations might prove particularly fruitful [106,107].
In this paper we have adapted the classic PBS model [36] in two key ways; allowing the fusion rate to evolve and subjecting the population to switching environments.In doing so, we have shown its capacity to parsimoniously capture the evolution of obligate binary cell fusion, obligate binary cell fission, and stressed induced binary cell fusion in unicellular organisms.These results offer particularly interesting implications for the evolution of binary cell-fusion as a precursor to sexual reproduction, as well as suggesting common mechanistic links between the evolution of binary cell fusion and multicellularity.Moreover, our analysis emphasises the importance of exploring the coevolutionary dynamics of a range of evolutionary parameters, and of developing computational and mathematical approaches to elucidate facultative sexual reproduction.
Supporting information S1 Text.Supporting information.pdf

Fig 1 .
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.

Fig 2
Fig 2 summarises the outcome of such evolutionary dynamics.Fig 2A shows that the population evolves towards an evolutionarily stable strategy (ESS) in m for a given environment.Should the environment suddenly become harsher (via an increase in β) the population evolves towards a new ESS, in which daughter cells are larger (i.e.

Fig
Fig 2B, we see that in the benign environment, α remains at zero, and the population evolves towards an ESS in m as in Fig 2A.However now when the population is introduced to a harsher environment, the evolutionary dynamics differ from those in Fig 2B (where α was held artificially at zero).Rather than cells evolving to be larger, we see a different response emerging; selection for binary cell fusion (α > 0).
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).

Fig 2 )
Fig 2).However if the population has been exposed to both the environments before

For
intermediate costs to cell fusion 0.39 ⪅ C ⪅ 0.86, there are two possible evolutionary outcomes.The outcome depends on the initial conditions (see Fig 3B).If the initial mass on the α = 0 boundary is small, selection acts to increase fusion rate and obligate fusion is the ESS.However, if the initial mass on the boundary is sufficiently large, the state of no cell fusion becomes the evolutionarily stable state.Finally when costs to fusion are extremely high (C ⪆ 0.86, see Fig 3C), selection for decreased fusion rate acts regardless of the initial value of m on the α = 0 boundary, and a state in which α = 0 (zero fusion rate) is the only evolutionary outcome.Under this high cost regime, dα/dt < 0 along the entire line α = 0 and so fusion rate is never selected to increase given any initial daughter cell mass.

39 and 1
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 × 10 7 growth cycles in panel A, 1.24 × 10 7 growth cycles in panel B and 10 7 growth cycles in panel C. Remaining system parameters are given in Supporting Information 7.

Fig 5 ,
we summarise the key results over the β 1 − β 2 parameter plane.Here we assume that the cost to cell fusion is intermediate (1 − e −2 > C > 1 − e −1/2 , i.e. 0.86 ⪆ C ⪆ 0.39) such that there are regions on the boundary α = 0 at which increased fusion rates are August 16, 2024 19/35

Fig 5 .
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, P 1 = 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.
Daughter cells are then introduced to a pool in which they can undergo binary cell fusion, with larger fused cells experiencing a survival advantage.Unlike in the classic PBS model, the fusion rate is also a trait subject to evolution.In the following sections we proceed to outline some insights developed from numerical