Evolution of networks for body plan patterning; interplay of modularity, robustness and evolvability

When comparing network architecture, we find that SF networks are simpler, with similar numbers of genes but significantly lower connectivity. With regards to the network's developmental output, we find that the two alternative strategies attain very similar overall fitness levels. However, SF type networks produce body plans with more segments then domains, whereas the SS type networks do exactly the opposite. In addition, the segments produced by SF networks are much more regularly sized than those produced by SS networks. Indeed, the developmental gene expression dynamics generated by the two network types differ significantly.

Figure 3 shows final evolved networks together with the generated intracellular gene expression dynamics, developmental space-time plot, and the final gene expression pattern for both an example SS (Figure 3A) and SF (Figure 3B) network.

We see that the evolved SS GRN is quite complex, containing 24 nodes and 72 connections (Figure 3A, top row). The network produces a complex time transient of gene expression (Figure 3A, bottom row) that upon passage of the maternal morphogen wavefront (gene type 0, arrow) is converted into a stable gene expression pattern. We furthermore observe that the gene types that become stably expressed at a location depend on the time when the wavefront passes. As a consequence, the complex time transient is translated into a temporally stable, but spatially diversified gene expression pattern (Video S1). The space-time plot (Figure 3A, top row) shows another representation of this process. We recognize the anterior to posterior progression of the morphogen wavefront as a distinct diagonal pattern, and see how it transforms the time varying gene expression into a stable spatial pattern (Video S2). If we look at the gene expression pattern at the end of development (Figure 3A, top row) we see that a spatially alternating expression of the segmentation gene (gene type 5) produces 7 body segments of different sizes. The combination of spatially varied expression of the identity genes (gene types 8 till 15) produces a total of 10 domains, also of varying sizes.

The SF network is indeed simpler, containing 23 genes and 57 connections (Figure 3B, top row). The networks produces a complex time transient of gene expression (Figure 3B, bottom row) in which a subset of genes (gene types 2, 5, 7, 10, 12, 13 and 15) display an oscillatory dynamics that we did not observe for the SS network. As for the SS network, the passing by of the morphogen wavefront converts the time-varying gene expression into a stable, spatially varied expression pattern (Video S3, Video S4). However, in this case the oscillatory dynamics of genes 2, 15, 5 and 7 are translated into a regular, alternating expression pattern, allowing gene type 5 (segmentation gene) to produce 12 regularly sized segments (Figure 3B, top row). This mode of producing segments resembles the process of somitogenesis in vertebrates. In addition, the non-oscillatory dynamics of genes 3, 6, 8 9 and 11 are converted to 4 continuous, staggered expression regions (Figure 3B, top row). This expression pattern resembles the typical expression pattern of Hox genes along the anterior posterior axis of bilaterian animals. As genes 8 till 15 all are identity genes, the combination of the alternating expression of gene 13 and 15 and the continuous staggered expression of genes 8, 9 and 11 produce a total of 7 different domains (if multiple regions express the same set of identity genes only the first counts as a domain).

Similar results were found for other SS and SF networks. Thus, while SS networks use a complex time transient to produce both segments and domains, SF networks use a similar complex time transient to produce domains, while using oscillatory dynamics to produce regularly sized segments. In later sections we discuss further details of these developmental dynamics in the context of network modularity.

We found that under the default parameter settings (Table S1 in Text S1) the SS strategy evolved more frequently than the SF strategy. Next, we investigated how the propensity of the two evolutionary strategies is affected by adding noise to our simulations. Previous research has shown that robustness evolves as a result of increased mutational or gene expression noise [46]. Here we thus assume that increased noise, independent of the type of noise, produces indirect selection for robustness. By assessing the frequency with which the different strategies evolve under increased noise we investigate which of the two strategies is more robust.

We performed 3 series of 50 simulations. In the first series mutation rate was increased by a factor 10. In the second the propagation speed of the maternal morphogen gradient was varied between individuals within a 30% range. In the third series 5% gene expression noise was incorporated. Table 2 shows the percentage of successful simulations and how often the different evolutionary trajectories were followed. Note that we did not observe any additional types of evolutionary trajectories, i.e. first evolving domains and then segments. We see that for all 3 additional series of simulations a shift occurred from SS as a dominant evolutionary strategy to SF as a dominant evolutionary strategy. Thus indirect selection for robustness favors the SF type networks, suggesting that these are more robust.

