Systems Analysis of Plant Functional, Transcriptional, Physical Interaction, and Metabolic Networks

GENOME-WIDE COFUNCTIONAL NETWORK INFERENCE AND ANALYSIS

Advances in high-throughput technologies have led to the generation of extensive gene expression data sets on a genome-wide scale. The observation that genes sharing a common function can be transcriptionally coordinated (Hughes et al., 2000; Kim et al., 2001; Wu et al., 2002; Stuart et al., 2003) has enabled these data sets to be analyzed such that putative functional associations between genes can be inferred. Coexpression approaches use statistical metrics to establish correlations between gene expression profiles across many samples based on the principle of guilt-by-association (Figure 1A) (Usadel et al., 2009). Given the abundance of existing gene expression data, the generation of coexpression networks is a feasible top-down approach to generate genome-wide cofunctional network models in plants.

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

Schematic Representation of Different Types of Networks.

(A) and (B) Cofunctional networks.

(A) A coexpression association where two genes are linked based on a common expression pattern.

(B) A cofunctional association where two genes are linked based on coexpression and/or other shared properties.

(C) A gene regulatory interaction where a TF directly binds the promoter to regulate the expression of a target gene.

(D) and (E) Two different approaches for uncovering GRNs using either a TF-centered approach (D) or a target-centered approach (E).

(F) Protein–protein interactions reflecting an experimentally determined physical interaction between two proteins.

(G) Metabolite interaction where two metabolites are linked through a common edge that represents a biochemical reaction converting one metabolite into the other.

Genes are colored in black, proteins are colored in gray, and metabolites are colored white.

Coexpression networks most commonly use the Pearson correlation coefficient to establish linear pairwise correlations between gene pairs in an adjacency matrix. Another associative metric that can be used is the Spearman correlation coefficient which enables nonlinear correlations between genes to be uncovered (Usadel et al., 2009). A modified graphical Gaussian model that takes into account partial correlations between genes after removing the effects of other adjacent genes has also been used (Ma et al., 2007). Following the establishment of gene associations, a cutoff threshold is then set, and pairwise interaction values exceeding this selected threshold are kept. The strength of the correlations between gene pairs can be considered as edge weights indicating the strength of coregulation between gene pairs. The end result of such an approach is a network consisting of nodes representing genes connected by edges representing a significant similarity in their common expression pattern. It is important to note that both positively and negatively acting components of a biological process can be coexpressed (Bassel et al., 2011b; Lee et al., 2011), so it is not possible to predict the function of putative coregulated gene pairs using this approach.

Functional modules consisting of subsets of highly interconnected nodes can be identified within networks. These cofunctional modules comprise strongly connected nodes, based on their collective interaction strength, with only weak interactions to other nodes in the network. In the case of coexpression networks, the modules represent groups of coexpressed genes or modular units of gene expression that may act together within a biological process with which they are associated. A variety of network clustering algorithms have been developed to identify modules in networks, including molecular complex detection (Bader and Hogue, 2003), Markov clustering algorithm (Enright et al., 2002), MCLUST (Fraley and Raftery, 2003), biclustering (Prelić et al., 2006), and the heuristic cluster chiseling algorithm (Mutwil et al., 2010). Some of these methods place genes into mutually exclusive clusters, which may not be realistic in a biological context where individual genes are capable of acting in multiple pathways, while other methods do not take into account interaction strength. A recently developed clustering algorithm named ClusterONE allows for overlap in module detection while considering edge weights (Nepusz et al., 2012). These algorithms produce clusters of varying gene members and sizes based on the subjective modification of their respective thresholds. Variability in module detection can also be compounded by the generation of different networks from the same data when using different methods (Aoki et al., 2007). Another source of variation is in the decision for placing the significance cutoff in the correlation coefficient affecting overall network size and topology. Too high a cutoff would not enable weakly coordinated processes to be captured, while too low a cutoff leads to large gene networks and a greater false positive rate. The use of known network topology to guide threshold cutoffs when generating cofunctional networks is also ambiguous, as conflicting reports have shown biological networks either to be scale free (follow a power law distribution; Barabasi and Albert, 1999) or not scale free (Khanin and Wit, 2006; Daudin et al., 2008).

Once a coexpression network has been generated, identifying modules by clustering can help extract biological meaning from the network. A module that contains genes that are annotated to overrepresented Gene Ontology (GO) categories (Ashburner et al., 2000) suggests this unit of genes may be involved in a given biological process. Uncharacterized genes within such a module can be candidates for participating in the same process. Similarly, genes directly connected to (or coexpressed with) known central regulators of a developmental process are candidates within this cofunctional framework. Identifying statistically overrepresented cis-elements in the promoters of the genes within a coexpression module may help infer regulatory interactions between transcription factors (TFs) and their targets. Given that members of a module in a coexpression network are highly transcriptionally coordinated across multiple samples, common upstream regulatory factors likely act through shared promoter elements. Shared cis-elements may be significantly overrepresented within a given module, providing clustered targets for TFs known to bind these sites. This information can be used to infer GRNs between TFs and their targets (see below).

