An explanatory evo-devo model for the developmental hourglass

The"developmental hourglass"describes a pattern of increasing morphological divergence towards earlier and later embryonic development, separated by a period of significant conservation across distant species (the"phylotypic stage"). Recent studies have found evidence in support of the hourglass effect at the genomic level. For instance, the phylotypic stage expresses the oldest and most conserved transcriptomes. However, the regulatory mechanism that causes the hourglass pattern remains an open question. Here, we use an evolutionary model of regulatory gene interactions during development to identify the conditions under which the hourglass effect can emerge in a general setting. The model focuses on the hierarchical gene regulatory network that controls the developmental process, and on the evolution of a population under random perturbations in the structure of that network. The model predicts, under fairly general assumptions, the emergence of an hourglass pattern in the structure of a temporal representation of the underlying gene regulatory network. The evolutionary age of the corresponding genes also follows an hourglass pattern, with the oldest genes concentrated at the hourglass waist. The key behind the hourglass effect is that developmental regulators should have an increasingly specific function as development progresses. Analysis of developmental gene expression profiles from Drosophila melanogaster and Arabidopsis thaliana provide consistent results with our theoretical predictions.


Introduction
The evolutionary mechanism of conservation during embryogenesis, and its connection to the gene regulatory networks that control development, are fundamental questions in systems biology [1][2][3] . Several models have been presented in the context of morphological, molecular, and genetic developmental patterns. The most widely discussed model is the "developmental hourglass", which places the strongest conservation across species in the "phylotypic stage". The first observations supporting the hourglass model go back to von Baer when he noticed that there exists a mid-developmental stage in which embryos of different animals look similar 4 . On the other hand, the "developmental funnel" model of conservation predicts increasing diversification as development progresses 5,6 .
Recently, the hourglass model has come under new light. Multiple studies have observed the hourglass pattern across diverse biological processes, including transcriptome divergence 7-12 , transcriptome age 7,13,14 , molecular interaction 15 , and evolutionary selective constraints 10,15,16 . Despite these observations the genomic basis and even the existence of the developmental hourglass effect have been the subject of an intense debate 1,6,14,[17][18][19][20][21][22] . More importantly, the underlying mechanism that can shape the developmental process in the hourglass or funnel forms is still unknown.
We aim to understand the conditions under which the hourglass effect can emerge in a general setting, based on an abstract model for the evolution of embryonic development. The model focuses on a hierarchical network that represents the temporal "execution" of the underlying Gene Regulatory Network (GRN) during development. Each layer of the network corresponds to a developmental stage. The nodes at each layer represent regulatory genes (i.e., genes encoding transcription factors or signaling molecules) that undergo significant activity change at that corresponding stage. The edges from genes at one layer to genes at the next layer represent regulatory interactions that cause those activity changes. We refer to this hierarchical network as Developmental Gene Execution Network (DGEN) to distinguish it from the corresponding GRN. A DGEN is subject to evolutionary perturbations (e.g., gene deletions, rewiring, duplication) that may be lethal, or that may impede development, for the corresponding organism.
The model predicts that the evolutionary process shapes the DGENs of a population in the form of an hourglass, under fairly general assumptions. Specifically, the number of genes at each developmental stage follows an hourglass pattern, with the smallest number at the "waist" of the hourglass. The main condition for the appearance of the hourglass pattern is that the DGEN should gradually get sparser as development progresses, with general-purpose regulatory genes at the earlier developmental stages and highly specialized regulatory genes at the later stages. This assumption is motivated from the well established patterns of increasing modularity as the embryo develops 23 . Under the previous condition, the model predicts that gene regulatory changes or rewiring in mid-development are more likely to cause cascades of removing non-essential genes from the DGEN, compared to early or late developmental stages. Another model prediction is that the evolutionary age of DGEN genes also follows an hourglass pattern, with the oldest genes concentrated at the waist.
We have examined the aforementioned predictions using transcriptome data from the development of Drosophila melanogaster and Arabidopsis thaliana. This data is insufficient to reconstruct the complete DGEN of these species but it allows to estimate the number of genes at each developmental stage, given an activity variation threshold. Under a wide range of this threshold, the inferred DGEN shape follows an hourglass pattern, the waist of that hourglass roughly coincides with the previously reported phylotypic stage for these species, and the age of the corresponding genes follows the predicted hourglass pattern.

Developmental gene execution networks
As a first-order approximation, a regulatory gene can be modeled in one of several discrete functional states 24 . 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

Amendments from Version 1
The main changes in Version-2 of this paper are: -New text in the Introduction and Model Description sections, providing connections between modularity (see review by G.P.Wagner et al.) and our assumption of increasing specificity as development progresses.
-Included additional context and insight in the Introduction about the results shown later in the paper.
-Clarified several points in the Model Description, especially about how we model the specificity function and the probability of regulatory failure.
-A new paragraph in the Discussion section, focusing on Raff's hypothesis in the context of this work.
-Provided a connection between this work (the computational results about lethality probability) and the work of Galis and Metz about perturbations in the phylotypic stage of rodents.
-Added several new references.
-Re-formatted all equations that were previously embedded in the main text ("inline") as separate lines.
-The quantitative results of the paper (simulations and data analysis) have not been modified in any way.

See referee reports
REVISED 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 25 , 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 23 .
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 Nexisting 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: 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: 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).
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.

Methods
Hourglass score H. Suppose that w(l) denotes the number of transitioning genes in stage l, also referred to as the width of stage l. Let b be the stage with the minimum number of such genes. We construct the sequences X = {w(l), l = 1, . . . b} and Y = {w(l), l = b, . . . L}. τ X and τ Y denote the normalized univariate Mann-Kendall statistic for monotonic trend in each sequence, respectively (τ is -1 for decreasing, +1 for increasing and almost 0 for random sequences) 26 .
The H score is defined as See Figure 3 for an illustration, and for the definition of a more robust version of H.

Statistical analyses.
All biological data processing and analyses were performed using custom scripts written in Java (JDK v1.6) [ Dataset 14].
The identification of transitioning genes follows the same method for both absolute and normalized expression levels. In the case of absolute expressions we define δ i,l = e i,l − e i,l-1 for each gene i and at each stage l = 2 . . . L, while in the case of normalized expressions, we similarly define δ i,l = e′ i,l − e′ i,l-1 . Gene i is considered "transitioning" at the stage-pair (l − 1, l) if where c is a given transition threshold. This condition is more robust to noise than a ratio-based rule (e′ i,l /e′ i,l-1 ) for the identification of transitioning genes. Note that a gene may be transitioning at more than one stage-pair, but it may also not be transitioning at any stage-pair.