Next we determined whether the two network types also differed in evolvability. It is frequently thought that a special selection regime is required for the evolution of evolvability [2], [13]–[16]. An often used approach is to impose indirect selection for evolvability by alternating between different selection regimes [44], [47]–[49]. Clearly, such a back and forth alternation between selection criteria is hardly realistic in a developmental context. However, it has been shown that robustness and evolvability of GRNs is strongly correlated [50], [51]. It is thus interesting to investigate whether the differences in robustness we observed between the two evolutionary strategies are correlated with differences in evolvability. Specifically, we tested for differences in the evolutionary potential of the two network types for evolving new segments and domains.

To do this, we first performed 20 simulations in which we selected for 6 segments and 6 domains (Figure 4). From these simulations we selected the successful ones. Next, we selected 3 SF and 3 SS simulations. From these 6 simulations we extracted the genome of a finally evolved, fit individual. Each of these 6 genomes were used as input for a series of 20 independent simulations in which now selection for 9 segments and 9 domains was imposed. Finally, we compared the success rates of these 6 series of simulations (Table 3) and whether these differed significantly (pairwise t-test) (Table 4).

We see that simulations started with SF type genomes have a considerably higher success rate than simulations started with SS type genomes (Table 3) and that these differences are significant (Table 4). In contrast, simulations started with different genomes but of the same strategy type have much more similar success rates (Table 3), differences being not or hardly significant (Table 4). Differences in success rate are thus not due to random differences between genomes from different simulations, but rather are due to the more fundamental differences between genomes evolved following SF versus SS type evolutionary trajectories. Clearly, genomes evolved in a SF trajectory have a higher evolvability for inventing new segments and domains. These results imply that increased network evolvability can occur as a byproduct of selection for robustness, rather than requiring selection for evolvability itself.

Note that it remains an interesting question for further research whether other types of evolvability have also increased. Particularly relevant would be whether the ease with which segmentation and differentiation patterns are maintained if embryo size changes, the ease with which celltypes within domains can be changed, or the ease with which segment and domain numbers can decrease are also increased.

Next, we determined the modularity scores for both SF and SS networks. Based on the higher robustness and evolvability of SF networks together with the fact that they use distinct expression dynamics to generate segments or domains, one would expect SF networks to be more modular. In contrast, independent of the method used we found for both the SF and SS networks an average modularity score of around 0.29 (see Table 1, Figure 5A and Tables S4 and S7 in Text S1).

To be able to interpret the meaning of these similar modularity scores, we also determined modularity scores of randomly generated networks, neutrally evolved networks (without a fitness target) and manually designed, architecturally modular networks (see Figure 5 and Text S1). Independent of the modularity algorithm used we found Q values of around 0.29 for random networks (Table S2 in Text S1). For modular networks we found Q scores of around 0.65 (Table S3 in Text S1). Interestingly, for neutrally evolved networks we obtained Q values of around 0.45 (see Text S1). This demonstrates that the mutational dynamics alone causes a significant bias towards architectural modularity, without any present functionality.

If we compare the modularity scores of our evolved networks to these data we see that they are only slightly higher than those of random networks and significantly lower than those of neutrally evolved networks. Thus, selection clearly does not increase the type of architectural modularity measured by the used methods in either the SF or SS networks. This result is further confirmed by the observation that during evolution no significant increases in Q values are observed (see Figure S8 in Text S1).

To determine whether non-functional and redundant network parts obscure an underlying architectural modularity we also determined Q values for the core networks derived from the evolved networks. Similar to the original networks, the SF core networks have significantly less connections than the SS core networks (Table 1, also compare Figure 6A and 6B, top rows). However, again Q values of around 0.3 were obtained for both SF and SS type networks (Table 1, Figure 5B, and Tables S4 and S7 of Text S1).