Another means of extracting biologically relevant information from coexpression networks is to focus on genes with the greatest number of edges, so-called hub nodes, which show the greatest amount of transcriptional coordination with other genes. The centrality hypothesis states that if these connections are functionally relevant, removal of this central hub would severely perturb the network, possibly leading to an altered phenotype. The exploration of hub nodes often results in pleiotropic phenotypes due to the genes of interest being involved in multiple processes and hence can be enriched in false positives (Han et al., 2004; Yu et al., 2008a).

Coexpression networks created using various combinations of these approaches have led to the identification of previously uncharacterized genes in Arabidopsis involved in cell wall synthesis (Persson et al., 2005), flower development (Usadel et al., 2009), seed germination (Bassel et al., 2011b), and plant growth (Mutwil et al., 2010). The conservation of coexpression networks has been investigated across diverse plant species within a consensus network named PlaNet (Mutwil et al., 2011). This work demonstrated that the net transcriptional output of plant cells is at least partially conserved across model and crop species.

Integrating multiple data types can enhance the confidence by which cofunctional relationships are inferred (Figure 1B). A drawback of cofunction networks utilizing coexpression alone is that not all biological processes are transcriptionally coordinated. In addition, the guilt-by-association principle has been reported to have limited applicability in the inference of gene function (Gillis and Pavlidis, 2012). This conclusion was reached based on the observation that functional information is not encoded throughout coexpression networks, but rather within a subset of key interactions. These observations limit the general applicability of predicting cofunctional relationships based solely on coexpression analysis.

A genome-wide cofunction network of Arabidopsis called AraNet has been constructed using a number of omics data sets from Arabidopsis in addition to functional linkage data from homologous gene pairs in other model organisms (Lee et al., 2010). These data were vetted against experimentally characterized genes in Arabidopsis through standardized GO annotations available from The Arabidopsis Information Resource (Lamesch et al., 2010, 2012) using a Bayesian probabilistic model. Several previously uncharacterized genes have been assigned a function based on the functions of their neighboring genes in AraNet and subsequent experimental validation. These genes include genes associated with lateral root formation (Lee et al., 2010) and drought sensing (Lee et al., 2010).

Another cofunction network named RiceNet has been generated by combining gene coexpression data from this species with diverse proteomic data sets (Lee et al., 2011). A subset of predicted cofunctional interactions was found to interact physically at the protein level, and three novel regulators were reported to function in planta as components of the biotic stress response in this key crop species. These examples show that combining multiple data types can enhance predictive capacity of cofunctional relationships. Web-based tools to generate cofunctional networks of genes of interest using multiple data types are described later in this review.

Many studies predicting cofunctional relationships have focused on using the greatest quantity of gene expression data available, generally coming from diverse plant tissues and cell types. While more data increase the statistical power of the analysis, a lack of context specificity within data sets will only capture processes that are common in the samples used. Examples of this are the identification of photosynthesis and flavonoid synthesis pathways conserved across PlaNet (Mutwil et al., 2011). However, different tissues are defined by distinct regulatory networks, as are the discrete cell types that comprise them. A condition-dependent approach using data from a single tissue (Bassel et al., 2011b) or from single cell types (Birnbaum et al., 2003; Brady et al., 2007; Dinneny et al., 2008) can increase the context specificity by which gene cofunctional interactions can be inferred. An example of this is SeedNet, which was inferred using data generated exclusively from whole Arabidopsis seeds. A core module of seed dormancy–related coexpressed genes was identified within a dormancy subdomain of the network. Characterizing this module led to the identification of 10 novel regulators controlling seed dormancy (Bassel et al., 2011b). These predicted gene interactions were absent in other condition-independent cofunctional networks, while the overall network topology was found to be conserved in seeds of the monocot crop species wheat (Triticum aestivum) when equivalent data sets were compared.

When using condition-independent data sets, strategies have been developed to uncover biologically significant cofunctional interactions. One approach to enrich for causal interactions is to remove from the network all edges that are conditionally dependent (i.e., those that cannot be explained by other variables). Such a selection of conditionally independent edges has been used successfully for small metabolic and developmental networks by Gaussian Graphical Modeling (Wille et al., 2004). However, computational complexity prevents genome-wide applications at this time.

The spatial and temporal modulation of gene networks during cellular development may be limited by the availability of a sufficient number of data sets capturing cell type–specific processes over time. AraNet is capable of detecting genes preferentially coexpressed within specific root cell types with greater accuracy than randomly generated networks, despite using diverse whole tissue data sets as input. This predictive capacity may be increased through the generation of condition-dependent cofunctional networks.

A novel computational approach to construct a cofunctional network using gene expression data has recently been presented (Bassel et al., 2011a), wherein a rule-based machine learning algorithm called BioHEL was used to associate genes functionally. This algorithm generates rules that predict the developmental fate of a biological system through the automated exploration of subsets of annotated gene expression data. This cofunctional association metric, termed “coprediction,” relies on the hypothesis that groups of genes that collectively are able to predict a developmental output have a greater probability of being functionally related. This approach was used to generate a seed coprediction network (SCoPNet), which was functionally validated with the identification of four novel regulators of seed germination. A key feature of coprediction using gene expression data is that it does not rely on establishing cofunctional relationships using common gene expression patterns. Additional computational approaches beyond correlation-based methods may therefore be used to associate genes functionally using large-scale datasets.

