An explanatory evo-devo model for the developmental hourglass

Developmental gene execution networks

As a first-order approximation, a regulatory gene can be modeled in one of several discrete functional states²⁴. In the simplest case, a regulatory gene can act as a binary switch (“on” and “off”) but in general a gene may have more than two functional states. The transition of a regulatory gene X from one functional state to another is often (but not always) caused by one or more upstream regulators of X that go through a functional state transition before X. We use the term transitioning gene to refer to a regulatory gene that goes through a functional state transition at a given developmental time anywhere at the developing embryo.

A DGEN is a directed and acyclic network; see Figure 1a for an abstract example. The vertical direction refers to developmental time, from the zygote at the top to the developed organism at the bottom. In the horizontal direction we can represent different spatial domains, even though this is not necessary and it is not done in our model. For instance, the zygote at the top of the DGEN would be a single domain, while the organism at the bottom of the DGEN would have the largest number of spatial domains.

Development is often approximated (conceptually and experimentally) as a succession of discrete developmental stages. The duration of a developmental stage can be thought of as the typical time that is required for a gene’s functional state transition, and it does not need to be the same for all stages. Each layer of a DGEN refers to a developmental stage, and it includes only the transitioning genes during that stage anywhere in the embryo. The same gene can appear in more than one stage if it goes through several functional state transitions during development. Additionally, a DGEN edge from a gene X at stage l to a gene Y at stage l + 1 implies that the functional transition of X caused the functional transition of Y at the next stage. If gene Y has more than one incoming edge, its functional state transition was caused by the coupled effect of more than one transitioning genes at the earlier stage. Any upstream regulators of Y that remained at the same functional state at stage l are not included in that stage of the DGEN.

The sequence of developmental stages is denoted by l = 1 . . . L. The set of transitioning genes at stage l is G(l). A gene g at stage l<L regulates a set of downstream genes at stage l+1 denoted by D(g) (outgoing edges from g). Similarly, a gene g at stage l>1 is regulated by a set of upstream genes U(g) at stage l-1 (incoming edges to g). The functional transitions at the first stage are assumed to be caused by regulatory maternal genes that initiate the developmental process.

Model description

The model captures certain aspects of both the developmental process, in the form of a given DGEN for each embryo, and of the evolutionary process, as random perturbations in the structure of individual DGENs in a population. The model does not need to capture the actual functional state transitions or the regulatory input function of each gene. It does capture however the dynamic and stochastic effect of structural network perturbations (gene deletion, rewiring and duplication) on the success of the developmental process, as explained in the following.

Similar to the Wright-Fisher model²⁵, we consider a population of N individuals, each represented by a DGEN. In each generation, individuals reproduce asexually, inheriting the DGEN of their parent. Various evolutionary events can cause structural changes in the DGEN of an individual that may result in “developmental failure”. Such individuals (and their DGENs) are replaced with developmentally successful individuals so that the population size remains constant.

The model is meant to be as general as possible and so the regulatory interactions between genes of successive stages are determined probabilistically, as follows. Each stage l is assigned a regulatory specificity, or simply specificity s(l) with 0 ≤ s(l) ≤ 1. A gene g at stage l acts as upstream regulator for a gene g′ at stage l + 1 with probability s′(l) = 1 − s(l). So, the specificity of a developmental stage determines how likely it is for regulatory genes of that stage to cause a state transition of the genes at the next stage; a higher specificity decreases that likelihood.

Our major assumption is that the regulatory specificity increases substantially as development progresses. In other words, the DGEN becomes gradually sparser along the developmental time axis, with s(1) << s(L). This assumption is plausible for the following reasons. First, as development progresses the embryo grows in size forming distinct spatial domains. So, extracellular gene regulation becomes more difficult, especially across different domains. Additionally, as development progresses the transitioning genes become more organ- or tissue-specific, implying that their downstream interactions become sparser. Unfortunately, a direct and empirical investigation of the increasing specificity assumption requires knowledge of the complete DGEN for a given species; this is currently not feasible for even the most well-studied model organisms. However, this assumption is plausible if we consider the well established patterns of increasing modularity as the embryo develops at the morphological, signaling, and genomic levels²³.

The DGEN structural changes we consider are gene deletions, gene duplications, and gene rewiring:

