BACKGROUND
biologists have long been fascinated by the green plant leaf and have tried to understand how leaves are born, live and die.
in the last decades, several new approaches to study the structure and function of leaves have emerged: molecular biology and molecular genetics have, for example, enabled identification of genes that regulate the primary function of the leaf - photosynthesis - and leaf development has been understood in much greater detail; high through-put transcriptomics has identified additional factors influencing leaf function, but traditional transcriptome analyses typically reduces the problem of finding key regulators to detecting differentially expressed genes or computing pair-wise similarity between targets and putative regulators .
in contrast, systems biology analysis of transcriptional programs treats genes as interacting rather than isolated entities.
thus these methods can begin to understand how so-called emergent properties such as complex phenotypes arise from interacting genes.
whether this can be seen as taking a holistic rather than a reductionistic approach to science has generated quite some debate  <cit> , but systems biology methods account for synergistic and competitive effects between regulators that individually could have low similarity to the target.
methods for reverseengineering the transcriptional network from collections of gene expression data have been pioneered on single-cell organisms, but have increasingly been applied to higher order organisms  <cit>  including plants  <cit>  where applications of systems biology methods are now emerging.
most systems biology studies have - not surprisingly - utilized using "the model plant" arabidopsis thaliana, where large transcriptomics programs have generated adequate quantities of high-quality data to enable systems analysis  <cit> .
for example, carerra et al.
<cit>  modeled the transcriptional network of arabidopsis and identified plant-specific properties such as high connectivity between genes involved in response and adaptation to changing environments.
however, not all aspects of plant biology can be studied in arabidopsis, which in many respects is a rather atypical plant.
indeed, it was not selected as a model system due to its physiological and ecological qualities, but rather for its suitability for genetic and genomic studies.
therefore, it is important to perform parallel studies in plants with other characteristics, as well as developing the methods to allow data from the arabidopsis system to inform studies in other organisms.
one rapidly emerging plant model system is populus  <cit> ; it's interesting biology  and the access to a sequenced genome  <cit>  represent an attractive combination.
correspondingly, more advanced data analyses approaches are now being applied in populus.
populus provides an attractive model system for studies of leaf biology.
for example, sjödin et al.
<cit>  exploited the fact that mature aspen  in boreal regions have the rather unique property that all leaves emerge simultaneously from overwintering buds.
this provides a synchronized system, resulting in a full temporal separation of the leaf developmental stages and subsequent acclimation that could be exploited using transcriptomics.
access to a centralized repository of much of the populus cdna microarray data  <cit>  and databases for the analysis of gene expression - and other - data  <cit>  substantially facilitates the ability to perform systems biology studies.
for example, grönlund et al.
<cit>  induced a co-expression network revealing modular architecture explaining gene function and tissue-specific expression; street et al.
<cit>  identified co-expression networks across a large collection of leaf transcriptomics data and found that some network hubs have existing functional evidence in arabidopsis; quesada et al.
<cit>  performed a comparative analysis of the transcriptomes of populus and arabidopsis, and found evidence of extensive remodeling of the transcriptional network, although some essential functions showed little divergence.
a few studies have also integrated promoter information to study regulatory control in populus.
shi et al.
<cit>  identified combinations of xylem-specific motifs in populus promoters.
another study inferred transcriptional networks in xylem, leaves, and roots, and showed that genes with conserved regulation across tissues are primarily cis-regulated, while genes with tissue-specific regulation are often trans-regulated  <cit> .
all these studies are essentially co-expression networks that visualize expression similarity between pairs of genes, but do not infer complex interactions.
network inference methods using expression data can be divided into those that aim to model the general influence that genes have on the expression of other genes   <cit>  and methods that aim to model the physical interaction between transcription factors and the regulated genes   <cit> .
both approaches employ common network inference methods , but those that infer gene regulatory networks also typically integrate motif finding and detection of transcriptional modules  <cit> .
approaches that describe how the regulatory genome orchestrates dynamic gene expression has developed from pilpel et al.
<cit> , who showed that yeast genes sharing pairs of binding sites in their promoters were significantly more likely to be co-expressed than genes sharing only single binding sites, to various machine learning methods that identify modules of co-expressed genes with common motif patterns in their promoters  <cit> .
here we apply a network inference method combining promoter information and expression data to describe the transcriptional network in populus leaves.
our aims were  to detect regulatory hubs in leaves,  to describe conservation of transcriptional regulation within populus and between populus and arabidopsis, and  to understand the regulatory complexity in leaves by comparing systems biology and traditional bioinformatics as methods for detecting target genes for further analysis.
this study goes beyond previous meta-analyses of populus transcriptome data by taking into account synergistic and competitive interactions between regulators, and by systematically integrating the regulatory genome and the transcriptome to infer networks.
we show that our network is robust, explains available gene function information and generalizes to new expression data in both populus and arabidopsis.
we identify the main regulators of primary processes in leaves, and show how some of these have regulatory partners orchestrating expression either in a synergistic or competitive manner.
such interactions are not considered by pair-wise similarity methods, and thus several of the regulators predicted here would not have been identified by traditional approaches.
RESULTS
we inferred the regulatory network of a collection of  <dig> leaf-specific populus genes with quantified transcription profiles across  <dig> samples in various experiments such as leaf primordial, budset, biotic infection and drought stress  <cit>  .
the approach employed two separate steps to construct the network : first, we discovered a set of representative transcriptional modules containing co-expressed genes with evidence of co-regulation in their promoters.
second, we inferred the most likely regulators  of each module based on gene expression predictability.
thus our model is based on the simple assumption that genes regulated by the same transcription factors should exhibit similar expression profiles across different condition and contain common sequence motifs in their promoters.
discovered transcriptional modules reflect important processes in leaves
putative modules were defined as co-expression genes that could be predicted from sequence motifs in promoters.
significant co-expression was required across all  <dig> conditions for at least five genes.
a large number of overlapping modules were initially induced to capture the rich dynamics of the system.
these were then set to compete against each other in an algorithm that produced a final representative library of  <dig> modules covering  <dig> genes.
figure  <dig> shows two examples of these transcriptional modules, while all  <dig> modules are displayed in additional file  <dig>  the first module  contains all genes with the two motifs cr~msa-like and ma <dig> _gamyb in their promoters.
these motifs were over-represented in co-expressed genes .
over-represented functional annotations indicate a role in drought stress and nucleosome assembly.
indeed, a high expression correlation can be observed for these genes in the drought stress experiment .
the second example module  exhibits high expression similarity in the leaf primordial experiment  and annotations indicate a role in photosynthesis.
interestingly, one of the two motifs  is a known abscisic acid  response element, with aba having a role in many plant developmental processes.
most modules were significantly co-expressed within developmental processes such as leaf primordial and budset, while only a few modules were co-expressed in stress responses such as biotic infection and elevated  .
since the expression data are measured by two-channel microarrays, where stress-experiments typically used normal conditions as reference, this indicates that these stress-conditions activate rather different regulatory responses than do development.
a notable exception is drought stress, where all but one module exhibit significant co-expression, indicating that drought affects leaf development through these same modules.
interestingly, all of the three modules with a role in nucleosome assembly  belong to the very small number of modules with a significant co-expression in stress.
the relationship between nucleosome organization and stress has also been reported by others  <cit>  and may indicate a role for epigenetic modifications in response to stress.
one of the goals of this study was to investigate regulatory complexity.
interesting, very few of the discovered modules are associated with only one sequence motif .
typically two or three motifs were required to find a significant correspondence between motifs and co-expression, indicating a complex relationship between observed expression and the regulatory genome.
to evaluate the biological significance of the discovered modules, and their suggested regulatory control, we used functional annotations from gene ontology and kegg.
in general, 71% of the modules had some evidence of biological relevance in terms of over-represented gene ontology annotations  and kegg annotations .
many of these were related to photosynthesis and ribosomal activity, and thus of relevance to leaf development .
since all genes in this study were leaf-specific with a corresponding over-representation of leaf-specific annotations  <cit> , one could argue that any division of these genes into modules would produce relevant annotations.
however, in our statistical tests we used only the leaf-specific genes, not the whole genome, as background to avoid that typical leaf-functions show up as significant just because of the bias in the dataset.
hence, the large fraction of significant modules indicates that our division into modules based on common motifs and co-expression is indeed relevant.
this was also confirmed by randomization experiments, which invariably resulted in modules with considerably lower significance than reported here.
regulatory network indicates complex regulations
a regulatory network was inferred by applying regression models to predict the expression of genes in the transcriptional modules from the expression of sets of possible regulators .
the regression models increasingly included more transcription factors until the prediction performance of the more complex model  did not significantly improve on the simpler model .
a network was then drawn based on the best regulators of each module .
the method allowed us to identify the regulatory hubs of the leaf transcriptional program.
as in most biological networks, we observe a few hubs regulating many modules while most transcription factors only regulated a few modules .
a particularly strong hub was the transcription factor with protein id  <dig>  the closest homolog in arabidopsis is asil <dig> .
this factor belongs to the trihelix family of plant-specific transcriptional activators.
in our network, it is predicted to be involved in the regulation of all  <dig> photosynthesis genes that are overrepresented in transcriptional modules .
table  <dig> contains a full list of transcription factors predicted to have a regulatory role in populus leaves.
populus v <dig>  protein id is given together with information on the closest homologue in arabidopsis.
the last column gives the number of modules  regulated by the factor in our systems biology-based network and in the co-expression network, respectively.
transcription factors in our systems biology-based network that are not in the co-expression network are marked in bold.
our method of increasingly evaluating more complex regulatory mechanism allowed us to quantify the complexity of the regulation in populus leaves.
the distribution of modules over the number of transcription factors in the predicted regulatory mechanism  roughly follows that of the number of motifs .
thus, the predictive power of the regulatory mechanisms of most modules benefit significantly from including more than one transcription factor.
both steps in our method predict expression of genes, however, while the module discovery approach finds sequence motifs predictive of gene expression clusters, the network inference approach finds transcription factors predictive of the gene expression in each module.
both approaches are guided by the principle of occam's razor, that is, that the simplest model explaining the data is the best, and both approaches, as we have seen, result in the same distribution for the number of regulators per module.
the regression models describe the expression profiles in modules using the expression profiles of transcription factors.
in the case of two regulators, the expression of a module m is represented as a weighted sum of the expression of the regulators tf1and tf <dig>  i.e.
m = β <dig> + β1tf <dig> + β2tf <dig> β12tf1tf <dig>  thus, after fitting this model to the available expression data, the values of β <dig> and β <dig> will reflect the importance of each individual regulator, while the value of β <dig>  will reflect the importance of the interaction between the two regulators.
if the cross term is close to zero, there is a linear relationship between the module and the regulators, and not necessarily an interaction between the regulators.
a positive value of the cross term indicates a synergistic relationship between the regulators, while a negative value indicates a competitive relationship  <cit> .
figure  <dig> shows that individual regulators have a strong preference towards positive regulation over negative .
we also see slightly more synergistic than competitive relationships between regulators .
seven modules are governed by statistically significant synergistic interactions, while four modules exhibit competitive regulation .
the network is fully connected except for a small sub-network of the three nucleosome assembly modules discussed earlier.
one of these modules is shown in figure 2a, and is predicted to be regulated by  <dig> .
this factor is a histone protein with a known role in nucleosome assembly .
the other two modules are predicted to be regulated by  <dig> in concert with  <dig> , also a histone protein with a known role in nucleosome assembly.
the protein  <dig> is itself a member of the example module from figure 2a.
the fact that we did not allow auto-regulations in our inference method might thus be the reason why this module only has one regulator .
the two modules associated with both factors are the two modules with the strongest competitive regulatory mechanisms in the network .
both these regulators have a significant individual influence on the expression of the modules, but they also have a highly significant negative cross-term indicating the competitive regulation.
intriguingly, these are the only two modules in the network with a significant co-expression during biotic infection, although they are also co-expressed in a number of other experiments.
regulatory network predicts expression in unseen experiments
bootstrap analysis is often used in computational studies to evaluate the statistical significance of models such as phylogenetic trees  <cit> .
a bootstrap dataset has the same number of genes and conditions as the original data, but with some conditions occurring several time and some conditions not occurring at all .
on average,  <dig> % of the conditions will not occur in the bootstrap dataset and we refer to this as the hold-out set.
our network was validated statistically by first inferring a number of networks from different bootstrap dataset, and then  assessing the agreement between these bootstrap networks and the original network  and  using the regression models from the bootstrap networks to predict expression values in the hold-out sets .
most predicted regulations in the network recurred in a majority of the bootstrap networks .
however, about every third regulation had low support  .
three hubs  were the sources of  <dig> of these  <dig> weak regulations .
they are predicted to co-regulate modules with other, stronger regulators, and typically do not regulate modules by themselves.
thus these predicted regulatory interactions are sensitive to data removal and may only be valid under some experimental conditions.
our populus network models show a remarkable ability to generalize to unseen conditions, although similar predictive capability has been demonstrated also for other organisms  <cit> .
since we use the expression of a set of transcription factors to predict one expression profile per module, the correlation between observed and predicted expression is limited by the degree of expression similarity of genes within modules.
still, all co-expressed genes in modules had a significant correlation between observed and predicted expression when using the bootstrap networks to predict the expression in the hold-out sets .
in fact, 90% of genes, and all the modules, obtained a correlation above  <dig>  .
we also held out entire experiments  and used the resulting networks to predict the expression values in the missing experiment .
since few modules have a significant expression similarity within modules in stress responses , we are naturally unable to predict the expression in these experiments.
however, the regulation of the developmental programs, in particular leaf primordia and budset, can be predicted from the other experiments .
this is also true for drought stress, indicating that regulation of drought response corresponds to the regulation of development in that there is a conserved relationship between regulating transcription factors and regulated gene modules.
a notable exception is the nucleosome assembly modules from figure 2a with a role in water deprivation response.
this role is confirmed by the fact that the expression profile of this module cannot be predicted without the drought stress dataset .
several regulatory mechanisms are conserved between populus and arabidopsis
the aim of comparative genomics is usually to investigate the conservation of sequence across different species.
however, while proteins have diverged surprisingly little between related species, regulatory networks are believed to evolve much faster  <cit> .
our predictive approach makes it possible to investigate to what degree regulatory mechanisms of modules inferred from populus are conserved in other plant systems.
we applied the regression models from our populus inferred network to predict the expression of closest homologues in arabidopsis using the atgenexpress developmental conditions  <cit> .
since we were predicting the expression of arabidopsis genes from the expression of arabidopsis transcription factors, we were not testing the co-expression of these genes between the two plants.
rather, we were testing whether the regulatory mechanism, i.e.
the relationship between transcription factors and genes, is conserved.
of the  <dig> modules with expressed homologues in arabidopsis, 50% showed conservation beyond what would be expected by chance .
these  <dig> conserved modules cluster in three distinct parts of the network with functional roles in  biosynthesis, protein metabolism and translation,  carbon fixation, and  nucleosome assembly .
on the other hand, the non-conserved modules are almost exclusively over-represented for photosynthesis genes, showing a clear functional distinction between modules with conserved regulation in arabidopsis and those without.
interestingly, the photosynthesis modules contain co-expressed genes also in arabidopsis, although less so than the modules with conserved regulation .
thus, the model does not predict that the photosynthesis modules themselves have diverged between populus and arabidopsis, but rather that the regulation of these genes has been rewired.
this predicted rewiring of photosynthesis could be explained by the divergence in expression of the hub asil <dig>  we also repeated the analysis for arabidopsis expression data observed under abiotic stress , biotic stress  and various light conditions .
thus, the observation that abiotic stress  perturbs our modules to a less degree than biotic stress  also extends to the arabidopsis data.
other studies have also investigated the conservation of gene expression across populus and arabidopsis.
quesada et al.
<cit>  reported evidence of extensive evolution of gene expression regulation.
street et al.
<cit>  identified hub-genes in leaf development, and quantified the fraction of conserved genes to about 60%.
our results seem to imply similar conclusions, although the present study directly identified conserved relationships between transcription factors and gene modules.
an interesting question not addressed here is to what degree evolution of gene expression can be explained by divergence in the regulatory regions  of the two species.
systems biology predicts new leaf regulators
our approach describes interactions between regulators by inferring sets of transcription factors that regulate modules in concert.
this systems biology approach differs from traditional analysis such as hierarchical clustering or co-expression networks that only consider pair-wise similarity between the regulator and the regulated genes.
to compare these two approaches, we also constructed a co-expression network where each module is regulated by the single transcription factor with the most similar expression to that module .
table  <dig> lists all transcription factors in our systems biology-based network, and compares these to the regulators in this reductionistic co-expression network.
while the co-expression network identifies  <dig> transcription factors as regulators of populus leaf transcription, our network includes  <dig> of the  <dig> transcription factors in the data.
from figure  <dig> it was apparent that most collaborative regulations in our network have a master regulator, and this is the regulator typically identified by the co-expression network.
thus, most of the new regulators in our network are due to the fact that collaboration between transcription factors explains more of the expression in modules than single factors.
thus, although the regulations in our network are considerably stronger in terms of prediction power, they rarely exclude the transcription factor found in the co-expression network.
however, two modules were predicted by the co-expression network to be regulated by transcription factor  <dig> even though this most-similar factor was excluded as a regulator in our network.
somewhat surprisingly, the proposed regulatory mechanisms in these two cases are a weighted sum of two and three transcription factors without a statistically significant synergistic or competitive interaction .
discussion
one of the aims of systems biology is to model the complex interactions in living cells, describing emerging properties not apparent from studying genes, proteins or metabolites individually.
still, most computational approaches just take pair-wise similarity, not interactions between genes, into account when inferring network from expression data.
the reason for this is at least two-fold.
first, exploring combinatorics is computationally expensive.
for example, there are over  <dig> transcription factors in populus giving rise to over  <dig> million pairs,  <dig>  billion triplets, etc.
second, more complex models  imply many more parameters that have to be estimated from data .
since we need more observations than model-parameters to avoid over-fitting the models, the number of required observations grows quadraticly with the number of regulators when considering pairs.
this curse of dimensionality represents a huge obstacle to studying interactions in biological systems.
here, we deal with these problems in several different ways.
first, we restrict our study to leaf-specific genes rendering far fewer combinations than an unfiltered whole-genome study.
second, rather than considering all regulators at once, we devised a method that starts with single regulators, and then moves to pairs and higher-order combinations.
this provides adequate observations to estimate parameters for each model , but because we test so many different models it comes with the risk of finding combinations that obtain high predictive power by change .
we deal with this problem by only increasing model complexity if a statistically significant boost in predictive power is observed on unseen data .
in the statistical test we used the highly-conservative bonferroni correction where the initial significance threshold  was divided by the number of transcription factor combinations tested.
for the leaf-specific genes studied here, the systems biology-based network mostly discovered co-regulators to the transcription factors also identified in the co-expression network.
that means that although  <dig> of  <dig> modules had regulatory mechanisms with a significant interaction term , these regulators also had significant individual contributions of which the strongest is detected by pair-wise similarity.
a situation where the cross-term is significant, while the individual contributions are non-significant, is not observed in this data.
an example of such a regulation is the logical xor, that is, the regulated module is up-regulated only if one of the regulators is up-regulated .
whether such regulations exist in populus leaves cannot be settled from this study considering the limited set of genes included.
interestingly, the interaction term was non-significant in both cases where the best individual regulator was not part of our regulatory mechanism, meaning that the single best regulator was outperformed by a linear combination of other regulators.
such examples demonstrate how systems biology approaches have a better power to dissect regulatory complexity of biological systems than traditional approaches  <cit> .
they also show that systems biology is able to better model the 'real world' as qtl analysis of quantitative traits typically identifies numerous genetic loci, suggesting the involvement of numerous genes.
a particularly appealing feature of regression-based networks is their ability to predict expression of genes based on the expression of transcription factors.
we have used this to quantify the stability and predictive power of the network, but also to study module-conservation between experiments in populus and in arabidopsis.
several interesting predictions were found when studying modules that are co-expressed and correctly predicted using bootstrap networks, but that lose their predictability in particular experiment when these are entirely removed before network inference.
we have already mentioned the nucleosome assembly modules that are predicted to be regulated by histone h2a proteins.
the drought response profiles in these modules cannot be predicted by networks not trained on drought stress data.
another module  was affected by the removal of the budset data and is predicted to be regulated by factor  <dig>  a known cell death regulator.
the module characterized by motif os~tggca looses predictability without the dynamics of leaf growth-dataset, and is predicted to be controlled by a cell growth regulator .
genes in this module are also over-represented for carbon fixation.
prediction is a central theme in this study, and we strongly believe that predictive models have a lot to offer experimental biology as hypotheses generators.
the complete and correct regulatory network of an organism cannot be reverse-engineered from a limited collection of gene expression data.
however, we believe that such models represent a powerful starting point for further analysis as both hypothesis generation and descriptive tools.
the hubs in our network  thus represent attractive candidates for chip-seq analysis, functional knock-down studies and regulon engineering.
the network we present here only reflects the best regulators of each module.
however, behind each module in this network there is a ranked list of regulatory mechanisms , and as we have seen through bootstrap analysis, the ranking of these lists is not written in stone.
in the future one might hope that additional, and higher quality, data  will enable creation of more robust network models that more accurately reflect the underlying biological truth.
obviously, even a perfect network inference method cannot be better than the data it is modeled on .
another route to more reliable networks lies in combining computational inference with experimental testing in an iterative modeling approach.
several studies have shown how systematic perturbation of critical pathway components can be used to refine network representations  <cit> .
in plants, the lignin systems-project is taking this approach to model the lignin biosynthesis pathway  <cit> http://www.ligninsystems.org.
other sources of information may also be integrated into the network, but were not considered here, including epigenetic signatures such as nucleosome positioning and methylation patterns  <cit> , predicted binding site strength and transcription factor binding site preference  <cit> , and mirna regulation  <cit> .
CONCLUSIONS
we have outlined a systems biology model of the regulatory network of populus leaves.
the approach goes beyond previous analyses of populus transcriptome data by systematically considering interactions between transcription factors, leading us to predict new regulators of leaf development not found by traditional genomics methods.
these regulators orchestrate the transcriptional program in a synergistic or competitive manner, and thus constitute non-obvious targets for further analysis.
the model is robust when applied to predict expression levels in new data, and reveals conserved and diverged regulation both in different conditions within populus and between populus and arabidopsis.

 biologists have long been fascinated by the green plant leaf and have tried to understand how leaves are born, live and die. in the last decades, several new approaches to study the structure and function of leaves have emerged: molecular biology and molecular genetics have, for example, enabled identification of genes that regulate the primary function of the leaf - photosynthesis - and leaf development has been understood in much greater detail; high through-put transcriptomics has identified additional factors influencing leaf function, but traditional transcriptome analyses typically reduces the problem of finding key regulators to detecting differentially expressed genes or computing pair-wise similarity between targets and putative regulators. in contrast, systems biology analysis of transcriptional programs treats genes as interacting rather than isolated entities. thus these methods can begin to understand how so-called emergent properties such as complex phenotypes arise from interacting genes. whether this can be seen as taking a holistic rather than a reductionistic approach to science has generated quite some debate, but systems biology methods account for synergistic and competitive effects between regulators that individually could have low similarity to the target. methods for reverseengineering the transcriptional network from collections of gene expression data have been pioneered on single-cell organisms, but have increasingly been applied to higher order organisms including plants where applications of systems biology methods are now emerging. most systems biology studies have - not surprisingly - utilized using "the model plant" arabidopsis thaliana, where large transcriptomics programs have generated adequate quantities of high-quality data to enable systems analysis. for example, carerra et al. modeled the transcriptional network of arabidopsis and identified plant-specific properties such as high connectivity between genes involved in response and adaptation to changing environments. however, not all aspects of plant biology can be studied in arabidopsis, which in many respects is a rather atypical plant. indeed, it was not selected as a model system due to its physiological and ecological qualities, but rather for its suitability for genetic and genomic studies. therefore, it is important to perform parallel studies in plants with other characteristics, as well as developing the methods to allow data from the arabidopsis system to inform studies in other organisms. one rapidly emerging plant model system is populus; it's interesting biology and the access to a sequenced genome represent an attractive combination. correspondingly, more advanced data analyses approaches are now being applied in populus. populus provides an attractive model system for studies of leaf biology. for example, sjödin et al. exploited the fact that mature aspen in boreal regions have the rather unique property that all leaves emerge simultaneously from overwintering buds. this provides a synchronized system, resulting in a full temporal separation of the leaf developmental stages and subsequent acclimation that could be exploited using transcriptomics. access to a centralized repository of much of the populus cdna microarray data and databases for the analysis of gene expression - and other - data substantially facilitates the ability to perform systems biology studies. for example, grönlund et al. induced a co-expression network revealing modular architecture explaining gene function and tissue-specific expression; street et al. identified co-expression networks across a large collection of leaf transcriptomics data and found that some network hubs have existing functional evidence in arabidopsis; quesada et al. performed a comparative analysis of the transcriptomes of populus and arabidopsis, and found evidence of extensive remodeling of the transcriptional network, although some essential functions showed little divergence. a few studies have also integrated promoter information to study regulatory control in populus. shi et al. identified combinations of xylem-specific motifs in populus promoters. another study inferred transcriptional networks in xylem, leaves, and roots, and showed that genes with conserved regulation across tissues are primarily cis-regulated, while genes with tissue-specific regulation are often trans-regulated. all these studies are essentially co-expression networks that visualize expression similarity between pairs of genes, but do not infer complex interactions. network inference methods using expression data can be divided into those that aim to model the general influence that genes have on the expression of other genes and methods that aim to model the physical interaction between transcription factors and the regulated genes. both approaches employ common network inference methods, but those that infer gene regulatory networks also typically integrate motif finding and detection of transcriptional modules. approaches that describe how the regulatory genome orchestrates dynamic gene expression has developed from pilpel et al., who showed that yeast genes sharing pairs of binding sites in their promoters were significantly more likely to be co-expressed than genes sharing only single binding sites, to various machine learning methods that identify modules of co-expressed genes with common motif patterns in their promoters. here we apply a network inference method combining promoter information and expression data to describe the transcriptional network in populus leaves. our aims were to detect regulatory hubs in leaves, to describe conservation of transcriptional regulation within populus and between populus and arabidopsis, and to understand the regulatory complexity in leaves by comparing systems biology and traditional bioinformatics as methods for detecting target genes for further analysis. this study goes beyond previous meta-analyses of populus transcriptome data by taking into account synergistic and competitive interactions between regulators, and by systematically integrating the regulatory genome and the transcriptome to infer networks. we show that our network is robust, explains available gene function information and generalizes to new expression data in both populus and arabidopsis. we identify the main regulators of primary processes in leaves, and show how some of these have regulatory partners orchestrating expression either in a synergistic or competitive manner. such interactions are not considered by pair-wise similarity methods, and thus several of the regulators predicted here would not have been identified by traditional approaches.  we inferred the regulatory network of a collection of leaf-specific populus genes with quantified transcription profiles across samples in various experiments such as leaf primordial, budset, biotic infection and drought stress. the approach employed two separate steps to construct the network: first, we discovered a set of representative transcriptional modules containing co-expressed genes with evidence of co-regulation in their promoters. second, we inferred the most likely regulators of each module based on gene expression predictability. thus our model is based on the simple assumption that genes regulated by the same transcription factors should exhibit similar expression profiles across different condition and contain common sequence motifs in their promoters. discovered transcriptional modules reflect important processes in leaves putative modules were defined as co-expression genes that could be predicted from sequence motifs in promoters. significant co-expression was required across all conditions for at least five genes. a large number of overlapping modules were initially induced to capture the rich dynamics of the system. these were then set to compete against each other in an algorithm that produced a final representative library of modules covering genes. figure shows two examples of these transcriptional modules, while all modules are displayed in additional file the first module contains all genes with the two motifs cr~msa-like and ma _gamyb in their promoters. these motifs were over-represented in co-expressed genes. over-represented functional annotations indicate a role in drought stress and nucleosome assembly. indeed, a high expression correlation can be observed for these genes in the drought stress experiment. the second example module exhibits high expression similarity in the leaf primordial experiment and annotations indicate a role in photosynthesis. interestingly, one of the two motifs is a known abscisic acid response element, with aba having a role in many plant developmental processes. most modules were significantly co-expressed within developmental processes such as leaf primordial and budset, while only a few modules were co-expressed in stress responses such as biotic infection and elevated. since the expression data are measured by two-channel microarrays, where stress-experiments typically used normal conditions as reference, this indicates that these stress-conditions activate rather different regulatory responses than do development. a notable exception is drought stress, where all but one module exhibit significant co-expression, indicating that drought affects leaf development through these same modules. interestingly, all of the three modules with a role in nucleosome assembly belong to the very small number of modules with a significant co-expression in stress. the relationship between nucleosome organization and stress has also been reported by others and may indicate a role for epigenetic modifications in response to stress. one of the goals of this study was to investigate regulatory complexity. interesting, very few of the discovered modules are associated with only one sequence motif. typically two or three motifs were required to find a significant correspondence between motifs and co-expression, indicating a complex relationship between observed expression and the regulatory genome. to evaluate the biological significance of the discovered modules, and their suggested regulatory control, we used functional annotations from gene ontology and kegg. in general, 71% of the modules had some evidence of biological relevance in terms of over-represented gene ontology annotations and kegg annotations. many of these were related to photosynthesis and ribosomal activity, and thus of relevance to leaf development. since all genes in this study were leaf-specific with a corresponding over-representation of leaf-specific annotations, one could argue that any division of these genes into modules would produce relevant annotations. however, in our statistical tests we used only the leaf-specific genes, not the whole genome, as background to avoid that typical leaf-functions show up as significant just because of the bias in the dataset. hence, the large fraction of significant modules indicates that our division into modules based on common motifs and co-expression is indeed relevant. this was also confirmed by randomization experiments, which invariably resulted in modules with considerably lower significance than reported here. regulatory network indicates complex regulations a regulatory network was inferred by applying regression models to predict the expression of genes in the transcriptional modules from the expression of sets of possible regulators. the regression models increasingly included more transcription factors until the prediction performance of the more complex model did not significantly improve on the simpler model. a network was then drawn based on the best regulators of each module. the method allowed us to identify the regulatory hubs of the leaf transcriptional program. as in most biological networks, we observe a few hubs regulating many modules while most transcription factors only regulated a few modules. a particularly strong hub was the transcription factor with protein id the closest homolog in arabidopsis is asil. this factor belongs to the trihelix family of plant-specific transcriptional activators. in our network, it is predicted to be involved in the regulation of all photosynthesis genes that are overrepresented in transcriptional modules. table contains a full list of transcription factors predicted to have a regulatory role in populus leaves. populus v protein id is given together with information on the closest homologue in arabidopsis. the last column gives the number of modules regulated by the factor in our systems biology-based network and in the co-expression network, respectively. transcription factors in our systems biology-based network that are not in the co-expression network are marked in bold. our method of increasingly evaluating more complex regulatory mechanism allowed us to quantify the complexity of the regulation in populus leaves. the distribution of modules over the number of transcription factors in the predicted regulatory mechanism roughly follows that of the number of motifs. thus, the predictive power of the regulatory mechanisms of most modules benefit significantly from including more than one transcription factor. both steps in our method predict expression of genes, however, while the module discovery approach finds sequence motifs predictive of gene expression clusters, the network inference approach finds transcription factors predictive of the gene expression in each module. both approaches are guided by the principle of occam's razor, that is, that the simplest model explaining the data is the best, and both approaches, as we have seen, result in the same distribution for the number of regulators per module. the regression models describe the expression profiles in modules using the expression profiles of transcription factors. in the case of two regulators, the expression of a module m is represented as a weighted sum of the expression of the regulators tf1and tf i.e. m = β + β1tf + β2tf β12tf1tf thus, after fitting this model to the available expression data, the values of β and β will reflect the importance of each individual regulator, while the value of β will reflect the importance of the interaction between the two regulators. if the cross term is close to zero, there is a linear relationship between the module and the regulators, and not necessarily an interaction between the regulators. a positive value of the cross term indicates a synergistic relationship between the regulators, while a negative value indicates a competitive relationship. figure shows that individual regulators have a strong preference towards positive regulation over negative. we also see slightly more synergistic than competitive relationships between regulators. seven modules are governed by statistically significant synergistic interactions, while four modules exhibit competitive regulation. the network is fully connected except for a small sub-network of the three nucleosome assembly modules discussed earlier. one of these modules is shown in figure 2a, and is predicted to be regulated by. this factor is a histone protein with a known role in nucleosome assembly. the other two modules are predicted to be regulated by in concert with, also a histone protein with a known role in nucleosome assembly. the protein is itself a member of the example module from figure 2a. the fact that we did not allow auto-regulations in our inference method might thus be the reason why this module only has one regulator. the two modules associated with both factors are the two modules with the strongest competitive regulatory mechanisms in the network. both these regulators have a significant individual influence on the expression of the modules, but they also have a highly significant negative cross-term indicating the competitive regulation. intriguingly, these are the only two modules in the network with a significant co-expression during biotic infection, although they are also co-expressed in a number of other experiments. regulatory network predicts expression in unseen experiments bootstrap analysis is often used in computational studies to evaluate the statistical significance of models such as phylogenetic trees. a bootstrap dataset has the same number of genes and conditions as the original data, but with some conditions occurring several time and some conditions not occurring at all. on average, % of the conditions will not occur in the bootstrap dataset and we refer to this as the hold-out set. our network was validated statistically by first inferring a number of networks from different bootstrap dataset, and then assessing the agreement between these bootstrap networks and the original network and using the regression models from the bootstrap networks to predict expression values in the hold-out sets. most predicted regulations in the network recurred in a majority of the bootstrap networks. however, about every third regulation had low support. three hubs were the sources of of these weak regulations. they are predicted to co-regulate modules with other, stronger regulators, and typically do not regulate modules by themselves. thus these predicted regulatory interactions are sensitive to data removal and may only be valid under some experimental conditions. our populus network models show a remarkable ability to generalize to unseen conditions, although similar predictive capability has been demonstrated also for other organisms. since we use the expression of a set of transcription factors to predict one expression profile per module, the correlation between observed and predicted expression is limited by the degree of expression similarity of genes within modules. still, all co-expressed genes in modules had a significant correlation between observed and predicted expression when using the bootstrap networks to predict the expression in the hold-out sets. in fact, 90% of genes, and all the modules, obtained a correlation above. we also held out entire experiments and used the resulting networks to predict the expression values in the missing experiment. since few modules have a significant expression similarity within modules in stress responses, we are naturally unable to predict the expression in these experiments. however, the regulation of the developmental programs, in particular leaf primordia and budset, can be predicted from the other experiments. this is also true for drought stress, indicating that regulation of drought response corresponds to the regulation of development in that there is a conserved relationship between regulating transcription factors and regulated gene modules. a notable exception is the nucleosome assembly modules from figure 2a with a role in water deprivation response. this role is confirmed by the fact that the expression profile of this module cannot be predicted without the drought stress dataset. several regulatory mechanisms are conserved between populus and arabidopsis the aim of comparative genomics is usually to investigate the conservation of sequence across different species. however, while proteins have diverged surprisingly little between related species, regulatory networks are believed to evolve much faster. our predictive approach makes it possible to investigate to what degree regulatory mechanisms of modules inferred from populus are conserved in other plant systems. we applied the regression models from our populus inferred network to predict the expression of closest homologues in arabidopsis using the atgenexpress developmental conditions. since we were predicting the expression of arabidopsis genes from the expression of arabidopsis transcription factors, we were not testing the co-expression of these genes between the two plants. rather, we were testing whether the regulatory mechanism, i.e. the relationship between transcription factors and genes, is conserved. of the modules with expressed homologues in arabidopsis, 50% showed conservation beyond what would be expected by chance. these conserved modules cluster in three distinct parts of the network with functional roles in biosynthesis, protein metabolism and translation, carbon fixation, and nucleosome assembly. on the other hand, the non-conserved modules are almost exclusively over-represented for photosynthesis genes, showing a clear functional distinction between modules with conserved regulation in arabidopsis and those without. interestingly, the photosynthesis modules contain co-expressed genes also in arabidopsis, although less so than the modules with conserved regulation. thus, the model does not predict that the photosynthesis modules themselves have diverged between populus and arabidopsis, but rather that the regulation of these genes has been rewired. this predicted rewiring of photosynthesis could be explained by the divergence in expression of the hub asil we also repeated the analysis for arabidopsis expression data observed under abiotic stress, biotic stress and various light conditions. thus, the observation that abiotic stress perturbs our modules to a less degree than biotic stress also extends to the arabidopsis data. other studies have also investigated the conservation of gene expression across populus and arabidopsis. quesada et al. reported evidence of extensive evolution of gene expression regulation. street et al. identified hub-genes in leaf development, and quantified the fraction of conserved genes to about 60%. our results seem to imply similar conclusions, although the present study directly identified conserved relationships between transcription factors and gene modules. an interesting question not addressed here is to what degree evolution of gene expression can be explained by divergence in the regulatory regions of the two species. systems biology predicts new leaf regulators our approach describes interactions between regulators by inferring sets of transcription factors that regulate modules in concert. this systems biology approach differs from traditional analysis such as hierarchical clustering or co-expression networks that only consider pair-wise similarity between the regulator and the regulated genes. to compare these two approaches, we also constructed a co-expression network where each module is regulated by the single transcription factor with the most similar expression to that module. table lists all transcription factors in our systems biology-based network, and compares these to the regulators in this reductionistic co-expression network. while the co-expression network identifies transcription factors as regulators of populus leaf transcription, our network includes of the transcription factors in the data. from figure it was apparent that most collaborative regulations in our network have a master regulator, and this is the regulator typically identified by the co-expression network. thus, most of the new regulators in our network are due to the fact that collaboration between transcription factors explains more of the expression in modules than single factors. thus, although the regulations in our network are considerably stronger in terms of prediction power, they rarely exclude the transcription factor found in the co-expression network. however, two modules were predicted by the co-expression network to be regulated by transcription factor even though this most-similar factor was excluded as a regulator in our network. somewhat surprisingly, the proposed regulatory mechanisms in these two cases are a weighted sum of two and three transcription factors without a statistically significant synergistic or competitive interaction. discussion one of the aims of systems biology is to model the complex interactions in living cells, describing emerging properties not apparent from studying genes, proteins or metabolites individually. still, most computational approaches just take pair-wise similarity, not interactions between genes, into account when inferring network from expression data. the reason for this is at least two-fold. first, exploring combinatorics is computationally expensive. for example, there are over transcription factors in populus giving rise to over million pairs, billion triplets, etc. second, more complex models imply many more parameters that have to be estimated from data. since we need more observations than model-parameters to avoid over-fitting the models, the number of required observations grows quadraticly with the number of regulators when considering pairs. this curse of dimensionality represents a huge obstacle to studying interactions in biological systems. here, we deal with these problems in several different ways. first, we restrict our study to leaf-specific genes rendering far fewer combinations than an unfiltered whole-genome study. second, rather than considering all regulators at once, we devised a method that starts with single regulators, and then moves to pairs and higher-order combinations. this provides adequate observations to estimate parameters for each model, but because we test so many different models it comes with the risk of finding combinations that obtain high predictive power by change. we deal with this problem by only increasing model complexity if a statistically significant boost in predictive power is observed on unseen data. in the statistical test we used the highly-conservative bonferroni correction where the initial significance threshold was divided by the number of transcription factor combinations tested. for the leaf-specific genes studied here, the systems biology-based network mostly discovered co-regulators to the transcription factors also identified in the co-expression network. that means that although of modules had regulatory mechanisms with a significant interaction term, these regulators also had significant individual contributions of which the strongest is detected by pair-wise similarity. a situation where the cross-term is significant, while the individual contributions are non-significant, is not observed in this data. an example of such a regulation is the logical xor, that is, the regulated module is up-regulated only if one of the regulators is up-regulated. whether such regulations exist in populus leaves cannot be settled from this study considering the limited set of genes included. interestingly, the interaction term was non-significant in both cases where the best individual regulator was not part of our regulatory mechanism, meaning that the single best regulator was outperformed by a linear combination of other regulators. such examples demonstrate how systems biology approaches have a better power to dissect regulatory complexity of biological systems than traditional approaches. they also show that systems biology is able to better model the'real world' as qtl analysis of quantitative traits typically identifies numerous genetic loci, suggesting the involvement of numerous genes. a particularly appealing feature of regression-based networks is their ability to predict expression of genes based on the expression of transcription factors. we have used this to quantify the stability and predictive power of the network, but also to study module-conservation between experiments in populus and in arabidopsis. several interesting predictions were found when studying modules that are co-expressed and correctly predicted using bootstrap networks, but that lose their predictability in particular experiment when these are entirely removed before network inference. we have already mentioned the nucleosome assembly modules that are predicted to be regulated by histone h2a proteins. the drought response profiles in these modules cannot be predicted by networks not trained on drought stress data. another module was affected by the removal of the budset data and is predicted to be regulated by factor a known cell death regulator. the module characterized by motif os~tggca looses predictability without the dynamics of leaf growth-dataset, and is predicted to be controlled by a cell growth regulator. genes in this module are also over-represented for carbon fixation. prediction is a central theme in this study, and we strongly believe that predictive models have a lot to offer experimental biology as hypotheses generators. the complete and correct regulatory network of an organism cannot be reverse-engineered from a limited collection of gene expression data. however, we believe that such models represent a powerful starting point for further analysis as both hypothesis generation and descriptive tools. the hubs in our network thus represent attractive candidates for chip-seq analysis, functional knock-down studies and regulon engineering. the network we present here only reflects the best regulators of each module. however, behind each module in this network there is a ranked list of regulatory mechanisms, and as we have seen through bootstrap analysis, the ranking of these lists is not written in stone. in the future one might hope that additional, and higher quality, data will enable creation of more robust network models that more accurately reflect the underlying biological truth. obviously, even a perfect network inference method cannot be better than the data it is modeled on. another route to more reliable networks lies in combining computational inference with experimental testing in an iterative modeling approach. several studies have shown how systematic perturbation of critical pathway components can be used to refine network representations. in plants, the lignin systems-project is taking this approach to model the lignin biosynthesis pathway http://www.ligninsystems.org. other sources of information may also be integrated into the network, but were not considered here, including epigenetic signatures such as nucleosome positioning and methylation patterns, predicted binding site strength and transcription factor binding site preference, and mirna regulation.  we have outlined a systems biology model of the regulatory network of populus leaves. the approach goes beyond previous analyses of populus transcriptome data by systematically considering interactions between transcription factors, leading us to predict new regulators of leaf development not found by traditional genomics methods. these regulators orchestrate the transcriptional program in a synergistic or competitive manner, and thus constitute non-obvious targets for further analysis. the model is robust when applied to predict expression levels in new data, and reveals conserved and diverged regulation both in different conditions within populus and between populus and arabidopsis. 