LARGE-SCALE TRANSCRIPTIONAL REGULATORY NETWORKS

Mapping and systems analysis of large scale regulatory networks allow for the discovery of emergent properties (Long et al., 2008). A GRN describes the interaction of regulatory factors and their targets; here, we focus on TFs and their regulation of target gene expression (Figure 1C). There are two complementary ways to map GRNs: TF centered and target centered (Figures 1D and 1E).

Chromatin immunoprecipitation (ChIP) coupled with microarrays (ChIP-chip) or next generation sequencing (ChIP-seq) is a TF-centered approach to map GRNs. This technique involves immunoprecipitation of a TF and its bound target DNA with subsequent identification of downstream targets. The ChIP approach results in a network that describes the in vivo binding targets of a single TF (Table 1). Yeast one-hybrid (Y1H) assays are a target-centered approach that identifies upstream TFs binding to a promoter or motif of interest. This assay requires the initial development of a TF collection or cDNA library of prey TFs, which are assayed against bait promoters of interest in yeast (Deplancke et al., 2004).

ChIP studies have focused on TFs that are well characterized and typically abundant. In plants, ChIP studies have been performed primarily in Arabidopsis, but occasionally also in rice (Oryza sativa) and maize (Zea mays) (Fornalé et al., 2010; Moreno-Risueno et al., 2010) (Table 1). A protein expressed in few cell types or developmental stages can be difficult to precipitate due to lack of abundance. Traditionally this has been overcome using overexpression and conditional induction of the gene of interest by fusion to the glucocorticoid receptor or to a hormone inducible promoter, such as the estrogen-inducible promoter. Although this provides the researcher with potential binding targets, it obscures the developmental or tissue/cell type–specific context of the transcriptional regulation identified. Additionally, high-quality antibodies against the TFs are needed to achieve specific and sufficient immunoprecipitation. Using epitope-tagged TFs or fusion with a protein to which high-quality antibodies have been developed (i.e., green fluorescent protein) has been a way to circumvent the development of protein-specific antibodies. However, the tagged protein must be verified to ensure that its expression reflects the wild type and that it can functionally complement a loss-of-function phenotype.

ChIP studies provide evidence for in vivo DNA binding but do not show the regulatory effect of the TF on its targets. Similarly, Y1H does not show the regulatory effect of the binding. Additionally, Y1H is limited by the heterologous nature of the assay and does not provide direct in planta evidence for TF-DNA binding. For example, a plant TF may not be able to form obligate heterodimers or higher order protein complexes to bind to promoter sequences. This will produce false negatives that cannot be detected in the yeast system. Therefore, this approach requires biological validation to be confident of the in vivo physical interactions through independent validation techniques, including directed ChIP experiments. However, Y1H can be adapted to high-throughput techniques to screen the promoters of multiple genes of interest simultaneously, leading to a typically larger and more comprehensive network (Brady et al., 2011; Gaudinier et al., 2011). The assay can be tailored to screen specific biologically relevant promoter sequences of interest, focusing the experimental approach and enriching for the identification of relevant targets when mapping GRNs.

Although previously only 1 to 10% of interactions identified by ChIP have been determined to be regulatory in nature (Farnham, 2009; Moreno-Risueno et al., 2010), more recently biological validation rates have been increasing (Kaufmann et al., 2011). Using Y1H assays to determine TF-DNA interactions, two of three interactions tested were reported to be regulatory in planta by ChIP and genetic analysis. Thirty-two of 66 additional interactions were confirmed through genetic approaches with the supposition that the same proportion above (two out of three) would bind in planta (Brady et al., 2011). Comprehensive biological validation of GRNs predicted by either of the above methods is a challenge. The ChIP approach provides in planta binding data and the Y1H assays show binding, at least in yeast, but neither method is sufficient to probe the regulatory nature of the interactions. Perturbing the network using loss-of-function and/or conditional induction lines of the TF (such as using glucocorticoid receptor fusions or estrogen-inducible promoters) and analyzing expression of the TF’s targets relative to the control will show if the TF acts as an activator, repressor, or has no effect. In vivo confirmation can also be performed using a transient expression system in protoplasts in which a TF is induced to see its effect on target gene expression (Pruneda-Paz et al., 2009; Lau et al., 2011), although the developmental and tissue-specific context of the interactions can be missed using this approach.