Transcriptome age index (TAI).
We collected the groups of orthologs for each gene in Drosophila using two databases, OrthoDB 29 and OrthoMCL 30 . The Eumetazoa data were taken from OrthoDB, while Fungi and Plants species were retrieved from OrthoMCL, and the two datasets were merged. Using these orthologs we then assigned an age index to each gene based on its absence and presence in a phylogenetic tree of 24 well-diverged species (see Figure 4) Figure 3. Illustration of the H score calculation. Let w(l) be the width of stage l. Let w b be the minimum width across all stages, and suppose that this minimum occurs at stage l = b; this is the waist of the network (ties are broken so that the waist is closer to ⎣L/2⎦). Consider the sequence X = {w(l ) }, l = 1, . . . b} and the sequence Y = {w(l)}, l = b, . . . L}. We denote the normalized univariate Mann-Kendall statistic for monotonic trend on the sequences X and Y as τ X and τ Y , respectively. The Mann-Kendall statistic varies between -1 (decreasing) and 1 (increasing), and it is approximately zero for a random sequence. We define H = (τ Y − τ X ) / 2 ; H is referred to as the hourglass score. H = 1 if the DGEN is structured as an hourglass, with a decreasing sequence of b stages followed by an increasing sequence of L − b stages. In the computational modeling results, we do not consider the width of the first stage because it can never decrease in Models-1,2,3. We also define a variation of the hourglass score in which we do not take into account adjacent stages in calculating the two Mann-Kendall statistics. That statistic is denoted by H and is referred to as the robust hourglass score. genes are assigned to stage-pair (l − 1, l). Denote by p i the phylogenetic rank (TAI value) of gene i. The age index assigned to that stage-pair is given by The same method is used when transitioning genes are identified based on absolute expression levels.

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.
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 15 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 where t 0 (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].
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.

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; 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 −4 . The graphs show the median (red lines) and the 10th, 25th, 75th, and 90 th percentiles (green boxes) of the location of the waist. , Γ = 100 genes at each stage initially, RF parameter z = 4, 500,000 generations, probability of RW event P RW = 10 −4 . 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. the corresponding figures for which are Figure 12 and Figure 13 [Dataset 8 and Dataset 9]. For Drosophila, we analyze Microarray 9 and RNA-Seq 27 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 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).

Figure 12a and
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 9 .
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 11 (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 13 , even though that study did not analyze transitioning genes.
The same approach was extended using transcriptome data 28 and TAI profiles 11 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.

Discussion
Early studies of the developmental hourglass effect mostly analyzed morphological and phenotypic similarities across species 2,34 . Recently, the focus has shifted towards genomic and molecular comparative studies 8,9,11,13,20 that investigate conservation of gene expression variation, sequence conservation, selective constraint on coding sequences, and evolutionary gene "age". These studies often report contradicting observations: some support strong conservation in earlier developmental stages 5,6,20 , while others support that strongest conservation occurs at a mid-developmental stage [7][8][9][10][11][13][14][15][16] . Nevertheless, the fact that the hourglass effect is observed in highly divergent species across deep phylogenetic scales (including fish, flies and plants), suggests that this observed pattern of conservation may stem from fundamental organization principles.
What these principles are has remained elusive. Earlier stages may be conserved because any changes therein could have large cascading effects in later stages 5,35,36 . Later stages may experience less constraint because as development progresses gene interactions become more modular, and so it is plausible that perturbations there have only local effects 1 . We refer to them as the "temporal" constraint model and the "spatial" constraint model, respectively, following Tian et al. 37 .
In this paper, we developed an evolutionary model of development that combines some aspects of the previous two models. Regulatory perturbations at a certain stage can cause cascades of regulatory failures at subsequent stages (temporal model), while the likelihood that a gene regulates genes at a subsequent stage decreases as development progresses (spatial model).
Our computational results lead to the following testable predictions: a) the number of transitioning genes during development follows an hourglass pattern, b) the evolutionary age of the transitioning genes also follows an hourglass pattern, with the oldest genes being at the waist of the hourglass, and c) the genes at the waist of that hourglass are the most essential, in the sense that their deletion maximizes the probability of developmental failure.
We have relied on developmental gene expression profiles of Drosophila melanogaster and Arabidopsis thaliana to examine the predictions of the model. The analysis of that data agrees with the  15 . That study has shown, based on the teratological literature for rodent development, that disturbances in the phylotypic stage lead to much higher mortality than in other stages. Further, such disturbances lead to the co-occurrence of several distinct anomalies in the developing embryo and so the increased mortality cannot be due to a single particularly vulnerable process that takes place at that stage. Further, our simulations confirm that the details of these regulatory perturbations, such as the probabilities of gene duplication and deletion or the parameter z in the regulatory failure probability, do not affect the results of the model, at least at the qualitative level.
It is interesting to examine our results in the context of Raff's hypothesis 1 . Raff proposed an interesting explanation for the developmental hourglass effect based on the notions of "interconnectivity between body elements" and "developmental flexibility". Specifically, his hypothesis states that the level of interaction (interconnectivity) between the elements of the developing embryo is maximized in mid-development. Early development is flexible because it governs robust and general global patterning processes, late development is also flexible because "signaling events within the primordia are little influenced by events in other primordia", while mid-development (phylotypic stage) is least flexible because of the "high interconnectivity between elements that will later come to represent separate modules". This hypothesis may be viewed initially as a contradiction with our increasing specificity assumption, which states that the density of regulatory interactions is maximized at the earliest stages of development. Note however that Raff's hypothesis was not stated in terms of gene regulatory interactions -he was referring more generally to "developmental flexibility". If we interpret the notion of "developmental flexibility" as the ability of an embryo to survive gene mutations and rewiring at different stages of the developmental process, then Raff's hypothesis is actually consistent with our results regarding the "lethality probability" at each developmental stage (see Figure 6-B, Figure 7-B and Figure 8-B). The lethality probability is maximized at the phylotypic stage, and it is significantly lower at early and late developmental stages, following exactly the same pattern with Raff's developmental flexibility.
The use of DGENs in this work was only as an abstract tool to study the effect of gene regulatory perturbations in the developmental process. In future work, it is important to infer the actual DGEN of model organisms. This will require information about gene regulatory interactions across time and space, but it should be possible for at least some developmentally well studied species 24 . Such DGENs would help to identify the specific genes that form the hourglass waist and their function. Additionally, an inferred DGEN would allow to directly test the increasing specificity assumption.
Finally, we note that the hourglass effect (sometimes referred to as the "bow-tie" effect) has also been observed in other complex biological and technological systems that exhibit hierarchical modularity and that are subject to evolutionary pressure or optimization tradeoffs 38-41 . Recently, Friedlander et al. have proposed a different mechanism for the emergence of hourglass-like patterns in evolving biological or technological networks, based on a linear system formulation 42 . They showed that if a system can be represented as a hierarchical linear transformation of an input vector to an output vector, and the desired transformation matrix is rank-deficient, then an evolutionary process that selects that particular transformation can, under certain conditions, converge to an hourglass-like structure. It is not clear yet how to adapt this linear model in the context of inherently non-linear systems, such as gene regulatory networks.
Another example of a hierarchical system that is structured as an hourglass is the Internet "protocol stack" 43 ; this pattern was not designed but it emerged through the competition between protocols that serve roughly the same function at each communication layer, during the last 30-40 years. In earlier work, we proposed an abstract model (EvoArch) that captures the evolution of protocol architectures and that predicts the emergence of an hourglass structure. Interestingly, both EvoArch and the model of this paper share the same principle: the underlying hierarchical networks that control both systems should be increasingly sparser as complexity increases, i.e., the specificity of each complexity stage (or layer) should be increasing. In the future, we will further investigate this common organization principle between biological and technological systems.

