BACKGROUND
functional genomics is a field of molecular biology that studies the functions and interactions between genes and proteins from large datasets generated by genome projects.
over recent years, the development of new high throughput technologies among different omics approaches has led to the availability of a large volume of transcriptomic, metabolomic, proteomic, physiological and phenotypic data.
in this sense, the multidisciplinary synergy between molecular biology, statistics and informatics emerge as necessary to support the use and interpretation of the results in the area of functional genomics in order to provide specific tools for the simultaneous analysis of many genes.
senescence is the final stage of leaf development and precedes cell death.
once the senescence program is triggered, important highly regulated gene expression changes occur, leading to important changes in the metabolism .
in sunflower, which is the fourth most important oil crop worldwide, a delay in leaf senescence has a strong impact on yield; this effect has also been shown in other crops.
maintaining the photosynthetic leaf area especially during the reproductive stage  has been shown to impact on the gap between potential and real yield observed, this is due to the incapacity of plants to maintain their green leaf area for longer periods  <cit> .
transcription factors  are key proteins involved in the regulation of gene expression and signal transduction networks, regulating different biological processes.
several transcription factor families have been associated with leaf senescence in a range of species  and their function is crucial for triggering and/or regulating the different aspect of this process.
various software tools have been developed with the aim of analyzing datasets to predict hub genes which could regulate a specific pathway or be involved in a particular biological process.
weighted gene correlation network analysis  is an r package method designed to find clusters  of highly correlated genes or metabolites.
wgcna calculates a pearson’s correlation matrix for all the genes and then transforms the correlation matrix into an adjacency matrix by raising all values to a soft threshold power; this has the effect of emphasizing the strong correlations and penalizing the weaker correlations.
the module eigengene  is the first principal component of a given module and can be considered as representative of the module’s gene expression profile.
modules often represent specific biological processes; highly connected hub genes within the module are often regulatory genes and represent candidate biomarker.
wgcna can be used to construct networks in which each node represents a gene or metabolite and the connecting lines  represent correlations between the nodes.
wgcna has been widely used as a method to cluster gene expression and metabolite data and identify hub biomarker genes in several non-model plant species such as tomato  <cit> , brassica  <cit> , petunia  <cit> , rice  <cit>  and ficus  <cit> .
biosignature discoverer  is a software application devised for identifying molecular signatures in different biological datasets, such as next generation sequencing, microarray data and metabolic profiles, in a statistically robust, computationally efficient, and user-friendly way.
a signature is defined as the minimal set of molecular quantities that collectively yield maximal predictive performance.
thus, a signature does not contain irrelevant or redundant quantities, given the selected ones.
in addition, there may be numerous equivalent signatures that lead to equally accurate predictions, employing different molecular quantities.
the core of this application uses a feature selection algorithm that belongs to the class of bayesian network, constraint-based learning  <cit>  and that is able to identify multiple, statistically equivalent signatures.
in this work we used these two complementary methods, wgcna and biosd, with the aim of identifying potential transcription factors and metabolites, as putative biomarkers associated with leaf senescence in sunflower.
methods
plant material and experimental conditions
the sunflower leaf senescence assays were conducted under field conditions.
the experiment was carried out at the inta balcarce experimental station  as previously described  <cit> .
the experiment was sowed during the 2010/ <dig> growing season.
the sunflower hybrid vdh  <dig>  was sown at a  <dig>  plants/m <dig> with three biological replicates , each one consisted of three randomly selected plants from each plot.
sunflower assay was conducted under control conditions, without limitations .
transcriptomic and metabolomics profiling was performed using the leaf  <dig>  at three different development stages labeled as t- <dig> , t- <dig>  and t- <dig>   <cit> .
the samples were immediately frozen in liquid nitrogen upon collection and saved at −80 °c until processing.
high quality total rna was isolated from 100 mg of frozen tissue using tripure, according to the manufacturer’s instructions .
the genomic dna was eliminated by treatment with dnase i for 20 min at room temperature .
the rna concentration was measured using a nanodrop nd- <dig> spectrophotometer .
the purity and integrity of total rna was determined by 260/280 nm ratio and by nanobioanalyzer rna- <dig> analysis .
data set
transcriptomic profiling was performed using a custom sunflower microarray  which has been previously described  <cit> .
background correction was performed using the rma algorithm from the backgroundcorrect function .
between arrays normalization was performed using the quantile method from the normalizebetweenarray function.
finally, gene expressions were transformed to log <dig> scale and information from technical replicates was incorporated by calculating the median parameter.
statistically significant probesets where identified using the lmfit and contrasts.fit functions from the r software limma bioconductor package  <cit>  a total of  <dig>  probes were shown to be statistically significant during leaf senescence as the plants develop .
the statistical parameters were p-value lower than  <dig>  and fold-change higher or lower than  <dig>  <cit> .
sunflower transcription factors  were identified by comparing approximately  <dig>  tfs sequences from arabidopsis lyrata, arabidopsis thaliana, oryza sativa, populus trichocarpa, vitis vinifera and zea mays available from the plant transcription factor database   <cit>  with sur v <dig>  database  <cit>  using blast software.
a total of  <dig> tfs were differentially expressed during leaf senescence.
metabolic profiling was performed using the gc-tof-ms system .
metabolite extraction was performed by promoting the extraction of lipophilic and polar compounds according to recently published protocols  <cit>  adapted for sunflower tissue samples  <cit> .
the chromatograms and spectra were evaluated using chromatof  and tagfinder  <cit> .
ion spectra were compared to the golm metabolome database .
metabolite levels were normalized to fresh weight using ribitol as the internal control.
changes in metabolite levels along leaf development were calculated as fold-change relative to the first sampled time .
weighted gene co-expression network analysis
wgcna was performed using the wgcna r package  as described by langfelder and horvath  <cit> .
the expression values for  <dig>  non-redundant genes and levels of  <dig> metabolites for the  <dig> samples were used to construct two separate networks.
samples were clustered using the function hclust to check that there were no outliers.
the picksoftthreshold function was used to select the soft threshold power used to construct a network based on the criterion of approximate scale-free topology.
the power value selected was the lowest power for which the scale-free topology fit index curve flattened out upon reaching a high value.
the next step was to transform the adjacency matrix into a topological overlap matrix , which summarizes the degree of shared connections between any two genes.
the tom matrix was then converted into a dissimilarity matrix.
genes were then clustered using the average linkage hierarchical clustering and the modules were identified in the resulting dendrogram using the dynamic hybrid tree cut method .
found modules were trimmed of genes whose correlation with the module eigengene  was less than minkmetostay.
if p-values of the higher correlations were smaller than those of the native module by the factor reassignthreshold, the gene was re-assigned to the closer module.
modules whose eigengenes were highly correlated were merged.
this was achieved by clustering module eigengenes using the dissimilarity given by one minus their correlation, cutting the dendrogram at the height mergecutheight and merging all modules on each branch below mergecutheight.
this was implemented using the blockwisemodules function in wgcna which performs the network construction and consensus module detection.
the following settings were used for the gene network; power =  <dig>  minmodulesize =  <dig>  mergecutheight =  <dig> , maxblocksize =  <dig>  deepsplit =  <dig>  reassignthreshold = 1e- <dig> and minkmetostay =  <dig> , networktype = “unsigned”, tomtype = “signed”.
the genes were divided into  <dig> modules; the module sizes ranged from  <dig> to  <dig>  genes,  <dig> genes were not assigned to any of the modules.
the following settings were used for the metabolite network; power =  <dig>  minmodulesize =  <dig>  mergecutheight =  <dig> , maxblocksize =  <dig>  deepsplit =  <dig>  reassignthreshold = 1e- <dig> and minkmetostay =  <dig> , networktype = “unsigned”, tomtype = “signed”.
the metabolites were divided into  <dig> modules; the module sizes ranged from  <dig> to  <dig> metabolites,  <dig> metabolites were not assigned to any of the modules.
module membership  is a quantitative measure of module membership; this is the correlation of the me and the gene expression profile for each gene in the module.
the gene with the highest connectivity in each module was returned using the function choosetophubineachmodule.
the networks were exported from wgcna using the function exportnetworktocytoscape.
the gene network showed approximate scale free topology, whereas the metabolic network did not satisfy scale free topology:  <dig> genes scalefreersquared  <dig>  slope − <dig>  and  <dig> metabolites scalefreersquared  <dig>  slope − <dig>  .
a total of  <dig> metabolites with a degree higher than  <dig> were then selected as metabolite hubs and correlated with the list of  <dig> tfs statistically significant during senescence  using pearson correlation.
correlations with p-value < <dig>  were selected and visualized using cytoscape  <cit> .
biosignature discoverer
the statistical analysis pipeline of biosd, a plugin of qiagen’s clc-bio workbench , has been used for identifying gene-expression signatures able to discriminate the age of sunflower plants.
particularly, the plugin has been applied on the expression values of  <dig> transcription factors that were found differentially expressed during leaf senescence, on the  <dig> metabolites measured in the study, and on the combination of these two data types.
additionally, the pipeline has been separately applied on the extended list of  <dig>  non-redundant genes.
biosd employs a complex machine-learning approach in order to discover such signatures and for quantifying their actual predictive power.
particularly, this tool employs novel, proprietary feature selection algorithms to identify the signatures, inspired by prior work on bayesian-network constraint-based learning  <cit> .
notably, the pipeline is also able to identify multiple equivalent signatures, whose predictive performances are statistically equivalent.
to avoid under fitting, the tool automatically tries multiple algorithms for feature selection and prediction modeling, as well as multiple values for tuning them.
the tool also employs sophisticated methods, based on an enhanced version of cross-validation  <cit> , for providing conservative estimations of the performances that can be achieved with the given signatures.
the biosd interface allows the user to modulate the trade-off between quality of the results and computational requirements through the ‘tuning effort’ parameter.
this parameter regulates the extent to which the tool will attempt tuning the feature selection and prediction modeling algorithms.
in all analyses, the ‘tuning effort’ was set to ‘extensive’, in order to favor result accuracy at the expense of higher computational time.
RESULTS
wgcna is a powerful method used to identify clusters of highly correlated genes which are potentially co-regulated.
all the connections for the metabolites with an edge weight > <dig>  produced a network of  <dig> nodes and  <dig> edges.
the connections between the  <dig> significant differentially expressed transcription factors were selected from the gene network with an edge weight > <dig> , producing a network of  <dig> nodes and  <dig> edges.
the networks were visualized in cytoscape by edge weight and degree .fig.
<dig> networks exported from wgcna visualized in cytoscape.
network were constructed using  <dig> transcription factors  and  <dig> metabolites  and were visualized in cytoscape  <cit>  with an edge weight higher than  <dig> .
the nodes represent genes and metabolites and the edges represent connections between them.
the node size and color is related to the number of connections, large red nodes represent highly connected hub genes and metabolites, small green nodes represent gene with few connections.
strong connections are visualized as wider lines
a total of  <dig> transcription factors and  <dig> metabolites showed the highest number of connections , over 80 % of the maximum number of connections observed , highlighting them as potentially co-regulated hub biomarkers during the senescence process in sunflower.
pearson correlation analysis between the  <dig> hub metabolites and the list of  <dig> tfs significant during senescence was performed in order to integrate both analyses.
the resulting network contained  <dig> nodes and  <dig> edges .fig.
<dig> integrated network of hubs metabolites and transcription factors.
pearson correlation analysis of  <dig> selected hubs metabolites and the list of  <dig> tfs statistically significant during senescence.
correlations with p-value <  <dig>  were selected and visualized in cytoscape  <cit>  by degree  and edge weight
the integration results showed that  <dig> metabolite hubs correspond to sugars and tca cycle metabolites and the nac tfs family, widely reported as senescence associated, showed  <dig> members  as hub genes.
on the other hand, biosd allows identification of the most informative genes and metabolites for discriminating between the various stages of the leaf senescence process.
table  <dig> lists all the transcription factors and metabolites whose values allow the predicting of senescence in sunflower plants.
all these quantities have equivalent discriminative power and 100 % stability, which mean that they are consistently selected even when other datasets are included.table  <dig> list of biosignatures detected in each dataset
each signature is composed of a single biomarker, and all signatures are expected to have the same predictive performance
table  <dig> reports the ‘in-sample’ and ‘out-of-sample’ performances of the selected transcripts and metabolites.
the in-sample values quantify the fitness of the predictive models, and poor results would indicate presence of outliers, unmeasured yet needed predictive factors, or poorly modeled trends in the data.
the out-of-sample values estimate the performances that are expected if the selected expression and metabolic values are used for predicting the age of new, independent samples of sunflower plants .
assuming y and ŷ are the observed and predicted age of n plants.
the r <dig> metric measures the proportion of y variance explained by the predictive model, r <dig> values close to  <dig> are indicative of good fit, while the mean absolute error  and mean squared error  statistics quantify predictions’ deviations from the actual age as mae=1n⋅∑absy−y^ and mse=1n⋅∑y−y^ <dig> table  <dig> predicted performances of the selected biosignatures
performances are reported in terms of the determination coefficient r <dig>  mean absolute error  and mean squared error .
the in-sample performances quantify the fitness of the predictive models on the training data, while the out-of-sample values estimate the expected performance on new data.
confidence intervals are calculated using a bootstrap approach
in the present application both in- and out-of-sample performances show quite good results, indicating  that the predictive models adequately fit the available samples and  that the selected transcript factors and metabolites have a high predictive power.
mae, a robust estimator of the standard deviation, is around  <dig>  days for the “out of sample” for both, metabolite and transcript data, meaning that it should be possible to predict the age of a new plant on the basis of transcriptional or metabolic information with an average error of approximately three days.
residual vs. predicted value plots did not showed any apparent pattern or trend, indicating that the models have fitted the data adequately.
the analysis of the complete list of gene  integrated with metabolic information identified twelve distinct signatures each composed of two elements.
particularly, each signature contains the gene hean_c_ <dig>  coupled either with the alanine metabolite or with one of the following eleven genes: hean_c_ <dig>  hean_c_ <dig>  hean_c_ <dig>  hean_c_ <dig>  hean_c_ <dig>  hean_c_ <dig>  hean_c_ <dig>  hean_s_ <dig>  hean_s_ <dig>  hean_s_ <dig>  hean_s_ <dig>  these signatures have slightly better performances than the single biomarkers reported in table  <dig> .
these results indicate that the alanine metabolite is interchangeable, in terms of predictive power with each of the other eleven genes.
thus, practitioners interested in deploying this signature in practice  can freely decide whether to build a predictive model using exclusively transcriptional information or a mixture of transcriptional and metabolic data.
this decision can depend on multiple factors, particularly considerations about the cost and technical practicability of measuring the quantities included in the signature.
in recent years, many studies have focused on the detection of different genes up or down-regulated along the progression of the leaf senescence process, identifying tfs that could be triggering the process .
hence, studying interactions between these tfs under different conditions, for a non-model organism, will open new insights into the details of gene regulation by identifying new interactions and comparing them to those already observed in model organisms  <cit> .
one of the most important aims in new generation biotech crops is to increase the yield and to improve tolerance to different stresses.
although sequencing-based approaches are becoming the method of choice for gene expression profiling compared to hybridization-based ones, for non-model species for which a reference sequence is not available, the hybridization approach can still be the more reliable method  <cit> .
this study reports on the use of two complementary software applications with the aim of producing a powerful network and predicting potential biosignatures enabling the identification of important tfs and metabolites involved during the leaf senescence process in sunflower.
using wgcna analysis, genes and metabolites are clustered in one or more modules.
genes in each module are highly correlated, and consequently they are likely to be co-regulated during the senescence process.
biosd identifies which genes are the most informative in order to discriminate between different stages of leaf development.
the selected genes, in this case, are not necessarily correlated among themselves.
in other words, wgcna attempts to identify the elements involved in the senescence process, and to gain insight of their interplay.
on the other hand, biosd tries to retain the minimum number of these elements that are needed for optimal prediction of the plant senescence process.
wgcna allowed the detection of metabolite hubs during leaf senescence, most of them correspond to sugars and tca cycle metabolites.
sugars are central elements in the sink-source relationships  <cit>  and have been reported as growth  <cit>  and photosynthesis regulators  <cit> .
sunflower is a crop plant with a strong demand for nutrients, especially sugar as a substrate in the oil synthesis during grain filling; this could potentially act as a trigger signal for the senescence process.
biosd allowed the detection of one metabolite and ten tfs as potential biomarkers.
the selected biosignatures are able to estimate the senescence of a plant with an error of approximately  <dig>  days.
when both tf and metabolite data are simultaneously analyzed, the expected error decreases to  <dig>  days, indicating that the integration of different omic data allows better predictions of sunflower senescence by biosd analysis.
transcription factors are major players and some of them constitute major hubs in signaling pathways.
in both analyses, the selected genes belong to transcription factor families widely reported as senescence associated in model species.
when we analyze the integration of metabolites and tfs  we found that three tfs  are detected as biomarkers independently in the two methods.
the expression profiles of the selected biomarkers during leaf senescence are displayed in the additional file  <dig>  hean_c_ <dig> presents high sequence similarity to aterf <dig>  which encodes a member of the erf .table  <dig> candidate biomarkers detected by the two complementary methods, wgcna and biosignature discoverer
hean_c_10939
at1g50640
ap2
hean_s_32569
at5g49620
myb
hean_c_11118
at1g69490
nac
hean_c_10939
at1g50640
ap2
hean_c_11118
at1g69490
nac
hean_s_32569
at5g49620
myb
biomarkers from wgcna were selected based on their degree  and biomarkers from biosignature discoverer correspond to the analysis of the integrated list.
biomarkers detected independently in the two methods are in bold
transgenic arabidopsis plants with enhanced expression of this tf showed precocious leaf senescence  <cit> .
additionally, hean_c_ <dig> has high sequence similarity to atnap/anac <dig> , a nac transcription factor which acts downstream of ein <dig> and ein <dig> genes.
ein <dig> positively regulates leaf senescence by activating ore <dig> and atnap.
genetic analysis suggest that both genes act in distinct and overlapping signaling pathways regulating leaf senescence in arabidopsis  <cit> .
hanac <dig>  a putative sunflower orthologous gene to ore <dig>  has previously been reported in sunflower as a candidate biomarker for leaf senescence  <cit> , which reinforces these genes as potential regulators of the senescence process in sunflower.
hean_s_ <dig> has high sequence similarity to myb <dig> .
this tf has not yet been directly associated with the senescence process, however, members of this family have been reported as associated to senescence and involved in response to several abiotic stresses  <cit> .
these results suggest that wgcna and the identification of molecular biosignatures are powerful tools for the detection of potential biomarkers.
however, it should be noted that these results are obtained on a relatively small experiment .
further studies would be advisable in order to better estimate the predictive performance of this approach and to achieve best fit for the list of candidate omics signatures.
the implementation of data integration analysis associated to different biotic and abiotic stresses is a powerful tool, especially in a non-model crop such as sunflower for which the complete genome sequence is not yet available.
these results open new strategies of analysis to explore and detect potential biomarkers associated with leaf senescence that will be useful for future molecular breeding programs.
CONCLUSIONS
the complementary analysis of transcriptomic and metabolomic data with wgcna and biosd emerge as a useful strategy to predict not only the age of the plants but also to detect and identify potential biomarkers associated with leaf senescence.
in the case of sunflower, a worldwide oil crop, the implementation of wgcna enabled the construction of a hypothetical network which was used to infer relationships between tfs and metabolites and the identification of hubs as potential candidate biomarkers involved in leaf senescence, whereas biosd selected transcripts and metabolites to discriminate the age of sunflower plants.
this methodological approach is a novel strategy for omics data integration, it can be used to elucidate and predict novel networks as well as identifying potential transcripts as putative biomarker hub genes at different developmental stages.
additional files
additional file 1: cluster dendrogram visualizing the modules from wgcna.
a genes, b metabolites were clustered using the average linkage hierarchical clustering and modules identified in the resulting dendrogram by the dynamic hybrid tree cut method.
additional file 2: scale-free topologies.
a genes, b metabolites.
the left panel shows a histogram of network connectivity.
the right panel shows a log-log plot of the histogram.
the approximate straight line relationship  shows approximate scale free topology.
additional file 3: list of  <dig> hubs metabolites and  <dig> transcription factors.
the list shows the identifiers and the normalized expression data of the different samples.
additional file 4: biosignature discoverer results on predicting senescence from integrated metabolomics and transcriptomics data in sunflower.
additional file 5: expression profiles of the selected biomarkers during leaf senescence.
the list shows the identifiers, the expression profile of candidate biomarkers in log <dig> scale and the standardization of this expression profile.
the figures show the genes and the metabolites in red and blue lines respectively.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
pf and sm conceived and designed the work, contributed to the implementation of wgcna for the data, tested software to visualize and understand the data analysis and wrote the first draft of the manuscript, jh implemented wgcna and tested different parameters to optimize the network design, jadr revised the manuscript and statistically supported software analysis, rh and np initiated the project and contributed to the work by the interpretation, discussion of the data and critically revised the manuscript.
all authors read and approved the final manuscript.

 functional genomics is a field of molecular biology that studies the functions and interactions between genes and proteins from large datasets generated by genome projects. over recent years, the development of new high throughput technologies among different omics approaches has led to the availability of a large volume of transcriptomic, metabolomic, proteomic, physiological and phenotypic data. in this sense, the multidisciplinary synergy between molecular biology, statistics and informatics emerge as necessary to support the use and interpretation of the results in the area of functional genomics in order to provide specific tools for the simultaneous analysis of many genes. senescence is the final stage of leaf development and precedes cell death. once the senescence program is triggered, important highly regulated gene expression changes occur, leading to important changes in the metabolism. in sunflower, which is the fourth most important oil crop worldwide, a delay in leaf senescence has a strong impact on yield; this effect has also been shown in other crops. maintaining the photosynthetic leaf area especially during the reproductive stage has been shown to impact on the gap between potential and real yield observed, this is due to the incapacity of plants to maintain their green leaf area for longer periods. transcription factors are key proteins involved in the regulation of gene expression and signal transduction networks, regulating different biological processes. several transcription factor families have been associated with leaf senescence in a range of species and their function is crucial for triggering and/or regulating the different aspect of this process. various software tools have been developed with the aim of analyzing datasets to predict hub genes which could regulate a specific pathway or be involved in a particular biological process. weighted gene correlation network analysis is an r package method designed to find clusters of highly correlated genes or metabolites. wgcna calculates a pearson’s correlation matrix for all the genes and then transforms the correlation matrix into an adjacency matrix by raising all values to a soft threshold power; this has the effect of emphasizing the strong correlations and penalizing the weaker correlations. the module eigengene is the first principal component of a given module and can be considered as representative of the module’s gene expression profile. modules often represent specific biological processes; highly connected hub genes within the module are often regulatory genes and represent candidate biomarker. wgcna can be used to construct networks in which each node represents a gene or metabolite and the connecting lines represent correlations between the nodes. wgcna has been widely used as a method to cluster gene expression and metabolite data and identify hub biomarker genes in several non-model plant species such as tomato, brassica, petunia, rice and ficus. biosignature discoverer is a software application devised for identifying molecular signatures in different biological datasets, such as next generation sequencing, microarray data and metabolic profiles, in a statistically robust, computationally efficient, and user-friendly way. a signature is defined as the minimal set of molecular quantities that collectively yield maximal predictive performance. thus, a signature does not contain irrelevant or redundant quantities, given the selected ones. in addition, there may be numerous equivalent signatures that lead to equally accurate predictions, employing different molecular quantities. the core of this application uses a feature selection algorithm that belongs to the class of bayesian network, constraint-based learning and that is able to identify multiple, statistically equivalent signatures. in this work we used these two complementary methods, wgcna and biosd, with the aim of identifying potential transcription factors and metabolites, as putative biomarkers associated with leaf senescence in sunflower. methods plant material and experimental conditions the sunflower leaf senescence assays were conducted under field conditions. the experiment was carried out at the inta balcarce experimental station as previously described. the experiment was sowed during the 2010/ growing season. the sunflower hybrid vdh was sown at a plants/m with three biological replicates, each one consisted of three randomly selected plants from each plot. sunflower assay was conducted under control conditions, without limitations. transcriptomic and metabolomics profiling was performed using the leaf at three different development stages labeled as t-, t- and t-. the samples were immediately frozen in liquid nitrogen upon collection and saved at −80 °c until processing. high quality total rna was isolated from 100 mg of frozen tissue using tripure, according to the manufacturer’s instructions. the genomic dna was eliminated by treatment with dnase i for 20 min at room temperature. the rna concentration was measured using a nanodrop nd- spectrophotometer. the purity and integrity of total rna was determined by 260/280 nm ratio and by nanobioanalyzer rna- analysis. data set transcriptomic profiling was performed using a custom sunflower microarray which has been previously described. background correction was performed using the rma algorithm from the backgroundcorrect function. between arrays normalization was performed using the quantile method from the normalizebetweenarray function. finally, gene expressions were transformed to log scale and information from technical replicates was incorporated by calculating the median parameter. statistically significant probesets where identified using the lmfit and contrasts.fit functions from the r software limma bioconductor package a total of probes were shown to be statistically significant during leaf senescence as the plants develop. the statistical parameters were p-value lower than and fold-change higher or lower than. sunflower transcription factors were identified by comparing approximately tfs sequences from arabidopsis lyrata, arabidopsis thaliana, oryza sativa, populus trichocarpa, vitis vinifera and zea mays available from the plant transcription factor database with sur v database using blast software. a total of tfs were differentially expressed during leaf senescence. metabolic profiling was performed using the gc-tof-ms system. metabolite extraction was performed by promoting the extraction of lipophilic and polar compounds according to recently published protocols adapted for sunflower tissue samples. the chromatograms and spectra were evaluated using chromatof and tagfinder. ion spectra were compared to the golm metabolome database. metabolite levels were normalized to fresh weight using ribitol as the internal control. changes in metabolite levels along leaf development were calculated as fold-change relative to the first sampled time. weighted gene co-expression network analysis wgcna was performed using the wgcna r package as described by langfelder and horvath. the expression values for non-redundant genes and levels of metabolites for the samples were used to construct two separate networks. samples were clustered using the function hclust to check that there were no outliers. the picksoftthreshold function was used to select the soft threshold power used to construct a network based on the criterion of approximate scale-free topology. the power value selected was the lowest power for which the scale-free topology fit index curve flattened out upon reaching a high value. the next step was to transform the adjacency matrix into a topological overlap matrix, which summarizes the degree of shared connections between any two genes. the tom matrix was then converted into a dissimilarity matrix. genes were then clustered using the average linkage hierarchical clustering and the modules were identified in the resulting dendrogram using the dynamic hybrid tree cut method. found modules were trimmed of genes whose correlation with the module eigengene was less than minkmetostay. if p-values of the higher correlations were smaller than those of the native module by the factor reassignthreshold, the gene was re-assigned to the closer module. modules whose eigengenes were highly correlated were merged. this was achieved by clustering module eigengenes using the dissimilarity given by one minus their correlation, cutting the dendrogram at the height mergecutheight and merging all modules on each branch below mergecutheight. this was implemented using the blockwisemodules function in wgcna which performs the network construction and consensus module detection. the following settings were used for the gene network; power =  minmodulesize =  mergecutheight = , maxblocksize =  deepsplit =  reassignthreshold = 1e- and minkmetostay = , networktype = “unsigned”, tomtype = “signed”. the genes were divided into modules; the module sizes ranged from to genes, genes were not assigned to any of the modules. the following settings were used for the metabolite network; power =  minmodulesize =  mergecutheight = , maxblocksize =  deepsplit =  reassignthreshold = 1e- and minkmetostay = , networktype = “unsigned”, tomtype = “signed”. the metabolites were divided into modules; the module sizes ranged from to metabolites, metabolites were not assigned to any of the modules. module membership is a quantitative measure of module membership; this is the correlation of the me and the gene expression profile for each gene in the module. the gene with the highest connectivity in each module was returned using the function choosetophubineachmodule. the networks were exported from wgcna using the function exportnetworktocytoscape. the gene network showed approximate scale free topology, whereas the metabolic network did not satisfy scale free topology: genes scalefreersquared slope − and metabolites scalefreersquared slope −. a total of metabolites with a degree higher than were then selected as metabolite hubs and correlated with the list of tfs statistically significant during senescence using pearson correlation. correlations with p-value < were selected and visualized using cytoscape. biosignature discoverer the statistical analysis pipeline of biosd, a plugin of qiagen’s clc-bio workbench, has been used for identifying gene-expression signatures able to discriminate the age of sunflower plants. particularly, the plugin has been applied on the expression values of transcription factors that were found differentially expressed during leaf senescence, on the metabolites measured in the study, and on the combination of these two data types. additionally, the pipeline has been separately applied on the extended list of non-redundant genes. biosd employs a complex machine-learning approach in order to discover such signatures and for quantifying their actual predictive power. particularly, this tool employs novel, proprietary feature selection algorithms to identify the signatures, inspired by prior work on bayesian-network constraint-based learning. notably, the pipeline is also able to identify multiple equivalent signatures, whose predictive performances are statistically equivalent. to avoid under fitting, the tool automatically tries multiple algorithms for feature selection and prediction modeling, as well as multiple values for tuning them. the tool also employs sophisticated methods, based on an enhanced version of cross-validation, for providing conservative estimations of the performances that can be achieved with the given signatures. the biosd interface allows the user to modulate the trade-off between quality of the results and computational requirements through the ‘tuning effort’ parameter. this parameter regulates the extent to which the tool will attempt tuning the feature selection and prediction modeling algorithms. in all analyses, the ‘tuning effort’ was set to ‘extensive’, in order to favor result accuracy at the expense of higher computational time.  wgcna is a powerful method used to identify clusters of highly correlated genes which are potentially co-regulated. all the connections for the metabolites with an edge weight > produced a network of nodes and edges. the connections between the significant differentially expressed transcription factors were selected from the gene network with an edge weight >, producing a network of nodes and edges. the networks were visualized in cytoscape by edge weight and degree.fig. networks exported from wgcna visualized in cytoscape. network were constructed using transcription factors and metabolites and were visualized in cytoscape with an edge weight higher than. the nodes represent genes and metabolites and the edges represent connections between them. the node size and color is related to the number of connections, large red nodes represent highly connected hub genes and metabolites, small green nodes represent gene with few connections. strong connections are visualized as wider lines a total of transcription factors and metabolites showed the highest number of connections, over 80 % of the maximum number of connections observed, highlighting them as potentially co-regulated hub biomarkers during the senescence process in sunflower. pearson correlation analysis between the hub metabolites and the list of tfs significant during senescence was performed in order to integrate both analyses. the resulting network contained nodes and edges.fig. integrated network of hubs metabolites and transcription factors. pearson correlation analysis of selected hubs metabolites and the list of tfs statistically significant during senescence. correlations with p-value <  were selected and visualized in cytoscape by degree and edge weight the integration results showed that metabolite hubs correspond to sugars and tca cycle metabolites and the nac tfs family, widely reported as senescence associated, showed members as hub genes. on the other hand, biosd allows identification of the most informative genes and metabolites for discriminating between the various stages of the leaf senescence process. table  lists all the transcription factors and metabolites whose values allow the predicting of senescence in sunflower plants. all these quantities have equivalent discriminative power and 100 % stability, which mean that they are consistently selected even when other datasets are included.table list of biosignatures detected in each dataset each signature is composed of a single biomarker, and all signatures are expected to have the same predictive performance table  reports the ‘in-sample’ and ‘out-of-sample’ performances of the selected transcripts and metabolites. the in-sample values quantify the fitness of the predictive models, and poor results would indicate presence of outliers, unmeasured yet needed predictive factors, or poorly modeled trends in the data. the out-of-sample values estimate the performances that are expected if the selected expression and metabolic values are used for predicting the age of new, independent samples of sunflower plants. assuming y and ŷ are the observed and predicted age of n plants. the r metric measures the proportion of y variance explained by the predictive model, r values close to are indicative of good fit, while the mean absolute error and mean squared error statistics quantify predictions’ deviations from the actual age as mae=1n⋅∑absy−y^ and mse=1n⋅∑y−y^ table predicted performances of the selected biosignatures performances are reported in terms of the determination coefficient r mean absolute error and mean squared error. the in-sample performances quantify the fitness of the predictive models on the training data, while the out-of-sample values estimate the expected performance on new data. confidence intervals are calculated using a bootstrap approach in the present application both in- and out-of-sample performances show quite good results, indicating that the predictive models adequately fit the available samples and that the selected transcript factors and metabolites have a high predictive power. mae, a robust estimator of the standard deviation, is around  days for the “out of sample” for both, metabolite and transcript data, meaning that it should be possible to predict the age of a new plant on the basis of transcriptional or metabolic information with an average error of approximately three days. residual vs. predicted value plots did not showed any apparent pattern or trend, indicating that the models have fitted the data adequately. the analysis of the complete list of gene integrated with metabolic information identified twelve distinct signatures each composed of two elements. particularly, each signature contains the gene hean_c_ coupled either with the alanine metabolite or with one of the following eleven genes: hean_c_ hean_c_ hean_c_ hean_c_ hean_c_ hean_c_ hean_c_ hean_s_ hean_s_ hean_s_ hean_s_ these signatures have slightly better performances than the single biomarkers reported in table . these results indicate that the alanine metabolite is interchangeable, in terms of predictive power with each of the other eleven genes. thus, practitioners interested in deploying this signature in practice can freely decide whether to build a predictive model using exclusively transcriptional information or a mixture of transcriptional and metabolic data. this decision can depend on multiple factors, particularly considerations about the cost and technical practicability of measuring the quantities included in the signature. in recent years, many studies have focused on the detection of different genes up or down-regulated along the progression of the leaf senescence process, identifying tfs that could be triggering the process. hence, studying interactions between these tfs under different conditions, for a non-model organism, will open new insights into the details of gene regulation by identifying new interactions and comparing them to those already observed in model organisms. one of the most important aims in new generation biotech crops is to increase the yield and to improve tolerance to different stresses. although sequencing-based approaches are becoming the method of choice for gene expression profiling compared to hybridization-based ones, for non-model species for which a reference sequence is not available, the hybridization approach can still be the more reliable method. this study reports on the use of two complementary software applications with the aim of producing a powerful network and predicting potential biosignatures enabling the identification of important tfs and metabolites involved during the leaf senescence process in sunflower. using wgcna analysis, genes and metabolites are clustered in one or more modules. genes in each module are highly correlated, and consequently they are likely to be co-regulated during the senescence process. biosd identifies which genes are the most informative in order to discriminate between different stages of leaf development. the selected genes, in this case, are not necessarily correlated among themselves. in other words, wgcna attempts to identify the elements involved in the senescence process, and to gain insight of their interplay. on the other hand, biosd tries to retain the minimum number of these elements that are needed for optimal prediction of the plant senescence process. wgcna allowed the detection of metabolite hubs during leaf senescence, most of them correspond to sugars and tca cycle metabolites. sugars are central elements in the sink-source relationships and have been reported as growth and photosynthesis regulators. sunflower is a crop plant with a strong demand for nutrients, especially sugar as a substrate in the oil synthesis during grain filling; this could potentially act as a trigger signal for the senescence process. biosd allowed the detection of one metabolite and ten tfs as potential biomarkers. the selected biosignatures are able to estimate the senescence of a plant with an error of approximately  days. when both tf and metabolite data are simultaneously analyzed, the expected error decreases to  days, indicating that the integration of different omic data allows better predictions of sunflower senescence by biosd analysis. transcription factors are major players and some of them constitute major hubs in signaling pathways. in both analyses, the selected genes belong to transcription factor families widely reported as senescence associated in model species. when we analyze the integration of metabolites and tfs we found that three tfs are detected as biomarkers independently in the two methods. the expression profiles of the selected biomarkers during leaf senescence are displayed in the additional file hean_c_ presents high sequence similarity to aterf which encodes a member of the erf.table candidate biomarkers detected by the two complementary methods, wgcna and biosignature discoverer hean_c_10939 at1g50640 ap2 hean_s_32569 at5g49620 myb hean_c_11118 at1g69490 nac hean_c_10939 at1g50640 ap2 hean_c_11118 at1g69490 nac hean_s_32569 at5g49620 myb biomarkers from wgcna were selected based on their degree and biomarkers from biosignature discoverer correspond to the analysis of the integrated list. biomarkers detected independently in the two methods are in bold transgenic arabidopsis plants with enhanced expression of this tf showed precocious leaf senescence. additionally, hean_c_ has high sequence similarity to atnap/anac, a nac transcription factor which acts downstream of ein and ein genes. ein positively regulates leaf senescence by activating ore and atnap. genetic analysis suggest that both genes act in distinct and overlapping signaling pathways regulating leaf senescence in arabidopsis. hanac a putative sunflower orthologous gene to ore has previously been reported in sunflower as a candidate biomarker for leaf senescence, which reinforces these genes as potential regulators of the senescence process in sunflower. hean_s_ has high sequence similarity to myb. this tf has not yet been directly associated with the senescence process, however, members of this family have been reported as associated to senescence and involved in response to several abiotic stresses. these results suggest that wgcna and the identification of molecular biosignatures are powerful tools for the detection of potential biomarkers. however, it should be noted that these results are obtained on a relatively small experiment. further studies would be advisable in order to better estimate the predictive performance of this approach and to achieve best fit for the list of candidate omics signatures. the implementation of data integration analysis associated to different biotic and abiotic stresses is a powerful tool, especially in a non-model crop such as sunflower for which the complete genome sequence is not yet available. these results open new strategies of analysis to explore and detect potential biomarkers associated with leaf senescence that will be useful for future molecular breeding programs.  the complementary analysis of transcriptomic and metabolomic data with wgcna and biosd emerge as a useful strategy to predict not only the age of the plants but also to detect and identify potential biomarkers associated with leaf senescence. in the case of sunflower, a worldwide oil crop, the implementation of wgcna enabled the construction of a hypothetical network which was used to infer relationships between tfs and metabolites and the identification of hubs as potential candidate biomarkers involved in leaf senescence, whereas biosd selected transcripts and metabolites to discriminate the age of sunflower plants. this methodological approach is a novel strategy for omics data integration, it can be used to elucidate and predict novel networks as well as identifying potential transcripts as putative biomarker hub genes at different developmental stages. additional files additional file 1: cluster dendrogram visualizing the modules from wgcna. a genes, b metabolites were clustered using the average linkage hierarchical clustering and modules identified in the resulting dendrogram by the dynamic hybrid tree cut method. additional file 2: scale-free topologies. a genes, b metabolites. the left panel shows a histogram of network connectivity. the right panel shows a log-log plot of the histogram. the approximate straight line relationship shows approximate scale free topology. additional file 3: list of hubs metabolites and transcription factors. the list shows the identifiers and the normalized expression data of the different samples. additional file 4: biosignature discoverer results on predicting senescence from integrated metabolomics and transcriptomics data in sunflower. additional file 5: expression profiles of the selected biomarkers during leaf senescence. the list shows the identifiers, the expression profile of candidate biomarkers in log scale and the standardization of this expression profile. the figures show the genes and the metabolites in red and blue lines respectively. competing interests the authors declare that they have no competing interests. authors’ contributions pf and sm conceived and designed the work, contributed to the implementation of wgcna for the data, tested software to visualize and understand the data analysis and wrote the first draft of the manuscript, jh implemented wgcna and tested different parameters to optimize the network design, jadr revised the manuscript and statistically supported software analysis, rh and np initiated the project and contributed to the work by the interpretation, discussion of the data and critically revised the manuscript. all authors read and approved the final manuscript. 