Both ChIP and Y1H led to genome-wide understanding of transcriptional regulatory networks. For example, ChIP-chip and ChIP-seq technologies have been used to identify binding targets of the general floral regulator SEPALLATA3 (SEP3) (Kaufmann et al., 2009) and of the A-function floral identity regulator APETALA1 (AP1). As many as 2298 genes were found as AP1 targets, with 10.8% of these genes as high-confidence targets that show binding and a moderate change in gene regulation by AP1 as determined by induction of AP1 and ChIP-seq analysis. AP1 bound to the SEP3 promoter and activated SEP3 expression. A comparison of AP1 and SEP3 targets (Kaufmann et al., 2010) showed overlap, suggesting that these proteins may function as a heterodimer to activate genes involved in early flower organ development. Using a combination of ChIP and network predictions, an opposing relationship between two proteins involved in the response to iron deprivation, POPEYE (PYE), a basic helix-loop-helix TF, and BRUTUS (BTS), a putative E3 ligase, has been described (Long et al., 2010). Using ChIP-chip of PYE, 70 targets were identified. Using the ATTED-II coexpression network tool with PYE as bait, BTS, which is involved in iron homeostasis, was identified. The pye-1 mutant plants are small and chlorotic when grown on iron-deficient media, while BTS plays an antagonistic role in negatively regulating the iron deprivation response as bts-1 plants are green and have an increased tolerance to iron deficiency. Y1H assays have been used at various scales to identify upstream transcriptional regulators. A directed Y1H screen, with a library of 186 Arabidopsis TFs that are circadian regulated, against the CCA1 promoter was used to identify CCA1 HIKING EXPEDITION (CHE), a TCP TF. Using tiled promoter fragments, the binding region was narrowed down to 171 bp (Pruneda-Paz et al., 2009). A transient expression assay in protoplasts determined that CHE acts as a repressor of CCA1, and in planta binding was confirmed by ChIP–quantitative PCR (qPCR). Yeast two-hybrid (Y2H) and coimmunopurification experiments additionally demonstrated a protein–protein interaction between CHE and TOC1 as well as between TOC1 and CCA1, indicating further putative combinatorial action between these TFs (Pruneda-Paz et al., 2009).

Finally, combining approaches will lead to more comprehensive networks. For example, the combination of Y1H and Y2H assays between 167 TFs and 93 TF and microRNA promoters as well as computationally and experimentally derived microRNA–mRNA interactions resulted in a more holistic Arabidopsis root stele–enriched GRN (Brady et al., 2011). A network of 103 interactions was mapped and interactions were tested in planta with a variety of biological validation techniques, including ChIP, mutant and conditional induction analyses to confirm 59% of the interactions. Only 16% of the GRN TF mutants screened showed a root phenotype but 65% showed a molecular phenotype. This demonstrates the robustness of the stele-enriched network as TFs regulating the same gene can compensate for a mutation in one of the regulators. Increasing this TF collection from 24% to 92% of stele-expressed TFs (Gaudinier et al., 2011) provides the ability to generate a more comprehensive network of transcriptional regulation across the root stele.

LARGE-SCALE PROTEIN–PROTEIN INTERACTIONS: THE INTERACTOME

To understand fully genotype-to-phenotype relationships at the systems level, comprehensive knowledge about the complex and dynamic protein–protein interactions that exist within an organism is required. Identifying binary protein–protein interactions and components of protein complexes is a crucial step toward elucidating associated biological activities (Figure 1F). Several assays have been developed to identify and predict protein interactions in yeast, in cell suspension, in planta, and in silico, and various public databases have compiled these interactions (Table 2).

As there is still a long way to full coverage of in planta protein–protein interaction networks, computational approaches, including coexpression, have been used to infer protein–interactions. Protein interactions can be inferred by homology to known interactions in other organisms (interaction ortholog or interlog) (Geisler-Lee et al., 2007), using indirect evidence or literature (Cui et al., 2008) or integrated methods (Xu et al., 2010). Interlogs can be filtered using functional association data to improve prediction reliability (De Bodt et al., 2009). Several databases of inferred protein interactions have been established (Geisler-Lee et al., 2007; Cui et al., 2008; Brandão et al., 2009; De Bodt et al., 2009; Lin et al., 2009, 2011a, 2011b) (Table 2). The predicted Arabidopsis interactome resource provides an interactome inferred from multiple indirect lines of evidence, including coexpression, colocalization, coevolution, annotation similarity, domain interaction, and homologous interactions in other species. The BAR Arabidopsis Interactions Viewer also provides similar functionality (Geisler-Lee et al., 2007). These computational approaches can guide the identification of potential in vivo and in planta protein–protein interactions. In this respect, coexpression is relevant because proteins generally should be expressed in the same cell at the same time for in vivo interaction. Care must be used when inferring protein expression from transcriptome data as transcript and protein abundance do not always correlate (Boruc et al., 2010; Van Leene et al., 2010; Arabidopsis Interactome Mapping Consortium, 2011; Petricka et al., 2012).

While inferring interactions using computational approaches can be useful, validation in a biological system is required to draw more meaningful conclusions. This is limited by only a few available in planta high-throughput assays. To gain insight into protein–protein interactions and to obtain interactome networks, one can identify binary interactions using Y2H, bimolecular fluorescence complementation (BiFC), split ubiquitin or split luciferase, or alternatively isolate components of a protein complex in planta using coimmunopurification or tandem affinity purification (TAP) (Van Leene et al., 2011). Uses of these different approaches are described below.

On a small, protein family scale, a Y2H approach was used to identify binary protein–protein interactions, such as AUXIN/INDOLE-3-ACETIC ACID (AUX/IAA)–AUXIN RESPONSE FACTOR (ARF) interactions (De Rybel et al., 2010; Vernoux et al., 2011), to map the dimerization network of MADS domain TFs (de Folter et al., 2005) and, in combination with BiFC assays, to reveal functional modules in core cell cycle binary protein–protein interactions (Boruc et al., 2010). To obtain a more global interactome network, a high-throughput Y2H approach has been used to map binary interactions between soluble proteins in plants (Arabidopsis Interactome Mapping Consortium, 2011). This represents the first large-scale experimentally validated protein interactome in plants and has also been used to determine a plant-pathogen immune network (Mukhtar et al., 2011).