Competing interests
No competing interests were disclosed.  Figure S1. We evaluated the agreement between the two Drosophila datasets in terms of the transitioning genes assigned to each stage-pair, considering only those genes that appear in both datasets. Because the appropriate transition threshold may be different at each dataset, we use a different threshold for each dataset, say c 1 and c 2 . For each pair (c 1 , c 2 ), we determine the transitioning genes at each stage-pair with the corresponding dataset., L − 1 pairs of gene sets), and then calculate the average Jaccard similarity across these L − 1 pairs. The Jaccard similarity maps show this average across all stage-pairs for various threshold pairs (c 1 , c 2 ). In graph (a) with normalized expression levels, when the two thresholds are roughly equal, the average Jaccard similarity is as high as 50%; this means that about 2/3 of the genes assigned to a certain stage-pair using one dataset are also assigned to the same stage-pair using the other dataset. Graph (b) shows a similar Jaccard similarity map for the case of absolute expression levels. We highly appreciate that the authors carefully addressed all our concerns and we now have no further reservations.

Supporting information
We also add a detailed response to the rebuttal letter:

Point#1 of our summarized comments to the first version of the manuscript: "The authors do cite Raff in the first sentence of their introduction, but should state more clearly in which manner their model was inspired by Raff's hypothesis."
Authors' response: The reviewers are right to highlight the connection between our "increasing ecificity" assumption and Raff's hypothesis. We have included a paragraph in the regulatory sp Discussion section to put Raff's hypothesis in the context of this work.
A first remark is that our assumption is consistent with Raff's basic premise that modularity increases as the embryo develops. Second, Raff's hypothesis (see chapter-6 of [1], and in particular Figures 6.6 and 6.7) also states that the "level of interaction" or the "interconnectivity between (body) elements" is maximized at the phylotypic stage. This may be viewed initially as a contradiction between Raff's hypothesis and our increasing specificity assumption, which states that the density of regulatory interactions is maximized at the earliest stages of development. Note however that Raff's hypothesis was not stated in terms of gene regulatory interactions --he was referring more generally to "developmental flexibility", arguing that early development is flexible because it governs robust and general global patterning processes, late development is also flexible because "signaling events within the primordia are little influenced by events in other primordia", while mid-development (phylotypic stage) is least flexible because of the "high interconnectivity between elements that will later come to represent separate modules." If we think of "developmental flexibility" as the ability of an embryo to survive gene mutations and rewiring at different stages of the developmental process, Raff's hypothesis is actually consistent with our results regarding the "lethality probability" at each developmental stage (see Figs 6-B, 7-B amd  8-B). The lethality probability is maximized at the phylotypic stage, and it is significantly lower at early and late developmental stages, following the same pattern with Raff's "developmental flexibility".
We appreciate that the authors now differentiate between Our comment to authors' response: We appreciate that the authors now differentiate between Our comment to authors' response: Raff's hypothesis concerning the modularization (interconnectivity) between body elements (modules) and their "increasing regulatory specificity" assumption. We agree that Raff's hypothesis is consistent with the "lethality probability" proposed by the authors in terms of the observed pattern of maximal interconnectivity and maximal probability of lethal effects during the phylotypic stage (mid embryogenesis). However, we would like to point out that a causal relationship between the observed patterns remains to be tested.
Point#2 of our summarized comments to the first version of the manuscript: "Correlate and reference major assumptions of the DGEN model that seem plausible with findings from experimental studies that validate this plausibility."

First, to the extent of our knowledge there is no prior work that Authors' response to Point #2:
validates the assumption of increasing regulatory specificity as formulated in our paper directly (i.e., with the specificity of a developmental stage defined as the topological density of regulatory edges in that stage). We are planning to examine the validity of this assumption in the future using various approaches: First we plan to analyze the aforementioned developmental GRNs for the sea-urchin that have been directly examined by Dr. Davidson's lab at Caltech.

Second, we are in the process of examining other genomic metrics, which can be viewed as "proxies" to regulatory specificity. For example, we have defined new measures to quantify how specific a gene's expression profile in a developmental stage is compared to other developmental stages, which we have tentatively named as "developmental stage specificity index (DSI)".
Similarly, we defined "tissue specificity index (TSI)", measuring the bias of a gene's expression in a specific tissue compared to many other tissues. Analyses of these measures using genomic data from Drosophila so far have revealed generally increasing trends, which are consistent with the idea of increasing specificity. However we are still working on generalizing these analytical tools to data sets generated by different methods and from different species. Importantly, our model focuses DGENs underlying development, and genomic data can obscure the true signals from genes constituting DGENs. Thus we are also working on examining these specificity measures using inferred DGENs. Additionally, we are developing other methods to approximate functional specificity. We hope to complete these analyses as a follow-up paper. The current paper provides theoretical motivation.

Third, as the reviewers also point out, there is a large body of prior work in developmental biology that confirms the increasing modularity in the developing embryo (for instance, see Wagner et al., ) as well as the increasing specificity of the developmental process at the signaling or 2007 genomic level. The connections between this pattern of increasing modularity and the structure of the underlying gene regulatory networks are still not well understood however. We think that it would be misleading if we had argued that these modularity patterns provide a direct validation for our "increasing regulatory specificity" assumption. Please note that we have cited the paper by Wagner et al. at the Introduction and Model Description sections, providing some connections between this work and earlier work on the role of modularity in developmental biology.
The authors point out that "the connections between this Our comment to authors' response: pattern of increasing modularity and the structure of the underlying gene regulatory network are still not well understood". This statement is valid for experimental studies that aim to address this presumptive causal relationship in a systematic manner and we appreciate that the authors share 3.

4.
presumptive causal relationship in a systematic manner and we appreciate that the authors share their future plans to introduce sufficient measures (DSI and TSI) to better understand the nature of this relationship.

Point#3 of our summarized comments to the first version of the manuscript: "Correlate the output of the DGEN simulation with existing findings of studies addressing similar questions of the causality of the developmental hourglass phenomenon and discuss potential simplifications of the model assumptions that do not fit experimental results or that are not derivable via experiments yet."
This is a very interesting idea and we thank the reviewers for this Authors' response to Point #3: suggestion. We are planning to pursue this investigation in a follow-up paper that will also include the previously mentioned analysis of sea-urchin GRNs as well as the DSI/TSI results. We believe that it would be distracting to try to include all these results in this first paper, given that the main focus of this work has been on the computational model and the computational results. We are also open to collaborate in this investigation with the reviewers or other researchers working in this area, if they are interested.
We are delighted to hear that the authors address this Our comment to authors' response: question in a follow up paper. We gladly offer our expertise to collaboratively investigate the evolutionary process on the regulatory level that might have led to a correlation between the maximized evolutionary age of expressed genes during mid-embryogenesis and the interconnectivity of (body) modules that are hypothesized to cause the morphological patterns during the phylotypic stage.

Points#4 and 5 of our summarized comments to the first version of the manuscript: a clear statement how the exact values of the stage specificity (regulatory specificity) s(l) have been derived. State more clearly that s(l) is the crucial parameter determining the observed phenomena."
The regulatory specificity functions we experimented Authors' response to Points #4 and 5: with are arguably the simplest one could think of. Specifically, we start (Model-1) with the simplest hypothesis that the regulatory specificity does not vary across developmental stages (in the paper we show results for s=0.5 but we have generated results for all values 0.1, 0.2, ... 0. 9 Figures-9 and 10 we consider a more general specificity function that increases in a non-linear manner, providing a simple way to control the stage at which it gives the mid-range value 0.5. This allows us to examine how the shape of the specificity function s(l) affects the location of the hourglass waist. We hope that follow-up work will reveal experimentally the actual shape of the regulatory specificity function s(l) and there is no doubt that it will not be identical to any of these simple mathematical functions.

