|Home | About | Journals | Submit | Contact Us | Français|
Advanced genomics tools enable powerful new strategies for understanding complex biological processes, including development. By extension, we should be able to use these methods in a comparative fashion to capture evolutionary mechanisms. This requires a capacity to go deep and broad, to analyze developmental gene regulatory networks in many organisms, especially nontraditional models. As we usher in a new era of next-generation GRN (gene regulatory network) analysis, it is important to ask how to evaluate the evolution of network interactions. Particularly problematic, as always, is defining “independence”: Are two character traits found together because they are functionally linked or because of historical accident? The same basic question applies to understanding developmental GRN evolution. However, the essential difference here is that a GRN defines a causal chain of events. An understanding of causal relations—how Genes A and B work in concert to drive expression of Genes C and D to create a new Territory E—gives hope for establishing “trait independence” in a way that purely correlative arguments—the association of the expression of Gene D in Territory E—never could. Insight into causality provides the key to interpretation, as seen in this simplified scenario. Real-world networks bring new degrees of complexity, but the elucidation of causal relations remains the same. Has the day arrived when a single laboratory has the wherewithal to conduct multiorganism gene network projects in parallel? No. However, we argue that day is closer than one might suppose. We describe how a speedboat GRN project in one’s favorite nonmodel organism(s) might look and provide a framework for comparative network analysis.
“A pattern corresponding in a set of compared traits is homologous, if after a common evolutionary origin, the pattern was maintained along diverging lineages by robust pattern transmission.” (Szucsich and Wirkner 2007).
Generally speaking, a pattern is a set of predictably recurring events or objects. Obviously, many developmental processes fall into this description and nowadays, developmental patterns of gene expression are commonly used for comparisons between species. Gene-expression patterns are employed to infer homologies and the shared evolutionary origin of cell types, structures, or body parts, but patterns can also be defined as a system of connections and interactions between different components (Szucsich and Wirkner 2007). In light of this interpretation, different kinds of interactions and connections can help characterize a particular pattern, and processes can be interpreted as spatiotemporal patterns (Szucsich and Wirkner 2007). Consequently, gene regulatory networks (GRNs) can be interpreted as highly dynamic patterns. This means GRNs are evolutionary characters that can be homologized. When interpreting GRNs as patterns, the genes or gene products and their interactions in the network are the components of the pattern (Szucsich and Wirkner 2007). They interact dynamically, spatially, and temporally by activation or repression.
Resolving the GRNs that underlie developmental processes in several species would enable comparison of these networks and identify similarities and differences in their components and interactions. Similarities in the GRN between two species may indicate that the pattern has been maintained along both lineages and thus has a common evolutionary origin. However, in order to determine whether any given GRN component, for instance, a group of coregulated genes, is the product of conservation or convergence, it is important to show that the investigated elements represents truly complex patterns. The more complex the correspondences between two or more components, the more plausible is a hypothesis of common origin. Only if the components are independent can it be said that the investigated pattern is complex, thus favoring the homology hypothesis. If the components of a pattern can be changed, they are independent of each other (Scholtz 2005). Thus, it has to be shown that individual genes or gene products and their interactions are not just functionally linked and dependent on each other, but that they are, or could be, altered.
A famous example is the hedgehog pathway, which is involved in several developmental processes across eumetazoans (e.g., Matus et al. 2008). The interactions between proteins of the hedgehog family with Patched and Smoothened and the signal transduction by Cubitus interruptus/Gli are highly conserved (Matus et al. 2008). The members of this pathway are functionally dependent on each other and cannot be altered. Thus, in an evolutionary comparison, the hedgehog pathway as a whole should be treated as one component of the network. Naturally, we must be open to the possibility of varying degrees of dependence and independence as opposed to speaking in terms of absolute dependence. Nonetheless, a high degree of dependence among factors might be a general feature of signaling systems.
At its most fundamental basis, developmental biology asks how certain genes are designated for expression, and others not, and at the heart of this central challenge lies the problem of scale. This is a familiar problem to any biologist no matter the system. With respect to GRNs, the factors making up the network may number in the hundreds, whereas the dynamic, i.e., multi-timescale, nature of biological regulation further confounds simple analyses. Although it may seem an academic point, to us this suggests the limited effectiveness for analyses done at any one scale, for instance to take extreme cases, a “systems” level approach that does not take into account the properties of the components of the system or, conversely, a “reductionist” view that never defines the components that make up the system. We would, instead, characterize our approach as one of “scale integration,” integrating data sets from multiple scales, namely, using our ability to capture the global network to define which genes make up the control system, then to focus progressively on the most relevant factors and interactions driving discrete outcomes.
There are three recurring themes to our scale-integration approach, which can only be touched on here. First, we rely on an ability to perform temporal modeling, as opposed to depending on a static, i.e., time-free, depiction of network relationships. The importance of incorporating time into modeling ultimately derives from the fact that biological regulation is itself highly dynamic. Many interactions operate on scales of milliseconds to seconds (think phosphorylation cascades), whiereas others, such as gene regulation, are measured in many minutes to hours, and various chromatin states are even more long-lived than that. Modeling that can capture high-order interactions better reflects the system in action.
Second, scale integration is fundamentally about balancing complementary prospective analyses. This can be thought of in terms of false positives and false negatives or, similarly, sensitivity (reduced false negatives at the expense of some false positives) and specificity (low false positives at the expense of false negatives). Balancing sensitivity and specificity over all assays and within particular phases represents a major part of dealing with the challenges of scale.
Third, we emphasize the utility of qualitative modeling techniques. In the past, there has been a tendency toward quantitative modeling for the purpose of fitting data, but a strict numerical straightjacket is often neither practical nor informative. With large-scale networks and very few of the kinetic parameters known, this is a particularly sterile exercise. Qualitative modeling on the other hand better captures many biological phenomena, such as “Factor A represses Gene B,” as well as higher order functions, such as bistability and switching behavior, a good example of which is found in Cosentino et al. (2012). A mathematical framework that allows us to incorporate both quantitative and qualitative data through the use of inequalities, plays an important role in integrating heterogeneous datasets. Integration of diverse, complementary datasets is a central feature of integration of scale.
What follows is a brief overview of the general scale integration procedure (more details can be found in Smith et al. 2011), then an outline of how comparative network analysis might proceed. Until recently, such an approach was not feasible. The good news is that neither experimental nor computational capacities are as limiting as before due to advances in the past decade. With many previous constraints removed, we try to restructure research objectives in line with the classic parameters of observation, and the generation and testing of hypotheses. Thus, we first take large-scale surveys to define the factors making up the control system. This is the “observational” phase: From transcriptome to interactome in the language of -omics. Then, more focused analyses based on controlled experimentation are used to resolve network topology, i.e., to generate falsifiable hypotheses (interactome to network). Finally, targeted cis-regulatory analysis and fine-scale kinetic modeling combine to test hypothesized network functions (network to subcircuit).
How can we identify the most relevant factors in a network? Indeed, how do we even define relevance? Advanced genomics platforms offer unprecedented powers of observation in this respect. First, we can simply ask which genes are expressed, in other words, define the transcriptome. This might reduce the field from, say, the ~30,000 genes encoded in a typical animal genome to perhaps 3000: A big step but still not enough, practically speaking, to provide a focus for systematic genetic perturbations. Using bioinformatics tools to screen for annotated “regulatory genes” might also help, although success with that approach is not guaranteed even for organisms with well-annotated genomes, never mind for new or emerging model systems. Again, this might help some, but not enough to get to the 50–150 or so genes we can expect using established systems as a numerical guide. Although this number still seems large, it is a small enough set to pursue the next phase: targeted genetic perturbations. However, transcriptome sequencing with annotation of genes is unlikely to achieve this step.
Far better would be a less biased method for defining genes encoding the active set of regulatory factors. One way that has been rather well explored, particularly in nondevelopmental systems, such as those investigated in human cancer research, is the use of various algorithms to measure the statistical dependencies of the expression level for each expressed gene with that of every other (e.g., Margolin et al. 2006). This process of “reverse-engineering” returns an interaction map, or interactome (Fig. 1), which is a rough guide to how genes might be functionally related. Being defined by correlative measures, the interactome is a nondirected graph, i.e., whether Gene A regulates Gene B or vice versa (or both regulate each other) is not known. Nor does an interactome provide enough information to resolve many types of interactions, for instance, when a gene product regulates another gene both directly and indirectly through other factors. This is a common and functionally important motif that we address below.
Nevertheless, an interactome can give us important insights. In particular, even a low-resolution interaction map should tell us which genes are more likely to act as regulators of many others. These will be the hubs, as in Fig. 1. The more “central” genes, literally the ones found in the center of the interactome graph, represent the genes connected to the most genes that are themselves connected to the most genes. In other words, these are the regulators of other regulatory genes, the most relevant components of the control system. Nonregulatory genes, in contrast, are those at the periphery, perhaps connected to few genes, which are in turn connected to zero or few genes.
Although some peripheral genes might be canonical regulatory genes, including those encoding transcription factors, an interactome placing them at the periphery indicates little relevance to the control system.
Although this method has worked well in other fields, how well will it work for developmental systems? First, there is always the possibility of combining approaches to improve detection. For instance, the interactome may indicate that a gene of particular interest in other organisms is not central, and therefore not relevant to the control system. At this stage, it would make sense to include this gene and any others suspected of being regulatory genes. The point here is to maximize assay sensitivity; to be inclusive and minimize false negatives, i.e., those genes not indicated as relevant that really are, as these are difficult to recover at later stages of inquiry (witness the endurance of Gene Xs and other unknown factors in gene-network models [see, e.g., Peter and Davidson 2011]). Thus, the specificity of our interactome assays, which describes our accuracy in making connections, is less important than in assay sensitivity. The interactome, by nature, overstates the degree of interactions, thus minimizing the chance of underrepresenting the set of relevant regulatory genes.
The second consideration for applying interactome reverse-engineering algorithms is technical. Statistical algorithms require a large sample size to achieve significant correlations. Large sampling is often easier in the otherwise more intractable human disease model systems, as there are existing collections of data on expression arrays from clinical samples. Moreover, development is, as mentioned above, a dynamic process, whereas human diseases are often a contrast in two near-equilibrium states: normal and pathologic, or premetastasis and postmetastasis (Margolin et al. 2006). The comparable data set for a developmental system might be several hundred biological replicates of single embryos at, say, the onset of gastrulation. Clearly, this is not practical. Recently, it has been shown that static reverse-engineering algorithms and dynamic models are highly complementary for defining interactomes (Greenfield et al. 2010; Madar et al. 2010; Yip et al. 2010). It will therefore, be important to identify the functional set of reverse-engineering methods suitable for developmental systems. This point can be illustrated by returning to how static reverse-engineering algorithms work: They capture the statistical dependency of each gene’s expression level with that of every other gene. It is easy to see how a highly resolved expression of time course, fast enough to capture transitions, would add a new dimension in dependency measures. Despite these areas of continuing challenges, it needs to be kept in mind that the interactome is necessarily a low-resolution map, acting only as a guide to narrow the field for the next step: resolving network interactions through controlled experimentation.
The perturbation of gene function represents an essential method for determining functional gene–gene interactions. Several perturbation techniques are available: The method of choice should be as cost-efficient and as fast as possible while providing the least amount of side effects or artifacts. This could be—depending on the organism—RNAi either by injection or soaking, electroporation or feeding of siRNA or dsRNA, or soaking or injection of morpholinos. After perturbation, the expression levels of all other regulatory genes are analyzed by quantitative RT–PCR, high-throughput RT–PCR, the nanostring n-counter system, or quantitative RNA-seq. It is important to ensure a broad coverage of potential target genes during the perturbation analyses, particularly early on as these may identify additional factors missed by the interactome screen. Once a specific upregulation or downregulation of a putative target gene has been detected it may then be subjected to cis-regulatory analysis to understand whether this is the result of a direct or an indirect interaction. However, it is important to note that although perturbation is an essential part of generating these testable predictions about networks, it entails systematic errors. In particular, certain network motifs cannot be systematically resolved by perturbations alone (Greenfield et al. 2010; Madar et al. 2010). For explanatory purposes, the basic problem can be reduced to a discussion of two common and functionally important motifs (Fig. 2). Even exhaustive perturbations cannot address the problems raised by these simple, yet important, relationships. In the first case of Fig. 2 (motif 1), whether gene A acts only indirectly on gene C via gene B, or both directly and indirectly, cannot be reliably determined, despite attempts to solve this very problem (Lefebvre et al. 2010). In the second case of Fig. 2 (motif 2), B* indicates a repressor expressed in a large domain that is then repressed in a small subset of this domain by A, thus allowing C to be expressed there; because the effects of perturbing A lead to a proportionally small effect on the total expression levels of B, perturbation analysis will not detect this interaction (although spatial expression analysis could resolve this, this is not feasible for the large-scale discovery of network relations). These might seem like minor problems, but there are important consequences in misinterpreting these relationships. Each motif in isolation would be expected to perform a low-pass filter function (Alon 2007), i.e., expression of C requires sustained expression of A, either while B is turned on in the first example, or while B mRNA and protein are degraded in the second. In practice, the situation is yet more complex as these motifs are never in isolation, but instead are functionally linked with other factors, possibly in (multiple) feedback loops and via cell–cell signaling. We have shown that these additional linkages not only hinder analysis further, they also play critical roles in governing outcomes (Smith et al. 2007; Smith and Davidson 2008, 2009). Thus, subcircuits based on these problematic topologies are common and execute essential information-processing functions. Indeed, it seems likely that most regulatory switches are knotted in highly recursive, and therefore problematic, relationships (Davidson and Erwin 2006), a possibility supported on theoretical grounds (Alon 2007). It is therefore essential to employ complementary approaches to resolve important functional relationships.
The identification of transcription-factor (TF)-binding sites would provide excellent complementary data in this respect. Computational identification of TF-binding sites in a genome is uninformative; as these are typically 6–8 bp, including some degeneracy, they occur randomly in a genome every 4 kb, or even more frequently. One strategy for identifying all binding sites for a given TF is ChIP-seq. However, ChIP-seq requires specific antibodies against each TF that should be analyzed. The production of specific antibodies is often time-consuming, difficult and expensive, and therefore not practical, to say the least, for many genes, particularly in a nontraditional model system for which community resources of this nature would not be available. Furthermore, ChIP-seq requires a large amount of starting material, which can be limiting. Therefore, this method is not the first choice in many cases.
Alternatively, the combination of two different assays might be able to provide datasets that are comparable with TF-specific ChIP-seq data. The first step is to restrict the area in the genome in which to search for a TF-binding site. This can be done with FAIRE-seq (sequencing of formaldehyde assisted isolation of regulatory elements) or DNase-seq. FAIRE-seq is a method for enriching nucleosome-depleted DNA using formaldehyde fixation followed by a phenol–chloroform extraction (Giresi and Lieb 2009). The idea is to crosslink chromatin with formaldehyde, shear the DNA by sonication, and extract it by phenol–chloroform. Most genomic DNA is crosslinked to nucleosomes, which are sequestered to the interphase during the extraction. DNA that is bound by nucleosomes cannot be accessed by TFs and other proteins and is thus, “inactive”; histone crosslinking is far more efficient than TF crosslinking. Therefore, the DNA that is recovered in the aqueous phase corresponds to nucleosome-depleted regions of the genome, including transcriptional start sites, enhancers, insulators, and active promoters, covering <10% of the genome.
DNase hypersensitivity assays are based on a similar mechanism and target nucleosome-depleted genomic regions. After isolating the cell nuclei, they are mildly treated with DNaseI. As the enzyme can only be active in nucleosome-depleted regions, these are the only regions that contain the breaking points. Deep sequencing of the resulting fragments reveals the active parts of the genome.
Although the first step covered by FAIRE-seq or DNase-seq may reduce the search space for TF-binding sites, perhaps to three or four regions of 500–1000 bp each (Song et al. 2011), these methods are not TF-specific. The second step “recovers” TF specificity by cross-referencing data on the active chromatin with known TF-binding preferences, which are relatively conserved throughout animal taxa. This data set represents a large set of possible direct interactions. On the other hand, the perturbation data represents possible direct and indirect interactions. Cross-referencing perturbation data and the direct data on TF interactions should therefore, correctly identify those regulatory interactions that are direct as opposed to only indirect. Although this analysis covers problematic motifs such as described by motif 1 in Fig. 2, this still leaves the interactions that perturbation data never pick up, as illusrated by motif 2 in Fig. 2, which falls below threshold levels of detection, a potentially problematic issue for highly spatially differentiated tissue as in postgastrulation embryos. However, the first iteration does return high-quality predictions for TF-binding sites within larger cis-active elements. These elements can in effect be reprobed for TF-binding sites even in the absence of perturbation data indicating an interaction, as this represents a high-quality, yet confined search space for sites that can be readily tested by cis-regulatory validation experiments in the next phase.
The network is now more resolved in terms of directionality and accuracy than is the interactome. It is also highly resolved in terms of the specific predictions of specific cis-regulatory target sites, thanks to careful selection of the data sets that were cross-referenced along with data on perturbations. These sites can now be rapidly tested by various reporter assays because candidate sites have been predicted by the combination of perturbation, FAIRE-seq/DNase-seq, and TF-binding preference data. An approach based on Bacterial Artificial Chromosome (BAC) reporters provides a means of rigorous analysis; though other large-insert genomic libraries may be more rapidly made and may be used instead. However, a more important aspect of cis-regulatory analysis than merely validating network predictions is its potential to shed light on the dynamics of the regulation of differential gene expression. We have developed an approach for testing network function by dominant-acting cis-regulatory mutations that we call cis-reengineering (Smith et al. 2007; Smith and Davidson 2008, 2009). This is distinct from the use of inert reporters, although it often includes reporter genes as markers and usually follows cis-regulatory analyses that identify functional TF-binding sites. Instead, we employ as transgenes a BAC carrying the wild-type coding sequence for a particular gene that contains a mutated TF-binding site(s) that we hypothesize will alter subcircuit function in a predictive manner, thus determining the relevance of any given site for overall subcircuit function. The cis-reengineering approach captures both the genomic context, as these disruptions may represent just a few nucleotide changes within 150+ kb of unaltered genomic sequence in the BAC clone; and it captures the context of the endogenous gene regulatory network, as these transgenes operate in an otherwise unaltered genetic background.
We have shown in several cases how cis-reengineering can be applied to test subcircuit function during sea urchins’ embryogenesis (Smith et al. 2007; Smith and Davidson 2008, 2009). In each example, these studies demonstrated, first, the importance of network relations in driving biological outcomes; although certain kinetic parameters were restricted, the highly reliable switching behavior induced by the subcircuits depends largely on the structure of the subcircuit. Second, and more important to the discussion here, because the “testing” phase of cis-analysis generated fresh questions and new lines of inquiry, our studies showed that network analysis is best seen as an iterative approach, and scale integration as a guide for getting the most mileage out of prospective analyses rather than a rigid hierarchical structure.
The decrease in sequencing costs and the growing number of accessible nonmodel organisms render the preceding research plan increasingly feasible, putting within reach the goal of resolving the developmental GRNs of more than one species within a reasonable timeframe. Table 1 shows our criteria for what makes a system amenable to network analysis according to the procedure we outlined above. Many organisms fill these needs as no unusual demands are made; nor are large community resources required, such as those used for genetic model species. This opens up the opportunity to compare GRNs and reconstruct their ancestral configuration.
The comparison of GRNs should follow the same guidelines as any other evolutionary study. The basis is a comparative approach to first identify homologous characters shared between two or more taxa. Then, the most parsimonious hypothesis about their origin should be preferred and a comparison with representatives of the outgroup should be performed to distinguish between old and new evolutionary traits.
In the case of comparisons across species networks, genes, gene products, and their interactions should be analyzed to evaluate their homology. Very often the comparison is restricted to subcircuits, although entire developmental GRNs also can be compared.
First it must be assessed how many genes or gene products are shared between the networks in question, followed by a comparison of their interactions. Are there “enough” shared elements to justify a comparison? Are there similar genes that are part of all the compared networks? Which are the genes that occur only in some of the networks and not in others? Which interactions are repeatedly found and which are unique to certain networks?
In order to perform a meaningful evolutionary comparison, a largely overlapping set of genes should be part of all the compared networks and only some genes or interactions should vary. However, it is a difficult decision and it is hard to predict how many elements in the comparison are “enough” and how many genes and interactions should be similar. In the early days of cross-species-network-analyses, it will be partially up to the researcher to evaluate the complexity of their networks in this respect. Experience will show the minimum number of genes that should be involved in such a comparison. For now, more than five to ten genes as part of each network should be a good estimate as a lower limit.
The cross-species-network-comparison will probably show that not all the genes are included in every network under study. Those genes should be pointed out and their function should be analyzed. The next step is to compare the interactions between genes and gene products. As numerous studies have shown, we must remain open to the possibility, and some might say likelihood, that more often than the genes themselves, the interactions between genes differ. The interaction of two genes can be changed from activating to repressing, or another gene is additionally required for the activation of a target gene.
Cases in which differences regarding the genes or their connections are discovered are of particular importance to prove the independence of the compared elements. If a comparison shows that all genes of a subcircuit and all of the interactions are similar between two species, it has to be considered that these subcircuits are composed of genes that are to some extent functionally dependent on each other and thus, may not represent a complex pattern but instead be treated collectively as one character.
The next step, after evaluating and confirming the homology of two or more patterns, is the reconstruction of the ground pattern of the last common ancestor. When incorporating GRNs into the ground pattern, the genes and gene products, as well as the interactions that are shared between all networks, should be collected. To evaluate whether an element in the ancestral network is an evolutionary novelty, a comparison with outgroup representatives has to be performed. Once a ground pattern is reconstructed, other GRNs from other species can be incorporated into the comparison to test whether the hypothesized ancestral network holds true or has to be refined.
However, in contrast to traditional evolutionary studies, cross-species-network-comparison allows even more than reconstructing the putative ground-pattern GRN. It also provides testable hypotheses about which genes or interactions are the causal agents resulting in taxon-specific differences. Cis-reengineering of constructs can then be used to test whether a particular cis-regulatory element (or new gene/gene-loss) is sufficient to recapitulate the situation as found in another species or as predicted in the ground pattern, similar to previous studies (e.g., Gompel et al. 2005). In many cases, the recapitulation of an ancestral state might be technically challenging, since numerous modifications of the gene expression could be required. However, the cross-species test of regulatory elements or gene expressions can be highly informative. The possibility of evaluating evolutionary hypotheses in an experimental background thus, adds a new arena for evo–devo research.
The author was supported by National Institutes of Health (5P30GM092374-02).