Deletions (DL): This event removes a gene from the DGEN, including its incoming and outgoing edges. There are many genetic mechanisms that may cause such events. A DL event deletes each gene of an individual and at each generation with probability P_DL.

Duplication (DP): This event creates an identical copy of a gene g with the same downstream and upstream regulators and at the same developmental stage as g. The two genes may have different fates if one of them is subject later to deletion or rewiring. Otherwise, the two genes are considered identical. A DP event duplicates each gene of an individual and at each generation with probability P_DP.

Rewiring (RW): This event changes the upstream and/or downstream regulators of a gene. Changes in the upstream versus downstream regulators may have different biological basis. The former occur, for instance, as a result of mutations in the transcription factor binding sites in a gene’s promoter or mutations in distal regulatory elements such as enhancers, while the latter may be mostly caused by coding sequence mutations. The details of the rewiring process do not affect the results qualitatively as long as the average density of edges in each stage remains consistent with the specificity of that stage. The specific rewiring mechanisms we use are presented next.

Suppose that a RW event affects gene g at stage l. The upstream regulators of g are recomputed based on the specificity of the previous stage, i.e., by choosing each distinct gene of stage l −1 with probability s′(l −1). For the downstream regulators of g, we randomly remove N_– existing outgoing edges of g, and then add N₊ outgoing edges to randomly chosen genes of stage l + 1 that g is not already connected to. If B(N, p) denotes a Binomial random variable with N trials and success probability p, the random variables N₋ and N₊ are independent and they both follow the B(|D(g)|, s′(l)) distribution (|X| denotes the cardinality of set X). Thus, the number of downstream edges of gene g after a Rewiring event becomes:

              |D_new(g)| = |D(g)| − N_–+ N₊,

which is at least 0 and at most 2 × |D(g)|. This captures that the downstream regulators of g are derived by incremental changes in D(g), instead of giving g a completely new network configuration (thereby, new regulatory function). The higher the regulatory specificity of a stage, the less likely these incremental changes are. An RW event rewires each gene of an individual and at each generation with probability P_RW.

A gene deletion or rewiring event at stage l can remove an upstream regulator from genes at stage l + 1. A loss of incoming edges may trigger the regulatory failure of a gene, as described next.

Regulatory failures (RF): A gene g may not be able to change functional state if some of its upstream regulators U(g) are lost due to DL or RW events (see Figure 1b). Even though regulatory networks are often robust to structural perturbations, even a partial gene loss in U(g) may disable g causing a regulatory failure. It is plausible that the probability of a regulatory failure increases with the fraction of lost upstream regulators. So, if U′(g) is the new set of upstream regulators and |U(g)| > |U′(g)| > 0, gene g is removed with probability:

Figure 1. An abstract DGEN.

(A) The circles denote state-transitioning genes, edges represent directed regulatory interactions, and colored boxes refer to spatial domains that form during development. If regulatory genes become increasingly function-specific as development progresses, the network gradually becomes sparser in that direction. (B) Evolutionary perturbations on a DGEN’s structure: Gene A is deleted (DL), while gene B is rewired (RW) losing an outgoing edge. This RW event causes the regulatory failure (RF) of gene C, which then causes a cascade of five more RF events. This cascade causes developmental failure (DF). Note that the removal of some upstream regulators does not always cause an RF event (e.g., genes regulated by A).

while if |U′(g)| = 0 we set P_RF (1) = 1. z is the RF parameter and it depends on the robustness of regulatory interactions to gene loss (see Figure 2).

Figure 2. Probability of Regulatory Failure (RF) for three values of the parameter z.

r is the fraction of upstream regulating edges that are lost due to a DL or RW event.

When a DL or RW event causes one or more RF events, the latter can trigger additional RF events in subsequent developmental stages, leading to cascades of regulatory failures. Such RF cascades may cause developmental failure, meaning that the developed embryo is unable to survive or reproduce.

Developmental failure (DF): The last stage of a DGEN represents the fully developed embryo. If that stage includes Γ transitioning genes at a successfully developed embryo, the simplest assumption is that an individual with less than Γ genes at stage-L has failed to develop properly. Such DGENs are removed from the population and they are replaced with randomly chosen but successfully developed DGENs. We have also experimented with two variations of the DF event: first, the individual is removed if its last stage has less than Γ − γ genes, where γ is small relative to Γ, and second, the probability of a DF event increases as the number of genes at stage-L decreases below Γ. The qualitative results, as described next, do not change with these two model variations.