and the conclusion remains that when the specificity is constant the network does not evolve into an hourglass shape). Then, we consider (Model-2) a linearly increasing specificity function. We do not argue that this is realistic or that there are experimental results that suggest this linearity or even the increasing trend. It is just the simplest and most parsimonious model that leads to the emergence of an hourglass pattern. Finally, in
Please note that we have revised the paragraphs about "Model-1" and "Model-2" in the Simulation section to clarify these points.
We thank the authors for clarifying the paragraphs about Our comment to authors' response: "Model-1" and "Model-2". 5.

"For the probabilities Specific reviewer comment on the first version of the manuscript P_{DL} and P_{DP} the explicit formula to compute the corresponding probabilities could also be included within the paragraph analogous to P_{RF}. This would enable a clearer reproducibility."
Actually there is no formula for the gene deletion and duplication Authors' response: probabilities. They are given probability values that do not depend on any other parameters. At each round of the simulation, a gene is duplicated with probability P_{DP}. If not rewired or duplicated, the gene is deleted with probability P_{DL}.
We thank the authors for clarifying this point for us. Our comment to authors' response: Figure 4, and mapped them to six different "phylostrata" (marked in Figure 4). Akin to the method described in the above references, we assigned each gene to a specific "phylostratum" based on the farthest detectable homologue of that gene, and this is referred to as "age index" in the paper. This made it possible to compare the age index based results we obtain in Drosophila to the phylostratigraphy based assignments in Arabidopsis.

"Explain why Point#7 of our summarized comments to the first version of the manuscript: merging OrthoDB results with OrthoMCL results might be plausible and why comparing it with results obtained from phylostratigraphy is reliable." The phylostratigraphy technique as employed by Domazet-Loso Authors' response to Point#7: et al. and Quint et al. involves the assignment of an evolutionary age to each gene in a given species, by tracing the most distant ancestral node containing at least one immediate daughter species with a detectable homologue. A phylogenetic tree of divergent species ranging from cellular organisms transitioning into simple and more complex eukaryota served as a reference to divide phyla radiating from consecutive ancestral points into distinct phylostratum layers. To determine the homologues, BLAST sequence similarity searches were done against complete genomes reliably annotated across the different phylostrata. Similarly, we obtained known homologues of each Drosophila gene from OrthoDB5 and OrthoMCL5 (that also use BLAST based searches) in the species depicted in
Homologues for Drosophila genes were sourced from these two databases, to increase coverage across the species depicted in the tree in Figure 4. Since these two databases use slightly different parameters to identify orthologs, we examined whether the usage of specific databases introduce bias by analyzing data from each database separately and observed that age index of each gene remained the same.
We are glad that the authors carefully addressed our initial Our comment to authors' response: concerns about the reliability of results obtained from merging methods of gene homology detection. Since the observed age indices are not influenced by the database used to assign gene homology relationships, we agree that merging OrthoDB5 and OrthoMCL5 results do not bias this specific analysis.

Reviewer summary
All remaining concerns were carefully revised and we thank the authors for investing the additional time to respond to our review in such detail.
No competing interests were disclosed.

Competing Interests:
We have read this submission. We believe that we have an appropriate level of expertise to confirm that it is of an acceptable scientific standard. The proposed model is based on a hierarchical network representation of the Gene Regulatory Networkthe so-called Developmental Gene Execution Network (DGEN) -and evolution proceeds through a series of stochastic perturbations involving gene duplication, re-wiring and deletion. The manuscript suggests that the key factor that reproduces the hourglass effect is the monotonicity in "time" (increasing) of a function they call specificity s(i) (a measure of likelihood of a gene regulating descendant genes at later times) that competes against perturbative effects such that there is a "waist" of the number of transitions genes w(l) at intermediate times. They then test predictions of the model in two datasets of developmental data and finding good qualitative agreement.

Drosophilia melanogaster Arabidopsis Thaliana,
The manuscript is properly motivated and well written (with some caveats, see comments below). The model proposed is quite intuitive and compelling and represents an excellent first attempt in beginning to uncover the mechanism behind this intriguing phenomenon. I anticipate that this will have high impact (considering the hourglass effect transcends biological phenomena and is also found in "designed" systems such as the Internet ) in multiple fields. Consequently I strongly endorse this protocol stack manuscript. 1.

2.
3. I do however have some concerns about some of the details about the model and therefore propose the following that may improve the readability of the manuscript and the plausibility and robustness of the model: The section describing the model is a bit dense, and can benefit from some rewriting. For example it is a bit hard to follow what the significance of a spatial domain is and its relevance to the model. In addition a bewildering plethora of parameters are introduced, that can be a bit overwhelming for the reader. Consequently I suggest putting all the important and relevant parameters in a table for easy access to the reader. The section on Rewiring (RW) is quite hard to read and will benefit from the introduction of explicit equations similar to that for P_{RF}.
Which brings me to my second point about the assumptions of the model. It seems to me that two key factors lead to the hourglass effect: a) the montonocity of s(l) with stage l and likewise b) of P_{RF} with r. While the choice of this dependence for the specificity is a simple linear dependence (which makes sense as a first pass), the choice of function for P_{RF} is highly non-trivial. Why not for example choose a simple linear dependence such as P_{RF} ~ r? The reason why I say this is that the authors would like to propose the most general and minimal model to explain this phenomenon. And it is not clear to me whether the specific choice of the functional dependences matter. For example with a non-linear choice of P_{RF} we see that a non-linear s(l) has the effect of shifting the waist. But what about a linear choice of P_{RF} and s(l)? Does that still preserve the hourglass effect? In my opinion a truly robust model will state that the main things that determine the hourglass are the monotonic dependencies in combination with non-linear or linear choices of the functions independent of their specific forms. This will allow the model to be applied to multiple settings.
At the moment I'm a bit concerned that the non-trivial choice of P_{RF}, as it stands in the manuscript, might be key to the observed effect. If it is indeed so, then one must motivate why one must make this particular choice. Same goes for the slightly peculiar choice for the non-linearity of s(l) introduced later in the model. The authors will do well to explain in detail what motivated them to make such a choice.
Additionally it would be helpful if the authors make it clearer (than they already have) that effects such as duplication and deletion are there only to make the model more realistic in terms of maintaining the number density of genes across different stages and play no part in the hourglass effect (or at least that's how I understand it).

Some minor nitpicks:
In the section Model Description, the authors refer to the Wright-Fisher model. I was not aware of what this was and had to look it up, so a reference for the reader should be provided.
The Hourglass score H is determined through Mann-Kendall statistics. A brief description of the method should be provided for the benefit of the reader (maybe in the supplementary material). Additionally, why this particular choice? Presumably any statistical test for monotonic trends should be robust to the parameters of the model. Would there be much of a difference is one used the Spearman rho for example?
Finally I think the authors will do well to highlight the importance of Specificity right at the 3. Finally I think the authors will do well to highlight the importance of Specificity right at the outset of the manuscript and move some simplified variant of the segment "What is the reason behind the hourglass shape of DGENS?" (page 7, line 11, second column) to somewhere in the introduction.
No competing interests were disclosed.

Competing Interests:
I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above.
Author Response 11 Dec 2014 , Georgia Institute of Technology, USA

Constantine Dovrolis
First, we would like to thank the author of this review for his very thorough reading of our paper. His comments and suggestions have been extremely helpful in revising the paper (please look at Version-2). We would also like to apologize for not revising the paper earlier.

Review: "The section describing the model is a bit dense, and can benefit from some rewriting. For example it is a bit hard to follow what the significance of a spatial domain is and its relevance to the model."
: We hope our extensive revision has mitigated this problem of the previous Response version of the paper.
We explain what we mean by "spatial domains" early in the paper (and in Fig-1) because it is an important and necessary concept when we justify later in the paper the assumption of increasing regulatory specificity (see paragraph that starts with "The major assumption is that the regulatory specificity…"). It is the formation of those distinct spatial domains (or " modules") that explains why the regulatory specificity is probably decreasing as development progresses.

Review: "The section on Rewiring (RW) is quite hard to read and will benefit from the introduction of explicit equations similar to that for P_{RF}."
: We have revised that paragraph to make it more clear. Response

Review: It seems to me that two key factors lead to the hourglass effect: a) the monotonicity of s(l) with stage l and likewise b) of P_{RF} with r. While the choice of this dependence for the specificity is a simple linear dependence (which makes sense as a first pass), the choice of function for P_{RF} is highly non-trivial. Why not for example choose a simple linear dependence such as P_{RF} ~ r? For example with a non-linear choice of P_{RF} we see that a non-linear s(l) has the effect of shifting the waist. But what about a linear choice of P_{RF} and s(l)? Does that still preserve the hourglass effect? In my opinion a truly robust model will state that the main things that determine the hourglass are the monotonic dependencies in combination with non-linear or linear choices of the functions independent of their specific forms."
: We chose the non-linear P_{RF} function shown in Fig-2 because it is quite Response flexible and it allows us to examine how the shape of this curve affects the behavior of the model by simply varying the parameter z. Please note that when z is equal to one this flexible and it allows us to examine how the shape of this curve affects the behavior of the model by simply varying the parameter z. Please note that when z is equal to one this function is almost linear. As discussed in Fig-10 however, as z decreases towards one the location of the hourglass waist also decreases. Consequently, we expect that a linearly increasing P_{RF} function would result in a shape that resembles a "funnel" rather than an hourglass.
We hope that new experimental work will provide some direct evidence for the shape of this function in the near future.