the wgcna tool allowed us to design and test a hypothetical network in order to infer relationships across selected transcription factor and metabolite candidate biomarkers involved in leaf senescence, whereas biosignature discoverer selected transcripts and metabolites which discriminate between different ages of sunflower plants. weighted gene correlation network analysis and biosignature discoverer were used to integrate transcriptomic and metabolomic data. wgcnabiosignature discovererleaf senescencetranscriptomicmetabolomicdata integrationsunflowerstatistical methods for omics data integration and analysis smodia 2014heraklion, crete, greece 10- november http://smodia com/issue-copyright-statement© the author 2016 transcription factors are key proteins in the regulation of gene expression, regulating different signaling pathways; their function is crucial for triggering and/or regulating different aspects of the leaf senescence process. wgcna allowed the detection of metabolites and tfs whereas biosd allowed the detection of metabolite and tfs as potential biomarkers. the methodology presented here would help to elucidate and predict novel networks and potential biomarkers of leaf senescence in sunflower. the study of tf interactions and their integration with metabolic profiles under different developmental conditions, especially for a non-model organism such as sunflower, will open new insights into the details of gene regulation of leaf senescence. leaf senescence is a complex mechanism controlled by multiple genetic and environmental variables, which has a strong impact on crop yield. the comparative analysis demonstrated that three transcription factors were detected through both methodologies, highlighting them as potentially robust biomarkers associated with leaf senescence in sunflower. the complementary use of network and biosignature discoverer analysis of transcriptomic and metabolomic data provided a useful tool for identifying candidate genes and metabolites which may have a role during the triggering and development of the leaf senescence process. in recent years, high throughput technologies have led to an increase of datasets from omics disciplines allowing the understanding of the complex regulatory networks associated with biological processes.
 in recent years, high throughput technologies have led to an increase of datasets from omics disciplines allowing the understanding of the complex regulatory networks associated with biological processes. leaf senescence is a complex mechanism controlled by multiple genetic and environmental variables, which has a strong impact on crop yield. transcription factors are key proteins in the regulation of gene expression, regulating different signaling pathways; their function is crucial for triggering and/or regulating different aspects of the leaf senescence process. the study of tf interactions and their integration with metabolic profiles under different developmental conditions, especially for a non-model organism such as sunflower, will open new insights into the details of gene regulation of leaf senescence.  weighted gene correlation network analysis and biosignature discoverer were used to integrate transcriptomic and metabolomic data. wgcna allowed the detection of metabolites and tfs whereas biosd allowed the detection of metabolite and tfs as potential biomarkers. the comparative analysis demonstrated that three transcription factors were detected through both methodologies, highlighting them as potentially robust biomarkers associated with leaf senescence in sunflower.  the complementary use of network and biosignature discoverer analysis of transcriptomic and metabolomic data provided a useful tool for identifying candidate genes and metabolites which may have a role during the triggering and development of the leaf senescence process. the wgcna tool allowed us to design and test a hypothetical network in order to infer relationships across selected transcription factor and metabolite candidate biomarkers involved in leaf senescence, whereas biosignature discoverer selected transcripts and metabolites which discriminate between different ages of sunflower plants. the methodology presented here would help to elucidate and predict novel networks and potential biomarkers of leaf senescence in sunflower. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords wgcnabiosignature discovererleaf senescencetranscriptomicmetabolomicdata integrationsunflowerstatistical methods for omics data integration and analysis smodia 2014heraklion, crete, greece 10- november http://smodia com/issue-copyright-statement© the author 2016 