Simulation

We simulate the presented model to examine the properties of the surviving DGENs as evolutionary time progresses. The initial population consists of N identical DGENs with Γ genes at each stage. The edges between genes are constructed probabilistically based on the specificity of each stage, as described previously. Simulating the complete model would not show the significance of individual mechanisms such as the increasing specificity assumption. For this reason we construct a sequence of four models with increasing complexity, presenting results separately for each of them:

Model-1: Constant specificity. Each stage has the same specificity, s(l) = 0.5 for l = 1 . . . L − 1. Further, this model does not include gene deletion and duplication. Gene rewiring can cause RF and DF events even if there are no DL or DP events. We have also experimented with other values of the specificity probability, and the results are qualitative the same (namely, when the specificity is constant there is no hourglass effect).

Model-2: Increasing specificity. The difference from Model-1 is that the specificity is gradually increasing across developmental stages. Unless noted otherwise, the specificity is linear, s(l) = l/L for l = 1 . . . (L − 1); a nonlinear specificity function is also considered, which we describe later. We do not claim that these particular specificity functions are realistic or that there are experimental results that suggest them. They are just the simplest and most parsimonious models that lead to the emergence of the hourglass pattern.

Model-3: With gene duplications. Model-3 adds DP events in Model-2. The duplication probability P_DP is set so that the average size of a DGEN, across the entire population, stays within a given range (70%–80% of L × Γ genes).

Model-4: With gene deletions. Model-4 adds DL events in Model-3 (complete model). The deletion probability P_DL is set so that the average size of a DGEN, across the entire population, stays within the same range as in Model-3.

In Model-1 and Model-2, genes can be only removed (due to RW events, potentially followed by RF cascades) and so the average DGEN size decreases as evolutionary time progresses, which is unrealistic. Model-3 and Model-4 are more realistic because they can maintain a roughly constant DGEN size in the long-term. However, as will be shown next, all aspects of the developmental hourglass effect can already be seen with Model-2 (but not with Model-1). This highlights that the increasing specificity assumption is sufficient to generate the hourglass effect. Further, the inclusion of additional biological mechanisms in the model, namely gene duplication and gene deletion, even though they make the model more realistic, they are not necessary for the emergence of the hourglass effect.

Hourglass shape. A first observation is that as evolutionary time progresses, DGENs acquire an “hourglass-like” shape in Models-2,3,4. This means that the width of each stage first decreases until a certain stage (referred to as the waist of the hourglass) and then gradually increases. The hourglass may not be symmetric with respect to the waist. To quantify this observation, we define an “hourglass score” H (see Methods and Figure 3) that is equal to 1 if the sequence of L stage widths consists of two segments: a decreasing sub-sequence of k ≥ 2 stages followed by an increasing sub-sequence of L − k ≥ 2 stages. Figure 6a shows the hourglass score for the population of DGENs in Model-2. Similar graphs for the three other models are shown in Figure 5a, Figure 7a, and Figure 8a. The H score quickly increases in the three models that have increasing specificity, and it fluctuates close to 1 afterwards.

Figure 5. Computational results with Model-1.

Parameters: 10 runs with different initial populations, N = 10 individuals, L = 10 stages, specificity function s(l) = 0.5 for all stages, Γ = 100 genes at each stage initially, RF parameter z = 4, 1,000,000 generations, probability of RW event P_RW = 10⁻⁴. The red line is the median and the green boxes are the 10th, 25th, 75th, and 90th percentiles, across all individuals and all runs. (a) The hourglass score H across evolutionary time. (b) Lethality probability at each stage. (c) Age of existing genes at the last generation.

Figure 6. Computational results with Model-2.

Parameters: 10 runs with different initial populations, N = 1000 individuals, L = 10 stages, specificity function s(l) = l/L (l = 1 . . . L − 1), Γ = 100 genes at each stage initially, RF parameter z = 4, 500,000 generations, probability of RW event P_RW = 10⁻⁴. The red line is the median and the green boxes are the 10th, 25th, 75th, and 90th percentiles, across all individuals and all runs. (A) The hourglass score H across evolutionary time. (B) Lethality probability at each stage. (C) Age of existing genes at the last generation.