Next, to check whether the parameter setting used causes a bias towards densely connected non modular networks we performed three series of additional simulations in which TFBS were increased. In the first two series we performed simulations that are the same as before, but with two times or five times higher TFBS deletion rates. This did not result in networks with significantly higher architectural modularity scores, independent of whether the original or core networks were evaluated (Table S10 of Text S1). In the final series, we started simulations with a previously evolved core network and continued its evolution under five times higher TFBS deletion rates. As a starting core network we took the core of the SF network shown in Figure 6B, as it has a relatively high Q value compared to average found values. Again, no significant increases in Q values were observed (Table S11 of Text S1). We conclude that frequently used, purely architectural methods to determine network modularity suggest that SS and SF networks are equally non-modular.

Summarizing, SS networks generate both segments and domains from a complex gene expression time transient, whereas SF networks use a complex time transient to generate domains and oscillating dynamics to generate segments. Furthermore, SF networks are more robust and more evolvable. Still, no differences in network modularity were found using frequently used, purely architectural methods. The question thus is whether SF networks indeed are not more modular than SS networks, or that the methods we used above perhaps fail to uncover certain types of modularity.

Recently, several alternative, more functionally oriented methods to asses network modularity have been suggested. Examples are the clustering of genes with similar expression in network attractors [36], or with similar knockout effects [40], or with a function in the same specific process [37]. Here we also took such a function based approach. We use the fact that networks were evolved to produce both segments and domains, and our observation that SF networks use different dynamics to generate segments or domains. We determine the minimum networks needed for either segmentation or differentiation alone to asses network modularity in an alternative manner (for details see Text S1).

Figure 6A shows the core, minimum segment and minimum domain networks derived from the example evolved SS network, together with the developmental dynamics and final gene expression patterns they generate. The core network has 21 genes and 64 regulatory connections (Figure 6A, top row), and the minimum segment network (Figure 6A, middle) still contains 16 genes and 38 connections. It produces a segmentation gene expression pattern that is shifted relative to the original pattern and capable of producing 6 of the original 7 segments. Furthermore, even though it is only required to produce segments, as a side effect it also produces 5 of the original 10 domains. The minimum domain network (Figure 6A, bottom row) consists of 17 genes and 36 connections. It generates an identity gene expression pattern that is very different from the original, and is capable of producing only 7 of the original 10 domains.

Summarizing, the minimum segment network produces a significant number of domains as a side effect of producing segments, and the minimum domain network performs rather poorly at reproducing the original domain pattern. We conclude that the evolved network is rather non-modular. Instead segments and domains are generated in a highly integrated manner. Indeed, if we compare the two minimum networks, we see that only 2 genes are unique for the minimum segment network and only 3 genes are unique for the minimum domain network (light blue), all other genes are used both for segmentation and domain formation.

Thus, to understand the mechanism behind body plan patterning we should look at the core network, which generates segments and domains in an integrated manner. The observation that a complex gene expression transient is translated into a spatial differentiation pattern suggests two things. First, the core network contains multiple attractors allowing for different stable cell types. Indeed, we see a total of 6 positive feedback loops, essential for attractor formation [52]–[54], in the core network (Figure 6A, top row). Second, the network produces complex and slow expression dynamics, allowing different times of wavefront passage to cause convergence to different attractors. In Text S1 we further explain this developmental mechanism and contrast it with the one described by Francois and Siggia in which a slow timer gene controls a linear sequence of gene activations [43]. Finally, to understand how segments arise as part of this process we study the regulation of the segmentation gene. We see that genes 14 and 15 activate and gene 8 represses gene 5 (Figure 6A, top row). Thus, the spatially alternating expression of gene 5 arises from integrating the inputs of these three genes. Each segment is thus generated by a different combination of regulatory inputs, in a very crude manner resembling Drosophila segmentation.

In Figure 6B we show the core, minimum segment and minimum domain networks derived from the example evolved SF network, combined with the developmental dynamics and final gene expression patterns they generate. We see that, in contrast to the SS network, the core network contains only 18 genes and 36 connections (Figure 6B, top row) and the SF minimum segment network contains only 7 genes and 10 connections (Figure 6B, middle row). The latter produces an oscillatory expression pattern that the passing wavefront transforms into a spatially alternating pattern, producing 11 of the 12 original segments. The segmentation network can be decomposed into a part responsible for generating bistability and a part responsible for producing oscillations, which in combination enable the translation of temporal oscillations into spatial stripes (see Figure S12 in Text S1). Also in contrast to before, the SF minimum segment network does not produce any domains as a side effect.

