Evolution of Complex Regulation for Cell-Cycle Control

Abstract Many questions remain about the interplay between adaptive and neutral processes leading to genome expansion and the evolution of cellular complexity. Genome size appears to be tightly linked to the size of the regulatory repertoire of cells (van Nimwegen E. 2003. Scaling laws in the functional content of genomes. Trends Gen. 19(9):479–484). In the context of gene regulation, we here study the interplay between adaptive and nonadaptive forces on genome and regulatory network in a computational model of cell-cycle adaptation to different environments. Starting from the well-known Caulobacter crescentus network, we report on ten replicate in silico evolution experiments where cells evolve cell-cycle control by adapting to increasingly harsh spatial habitats. We find adaptive expansion of the regulatory repertoire of cells. Having a large genome is inherently costly, but also allows for improved cell-cycle behavior. Replicates traverse different evolutionary trajectories leading to distinct eco-evolutionary strategies. In four replicates, cells evolve a generalist strategy to cope with a variety of nutrient levels; in two replicates, different specialist cells evolve for specific nutrient levels; in the remaining four replicates, an intermediate strategy evolves. These diverse evolutionary outcomes reveal the role of contingency in a system under strong selective forces. This study shows that functionality of cells depends on the combination of regulatory network topology and genome organization. For example, the positions of dosage-sensitive genes are exploited to signal to the regulatory network when replication is completed, forming a de novo evolved cell cycle checkpoint. Our results underline the importance of the integration of multiple organizational levels to understand complex gene regulation and the evolution thereof.

revealed the same network topology. In some cases, and in fact in most genomes of each population, more genes appeared, although without a substantial functional role.
One crucial improvement is the appearance of baseline activity of g1 (θ < 0). This abolishes the second non-functional attractor that was present in the network of Quiñones-Valles et al. (2014) and which has become a liability in our model due to the stochasticity in regulatory dynamics. Another notable change is the appearance of the activation of g3 by g5. Although this change likely did not alter the cell-cycle in this simple initial experiment, it became an important subject for evolution in the ten replicate evolution experiments (see main text, Fig. 5).