Figure 7. Computational results with Model-3.

Parameters: 10 runs with different initial populations, N = 10 individuals, L = 10 stages, specificity function s(l) = l/L (l = 1 . . . L − 1), Γ = 100 genes at each stage initially, RF parameter z = 4, 1,000,000 generations, probability of RW event P_RW = 10⁻⁴. The probability of gene duplication P_DP is adjusted dynamically so that the average DGEN size stays between 700 and 800 genes. The red line is the median and the green boxes are the 10th, 25th, 75th, and 90th percentiles, across all individuals and all runs. The hourglass score H across evolutionary time. (b) Lethality probability at each stage. (c) Age of existing genes at the last generation.

Figure 8. Computational results with Model-4.

Parameters: 10 runs with different initial populations, N = 10 individuals, L = 10 stages, specificity function s(l) = l/L (l = 1 . . . L − 1), Γ = 100 genes at each stage initially, RF parameter z = 4, 1,000,000 generations, probability of RW event P_RW = 10⁻⁴. The probability of gene duplication P_DP is adjusted dynamically so that the average DGEN size stays between 700 and 800 genes. The probability of gene deletion (DL) is P_DL = 10⁻⁶. The red line is the median and the green boxes are the 10th, 25th, 75th, and 90th percentiles, across all individuals and all runs. (a) The hourglass score H across evolutionary time. (b) Lethality probability at each stage. (c) Age of existing genes at the last generation.

What is the reason behind the hourglass shape of DGENs? When a gene g is rewired at stage l, it may trigger RF events in stage l + 1 depending on the number of its lost outgoing edges. In the first few stages, where specificity is low, it is unlikely that a gene would lose a large fraction of its (typically many) incoming edges. In the last few stages, where specificity is high, edges are unlikely to get rewired in the first place. In the mid-stages however, where the specificity is close to 50%, there is higher variability in the number of outgoing edges lost or gained due to RW events. The loss of several outgoing edges due to an RW event at stage l can trigger several RF events and gene removals in the subsequent stage. Thus, the probability of RF events in mid-stages is higher than in early/late stages, making the removal of genes more likely in the former.

The constant specificity of Model-1 does not result in an hourglass pattern [Dataset 1] (see Figure 5a) for the following reason. RW events at stage l can cause RF events at the next stage with the same probability, independent of l. However, after the occurrence of an RF event, the size of the potential cascade increases as l decreases simply because there are more subsequent stages to affect. This gives DGENs a “funnel-like” shape with a gradually increasing number of genes after stage-1; H fluctuates around 0.5, as expected for an increasing sequence.

Stage lethality. Another aspect of the developmental hourglass is in terms of the significance of each stage for the survival of the embryo. We define lethality of stage l as the probability that a RW or DL event at stage l starts a RF cascade that eventually leads to a DF event. We estimate this probability at generation i as the fraction of RW and DL events, during the past i generations, that occurred at stage l and led to a DF.

In Model-1, there is no clear trend for the stage lethality probability (see Figure 5b); with the exception of the last stage (in which RW events cannot result in gene loss), the lethality probability is roughly the same at all stages. In the three models with increasing specificity, however, we observe a clear pattern: the lethality gradually increases until the waist of the hourglass, and then it decreases [Dataset 2, Dataset 3 and Dataset 4] (see Figure 6b, Figure 7b, and Figure 8b). The reason, as explained earlier, is that the probability of RF events in mid-stages is higher that in early/late stages. Additionally, after the formation of the hourglass shape the mid-stages have relatively few genes and so an RF event in those genes is more likely to initiate a lethal RF cascade.

These computational results for the lethality probability across development are consistent with the empirical observations of Galis and Metz¹⁵ about the increased mortality of rodents due to perturbations in the phylotypic stage (see also Discussion section).

Age of genes. A third aspect of the developmental hourglass effect is related to the evolutionary age of genes. The age of a gene g at generation i is defined as

          A(g) = i − t₀(g),

where t₀(g) is the generation at which g was most recently rewired (and 0, if it was not rewired earlier). The rationale behind this definition is that a rewiring event may give that gene a new function, at least in terms of its upstream and downstream regulators^32,33.