However, Y2H is mainly suitable for soluble proteins that can traffic to the nucleus, ignoring membrane proteins, which often play essential roles in fundamental biological processes, including signaling, homeostasis, nutrient acquisition, and metabolism. Despite their importance, we know little about the functions of most membrane proteins. To systematically elucidate a map of membrane protein interactions, a mating-based split ubiquitin system (mbSUS) was developed (Obrdlik et al., 2004; Miller et al., 2005). The split ubiquitin system is similar to the classical Y2H as it uses yeast as a heterologous system and has a similar readout, but it specifically allows the detection of interactions of full-length membrane proteins. The concept of mbSUS relies on the release of a TF from a membrane protein if it interacts with another membrane (or soluble) protein. Similar to other detection systems based on the reconstitution of two halves of a protein, mbSUS uses an ubiquitin molecule split into two halves. The use of a mutated N-terminal domain of ubiquitin (Nub) and C-terminal half (Cub) fused to a TF to test for physical interaction between two fused proteins is the basis of mbSUS, whereby a functional ubiquitin is reconstituted only when brought into vicinity via interaction between two fusion partners (Johnsson and Varshavsky, 1994). Cleaving the reunited ubiquitin molecule by an endogenous ubiquitin-specific protease releases the TF to move into the nucleus and activate reporter genes. The mbSUS was used to test for interactions of the translocon complex at the outer chloroplast membrane (Rahim et al., 2009) and to screen for potential interactions among 490 Arabidopsis membrane and signaling proteins (Lalonde et al., 2010). However, the currently available membrane interaction networks for Arabidopsis cover only a small portion of the membrane-bound proteome. A recent study used a computational classification system to predict putative membrane protein interactions and decrease the sample space of interactions to be proved (Chen et al., 2012). These predictions were validated through the verification of 541 interactions between 239 membrane proteins, enriched in transporters.

The above-described yeast-based systems are informative, though high-throughput proteome-wide in planta interaction data are more biologically relevant. An alternative approach is TAP, which allows investigation of protein interactions through in situ purification of protein complexes. In a purification experiment, complexes that are associated with a protein of interest at the time of extraction can be isolated. The TAP approach is based on the expression of a bait protein fused to a double affinity tag and a two-step purification process followed by mass spectrometry of protein complex members (Van Leene et al., 2008). A targeted mapping of the cell cycle interactome, using TAP, provided a first draft of the basic cell cycle complex machinery in Arabidopsis (Van Leene et al., 2010). The TAP approach has also been employed to isolate the core jasmonate signaling module and to find new interactors of JASMONATE ZIM-DOMAIN repressor proteins (Pauwels et al., 2010; Fernández-Calvo et al., 2011) to determine whether PROHIBITINS are present within a multimeric complex (Van Aken et al., 2007), to purify 14-3-3 complexes (Chang et al., 2009), and to identify components of complexes, including the COP9 signalosome 3 complex (Rubio et al., 2005), novel WPP2-interacting proteins (Zhao et al., 2008), the BRASSINAZOLE-RESISTANT1–dephosphorylating phosphatase (Tang et al., 2011), and E2F TARGET GENE1–associated proteins (Takahashi et al., 2008). However, current studies using TAP rely on the use of cell suspension cultures. Implementing this technology in planta within the context of multicellular tissues will provide further biological insights. Furthermore, improving the methodology to identify membrane-associated complexes awaits development.

For the cell cycle interactome study, interaction between 17 protein pairs was validated using split luciferase assays. Here, the firefly luciferase (LUC) protein is split into two halves and fused to two different proteins. LUC activity is only reconstituted when the N- and C-terminal LUC moieties are brought together by the two interacting proteins. The validation rate was 41%, which is a great improvement over a similar type of split protein reconstitution performed following TAP tagging in yeast (Yu et al., 2008a; Van Leene et al., 2010). The actual false positive rate is expected to be much smaller because in a TAP tagging experiment proteins may interact with each other indirectly and therefore may not be confirmed when validated by other assays. Furthermore, these interactions likely occur in specific cell types that may require additional temporal cues not present in transient systems. The LUC-based approach was also used to detect SNARE–SNARE interactions in Arabidopsis protoplasts (Kato et al., 2010). BiFC, which uses a similar strategy of bringing together two parts of a protein, in this case yellow fluorescent protein, was also recently used to validate 78 G-protein interactions in planta that were previously identified by Y2H analysis. Here, the validation rate was a striking 95%, which is the highest validation rate reported for any organism (Klopffleisch et al., 2011). However, this bait set of proteins was biased in their original selection by their common participation in a biological process, which may partially explain the high rate of validation relative to other studies. In addition, the high sensitivity of the stable reconstituted fluorescent protein can lead to unspecific reconstitution, especially upon overexpression of the fusion proteins of interest (Lalonde et al., 2008).

