Metabolite-Centric Reporter Pathway and Tripartite Network Analysis of Arabidopsis Under Cold Stress

The study of plant resistance to cold stress and the metabolic processes underlying its molecular mechanisms benefit crop improvement programs. Here we investigate the effects of cold stress on the metabolic pathways of Arabidopsis when directly inferred at system level from transcriptome data. A metabolite-centric reporter pathway analysis approach enabled the computation of metabolites associated with transcripts at four time points of cold treatment. Tripartite networks of gene-metabolite-pathway connectivity outlined the response of metabolites and pathways to cold stress. Our metabolome-independent analysis revealed stress-associated metabolites in pathway routes of the cold stress response, including amino acid, carbohydrate, lipid, hormone, energy, photosynthesis, and signaling pathways. Cold stress first triggered the mobilization of energy from glycolysis and ethanol degradation to enhance TCA cycle activity via acetyl-CoA. Interestingly, tripartite networks lacked power law behavior and scale free connectivity, favoring modularity. Network rewiring explicitly involved energetics, signal, carbon and redox metabolisms and membrane remodeling.

Introduction

Environmental stresses constrain the growth and production of plants. Abiotic stresses include drought, cold, salt, and osmotic shock and biotic stresses include wounding and pathogen attack. Both kinds of stresses can form reactive oxygen species (ROS) leading to oxidative damage of proteins, nucleic acids, and lipids (Mittler, 2002). Stress tolerance and physiological adaptations mitigate the effects of stress. These processes are mediated by differential gene expression regulating the activity of a number of metabolic pathways and their associated metabolites (Lakshmanan et al., 2015).

Metabolic networks consist of collections of gene-encoded enzyme-catalyzed reactions that convert metabolites that transfer energy and matter along reaction pathways. The exploration of metabolic changes in plants that occur at system level in response to environmental perturbations can be modeled with systems biology approaches. Currently, modeling all components and interactions of a biological system is not realistic due to the large number of parameters that need to be considered (Cakir et al., 2006). One shortcut is to model metabolic activities from transcriptome data to reveal plant responses to abiotic stress (Liu et al., 2016). Many genome-scale models have been developed for model plants and plants of agronomic importance, including Arabidopsis (de Oliveira Dal'Molin et al., 2010) and Chlamydomonas reinhardtii (Dal'Molin et al., 2011), rice (Lakshmanan et al., 2015), and maize (Saha et al., 2011), facilitating global and realistic documentation of perturbation effects in plants.

Omics scale data enables global metabolic overviews. However, understanding cellular physiology requires gaining information of the metabolome and the complexity and high connectivity of metabolic reactions (Nielsen, 2003; Cakir et al., 2006). Metabolic network, sub-network and reporter metabolites have been identified using well-annotated biological data (Patil and Nielsen, 2005). The biological knowledge obtained from these genomic scale approaches is translated into many potential applications, such as engineering agronomic traits and designing stress response gene circuits in crops. Understanding plant response mechanisms requires comprehensive integration of all the outputs of multi-omics approaches (Tatli et al., 2015). In recent years, there have been many profiling studies of the transcriptome, proteome, and metabolome focused on how plants interact with the environment when subjected to stress (Rensink et al., 2005; Budak et al., 2013; Kim et al., 2015; Koç et al., 2015c). However, uncovering genome scale metabolic networks comprising biochemical reactions from qualitative experimental proteomic and metabolomics data is difficult due to the limited sensitivity of current technologies (Yue et al., 2014).

Cold stress affects agricultural productivity in cold regions of the world. Numerous physiological and molecular changes occur during cold acclimation. Physiological changes include chilling injury and growth retardation, and molecular changes include gene and miRNA regulation (Koç et al., 2015a). Here we study the cold stress condition at different time points. We use transcriptome data to gain insight into changes of key metabolites and pathways of Arabidopsis. The Gene Set Enrichment (GSE) method is one of several approaches capable of revealing significant operating pathways from microarray and RNAseq data. These methods identify statistically significant gene sets that are differentially expressed in the gene or pathway lists. While they are often sufficient, there are some biological cases where only few differentially expressed genes are found amid large perturbations. Examples include diagnosis-relapse and tumor tissues (Staal et al., 2003; Yang et al., 2008; Tegge et al., 2012; Akkiprik et al., 2015). Although GSE provides pathway-level views, the general approach still evaluates each gene individually (Tegge et al., 2012). The approach does not describe how genes coordinate their activities within the pathways. Here we take advantage of an algorithm that maps metabolic gene expression by identifying the most significantly changed reporter metabolites using neighboring genes and gene-metabolite association. This reporter metabolite analysis, which has been successfully used in organisms such as yeast, Arabidopsis and rice (Çakir, 2015; Kim et al., 2015; Lakshmanan et al., 2015), elucidates the metabolic pathway regulation of cold stress in Arabidopsis.

Materials and Methods

Data

The gene expression data of Arabidopsis under different abiotic stresses was retrieved from the GEO database of NCBI (ncbi.nlm.nih.gov/geo/). The data comprises cold stress ID series GSE5620 (Control plants) and GSE5621 (Cold stress). These datasets were from experiments carried out using Affymetrix ATH1 microarray chip technology, which included an array of over 22,000 genes (Kilian et al., 2007). In these samples, seeds of A. thaliana Wild Type (col-0) were sown on Magenta boxes containing MS-Agar media and grown as previously described (Kilian et al., 2007). Stress treatment was applied at day 16. In the present study, we investigated the stress response by selecting samples from a transcriptome data set that was sufficiently complete and time-resolved (3, 6, 12, and 24 h of cold treatment). The data correlation of the biological replicates for the cold stress experiment was over 0.95 (Kilian et al., 2007). Initially, data were quantile normalized across all samples, log₂ transformed based on each time point of control, and processed with Principal Component Analysis (PCA) to reduce data dimensionality and help identify outliers. We extracted metabolic genes from AraCyc and then structured a custom-based Excel query to identify metabolic genes of interest within the Affymetrix ATH1 gene information file. The resultant metabolic gene set included 4,730 genes.

Metabolite-Centric Reporter Pathway Analysis (RPA^m)