the network is reverse-engineered from promoter information and expression profiles of leaf-specific genes measured over a large set of conditions related to stress and developmental. the network model incorporates interactions between regulators, such as synergistic and competitive relationships, by evaluating increasingly more complex regulatory mechanisms, and is therefore able to identify new regulators of leaf development not found by traditional genomics methods based on pair-wise expression similarity. we present a systems biology model of the regulatory network of populus leaves. although systems biology models should be used with care considering the complexity of regulatory programs and the limitations of current genomics data, methods describing interactions can provide hypotheses about the underlying cause of emergent properties and are needed if we are to identify target genes other than those constituting the "low hanging fruit" of genomic analysis. today, comprehensive databases of gene expression data enable us to apply increasingly more advanced computational methods for reverse-engineering the regulatory network of leaves, and to begin to understand the gene interactions underlying complex emergent properties related to stress-response and development. we outline a computationally inferred model of the regulatory network of populus leaves, and show how treating genes as interacting, rather than individual, entities identifies new regulators compared to traditional genomics analysis. the approach is shown to explain available gene function information and to provide robust prediction of expression levels in new data.
 green plant leaves have always fascinated biologists as hosts for photosynthesis and providers of basic energy to many food webs. today, comprehensive databases of gene expression data enable us to apply increasingly more advanced computational methods for reverse-engineering the regulatory network of leaves, and to begin to understand the gene interactions underlying complex emergent properties related to stress-response and development. these new systems biology methods are now also being applied to organisms such as populus, a woody perennial tree, in order to understand the specific characteristics of these species.  we present a systems biology model of the regulatory network of populus leaves. the network is reverse-engineered from promoter information and expression profiles of leaf-specific genes measured over a large set of conditions related to stress and developmental. the network model incorporates interactions between regulators, such as synergistic and competitive relationships, by evaluating increasingly more complex regulatory mechanisms, and is therefore able to identify new regulators of leaf development not found by traditional genomics methods based on pair-wise expression similarity. the approach is shown to explain available gene function information and to provide robust prediction of expression levels in new data. we also use the predictive capability of the model to identify condition-specific regulation as well as conserved regulation between populus and arabidopsis.  we outline a computationally inferred model of the regulatory network of populus leaves, and show how treating genes as interacting, rather than individual, entities identifies new regulators compared to traditional genomics analysis. although systems biology models should be used with care considering the complexity of regulatory programs and the limitations of current genomics data, methods describing interactions can provide hypotheses about the underlying cause of emergent properties and are needed if we are to identify target genes other than those constituting the "low hanging fruit" of genomic analysis. 