Another approach to probe protein interactions over a large scale is to use protein microarrays (Popescu et al., 2007, 2009). This involves the generation of recombinant proteins in planta through transient expression and spotting the purified proteins onto microchips. These chips may then be hybridized with fluorescently tagged proteins of interest or subjected to enzymatic modification with a labeled substrate. Imaging of the arrays following incubation can identify protein–protein interactions and substrates of posttranslational modification. Using this approach, large-scale calmodulin interaction (Popescu et al., 2007) and mitogen-activated protein kinase phosphorylation networks have been indentified (Popescu et al., 2009).

In addition to these approaches, direct in planta immunoprecipitation of tagged proteins of interest, followed by mass spectrometry of coimmunoprecipitated proteins has been successful in identifying protein complex members for several individual proteins in Arabidopsis, such as SOMATIC EMBRYOGENESIS RECEPTOR-LIKE KINASE1 (Karlova et al., 2006), CELLULOSE SYNTHASE CATALYTIC SUBUNITs (Desprez et al., 2007), MADS domain TFs (Smaczniak et al., 2012), and the ESCRT-III component VPS2.2 (Ibl et al., 2012). While this is a powerful approach, uncovering rare, dynamic complexes in specific plant cells and tissues remains challenging.

A combination of the above approaches led to the description of a number of interactomes in Arabidopsis, not only for the cell cycle (Boruc et al., 2010; Van Leene et al., 2010, 2011), but also for the chloroplast (Yu et al., 2008b) and for ESCRT, TOPLESS (TPL), TPL-related (TPR), and G-proteins (Klopffleisch et al., 2011; Shahriari et al., 2011; Causier et al., 2012; Ibl et al., 2012). For example, the TPL and TPR interactomes revealed that TPL proteins have been co-opted multiple times in evolution to cause transcriptional repression. In any systems-level analysis of protein–protein interactions, in planta validation of the entire spectrum of interactions that were observed or predicted is unlikely, but the tools above can be used to confirm a subset of these interactions. In addition, the cell cycle interactome showed that various techniques may be complementary as the overlap was rather small between Y2H and BiFC (Boruc et al., 2010).

Genetic analysis is an additional method of in planta validation. In an experiment mapping transcriptional regulatory interactions in the Arabidopsis stele, physically interacting TFs were identified using Y2H. TFs that interact with each other previously have been shown to influence each other’s expression, either directly or indirectly (Cui et al., 2007). For each TF of an interacting pair, expression of the interactor was measured in each respective TF loss-of-function background. The majority of interacting TFs was found to have a regulatory influence on each other’s expression (Brady et al., 2011). In the Y2H interactome study to determine a plant-pathogen immune network, interaction between Pseudomonas syringae and Hyaloperonspora arabidopsidis effectors with Arabidopsis proteins was determined (Mukhtar et al., 2011). A prediction was that Arabidopsis proteins that were targeted by effectors from both pathogens could display pathogen resistance or susceptibility. Fifteen of the 17 proteins targeted by effectors from both pathogens showed mutant phenotypes consistent with immune system functions (Mukhtar et al., 2011).

METABOLIC NETWORKS

Genome-scale metabolic network models have been built and applied to study a wide variety of topics, including adaptive evolution, network properties, and metabolic engineering (Figure 1G) (Feist and Palsson, 2008). Genome-scale metabolic models in metabolic engineering have been used to couple growth rate with overproduction of lactate (Fong et al., 2005) and l-Val (Park et al., 2011) in Escherichia coli and prioritization of alternative engineering strategies for overproducing ethanol in yeast (Bro et al., 2006). Although most of the work in metabolic network reconstruction and application has been performed in microbes, several advances in plant metabolic network reconstruction and modeling have been made recently, paving the way for rational plant metabolic engineering.

Two types of plant metabolic networks have been developed: descriptive and predictive. Genome-wide descriptive metabolic networks have been predicted from metabolic pathway databases that have been curated from experimental data in the literature (Caspi et al., 2012; Kanehisa et al., 2012). These genome-wide metabolic network representations are now available for several plant species, such as Arabidopsis (Mueller et al., 2003; Zhang et al., 2010), Populus trichocarpa (Zhang et al., 2010), Chlamydomonas reinhardtii (May et al., 2009), Medicago truncatula (Urbanczyk-Wochniak and Sumner, 2007), grasses (Youens-Clark et al., 2011), and shade plants (Bombarely et al., 2011). Recently, the plant metabolic network group published genome-wide metabolic networks for maize, soybean (Glycine max), cassava (Manihot esculenta), and wine grape (Vitis vinifera) online (www.plantcyc.org; Zhang et al., 2010). Most of these networks (except for Arabidopsis) are not associated with extensive validation and curation with experimental data from the literature, and caution should be taken when interpreting the data. These networks have been used for a variety of studies, including identification of novel genes in glucosinolate metabolism (Chan et al., 2011), the role of duplicated enzymes (Hanada et al., 2011), and metabolomic responses to UV-B (Kusano et al., 2011). While the descriptive, qualitative networks have been useful in many studies, they are not predictive and are therefore limited in facilitating metabolic engineering.