Review: "The Hourglass score H is determined through Mann-Kendall statistics. A brief description of the method should be provided for the benefit of the reader (maybe in the supplementary material). Additionally, why this particular choice? Presumably any statistical test for monotonic trends should be robust to the parameters of the model. Would there be much of a difference is one used the Spearman rho for example?"
: We have added a reference for this statistic: : The reviewer is right. We revised the fourth paragraph of the introduction to Response include a "hint" about the main reason behind the hourglass effect. Obviously, we cannot write much more at that early point of the paper because we still haven't defined what we mean by DGEN, RFs, cascades, etc.
No competing interests were disclosed. The article proposes and analyzes a network model that is capable of reproducing the hourglass effect in the number of genes that undergo a functional state transition as a function of the developmental stage. The hourglass effect refers to higher morphological divergence at earlier or later stages of embryonic developmental, compared to medium stages. The main idea of the model is to explain the waist of this hourglass (the number of "transitioning genes" w(l) having a minimum at medium developmental stages l) by an interplay between the impact of random rewiring events in the developmental gene execution network (DGEN) and stage specificity. The DGEN is a network representing regulatory interactions between genes at different developmental stages, while stage specificity s(l) defines the probability 1-s(l) with which a gene at stage l regulates a gene at stage l+1. One of the key assumptions in the model is that s(l) is an increasing function of l. As a result of the interplay between a decreasing impact of DGEN perturbations as a function of l, and increasing s(l), the lethality of perturbations is maximized at mid l's, leading to evolutionary incentives to minimize w(l) at those mid stages. The article consists of three parts. The first part describes the model in detail. The second part discusses extensive simulation results of the model. The third part presents an extensive analysis of existing developmental data on Drosophila and in the model context.

melanogaster Arabidopsis thaliana
The most interesting aspect of the article is that provides a possible intriguing explanation of the hourglass effect in developmental biology, which may foster future creative thinking and research in this important direction. The most obvious reservation one can express about the article is that the model cannot be currently refuted since it is formulated at the DGEN level, and there is currently no data from which a reliable DGEN reconstruction would be possible. The fact that some outcomes of model simulations qualitatively agree with available data does not directly validate either the model, or its main assumptions. The article does not contain such claims, however, that would not be supported by the data.
My minor comments that may help to improve the article are: The main idea of the model is buried somewhere on page 6 in the middle of a paragraph after sentence "What is the reason behind the hourglass shape of DGENs?" Why not explain this as early as possible, certainly in the Introduction, if not in the abstract?
The third paragraph in the Introduction reads somewhat unclear and imprecise, as well as the "Developmental gene execution networks" section.
When D(g) and Gamma notations first appear, it's not spelled out what they are.

DGEN removals upon DF events have a clear meaning. But what do replacements correspond to in reality?
The Mann-Kendall statistics could be briefly described, or at least a reference could be provided.
No competing interests were disclosed.

Competing Interests:
I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.
Author Response 11 Dec 2014 , Georgia Institute of Technology, USA