In the case of Model-2, Figure 6c shows the median age of the genes at each stage, considering the population of all genes across all individuals at a given generation. See Figure 5c, Figure 7c, and Figure 8c for the same results with the three other models. The evolutionary age at stage l follows the same pattern as the lethality probability: it gradually increases until we reach the waist of the hourglass, and then it gradually decreases. Genes at intermediate stages tend to be older because, as explained earlier, they are fewer and their rewiring is more likely to be lethal. When one of those genes g is rewired or deleted from a DGEN, the corresponding individual is often replaced (DF event) by another individual that has the same gene g. So, the genes at the waist of a DGEN tend to be more conserved than genes at earlier or later stages.

Location of waist. What controls the location of the hourglass waist in the developmental process?

The location of the waist is mostly affected by two parameters of the model: the shape of the specificity function and the RF parameter z. To examine the effect of the former we use the nonlinear function shown in Figure 9. γ is the stage at which the specificity is 50%, and so that stage has the maximum variance in the number of outgoing regulatory edges. RW events at this stage can cause the largest extent of rewiring and so, the highest likelihood of RFs in genes of the next stage. Figure 10a shows that the location of the hourglass waist is ”pushed” towards stage γ, even though it does not coincide always with that stage. Parameter z controls the shape of the RF probability: increasing z makes RF events more likely, also increasing the likelihood of lethal RF cascades. Figure 10b shows that as we increase z the hourglass waist moves towards later developmental stages [Dataset 5].

Figure 9. A nonlinear specificity function, s(l) = 0.9 − 0.8/(1 + e^γ−l), for three values of the parameter γ.

This function allows us to control the stage γ at which the specificity is 50%.

Figure 10. We examine the effect of the two model parameters that affect the location of the DGEN hourglass waist.

The first is the specificity function. To examine its effect, we use a sigmoid-like mathematical function that controls the stage γ at which the specificity is 50% (see Figure 9). This is the stage with the maximum variance in the number of outgoing regulatory edges. RW events at this stage can cause the largest extent of rewiring and so, the highest likelihood of RFs in genes of the next stage. Graph (a) shows that the location of the hourglass waist is “pushed” towards stage γ, even though it is not always exactly at that stage. The second way to affect the location of the hourglass waist is the parameter z that controls the shape of the RF probability. Increasing z makes RF events more likely, also increasing the likelihood of lethal RF cascades. Graph (b) shows that as we increase z the hourglass waist moves towards later developmental stages. These results are obtained using Model-2. Parameters: 10 runs with different initial populations, N = 1000 individuals, L = 10 stages, specificity function s(l) = l/L (l = 1 . . . L − 1), Γ = 100 genes at each stage initially, RF parameter z = 4, 500,000 generations, probability of RW event P_RW = 10⁻⁴. The graphs show the median (red lines) and the 10th, 25th, 75th, and 90^th percentiles (green boxes) of the location of the waist.

Gene prevalence. We introduce a “gene prevalence” metric for gene g at time t as the fraction of individuals that include g at evolutionary time t [Dataset 6]. Figure 11a shows the gene prevalence metric across developmental stages after 500,000 generations, while Figure 11b shows the relation between gene prevalence and gene age. The genes at the waist of the developmental hourglass are not only the oldest but also the most prevalent across the population. The implication of this simulation result is that we can expect that those genes that are transitioning near the waist of the developmental hourglass will be the most conserved genes in a population.

Figure 11. The prevalence of a gene g in a population of N individuals is the fraction of individuals in which gene g appears.

These results are obtained using Model-2. Parameters: 10 runs with different initial populations, N = 1000 individuals, L = 10 stages, specificity function s(l) = l/L (l = 1 . . . L −1), Γ = 100 genes at each stage initially, RF parameter z = 4, 500,000 generations, probability of RW event P_RW = 10⁻⁴. The graphs show the median (red lines) and the 10th, 25th, 75th, and 90th percentiles (green boxes) for: (a) prevalence of genes in each stage after 500,000 generations, and (b) gene age as a function of gene prevalence. As expected, older genes tend to be more prevalent in the population.

Data analysis