Predictive metabolic modeling approaches can be broadly grouped into kinetic and stoichiometric modeling (Sweetlove et al., 2008). Kinetic modeling uses enzyme kinetics numerically to simulate and test metabolic fluxes and can explain mechanism of flux changes. However, the difficulty of determining in vivo enzyme kinetics has limited this modeling to a small number of pathways. The most widely adapted modeling approach used for genome-wide metabolism is constraints-based modeling that includes stoichiometric, thermodynamic, and flux capacity constraints to model the fluxes of metabolites (Thiele and Palsson, 2010). This approach has been used to build predictive models of metabolism for Arabidopsis (Poolman et al., 2009; de Oliveira Dal’Molin et al., 2010), maize (Saha et al., 2011), and C. reinhardtii (Chang et al., 2011). Most of these models have not been validated extensively using flux measurements, though advances in metabolic flux analysis using ¹³C-labeling and metabolomics approaches hold promise (Schwender, 2008; Sweetlove et al., 2008; Allen et al., 2009).These predictive models have been applied in studies in many fields, such as metabolic engineering, drug discovery, drug target discovery, identification of novel gene function, evolutionary processes, network behaviors, and interpretations of mutant phenotypes (Feist and Palsson, 2008). The most common algorithm used in such studies has been flux balance analysis (FBA), which attempts to balance the stoichiometry of the metabolites within the metabolic network system with a goal (objective function) of maximal growth or maximal biomass accumulation. While prediction of fluxes using FBA matches well with experimental data (Burgard and Maranas, 2003), its assumptions may not always hold true, especially for engineered mutant lines. Several algorithms that have the goal of minimizing the change in the metabolic network upon perturbation have been developed and appear to perform better than FBA in explaining fluxes of mutants (Segrè et al., 2002; Shlomi et al., 2005; Herrgård et al., 2006).

The field of metabolic network analysis in plants is still young, and there are several challenges ahead. First, the process of generating high-quality networks is still time consuming and requires much manual input. More advances in automating these steps, such as the Model SEED system, are needed (Henry et al., 2010). Second, compartmentalization at the subcellular and tissue level must be considered for more accurate modeling of plant metabolic networks. Recent advances, such as the large-scale C4 photosynthesis model that incorporates mesophyll and bundle sheath–specific reactions (Dal’Molin et al., 2010) and Arabidopsis tissue-specific network models (Mintz-Oron et al., 2012), are early examples of the next-generation metabolic network models. Third, metabolic network models need to be integrated with regulatory and signaling models to allow more accurate and comprehensive understanding of plant metabolism. For example, large-scale experimental validations showed that the metabolic network in E. coli can explain fluxes in ∼50% of the conditions (Covert et al., 2004). Adding transcriptional regulatory constraints increased predictability to ∼70%. Transcriptional regulatory constraints have been incorporated into FBA analysis using the steady state regulatory FBA approach (Shlomi et al., 2007). More recently, an attempt has been made to develop a system to integrate signaling, metabolic, and regulatory networks into a common framework (Lee et al., 2008). Finally, metabolomics (comprehensive measurement of small molecules) in plants promises to be a useful source for validating the structure of, identifying missing links in, and understanding the dynamics and evolution of, metabolic networks (Fiehn et al., 2000; Kliebenstein, 2009; Nakabayashi et al., 2010; Tohge et al., 2011). However, several issues remain before these data can be integrated into metabolic network analysis, the biggest of which may be the high portion of unnamed metabolites in these studies. In addition, the mechanistic relationships among the levels of metabolites, transcripts, proteins, and enzyme activities are as yet unclear. An elegant study examining diurnal cycles of Arabidopsis suggests that metabolites regulate transcripts and not vice versa (Gibon et al., 2006). Despite these challenges, the advances in the omics and visualization technologies will help move this field forward in both expected and unexpected ways.

TOOLS AND TOYS FOR NETWORK ANALYSES

Any systems biology study depends on suitable software and data access, and Table 2 provides a selection of some useful tools. Much can be learned from available information, and online databases and data warehouses provide access to this body of prior knowledge. Repository databases focus on data storage and retrieval, while data warehouses contain extensively curated data and often include analysis tools. Many of the popular data warehouses in plant biology provide access to transcriptome data. Tools such as Genevestigator, the Bio-Array Resource for Plant Biology, and ATTED II are useful to search for genes or proteins with particular expression profiles or coexpression with bait genes (Toufighi et al., 2005; Hruz et al., 2008; Obayashi et al., 2011). Functional modules in biological networks are often composed of coexpressed and colocalized proteins that are highly connected in protein interaction graphs. A large set of databases is dedicated to protein–protein interaction data, some of which include information on coexpression. Several initiatives, such as IntAct (Kerrien et al., 2006), the Molecular Interaction database (Chatr-aryamontri et al., 2007), the Database of Interacting Proteins (Salwinski et al., 2004), the Biomolecular Interaction Network Database (Bader et al., 2001, 2003; Alfarano et al., 2005), and BioGRID (Stark et al., 2006), have been established to systematically collect and organize the interaction data reported by proteome-scale high throughput experiments and low-throughput studies focusing on individual proteins and pathways. In this context, the Arabidopsis network analysis pipeline is a useful tool because it integrates 11 Arabidopsis protein interaction databases (Wang et al., 2012). One of the most comprehensive resources for protein localization is the subcellular location database for Arabidopsis proteins, which collects both experimental and predicted subcellular protein localizations in Arabidopsis (Heazlewood et al., 2007).

