Amphibian segmentation clock models suggest mechanisms of slowed development across increasing genome size and nuclear volume

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 ten-fold in genome size: Xenopus laevis (3.2 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, which remain relatively unexplored, are significant drivers of slowed developmental tempo.


Introduction
Across the tree of life, genome size and cell volumes span a remarkable range, and a positive 18 correlation has been observed between increases in genome size, cell volume, and nuclear volume (Gregory 2001, Malerba and Marshall 2021, Sessions 2008. However, the mechanisms underlying these relationships and the implications associated with such increases remain areas of ongoing research. For example, evolutionary increases in genome, cell, and nuclear size have been found to slow developmental processes (Jockusch 1997, Sessions and Larson 1987, Wyngaard et al. 2005), but the driving mechanisms are poorly understood. Development emerges from the progression 24 and interaction of a wide range of processes taking place at the single-cell level, where increasing genome, cell, and nuclear size impact transcription dynamics, intra-cellular distances, surface area to volume ratios, and other fundamental characteristics (Cadart et al. 2023, Sessions and 27 Wake 2021). We therefore consider how alterations in single-cell processes might translate to slowed developmental tempo as increasing genome, cell, and nuclear sizes change cell structure and functionality. 30 Vertebrates comprise a large portion of the overall range in genome and cell size across the tree of life, and despite variation in genome size, cell volume, and developmental tempo, many developmental processes remain conserved. Somitogenesis is one such process that is relatively 33 well understood; it is the process through which bilateral pairs of somites, blocks of presomitic mesoderm (PSM) tissue, are patterned along the head-to-tail axis in vertebrate embryos. Somitogenesis is typically described as operating through a clock and wavefront mechanism in which 36 cell-autonomous oscillations of a somitogenesis gene (i.e. the "segmentation clock") interact with Notch, Wnt, FGF, and retinoic acid pathways across the PSM tissue to coordinate proper timing of segmentation of groups of neighboring cells into bilateral pairs of somites (Aulehla and Pourquié 39 2010, Cooke and Zeeman 1976, Gibb et al. 2010, Klepstad and Marcon 2023. The segmentation clock operates via oscillatory gene expression at the single-cell level, while the wavefront takes place at the intercellular level across the PSM. Vertebrates exhibit species-specific segmentation 42 clocks directly related to oscillatory expression of an autoregulated somitogenesis gene at the single cell level (Diaz-Cuadros et al. 2023, Lázaro et al. 2023, Matsuda et al. 2020. One period of oscillation, the time required for one cycle of expression, at the single cell level determines 45 the time needed for a group of neighboring cells to segment off from the larger block of unsegmented PSM tissue, creating a bilateral pair of somites. As a conserved phenomenon that drives developmental tempo while operating at the single-cell level, the segmentation clock is an 48 appropriate lens through which to examine single-cell processes as an underlying link between increasing genome and cell size and slowed development. While species-specific segmentation rates have been linked to biochemical differences at the intracellular level  2023, Lázaro et al. 2023, Matsuda et al. 2020, Rayon et al. 2020, the role of genome and cell size as potential mediators of species-specific gene expression kinetics and therefore developmental tempo have not been explicitly examined. 54 The segmentation clock is well modeled by the following system of delayed differential equations, first proposed by Lewis [2003], that describe coupled oscillatory expression of mRNA and protein, 57 dp dt = am(t − T p ) − bp(t) (1) Similarly, equation (2) models the change in mRNA expression over time, dm dt . cm(t) describes the rate at which mRNA molecules degrade while f (p), given below, is a Hill function that describes 72 the rate at which mRNA molecules emerge from the nucleus into the cytoplasm, where k is a rate constant for mRNA synthesis (mRNA/min/cell) in the absense of repres-75 sion, p crit is the number of protein molecules in the nucleus needed to yield an assumed critical concentration of 10 −9 M within the nucleus (from Lewis [2003]) associated with transcriptional repression, and n is the so called "Hill coefficient". The rate of mRNA emergence is dependent 78 upon and inversely related to protein quantity, indicating the presence of an auto-repressive mechanism (i.e. a negative feedback loop). Like Lewis [2003], we assume that these repressive proteins act as dimers and let n = 2.

81
As touched on above, the parameters T p and T m account for delays associated with protein and mRNA production, respectively. Delay parameters reflect the reality that biological processes are often a non-instantaneous affair. Mathematically, they distinguish the DDE system 84 above from an ordinary differential equation system and are necessary to generate the sustained oscillations that yield species-specific periods of gene expression (Lewis 2003), corresponding to species-specific rates of somite segmentation.

87
Here, we adapt Lewis' model to assess the impact of increasing genome and cell size on the segmentation clock periods associated with two amphibian species that exhibit a 10-fold difference in genome size: Xenopus laevis and Ambystoma mexicanum, the model frog and salamander, 90 respectively. We break delay parameters down into specific transcription, post-transcription, and translation processes, and we consider the potential impact of increasing genome and/or cell size on each individual component. We also adjust critical protein threshold values to reflect species-93 specific nuclear volume estimates. Finally, we consider additional potential roles for mRNA and protein stability, which are not directly related to genome or cell size, in the mediation of developmental tempo. We simulate the Lewis model under all of these different scenarios to test 96 whether we can reproduce the observed periodicity of the somitogenesis clock. We buttress our simulations with an analytical derivation of the minimal conditions for oscillations in the Lewis model. With our approach, we are able to establish direct links between increases in genome size 99 and nuclear volume and their specific impacts on gene expression (i.e. increases in network delays and threshold value) that yield slower developmental tempo, and we are also able to assess potential indirect roles for gene product stability in the mediation of developmental tempo.

Taxon selection
We choose to adapt the Lewis model to X. laevis and A. mexicanum. Their recorded genome sizes 105 are ∼3.1 Gb (X. laevis) and ∼32 Gb (A. mexicanum) (Hellsten et al. 2010, Smith et al. 2019, and, although volumes vary across type, A. mexicanum cell volumes are typically larger than their X. laevis counterparts. A. mexicanum nerve cells, for example, are ∼2-times larger than X. laevis 108 nerve cells (Roth and Walkowiak 2015), and their red blood cells (RBCs) are ∼10-times than in X.
laevis (Gregory 2023, using cell volume measurements for Ambystmoma tigrinum whose average reported genome size is also ∼32 Gb. Note also that amphibian RBCs are nucleated). Meanwhile, 111 there is about a 3-fold difference in the rate of somite segmentation. In X. laevis, bilateral pairs of somites are segmented off from the PSM every 50 minutes (extrapolated by Curran et al. [2014] from Faber and Nieuwkoop [1994] and Hamilton [1969]); in A. mexicanum, somite segmentation 114 occurs every ∼155 minutes (Armstrong and Graveson 1988).

Generating species-specific parameter values
Our first goal was to test if parameter changes directly related to increasing genome, cell, and 117 nucleus size are sufficient to recapitulate the slowed rate of somite segmentation in A. mexicanum relative to X. laevis. To this end, we start by generating species-specific delay time and critical parameters (rates of production and degradation) at constant values or ranges between species.
We derive protein and mRNA-production delays, T p and T m , by applying estimation methods from Lewis to somitogenesis gene candidates for X. laevis and A. mexicanum. Vertebrate clock 123 genes are members of the Hairy and enhancer of Split (Hes/Her) family of basic helix-loop-helix (bHLH) genes (Kageyama et al. 2007). Hes/Her gene family size varies across vertebrates, and the individual cycling members of the clock network also vary (Eckalbar et al. 2012, Krol et al. 126 2011). In zebrafish and mice, hes7 has been shown to be the central component of the oscillator (Bessho et al. 2001, Hirata et al. 2004, Holley et al. 2002, Oates and Ho 2002, and its oscillatory pattern has also been shown in the PSM in the lizard Anolis carolinensis (Eckalbar et al. 2012).
We therefore infer that hes7 is the ancestral clock for vertebrates and that it retains clock function in A. mexicanum. In Xenopus laevis, in contrast, the Hes/Her gene family is greatly expanded to 37 copies, following both ancient allotetraploidization and tandem duplication ( Kuretani et al. 132 2021, Watanabe et al. 2017). The hes7 orthologs in X. laevis are not cyclically expressed in the PSM and, therefore, cannot act as the clock (Davis et al. 2001, Jen et al. 1999, Shinga et al. 2001).
In contrast, in Xenopus laevis, three Hes/Her family genes are known to oscillate in the PSM: 135 hes5.3, hes5.5, and hes5.7 (Blewitt 2009, Li et al. 2003. Of these, the strongest candidate is hes5.7L based on experimental data showing that de novo protein synthesis is required to repress hes5.7L transcription during somite formation, and that hes5.7L RNA instability is part of the mechanism 138 underlying its cyclic expression (Davis et al. 2001, Li et al. 2003. Hes5.7 is absent from the closely related X. tropicalis and is thus inferred to be specific to X. laevis (Kuretani et al. [2021], Watanabe et al. [2017]), suggesting a case of developmental system drift (Haag and True 2021).

141
In X. laevis, hes5.7L has a primary sequence length of 1,604 nt; it is made up of 3 exons and 2 introns (lengths: 166 and 113 nt), and its coding sequence is 465 nt (Sayers et al. 2022 Critical protein threshold, p crit We calculate the number of protein molecules needed to achieve a critical concentration of 10 −9 M in the nucleus, based on species-specific nuclear volumes. We assume a spherical nucleus (V = 4 3 πr 3 ) and estimate species-specific radii of 4 and 5.5 Protein production delay, T p We assume that the delay associated with protein production is 153 equal to translation delay. We estimate species-specific translation delays by applying a translation rate of 6 nucleotides per second (Lewis 2003) to the reported coding sequence lengths. mRNA production delay, T m We consider mRNA-production delay to be a cumulative sum 156 of transcription T tx , intron-splicing T in , and mRNA export T exp delays. We estimate speciesspecific T tx values by applying a transcription rate of 20 nucleotides per second (Lewis 2003) to the reported primary sequence lengths. Hoyle and Ish-Horowicz [2013] find that in vivo intron-159 splicing delay constitutes a relatively constant proportion of ∼8.3 % of the overall segmentation clock period in mice, chick, and zebrafish embryos. We apply this proportion to the reported clock periods in X. laevis (∼50 min) and A. mexicanum (∼155 min) to get species-specific T in 162 values. We estimated species-specific mRNA export delays, T exp , using simulations of particle diffusion within a sphere, described in detail below.
Simulation of mRNA export time parameter, T exp Before they are released into the cytoplasm, To generate species-specific estimates for nuclear export, we simulate the (3D) random walk of a diffusing particle within spheres of radii between 3 and 6 µm. This range of radii captures measurements for X. laevis and A. mexicanum PSM cell nuclei, 4 and 5.5 µm, respectively. We simulate mRNA transcript trajectories and record the number of steps needed for our simulated mRNA molecule to first cross the nuclear periphery from the nuclear center. 10,000 trajectories are simulated for each sphere in the specified range (radius of 3 to 6 with intervals of 0.5).

186
The transcript trajectory for normal diffusion is simulated by generating x, y, and z position vectors as cumulative sums of increments (i.e. step sizes) chosen from a normal random distribution. The transcript trajectory for obstructed diffusion is simulated similarly, with the x, y, 189 and z position vectors generated directly by a fractional Brownian Motion function in MATLAB, wfbm, with a Hurst parameter of 0.25. For both simulations, the number of steps required by a particle trajectory to first exit the domain is retrieved and averaged. We assume that each step 192 takes one second, and we scale the particle trajectory such that the estimated mean first exit times agree reasonably well with mRNA exit times observed for the somitogenesis gene in Danio rerio, or zebrafish. It has been reported that nuclear export of her1(hes7), the somitogenesis clock We assume the nuclear center as the initial mRNA position based on known chromosomal 201 territories associated with the somitogenesis gene in humans and mice coupled with patterns of synteny observed across vertebrates. That is, in both mice and humans, the somitogenesis gene, hes7, is found on gene-rich chromosomes 11 and 17 (which are homologs), respectively; both 204 chromosomes localize in their respective nuclear centers Boyle et al. 2001, Kile et al. 2003, Mayer et al. 2005, Zody et al. 2006. Broad patterns of synteny conservation and topologically associated domain conservation have been observed across vertebrates (Schloissnig et al. 2021, Smith et al. 207 2019), suggesting that this pattern extends beyond humans and mice.

Parameters held constant
For the first set of analyses, we hold the rates of mRNA production (in the absence of inhibition), protein production, and mRNA degradation constant, at k = 33 210 mRNA/min, a = 4.5 protein/mRNA/min, and c = ln(2)/3 mRNA/min corresponding to a halflife h m = 3 minutes (Lewis 2003); and we consider a set range of protein stability/degradation ln(2)/23 < b < ln(2)/3 protein/min, corresponding to a half-life of 3 < h p < 23 minutes. This 213 range is chosen based on typical reported and estimated protein stability of the somitogenesis gene across model vertebrates, namely zebrafish and mice; we test the model across a range of protein half-lives due to its important role in mediating the period of gene expression (Hirata 216 et al. 2004, Lázaro et al. 2023, Lewis 2003, Matsuda et al. 2020, Takashima et al. 2011) (see also Extrapolation of periodicity to compare with known somite segmentation rates 219 For our first set of analyses, we plug the parameter values described above back into the Lewis model, and we use the DDE solver ddesd in Matlab to generate solutions across identical ranges of protein stability h p and species-and diffusion-specific ranges of total delay time, T m + T p . We 222 assess the period of gene expression that emerges for each set of solutions, and we compare it to the known somite segmentation rate of the corresponding species. In doing so, we aim to first verify that our parameter selection does in fact yield the correct period of oscillation for X. laevis, 225 and to then determine if parameter changes directly driven by genome and cell size differences are sufficient to capture slowed developmental rate in A. mexicanum relative to X. laevis. To assess periodicity, we create vectors to store the local extrema (i.e. local minimum and 228 maximum values) and corresponding time stamps for each solution. Gene expression tends to spike in the first 4 to 5 cycles of oscillation before settling into a long-term pattern, so we remove the first 5 cycles of oscillation from our data to avoid skewing. The period of oscillatory 231 gene expression is calculated by taking the average difference between successive time stamps associated with local minima (using local maxima would yield the same results). We do this across a time span of 0 to 3,100 minutes. The time span, 3,100 minutes, is chosen based on  (Dali et al. 2002), each taking ∼50 minutes to form (Li et al. 237 2003). For both species, oscillations must remain robust throughout this time span. We define robust oscillations as having an amplitude (the height of an oscillation, or difference between local minimum and maximum) of no fewer than 10 molecules throughout the time span. If this 240 requirement is not met, the oscillations are considered damped and we define the periodicity as Inf (infinite). We choose 10 molecules as a conservative finite cut-off based on observed average RNA transcript amplitudes of zebrafish segmentation clock genes her1 and her7 which are ∼41 243 and ∼49 molecules, respectively (Keskin et al. 2018).
The assessment of periodicity described above is done for both mRNA and protein counts, and the difference in periodicity is always within 0.2 minutes (see Supplemental Material 1).

246
In other words, results are relatively similar for both sets of oscillations, so we choose mRNA periodicity to represent overall system behavior.
Calculating average amplitude of expression 249 We also calculate the average amplitude of mRNA expression across the parameter combinations in each model. To do this, we use the extrema vectors described above to create a vector that stores the difference between local minima and maxima corresponding to every complete cycle of 252 oscillation, excluding data from the first 5 oscillation cycles that were cut out. The resulting vector gives the amplitude associated with each complete oscillation cycle, and we take the average of all vector entries to get an average amplitude associated with a particular parameter combination.

255
For all combinations with periodicity defined as Inf (infinite), we set amplitude equal to 0, to keep results consistent with each other.

258
To assess the impact of individual changes on the period of gene expression, we increase and decrease each individual parameter (assessing total delay as a parameter as opposed to taking T p and T m individually) by 50% while holding all other parameters constant. We first extrap-261 olate the resulting period and amplitude of gene expression for an original set of parameters (corresponding to an A. mexicanum Brownian Motion model with protein half-life arbitrarily set at h p = 15 minutes) and then for each parameter change using the methods described above.

264
Testing whether scaling mRNA stability with export time yields the rate of somite segmentation in A. mexicanum Our mRNA export simulations suggest that transcripts take much longer to leave larger nuclei.

267
We therefore assess the impact of increasing mRNA stability with estimated nuclear export time on the A. mexicanum segmentation clock period. To this end, we re-considered both A. mexicanum models, normal and fractional Brownian Motion corresponding to normal and obstructed diffu-  Generalizing mRNA export simulations 276 We generalize our mRNA export simulations by running the simulations described above across a range that encompasses what has been observed across the tree of life, that is nuclei of radius between ∼0.5 and 13 µm, based on the minimum and maximum nuclear volumes reported in while assuming a spherical volume V = 4 3 πr 3 . We also run simulations for which the initial position is selected from a uniform distribution as opposed to always being set at the origin 282 (assumed to be at the nuclear center). The uniform distribution we draw our initial x, y, and z positions from encompasses a domain that is 3 4 of each radius. In doing so, we allow for initial positions to be drawn from positions throughout our theoretical nuclei, excluding the periphery. 285 We choose to exclude the nuclear periphery because this is where heterochromatin is spatially concentrated and reduced transcriptional activity has been observed (Bizhanova and Kaufman

309
Parameter changes directly linked to increasing genome and cell size can mathematically recapitulate slowed developmental tempo Using simulated export times, we generate species-and diffusion-specific general delay times 312 associated with mRNA production, a sum of transcription, intron splicing, and nuclear export delays. All other parameters are held constant either across species (a, k, b, c) or at species-specific values across diffusion models (p crit , T p ). All of these values are shown in Table 2.

315
The resulting periods of oscillation for each model are shown in figure 1 (with non-oscillatory combinations/regions, assigned period Inf, shown in dark purple). We test across a range of parameter combinations: on the x-axis, we have values of total delay, T m + T p , that fall within 318 ±5 minutes of our species-and diffusion-specific values given in Table 2; on the y-axis, we have a range of protein stability corresponding to 3 < h p < 23 minutes. The observed rate of somite segmentation in X. laevis (∼50 minutes) is captured by a subset of total delay and protein stability is outlined in the plot by a dashed-line in figures 1A and 1B. In confirming that we can achieve the correct period of oscillation for X. laevis, the results in figures 1A and 1B provide support for our 324 methods of parameter estimation; these results also act as a plausible baseline against which the A. mexicanum models can be compared. Meanwhile, the observed rate of somite segmentation in A. mexicanum (∼155 min) is only captured under sub-diffusive conditions. Genome and nucleus 327 size-driven increases in delay time and concentration threshold are sufficient to fully recapitulate slowed development in A. mexicanum when nucleoplasmic movement of transcripts is assumed to be sub-diffusive ( figure 1D), but are insufficient when normal diffusion is assumed ( figure 1C).

330
When normal diffusion is assumed, genome and cell size-driven parameter changes can slow the period of oscillatory gene expression down to ∼125 minutes, 30 minutes faster than the known rate of somite segmentation. The additional delay introduced by assuming sub-diffusion, about 333 14 minutes longer, is needed to produce a set of total delay and protein stability combinations that yield a ∼155 min period of oscillation.
mRNA export is 2-to 3-fold slower in A. mexicanum relative to X. laevis 336 mRNA export time roughly doubles in A. mexicanum relative to X. laevis when normal diffusion is assumed, and roughly triples when obstructed diffusion is assumed (see T exp in Table 2). The impact of obstructed diffusion on export time becomes more pronounced (i.e. the gap between 339 export times for the two diffusion types becomes wider) as nuclear radius and therefore volume increase (see figure 3).
Scaling mRNA stability with nuclear export time yields biologically plausible 342 recapitulation of slowed developmental tempo Although the fBM model shown in figure 1C yields the known rate of somite segmentation in A.
mexicanum based solely on genome and cell size differences in parameter values, an additional 345 increase in mRNA stability seems logically necessary given our simulated nuclear export times. mRNA half-life is held constant at: A, B diffusion-specific estimates for mRNA export delay; C, D half of estimated mRNA export delays; E, F a quarter of estimated mRNA export delays. * Note: Color of star/outline is chosen for contrast and has no additional meaning.
The models shown in figure 1 assume an mRNA degradation rate associated with a half-life of 3 minutes. However, we are working with simulated mean export times of ∼12 and ∼26 minutes in A. mexicanum PSM nuclei under normal and sub-diffusive conditions, respectively. Under these assumptions, a vast majority of mRNA molecules are expected to degrade long before ever leaving the nucleus. 351 We therefore test if the A. mexicanum segmentation clock can also be recovered when mRNA stability is increased relative to export time such that a greater proportion of transcripts are able to exit the nucleus before degrading. Although it is established that some fraction of RNA 354 transcripts will degrade in the nucleus before exiting into the cytoplasm, the extent of degradation seems to differ across transcript types and remains relatively understudied (Smalec et al. 2022 To account for the wide range of empirical and theoretical estimates of mRNA degradation in the nucleus, we explore 3 potential scaling scenarios for mRNA stability and nuclear export time 369 in A. mexicanum. We first set mRNA half-life equal to diffusion-specific export times, h m = T exp , for which 50% of mRNA transcripts degrade before leaving the nucleus (and 50% leave the nucleus before degrading); then we set mRNA half-life equal to 50% of the simulated export 372 times, h m = 1 2 T exp , for which 75% of mRNA transcripts degrade before leaving the nucleus (and 25% leave before degrading); finally, we set mRNA half-life equal to 25% of the simulated export times, h m = 1 2 T exp , for which 87.5% of mRNA transcripts degrade before leaving the nucleus (and 375 12.5% leave before degrading).
In Figure

384
Previous research on developmental timing has pointed towards intron length and differences in biochemical characteristics, like degradation rates and network delays, as sources of speciesspecific tempo (Lázaro et al. 2023, Matsuda et al. 2020, Rayon et al. 2020 2008, . In our study, we not only reconsider these points while contextualizing them within the scope of increasing genome and cell size, but we also make spatially explicit considerations that set our study apart from others.

390
In Table 2, we outline how species-specific parameters differ between X. laevis and A. mexicanum, and we note whether each change (based on our methods of estimation) captures increases in genome size or nuclear volume, or neither. Using these species-and diffusion-specific 393 parameters, we first verified that our parameter estimates for X. laevis were able to recapitulate the ∼50 minute segmentation clock. This indicated that our methods of parameter estimation were reasonable, and we moved forward using the same methods for our A. mexicanum models, 396 testing for parameter changes and combinations that recapitulate the ∼155 minute segmentation clock. The results of our models reveal the following potential mechanisms through which developmental tempo slows with increased genome and nuclear size.

399
Increasing intron length impacts developmental tempo through transcriptional delays Increasing the sum of delays in the gene regulatory network has the most significant impact on 402 the period of gene expression (see Table 3: Sensitivity Analysis). When we break total delay down into individual components, we see that while transcriptional delay, T tx , is not the most significant contributor to total delay, T m + T p , it confirms an intuitive link between increasing 405 genome size and the alteration of oscillatory gene expression kinetics. This is because transcriptional delay includes the time needed to transcribe the entire primary gene sequence, including intronic regions, which scale positively with genome size. Indeed, increasing intron size (in part 408 due to transposable element insertions) is a major driver of genomic expansion in A. mexicanum (Nowoshilow et al. 2018), and while intronic regions only constitute ∼17% of the hes5.7L primary sequence in X. laevis (279 bp of intronic sequence), they make up ∼76% of the hes7 primary se-411 quence in A. mexicanum (6,307 bp of intronic sequence). This is a general pattern observed across orthologous genes in A. mexicanum and X. laevis (Nowoshilow et al. 2018). Interestingly, this pattern of intronic expansion seems to be constrained in developmental genes, where an ∼11-414 fold increase in average intron length is observed between X. laevis and A. mexicanum in contrast with an almost 20-fold increase seen in non-developmental genes (Nowoshilow et al. 2018). In Nowoshilow et al. [2018], the role of transcription in developmental timing is suggested as an 417 explanation for this pattern. The time required to transcribe intronic regions of genes (i.e. "intron delay") has not only been identified as a potential mediator of developmental timing via its role in gene regulatory networks , but has also been suggested to link 420 genomic gigantism in particular with slowed development and regenerative abilities through its impact on gene expression kinetics (Sessions and Wake 2021).
However, according to our computational results, increased transcriptional delay alone cannot 423 drive the vast differences in overall delay time between X. laevis and A. mexicanum. Our model therefore also implies that intron delay alone cannot drive slowed developmental tempo observed in the A. mexicanum segmentation clock relative to that of X. laevis. Instead, we see in Table 2 that 426 the delay associated with increasing nuclear volume, T exp is the most significant contributor to overall delay time across all species-and diffusion-specific models. Although the delay associated with intron-splicing, T in is an important contributor to total delay, and the difference between 429 species-specific values for T in is in fact larger than that between species-specific values for T tx , inter-species differences in splicing kinetics cannot be clearly linked with differences in intron or, by extension, genome size (Khodor et al. 2012).

432
mRNA export time increases with nuclear volume and has a pronounced impact on developmental tempo Increasing nuclear volume between X. laevis and A. mexicanum is captured by the critical protein 435 threshold value, p crit , and mRNA export delay, T exp , parameters. The increase in critical protein threshold accounts for the fact that in a larger nucleus, more protein molecules are needed to reach the same critical concentration of 10 −9 M required for transcriptional repression to act "in 438 earnest" (Lewis 2003). However, changes to p crit have a very limited impact on the period of gene expression relative to other parameter changes. Meanwhile, the increase in mRNA export delay in A. mexicanum relative to X. laevis has a relatively pronounced impact on the period of gene 441 expression.
Although the increase in the estimated radius of PSM nuclei in A. mexicanum relative to X.
laevis may seem minimal, 5.5µm compared to 4µm, mRNA export delay still nearly doubles in 444 A. mexicanum when normal diffusion is assumed and nearly triples when obstructed diffusion is assumed. The impact of increasing radius is greater when obstructed diffusion is assumed, and differences in estimated export times (between lengths and diffusion types) become more pro-447 nounced as radius increases. As a result, we would expect mRNA export delays to become more pronounced in gene regulatory networks across increasing nuclear volume, and we might also expect mRNA export to be the largest contributor to total delay in a gene regulatory network, ). An initial position at the nuclear center is assumed. Trajectories are scaled such that a radius of 3µm corresponds to a mean export time of ∼3.36 minutes (shown by the red arrow) to match the reported export time of her1 (hes7) in zebrafish Hoyle and Ish-Horowicz 2013. * Note: There is a limit to how closely 3µm can be scaled with ∼3.36 minutes. As a result, obstructed diffusion mean export times start off slightly faster than normal diffusion when radius r < 3.5µm, but this is not necessarily biologically meaningful, and mean export times are quick to converge back to expectations. especially in larger nuclei. This can be seen in figure 3, where we plot simulated mean export times under both normal and obstructed diffusion for radii between 3 and 6 µm. If we expand the range of radii past 6 µm, the impact of increasing radius and obstructed diffusion on ex-453 port time becomes even more pronounced, suggesting that the nuclear center in species with the largest known genome sizes may have become "uninhabitable" for genes with dynamic expression patterns. If we introduce transcripts whose positions are drawn from a uniform distribution 456 within the nucleus, as opposed to always starting from the origin/nuclear center, the impacts of increasing radius and obstructed diffusion on mean export time are reduced. The distribution of export times becomes more skewed towards lower values (because most transcript trajectories 459 are starting closer to the periphery) (see Supplemental Material 3). In very large nuclei, we hypothesize that genes -especially those with dynamic expression patterns -may be constrained to occupy these locations away from the nuclear center, suggesting an overall effect of genome 462 expansion on the organization of chromosomal territories.
Regardless of model type, assuming an intial position at the nuclear center, our simulations position mRNA export as the largest contributor to total delay time. This is consistent with 465 observations that find mRNA export delays to be longer than both transcriptional and intron splicing delays in mice, chick, and zebrafish segmentation clocks (Hoyle and Ish-Horowicz 2013).
As the largest contributor to total delay time, mRNA export is also the largest contributor to 468 overall differences in X. laevis and A. mexicanum segmentation clocks, implying that spatially induced delays in gene regulatory networks play a significant role in the slowing of cellular and developmental processes. This suggests an important albeit underexplored spatial component in With these empirical possibilities in mind, we turn towards our theoretical results. In figures

486
1A and 1B, we see that both diffusion models are able to fully recapitulate the X. laevis segmentation clock, while obstructed diffusion (modeled by fractional Brownian Motion, shown in figure   1D) must be assumed to fully recapitulate the slowed A. mexicanum segmentation clock. Simi-489 larly, in figure 2, we have that obstructed diffusion still must be assumed, regardless of mRNA gene expression period of 155 minutes. This is in agreement with previous studies that point to gene product stability as a mediator of species-specific segmentation clock periods (Hirata et al. 2004, Lázaro et al. 2023, Matsuda et al. 2020, Rayon et al. 2020). However, some mRNA stability 507 scenarios yield more reasonable results than others. In figures 2B and 2D, with mRNA half-lives of 26.27 and 13.14 minutes, the 155 minute period of oscillation is captured for protein half-lives that range from 3 to 8 and 3 to 13 minutes, respectively. As a result, we would have to assume that 510 protein is as stable or less stable than mRNA, which is inconsistent with reported patterns for Hes7 in mice and humans (Matsuda et al. 2020). In figure 2F, for which mRNA half-life is set to 6.54 minutes, we get a larger subset of parameter combinations that yield a period of 155 minutes.

513
Across this subset, protein half-lives range from being only slightly more stable than mRNA at ∼8 minutes to almost 4-times as stable as mRNA at ∼23 minutes. The resulting relationship between mRNA and protein stability, that protein molecules exhibit increased stability relative to 516 their transcripts, is more consistent with the patterns reported in Matsuda et al. [2020]. However, setting h m = 1 4 T exp implies that a vast majority of mRNA transcripts, 87.5%, will degrade before exiting the nucleus. We might expect amplitude (i.e. number of transcripts in the cell) to be 519 extremely low as a result, but this is not the case. In fact, mRNA degradation is related in complicated ways to not only the amplitude of mRNA expression, but to the amplitude of protein expression as well, resulting in counter-intuitive patterns of expression across different mRNA 522 stability scenarios (see Supplemental Material 5).
Considering again the period of gene expression, an interesting tradeoff exists between the possible mRNA stability scenarios described above: in order to mathematically recapitulate the 525 slowed A. mexicanum segmentation clock, we must either assume that a vast majority of mRNA transcripts degrade within the nucleus, or that the protein molecules in this system are less stable relative to their transcripts. The assumption of obstructed diffusion necessitates some kind of 528 deviation -either high rates of nucleoplasmic degradation, or low stability of protein molecules relative to their transcripts -from conventional gene expression kinetics as genome and nuclear size increase between X. laevis and A. mexicanum. Based on reported patterns for Hes7 in mice 531 and humans (Matsuda et al. 2020), we propose that widespread nucleoplasmic degradation of the somitogenesis transcript is more plausible than transcripts that are more stable than the corresponding protein molecules.

534
What remains the same across genome size and nuclear volume is that, as shown in previous implementations of the Lewis model (Hirata et al. 2004, Matsuda et al. 2020), the relative stabilities of different gene products interact to shape clock tempo. This is demonstrated by to mediate developmental tempo, our model results suggest that across increasing genome size and nuclear volume, gene product stability does not act alone to slow developmental tempo.

Conclusion
Our results suggest that spatial differences must be accounted for to recapitulate slowed development as a result of evolutionary increases in genome and cell size. In agreement with previous 549 studies, our results also suggest an important role for gene product stability when it comes to mediating species-specific rates of development. Taken together, we show that the physical and spatial delays predicted by increased intron length and nuclear size, coupled with alterations to 552 gene product stability that ensure the products persist for long enough to fulfill their molecular function, mathematically recapitulate the slow developmental tempo found in species with large genomes. togenesis in the anole lizard and alligator reveals evolutionary convergence and divergence in