We have examined the predictions of the previous model using transcriptome data for Drosophila melanogaster and Arabidopsis thaliana [Dataset 14]. We summarize results from both species here; the corresponding figures for which are Figure 12 and Figure 13 [Dataset 8 and Dataset 9]. For Drosophila, we analyze Microarray⁹ and RNA-Seq²⁷ temporal expression profiles during the first 20 hours of development, taken at 10 stages of 2-hr intervals. We examine whether a) the number of transitioning genes follows an hourglass pattern, b) the waist of that hourglass coincides with the Drosophila phylotypic stage, and c) the evolutionary age of the transitioning genes follows a similar hourglass pattern. The two datasets are described in more detail in the Methods section. With such limited data, we cannot infer the regulatory edges between transitioning genes and we cannot reconstruct the underlying DGEN. However, we can identify the transitioning genes at each developmental stage given a “transition threshold” c (see Methods). Even though the correct value of this threshold is not known, the following results are robust in a wide interval of c, which includes most of the expression variation range across successive developmental stages (see Figure S3 for the CDFs of expression level variations across successive stages [Dataset 12]).

Figure 12. Drosophila results using normalized expression levels.

Graphs (A) and (B) show the hourglass score (normal and robust) as a function of the transition threshold c for the two datasets. Graphs (C) and (D) show the location of the hourglass waist (stage-pair) as a function of the transition threshold c for the two datasets. Graphs (E) and (F) show the Transcriptome Age Index of transitioning genes for three different values of c (chosen so that the number of genes with known age index assigned to each stage is at least 25) for the two datasets. Graph (G) shows the transitioning genes for the Microarray dataset with c = 0.0005. The transitioning genes constitute 11% of all genes in that dataset. 53% of those genes transition in a single stage-pair. Of the remaining, 64% transition only in consecutive stage-pairs. Note that if a gene transitions n times, it is counted in n stage-pairs. Similarly, graph (H) shows the transitioning genes for the RNA-Seq dataset with c = 0.00025. The transitioning genes constitute 5% of all genes in that dataset. 45% of those genes transition in a single stage-pair. Of the remaining, 52% transition only in consecutive stage-pairs.

Figure 13. Arabidopsis thaliana results using normalized expression levels.

Graph (a) shows the hourglass score (normal and robust) as a function of the transition threshold c. Graph (b) shows the location of the hourglass waist (stage-pair) as a function of the transition threshold c. Graph (c) shows the Transcriptome Age Index of transitioning genes for five different values of c (chosen so that the number of genes with known age index assigned to each stage-pair is at least 290). Graph (d) shows the CDFs of the expression level absolute variations |δ| across successive stage-pairs. Note that in this case the hourglass waist (in terms of number of transitioning genes) appears in stage-pair (3,4), while the oldest genes appear in the next stage-pair. Graph (e) shows the transitioning genes with c = 0.0001. The transitioning genes constitute 7% of all genes in that dataset. 49% of those genes transition in a single stage-pair. Of the remaining, 36% transition only in consecutive stage-pairs.

Figure 12a and Figure 12b show the hourglass resemblance score H (and its more robust variant) as function of c. Note that the H score is close to 1 for a wide range of c, confirming the presence of an hourglass-like structure in terms of the number of transitioning genes. Figure 12g and Figure 12h exhibit this pattern more clearly in the number of transitioning genes for a specific value of c. The two datasets also show reasonable agreement in terms of the assignment of transitioning genes to stage-pairs [Dataset 10] (see Figure S1).

Second, the location of the waist in this hourglass pattern, shown in Figure 12c and Figure 12d, occurs at the stage-pair (3,4) or (4,5), depending on c. This is roughly 8 hours after the formation of the zygote, and it includes the phylotypic stage for Drosophila melanogaster⁹.

We have also estimated the evolutionary age of most of the transitioning genes at each developmental stage-pair using the Transcriptome Age Index (TAI) metric¹¹ (see Methods). TAI is lower for older genes. Figure 12e and Figure 12f show the average TAI for transitioning genes, weighted by the expression level of each gene, at each stage-pair and for each dataset using three values of c. The TAI index follows the pattern that the model predicts, with older genes (lower TAI values) close to the waist of the hourglass. This result appears consistent with the main observation of Domazet-Loso and Tautz¹³, even though that study did not analyze transitioning genes.

The same approach was extended using transcriptome data²⁸ and TAI profiles¹¹ from Arabidopsis. The results obtained from the analysis (Figure 13) were consistent with the predictions of the model, as well as the results obtained from Drosophila data.