The initial metabolites resulting from sphingolipid biosynthesis include the sphingoid bases dihydrosphingosine (DHS) and phytosphingosine (PHS), which in turn serve as metabolic precursors for synthesis of an array of chemically diverse species including ceramides and sphingoid base phosphates (). Though the homologous mammalian lipid sphingosine-1-phosphate activates a variety of signaling events through receptor-mediated modulation of characterized signaling pathways (Alvarez et al, 2007
), specific cell functions for sphingoid base phosphates in yeast remain unknown. PHS1P was shown earlier to transiently increase during heat stress (Skrzypek et al, 1999
), suggesting a potential role for this lipid in the heat stress response. Moreover, mutant yeast strains that accumulate PHS1P exhibited poor growth and heat stress resistance (Skrzypek et al, 1999
; Kim et al, 2000
), though a specific role for this lipid in the heat stress response remains unidentified. We hypothesized that PHS1P may mediate subprograms of the heat stress response, and, therefore, we designed experiments to perturb PHS1P metabolism using a conventional gene deletion approach. Deletion of LCB4 and LCB5 genes, encoding the yeast sphingoid base kinases (Nagiec et al, 1998
), attenuates PHS1P production whereas deletion of DPL1, encoding the sphingoid base phosphate lyase (Saba et al, 1997
), blocks its degradation and results in accumulation of PHS1P (). Thus, mutating the above genes allows a ‘clamp' on PHS1P in cells at either low or high levels, respectively. Log-phase cultures of the lcb4
Δ, the dpl1
Δ, and wt
strains were subjected to heat stress, and samples were collected at 5-min intervals over a time course of 30 min in two duplicate experiments. To collect multiple ‘-omics' data reflecting cellular signaling systems, each sample was divided into two portions, thus allowing mRNA and sphingolipid extraction from each sample for transcriptomic and lipidomic analyses (Bielawski et al, 2006
), generating a total of 42 transcriptomics measurements paired with 42 corresponding lipidomics measurements.
Summary of major sphingolipid biosynthetic pathways in Saccharomyces cerevisiae.
Transcriptomic responses to heat stress
As previously reported (Gasch et al, 2000
; Gasch and Werner-Washburne, 2002
; Cowart et al, 2003
), heat stress significantly induced or repressed the expression of over a thousand genes in the wt
strain, and a significant portion of this program was conserved in both the lcb4
Δ and the dpl1
Δ strains. However, we further identified the genes that are differentially expressed among the strains under comparable conditions. This group included 687 probe sets (corresponding to 441 genes with gene names) that were differentially expressed after heat stress in the lcb4
Δ mutant when compared to wt
strains. Analysis of this group using GOStat (Beissbarth and Speed, 2004
) revealed a broad range of functional categories including DNA repair, translational regulation, post-translational modification of proteins, cell wall organization and biogenesis, and others ().
Figure 3 Effects of mutations in specific sphingolipid metabolic genes on gene expression and total sphingolipid profiles. (A) Overrepresented Gene Ontology annotations for genes aberrantly regulated during heat stress in the lcb4Δ/lcb5Δ mutant (more ...)
Lipidomics responses to heat stress
Recognizing the connectedness of sphingolipid metabolism (Alvarez-Vasquez et al, 2005
; Hannun and Obeid, 2008a
), we hypothesized that the mutations might cause broad changes in sphingolipid metabolism and that only a subset of these genes spanning broad functional categories actually resulted from the lack of PHS1P, whereas many of these genes' aberrant regulation might result from modulation of other sphingolipids in this mutant. Therefore, we further evaluated the lipidomic data to determine the impact of heat stress on lipid species and to determine whether the effects of mutations generate a metabolic ‘ripple effect' leading to changes in other metabolites. Indeed, out of 40 distinct sphingolipid species measured, 28 demonstrated changes at least at one time point of heat stress in wild-type cells (). Heat stress induced an increase in C18 PHS1P from nearly undetectable levels to 0.02 pmol/nmol phosphate (see Supplementary Table 1
), peaking around 20 min and returning to basal levels by 30 min. These kinetics were similar to previously published data that demonstrated an eight-fold increase in PHS1P that peaked around 10 min and returned to basal levels around 20 min (Skrzypek et al, 1999
). As expected, these changes did not occur in the lcb4
Δ mutant strain. On the other hand, the dpl1
Δ strain showed constitutively elevated PHS1P, which further increased during heat stress (, bottom row). More importantly, the data revealed widespread differences in lipid profiles between the three strains, both basal and under heat stress. For example, the lcb4
Δ strain demonstrated significant elevation of PHS, increased α-hydroxylated ceramides of 20 and 22-carbon N
-acyl chain length, increased phytoceramides of 24 carbon chain length, and decreased dihydroceramide of 26 carbon chain length. In fact, out of the 28 lipid species showing changes in the wild type over the time course, at least 19 species showed differences between wt
and the lcb4
Δ strain at least one time point.
Lipidomics profiles of the dpl1Δ strain also demonstrated widespread differences as compared to the parental strain (, right panel). With the exception of the expected accumulation of PHS1P, lipid measurements in this strain revealed fewer and more subtle variation from the wild-type strain at basal temperature; however, time course measurements indicated that deletion of DPL1 significantly altered the heat stress sphingolipid response in that at least 20 of the 28 measured lipids exhibited differences from the parental strain in at least one point of the time course. Moreover, changes in this mutant were partially distinct from changes observed in the lcb4Δ/lcb5Δ mutant strain.
In summary, these mutations not only had the expected effects on PHS and PHS1P (substrate and product), but also caused widespread changes in many sphingolipid species. Therefore, any changes in gene regulation observed in the lcb4Δ/lcb5Δ or dpl1Δ mutant strains could not be readily attributed to PHS1P (or any single lipid species). To circumvent the limitations of ‘gene-centric' approach, an alternative approach was devised to integrate information from lipidomic and transcriptomic data in a manner that allows inferring lipid-mediated events.
Indentify potential sphingolipid-regulated genes
To identify and quantify the information connecting lipidomic and transcriptomic changes, we performed covariance analysis (DeGroot and Schervish, 2002
) between all lipid-probe pairs, which led to a genes-
matrix in which an element contained a correlation coefficient (ρ) for a lipid–gene pair if the ρ is statistically significant or a zero otherwise. It was noted that a row of the matrix constituted a correlation (information) pattern demonstrated by a gene with respect to all lipids, and similarly a column encoded an information pattern demonstrated by a lipid with respect to all genes. We sought to identify the shared information patterns among genes and lipids by applying a double-sided hierarchical clustering analysis, which grouped genes (and lipids) sharing similar information patterns into clusters. Selected results focusing on key sphingolipids are shown in as a heat map (also see Supplementary Table 2
). The results indicate that lipid species demonstrated distinct correlation patterns with respect to modules of genes, suggesting potential regulatory roles of the lipids on the modules. The double-sided clustering revealed clusters of genes sharing similar information with respect to lipid (blocks across rows) and clusters of lipids sharing similar information with respect to gene expression data (columns with similar correlation coefficient patterns). For example, DHS and PHS, two closely related metabolites in the metabolic network, were grouped together because of their shared information with respect to clusters of genes at the top region , but they also showed distinct information with other gene clusters.
Importantly, this procedure identified subsets of genes that were significantly correlated to PHS1P, with 44 positively correlated probe sets (mapped to 23 named genes) from the microarrays and 61 negatively correlated ones (mapped to named 54 genes), which was significantly less than the genes identified in the microarray analysis of the lcb4Δ/lcb5Δ mutant (441 genes). Among these 77 PHS1P-sensitive genes, 40 genes were also deemed differentially expressed according to differential expression analysis, whereas the other 33 genes were not. The results indicate a tentative statistical advantage of correlation analysis over the differential expression analysis in that the former can use all samples (42 samples for a lipid-versus-gene pair) whereas the latter can only use samples from specific conditions (12 lcb4Δ/lcb5Δ microarrays versus 12 wt microarrays after heat stress).
We applied a method referred to as GO Steiner Tree (GOSteiner, 2009
) to analyze the functional coherence of the gene sets and to visualize their functional relationships. The method represents the genes and their associated Gene Ontology terms as a graph, finds a subgraph (a Steiner tree) connecting all genes and their annotations with a shortest total functional semantic distance, and finally evaluates the statistical properties of the tree as metrics of functional coherence. Supplementary Figure 1
shows the GO Steiner tree for genes that demonstrated failure to induce on heat stress in the lcb4
Δ strain when compared to wt
. Although the analysis showed that the lcb4
Δ-sensitive genes had diverse functions, the visualization of the GO Steiner tree revealed clusters of genes with coherently related functions, including genes involved in mitochondrial metabolism, for example, COX4, INH1, and ATP17, and vesicular transport, for example, APS2, SVP26, TPM2, and TVP23. The results indicate that the deletion of LCB4 and LCB5 led to aberrant regulation of several distinct groups of genes, among which genes showing high correlation to PHS1P represented only a subset of those genes (highlighted in Supplementary Figure 1
). Importantly, this subset included genes in specific sub-categories (e.g. the mitochondrial metabolism), but not others (e.g. vesicular transport) of the LCB4/LCB5-regulated genes.
The above results led to the hypothesis that, although both sets of genes were lcb4Δ/lcb5Δ sensitive, only the genes showing strong correlation to PHS1P levels were truly PHS1P sensitive, whereas dysregulation of other genes in this mutant strain could result from confounding changes in other sphingolipids which had become apparent from the lipidomics analysis (). We tested this hypothesis by treating wild type cells (wt) with exogenous PHS1P in at non-heat stress temperature and monitored the expression of sample genes from the putative PHS1P-dependent set involved in the mitochondrial metabolism (COX4, INH1, and ATP17) and a putative non-PHS1P-dependent set involved in vesicular transport (APS2, SAR1, SVP26, TPM2, and TVP23). Indeed, expression of the genes involved in vesicular transport failed to be induced by PHS1P (), whereas the mitochondrial metabolism genes demonstrated ~2.4–3-fold increases on treatment (). To further test whether the effects were specific for PHS1P, wt and lcb4Δ/lcb5Δ mutant cells were treated with PHS, the metabolic precursor of PHS1P, which undergoes conversion to PHS1P in the wild type but not in the lcb4Δ/lcb5Δ mutant. Indeed, in wild-type cells, PHS produced a significant upregulation of the genes in the respiration set, which was totally lacking in the lcb4Δ/lcb5Δ mutant (). Thus, the integrated analysis enabled the identification of the PHS1P-dependent subset of genes within the larger set of genes that showed failed regulation in the lcb4Δ/lcb5Δ mutant.
Figure 4 PHS1P-mediated gene expression regulation. (A) Treatment with PHS1P in the absence of heat stress did not induce all genes aberrantly regulated in the lcb4Δ/lcb5Δ strain. (B) Treatment with PHS1P in the absence of heat stress induces gene (more ...)
Identifying candidate TFs
Defining specific lipid-mediated responses, beyond the overall gene-mediated responses manifested by mutations, is important in that it allows studying the distinct roles and mechanisms of specific lipids—the putative functional mediators—that contribute to the overall responses. Highly selective PHS1P-mediated expression of specific genes led to the hypothesis that PHS1P regulates a specific pathway by activating some downstream TFs, which then mediate the regulation of these genes in response to PHS1P. To identify such putative TFs, a transcription factor binding site (TFBS) matrix was constructed, in which an element contains a binary variable indicating if a gene (g
) can be bound by a TF (t
). AN element, btg
, of the matrix is set to 1 if the analysis of the chromatin immunoprecipitation experiments (Lee et al, 2002
; MacIsaac et al, 2006
) indicates that the TF t
is capable of binding to the promoter sequence of gene g
, or if there is documentations of the interaction between the TF and the gene according to the yeast TF database, YEASTRACT (Monteiro et al, 2008
). The knowledge of TFBSs enabled us to evaluate whether any TFBS was significantly enriched in the promoters of a gene set. Using a hypergeometric-distribution-based model, we assessed enrichment of TFBSs in the promoters of the genes that demonstrated positive correlation to PHS1P and identified 22 TFBSs as significantly ‘enriched' (see Supplementary information
). This led to the hypothesis that the activation states of some of the candidate TFs were sensitive to changes in PHS1P, thus transmitting the signal from PHS1P to genes. To test the hypothesis, we further developed a two-stage Bayesian information integration model to reveal information flow from bioactive lipids → TFs → gene expression.
Inferring activation states of TFs
We developed a novel Bayesian model to infer the activation states of the TFs under each experimental condition. This model extends our previous published method (Lu et al, 2004
) such that the genomic data of TFBSs can be integrated with transcriptomic data to infer the activation state of TFs. In this model, the activation state of a TF (t
) under a specific condition (a
) is represented as a binary variable, sat
, such that sat
=1 indicates the TF is active and sat
=0 otherwise. Representing activation states of TFs as binary variables has two advantages: (1) it provides a more intuitive representation of activation (active/inactive) state of a TF, in comparison to possible negative activation state allowed by some other models (Lee and Batzoglou, 2003
; Liao et al, 2003
; Gao et al, 2004
; Ochs et al, 2004
; Battle et al, 2005
; Sun et al, 2006
); (2) it renders the mathematical convenience to model the relationship between lipids and activation states of TFs using logistical regression, in which the sigmoid function mimics dose–response curves commonly observed in signal transduction pathways.
The model specifies that the expression value of a gene from a specific experiment (ega
) is influenced by three factors: (1) TFs that bind to its promoter, indicated by bgt
}; (2) the states of each of the TFs under this specific condition, represented by sat
}; (3) and the strength and direction (induction or repression) of an activated TF on its expression, represented by wgt
}. We define the probabilistic relationship between the above parent variables and the gene expression value as follows:
and τ represent the noise of the system and N
stands for the Gaussian distribution. It is of interest to note that the product of two binary variables, bgtsat
, encodes a logic AND relationship between the two variables in the equation, such that the equation can be interpreted as follows: TF t
influences the expression value of gene g
under the condition a
if and only if it has a binding site in the promoter of the gene (bgt
=1) AND it is activated under the condition (sat
=1). The equation also reflects coordinated influences of multiple activated TFs on a gene's expression in a linear form (usually in logarithmic scales), an assumption widely used in modeling of expression systems (Liao et al, 2003
; Gao et al, 2004
; Lu et al, 2004
; Battle et al, 2005
). Given TFBS matrix and expression data, the model probabilistically infers the state of each TF under a specific condition using a variational Bayes technique (Lu et al, 2004
) (see Supplementary information
for detailed description). shows that many TFs switch states during heat stress, among which many are well-documented stress-responding TFs.
Modeling the role of lipids in TF activation
On the basis of the estimated binary TF states from the previous section, a Bayesian logistic regression model was applied to investigate the relationship between the levels of sphingolipids and the putative TF activation states. In this model, the states of a TF under each experimental condition were modeled as a sigmoid function of the concentrations of sphingolipids, where the probability that a TF (t
) is active, conditioning on observed lipid profiles (l1,
) under an experiment condition (a
), is defined as follows:
. In this model, a parameter βlt
reflects the strength and direction of influence that the sphingolipid l
has on the activation state of the TF t
. Using a Gibb's sampling-based Bayesian logistic regression, we identified all statistically significant parameters, and we interpret a significant parameter βlt
with respect to a TF t
as an indication that the TF is regulated by lipid l
, and the selected results are shown in . A total of 13 TFs showed with significant coefficient with respect to PHS1P, with 11 positively influenced and 2 negatively influenced. Among the positive TFs, Hap2p, Hap4p, NRG1 and MGA1 were also among the 22 TFs deemed significantly enriched in the PHS1P postively correlated gene set. Hap2p and Hap4p are members of the HAP complex (Chodosh et al, 1988
; Buschlen et al, 2003
), whose binding sites were dominant in the PHS1P postively correlated gene set, particularly in the subset involved in the cellular respiration. Thus, modeling of the relationship between sphingolipids and TF states, in combination with the static information inferred from promoter analysis, provided evidence at a mechanistic level that gave rises to the hypothesis that PHS1P regulates expression of respiratory genes through modulating the activity of the HAP complex.
To test this hypothesis, we evaluated the ability of PHS1P to regulate PHS1P-dependent genes in a mutant strain deleted for the gene encoding Hap4p (hap4Δ cells), an essential component of the HAP complex. Wt cells or hap4Δ cells were treated with exogenous PHS1P added into the culture media. RNA was isolated from each culture, and gene expression was determined by real time RT–PCR. Indeed, as compared with the parental background strain, the hap4Δ strain demonstrated complete loss of PHS1P-mediated induction of the target genes (). The result indicated that regulation of these genes by PHS1P required a functional HAP complex.