Gene set analysis helps characterizing user-defined gene lists. Gene set analysis searches for distinguishing properties of gene lists, for instance, the overrepresentation of certain GO categories or functional bins for metabolic and regulatory processes, such as those defined in the MapMan software (Usadel et al., 2005). Gene set analysis should be supported by statistical evidence that frequencies of called categories are different from what is expected by chance. The hypergeometric test is an appropriate statistical framework to evaluate overrepresentation in small gene lists. Note that testing of multiple categories requires control of type I (false positive) error rates (e.g., the Bonferroni correction) (Dudoit and VanderLaan, 2008). The reference population used to construct the test’s null hypothesis must be chosen carefully to include only elements from the sampled space and not from the entire genome. Gene set analysis together with cofunction networks can be used to infer gene function. Gene associations solely based on coexpression are not always informative, but predictive power increases when associations include other type of data, such as cocitation or protein interactions. AraNet, GeneMania, and STRING are three such tools that are user friendly (Warde-Farley et al., 2010; Hwang et al., 2011; Szklarczyk et al., 2011).

Traditionally, systems biology often focuses on single-model species. Nevertheless, much can be learned by exploiting comparative genomics tools, which can help assign functions to multigene family members by defining ortholog sets, identify relevant gene regulatory motifs, or study how evolution or breeding shaped biological systems. The Plaza database, for instance, makes genomic data from 25 plant species available (Van Bel et al., 2012). Deciphering transcriptional regulation involves identification of both cis-regulatory elements, for instance, by motif samplers, and of regulatory trans-acting TFs. Several algorithms and tools are available for such analyses (Ladunga, 2010), but the currently fragmented knowledge about plant TF binding sites makes reverse engineering of promoter functions a nontrivial task.

Data integration and network inference are two major needs in systems biology. LeMoNe and ENIGMA, for instance, are software packages to extract module networks from gene expression data (Michoel et al., 2007; Maere et al., 2008). Coexpression-based correlation networks can also be constructed with TMEV (Saeed et al., 2006). Similarly, CORNET constructs coexpression networks but also provides workflows to integrate such networks with protein–protein interaction and localization data (De Bodt et al., 2010) (Figure 2). Network inference is an active field in machine learning and statistics communities, and powerful new algorithms are continuously being developed. To unleash the full potential of available algorithms, software environments with more open access are needed. Here, the general computation environments R (www.R-project.org/) and MatLab (www.mathworks.co.uk/products/matlab/index.html), with their large package libraries for specialized tasks, and the Systems Biology Workbench, with its rich model library (Sauro et al., 2003), are powerful alternatives. These environments lack much of the convenient graphical user interfaces characteristic of most of the previously mentioned tools, but their inherent scripting design standardizes reproducing and transferring analysis workflows within or between laboratories. In particular, the Systems Biology Workbench offers great modeling options for biological systems with minimal scripting needs and is well suited to analyze biochemical reaction models (Sauro et al., 2003).

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

Network Examples for Protein–Protein and Regulatory Transcriptional Interactions.

Examples of networks starting from a few lateral root initiation components with proven protein–DNA and protein–protein interactions. We selected ARF7/NPH4 (AT5G20730), ARF19 (AT1G19220), SLR (AT4G14550), LBD18 (AT2G45420), LBD33 (AT5G06080), E2FA (AT2G36010), and DPA (AT5G02470) to illustrate the complexity of the networks at various levels. A coexpression network was tested for protein–protein interactions (A) and for regulatory interactions (B). We used standard settings (0.6 to 1.0) for the Pearson’s correlation coefficient, and all the available data sets in the CORNET database (https://cornet.psb.ugent.be/) (with at least meeting the conditions in three data sets) were used.

For identification of protein–DNA interactions, ChIP-chip and ChIP-seq data analysis is a multistep procedure. In the case of ChIP-seq, reads first need to be aligned to a reference genome with a read-aligner, such as bowtie (Langmead et al., 2009), to generate genome coverage data. Then, peaks must be detected from coverage data or microarray signals. Finally, results should be visualized as an important validation step to test whether parameters and program settings were adequate. The experienced user will perform many of these tasks in R using packages such as CSAR (Muiño et al., 2011b). Others can choose from a number of software suites and analysis pipelines, such as the stand-alone program CisGenome (Ji et al., 2008), the Web tool PRI-CAT (Muiño et al., 2011a), or the commercial CLC Genomics Workbench (www.clcbio.com). More thorough descriptions of procedures to analyze ChIP-chip and ChIP-seq data can be found elsewhere (Pepke et al., 2009; Wilbanks and Facciotti, 2010). For visualization, various genome browsers are available. We found the integrated genome browser (Nicol et al., 2009) user friendly and robust. Like other such browsers, it allows sharing of results via network access. Researchers need only to upload their results to a server and make the URL publicly available to allow the community to browse through their findings.