Constantine Dovrolis
First, we would like to thank the author of this review for his thorough reading of our work. His comments have been very helpful in revising this paper (please see Version 2). We would also like to apologize for not revising the paper earlier.
Review: "The most obvious reservation one can express about the article is that the model cannot be currently refuted since it is formulated at the DGEN level, and there is currently no data from which a reliable DGEN reconstruction would be possible. The fact that some outcomes of model simulations qualitatively agree with available data does not directly validate either the model, or its main assumptions. The article does not contain such claims, however, that would not be supported by the data." : It is true that a direct validation of the model is not possible with any currently Response available data (at least to the extent of our knowledge). Such a direct validation would require the reconstruction of the complete DGEN for specific model organisms, i.e., knowing not only the expression profile of each gene as a function of developmental time and at different tissues of the embryo but also knowing the regulatory inputs for each gene during development. However, there are research groups that are working in that direction (e.g., Prof. Davidson's group at Caltech is gradually "reverse-engineering" the DGEN of the sea-urchin). Our hope is that a direct validation of our model will be possible within the next few years.
Having said that however, even though our model cannot be directly validated, we are providing indirect evidence that certain predictions of the model are true in Drosophila and Arabidopsis (see Results -Data Analysis section). Consequently, the model should not be viewed as completely hypothetical either.
Review: "The main idea of the model is buried somewhere on page 6 in the middle of a paragraph after sentence "What is the reason behind the hourglass shape of DGENs?" Why not explain this as early as possible, certainly in the Introduction, if not in the abstract?" : The reviewer is right. We revised the fourth paragraph of the introduction to Response include a "hint" about the main reason behind the hourglass effect. Obviously, we cannot write much more at that early point of the paper because we still haven't defined what we mean by DGEN, RFs, cascades, etc. The revised paragraph is: "The model predicts that the evolutionary process shapes the DGENs of a population in the form of an hourglass, under fairly general assumptions. Specifically, the number of genes at each developmental stage follows an hourglass pattern, with the smallest number at the "waist" of the hourglass. The main condition for the appearance of the hourglass pattern is that the DGEN should gradually get sparser as development progresses, with general-purpose regulatory genes at the earlier developmental stages and highly specialized regulatory genes at the later stages. Under this assumption, the model predicts that gene regulatory changes or rewiring in mid-development are more likely to cause cascades of removing non-essential genes from the DGEN, compared to early or late developmental stages. Another model prediction is that the evolutionary age of DGEN genes also follows an hourglass pattern, with the oldest genes concentrated at the waist." Review: "When D(g) and Gamma notations first appear, it's not spelled out what they are."

Review: "When D(g) and Gamma notations first appear, it's not spelled out what they are."
: We checked this. The notation D(g) is introduced at the end of the Response "Developmental gene execution networks" section, and that is the first time it appears in the paper. The notation Gamma first appears in the "Developmental Failure (DF)" paragraph, and that is also where it is defined.

Review: "DGEN removals upon DF events have a clear meaning. But what do replacements correspond to in reality?"
: DGEN replacements model the effect of selection: an embryo that did not Response develop properly (DF) dies and its DGEN is removed from the population, while healthy embryos develop and eventually give birth to new a embryo with the same genotype as its parent (asexual reproduction).
Review: "The Mann-Kendall statistics could be briefly described, or at least a reference could be provided." For this purpose their evolutionary model focuses on hierarchical gene regulatory networks that control the corresponding developmental processes. Based on the output of the model the authors predict the emergence of an hourglass pattern in the structure of a temporal representation of the underlying gene regulatory network and correlate this effect with the evolution of protocol architecture found for the Internet's "protocol stack".
Based on this universal finding the authors speculate that the developmental hourglass pattern might be a causal effect of a common organization principle during the establishment of complexity and is based on the underlying hierarchical network structure.
the underlying hierarchical network structure.
The paper is very interesting, well written, and will potentially help the community to better understand the developmental hourglass phenomenon. However, we see some issues with the model assumptions and strongly suggest including the following aspects in a revised version. In general, most of our criticisms concern the model description and the methods section. Both are integral to the validity of the obtained results, which are nicely reported and discussed.

Introduction
Modeling the evolution of embryonic development as directed acyclic graph (DAG) in which the nodes correspond to state-transitioning genes and the edges model directed regulatory interactions causing significant activity change at the corresponding stage seems plausible. The authors refer to this model as DGEN. A testable hypothesis proposed by the authors is that if regulatory genes become increasingly function-specific as development progresses, the network gradually should become sparser in that direction and evolutionary older genes should be concentrated at the waist of the hourglass (Figure 1). This hypothesis is essential to detect possible regulatory mechanisms that cause the developmental hourglass pattern. We therefore suggest including references to existing biological studies investigating the phenomenon of morphological complexity being correlated with regulatory specialization to support the validity and plausibility of their DGEN model assumptions.
A biological motivation of this hypothesis comes from Rudolf A. Raff who proposed that early embryogenesis is governed by global pattern formation processes and that during later development, the embryo is being organized into developmental modules, allowing each module to differentiate autonomously. Raff furthermore proposed that the phylotypic stage marks the transition from early global development to the modular mode later in development . 1 development to the modular mode later in development .
The authors do cite Raff in the first sentence of their introduction, but should state more clearly in which manner their model was inspired by Raff's hypothesis.
Numerous studies investigated the phenomenon that increasing morphological complexity is correlated with increasing modularity and specialization of the underlying regulatory apparatus . Here the authors predict that the main condition leading to the appearance of the hourglass pattern is that the DGEN should gradually get sparser as development progresses, with general-purpose regulatory genes at the earlier developmental stages and highly specialized regulatory genes at the later stages, expressing the oldest genes during mid development.
We suggest to also include findings investigating the functionality and age of the genes that have been shown to be active during mid embryogenesis to correlate potential outcomes of the DGEN prediction that evolutionary older genes should be concentrated at the waist of the hourglass with results obtained by recent studies . These aforementioned studies were able to show that during mid embryogenesis evolutionary old genes or evolutionary old processes that are shared within and between phyla are most active.
In our opinion it would be very interesting to test the degree of intersection between genes predicted to be causal for the waist of the DGEN hourglass and evolutionary conserved genes that have been found to be expressed during mid embryogenesis. We are aware that the model uses the number of genes at each developmental stage combined with the age of the corresponding genes to observe constraints during development. Hence, only the authors can estimate how much effort would be needed to identify such intersecting genes.
If the search for intersecting genes returned by the DGEN model during mid development and genes found by previous studies should be possible, the authors could integrate a motivation section to find these intersecting genes within the Introduction section.

Model description
The model description is quite intuitive. However, as the validity of the whole study depends on the model parameters chosen, several points need clarification including the rationale and motivation to choose exactly these parameters.
One of the most crucial parameters of the model is the regulatory specificity , with 0 <= <= 1. We s(l) s(l) think, that adding a more detailed paragraph discussing how is being estimated, would be highly s(l) beneficial due to the importance for the parameter for the entire model and predictions. In Model-1-4 s(l) the authors choose three different values for : In Model-1: = 0.5 and in Model-2-4: = l/L, and in s(l) s(l) s(l) Figure 9 a non-linear function = 0.9 -(0.8 / 1 + exp(\gamma -l) ). A gene g at stage l is modeled to act s(l) as upstream regulator for a gene g' at stage l + 1 with probability = 1 -. From this follows that each s'(l) s(l) potential upstream regulator in stage l has the same probability to regulate genes g' in stage l + 1. In s'(l) other words, Model-1 assigns each potential upstream regulator with the same probability = 1 -0.5 = s'(l) 0.5 to regulate genes g' in stage l + 1, in Model-2: '(l) = 1 -l/L = const., and for the non-linear function s s(l) : = 1 -(0.9 -(0.8 / 1 + exp(\gamma -l) )) = const.

s'(l)
It would be advantageous, if the authors could include motivations for the three cut-offs, especially how they derived the sigmoid-like mathematical function ( Figure 10, legend caption) they later denote as non-linear function . We do see that modeling the stage specificity with a constant value, an increasing s(l) . We do see that modeling the stage specificity with a constant value, an increasing s(l) value, and a non-linear function is plausible, but a clear motivation or derivation of the exact values would be beneficial for a better understanding of the modeling process.
For the probabilities and the explicit formula to compute the corresponding probabilities could P P also be included within the paragraph analogous to P . This would enable a clearer reproducibility.
The authors discuss that their major assumption is that the regulatory specificity increases substantially as development progresses, hence the DGEN becomes gradually sparser along the developmental time axis starting with ~ 0, … , ~ 1. Here the authors should include references to experimental s (1) s(L) observations that support their assumption. Since Models 2-4 rely on this assumption, published experimental studies in line with this assumption need to be referenced. Otherwise, the models applied would seem largely academic.

Hourglass score H
Here is defined as the number of transitioning genes in stage l. In Figure 3 and in section w(l) 'Hourglass is referred to as . In case the definition of is equivalent to the , it shape' w(l) stage width w(l) stage width would be beneficial for the reader to see the term in the same sentence as the initial definition stage width of . A correct understanding of the is important to understand the measurement of . w(l) stage width H Furthermore, the univariate Mann-Kendall statistic should be referenced to allow a non-statistical community to better understand this way of statistical modeling.

Transcriptome age index (TAI)
In this section the authors describe and perform two different methods to assign each gene a corresponding evolutionary age. For they collected groups of orthologs from OrthoDB D. melanogaster and OrthoMCL, whereas for they obtained phylostratum assignments from Quint et al. that A. thaliana was derived by phylostratigraphy[ref- 13.]We do not understand the motivation of merging orthologous gene groups obtained from OrthoDB and OrthoMCL, as both approaches apply different parameters and heuristics to determine orthologous genes. Was this done to maximize the number of potential orthologs to be included in the study? In addition, why did the authors use a different approach for the plant gene set? This seems somewhat inconsistent and statement concerning the plausibility to compare results returned by the DGEN model based on age assignments obtained from phylostratigraphy ( ) A. thaliana with OrthoDB and OrthoMCL merged orthologs ( ) would be appreciated.

D. melanogaster
Furthermore, the database specific parameters they used in OrthoDB and OrthoMCL should be mentioned as well as the database version they used, e.g. OrthoDB2, or OrthoDB3, … , or OrthoDB6.
Age index for each stage-pair There seems to be a typo within the TAI formula. As defined by Domazet-Loso and Tautz and adapted to the author's notation the formula of the TAI should begin not TAI(l) TAI (1).

Results
The authors show based on their findings returned by Model-2-4 that the increasing specificity The reviewers are right to highlight the connection between our Response to Point #1: "increasing regulatory specificity" assumption and Raff's hypothesis. We have included a paragraph in the Discussion section to put Raff's hypothesis in the context of this work. A first remark is that our assumption is consistent with Raff's basic premise that modularity increases as the embryo develops. Second, Raff's hypothesis (see chapter-6 of [1], and in particular Figures 6.6 and 6.7) also states that the "level of interaction" or the "interconnectivity between (body) elements" is maximized at the phylotypic stage. This may be viewed initially as a contradiction between Raff's hypothesis and our increasing specificity assumption, which states that the density of regulatory interactions is maximized at the earliest stages of development. Note however that Raff's hypothesis was not stated in terms of gene regulatory interactions --he was referring more generally to "developmental flexibility", arguing that early development is flexible because it governs robust and general global patterning processes, late development is also flexible because "signaling events within the primordia are little influenced by events in other primordia", while mid-development (phylotypic stage) is least flexible because of the "high interconnectivity between elements that will later come to represent separate modules." If we think of "developmental flexibility" as the ability of an embryo to survive gene mutations and rewiring at different stages of the developmental process, Raff's hypothesis is actually consistent with our results regarding the "lethality probability" at each developmental stage (see Figs 6-B, 7-B amd 8-B). The lethality probability is maximized at the phylotypic stage, and it is significantly lower at early and late developmental stages, following the same pattern with Raff's "developmental flexibility".
First, to the extent of our knowledge there is no prior work that Response to Point #2: validates the assumption of increasing regulatory specificity as formulated in our directly paper (i.e., with the specificity of a developmental stage defined as the topological density of regulatory edges in that stage). We are planning to examine the validity of this assumption in the future using various approaches: First we plan to analyze the aforementioned developmental GRNs for the sea-urchin that have been directly examined by Dr. Davidson's lab at Caltech.
Second, we are in the process of examining other genomic metrics, which can be viewed as "proxies" to regulatory specificity. For example, we have defined new measures to quantify how specific a gene's expression profile in a developmental stage is compared to other developmental stages, which we have tentatively named as "developmental stage specificity index (DSI)". Similarly, we defined "tissue specificity index (TSI)", measuring the specificity index (DSI)". Similarly, we defined "tissue specificity index (TSI)", measuring the bias of a gene's expression in a specific tissue compared to many other tissues. Analyses of these measures using genomic data from Drosophila so far have revealed generally increasing trends, which are consistent with the idea of increasing specificity. However we are still working on generalizing these analytical tools to data sets generated by different methods and from different species. Importantly, our model focuses DGENs underlying development, and genomic data can obscure the true signals from genes constituting DGENs. Thus we are also working on examining these specificity measures using inferred DGENs. Additionally, we are developing other methods to approximate functional specificity. We hope to complete these analyses as a follow-up paper. The current paper provides theoretical motivation.
Third, as the reviewers also point out, there is a large body of prior work in developmental biology that confirms the increasing modularity in the developing embryo (for instance, see ) as well as the increasing specificity of the developmental process at Wagner ., 2007 et al the signaling or genomic level. The connections between this pattern of increasing modularity and the structure of the underlying gene regulatory networks are still not well understood however. We think that it would be misleading if we had argued that these modularity patterns provide a direct validation for our "increasing regulatory specificity" assumption. Please note that we have cited the paper by Wagner et al. at the Introduction and Model Description sections, providing some connections between this work and earlier work on the role of modularity in developmental biology.
Review: "Numerous studies investigated the phenomenon that increasing morphological complexity is correlated with increasing modularity and specialization of the underlying regulatory apparatus . Here the authors predict that the main condition leading to the appearance of the hourglass pattern is that the DGEN should gradually get sparser as development progresses, with general-purpose regulatory genes at the earlier developmental stages and highly specialized regulatory genes at the later stages, expressing the oldest genes during mid development. We suggest to also include findings investigating the functionality and age of the genes that have been shown to be active during mid embryogenesis to correlate potential outcomes of the DGEN prediction that evolutionary older genes should be concentrated at the waist of the hourglass with results obtained by recent studies . These aforementioned studies were able to show that during mid embryogenesis evolutionary old genes or evolutionary old processes that are shared within and between phyla are most active. In our opinion it would be very interesting to test the degree of intersection between genes predicted to be causal for the waist of the DGEN hourglass and evolutionary conserved genes that have been found to be expressed during mid embryogenesis. We  Response: This is a very interesting idea and we thank the reviewers for this suggestion. We are planning to pursue this investigation in a follow-up paper that will also include the previously mentioned analysis of sea-urchin GRNs as well as the DSI/TSI results. We believe that it would be distracting to try to include all these results in this first paper, given that the main focus of this work has been on the computational model and the computational results. We are also open to collaborate in this investigation with the reviewers or other researchers working in this area, if they are interested.
Review: "The model description is quite intuitive. However, as the validity of the whole study depends on the model parameters chosen, several points need clarification including the rationale and motivation to choose exactly these parameters. One of the most crucial parameters of the model is the regulatory specificity s(l), with 0 <= s(l) <= 1. We think, that adding a more detailed paragraph discussing how s(l) is being estimated, would be highly beneficial due to the importance for the parameter s(l) for the entire model and predictions.  Figure 9 a non-linear function s(l) = 0.9 -(0.8 / 1 + exp(\gammal) ). A gene g at stage l is modeled to act as upstream regulator for a gene g' at stage l + 1 with probability s'(l) = 1 -s(l). From this follows that each potential upstream regulator in stage l has the same probability s'(l) to regulate genes g' in stage l + 1. In other words, Model-1 assigns each potential upstream regulator with the same probability s'(l) = 1 -0.5 = 0.5 to regulate genes g' in stage l + 1, in Model-2: s'(l) = 1 -l/L = const., and for the non-linear function s(l): s'(l) = 1 -(0.9 -(0.8 / 1 + exp(\gamma -l) )) = const. It would be advantageous, if the authors could include motivations for the three cut-offs, especially how they derived the sigmoid-like mathematical function (Figure 10, legend caption) they later denote as non-linear function s(l). We do see that modeling the stage specificity with a constant value, an increasing value, and a non-linear function is plausible, but a clear motivation or derivation of the exact values would be beneficial for a better understanding of the modeling process.
The authors discuss that their major assumption is that the regulatory specificity increases substantially as development progresses, hence the DGEN becomes gradually sparser along the developmental time axis starting with s(1) ~ 0, … , s(L) ~ 1. Here the authors should include references to experimental observations that support their assumption. Since Models 2-4 rely on this assumption, published experimental studies in line with this assumption need to be referenced. Otherwise, the models applied would seem largely academic.
Include a clear statement how the exact values of the stage specificity (regulatory specificity) s(l) have been derived. State more clearly that s(l) is the crucial parameter determining the observed phenomena."

Response:
The regulatory specificity functions we experimented with are arguably the simplest one could think of. Specifically, we start (Model-1) with the simplest hypothesis that the regulatory specificity does not vary across developmental stages (in the paper we show results for s=0.5 but we have generated results for all values 0.1, 0.2, ... 0.9 and the conclusion remains that when the specificity is constant the network does not evolve into an hourglass shape). Then, we consider (Model-2) a linearly increasing specificity function. We do not argue that this is realistic or that there are experimental results that suggest this do not argue that this is realistic or that there are experimental results that suggest this linearity or even the increasing trend. It is just the simplest and most parsimonious model that leads to the emergence of an hourglass pattern. Finally, in Figures-9 and 10 we consider a more general specificity function that increases in a non-linear manner, providing a simple way to control the stage at which it gives the mid-range value 0.5. This allows us to examine how the shape of the specificity function s(l) affects the location of the hourglass waist. We hope that follow-up work will reveal experimentally the actual shape of the regulatory specificity function s(l) and there is no doubt that it will not be identical to any of these simple mathematical functions.
Please note that we have revised the paragraphs about "Model-1" and "Model-2" in the Simulation section to clarify these points.
Review: "For the probabilities P_{DL} and P_{DP} the explicit formula to compute the corresponding probabilities could also be included within the paragraph analogous to P_{RF}. This would enable a clearer reproducibility." : Actually there is no formula for the gene deletion and duplication probabilities.

Response
They are given probability values that do not depend on any other parameters. At each round of the simulation, a gene is duplicated with probability P_{DP}. If not rewired or duplicated, the gene is deleted with probability P_{DL}.
Review: "Hourglass score H: Here w(l) is defined as the number of transitioning genes in stage l. In Figure 3 and in section 'Hourglass shape' w(l) is referred to as stage width. In case the definition of w(l) is equivalent to the stage width, it would be beneficial for the reader to see the term stage width in the same sentence as the initial definition of w(l). A correct understanding of the stage width is important to understand the measurement of H." : The reviewer is right. We will use the "stage width" term when we first introduce Response the notation w(l).
Review: "Furthermore, the univariate Mann-Kendall statistic should be referenced to allow a non-statistical community to better understand this way of statistical modeling." : Please see how we addressed a similar comment (comment #5) in Dmitri Response Krioukov's review. : The phylostratigraphy technique as employed by Domazet-Loso et al. and Response Quint et al. involves the assignment of an evolutionary age to each gene in a given species, by tracing the most distant ancestral node containing at least one immediate daughter species with a detectable homologue. A phylogenetic tree of divergent species ranging from cellular organisms transitioning into simple and more complex eukaryota served as a reference to divide phyla radiating from consecutive ancestral points into distinct phylostratum layers. To determine the homologues, BLAST sequence similarity searches were done against complete genomes reliably annotated across the different phylostrata. Similarly, we obtained known homologues of each Drosophila gene from OrthoDB5 and OrthoMCL5 (that also use BLAST based searches) in the species depicted in Figure 4, and mapped them to six different "phylostrata" (marked in Figure 4). Akin to the method described in the above references, we assigned each gene to a specific "phylostratum" based on the farthest detectable homologue of that gene, and this is referred to as "age index" in the paper. This made it possible to compare the age index based results we obtain in Drosophila to the phylostratigraphy based assignments in Arabidopsis.

Review
Homologues for Drosophila genes were sourced from these two databases, to increase coverage across the species depicted in the tree in Figure 4. Since these two databases use slightly different parameters to identify orthologs, we examined whether the usage of specific databases introduce bias by analyzing data from each database separately and observed that age index of each gene remained the same.

Review: "
There seems to be a typo within the TAI formula. Age index for each stage-pair: As defined by Domazet-Loso and Tautz and adapted to the author's notation the formula of the TAI should begin TAI(l) not TAI (1)." : We have fixed this problem. Thank you. Response Review: "The authors show based on their findings returned by Model-2-4 that the increasing specificity assumption is the key property behind the developmental hourglass effect. We think that this statement is not clear enough. The authors could phrase more clearly what they mean by "key property behind the hourglass effect" inferred from Model-2-4." : The reviewer is right. We revised that paragraph as follows: Response "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 6 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." Review: "Age of genes: Here the authors motivate that defining the age of a gene as A(g) = i -t0(g) is plausible, because a rewiring event may give a gene a new function, at least in terms of its upstream and downstream regulators. This statement needs validation from the literature and from experimental studies that already tested this hypothesis for specific rewiring events." : We have added two references that provide some evidence to the previous Response hypothesis: ) and . Guet  (Figure 8c) that the evolutionary age at stage l follows the same pattern as the lethality probability: gradually increasing until the waist of the hourglass, and then gradually decreasing. This finding should be correlated with the findings by Galis and Metz and should be discussed more clearly in the Discussion." : The reviewer is right. The findings of Galis and Metz are highly relevant and Response consistent with our computational results regarding the lethality probability. We modified the paper as follows: In the Results section: "These computational results for the lethality probability across development are consistent with the empirical observations of Galis and Metz \cite{galis2001testing} about the increased mortality of rodents due to perturbations in the phylotypic stage (see Discussion section)." In the Discussion section: "The third prediction of the model is in direct agreement with the empirical observations of Galis and Metz regarding the increased mortality caused by perturbations during the phylotypic stage \cite{galis2001testing}. That study has shown, based on the teratological literature for rodent development, that disturbances in the phylotypic stage lead to much higher mortality than in other stages. Further, such disturbances lead to the co-occurrence of several distinct anomalies in the developing embryo and so the increased mortality cannot be due to a single particularly vulnerable process that takes place at that stage." Review: "In the first sentence of the last paragraph before the Discussion, "… and TAI profiles from Arabidopsis" should be referenced to reference 11 and not 12." : Thank you, we have fixed this problem. Response : Review "Discussion: The authors discuss that the observed pattern of conservation (developmental hourglass) may stem from fundamental organization principles and that the exact origin of these principles remains elusive. Here the authors should include existing studies that also aim at predicting a universal organization principle during development (modularization, etc.) and what exact problems are preventing studies to elucidate this universal phenomenon. We do note, that the authors cite studies that aim to explain the effects causing the early phase of the developmental hourglass and studies aiming 3 14 15 16 effects causing the early phase of the developmental hourglass and studies aiming to explain effects causing the late phase of the developmental hourglass , but we think that it should be explicitly stated that the model is designed to test aspects of Raff's initial hypothesis and that also more recent studies have been investigating the underlying reasons for the pattern in the early phase of the hourglass ." : The reviewer is right. We have revised the Discussion section to also include a Response paragraph about Raff's hypothesis (see our response to comment 1). We have also added another paragraph (see below) to discuss a very recent paper that proposes a different mechanism for the emergence of the hourglass effect.
"Recently, Friedlander have proposed a different mechanism for the emergence of et al. hourglass-like patterns in evolving biological or technological networks, based on a linear system formulation . They showed that if a system can be \cite{friedlander2014evolution} represented as a hierarchical and layered linear transformation of an input vector to an output vector, and the desired transformation matrix is rank-deficient, then an evolutionary process that selects that particular transformation can, under certain conditions, converge to an hourglass-like structure. It is not clear yet how to adapt this linear model in the context of inherently non-linear systems, such as gene regulatory networks." No competing interests were disclosed. Competing Interests: 14