First, we constructed a draft consensus of an Arabidopsis genome scale network by collecting annotated metabolic genes and their biochemical reactions from AraCyc v.14 (Mueller et al., 2003). Although the AraCyc model, which contains 3,225 reactions, 5,276 enzymes, 2,802 metabolites, and 542 pathways, the draft network was manually re-checked for the presence of reactions in Arabidopsis using KEGG (Kanehisa and Goto, 2000), wikipathways (https://www.wikipathways.org), and especially UniProt (https://www.uniprot.org). In particular, we focused on significantly changed metabolites and pathways for accuracy. Differentially expressed genes between control and stress samples were identified in the Matlab environment using the ttest2 command based on two-sample t-test. Subsequently, P-values were calculated for each gene.

We constructed three files, a P-value data file, a gene metabolite association file, and a pathway metabolite association file for MATLAB. In RPA^m, a metabolite-score must be computed first by the Reporter Metabolic Centric algorithm (Oliveira et al., 2008). The approach depends on the underlying structure of the data and the basic assumptions of the gene expression analysis technology. Since a gene is usually represented by several sets of probes in for example microarray analysis, the significance of differentially expressed genes was assessed by assigning a P-value to each probe set and taking into account the variation within the groups being compared; the P-value for a differentially expressed gene returns the value from the top probe set following established probe set ranks. For reporter metabolite analysis, the genes were linked to the corresponding metabolites via controlled reactions with an algorithm that takes as input P-values of gene expression data and the topology of bio-molecular interaction networks represented as a graph (Ideker et al., 2002). P-values were then assigned to relevant enzymes and converted to Z scores with an inverse cumulative distribution function (θ⁻¹). The use of this distribution function converts uniformly distributed P-values for each probe set into normally distributed standard Z score random variables. For isozymes, the lowest P-values of isozymes were selected.

Once each enzyme was scored with a Z-value, a Z_metabolite score was calculated for each metabolite as an aggregate Z score of its k neighboring genes in the metabolic network scale. In this step, the Z scores of genes connected to a metabolite in the metabolic model were averaged to obtain a final Z score for the metabolite (Cakir et al., 2006).

$\begin{array}{l}{Z}_{metabolite}=\frac{1}{\sqrt{k}}{\displaystyle \sum _{i=1}^k}{Z}_i\end{array}$

Subnetworks of all sizes are comparable under this scoring system. If the Z_gene scores are independently drawn from a standard normal distribution, the aggregate Z_metabolite score will also be distributed according to a standard normal distribution, independent of k (Ideker et al., 2002). A high Z_metabolite score indicates a biologically active subnetwork.

In order to evaluate the significance of metabolite scores, the score of each metabolite was corrected by calculating the mean (μ_k) and standard deviation (σ_k) of a background Z score distribution resulting from 10,000 rounds of random sampling of k metabolites of metabolic networks (Väremo et al., 2013). The number of random samples of the distribution was determined by checking the sensitivity of the resulting background scores relative to the number of random samples. The Z score can be considered a parameter of statistical significance (typically Z = 1.96 corresponds to P = 0.05). Furthermore, the reporter metabolite analysis is similar to the method of Stouffer. However, the statistical parameters of the gene set were corrected for the background distribution prior to the estimation of significance (for details see Oliveira et al., 2008; Väremo et al., 2013).

$\begin{array}{l}{Z^{\prime }}_{metabolite}=\frac{Z_{metabolite}{-u}_k}{{\sigma }_k}\end{array}$

Here, the analysis of a “corrected” Z-value, ${Z^{\prime }}_{metabolite}$ is the reporter metabolite analysis (Patil and Nielsen, 2005). The final Z score was then transformed into a P value using the normal cumulative distribution and assigned to each metabolite (Figure 1) in order to detect those that are significant at P ≤ 0.05 level, which represent reporter metabolites.

FIGURE 1

Figure 1. Strategy and schematic representation of the metabolic-centric reporter pathway analysis (RPA^m) and its visualization with a tripartite network. (A) The flow diagram describes the computational strategy of the reporter pathway analysis and the integration of transcriptome data. P-values for differentially expressed genes are mapped onto metabolites associated genes. Each enzyme is assigned a score based on the average of the P-values of the probe sets representing the corresponding gene. Minimum P-values are chosen for a reaction catalyzed by an enzyme complex or a set of isoenzymes. P-values are then converted into Z scores and corrected with a background Z score distribution. The resulting scores are linked to related pathways, which are represented with a tripartite network. (B) Open and closed tripartite networks always connect nodes of different sets, with or without restrictions in how sets are connected with each other, respectively. (C) A pathway (yellow octagon) defines a set of reactions (red circles labeled with letter R) connecting metabolites (blue squares labeled with M). Genes (purple triangles labeled with G) can be associated to specific metabolites (gray arrows) when these are collectively significantly expressed at some P-value (labeled with p) over some threshold. (D) The resulting open tripartite network of genes, metabolites and pathways derived from the associations of (C).

Note that we removed some metabolites from the reporter metabolite analysis, including proton, water, oxygen molecule, NADH, NAD⁺, NADPH, NADP⁺, ATP, ADP, AMP, GTP, GDP, UDP, CoA, FAD, diphosphate, carbon dioxide, carbon monoxide, phosphate, ammonia, hydrogen peroxide, oxidized electron acceptor, and reduced electron acceptor. Their role in many reactions bias metabolite measurements (Kim et al., 2015).

Each metabolic pathway is similarly scored based on its metabolites. Scoring of significant pathways is based on P-values and Z scores of the corresponding metabolites following 10,000 iterations of the sampling algorithm.

$\begin{array}{l}{Z}_{pathway}=\frac{1}{\sqrt{n}}{\displaystyle \sum _{metabolite=1}^n}{Z}_{metabolite}\end{array}$

n is the number of metabolites associated with the pathway. A significance level of P ≤ 0.05 was used as a threshold to identify reporter pathways. Note that the Z_pathway score is corrected for the background distribution as calculated in ${Z^{\prime }}_{metabolite}$. Pathways with significant ${Z^{\prime }}_{metabolite}$ scores are called reporter pathways (Çakir, 2015).

We used the genome scale connection of genes and metabolites that are embedded in pathways to construct a tripartite network. They are represented in the graph of Figure 1, which illustrates these connections. This atypical network has vertices describing significantly changed genes, metabolites, and pathways and edges representing the connectivity of significantly changed vertices. The network was used to study time resolved responses of Arabidopsis to cold stress.

Network Analysis

The biological networks of genes, metabolites and pathways were visualized with Cytoscape 3.2 (Shannon et al., 2003). Clustering analysis was performed with the ClusterViz plugin using the MCODE algorithm, detecting clusters with default cutoffs. Scale free network behavior was examined using log-log plots, i.e., log of P(k) vs. log of k. The exponent γ of the power law [P(k) ~ k^−γ] and the determination coefficient (R²) were derived from linear regression models. γ is the negative of the slope of the log linear model. R² explains how much of the data fits the linear model. The scale free behavior of the network should be considered strongly supported when both γ and R² are high. The Kolmogorov-Smirnov (KS) fit statistic compares the fitted distribution with the input degree vector and the KS P-value (Newman, 2005; Clauset et al., 2009). The null hypothesis corresponds to data drawn from the power law distribution. We determined the maximum log likelihood of the fitted parameters. We also determined if the power law fit pattern was continuous. We created reference networks with Barabási methods using the R's igraph package (Csardi and Nepusz, 2006) to simulate a power law for the statistical analysis of networks. Modularity analyses were performed by using Fast Greedy Community (FGC) and Newman-Girvan (NG) algorithms (Clauset et al., 2004; Newman and Girvan, 2004). The NG algorithm iteratively computes the “betweenness” of edges, i.e., the number of shortest paths that run through an edge; the network is decomposed into communities by systematically removing edges and recalculating betweenness measures after each removal (Newman and Girvan, 2004). The FGC hierarchical agglomeration algorithm runs in time O(md log n) ~ O(n log² n), with m, n, and d representing edges, vertices and the depth of the dendrogram describing the structure of the community, respectively (Clauset et al., 2004). In addition, time resolved tripartite networks were overlaid and visualized using CompNet software (Kuntal et al., 2016).

Results and Discussion

Reporter Pathway and Tripartite Network Analysis

We explored the effect of cold stress on the metabolic networks of Arabidopsis with the RPA^m approach by: (i) analyzing differentially expressed genes in two transcriptome datasets (GSE5620 and GSE5621), (ii) associating these genes to metabolites, metabolic reactions, and metabolic pathways, and (iii) visualizing genes, metabolites and pathways using a tripartite graph representation (Figure 1A). The central data that is mined is the transcriptome. We assume that the differential expression of metabolic genes at a biological endpoint indicates their active involvement in core metabolic functions and/or their regulation.

Tripartite graphs are a special case of k-partite graphs with k = 3. k-partite graphs have nodes (graph vertices) that can be divided (partitioned or colored) into k disjoint sets (partitions or colors) and connections (graph edges) that always connect nodes belonging to different sets. Closed tripartite graphs do not impose restrictions on the tripartite structure of connected nodes (all sets can connect to each other) (Figure 1B). In contrast, open tripartite graphs do not allow a circular connecting structure. A gene-metabolite-pathway tripartite network describes the association of genes to metabolites at some minimum p-value and the association of metabolites to metabolic pathways through reactions (Figure 1C). We describe these associations with an “open” tripartite graph representation of genes (G_i), metabolites (M_i), and pathways (P_i), in which allowed connections are only between G_i and M_i nodes and between M_i and P_i nodes (Figure 1D). One benefit of open tripartite graphs is that they can be decomposed into one-mode and two-mode (bipartite) graph projections to improve visualization. For example, our tripartite graph can be transformed into a more coarse-grained bipartite graph of G_i and M_i nodes by excluding pathway P_i node information, i.e., how metabolites connect many different operating pathways. Alternatively, exclusion of gene G_i information produces a bipartite graph of M_i and P_i nodes.

An initial exploration of the differential expression of over 22,000 genes with PCA analysis showed samples clustering into groups and a time resolved response to cold stress occurring at 4 time points, 3, 6, 12, and 24 h of cold treatment (Figure S1). We also conducted PCA analysis using only metabolic genes. The projections of metabolic genes were consistently similar to those of all genes. These findings supported the use of the metabolic gene subset to study metabolic mechanisms behind stress perturbations in Arabidopsis.

Note that a reporter metabolite approach of our metabolite-centric RPA^m better represents the effects of environmental perturbations at network-level since it considers all connecting reactions that consume or produce the metabolites of the pathway (Çakir, 2015). Notably, the RPA^m algorithm does not require a priori knowledge of changes that are significant at the level of each node (e.g., the transcript of a gene) (Oliveira et al., 2008). In contrast, a reaction-centric pathway analysis takes into account only significantly changed gene transcripts which are involved in associated reactions. Note that genes that are not differentially expressed at significant levels can operate in coordination with others (Mootha et al., 2003). Thus, it is important to consider both the magnitude of the significance of their differential expression and their connectivity in the network (Çakir, 2015). Using the reporter metabolite algorithm, we detected reporter metabolites in Arabidopsis and found that 64, 68, 71, and 102 reporter metabolites were significantly (P ≤ 0.05) regulated in the 3, 6, 12, and 24 h cold acclimated plants, respectively (Table 1, Supplementary File 1). In turn, 79, 79, 103, and 94 reporter pathways were significantly regulated in the corresponding cold treatments (Table 2, Supplementary File 1). We identified several metabolic pathways regulated under cold stress, including glutathione redox reactions, sucrose biosynthesis, and the biosynthesis and degradation of starch, mannitol, xylan, sucrose, triacylglycerol, cellulose, 3-phosphoinositide, L-ascorbate, and trehalose. In contrast, GSE analysis led to the identification of 45, 58, 89, and 47 pathways in the 3, 6, 12, and 24 h cold acclimated plants, respectively (data not shown). These pathways were mainly related to redox, sugar, lipid, amino acid metabolisms, which are similar to those of RPA^m. Remarkably, GSE analysis did not capture the TCA cycle at any time point, which is a known effect of cold stress (Kaplan et al., 2004).

TABLE 1

Table 1. Top 20 reporter metabolites identified for each time point of cold stress are labeled with matching superscript letters (a, b, c, and d for 3h, 6h, 12h and 24h, respectively) and listed with associated P-values.

TABLE 2

Table 2. Top 20 reporter pathways identified for each time point of cold stress are labeled with matching superscript letters (a, b, c and d for 3h, 6h, 12h and 24h, respectively) and listed with associated P-values.

In our studies, regulated pathways are made explicit through the gene, metabolite and pathway associations of the tripartite gene-metabolite-pathway network that is illustrated in Figure 2. The network and its projections consists of three hierarchical sets of entities (nodes) describing genes (triangles), metabolites (squares), and pathways (octagons), which are connected to each other when significant changes are detected. The size of vertices was scaled to vertex connectivity and colored by P-value within the scale 0–0.05 making hub-like behavior and the metabolite-centric approach explicit in the tripartite network.

FIGURE 2

Figure 2. A gene-metabolite-pathway tripartite network of significant changes (P ≤ 0.05) for 3 h (A), 6 h (B), 12 h (C), and 24 h (D) cold acclimated Arabidopsis. Triangles represent genes, rectangles represent metabolites and octagons represent pathways. Red color range shows most significant genes, metabolites and pathways. Node sizes were scaled to number of neighbors. The networks were visualized by Force directed layout. The white squares represent snapshot of the network. Cytoscape files for network visualization are deposited in https://github.com/gcalab/files.

RPA^m led to identification of the most significant effects of cold stress. For example, metabolisms associated with energy generation were significantly detected. The TCA cycle and ethanol degradation were significant pathways which exhibited high network connectivity (see below). Moreover, the glycolytic pathway that is linked to perturbation of energy was significantly identified. On the other hand, ethanol (P = 0.01) was a significant metabolite which leads to membrane lipid composition for stabilization (Kaplan et al., 2004). The RPA^m identified several common pathways that were significant for all time points, including glutathione redox reactions, sucrose biosynthesis, cellulose biosynthesis, xylan biosynthesis, and D-myo-inositol (1,4,5)-trisphosphate biosynthesis. All of these pathways suggest stress related changes occur in cell membrane lipid composition at all time points through signal metabolism. Moreover, the regulatory role of sucrose in cold acclimation was also reported (Rekarte-Cowie et al., 2008). The following sections describe in detail the most significant effects of cold stress on metabolic networks.

Carbohydrate and Amino Acid Metabolic Activity Under Cold Stress

Plants subjected to cold stress often show water stress symptoms resulting from cold induced water stress (Singh et al., 2012). Adverse environmental stresses including cold trigger rapid regulatory responses of carbohydrate metabolism. For example, starch and sucrose are osmolytes that modulate regulation of water content and cellular dehydration during cold-induced ice formation. While cold stress induces the accumulation of proline, which is a well-known osmo-protectant, we found it was not significantly regulated in our study (P = 0.1 for 24 h cold stress). In contrast, sucrose biosynthesis was significantly regulated at all time points, indicating a key role of sucrose in cold perturbation (Table 2).

Carbohydrate metabolic changes were also found associated with the α-D-glucose-1- phosphate (Glc-1-P) metabolite (P = 0.02 for 24 h cold acclimation; Figures S2A–C), which is an important intermediate in several major carbon fluxes, including those involved in starch, sucrose, and cellulose biosynthesis (Figure 2D). Several carbohydrate degradation pathways were also regulated at different cold acclimation time points, including D-mannose degradation, (1,4)-beta-xylan degradation, sucrose degradation II, galactose degradation I (Leloir pathway), galactose degradation III, and starch degradation II. Taking into account the regulation of these pathways, our results suggest that the breakdown of carbohydrate reserves supported the energetics necessary to survive under cold stress conditions. In previous studies, the whole-plant metabolome of Arabidopsis indicated that cold stress increased many metabolites, including amines (asparagine, aspartate, glycine, serine, tryptophan, glutamine and proline, glycine betaine), organic acids (ascorbate, gluconate, malate, and α-ketoglutarate) and carbohydrates (sucrose, maltose, inositol, glucose, fructose, and trehalose) (Cook et al., 2004; Kaplan et al., 2004). In agreement, our method captured glucose, sucrose, maltose, amylopectin, myo-inositol, ascorbate, glycine betaine, L-allo-threonine, L-cysteine, valine, prephenate (precursor of phenylalanine) and other metabolites as significantly regulated (Supplementary File 1, Figure S2). Expectedly, several amino acid metabolic networks were significant in this study, including L-glutamine biosynthesis II, phenylalanine biosynthesis II, tyrosine biosynthesis II, L- glutamine biosynthesis II, and cysteine biosynthesis. Results are consistent with the hypothesis that amino acids and many carbohydrate metabolites were reconfigured under cold stress as a consequence of changes in metabolic activities. The significantly changed metabolites were related to various metabolic pathways including amino acid metabolism, carbohydrate metabolism, cell wall metabolism, fatty acid metabolism, redox metabolism, secondary metabolism, and signal transduction (Figure S2). These cold-induced metabolites may play key roles in signal transduction for gene expression. Therefore, the metabolite-centric approach may reveal whether pathways are hierarchically or metabolically regulated in conjunction with pathway cross-talk (Cakir et al., 2006; Çakir, 2015).

Energy, Protective, and Signal Metabolic Activity Under Cold Stress

We identified key pathways that consecutively build protective and signal metabolisms under cold stress. These pathways included callose, ascorbate, D-myo-inositol-5-phosphate, 3-phosphoinositide, xylogalacturonan, xylan, cellulose, homogalacturonan, indole-3-acetyl-amide conjugate biosynthesis, anthocyanidin modification, and phosphatidate metabolism (as a signaling molecule). The 1-phosphatidyl-1D-myo-inositol 4,5-bisphosphate (PIP2) metabolite (P = 0.005 for 3 h) is an important signal molecule of the D-myo-inositol-5-phosphate biosynthesis pathway, which is involved in abiotic stress, pollen tube growth, ion channel activity, and guard cell movement (Kost et al., 1999; DeWald et al., 2001; Jung et al., 2002; Liu et al., 2005). Our reporter pathway analysis showed that the gluconeogenesis pathway also scored significantly. This pathway makes glucose from lipids and proteins. Thus, it seems critical to facilitate glucose storage during cold stress when breeding for stress resistance. In addition to providing energy supplements, glucose also plays a key role in membrane stabilization. In addition, the analysis showed glycolysis I, glycolysis II, glycolysis IV, 2-oxoglutarate decarboxylation to succinyl-CoA, pyruvate decarboxylation to acetyl CoA, TCA cycle II (plants and fungi) and, TCA cycle V pathways were significant pathways. In previous studies, TCA cycle intermediates were increased in response to temperature and drought (Kaplan et al., 2004; Urano et al., 2009). When we explored the expression of TCA cycle genes, they were perturbed during early time points of cold stress and then showed only slight increases (Figure 3). Since ATP production needed to maintain adequate ATP/ADP ratios is not decreased in cold stress or cold acclimated Arabidopsis (Strand et al., 1997; Hurry et al., 2000), our finding supports the notion that cold stress first induces glycolysis and then adjusts the TCA cycle for ATP production.

FIGURE 3

Figure 3. The heatmap diagram of expression profile of TCA cycle genes. Red color hues indicates gene up-regulation while green color hues indicates down-regulation. The green-to-red scale below the heatmap describes expression values.

In addition, low-temperature stress causes an accumulation of phosphorylated metabolites (Savitch et al., 2001). This is consistent with ATP production from glycolysis, sucrose degradation, and D-mannose degradation pathways, which support the glycolysis pathway, being also affected in our study. The increase in glycolytic flux requires a coordinated increase in the flux of the pentose phosphate pathway since they are both parallel pathways of the metabolic network that produces ribose 5-phosphate (P = 0.017 for 12 h), a precursor for the synthesis of nucleotides. Significantly changed activities of the pentose phosphate (oxidative branch) I and pentose phosphate (non-oxidative branch) pathways are also the source of carbon skeletons for the synthesis of aromatic amino acids, phenylpropanoids, and their derivatives (Herrmann and Weaver, 1999). AT1G20950, AT4G04040, AT1G76550, and AT1G12000 encode fructose-6- phosphate 1-phosphotransferase (EC 2.7.1.90). This enzyme catalyzes the initial reaction of the glycolysis IV (plant cytosol) pathway, which converts β-D-fructofuranose 6-phosphate into fructose-1,6-bisphosphate. In the last reaction of glycolysis IV, AT3G22960, AT5G52920 and, AT3G55650 (pyruvate kinase) encode for enzymes that catalyze the conversion of phosphoenolpyruvate into pyruvate. AT2G29350 encodes the alcohol dehydrogenase (EC 1.1.1.1) enzyme that converts ethanol into acetaldehyde. In turn, AT3G48000, AT1G44170, and AT3G24503 (aldehyde dehydrogenase; EC 1.2.1.3) catalyze the conversion of acetaldehyde into acetate, and AT5G36880, AT3G16910 and, AT1G55320 (acetate-CoA ligase; EC 6.2.1.1) catalyze the conversion of acetate into acetyl-CoA. Previous studies have reported a pyruvate dehydrogenase bypass that produces ethanol to support the TCA cycle and lipid biosynthesis (Mellema et al., 2002; Gass et al., 2005; Yue et al., 2014). Additionally, the antimycin A inhibition of the mitochondrial electron transport chain triggered the immediate production of ethanol and a fast re-arrangement of metabolic pathways. Ethanol availability and its fermentation compensated for reduced ATP production (Obermeyer et al., 2013). Our study shows that glycolysis and ethanol fermentation (P = 0.01 for 6 h, Table 1, Figure S2B) under cold stress diverted energy into the TCA cycle, turning it into an alternative energy producing pathway of Arabidopsis (Figure 4).

FIGURE 4

Figure 4. The TCA and ethanol degradation II pathways of Arabidopsis. The Aracyc diagram illustrates the enzymes with corresponding EC numbers and substrates that makeup these metabolic routes.

Interestingly, the Rubisco shunt was also significantly changed in Arabidopsis. In general, the glycolysis I (from glucose 6-phosphate) and glycolysis IV (plant cytosol) pathways facilitate the conversion of carbohydrates to oil. Pyruvate, which is the end product of glycolysis, is further decarboxylated into acetyl-CoA, which is a precursor for fatty acid and oil biosynthesis (Neuhaus and Emes, 2000). However, this process causes the loss of one-third of carbon as CO₂ (Schwender et al., 2004). This is not metabolically advantageous. When compared to glycolysis, the Rubisco shunt increases the efficiency of carbon conversion, producing 20% more acetyl-CoA with 40% less carbon loss. Rubisco operates a metabolic route between carbohydrate and lipid metabolism in coordination with Acetyl-CoA derivatives and CO₂ producing pathways (2-oxoglutarate decarboxylation to succinyl-CoA and pyruvate decarboxylation to acetyl-CoA), which were reprogrammed for cold acclimation. This occurs in three stages: (1) the conversion of hexose phosphates to ribulose-1,5-bisphosphate by phosphoribulokinase (the non-oxidative reaction); (2) the conversion of ribulose-1,5-bisphosphate and CO₂ (mostly produced by PDH3) to PGA by Rubisco (ribulose 1,5-bisphosphate carboxylase/oxygenase); and (3) the metabolism of PGA to pyruvate and then to fatty acids (Figure 5). We hypothesized that the Rubisco shunt balances carbohydrate and oil production in order to make efficient use of the carbon pool under cold stress.

FIGURE 5

Figure 5. Representation of the Rubisco shunt in Arabidopsis with consecutive enzymes. The pathway was taken from Aracyc.

Lipid Metabolic Activity Under Cold Stress

Several fatty acid pathways were impacted by cold stress, including very long chain fatty acid biosynthesis I and II, choline biosynthesis III, triacylglycerol biosynthesis, diacylglycerol biosynthesis (polyunsaturated fatty acids), poly-hydroxy fatty acids biosynthesis, phospholipid remodeling, phospholipid desaturation, glycolipid desaturation, phospholipases, and phosphatidylcholine biosynthesis I and II. Cold, heat, and osmotic stress change the physical properties of lipids in cell membranes (Los and Murata, 2004). It is well-documented that the cell membranes are a primary site of cold-induced damage (Thomashow, 2001; Matteucci et al., 2011). In these studies, cold acclimation activates mechanisms that protect the fluidity of membranes by ensuring that associated enzymes have optimal activity. Cold responsive plant species can change membrane lipid fluidity by increasing the levels of unsaturated fatty acids (Matteucci et al., 2011). These molecules are crucial to cold tolerance and photosynthesis (Meyerowitz and Somerville, 1994). In one study, profiling membrane lipids of Arabidopsis subjected to 3 days of cold acclimation revealed an increase in desaturation in all measured phospholipids (Welti et al., 2002). Similarly, low temperature resulted in an increase in the degree of fatty acid unsaturation in chickpea (Bakht et al., 2006, p. 2). In accordance with these lipid profiling studies, our reporter pathway analysis showed changes of fatty acid metabolism necessary to alter membrane composition in Arabidopsis under cold stress.

Cell Wall Metabolism Under Cold Stress

Exposure to chilling temperature causes alterations in the composition of the cell wall and in the activities of cell wall-modifying enzymes (Le Gall et al., 2015). For example, cell wall strength was increased in grape suspension culture and apple under cold acclimation (Rajashekar and Lafta, 1996). We identified several pathways related to cell wall composition as significantly changed. Among them, cellulose biosynthesis, xylan biosynthesis, homogalacturonan biosynthesis, xylogalacturonan biosynthesis, xyloglucan biosynthesis, cuticular wax biosynthesis, cutin biosynthesis were significant. Xylan biosynthesis comprises the synthesis of several components of the cell wall such as xylan and arabinoxylan. In addition, homogalacturonan (P < 0.05 for 3, 6, and 12 h) accounts for 60% of the total pectin in primary cell walls of plants (Mohnen et al., 1996). Its content was reported to increase in oilseed rape leaves in response to cold exposure (Solecka et al., 2008). Thus, pectin (Table 1, Figures S2A,D), which was a significant metabolite of our study, appears a key cell wall component affected by the plant response to cold stress.

Methanol Production Under Cold Stress

We found that pectin-interacting methanol was a crucial metabolite in the 3 and 24 h cold treatments (Figures S2A,D). Most plants produce and extrude methanol, especially because of pectin demethylation during the early stages of leaf expansion (Fall and Benson, 1996). The degree of methylation plays a key role in properties of the cell wall and cell adhesion (Liu et al., 2014). PEG treatment decreased the degree of methylation of pectin in root tips of Phaseolus vulgaris. This metabolite is incorporated into the methyl groups of numerous metabolites in higher plants, including methionine, serine and phosphatidylcholine. Additionally, methanol induces the de novo synthesis of methyl-β-d-glucopyranoside (Gout et al., 2000). In this regard, overwintering leaves of alpine Rosaceae plant Dryas octopetala produced the highest levels of methyl-β-d-glucopyranoside (Aubert, 2004). Therefore, methanol appears involved in cold stress and its role merits further study.

Redox Metabolism Under Cold Stress

Our reporter pathway analysis identified several redox metabolisms, including L-ascorbate biosynthesis I (L-galactose pathway), L-ascorbate biosynthesis V, ascorbate glutathione cycle, glutathione redox reactions I and II, gamma-glutamyl cycle, formaldehyde oxidation II (glutathione-dependent), sulfate reduction II, formaldehyde oxidation II and thioredoxin pathway in cold stress (Supplementary File 1, Figure 2). L-ascorbic acid, popularly known as vitamin C, is one of the most abundant water-soluble antioxidants. Its ability to scavenge free radicals and ROS produced by drought, salt, and high or low temperature stress in plants demonstrates its important role in plant growth, development and abiotic stress tolerance (Zhang et al., 2015). Additionally, we found that glutathione (P < 0.05 for 3, 12, and 24 h) was a key metabolite in redox metabolism (Figure S2). When coupled with the ascorbate/glutathione cycle and the redox system, ascorbate can effectively regulate cellular hydrogen peroxide levels in plants (Ishikawa et al., 2008). Most genes encoding enzymes in this pathway, including GDP-mannose pyrophosphorylase (GMP) (Badejo et al., 2008), GDP-mannose-3′,5′-epimerase (Zhang et al., 2011), GDP-L-galactose phosphorylase (Bulley et al., 2012), L-galactose-1-phosphate phosphatase (Zhou et al., 2012), L-galactose dehydrogenase (Gatzek et al., 2002), and L-galactono-1,4-lactone dehydrogenase (Liu et al., 2011), have been shown to modulate ascorbic acid levels in plants. Ascorbate directly counteracts ROS and participates in the regeneration of α-tocopherol, which scavenges both ROS and lipid peroxyl radicals and zeaxanthin, which plays roles in the photoprotective xanthophyll cycle (Bielen et al., 2013). The glutathione redox pathway maintains the reduction of various molecules such as peroxides, glutaredoxin, or protein-disulfides that are formed from two molecules of glutathione (GSH) and are catalyzed into the oxidized form of glutathione disulfide (Milla et al., 2003). A recent study showed that the accumulation of both reduced and oxidized ascorbate coupled with oxidized and reduced glutathione were induced by cold and salt stresses (Wang et al., 2016). We also identified zeaxanthin (P = 0.04) and lutein metabolism (carotenoid metabolisms; P = 0.02), functioning as antioxidants. Similarly, Coffea canephora showed higher contents of the zeaxanthin and lutein in cold perturbation (Partelli et al., 2009). Aerobic respiration III (alternative oxidase pathway) was found regulated in cold stress. This pathway is not associated with ATP production. Nevertheless, alternative oxidase activity in this pathway decreases the formation of ROS by hindering the over-reduction of the ETC (Møller, 2001). Therefore, pathways linked to redox metabolism are expected to be active under cold stress for rebalancing the impaired redox state of the cellular machinery.

Photosynthesis Pathway Under Cold Stress

Photosynthesis is undoubtedly a crucial sensor of stress in plants. Photosystem I (PS I) and II (PS II) and Rubisco are considered to act as primary stress sensors of the chloroplast (Koç et al., 2015c). Stress perturbation of these sensors generate signals such as ROS production, change of sugar levels, redox reactions of the photosynthetic electron transport system, and energy imbalance, which cause metabolic/molecular reconfiguration of stress adaptation (Biswal et al., 2011; Koç et al., 2015b). Abiotic stresses can also affect the rates of primary photochemical reactions and reduce the activities of the enzymes of the Calvin cycle (Biswal et al., 2011). We identified that the Calvin-Benson-Bassham cycle (Calvin cycle) and photosynthesis light reactions were significantly regulated in Arabidopsis during cold acclimation. The Calvin cycle determines the net fixation of CO₂ in plants. Collectively, the reporter metabolites of the Calvin cycle and Rubisco shunt indicate that cold reduces the efficiency of carbon fixation. However, Rubisco participates in a previously described metabolic route between carbohydrate and lipid metabolism mediated by increases of acetyl-CoA levels. This strategy developed by plants can be thought as efficient metabolic flux to carbon storage (Schwender et al., 2004). It is well known that plants reduce the oxidative damage caused by low temperatures by down-regulating light reactions known to increase the oxidative damage. In this respect, we also identified that the chlorophyll a degradation I pathway was regulated under cold acclimation. Several studies reported that cold stress causes degradation of chlorophyll (Tewari and Tripathy, 1998; Liu et al., 2012; Singh et al., 2012). To this end, the chlorophyll biosynthetic pathway, which plays significant roles in photosynthesis and in avoiding the accumulation of phototoxic chlorophyll intermediates, must be tightly regulated (op den Camp et al., 2003; Rodríguez et al., 2013).

Hormone Metabolic Activity Under Cold Stress

Increasing evidence indicates that synergistic or antagonistic hormone cross-talk play key roles during abiotic stress (Santner and Estelle, 2009). Plant hormones often regulate metabolism by rapidly modulating transcriptional factor activity (Peleg and Blumwald, 2011). Our reporter pathway analysis identified phytohormone metabolism of auxin (IAA) biosynthesis VII, indole- 3-acetyl-amide conjugate biosynthesis, ethylene biosynthesis I, cytokinins 7-N-glucoside biosynthesis, and cytokinin-O-glucosides biosynthesis as significantly changed during cold acclimation (Figure 2). Auxin is involved in many physiological processes, including plant growth, development and the abiotic stress response (Jain and Khurana, 2009). In Arabidopsis, cold stress is known to affect the auxin transport pathway by inhibiting various intracellular proteins, including auxin efflux carriers (Shibasaki et al., 2009). Similarly, the expression profile of the ARF and Aux/IAA gene family members of Arabidopsis changed during cold acclimation (Hannah et al., 2005). Transcriptome analysis in rice showed that over 154 genes were induced and 50 were repressed by auxin under desiccation, salt, and cold stress (Jain and Khurana, 2009). Most IAAs are considered to be conjugated to a number of amino acids (Staswick et al., 2005). Some of these amino acid modifications can be reversed through the activity of amido hydrolases (Davies et al., 1999), suggesting that compounds such as IAA-alanine and IAA-leucine are storage forms of auxin that can be further metabolized to feed into the free auxin pool (Staswick et al., 2005). Indole-3-butyrate (IBA) was also a significantly changed metabolite (P = 0.01 for 24 h cold stress), suggesting it functions as an endogenous IAA precursor (Ljung, 2013). Some investigations of IBA biosynthesis show that its level is likely regulated by plant hormones and various stresses (Ludwig-Müller et al., 1995, 1997). Another phytohormone, ethylene, regulates germination, fruit ripening, organ abscission, and pathogen, senescence, and biotic/abiotic stresses (Chen et al., 2005). Increased ethylene levels confer enhanced freezing tolerance in Arabidopsis by decreasing freezing tolerance and inhibiting ethylene biosynthesis/signaling (Shi et al., 2012). Ethylene synthesis is sustained by inhibiting 1-aminocyclopropane-1-carboxylate oxidase (Vriezen et al., 1999). Since AT1G12010 and AT1G05010 (1-aminocyclopropane-1-carboxylate oxidase), which convert 1-aminocyclopropane-1-carboxylate to ethylene, were down-regulated, we find evidence that the ethylene pathway was perturbed under cold stress (Figure S3).

The signaling and response to cold stress has been associated with cytokinin signal transduction and a subset of a two-component signaling (TCS) system (Jeon et al., 2010; Jeon and Kim, 2012). The ARABIDOPSIS RESPONSE REGULATOR1 mediates signaling through the activity of AHP2, AHP3, or AHP5 to induce a subset of type-A ARRs (Jeon et al., 2010; Jeon and Kim, 2012). Together with cytokinin, this regulates tolerance and the cold stress response. However, the exact role of cytokinins on the cold stress response is unclear and needs to be further explored. In addition, our reporter metabolite analysis also identified some other phytohormone-related pathways of brassinosteroid (superpathway of C28 brassinosteroid biosynthesis), jasmonate (jasmonic acid biosynthesis), and gibberellin (GA12 biosynthesis) hormones (Table 2, Supplementary File 1). Our findings suggest plant hormones may coordinate responses to increase the plant adaptability to cold stress.

The Structure of Tripartite Networks

Stress responses have been widely studied in plants. Here we examine stress-induced patterns by visualizing changes in metabolites and pathways using reporter metabolite analysis in absence of metabolome data. The tripartite networks quantify the dynamic aspects of stress-induced change by making connectivity and modularity explicit in the network structures that describe the time resolved cold response. All tripartite networks resulted in disconnected undirected graphs. Only few small intra-connected components were detected to be isolated from the large connected component. As time progresses, the appearance of nodes and links emerge in the tripartite networks with increasing cold acclimation. The emergent networks had 218, 246, 306, and 320 nodes with a network density (actual/possible number of edges) of 0.023, 0.021, 0.026, and 0.028 at 3, 6, 12, and 24 h of cold acclimation, respectively (Table 3). An analysis of network topology showed that networks had diameters of 12, 9, 10, and 9, respectively, and exhibited average numbers of neighboring nodes of 5.1, 4.3, 3.9, and 3.7, respectively (Table 3). These results suggest that tripartite networks become more compact and gain decentralized infrastructure with increasing cold acclimation.

TABLE 3

Table 3. Graph properties of the tripartite networks at different times of cold acclimation.

Clustering analysis revealed that networks were well-structured. We see more connected rewiring patterns in response to stress over time (Figures 2, 6). Remarkably, the nodes that are common between all the time points of cold treatment (centrally located in Figure 6) exhibit an increasing connectivity trend that reaches a plateau at 12 h of cold acclimation (Figure 7A). For example, sucrose biosynthesis II showed the most significant pattern of increase, suggesting that these networks are highly dynamic in their ability to rewire over time of cold acclimation. When we compared the four time points by overlaying networks on each other, patterns of connectivity over the course of stress were visible, including the existence of mutually exclusive nodes and edges that decompose graphs into maximal non-separable subgraph blocks that are nearly non-overlappable. The number of exclusive nodes increased with time but reached a plateau at 12 h of acclimation. Similarly, the number of exclusive edges increased with time throughout the timeline (Table 3). These trends fit patterns of network clustering that suggest the network is modular (see below).

FIGURE 6

Figure 6. The union of tripartite networks. The pie-node representation enables to identify presence/absence of individual nodes across four time dependent networks.

FIGURE 7

Figure 7. Analysis of network structure. (A) Connectivity of common nodes for all time points. (B) Analysis of network modularity using the clustering algorithm and some statistical descriptor of power law behavior of tripartite networks over time. (C) Boxplots show the span of connectivity measured by node degree for each time point. (D) Analysis of modularity with the FGC (closed circles) and NG (open circles) algorithms.

The connectivity of a network can be characterized by the probability P(k) that a node has k connections. For inhomogeneous networks, highly connected nodes get more connections. In these random networks, the probability P(k) of a node being connected to k other nodes decays as a power law, P(k) ~ k^−γ, without a characteristic scale (Jeong et al., 2000). For example, metabolic networks have a tendency to follow scale free network behavior with γ = 2.1–2.4 exponents (Jeong et al., 2000). In order to test if the tripartite networks at different times of cold acclimation tend to follow the scale-free distribution we measured power law behavior with network statistics (Table 3). We plotted log P(k) vs. log k for all tripartite networks (Figure S4). Log-log plots showed a limited correlation fit [Kolmogorov–Smirnov (KS) fit ranging 0.4–0.46] to power law behavior (γ ranging from 1.06 to 1.37). The γ for the tripartite networks is lower than the γ reported for metabolic networks (Jeong et al., 2000). Determination coefficients (R²) of ~ 70% support the linear model. KS fit statistic higher than >0.10 and low P-values of the KS test (< 0.05) (Clauset et al., 2009) rejected network data being drawn from the fitted power-law distribution for all cold acclimation time points. Similarly, the log-likelihoods of the fitted power law parameters were much lower than zero for all networks. This suggests the power law distribution is unlikely. Overall, the statistical analysis of tripartite networks over time rejected power law behavior (Figure 7B, Table 3). This reveals that tripartite networks were not dictated by hubs. Instead they harbored large components of highly rewired vertices.

Boxplot graphs show connectivity variation measured by the number of incident connections of a node (its degree) (Figure 7C). Variation was minimal at 3 h cold acclimation and maximal at 24 h. Results explicitly show that stress propels the dynamic behavior of tripartite networks. Networks can show modular behavior. Modularity can be estimated with the average clustering coefficient (C), which describes the ratio of connected triplets in the graph normalized over all vertices and ignoring the direction and weights of the edges (Aziz et al., 2016). From a topological perspective, the tripartite networks of Arabidopsis have high clustering coefficients, with values that increased from 0.69 to 0.77 along the cold acclimation sequence (Table 3). This property is suggestive of modular organization. This trend corresponds to the KS fit statistics which in turn reject power law behavior. The modularity of tripartite networks was assessed using the FGC and NG algorithms that reveals scores of communities in the networks. In particular, FGC provides a hierarchical perspective of agglomerative community structure. Positive index values indicate the significant existence of modules in the connectivity structure of the network. Notably, FGC and NG similarly showed high modularity upon stress application (Figure 7D). Indices decreased at 12 and 24 h post-acclimation but remained positive. Analysis of modularity shows that connectivity increases (exclusive nodes and edges, densely interconnected functional modules) along with high cluster coefficients (Figures 6, 7, Table 3). High indices indicate that the connectivity of vertices belonging to modules is higher than those between modules. In turn, low indices indicate high community structure. High modularity appears in the networks as Arabidopsis progressed toward new programming of metabolic diversity and increased modularity upon perturbation. This suggests that tripartite networks under cold stress are characterized by a high intrinsic potential modularity. Apparently, this pattern was dictated by a development of network clusters within modules outside the context of scale free organization. Interestingly, shared vertices (commonly present in four networks) were highly connected and located in the center of the network while individual vertices were located toward the periphery (Figure 6). Among the common pathways for all time points, carbohydrate, lipid, cellulose, signal and redox metabolisms were explicit. This observation supports stress-induced trends toward redox metabolism, signal metabolism, carbon metabolism and cell wall metabolism (Supplementary File 2). This feature also uncovers the rewiring change of network with stable components of physiology. On the other hand, we see dynamic change of metabolites and pathways as a response to stress application, which evolves modular structure of the tripartite networks.

Network Cluster Analysis

In addition to modules, extracting clusters can identify highly connected subgraphs in which nodes have more interactions than the rest. The association between components of a cluster can reduce network complexity. We used cluster analysis to highlight key metabolic pathways and delineate central associations between pathways. Using cluster analysis, we identified >20 clusters of sizes from 3 to 27 (data not shown). For example, extraction of a sub-network from the 3 h cold acclimation network showed redox metabolisms linked to each other (Figure 8A). Interestingly, the jasmonic acid pathway was highly connected to the TCA cycle, acetyl–CoA biosynthesis, pyruvate decarboxylation to acetyl-CoA, and ethanol degradation pathways (Figure 8B). Jasmonate targets the mitochondria of cancer cells by perturbing them (Rotem et al., 2005). Similarly, jasmonate may induce perturbation of mitochondrial function in plants under stress conditions. It would be interesting to validate this hypothesis in future studies. Furthermore, plants seem to increase the glycolysis and ethanol degradation pathways under cold stress since the bio-energetic machinery is impaired. The resulting sub-network again supported the notion that the TCA cycle works alongside ethanol degradation to produce excess of acetyl-CoA. Focusing on the Rubisco shunt and Calvin cycle interaction, we extracted a sub-network from the 12 h cold acclimation network. The resulting sub-network also showed main carbon metabolic pathways were highly connected to each other (Figure 8C). Plants confront many imbalances in carbon metabolism under perturbation. As expected, photosynthesis is significantly down-regulated whereas starch, sucrose and amino acid metabolisms are up-regulated in such stress conditions. This suggested that when photosynthesis is perturbed, plants may increase the activity of alternative pathways to modulate CO₂ fixation, including activating the Rubisco shunt. This involves maintaining a steady-state of energy metabolism (Figure 8C).

FIGURE 8

Figure 8. Sub-network association of pathways at selected time points. (A) Sub-network of redox metabolism at 3 h cold acclimation. (B) Sub-network of energy metabolism at 3 h cold acclimation. (C) Sub-network of carbon metabolism at 12 h cold acclimation. Red color range shows most significant genes, metabolites, and pathways.

Conclusion

A metabolite-centric reporter pathway analysis of transcriptomes uncovered a number of metabolites and pathways that differentially modulate the plant response to cold stress in Arabidopsis. Notably, cold stress levies great demands on energy production and biochemical requirements. Initially, plants mobilize glycolysis to rapidly produce ATP and divert the main carbon flux routes to amino acid and lipid metabolic pathways. In this process, many metabolic pathways are activated to produce cellular materials needed to sustain stress. We find that energy rich metabolites, such as ethanol, may help the functioning of the TCA cycle under perturbation. Signal pathways and structural integrity related pathways were also activated. Clearly, our metabolite-centric reporter pathway analysis provides a new framework with which to help scientists achieve a deeper understanding of the co-regulation of metabolites and metabolic pathways. Our analysis also revealed modularity and rejection of the power law behavior in the structure of tripartite networks of gene-metabolite-pathway connectivity.

Author Contributions

IK and IY performed experiments. IK, IY, and GC-A analyzed the data. IK and GC-A wrote the manuscript.

Funding

Research at Illinois is supported by USDA National Institute of Food and Agriculture, Hatch project 1014249 and a Blue Waters supercomputer allocation to GC-A.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe.2018.00121/full#supplementary-material

Figure S1. Principal component analysis (PCA) of control and cold acclimated samples derived from all genes from the Affymetrix ATH1 microarray chip (A) and the subset of metabolic genes (B)

Figure S2. Subgraph of tripartite network showing significantly changed metabolites (P ≤ 0.05) for 3 h (A), 6 h (B), 12 h (C), and 24 h (D) cold acclimated Arabidopsis. Nodes represent metabolites. Their sizes were scaled to number of neighbors.

Figure S3. Mapping gene expression onto the ethylene pathway of the KEGG database (red circle). The graph was produced with the R pathview library. Rectangles represent enzyme-encoding genes labeled with Enzyme Commission (EC) number. Red color hues indicate gene up-regulation while green color hues indicate down-regulation. EC:1.14.17.4, aminocyclopropanecarboxylate oxidase.

Figure S4. Log-log plots show a scatter plot of natural log of k against P(k) of the tripartite networks for each time point.

Abbreviations

RPA^m, Metabolic centric reporter pathway analysis.

References

Akkiprik, M., Peker, I., Özmen, T., Amuran, G., Güllüoglu, B., Kaya, H., et al. (2015). Identification of differentially expressed IGFBP5-related genes in breast cancer tumor tissues using cDNA microarray experiments. Genes 6, 1201–1214. doi: 10.3390/genes6041201

PubMed Abstract | CrossRef Full Text | Google Scholar

Aubert, S. (2004). Methyl- -D-glucopyranoside in higher plants: accumulation and intracellular localization in Geum montanum L. leaves and in model systems studied by 13C nuclear magnetic resonance. J. Exp. Bot. 55, 2179–2189. doi: 10.1093/jxb/erh235

PubMed Abstract | CrossRef Full Text | Google Scholar

Aziz, M. F., Caetano-Anollés, K., and Caetano-Anollés, G. (2016). The early history and emergence of molecular functions and modular scale-free network behavior. Sci. Rep. 6:25058. doi: 10.1038/srep25058

PubMed Abstract | CrossRef Full Text | Google Scholar

Badejo, A. A., Tanaka, N., and Esaka, M. (2008). Analysis of GDP-D-mannose pyrophosphorylase gene promoter from Acerola (Malpighia glabra) and Increase in ascorbate content of transgenic tobacco expressing the Acerola gene. Plant Cell Physiol. 49, 126–132. doi: 10.1093/pcp/pcm164

PubMed Abstract | CrossRef Full Text | Google Scholar

Bakht, J., Bano, A., and Dominy, P. (2006). The role of abscisic acid and low temperature in chickpea (Cicer arietinum) cold tolerance. II. Effects on plasma membrane structure and function. J. Exp. Bot. 57, 3707–3715. doi: 10.1093/jxb/erl120

PubMed Abstract | CrossRef Full Text | Google Scholar

Bielen, A., Remans, T., Vangronsveld, J., and Cuypers, A. (2013). The influence of metal stress on the availability and redox state of ascorbate, and possible interference with its cellular functions. Int. J. Mol. Sci. 14, 6382–6413. doi: 10.3390/ijms14036382

PubMed Abstract | CrossRef Full Text | Google Scholar

Biswal, B., Joshi, P. N., Raval, M. K., and Biswal, U. C. (2011). Photosynthesis, a global sensor of environmental stress in green plants: stress signalling and adaptation. Curr. Sci. 101, 47–56. Available online at: http://www.jstor.org/stable/24077862

Google Scholar

Budak, H., Akpinar, B. A., Unver, T., and Turktas, M. (2013). Proteome changes in wild and modern wheat leaves upon drought stress by two-dimensional electrophoresis and nanoLC-ESI–MS/MS. Plant Mol. Biol. 83, 89–103. doi: 10.1007/s11103-013-0024-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Bulley, S., Wright, M., Rommens, C., Yan, H., Rassam, M., Lin-Wang, K., et al. (2012). Enhancing ascorbate in fruits and tubers through over-expression of the l-galactose pathway gene GDP-l-galactose phosphorylase. Plant Biotechnol. J. 10, 390–397. doi: 10.1111/j.1467-7652.2011.00668.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Çakir, T. (2015). Reporter pathway analysis from transcriptome data: metabolite-centric versus reaction-centric approach. Sci. Rep. 5:14563. doi: 10.1038/srep14563

PubMed Abstract | CrossRef Full Text | Google Scholar

Cakir, T., Patil, K. R., Önsan, Z. I., Ülgen, K. Ö., Kirdar, B., and Nielsen, J. (2006). Integration of metabolome data with metabolic networks reveals reporter reactions. Mol. Syst. Biol. 2:50. doi: 10.1038/msb4100085

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y.-F., Etheridge, N., and Schaller, G. E. (2005). Ethylene signal transduction. Ann. Bot. 95, 901–915. doi: 10.1093/aob/mci100

PubMed Abstract | CrossRef Full Text | Google Scholar

Clauset, A., Newman, M. E., and Moore, C. (2004). Finding community structure in very large networks. Phys. Rev. E70:066111. doi: 10.1103/PhysRevE.70.066111

CrossRef Full Text | Google Scholar

Clauset, A., Shalizi, C. R., and Newman, M. E. (2009). Power-law distributions in empirical data. SIAM Rev. 51, 661–703. doi: 10.1137/070710111

CrossRef Full Text | Google Scholar

Cook, D., Fowler, S., Fiehn, O., and Thomashow, M. F. (2004). A prominent role for the CBF cold response pathway in configuring the low-temperature metabolome of Arabidopsis. Proc. Natl. Acad. Sci. U.S.A. 101, 15243–15248. doi: 10.1073/pnas.0406069101

PubMed Abstract | CrossRef Full Text | Google Scholar

Csardi, G., and Nepusz, T. (2006). The igraph software package for complex network research. Inter J. Complex Syst. 1695, 1–9. Available online at: http://igraph.sf.net

Google Scholar

Dal'Molin, C. G., Quek, L. E., Palfreyman, R. W., and Nielsen, L. K. (2011). AlgaGEM–a genome-scale metabolic reconstruction of algae based on the Chlamydomonas reinhardtii genomme. BMC Genomics 12(Suppl. 4):S5. doi: 10.1186/1471-2164-12-S4-S5

PubMed Abstract | CrossRef Full Text | Google Scholar

Davies, R. T., Goetz, D. H., Lasswell, J., Anderson, M. N., and Bartel, B. (1999). IAR3 encodes an auxin conjugate hydrolase from Arabidopsis. Plant Cell 11, 365–376. doi: 10.1105/tpc.11.3.365

PubMed Abstract | CrossRef Full Text | Google Scholar

de Oliveira Dal'Molin, C. G., Quek, L.-E., Palfreyman, R. W., Brumbley, S. M., and Nielsen, L. K. (2010). AraGEM, a genome-scale reconstruction of the primary metabolic network in Arabidopsis. Plant Physiol. 152, 579–589. doi: 10.1104/pp.109.148817

PubMed Abstract | CrossRef Full Text | Google Scholar

DeWald, D. B., Torabinejad, J., Jones, C. A., Shope, J. C., Cangelosi, A. R., Thompson, J. E., et al. (2001). Rapid accumulation of phosphatidylinositol 4, 5-bisphosphate and inositol 1, 4, 5-trisphosphate correlates with calcium mobilization in salt-stressed Arabidopsis. Plant Physiol. 126, 759–769. doi: 10.1104/pp.126.2.759

PubMed Abstract | CrossRef Full Text | Google Scholar

Fall, R., and Benson, A. A. (1996). Leaf methanol—the simplest natural product from plants. Trends Plant Sci. 1, 296–301. doi: 10.1016/S1360-1385(96)88175-0

CrossRef Full Text | Google Scholar

Gass, N., Glagotskaia, T., Mellema, S., Stuurman, J., Barone, M., Mandel, T., et al. (2005). Pyruvate decarboxylase provides growing pollen tubes with a competitive advantage in petunia. Plant Cell 17, 2355–2368. doi: 10.1105/tpc.105.033290

PubMed Abstract | CrossRef Full Text | Google Scholar

Gatzek, S., Wheeler, G. L., and Smirnoff, N. (2002). Antisense suppression of l-galactose dehydrogenase in Arabidopsis thaliana provides evidence for its role in ascorbate synthesis and reveals light modulated l-galactose synthesis. Plant J. 30, 541–553. doi: 10.1046/j.1365-313X.2002.01315.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Gout, E., Aubert, S., Bligny, R., Rébeillé, F., Nonomura, A. R., Benson, A. A., et al. (2000). Metabolism of methanol in plant cells. Carbon-13 nuclear magnetic resonance studies. Plant Physiol. 123, 287–296. doi: 10.1104/pp.123.1.287

PubMed Abstract | CrossRef Full Text | Google Scholar

Hannah, M. A., Heyer, A. G., and Hincha, D. K. (2005). A global survey of gene regulation during cold acclimation in Arabidopsis thaliana. PLoS Genet. 1:e26. doi: 10.1371/journal.pgen.0010026

PubMed Abstract | CrossRef Full Text | Google Scholar

Herrmann, K. M., and Weaver, L. M. (1999). The shikimate pathway. Annu. Rev. Plant Biol. 50, 473–503. doi: 10.1146/annurev.arplant.50.1.473

PubMed Abstract | CrossRef Full Text | Google Scholar

Hurry, V., Strand, A., Furbank, R., and Stitt, M. (2000). The role of inorganic phosphate in the development of freezing tolerance and the acclimatization of photosynthesis to low temperature is revealed by the pho mutants of Arabidopsis thaliana. Plant J. 24, 383–396. doi: 10.1046/j.1365-313x.2000.00888.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Ideker, T., Ozier, O., Schwikowski, B., and Siegel, A. F. (2002). Discovering regulatory and signalling circuits in molecular interaction networks. Bioinformatics 18, S233–S240. doi: 10.1093/bioinformatics/18.suppl_1.S233

PubMed Abstract | CrossRef Full Text | Google Scholar

Ishikawa, T., Nishikawa, H., Gao, Y., Sawa, Y., Shibata, H., Yabuta, Y., et al. (2008). The pathway via D-galacturonate/L-galactonate is significant for ascorbate biosynthesis in Euglena gracilis Identification and functional characterization of aldonolactonase. J. Biol. Chem. 283, 31133–31141. doi: 10.1074/jbc.M803930200

PubMed Abstract | CrossRef Full Text | Google Scholar

Jain, M., and Khurana, J. P. (2009). Transcript profiling reveals diverse roles of auxin-responsive genes during reproductive development and abiotic stress in rice. FEBS J. 276, 3148–3162. doi: 10.1111/j.1742-4658.2009.07033.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Jeon, J., and Kim, J. (2012). Arabidopsis response regulator 1 (ARR1) and Arabidopsis histidine phosphotransfer protein 2 (AHP2), AHP3, and AHP5 function in cold signaling. Plant Physiol. 112:207621. doi: 10.1104/pp.112.207621

CrossRef Full Text | Google Scholar

Jeon, J., Kim, N. Y., Kim, S., Kang, N. Y., Novák, O., Ku, S.-J., et al. (2010). A subset of cytokinin two-component signaling system plays a role in cold temperature stress response in Arabidopsis. J. Biol. Chem. 285, 23371–23386. doi: 10.1074/jbc.M109.096644

PubMed Abstract | CrossRef Full Text | Google Scholar

Jeong, H., Tombor, B., Albert, R., Oltvai, Z. N., and Barabási, A.-L. (2000). The large-scale organization of metabolic networks. Nature 407, 651–654. doi: 10.1038/35036627

PubMed Abstract | CrossRef Full Text | Google Scholar

Jung, J.-Y., Kim, Y.-W., Kwak, J. M., Hwang, J.-U., Young, J., Schroeder, J. I., et al. (2002). Phosphatidylinositol 3-and 4-phosphate are required for normal stomatal movements. Plant Cell 14, 2399–2412. doi: 10.1105/tpc.004143

PubMed Abstract | CrossRef Full Text | Google Scholar

Kanehisa, M., and Goto, S. (2000). KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 28, 27–30. doi: 10.1093/nar/28.1.27

PubMed Abstract | CrossRef Full Text | Google Scholar

Kaplan, F., Kopka, J., Haskell, D. W., Zhao, W., Schiller, K. C., Gatzke, N., et al. (2004). Exploring the temperature-stress metabolome of Arabidopsis. Plant Physiol. 136, 4159–4168. doi: 10.1104/pp.104.052142

PubMed Abstract | CrossRef Full Text | Google Scholar

Kilian, J., Whitehead, D., Horak, J., Wanke, D., Weinl, S., Batistic, O., et al. (2007). The AtGenExpress global stress expression data set: protocols, evaluation and model data analysis of UV-B light, drought and cold stress responses. Plant J. 50, 347–363. doi: 10.1111/j.1365-313X.2007.03052.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, T., Dreher, K., Nilo-Poyanco, R., Lee, I., Fiehn, O., Lange, B. M., et al. (2015). Patterns of metabolite changes identified from large-scale gene perturbations in Arabidopsis using a genome-scale metabolic network. Plant Physiol. 167, 1685–1698. doi: 10.1104/pp.114.252361

PubMed Abstract | CrossRef Full Text | Google Scholar

Koç, I., Çalisşkaner, Z. O., and Filiz, E. (2015c). Application of data analysis in cold stress: a case study of Nicotiana benthamiana. Turk. J. Bot. 39, 1021–1032. doi: 10.3906/bot-1503-25

CrossRef Full Text | Google Scholar

Koç, I., Filiz, E., and Tombuloglu, H. (2015a). Assessment of miRNA expression profile and differential expression pattern of target genes in cold-tolerant and cold-sensitive tomato cultivars. Biotechnol. Biotechnol. Equip. 29, 851–860. doi: 10.1080/13102818.2015.1061447

CrossRef Full Text | Google Scholar

Koç, I., Vatansever, R., Ozyigit, I. I., and Filiz, E. (2015b). Identification of differentially expressed genes in chilling-induced potato (Solanum tuberosum L.); a data analysis study. Appl. Biochem. Biotechnol. 177, 792–811. doi: 10.1007/s12010-015-1778-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Kost, B., Lemichez, E., Spielhofer, P., Hong, Y., Tolias, K., Carpenter, C., et al. (1999). Rac homologues and compartmentalized phosphatidylinositol 4, 5-bisphosphate act in a common pathway to regulate polar pollen tube growth. J. Cell Biol. 145, 317–330. doi: 10.1083/jcb.145.2.317

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuntal, B. K., Dutta, A., and Mande, S. S. (2016). CompNet: a GUI based tool for comparison of multiple biological interaction networks. BMC Bioinformatics 17:185. doi: 10.1186/s12859-016-1013-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Lakshmanan, M., Lim, S.-H., Mohanty, B., Kim, J. K., Ha, S.-H., and Lee, D.-L. (2015). Unraveling the light-specific metabolic and regulatory signatures of rice through combined in silico modeling and multi-omics analysis. Plant Physiol. 169, 3002–3020. doi: 10.1104/pp.15.01379

CrossRef Full Text | Google Scholar

Le Gall, H., Philippe, F., Domon, J.-M., Gillet, F., Pelloux, J., and Rayon, C. (2015). Cell wall metabolism in response to abiotic stress. Plants 4, 112–166. doi: 10.3390/plants4010112

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, H., Ma, Y., Chen, N. A., Guo, S., Liu, H., Guo, X., et al. (2014). Overexpression of stress-inducible OsBURP16, the β subunit of polygalacturonase 1, decreases pectin content and cell adhesion and increases abiotic stress sensitivity in rice. Plant Cell Environ. 37, 1144–1158. doi: 10.1111/pce.12223

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, K., Li, L., and Luan, S. (2005). An essential function of phosphatidylinositol phosphates in activation of plant shaker-type K+ channels. Plant J. 42, 433–443. doi: 10.1111/j.1365-313X.2005.02384.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, L., Shen, F., Xin, C., and Wang, Z. (2016). Multi-scale modeling of Arabidopsis thaliana response to different CO2 conditions: from gene expression to metabolic flux. J. Integr. Plant Biol. 58, 2–11. doi: 10.1111/jipb.12370

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, X.-G., Xu, H., Zhang, J.-Y., Liang, G.-W., Liu, Y.-T., and Guo, A.-G. (2012). Effect of low temperature on chlorophyll biosynthesis in albinism line of wheat (Triticum aestivum) FA85. Physiol. Plant. 145, 384–394. doi: 10.1111/j.1399-3054.2012.01604.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Y., Yu, L., and Wang, R. (2011). Level of ascorbic acid in transgenic rice for L-galactono-1, 4-lactone dehydrogenase overexpressing or suppressed is associated with plant growth and seed set. Acta Physiol. Plant. 33, 1353–1363. doi: 10.1007/s11738-010-0669-5

CrossRef Full Text | Google Scholar

Ljung, K. (2013). Auxin metabolism and homeostasis during plant development. Development 140, 943–950. doi: 10.1242/dev.086363

PubMed Abstract | CrossRef Full Text | Google Scholar

Los, D. A., and Murata, N. (2004). Membrane fluidity and its roles in the perception of environmental signals. Biochim. Biophys. Acta BBA Biomembr. 1666, 142–157. doi: 10.1016/j.bbamem.2004.08.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Ludwig-Müller, J., Kaldorf, M., Sutter, E. G., and Epstein, E. (1997). Indole-3-butyric acid (IBA) is enhanced in young maize (Zea mays L.) roots colonized with the arbuscular mycorrhizal fungus Glomus intraradices. Plant Sci. 125, 153–162. doi: 10.1016/S0168-9452(97)00064-2

CrossRef Full Text | Google Scholar

Ludwig-Müller, J., Schubert, B., and Pieper, K. (1995). Regulation of IBA synthetase from maize (Zea mays L.) by drought stress and ABA. J. Exp. Bot. 46, 423–432. doi: 10.1093/jxb/46.4.423

CrossRef Full Text | Google Scholar

Matteucci, M., D'angeli, S., Errico, S., Lamanna, R., Perrotta, G., and Altamura, M. M. (2011). Cold affects the transcription of fatty acid desaturases and oil quality in the fruit of Olea europaea L. genotypes with different cold hardiness. J. Exp. Bot. 62, 3403–3420. doi: 10.1093/jxb/err013

PubMed Abstract | CrossRef Full Text | Google Scholar

Mellema, S., Eichenberger, W., Rawyler, A., Suter, M., Tadege, M., and Kuhlemeier, C. (2002). The ethanolic fermentation pathway supports respiration and lipid biosynthesis in tobacco pollen. Plant J. 30, 329–336. doi: 10.1046/j.1365-313X.2002.01293.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Meyerowitz, E. M., and Somerville, C. R. (1994). Arabidopsis. Cold Spring Harbor NY: Cold Spring Harbor Laboratory Press.

Google Scholar

Milla, M. A. R., Maurer, A., Huete, A. R., and Gustafson, J. P. (2003). Glutathione peroxidase genes in Arabidopsis are ubiquitous and regulated by abiotic stresses through diverse signaling pathways. Plant J. 36, 602–615. doi: 10.1046/j.1365-313X.2003.01901.x

CrossRef Full Text | Google Scholar

Mittler, R. (2002). Oxidative stress, antioxidants and stress tolerance. Trends Plant Sci. 7, 405–410. doi: 10.1016/S1360-1385(02)02312-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Mohnen, D., Doong, R. L., Liljebjelke, K., Fralish, G., and Chan, J. (1996). Cell free synthesis of the pectic polysaccharide homogalacturonan. Prog. Biotechnol. 14, 109–126. doi: 10.1016/S0921-0423(96)80250-4

CrossRef Full Text | Google Scholar

Møller, I. M. (2001). Plant mitochondria and oxidative stress: electron transport, NADPH turnover, and metabolism of reactive oxygen species. Annu. Rev. Plant Biol. 52, 561–591. doi: 10.1146/annurev.arplant.52.1.561

PubMed Abstract | CrossRef Full Text | Google Scholar

Mootha, V. K., Lindgren, C. M., Eriksson, K.-F., Subramanian, A., Sihag, S., Lehar, J., et al. (2003). PGC-1α-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes. Nat. Genet. 34, 267–273. doi: 10.1038/ng1180

PubMed Abstract | CrossRef Full Text | Google Scholar

Mueller, L. A., Zhang, P., and Rhee, S. Y. (2003). AraCyc: a biochemical pathway database for Arabidopsis. Plant Physiol. 132, 453–460. doi: 10.1104/pp.102.017236

PubMed Abstract | CrossRef Full Text | Google Scholar

Neuhaus, H. E., and Emes, M. J. (2000). Nonphotosynthetic metabolism in plastids. Annu. Rev. Plant Biol. 51, 111–140. doi: 10.1146/annurev.arplant.51.1.111

PubMed Abstract | CrossRef Full Text | Google Scholar

Newman, M. E. (2005). Power laws, Pareto distributions and Zipf's law. Contemp. Phys. 46, 323–351. doi: 10.1016/j.cities.2012.03.001

CrossRef Full Text | Google Scholar

Newman, M. E., and Girvan, M. (2004). Finding and evaluating community structure in networks. Phys. Rev. E69:026113. doi: 10.1103/PhysRevE.69.026113

CrossRef Full Text | Google Scholar

Nielsen, J. (2003). It is all about metabolic fluxes. J. Bacteriol. 185, 7031–7035. doi: 10.1128/JB.1805.24.7031-7035.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

Obermeyer, G., Fragner, L., Lang, V., and Weckwerth, W. (2013). Dynamic adaption of metabolic pathways during germination and growth of lily pollen tubes after inhibition of the electron transport chain. Plant Physiol. 162, 1822–1833. doi: 10.1104/pp.113.219857

PubMed Abstract | CrossRef Full Text | Google Scholar

Oliveira, A. P., Patil, K. R., and Nielsen, J. (2008). Architecture of transcriptional regulatory circuits is knitted over the topology of bio-molecular interaction networks. BMC Syst. Biol. 2:17. doi: 10.1186/1752-0509-2-17

PubMed Abstract | CrossRef Full Text | Google Scholar

op den Camp, R. G., Przybyla, D., Ochsenbein, C., Laloi, C., Kim, C., Danon, A., et al. (2003). Rapid induction of distinct stress responses after the release of singlet oxygen in Arabidopsis. Plant Cell 15, 2320–2332. doi: 10.1105/tpc.014662

PubMed Abstract | CrossRef Full Text | Google Scholar

Partelli, F. L., Vieira, H. D., Viana, A. P., Batista-Santos, P., Rodrigues, A. P., Leitão, A. E., et al. (2009). Low temperature impact on photosynthetic parameters of coffee genotypes. Pesqui. Agropecuária Bras. 44, 1404–1415. doi: 10.1590/S0100-204X2009001100006

CrossRef Full Text | Google Scholar

Patil, K. R., and Nielsen, J. (2005). Uncovering transcriptional regulation of metabolism by using metabolic network topology. Proc. Natl. Acad. Sci. U.S.A. 102, 2685–2689. doi: 10.1073/pnas.0406811102

PubMed Abstract | CrossRef Full Text | Google Scholar

Peleg, Z., and Blumwald, E. (2011). Hormone balance and abiotic stress tolerance in crop plants. Curr. Opin. Plant Biol. 14, 290–295. doi: 10.1016/j.pbi.2011.02.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Rajashekar, C. B., and Lafta, A. (1996). Cell-wall changes and cell tension in response to cold acclimation and exogenous abscisic acid in leaves and cell cultures. Plant Physiol. 111, 605–612. doi: 10.1104/pp.111.2.605

PubMed Abstract | CrossRef Full Text | Google Scholar

Rekarte-Cowie, I., Ebshish, O. S., Mohamed, K. S., and Pearce, R. S. (2008). Sucrose helps regulate cold acclimation of Arabidopsis thaliana. J. Exp. Bot. 59, 4205–4217. doi: 10.1093/jxb/ern262

PubMed Abstract | CrossRef Full Text | Google Scholar

Rensink, W. A., Iobst, S., Hart, A., Stegalkina, S., Liu, J., and Buell, C. R. (2005). Gene expression profiling of potato responses to cold, heat, and salt stress. Funct. Integr. Genomics 5, 201–207. doi: 10.1007/s10142-005-0141-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Rodríguez, V. M., Velasco, P., Garrido, J. L., Revilla, P., Ordás, A., and Butrón, A. (2013). Genetic regulation of cold-induced albinism in the maize inbred line A661. J. Exp. Bot. 64, 3657–3667. doi: 10.1093/jxb/ert189

PubMed Abstract | CrossRef Full Text | Google Scholar

Rotem, R., Heyfets, A., Fingrut, O., Blickstein, D., Shaklai, M., and Flescher, E. (2005). Jasmonates: novel anticancer agents acting directly and selectively on human cancer cell mitochondria. Cancer Res. 65, 1984–1993. doi: 10.1158/0008-5472.CAN-04-3091

PubMed Abstract | CrossRef Full Text | Google Scholar

Saha, R., Suthers, P. F., and Maranas, C. D. (2011). Zea mays iRS1563: a comprehensive genome-scale metabolic reconstruction of maize metabolism. PLoS ONE 6:e21784. doi: 10.1371/journal.pone.0021784

PubMed Abstract | CrossRef Full Text | Google Scholar

Santner, A., and Estelle, M. (2009). Recent advances and emerging trends in plant hormone signalling. Nature 459, 1071–1078. doi: 10.1038/nature08122

PubMed Abstract | CrossRef Full Text | Google Scholar

Savitch, L. V., Barker-Astrom, J., Ivanov, A. G., Hurry, V., Öquist, G., Huner, N. P., et al. (2001). Cold acclimation of Arabidopsis thaliana results in incomplete recovery of photosynthetic capacity, associated with an increased reduction of the chloroplast stroma. Planta 214, 295–303. doi: 10.1007/s004250100622

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwender, J., Goffman, F., Ohlrogge, J. B., and Shachar-Hill, Y. (2004). Rubisco without the Calvin cycle improves the carbon efficiency of developing green seeds. Nature 432, 779–782. doi: 10.1038/nature03145

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, Y., Tian, S., Hou, L., Huang, X., Zhang, X., Guo, H., et al. (2012). Ethylene signaling negatively regulates freezing tolerance by repressing expression of CBF and type-A ARR genes in Arabidopsis. Plant Cell 24, 2578–2595. doi: 10.1105/tpc.112.098640

PubMed Abstract | CrossRef Full Text | Google Scholar

Shibasaki, K., Uemura, M., Tsurumi, S., and Rahman, A. (2009). Auxin response in Arabidopsis under cold stress: underlying molecular mechanisms. Plant Cell 21, 3823–3838. doi: 10.1105/tpc.109.069906

PubMed Abstract | CrossRef Full Text | Google Scholar

Singh, I., Kumar, U., Singh, S. K., Gupta, C., Singh, M., and Kushwaha, S. R. (2012). Physiological and biochemical effect of 24-epibrassinoslide on cold tolerance in maize seedlings. Physiol. Mol. Biol. Plants 18, 229–236. doi: 10.1007/s12298-012-0122-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Solecka, D., Zebrowski, J., and Kacperska, A. (2008). Are pectins involved in cold acclimation and de-acclimation of winter oil-seed rape plants? Ann. Bot. 101, 521–530. doi: 10.1093/aob/mcm329

PubMed Abstract | CrossRef Full Text | Google Scholar

Staal, F. J. T., van der Burg, M., Wessels, L. F. A., Barendregt, B. H., Baert, M. R. M., van den Burg, C. M. M., et al. (2003). DNA microarrays for comparison of gene expression profiles between diagnosis and relapse in precursor-B acute lymphoblastic leukemia: choice of technique and purification influence the identification of potential diagnostic markers. Leukemia 17, 1324–1732. doi: 10.1038/sj.leu.2402974

PubMed Abstract | CrossRef Full Text | Google Scholar

Staswick, P. E., Serban, B., Rowe, M., Tiryaki, I., Maldonado, M. T., Maldonado, M. C., et al. (2005). Characterization of an Arabidopsis enzyme family that conjugates amino acids to indole-3-acetic acid. Plant Cell 17, 616–627. doi: 10.1105/tpc.104.026690

PubMed Abstract | CrossRef Full Text | Google Scholar

Strand, A., Hurry, V., Gustafsson, P., and Gardeström, P. (1997). Development of Arabidopsis thaliana leaves at low temperatures releases the suppression of photosynthesis and photosynthetic gene expression despite the accumulation of soluble carbohydrates. Plant J. 12, 605–614. doi: 10.1046/j.1365-313X.1997.00605.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Tatli, Ö., Nikerel, I. E., and Özdemir, B. S. (2015). Evaluation of metabolite extraction protocols and determination ofphysiological response to drought stress via reporter metabolites in model plant Brachypodium distachyon. Turk. J. Bot. 39, 1042–1050. doi: 10.3906/bot-1505-28

CrossRef Full Text | Google Scholar

Tegge, A. N., Caldwell, C. W., and Xu, D. (2012). Pathway correlation profile of gene-gene co-expression for identifying pathway perturbation. PLoS ONE 7:e52127. doi: 10.1371/journal.pone.0052127

PubMed Abstract | CrossRef Full Text | Google Scholar

Tewari, A. K., and Tripathy, B. C. (1998). Temperature-stress-induced impairment of chlorophyll biosynthetic reactions in cucumber and wheat. Plant Physiol. 117, 851–858. doi: 10.1104/pp.117.3.851

CrossRef Full Text | Google Scholar

Thomashow, M. F. (2001). So what's new in the field of plant cold acclimation? Lots! Plant Physiol. 125, 89–93. doi: 10.1104/pp.125.1.89

PubMed Abstract | CrossRef Full Text | Google Scholar

Urano, K., Maruyama, K., Ogata, Y., Morishita, Y., Takeda, M., Sakurai, N., et al. (2009). Characterization of the ABA-regulated global responses to dehydration in Arabidopsis by metabolomics. Plant J. 57, 1065–1078. doi: 10.1111/j.1365-313X.2008.03748.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Väremo, L., Nielsen, J., and Nookaew, I. (2013). Enriching the gene set analysis of genome-wide data by incorporating directionality of gene expression and combining statistical hypotheses and methods. Nucleic Acids Res. 41, 4378–4391. doi: 10.1093/nar/gkt111

PubMed Abstract | CrossRef Full Text | Google Scholar

Vriezen, W. H., Hulzink, R., Mariani, C., and Voesenek, L. A. (1999). 1-Aminocyclopropane-1-carboxylate oxidase activity limits ethylene biosynthesis in Rumex palustris during submergence. Plant Physiol. 121, 189–196. doi: 10.1104/pp.121.1.189

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Q.-J., Sun, H., Dong, Q.-L., Sun, T.-Y., Jin, Z.-X., Hao, Y.-J., et al. (2016). The enhancement of tolerance to salt and cold stresses by modifying the redox state and salicylic acid content via the cytosolic malate dehydrogenase gene in transgenic apple plants. Plant Biotechnol. J. 14, 1986–1997. doi: 10.1111/pbi.12556

PubMed Abstract | CrossRef Full Text | Google Scholar

Welti, R., Li, W., Li, M., Sang, Y., Biesiada, H., Zhou, H.-E., et al. (2002). Profiling membrane lipids in plant stress responses role of phospholipase Dα in freezing-induced lipid changes in Arabidopsis. J. Biol. Chem. 277, 31994–32002. doi: 10.1074/jbc.M205375200

CrossRef Full Text | Google Scholar

Yang, J. J., Bhojwani, D., Yang, W., Cai, X., Stocco, G., Crews, K., et al. (2008). Genome-wide copy number profiling reveals molecular evolution from diagnosis to relapse in childhood acute lymphoblastic leukemia. Blood 112, 4178–4183. doi: 10.1182/blood-2008-06-165027

PubMed Abstract | CrossRef Full Text | Google Scholar

Yue, X., Gao, X.-Q., Wang, F., Dong, Y., Li, X., and Zhang, X. S. (2014). Transcriptional evidence for inferred pattern of pollen tube-stigma metabolic coupling during pollination. PLoS ONE 9:e107046. doi: 10.1371/journal.pone.0107046

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, C., Liu, J., Zhang, Y., Cai, X., Gong, P., Zhang, J., et al. (2011). Overexpression of SlGMEs leads to ascorbate accumulation with enhanced oxidative stress, cold, and salt tolerance in tomato. Plant Cell Rep. 30, 389–398. doi: 10.1007/s00299-010-0939-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, G.-Y., Liu, R.-R., Zhang, C.-Q., Tang, K.-X., Sun, M.-F., Yan, G.-H., et al. (2015). Manipulation of the rice L-galactose pathway: evaluation of the effects of transgene overexpression on ascorbate accumulation and abiotic stress tolerance. PLoS ONE 10:e0125870. doi: 10.1371/journal.pone.0125870

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, Y., Tao, Q. C., Wang, Z. N., Fan, R., Li, Y., Sun, X. F., et al. (2012). Engineering ascorbic acid biosynthetic pathway in Arabidopsis leaves by single and double gene transformation. Biol. Plant. 56, 451–457. doi: 10.1007/s10535-012-0119-x

CrossRef Full Text | Google Scholar