1 Department of Biochemistry, University of Cambridge, Cambridge CB2 1QW, United Kingdom
2 MRC Clinical Sciences Centre, Faculty of Medicine, Imperial College London, Hammersmith Hospital, London W12 ONN, United Kingdom
3 Metabolic Research Laboratory, Oxford Centre for Diabetes, Endocrinology and Metabolism, University of Oxford, Churchill Hospital, Headington, Oxford OX3 9JL, United Kingdom
4 Biomedical Sciences, Faculty of Medicine, Imperial College London, Exhibition Road, South Kensington, London, United Kingdom
5 Genetics and Genomics Research Institute, Faculty of Medicine, Imperial College London, South Kensington, Armstrong Road, London SW7 2AZ, United Kingdom
![]() |
ABSTRACT |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
transcriptomics; metabolomics/metabonomics; DNA microarray; system biology; bioinformatics
![]() |
INTRODUCTION |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
The analysis of multivariate data sets that are generated within these different tiers of biological organization may potentially yield critical information about the mechanisms involved in biological transitions. However, enhanced understanding of complex biological systems requires the development of strategies for integration of the different tiers, which link information derived at each level in a systematic way (24, 34). In particular in complex biological systems where large data sets will be the routine, there is a need to integrate the different functional genomic approaches, taking into consideration the scaling and correlation problems associated with this integration.
In this study a simple perturbation of a biological system of relevance to our sphere of interest is investigated in terms of an integrative approach to explain the linkage between changes in gene expression and associated changes in metabolic pathways identified by NMR. Profound fatty liver can be induced in certain vertebrates by the feeding of orotic acid, a naturally occurring biological precursor of uridine (19, 25). This phenotype is maintained without obvious deleterious effect to the animal throughout extended periods of orotic acid feeding and is manifest over a period of several days in liver cells taken from orotic acid fed rodents and grown in primary culture in the absence of orotic acid (15, 20). Thus the hepatic phenotype should be tractable to functional genomic analysis.
We have used a two-tier strategy to define the metabolic and gene regulatory responses, by metabolic profiling (metabolomics/metabonomics) and DNA microarrays, created by administrating orotic acid in an inbred (Kyoto) and outbred (Wistar) strain of rat. As such, the system consisted of a simple perturbation but was subject to polygenic variation. The integration of the data from this two-tier system describes, in an understandable and biologically coherent direction, a series of relationships between nucleotide, fatty acid, glycerolipid, and cholesteryl ester biosynthesis, and of stress responses during orotic acid feeding. The approach also identified unexpected links, including an inverse correlation, between stearoyl-CoA desaturase 1 (SCD 1) mRNA levels and unsaturated fatty acids. We conclude that the integration of 1H-NMR and gene expression data can provide a robust approach to identify specific pathways and cellular responses underlying a changing phenotype and, in addition, yield new insights into the mechanisms involved in the biological transition.
![]() |
MATERIALS AND METHODS |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
Metabolic profiles derived via 1H-NMR spectroscopy.
Solution-state 1H-NMR spectroscopy was performed on hepatic tissue extracts and blood plasma at 600.2 MHz using a Bruker AVANCE spectrometer interfaced with a 14.1-Tesla superconducting magnet and high-resolution inverse geometry 1H-NMR probe. Extracts of hepatic tissue (n = 3 for each time point and both strains) were prepared using an acetonitrile/water mixture as previously described (13) and reconstituted in 2H2O containing 4 mM trisilylpropionic acid (TSP). For the analysis of blood, plasma (n = 3 for each time point and both strains) was diluted twofold in a phosphate buffer containing 2H2O and 4 mM TSP. Solvent-suppressed spectra were acquired for tissue extracts and blood plasma into 32 k data points, averaged over 128 acquisitions. Spectral assignments were made with reference to previously published literature (2) and confirmed using two-dimensional pulse sequences.
Intact liver tissue (10 mg; n = 3 for each time point and both strains) was examined using a high-resolution magic angle spinning (HRMAS) probe placed in the spectrometer described above. Solid-state spectra were acquired by spinning the samples at 5,000 Hz using a solvent-suppression pulse sequence to reduce the spectral intensity of the water resonance. Spectra were averaged over 256 acquisitions.
Spectra were processed using XWINNMR software, version 3.1 (Bruker, Germany). Following multiplication with a preexponential factor of 1 Hz and Fourier transformation, spectra were integrated across 0.04 ppm spectral regions between 0.4 and either 4.2 or 9.4 ppm for solid- and solution-state NMR, respectively, using the AMIX software package (Bruker). The different integral regions used were a result of the rolling baseline produced in some solid-state spectra caused by unresolved protein resonances. For solution-state spectra the region around the water resonance was excluded to minimize the effect of variable water suppression in the subsequent pattern recognition. Output vectors representing each spectrum were normalized across the integral regions so that each individual 0.04-ppm region was represented as a ratio to the entire spectral intensity, excluding the water resonance (13). While this produces data sets with inherent redundancy due to metabolites with multiple integral regions at different chemical shifts, this redundancy confirms spectral assignments. Metabolite ratios were derived by manual integration of the resonances specified.
Lipid profiling by gas chromatography.
Livers (n = 3 for each time point and both strains) were homogenized using 5 vol (wt/vol) of 0.25 M sucrose. Following solvent extraction of the total lipid fraction, the triglyceride and cholesterol contents were determined using commercial kits (Boehringer-Mannheim, Lewes, UK). Fatty acid profiles of the total hepatic lipid content were obtained by gas chromatography following transesterification using a mixture of toluene/sulfuric acid in methanol (11).
Gene expression profiling using microarrays.
Total RNA was extracted by an RNA isolation kit (Stratagene) from the livers of Wistar and Kyoto rats at day 0 (Wistar n = 3, Kyoto n = 2), days 1 and 3 (Wistar n = 2, Kyoto n = 2), and day 14 (Wistar n = 3, Kyoto n = 2). cDNA was generated from total RNA using the SuperScript double-stranded cDNA synthesis kit (Invitrogen, Life Technologies), and cRNA was prepared by in vitro transcription (BioArray High Yield RNA transcript labeling kit; Affymetrix UK) and purified on an affinity resin column (RNeasy, Qiagen, UK). Biotinylated cRNA was fragmented and hybridized to RG-U34A Affymetrix GeneChips (Affymetrix, Santa Clara, CA; http://www.affymetrix.com/support/technical/datasheets/rgu34_datasheet.pdf). Data were analyzed with the Affymetrix Microarray Suite 5.0 algorithm to generate a "signal" value and a "detection" label. Prior to normalization, probes that generated an absent call in all cells were removed. Normalization was calculated using a global scaling method based on the entire data set. Globally scaled values <10 were rounded up to 10. The data set is deposited as a series at the Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo) and has the accession number GSE665.
Pattern recognition and data processing.
All data are represented as means ± SD unless otherwise specified. The 1H-NMR spectra and transcriptional data sets were imported into the SIMCA package (Umetrics, Umea, Sweden) and preprocessed by three different scaling processes (9): 1) mean centering by measuring the variance of each integral region about the mean; 2) univariate using mean centering followed by scaling the variance by 1/sk, where s is the standard deviation of the variable (integral region) k; and 3) Pareto scaled to (1/sk)1/2. Pareto-scaled data sets produced the most robust models and are reported in this paper. Initially each data set was examined using principal component analysis (PCA) to examine trends and clusterings. Time progression of the data sets were interrogated by prediction to latent structures through partial least squares (PLS) using each integral region as an X-variable and time as the Y-variable. PLS calculates the linear relationship between a matrix Y (in this case time) and a matrix X (integral regions representing metabolites or transcripts). This may be represented as Y = f(X) + E, where E represents the residuals not modeled.
Three types of preprocessing were applied to the X-matrix of integral values. Correlations were tested as part of the routine by sequentially leaving out every seventh sample and predicting its class membership. This goodness of fit algorithm was used to determine whether a correlation was significant (Q2 > 0.097) and was derived from the predicted error sum of squares between the actual and predicted Y-matrix (PRESS) and the sum of squares proportion of the X-matrix explained by the model (SS) according to: Q2 = 1 PRESS/SS. While PLS models two data matrices in terms of linear regression, nonlinear responses may be modeled using a polynomial relationship or produces curvature in a linear model between PLS components representing the X- and Y-block data sets.
Transcriptional data was also ordered using hierarchical clustering with a variety of algorithms, including Euclidean distance, Standardized Euclidean distance (each coordinate in the sum of squares is inversely weighted by the sample variance of that coordinate), Mahalanobis distance, City Block/Manhattan distance, or Minkowski distance according to temporal and strain profiles within the program Matlab (http://www.mathworks.com). Hierarchical clustering plots were cross-checked using the cophenetic correlation coefficient between the original distance data and the clusterings formed as described in the statistical toolbox of Matlab.
Genomic and metabolic correlations.
To correlate differentially expressed transcripts with changes in metabolites, two matrices were formed for PLS regression within SIMCA (Umetrics) with 1H-NMR spectra forming the X-block and transcriptional changes the Y-block. Both strains of rats were included in the same model, but as separate genomic/metabolic inputs. The PLS model mapped each observation in terms of metabolic and transcriptional variation to identify those animals where these changes were maximal and whether there was a trend in these changes in terms of time or rat strain. The variables most significant to the regression model were identified using the loadings plot derived from the PLS model, and these were used to identify transcriptional and metabolic co-responses. For tissue taken from animals at days 1, 3, and 14 with no corresponding gene array data, an average was taken for that rat strain and exposure period.
![]() |
RESULTS |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
|
To substantiate the differences in the hepatic lipid profiles of rats fed orotic acid, we used gas chromatography analysis on liver tissue extracts. The total amount and proportion of lipid associated with 16:0, 16:1, 18:0, and 18:1 fatty acids increased with exposure to orotic acid. Importantly, there was a more marked increase in 16:1 and 18:1 fatty acid lipids, along with a proportional decreases in 18:0 fatty acids (Table 1). Thus the increased ratio of unsaturated to saturated lipids in the cytoplasm detected by HRMAS 1H-NMR spectroscopy was also reflected in the total lipid content.
|
Orotic acid feeding also induced systemic changes in blood metabolite concentrations in both strains (Wistar rats are shown Fig. 1E). By TSP there were marked increases in plasma concentrations of ß-hydroxybutyrate, which was undetectable in plasma from control rats. Lactate, acetate (not shown), and glucose (not shown) were also increased. There was also the expected decrease in the intensity of the LDL and VLDL signals (compound resonances at 1.0 and 1.3 ppm for both CH2 and CH3) and a decrease in phosphatidylcholine.
A feature of the orotic acid fed rat is that the phenotype completely corrects with concomitant adenosine feeding (19). In the present study, the metabolic profiles of intact liver, liver extracts, and plasma from control animals and Wistar rats fed orotic acid plus adenosine for 14 days were indistinguishable (data not shown).
Overall orotic acid feeding increased hepatic cytoplasmic triglyceride accumulation, and this was more marked in the predisposed Kyoto rat strain than the Wistar strain. In both strains there were also prominent increases in total cholesterol, phosphocholine, phosphatidylcholine, and unsaturated fatty acids compared with saturated fatty acids in terms of the total lipid content of the tissue. The methyl donors, betaine and TMAO, and pyrimidine nucleotides were increased in both strains and were particularly marked in the predisposed Kyoto strain.
Gene expression profiling during orotic acid feeding.
The gene expression profiles of the two rat strains at 1, 3, and 14 days were assessed relative to control animals using the Affymetrix RG-U34A GeneChip, containing 7,000 known transcripts and 1,000 expressed sequence tags (ESTs). Of these 8,000 transcripts, 4,132 expression elements were expressed in the microarray experiment.
We further investigated the gene microarray data set using a simple statistical filter to identify the principal genes involved in the orotic acid feeding perturbation in both strains of rats compared with the pre-dose animals. We argued that by combining the transcription data sets that genes changes common to both strains compared with control animals would represent those most affected by the orotic acid treatment. Although these transcript changes were at different levels in the different strains, combining the gene expression data from both strains increased the power for detecting genes affected by orotic acid feeding, while minimizing the number of false possibilities in the data set.
Using a threefold change in conjunction with statistical support of P < 0.05 for a change in transcription as a filter for the data, we identified the expression of 48 gene sequences as being altered by orotic acid feeding over the time course (Fig. 2A) (17). At a P value of <0.001, there was also support for orotic acid having an effect (i.e., between 1.5-fold and 3-fold expression change) on the expression of 12 further sequences. This analysis demonstrated that many of the transcripts identified by multivariate analysis did not pass through this more stringent reproducibility test, a problem associated with the much larger transcriptional data sets compared with those derived from 1H-NMR spectra. Also, in the case of ATP citrate lyase only one of the four transcripts represented on the microarray GeneChip passed the filter as being significantly increased, highlighting the possibility of identifying false positives in microarray experiments. However, as microarray experiments are also affected by false negatives, ATP citrate lyase was kept in the data set for subsequent analyses. ATP citrate lyase was the only gene highlighted in this manner. This transcriptional data set of 60 transcripts was used in subsequent correlation analyses, but the different strains were considered independently to examine how these common transcriptional changes differed in the two rat strains.
|
Hierarchical clustering was used to organize the 60 differentially expressed gene transcripts according to temporal and strain profiles (Fig. 2B for Manhattan distance method). All of the plots were robust, as judged by cophenetic coefficient values >0.88, and as expected, transcripts from the same gene, or genes with similar function, clustered together. However, the majority the 60 transcripts could not be separated into distinctly regulated subgroups and thus provided minimal information. Despite this, the transcripts consistently perturbed in the two strains indicate the areas of metabolism affected by orotic acid and proteins on the metabolic pathways involved. These results are highly plausible biologically and consistent in their direction of change with the changes in metabolites. This is with the exception of the increase in unsaturated fatty acids induced by orotic acid and the decrease in the mRNA for SCD 1, which appear to move in counterintuitive directions. The metabolite and transcript changes are illustrated on Fig. 3.
|
|
|
Correlation of transcripts to metabolic changes using bootstrapping.
To further investigate the regression models of transcriptional and metabolic changes and test their robustness, correlations between transcripts and key metabolic changes were investigated using a boot strapping routine (Matlab, Mathworks). Correlation coefficients were calculated for all 60 transcripts and the concentration of TMAO, the metabolite that demonstrated by far the largest change of all in the aqueous extracts, using bootstrapping with 10,000 iterations. This sampling with replacement routine allowed quantification of the reproducibility of these correlation coefficients for the transcripts.
The most positively correlated transcript with TMAO was GRP94, while the most negatively correlated was SCD 1 (Fig. 5). SCD 1 was also negatively correlated with CH2CH2CH2 and CHCH lipids, mitochondrial very-long-chain acyl CoA thioesterase, and GPAT transcripts, and positively correlated with apoCIII. Accordingly, a negative correlation was observed between apoCIII and GPAT, which was positively correlated with fatty acid CoA ligase 4.
|
![]() |
DISCUSSION |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
The metabolic phenotypes derived by NMR spectroscopic analysis were used to define "normality" in the outbred genetically heterogeneous Wistar rat relative to the inbred Kyoto rat. The NMR hepatic lipid profile of Kyoto rats relative to Wistar rats prior to orotic acid exposure suggested that the strain might be predisposed to fatty liver. These differences were borne out in the dynamic changes induced by orotic acid in Kyoto rats, with a much more marked increase in mobile lipids. Thus the increased susceptibility of Kyoto rats to fatty liver disease becomes overt following the challenge imposed by orotic acid. The reduced levels of glucose and glycogen in the livers of the Kyoto rats also suggest that this strain may be susceptible to insulin resistance, sharing some of the genetic features of the insulin-resistant Kyoto-derived spontaneously hypertensive rat (1). This reverse functional genomic approach using NMR spectra may have broad implications for the identification of genetic predispositions to other pathological states in pharmacogenomics.
In light of the observation that adenosine supplementation in conjunction with orotic acid prevents and reverses fatty liver (19), a crucial finding is the increased NMR signal for phosphatidylcholine, the methyl donor betaine, and the choline oxidation product TMAO. Phosphatidylcholine can be produced by CTP- and ATP-dependent pathways. The ATP-dependent pathway involves the conversion of phosphatidylserine to phosphatidylethanolamine and in turn to phosphatidylcholine and may proceed through the activation of L-methionine with ATP via the activity of L-methionine adenosyltransferase. Methionine for this process is regenerated from homocysteine by methyl donation from betaine, suggesting a link between the changes detected in methyl donors and phospholipid fatty acid metabolism of rats fed orotic acid. These changes may also link into the reduced secretion of hepatic apoB-containing lipoproteins from these animals, since a previous study has reported that the overexpression of betaine homocysteine S-methyltransferase, which catalyzes the donation of methyl groups from betaine to homocysteine, enhances the secretion of apoB-containing lipoproteins from the liver cell line, McArdle RH-7777 (29). Conversely, the secretion of VLDL triglyceride is reduced in mice deficient for phosphatidylethanolamine N-methyltransferase (22). Furthermore, we also note betaine and TMAO are oxidation products of choline, and choline deficiency induces fatty liver disease through deficiency of methyl donors (6, 28). In addition, methionine adenosyltransferase 1A knockout mice are more susceptible than wild-type animals to choline-deficient diet-induced fatty liver disease and express many acute phase-response and inflammatory markers, including metallothionein and growth-related genes (18) as seen in this study.
Several genes in the lipogenic pathway are regulated by sterol regulatory element binding protein-1c (SREBP-1c) including FAS, GPAT, and SCD 1. The decrease in the expression of SCD 1 by orotic acid suggests decreased transcriptional activation by SREBP-1c and is also consistent with previous reports demonstrating reduced de novo fatty acid synthesis in rats fed orotic acid (30, 31). SREBP-1c-mediated downregulation of SCD 1 and of lipogenic flux could relate to the large increase in the concentration of 16:1 and especially 18:1 fatty acids in the liver following orotic acid feeding (16). Unsaturated fatty acids such as these have been shown to decrease SREBP-1c mRNA levels by two mechanisms. First, liver X-receptor- (LXR
) is a transcriptional regulator of SREBP-1c (3, 16. 23), and activation of LXR
by its natural ligands is antagonized by unsaturated fatty acids (32). Second, unsaturated fatty acids accelerate SREBP-1c transcript decay (32). Although these SREBP-1c-related mechanisms would explain the decreased SCD 1 expression and decreased lipogenic flux following orotic acid feeding, our data indicate that the expression of other SREBP-1c target genes such as FAS and GPAT are not decreased under these conditions. This is a paradoxical finding that suggests an orotic acid-mediated uncoupling of the normal coordination of the expression of genes involved in regulation of the lipogenic pathway.
An early and dramatic response to orotic acid feeding is the increased transcription of chaperone and stress response genes and of genes involved in growth and inflammation control. GRP94 and the PDI-related protein are endoplasmic reticulum (ER)-resident proteins involved in the folding of apoB during the assembly of the triglyceride-rich lipoproteins in VLDL (21, 33). Surprisingly, the apoB gene is consistently induced, perhaps as a response to the accumulation of lipid in ER vesicles with orotic acid feeding. The metallothioneins and fibronectin are important early response general stress genes, which respond to inflammation and oxidative stress (4, 7, 8). These changes are accompanied by decreased activity of growth control genes such as TGFß inducible early growth response, early growth response 1, cyclin D3 and inflammation regulating genes like B lymphocyte chemoattractant precursor (weakly similar), and p38 mitogen activated protein kinase.
One caveat with interpreting the potential pathophysiological implications of the data generated is that the NMR and microarray analyses were done on whole liver extracts or whole tissue samples. Thus it was not possible to determine which cell types, for example, hepatocytes, macrophages, or stellate cells, were being influenced directly. Previous histopathological studies indicate that orotic acid largely acts directly on hepatocytes (19, 20, 25). However, improvements both in NMR spectroscopy miniaturization and amplification techniques for mRNA transcripts suggest that in the future such approaches may be coupled to laser-capture dissection of specific cell types.
The reverse functional genomics strategy described in the present paper demonstrates that metabolic phenotypes can be correlated with transcriptional data from DNA microarray experiments and used to define critical metabolic pathways, and importantly this approach can accommodate information from different animal strains across a time series. Although in the present study we have not evaluated posttranslational regulation, allosteric controls, or flux through pathways, we have been able to obtain important new information on the orotic acid-induced fatty liver phenotype through the integration of transcriptional and metabolic changes. Multivariate analysis of the combined transcription and NMR spectroscopy data sets were well suited to the large matrices produced and could be used as a strategy for investigating toxicological insults and the implications of pharmacogenomic predisposition to therapeutic interventions.
![]() |
GRANTS |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
![]() |
ACKNOWLEDGMENTS |
---|
Present address of C. Mann: Genetix, Queensway, New Milton, Hampshire, BH25 5NN, UK.
![]() |
FOOTNOTES |
---|
Addresses for reprint requests and other correspondence: J. L. Griffin, Department of Biochemistry, University of Cambridge, Cambridge CB2 1QW, UK (E-mail: jlg40{at}mole.bio.cam.ac.uk); and J. Scott, Genetics and Genomics Research Institute, Faculty of Medicine, Imperial College London, South Kensington, Armstrong Road, London SW7 2AZ, UK (E-mail: j.scott{at}imperial.ac.uk).
10.1152/physiolgenomics.00158.2003.
![]() |
REFERENCES |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|