The minimum domain network (Figure 6B, bottom row) contains 13 genes and 15 regulatory interactions. It produces a complex gene expression transient that generates 4 continuous staggered expression domains, very similar to the Hox-like domains produced by the original network. The SF minimum domain network uses the same developmental mechanism as we discussed before for the core SS network to generate different stable expression domains. In this case, the network contains 3 positive feedback loops: a loop consisting of genes 3, 6 and 11, and positive autoregulation of genes 8 and 9. Further details of this developmental mechanism can be found in Text S1.

However, we also see that the spatially alternating expression of identity genes is not reproduced by the minimum domain network (compare Figure 6B top and bottom row), causing 4 rather than 7 domains to be formed. This shortcoming is due to the standard removal of the segmentation gene from the minimum domain networks (see Methods and Text S1). In the original network the segmentation gene causes genes 13 and 15 to have an alternating expression pattern that contributes to the number of domains. Note however that in contrast to the SS network, the subset of domains that is generated by the minimum domain network corresponds well to those generated by the original network, rather than being shifted in position or expressing different gene combinations.

In contrast to the SS network, the SF minimum networks are thus well capable of generating either segments or domains autonomously and independently. Indeed, the dynamics and expression patterns generated autonomously by the minimum segment and domain networks to a large extent add up to the behavior of the original network. The only clear exception is formed by a subset of identity gene expression domains that are dependent on the segmentation process (see above). However the correspondence is not perfect. For example, the minimum segment number generates a first segment that is too wide and a total of 11 rather than 12 segments (compare Figure 6B top and middle row). In addition, the expression patterns of genes 3, 8, 9, 6 and 11 produced by the minimum domain network are somewhat different than those produced by the original network (compare Figures 6B top and bottom row) (for more details see Text S1). Apparently some of the network parts present in the core network but not in the minimum segment and domain networks are needed both for some segmentation dependent domains and for some additional fine tuning of the segmentation and differentiation processes.

Also in contrast to the SS minimum networks, the two SF minimum networks together contain 17 unique genes, of which only 3 (colored yellow) are shared between the two networks. Together these observations demonstrate that the SF minimum segment and minimum domain networks are modules that are largely independently capable of segmenting and differentiating the body plan. We conclude that the SF network is significantly more modular than the SS network. As discussed above, the SF network is not completely modular: some domains are segmentation gene dependent, some fine tuning between segmentation and differentiation is needed, and a few connections and genes are shared between the minimum networks.

Given the observed modularity of the SF minimum segment and domain networks, we next investigated whether the earlier used purely architectural modularity methods are capable of retrieving this modularity. Put differently, if we sum the minimum segment and domain networks into a single network, do the Q value methods retrieve these modules and assign the summed minimum network a high Q value? Perhaps surprisingly, modularity scores for the summed minimum networks are still lower than those of neutrally evolved networks ( and , see Text S1). In addition, found architectural modules are inconsistent between the two used methods and unrelated to the above discussed segmentation and differentiation modules (for details see Figure S11 in Text S1).

We suspect that apart from not taking functional aspects into account, an important problem of the architectural modularity algorithms is that even a limited amount of overlap in genes used between functional modules causes them to not be recognized as architectural modules. In contrast, with our alternative method we simply classify a network as being more modular if fewer overlaps between minimum segment and domain networks are found.

As a final part, we investigated whether the differences between the SS and SF evolutionary and developmental strategies are reflected in further differences between their evolutionary dynamics.

Figure 7A shows the evolutionary dynamics of segment and domain numbers, attractor numbers and genome size along the line of ancestry leading to fit individuals at the end of the example SS simulation. The arrow indicates the position along this ancestral line of the individual we have analyzed in detail in Figures 3A and 6A. The inset shows the initial phase of evolution. As expected, we see a gradual and simultaneous increase of segment and domain numbers. In addition, the increase in domain numbers appears correlated with an evolutionary increase in GRN attractor numbers. In contrast, we observe no clear correlation between genome size and increases in segment and domain numbers. Instead, genome size shows intermittent periods of expansion and contraction. Interestingly, core genome size is only slightly smaller and has a similar dynamics. This suggests that information is stored in a diffuse, distributed manner, so that when the amount of encoded information (number of segments and domains) increases, the size of the core genome does not change so much.