Evolution of networks for body plan patterning; interplay of modularity, robustness and evolvability
A major goal of evolutionary developmental biology (evo-devo) is to understand how multicellular body plans of increasing complexity have evolved, and how the corresponding developmental programs are genetically encoded. It has been repeatedly argued that key to the evolution of increased body plan complexity is the modularity of the underlying developmental gene regulatory networks (GRNs). This modularity is considered essential for network robustness and evolvability. In our opinion, these ideas, appealing as they may sound, have not been sufficiently tested. Here we use computer simulations to study the evolution of GRNs' underlying body plan patterning. We select for body plan segmentation and differentiation, as these are considered to be major innovations in metazoan evolution. To allow modular networks to evolve, we independently select for segmentation and differentiation. We study both the occurrence and relation of robustness, evolvability and modularity of evolved networks. Interestingly, we observed two distinct evolutionary strategies to evolve a segmented, differentiated body plan. In the first strategy, first segments and then differentiation domains evolve (SF strategy). In the second scenario segments and domains evolve simultaneously (SS strategy). We demonstrate that under indirect selection for robustness the SF strategy becomes dominant. In addition, as a byproduct of this larger robustness, the SF strategy is also more evolvable. Finally, using a combined functional and architectural approach, we determine network modularity. We find that while SS networks generate segments and domains in an integrated manner, SF networks use largely independent modules to produce segments and domains.
Surprisingly, we find that widely used, purely architectural methods for determining network modularity completely fail to establish this higher modularity of SF networks. Finally, we observe that, as a free side effect of evolving segmentation and differentiation in combination, we obtained in-silico developmental mechanisms resembling mechanisms used in vertebrate development.
A major goal of evolutionary developmental biology (evo-devo) is to understand how multicellular body plans of increasing complexity have evolved, and how the underlying developmental programs are encoded in the genome and gene regulatory network (GRN).
Modern evo-devo research shows more and more that a shared developmental toolkit of signaling, adhesion and transcription factor genes are essential for the development of organisms ranging in body plan complexity from cniderians to arthropods and vertebrates [1]–[3]. Therefore the current paradigm is that body plans of increasing complexity are the result of increases in the complexity of regulation of this similar set of genes [1], [2], [4]–[9] combined with increases in the number of variants of certain developmental toolkit genes [10]–[12]. As a consequence, a strong focus in current evo-devo research is on changes in spatio-temporal gene expression patterns and the differences in architecture of the developmental networks generating them.
Network characteristics that are considered key for the evolution of increasingly complex body plans are modularity, robustness and evolvability. It is frequently argued that developmental GRNs are typically modular, i.e. that different functions are performed by largely independent network parts [2], [13]–[16], and that this is the key property responsible for both network robustness and evolvability. The idea is that modularity reduces pleiotropy, allowing for the malfunctioning of or tinkering with network parts involved in one function without producing failure in other functions [2], [13]–[16]. Although this reasoning sounds appealing and intuitively correct, little research has been done to explicitly test the roles and relationships of developmental network modularity, robustness and evolvability in the evolution of complex body plans. Indeed, we argue that it is currently unclear how modular developmental networks are, how such modularity evolves, and how this modularity looks.
Today, only a limited number of developmental GRNs have been studied in considerable detail. These studied networks are mostly involved in the patterning of a single organ or developmental phase, without detailed knowledge on their relationships with the rest of the developmental network [17]–[22]. As a consequence, although these networks have often been claimed to be modular, it is currently hard to fully assess the modularity of developmental networks.
Based on theoretical studies it has been argued that evolution should neither be expected to produce nor to preserve architectural modular networks. This follows from the fact that modular networks form only a small subset of the possible network architectures capable of performing a particular function [23]. Indeed, theoretical studies aimed at investigating the evolution of architecturally modular networks have had to use quite specific fitness targets to obtain modular networks [24]–[27]. On the other hand, it has previously been shown for other genome [28], [29] and network [30] properties that these may arise as a neutral side effect of the mutational dynamics rather than requiring an adaptive explanation. Similar suggestions have been made for network modularity [31], [32].
With regards to the appearance of modularity, note that in its most general sense network modularity is defined fairly functional -different functions are performed by largely independent network parts- but is currently most frequently measured entirely architectural -different modules of genes that are more densely connected with genes within the module than genes in different modules [33]–[35]. However, it is recently being suggested that functional or dynamic rather than architectural network modularity may be most relevant for network functioning and evolution [36]–[40]. Note that architectural and functional modularity do not necessarily overlap. This might mean that different, more functionally oriented methods to measure modularity are needed [36], [37], [40]. Recently, several such methods have been proposed, among which 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].
Thus, currently both the extent and shape of developmental network modularity remain unclear. In addition, it is not well known whether evolution of this modularity requires selection for robustness or evolvability or arises neutrally. The goal of the current study is to use computer simulations to investigate what type of network architecture and properties evolve during the evolution of complex body plan patterning. This will allow us to check to what extent evolved developmental networks are modular, whether network modularity is related to increased robustness and evolvability, and what exactly network modularity looks like. In our simulations we select for segmented and differentiated body plans. Segmentation and extensive anterior posterior domain differentiation are considered key innovations of the bilaterian clade, and have been extensively studied both experimentally and theoretically. This will allow us to compare our in-silico evolved developmental networks with actual biological patterning networks and results of previous simulation studies. Furthermore, by independently selecting for segmentation and domain formation we enhance the chances for modular networks to evolve.
We study the different types of evolutionary trajectories that arise, and compare them with respect to network robustness, evolvability and modularity and the type of developmental mechanism they produce. Quite interestingly, we find that there are only two distinct evolutionary strategies to evolve a segmented and differentiated body plan, each resulting in a distinct developmental mechanism. In one strategy, first most segments and only then domains evolve (SF strategy), while in the other segments and domains evolve more or less simultaneously (SS strategy). In addition, we show that in the SF strategy, a complex time transient is responsible for domain differentiation, while a genetic oscillator produces regular body segments. In contrast, in the SS strategy, a complex time transient generates both the body segments and domains. We find that imposed indirect selection for robustness causes the SF strategy to evolve much more frequently than the SS strategy. Furthermore, the SF strategy was also found to be more evolvable.
The different types of expression dynamics involved in segmentation and domain formation, together with the larger robustness and evolvability of SF networks suggests that they may also be more modular. However, frequently used, purely architectural modularity scores suggest that the two network types are equally non-modular. Pruning of non necessary network parts that potentially obscure architectural modularity did not change these results. Furthermore, changing model parameters such that less densely connected networks evolve also did not produce architecturally modular networks. Therefore, we also used a more functionally oriented method. Specifically, we take into account the fact that the networks generate both segments and domains and investigate whether or not there are relatively independent network parts responsible for these two processes. Using this approach we could demonstrate that while SS networks generate segments and domains in an highly integrated manner, SF networks generate segments and domains in a more modular manner.
Our results show that evolved developmental networks are not necessarily highly modular, robust or evolvable. However, upon significant selection for robustness, networks that are more modular, robust and evolvable will dominate. Our results thus confirm the relationship between modularity, robustness and evolvability. Our results also show that the type of modularity that evolved could not be detected by frequently used, automated, purely architectural algorithms, but required a more functionally oriented method. Beslon recently reported similar results [40]. Importantly, these results suggest that for the detection of biologically meaningful modularity purely architectural methods are less suitable and approaches (also) taking into account network dynamics and function should be preferred.
Intriguingly, we find that the patterning mechanism employed by our SF networks shares key characteristics with vertebrate somitogenesis and axial patterning, even though this was not a specific aim of our study or explicit part of our model design.
Here we provide a succinct description of the methods used, for a more elaborate description we refer to Text S1.
Briefly, we use an individual based, spatially embedded model of a population of evolving embryo-organisms (Figure 1). The organisms consist of a one dimensional row of 100 cells, similar to the approach followed in [41]–[43]. The organisms have a genome that contains genes coding for transcription factors (TFs) and upstream regulatory regions with transcription factor binding sites (TFBS) [44], [45].
Genes have a certain type, indicated with a number ranging from 0 till 15. There can be multiple genes of the same type. The gene types can be subdivided into a few functional categories. Gene type 0 is a maternal gene. Its expression is not controlled by the organism, but instead is imposed to give rise to a morphogen wavefront. This wavefront moves from the anterior to the posterior of the embryo, switching the expression from gene type 0 from a level of 100 to 0. Gene type 5 is a gene that the organisms can use to indicate the boundaries of body segments. Differential expression of gene types 8 till 15 can be used to subdivide the body into functionally different regions (domains). Finally, gene types 1 till 4, 6 and 7 are general transcription factors. By assigning gene type 5 to segmentation and gene types 8 till 15 to differentation the evolving segmentation and differentiation processes are not forced to be coordinated but can in principle use completely disjoint sets of genes.
The genome codes for a gene regulatory network, with genes corresponding to nodes, and TFBS defining the activating and repressing influence of genes on each other. These regulatory links have a non-evolving impact strength of +1 or −1, respectively.
At the beginning of development, gene expression in each cell of an organism is initialized with gene types 1 to 4 having an expression level of 100 and all other genes having an expression level of 0. Subsequent gene expression dynamics and protein levels are governed by the GRN and are modeled with ordinary differential equations using a similar approach as in [41].
The gene expression pattern present at the end of development is used to determine the number of segments and domains the body is patterned in. A segment boundary is defined as a position in space where the expression of the segmentation gene switches from a high to a low level or vice versa. A domain is defined as a region in space where cells express the same combination of differentiation genes at a high level. The minimum length for a segment and domain is 7 cells, allowing for a maximum number of 14 segments and domains. To ensure stable differentiation, we compare gene expression at the end of development with that 20 time steps before. For each cell that has different gene expression levels at these two time points a fitness penalty is applied. In addition, to prevent excessive genome growth small fitness penalties are applied for each gene and TFBS present in the genome (See Table S1 in Text S1).
At the start of evolution the population is initialized with a group of 50 identical organisms in a field of size 30×30. These organisms have a genome containing a single copy of each gene type in a randomized order and with an average of 2 TFBS, randomly drawn from the possible types of TFBS, upstream of each gene. Evolution occurs through mutations on the genome and fitness dependent reproduction. We apply gene duplications and deletions, TFBS duplications and deletions, and changes in the type and weight (activating or repressing) of TFBS. Note that in contrast to some previous studies [41]–[43] we do not evolve gene expression rates, protein decay rates, or TF activation and inactivation threshold levels here. Tournament selection is used to determine which organisms may reproduce. Death occurs with a constant probability of 0.5. After an initial transient population sizes plateau at around 600 individuals.
As explained, we are interested in the robustness, evolvability and modularity of the evolved developmental GRNs. To give evolution the freedom to evolve networks producing segments and differentiation domains either in a modular or integrated manner, we choose our fitness function such the number of segments and domains contribute independently to fitness (i.e. we use rather than e.g. ). As a side effect of this choice, evolution is also free to evolve only segments or domains, rather than both. For our analysis we select those simulations that were successful in evolving both a significant number (≥7) of segments and domains.
To determine the evolutionary history of a developmental mechanism and its underlying GRN we traced the ancestry of the final fit evolved individuals.
We performed a total of 50 simulations using the default parameter settings of our model (see Table S1 in Text S1). We analyzed the networks that successfully evolved segmentation and differentiation in terms of evolutionary strategy followed (whether segments and domains evolve sequentially, simultaneously, or something in between), network size and architecture (number of genes and connections, positive feedback loops, attractors) and generated developmental dynamics (type of spatiotemporal gene expression patterns and how this generates segments and domains).
Furthermore we evaluated the robustness, evolvability and modularity of different evolved network types. First, to determine robustness of different evolved network types we performed three additional series of 50 simulations. We increased mutation rate, added gene expression noise, or added variability in morphogen wavefront speed (see Text S1). From the frequency with which the different evolutionary strategies (SS or SF) occur we determine their relative robustness. Second, we performed a total of 140 simulations to find how network types differ with respect to evolvability. Here, we first performed 20 simulations with a fitness target of 6 segments and 6 domains. From these we selected 6 successful networks that differed in type (SS or SF). These were each used as a starting point for 20 independent simulations with a fitness target of 9 segments and 9 domains (see Results). From differences in rates of success of evolving to this second target we determine the relative evolvability of the different network types.
Finally, we determined the modularity of the different network types. Here we used a range of approaches. First, we determined the architectural modularity of the evolved networks using algorithms that try to find the optimal modularity score or Q value for a network. To ensure that our results were not biased by the particular details of the algorithm used, we used two different methods applying different heuristics. The first uses Newman's leading eigenvector method to determine optimal modularity [33], [34], the second method uses a random walk approach to determine Q values [35]. Furthermore, to allow interpretation of the thus found Q values, we determined Q values for not only random and architecturally modular networks, but also for neutrally evolved networks. These neutrally evolved networks serve as a benchmark against which to test whether there is selection for architectural modularity in our simulations.
However, architectural network modularity can easily be obscured by the presence of non-functional or redundant genes and regulatory interactions. Therefore, we pruned the original evolved networks to a minimal essential core network (see Text S1) and also determined Q values for these core networks. Furthermore, architectural modularity may be obscured by the particular model parameter setting used, when these tend to cause the evolution of densely connected networks. To determine whether this was the case, we performed 3 additional series of simulations in which the impact of TFBS deletion rates on modularity was tested. In the first two series, TFBS deletion rates were increased either twofold or fivefold, while all else was kept the same as in the default simulations. In the last series of experiments, a core network with a relatively high Q value was selected from the set of default simulations. This core network was subsequently taken as a starting point for continued evolution simulations with a fivefold higher than normal TFBS deletion rate.
Finally, as an alternative to these automatic, purely architectural methods of determining network modularity, we also assessed modularity in an alternative way. Here we used the core networks as a starting point to determine the minimal networks needed for either segmentation or differentiation alone. To determine how modular a network is we subsequently asses three points. First, we check how well the minimum networks are capable of autonomously reproducing the original segment or domain pattern. Second, we determine how well they can produce one thing (segments) without as a side effect also accidentally producing the other thing (domains). Finally, we assess the amount of overlap between the two minimum networks. Thus, we assess how functionally autonomous and how functionally and architecturally independent these network parts are. The method thus takes into account prior knowledge about network function (they generate both segments and domains) and considers both functional and architectural aspects of modularity. If the minimum segment and domain networks function are good at reproducing either only the original segment or the original domain pattern and contain only a few overlapping genes and connections, we will classify the network as modular. In contrast, Q value based algorithms may fail to detect modularity if modules share not only connections but also a few genes.
Figure 2A schematically shows the phase space of possible evolutionary trajectories of evolving both segments and domains. In it we show 3 theoretically possible extreme trajectories: 1) all segments evolve before domains evolve; 2) the opposite, all domains evolve before segments evolve, 3) the intermediate, segments and domains evolve simultaneously.
In our analysis we focus on those 30 simulations (out of the total of 50) in which ≥7 segments and ≥7 domains evolved. We find that in 10 of these simulations (33%) first most segments and then domains evolved. In Figure 2B the evolutionary trajectory of 5 of these simulations is shown. Each trajectory shows the maximum number of segments and domains in the population as a function of evolutionary time. In the 20 other simulations (67%) segments and domain numbers increased more or less simultaneous over evolutionary time. Figure 2C shows the trajectories of 5 of these simulations. None of the simulations first evolved most domains and then segments.
A detailed overview of the results of the 10 SF simulations and 20 SS simulations can be found in Tables S4–S9 of Text S1. These results are summarized in Table 1.