Amphibian Segmentation Clock Models Suggest How Large Genome and Cell Sizes Slow Developmental Rate

Synopsis Evolutionary increases in genome size, cell volume, and nuclear volume have been observed across the tree of life, with positive correlations documented between all three traits. Developmental tempo slows as genomes, nuclei, and cells increase in size, yet the driving mechanisms are poorly understood. To bridge this gap, we use a mathematical model of the somitogenesis clock to link slowed developmental tempo with changes in intra-cellular gene expression kinetics induced by increasing genome size and nuclear volume. We adapt a well-known somitogenesis clock model to two model amphibian species that vary 10-fold in genome size: Xenopus laevis (3.1 Gb) and Ambystoma mexicanum (32 Gb). Based on simulations and backed by analytical derivations, we identify parameter changes originating from increased genome and nuclear size that slow gene expression kinetics. We simulate biological scenarios for which these parameter changes mathematically recapitulate slowed gene expression in A. mexicanum relative to X. laevis, and we consider scenarios for which additional alterations in gene product stability and chromatin packing are necessary. Results suggest that slowed degradation rates as well as changes induced by increasing nuclear volume and intron length, which remain relatively unexplored, are significant drivers of slowed developmental tempo.


Supplement 2: Analytical conditions for the emergence of oscillations
Following the derivation given in Verdugo and Rand [2008], we derive formulae that give analytical conditions for the emergence of oscillations for the system of delayed differential equations (DDE) modelling the segmentation clock as proposed by Lewis [2003].
We begin with Lewis' DDE system, using ṗ and ṁ derivative notation in place of dp dt and dm dt , respectively, for convenience: and we start by rescaling the system.

Rescaling
We rescale, following Verdugo and Rand, and define the following variables: We also use the following notation: m d = m(t − T p ) and p d = p(t − T m ), and we define x d and y d similarly.From the variables defined in Equation (3), we have: Steady state solution At the steady state (x * , y * ) we have x d = x * and y d = y * .Setting ẋ = 0 and ẏ = 0 yields the following steady state equations: x * − by * = 0 (7) Eqn. ( 7) implies x * = by * .Upon substituting in Eqn.(8) we get the following cubic equation for This equation can be solved by Mathematica.It has one real root given by, where, Linearize about a fixed point To linearize about a fixed point we define the following deviations from the steady state (x * , y * ), and we use the subscript d to signify the lagged variable such that and, We now expand Eqn. 12 for small η d .To linear order we get, Note that the second and third terms on the right hand side of Eqn. ( 14) sum to zero.
In order to compare the results in Eqn. ( 14) to Verdogu and Rand, we rewrite the coefficient of η d in terms of β = y * /Y 0 and define: We therefore end up with the following linearized equations: which are analogous to those presented in Verdugo and Rand Oscillatory solution In a supercritical Hopf bifurication, a stable fixed point becomes unstable and is surrounded by a stable limit cycle.Close to the bifurcation, the amplitude of the limit cycle is very small and can be approximated by cosine functions.Therefore, we assume that Eqns.16 and 17 have solutions given by, We can now substitute Eqns.( 18) and ( 19) into the linearized equations, Eqns.( 16) and ( 17), above to solve for ω in terms of model parameters.The condition for oscillation is thus that ω is a real, positive, non-zero number.After some tedious trigonometry and algebra, we find that this condition is satisfied when K > bc.In other words the geometric mean of the degradation constants has an upper bound given by the following.
The frequency of oscillation ω is given by: We may also ask if there exists a minimum total delay, T = T p + T m , required for the emergence of oscillations.Solving for such a T, we get: and we define this T to be a critical total delay, T crit , because it is the time delay required for an oscillatory solution just at the bifurcation.It follows that a second condition for oscillation can be defined, that is: Note that even when the conditions given in ( 20) and ( 23) are met, a finite cut-off in numerical simulations may lead a solution to be classified as non-oscillatory.
In figures S1 and S2, we plot h p against T crit for every model shown in figures 1 and 2 in the main text, respectively (note that the axes are swapped relative to their corresponding plots in the main text, to clearly demonstrate the positive relationship between protein half-life and critical total delay).T crit is dependent on protein and mRNA degradation rates, b and c, respectively, as well as K. K is dependent on β and y * , both of which are dependent on Y 0 which is defined as As a result, we can say that T crit is dependent on only three parameters from the Lewis model: b, c, and p crit .Therefore, we need only generate multiple species-specific plots when h m is changing between models (recall that b = ln(2)/h p and c = ln(2)/h m ) because the range of h p and value of p crit will remain the same across all species-specific models, regardless of diffusion type.This is why there are only two plots in figure S1 (corresponding to species-and diffusion-specific for which mRNA half-life is held constant at h m = 3), but there are six plots in

87
There is a clear and consistent trend across all plots in both figures.That is, as protein half-life increases, the total delay time required for oscillations to emerge, T crit , also increases.Analytical analysis therefore confirms that total delay time must increase with protein stability to ensure the emergence of oscillations  In figure S4, we plot the distribution of nuclear export times for all four diffusion models shown in figure S3 (BM with initial position at the origin, BM with initial positions drawn from a uniform distribution, fBM with initial position at the origin, and fBM with initial positions drawn from a uniform distribution), and for three different nuclear radii (3.5, 6, and 13 µm; most species have nuclear radii well below 6 µm).We can see that nuclear export distributions are skewed more towards the left (i.e.towards smaller times) when initial positions are drawn from a uniform distribution relative to when initial position is always at the origin or nuclear center.
This pattern holds across nuclear size.segmentation clock down, we "break" it altogether by damping oscillations.Additionally, results in figure S5B and S5C suggest that to achieve a period close to 154 minutes protein molecules are either only slightly more stable than their transcripts (S5B), or less stable than their transcripts (S5C).This is an issue discussed in the main text.
Resulting periods of mRNA and protein expression given parameter values corresponding to X. laevis and A. mexicanum Brownian Motion models minutes) a = 4.5, k = 33, p crit = 161, h m = 3, h p = 3

Figure S1 :
Figure S1: T crit plotted across increasing protein stability corresponding to a range of half-lives

Figure S2 :
Figure S2: T crit plotted across increasing protein stability corresponding to a range of half-lives between 3 and 23 minutes.A, C, E A. mexicanum BM models with h m = T exp , h m = 1 2 T exp , and h m = T exp , respectively.B, D, F A. mexicanum fBM models with h m = T exp , h m = 1 2 T exp , and h m = T exp , respectively.All panels correspond with their counterparts from figure 2.

Figure S3 :
Figure S3: Nuclear export simulation results.A Mean nuclear export time across nuclear radii

Figure S5 :
Figure S5: Further increasing gene product stability in the A. mexicanum BM model.A We keep

Figure S6 :
Figure S6: Resulting amplitudes of mRNA expression for A. mexicanum models corresponding to