Gene family analysis
To assess genotypic variation in populations, and study specific patterns in that variation (i.e. Fig. S1), we performed a gene family analysis (Fig. S2, S3). We defined gene orthologies per population and assigned to these the genes from all genomes in that population. For each population, we find a small set of core gene families (mean copy number per cell > 0.95 and present in > 95% of unique genomes), which always include g1-5 (Fig. S2a). Out of the non-core families, the large majority is only found in a few 1 Figure S1: The architecture of gene regulatory networks varies between replicate evolution experiments but is comparable within each experiment. While some variation across the gradient is functionally relevant (e.g. the extra gene g8 in sector 9 of R8), most variation is likely neutral (e.g. the two extra genes g7 and g11 in sector 7 of R8, neither of which was observed to be expressed during phenotypic profiling for over 1500 timesteps).
genomes (mean copy number per cell < 0.1, present in < 10% of unique genomes), which accounts for less than 2% of gene content on average. Two interesting exceptions are g8 in R8 and g9 in R2, which are ubiquitous in one half of the gradient and completely absent in the other (mean copy number ≈ 0.5, see Fig. S2a). Core gene families do not show such sharp presence-absence variation across the gradient, but some of them vary in copy numbers. By comparing variation between sectors to variation within sectors (using ANOVA's η 2 = SSgroups SSgroups+SSerrors ), we can quantify several types of specialization to the gradient. Besides copy number specialization, we obtain η 2 for 9 other properties of the 69 core gene families: activation threshold θ, regulatory weight w, position on the genome χ, number of binding sites n b , regulatory input (k inS and k inW , for edges with strong and weak affinity, respectively), regulatory output (k outS and k outW ) and binding sequence B (see Fig. S3). Figure S2c,d shows the degree of specialization (η 2 (c) and η 2 ) in gene families from all populations subdivided by strategy. Gene families from specialists tend to be more specialized to the gradient than families from generalists ( Fig. S2). In Figure S3, we show examples of how specialists adapt to the gradient through genomic differentiation.
This analysis shows that the genotypic structure of a population recapitulates its eco- evolutionary strategy (see main text).
Several more patterns can be found in the genotypic data. As expected, the number of binding sites upstream of a gene correlates with the number of regulatory inputs (n b ×k inS : Similarly, the sector-specific variation of these properties is also correlated (η 2 n × η 2 k : r = 0.7, R 2 = 0.49, p < 0.0005, N = 247). Interestingly, specialization of a gene's binding sequence is negatively correlated with the number of binding sites that it targets (n b × k inS : r = −0.66, R 2 = 0.44, p < 0.0005, N = 247). Thus, in line with literature, the binding domain of a transcription factor is more constrained by purifying selection if it coevolves with more targets.

Model testing
To assess how different features of our model impact its dynamics, we performed evolution experiments with modifications of the main model (Fig. S4).

Replication cost
First, the constraint on genome size is removed by making the replication cost independent of genome size: the nutrient level then defines the proportion of the genome rather than absolute number of beads that is replicated in one timestep. Under such proportional replication, massive and rapid genome expansion occurs (L > 120 by t = 0.25 · 10 6 ).
Most added beads are likely not functional, simply being retained due to reduced purifying selection on genome size. Still, populations also grow rapidly, showing that, as in the standard model, genome expansion brings about adaptation. The 10 replicates with proportional replication were terminated early because too large genomes become a computational burden.

Stochastic binding
Genome expansion, and thereby adaptation, is entirely absent when stochasticity is removed by setting a binding affinity threshold at H = 17-i.e. the effects of strongly interacting genes at each binding site are averaged, while the rest is ignored. Instead, genomes are streamlined very early from L = 64 to L ≤ 60 so that the short ancestral cell-cycle-which, in the absence of stochasticity, is at the theoretical optimum for high nutrient levels-can be employed in a slightly wider niche range. In 3 out of 10 replicates, cells expanded their range even further by loss of 8 household genes (obtaining L = 52) despite large fitness penalties.

Duplication and deletion rate
When we increase the rate of duplication and deletion 5-fold, genome expansion becomes more prevalent, but the effect is more subtle than in the two previous experiments. The energetic cost of replication still dominates, so the emergence of genomes larger than L ≈ 80 is generally prevented. Still, the greater influx of additional regulatory elements helps to reorganize cell-cycle regulation and thus favours adaptation, as seen in the greater population sizes than the standard model during the entire experiment. There actually appears to exist an optimum for the rate of duplication and deletion, as another 5-fold increase (25fold compared to the standard model) yields smaller populations than the standard model.

Innovation rate
In the main text, we firmly established that genome expansion was linked to adapta-tion. Not only did we show that genome size was correlated with final population size across replicate experiments, but we also revealed the mechanisms behind these adaptations (e.g. rewiring of the network). Still, the implementation of per-genome innovation rates for genes and binding sites could result in a mutational bias towards larger genomes.
Unlike the per-gene duplication rate, which is opposed by an equal per-gene deletion rate, the innovation rate is not balanced. Thus, to test whether genome expansion depends on unopposed innovation rates, we performed an additional evolution experiment without innovation (µ b,in = µ g,in = 0, Fig. S4). Although innovations can be helpful for evolving cells by providing them with new material, 3 out of the 10 replicates of this additional experiment still manage to achieve substantial genome expansion and corresponding range expansion. This supports our previous results (see "Model testing") that mutation rates are not driving genome expansion directly; they are just important suppliers of raw material to selection. Figure S5: The number of household genes is minimised by selection, even if genomes with more than 50 household genes are not directly penalized as in the Standard model. Lines show number of household genes along the ancestral lineage for a single cell from the middle of the gradient (sector 6) at the end of the experiment (t = 2 · 10 6 ). For each of the ten replicates of both models one such random lineage is plotted (Standard and Redundant household; one replicate of each model is shown in bold).

Household gene limit
Finally, the strength of selection can also be demonstrated by the evolution of the number of household genes. Household genes are simple in that they do not have a role in the regulatory network, yet they are a burden for replication just like binding sites and regulatory genes. In the main model, cells that do not have exactly 50 household genes incur a penalty at division. Yet, if we perform experiments with cells that are only penalized for having fewer than 50 household genes-and not for having more than 50-the number of household genes barely increases, from 49.9 to 50.7 (means at t = 2 · 10 6 across ten replicates). As with higher rates of duplication, deletion and innovation, allowing cells to encode redundant household genes grants them more flexibility to change their genome size during evolution (Fig. S5). Still, selection against genome expansion is so strong that the behaviour of the modified model is almost indistinguishable from the main model (Fig.   S4).

Appendix 2: Calculating Fitness
To quantify cell-cycle behaviour, we measured the cell-cycle duration τ and the division rate ρ (fraction of cell-cycles that yield proper division; 1 − ρ being the fraction of cellcycles that yield death). With these two parameters, a simple ODE model representing the cell dynamics of the full model can be constructed. For simplicity, we assume a fixed nutrient density, in line with how τ and ρ of single cells is measured.
The first term on the left hand side represents cell birth, which happens in proportion to the birth rate ( ρ τ ) and the current population size (N ). Furthermore, cell birth only leads to increase in population size if an empty site is chosen for the new cell (1 − N ).
Instead, if the new cell overgrows an existing cell, no change in the population size occurs.
Although overgrowth would affect the actual birth rate by "resetting" cells that may be halfway through their cell-cycle, we ignore this more subtle effect. The second term represents failed cell-cycles (note the 1 − ρ instead of ρ as in the first term). The third term denotes instantaneous cell death by the fixed parameter δ. From the ODE model, the fitness measure R 0 can be calculated as: Figure S7: Evolved cell-cycle behaviour, network topology and genome organization of a single representative of all ten replicate experiments (see main text, Fig. 5). State-space shows only the most frequently visited states. Networks were manually derived from the regulatory dynamics. In terms of cell-cycle behaviour, roughly two clusters appear coinciding with the level of genome expansion; the five replicates on the left retain cell-cycle behaviour similar to the ancestor, whereas the five replicates on the right have greatly enhanced S-loops. At the network level, we can also point out some similarities. R1 and R7 have identical topology (gene names beyond g1-5 are arbitrary). R9 and R8 are identical except that the three inhibitory effects of g7 in R9 are performed by g5 in R8. Finally, R10, R4 and R6 have very similar dynamics, all based on the inactivation of g3 to enter G2. With respect to R10, R6 has only added a single link (g6) upstream of the sensing motif. In R4, as described in the main text, g6 has evolved to directly inhibit g3 (rather than merely hindering its activation by g5) and to be controlled through an AND-gate by inputs from g1 and g4 (see "Evolution of a cell-cycle checkpoint").
9 Figure S8: Adaptation becomes constrained very early in evolution, as shown by population growth of the ten replicate evolutionary experiments plus two alternative trajectories starting at t = 2.5 × 10 5 time steps (lines with same color as original experiment). Alternative trajectories generally end up with similar population size (blue lines high, red lines low), showing that adaptive potential is already constrained by t = 2.5 × 10 5 .