Global temporal changes in hepatic transcriptome during the 16-week time-course

To assess temporal changes in hepatic gene expression over the 16-week period under the control and two high-fat diets, each of the time points per diet was compared to time-point day 0 in a pairwise fashion using limma statistical package [27]. Applying statistical cutoff of false discovery rate (FDR) <0.1, we identified 839, 3027 and 3316 genes differentially expressed by chow, HFBT and HFP feeding, respectively, across any of eight time-points (Figure S1). This showed that ageing of animals from 14 weeks (at the beginning of the experiment) to 30 weeks (at its end) affected expression of a portion of genes, also in the control condition. Nevertheless, the combined effects of ageing and high-fat feeding observed in HFBT and HFP conditions was substantially larger compared to solely aging effects.

In addition to overall temporal effects, comparing each of the time points to day 0 allowed assessment of dynamics of transcriptional response by detecting the magnitude of the gene expression changes at each time point compared to the starting condition. In both HF conditions we observed three phases of hepatic transcriptional response characterized by local peaks in the number of differentially expressed genes: early (day 1 to week 1), mid (week 2 and week 4) and late (week 8 to week 16) (Figure S1). In addition to the phasic trend, the number of differentially expressed genes in HF conditions increased gradually during the course of the experiment. This may be expected considering that using day 0 as a common reference leads to combined contribution of both high-fat intervention and ageing to the observed effects. In chow condition, temporal trends identified for HF diets were absent and the numbers of changed genes in each of the time-points remained below the ones identified in HF conditions (Figure S1).

Hepatic transcription response to high-fat diets

After determination of the temporal effects of all three diets across 16-week time-course, we focused on dissecting the HF-specific effects on the hepatic transcriptome. To analyse the effects of HF diets independently of changes occurring due to the animal age, we compared the gene expression of mice fed HFBT and HFP diets to those of mice fed the chow diet in each of the corresponding time-points. The pairwise comparisons were performed using the limma statistical package, and FDR q-value<0.1 was used as a threshold for significance. The numbers of identified differentially expressed genes in each diet and time-point are shown in Table 1 and Table S1. Per time-point, the number of up- and down-reglated genes was largely balanced, and genes changed expression in the same direction in both HF conditions (Table S1). In total, during the 16-week time-course the HFBT diet significantly affected the expression of 1263 genes, while HFP diet affected 836 genes (8.3% and 5.5% of the monitored genome, respectively). The statistically significant differentially expressed genes (DEG) under the two HF diets largely overlap (436 genes, p-value 1.7E-265) (Figure 1B).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Number of statistically significant differentially expressed genes per time-point (FDR q-value<0.1).

https://doi.org/10.1371/journal.pone.0006646.t001

To compare gene expression patterns of identified HFBT and HFP DEGs across all conditions, expression profiles of genes that are changing in either of HF diets compared to chow diet per time-point (1663 genes in total) were hierarchically clustered (Figure 1C). The close proximity of condition tree branches corresponding to two diets in each of the time-points revealed a striking similarity between the effects of HFBT and HFP diets on the expression of all 1663 genes. Therefore, not only 436 overlapping genes, but also the remaining genes unique to each condition changed similarly under the HFBT and HFP diets, although they did not pass significance threshold in both conditions. Due to this highly comparable gene expression response, the union of differentially expressed genes under the HFBT and HFP diets was further considered as a total set of 1663 HF-responsive genes (Table S1, Table S2). This set of HF-responsive genes was used for all downstream functional analyses that required a limited set of genes of interest as input. Nevertheless, to be able to assess subtle differences between the effects of two HF diets caused by different fat origin and/or the presence or absence of cholesterol in a diet, gene expression changes are always displayed separately for each of two diets.

Three distinct phases (early, mid and late) of hepatic transcriptional response to HF diets were observed by comparing each time-point to day 0 (Figure S1). These phases are also revealed by hierarchical clustering of the experimental conditions in Figure 1C, where gene expression of mice fed HFBT and HFP diets was compared to those of mice fed the chow diet at each of the corresponding time-points. To further characterize the expression profiles of the 1663 HF-responsive genes, Smoothing Spline Clustering was performed [28] (Figure S2, S3, S4 and S5). The identified gene expression profiles typically show peaks at early (day 3), mid (week 2) and/or late (week 8/12) phase of the time-course. Especially prominent expression changes were observed at the early and the late phase, as determined by the intensity of the gene expression changes and the abundance of the significant differentially expressed genes (Table 1, Figure 1C, Table S1, Table S2). Interestingly, many genes change the direction of their regulation between the early and the late phase of the time-course (Figure 1C, Figure S2, Figure S3, Figure S4 and Figure S5). Similarly to the observation obtained by comparing each time-point to day 0, the peak-intervening time-points are characterized by a more modest transcription response and may reflect resumptions of the local homeostasis resulting from the transient adaptation to excess dietary fat.

To identify which cellular processes are most affected by the hepatic exposure to excess dietary fat over the entire time-course, we first analyzed the overrepresented functional categories among the 1663 HF-responsive genes [29]. The most prominent significantly enriched functional clusters are related to lipid, cholesterol and oxido-reductive metabolism, as well as inflammation, immune response, apoptosis, cell cycle, protein folding and the regulatory pathways controlling these processes (Figure 2, Table S3 (“All HF-responsive genes”)). The expression changes per diet and time-point of the selected representative genes for each functional cluster are also shown in Figure 2. To assess if some of the categories are specifically enriched in either up- or down-regulated genes, we additionally performed the equivalent type of analysis using separately lists of HF-responsive up- and downregulated genes (Table S3 (“Upregulated HF-responsive genes” and “Downregulated HF-responsive genes”)). As a number of genes changed the direction of their regulation between the early and the late phase of the time-course, many categories were represented in both lists. Specifically, the enrichment of lipid metabolism and inflammatory processes were highly significant among the upregulated HF-responsive genes, and glutathione metabolism and cholesterol biosynthesis were most significant among the downregulated HF-responsive genes.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 2. Functional characterization of the high-fat responsive genes.

Representative overrepresented functional categories in the set of 1663 high-fat responsive genes are grouped according to their biological function: (a) lipid and cholesterol metabolism, (b) oxidative and metabolic processes, (c) inflammatory and immune response, (d) apoptosis and protein folding, (e) cell growth and cell cycle and (f) transcription regulation and signal transduction. For each functional group, representative genes are listed and their expression profiles (average HFBT vs. chow and HFP vs. chow expression ratios per time-point) are shown in the adjacent diagrams.

https://doi.org/10.1371/journal.pone.0006646.g002

Integrative temporal and functional characterization of the hepatic transcription response to high-fat diets

To include prior biological knowledge in pathway analysis, both public and proprietary gene sets were used for Gene Set Enrichment Analysis (GSEA) [30] of each of the HFBT and HFP versus chow per time-point comparisons. In addition, to facilitate the dynamic interpretation of the identified biological functions, statistically significant GSEA results were hierarchically clustered across all conditions, thus integrating temporal and functional information into single visual output (Materials and methods). Hierarchical clustering of the GSEA-calculated normalized enrichment scores (NES) resulted in an aggregation of the similarly regulated gene sets, facilitating the visualization of pathway activities during the 16-week response to HF diets. The resulting cluster of 314 gene sets, significant (FDR q-value<0.1) in at least one of the HF conditions, can be visually divided into five temporal modules (Figure 3).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 3. Temporal modules of pathway activities during the hepatic high-fat response.

Hierarchical clustering (Pearson correlation (uncentered), average linkage) of Normalized Enrichment Scores (NES), as determined by Gene Set Enrichment Analysis. The NES values of 314 gene sets, significant (FDR q-value<0.1) in at least one of HFBT vs. chow and HFP vs. chow per time-point comparisons are used as an input for hierarchical clustering. The NES scores, represented by the color gradient, correspond to the relative up- (red) and down- (blue) regulation of the gene sets under each of experimental conditions. The cluster diagram can be divided into five main temporal themes (depicted as modules 1 to 5), highlighting the main trends in temporal pathway activities: (1) early activation/late repression; (2) constant repression; (4) constant activation and (5) early repression/mid and late activation. Module (3) includes fuzzy pathway profiles that bring to light differences in transcription response to beef tallow- (HFBT) and palm oil-based (HFP) high-fat diets.

https://doi.org/10.1371/journal.pone.0006646.g003

An examination of the cluster heatmap reveals that the vast majority of the gene sets change the direction of their transcriptional regulation throughout the time-course. The first and largest module contains gene sets that are generally significantly upregulated during the early phase (day 3) and downregulated at the late phase of the HF-response (week 8, 12 and/or 16). The most characteristic gene sets in the module 1 are associated with inflammation and immune response and their regulation, such as Interleukin (IL)-1, IL-2, IL-3, IL-4, IL-5, IL-7 and IL-9 pathways, CD40 pathway, antigen processing and presentation, T and B cell receptors signaling, natural killer (NK) cell cytotoxicity, leukocyte migration and Tumor necrosis factor (TNFalpha), Nuclear factor-kappa B (NF-κB) and Toll-like Receptor (TLR) signaling pathways. The second representative functional theme in module 1 is related to cell growth, proliferation and differentiation. Examples of these gene-sets include cyclins and cell cycle regulators, G1 to S and G2 to M checkpoint controllers, DNA replication reactome and Mitogen-activated protein (MAP) kinases, Epidermal growth factor (EGF) and Transforming growth factor beta (TGF-β) signaling pathways. Additionally, various gene sets related to cancer and development also follow the early induction/late repression expression pattern.

In contrast to module 1, gene sets in the last module (5) are generally repressed during the early phase and significantly upregulated during mid and late phases. The most prominent functional characteristic of module 5 is the presence of the PPAR signaling pathway as well as many PPARs-regulated gene sets, including those associated with adipocytes differentiation (IDX_TSA_UP_CLUSTER6, ADIP_VS_FIBRO_UP, ADIP_DIFF_UP, NADLER_OBESITY_DN), fatty acid oxidation and lipogenesis. The transcriptional activation of these gene sets during the mid phase and their amplification during the late phase implies an important role for PPARs in regulating the transition from short- to long-term effects of hepatic exposure to excess dietary fat. Additionally, hepatic activation of the genes involved in adipocytes differentiation and lipogenesis suggests that fat accumulation and adipogenic transformation likely take place in the liver after long-term exposure to a high-fat diet. Other aspects of the metabolic control, such as amino acid metabolism and tricarboxylic acid cycle, are also upregulated during mid and late phases of the HF-response.

Despite the large overall similarity in transcription response to HFBT and HFP diets, a comparison of pathway activities reveals specific differences between the two high-fat conditions. This is particularly evident in the gene sets clustered in the lowest part of the module 3. Notably, the regulation of gene sets involved in energy metabolism (glycolysis and gluconeogenesis, oxidative phosphorylation, pentose phosphate pathway, electron transport chain) is sensitive to variations in fat origin and/or specific compositions of HF diets. Specifically, a palm oil-based HF diet (HFP) causes a transient induction of the gene sets involved in energy metabolism at day 1 but shows attenuated induction of these gene sets in the late phase compared to a beef tallow-based HF diet (HFBT). Similar deviation in the gene set activity patterns between HFP and HFBT conditions is also visible in the module 5.

Finally, a small fraction of all the represented gene sets retains a constant transcriptional pattern throughout the time-course. These are represented in modules 2 (constantly downregulated gene sets, e.g. cholesterol biosynthesis, glutathione metabolism and metabolism of xenobiotics by cytochrome P450) and module 4 (constantly upregulated gene sets, such as lipoprotein metabolism).

The reciprocal transcriptional profiles of the pathways represented in modules 1 and 5 emerge as the principal signature of the transition from early to late hepatic transcription response to excess dietary fat. The coincidence of the repression of inflammatory, immune and cell proliferation pathways and the induction of metabolic, lipogenic and adipogenic pathways prompted the hypothesis that these events may be interdependent and their swap relevant for driving the transition from a stressed to pathological hepatic state. The relationships between inflammatory and metabolic processes and the key regulators controlling them were further investigated using biological network analysis.

Network analysis of the HF-responsive genes: interplay between PPARγ and NF-κB regulatory modules

To further explore the control of and the biological connectivity between the HF-responsive genes, the set of 1663 genes (Figure 1B, 1C, Table S1, Table S2) was used as an input for the network analysis within the Ingenuity Pathway Analysis suite [31]. The networks with the highest significance score (network score equal to or higher than 35) and their associated biological functions are listed in the Table 2. To focus on the interactions between the processes identified as crucial for the transition from early to late hepatic response to excess dietary fat, networks related to immune response, lipid metabolism and hepatic steatosis (networks 1, 2 and 4) were merged for further examination (Figure 4, 5). The network number was limited to three to restrict the size and facilitate the clarity of the resulting network.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 4. Molecular network underlying hepatic response to high-fat diets (day 3).

The global molecular associations of the high-fat responsive genes are functionally characterized and divided into networks based on the functions and/or diseases that are most significant to the network objects (Ingenuity Pathway Analysis). Depicted is the result of merging the network 1 (Immune Response, Tissue Development, Skeletal and Muscular System Development and Function), network 2 (Cellular Development, Connective Tissue Development and Function, Lipid Metabolism) and network 4 (Hepatic System Disease, Liver Steatosis, Cancer). The overrepresented “Function and disease” (Fx) categories “immune response” and “hepatic steatosis” are overlaid onto resulting network, showing which genes (nodes) are directly involved in these processes. The interactions between nodes that are directly connected to both processes are highlighted in pink. Color coding of the nodes corresponds to the direction of gene expression changes at day 3 in HFBT vs. chow diet comparison (upregulated genes are shown in red and downregulated in green).

https://doi.org/10.1371/journal.pone.0006646.g004

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 5. Molecular network underlying hepatic response to high-fat diets (week 12).

The equivalent network to that in Figure 4, except that the color coding of the nodes corresponds to the direction of gene expression changes at week 12 in HFBT vs. chow diet comparison.

https://doi.org/10.1371/journal.pone.0006646.g005

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Top scoring networks among the high-fat responsive genes (Ingenuity Pathway Analysis).

https://doi.org/10.1371/journal.pone.0006646.t002

The analysis identified several major regulators of the cellular response to high-fat diets appearing as the hubs in the resulting network: PPARγ, SREBF1 and SREBF2 - regulators of lipid, fatty acid and cholesterol metabolism, NF-κB – regulator of immune response and Akt – regulator of cell growth, proliferation and differentiation. Aside from the constant transcriptional shutdown of the SREBF2 local subnetwork that regulates cholesterol biosynthesis, the majority of network components show a characteristic swap in transcription response between early (day 3) and late (week 12) phase of the time-course (Figure 4, 5).

To further investigate relations between the network components that show the observed swap in transcription response and to characterize their associated functions in more detail, we examined “Function and disease” categories that are overrepresented in the network. By overlaying the most significant categories over the network, the categories “hepatic steatosis” (p-value 1.41E-08) and “immune response” (p-value 5.85E-07) were identified as the best representatives of processes that are reciprocally regulated during the early and the late phase of high-fat response. In addition to the largely reciprocal regulation, these processes are interconnected via few key network components (pink lines, Figure 4, 5). The gene expression of nearly all the network components involved in the promotion of inflammation and immune response is induced during the early phase and repressed during the late phase of the HF response, including NF-κB regulatory factors, interferon, cytokine and chemokine signaling molecules, acute phase response reactants and complement components. In contrast, network members promoting development of hepatic steatosis, such as PPARγ, SCD, SREBF1, ACOX1, CIDEC and CFD (adipsin) are repressed during early phase and induced during the late phase of hepatic response to HF diets. The exception to the observed global repression of inflammatory and immune response in the late phase are interleukin-1 pathway components and, somewhat, the acute phase reactants that show statistically insignificant, low-grade re-induction at week 12 and week 8, respectively. This indicates that a modest fraction of the inflammatory response, likely mediated by JNK/AP-1 pathway parallel to NF-κB signaling escapes the global repression [32], [33].

Particularly interesting are the network members that are associated with both inflammatory and steatotic transcriptional modules, namely PPARγ and IKBKG (inhibitor of kappaB kinase gamma). It has been previously shown that PPARγ antagonizes inflammatory responses by a transrepression of NF-κB regulators and that its hepatic activation leads to the development of liver steatosis [11], [12], [34]–[38]. In contrast, IKBKG, also known as NF-κB essential modulator (NEMO), is required for the activation of the NF-κB complex by proinflammatory stimuli, and it has recently been recognized as suppressor of hepatic steatosis, possibly through the negative interaction with PPARs [39]–[42]. The gene expression data show significant repression of PPARγ during the early phase of the high-fat response. In contrast, at the mid and/or late phases, PPARγ and its target genes are induced, while NEMO is simultaneously significantly repressed (Figure 6, Table S2). This reciprocal transcriptional activity proposes an appealing model where the direct trans-repression between PPARγ and NEMO/NF-κB regulators may occur, controlling the transition from early to late hepatic response to HF diets. Evidently, further biochemical studies are required to confirm the suggested temporal transrepression.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 6. Reciprocal activation of regulators of an inflammatory and steatotic transcriptional programs during high-fat feeding time-course.

The average gene expression profiles of NF-kB regulators (RelB, IKBKG, IKBKE, IKBKAP) (red line) and PPARγ/hepatic steatosis-associated genes (PPARγ, SREBF1, SCD1, ACOX1, CIDEC, CFD) (yellow line) during the 16-week high-fat feeding time-course. (A) HFBT vs. chow diet. (B) HFP vs. chow diet.

https://doi.org/10.1371/journal.pone.0006646.g006

Transition from hepatic inflammation to steatosis is reflected by inflammatory plasma proteins and hepatic triglyceride levels

The transition from an inflammatory (early) to steatotic (late) transcriptional program in livers of ApoE3L mice observed by gene expression profiling is supported by expression of plasma proteins and liver triglyceride content. A series of inflammatory plasma proteins were quantified by multiplex immunoassay during the complete high-fat feeding time-course (Materials and methods). The transient activation of many inflammation-related proteins has been observed during the early phase of the time-course (data not shown). The levels of plasma proteins likely reflect a systems response of multiple organs to excess dietary fat. Nevertheless, the trend of early activation coupled with late repression, which has been observed in the expression of many hepatic inflammatory genes when HF diets are compared to chow diet, can also be found in the expression of NF-κB-related plasma inflammatory proteins. Examples of such proteins are NF-κB activating protein Immunoglobulin A (Figure 7A) and NF-κB activated proteins Beta-2 Microglobulin, Interleukin-18 and Macrophage-Derived Chemokine (CCL22) (Figure 7B-D) [43]–[47]. The upregulation of Interleukin-18 at week 4 and week 8 under all three dietary conditions is followed by HF-specific repression at the two last time-points (Figure 7C).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 7. Expression changes of plasma proteins caused by HFBT and HFP high-fat diets.

Expression changes of a subset of inflammatory plasma proteins that show trend of early activation coupled with late repression compared to the control condition during the 16-week high-fat (HF) feeding time-course, as measured by multiplex immunoassay. Plotted are average protein expression levels per time point in mice fed chow, HFBT and HFP diets. Statistically significant changes in protein expression of HFBT and HFP fed mice compared to chow fed mice per time-point are marked with asterisk and hash symbol, respectively (p value<0.05). All four proteins are associated with NF-κB activation. (A) Immunoglobulin A, protein that activates NF-κB. (B–D) Beta-2 microglobulin, Interleukin-18 and Macrophage-derived chemokine (CCL22), proteins activated by NF-κB.

https://doi.org/10.1371/journal.pone.0006646.g007

The activation of the steatotic transcriptional program during the late phase of the high-fat feeding time-course observed by transcriptome analysis harmonizes with the increased hepatic triglyceride content at the late time-points (Figure 8A). To identify genes whose expression is most predictive of hepatic fat accumulation under high-fat dietary conditions, liver triglyceride content and expression of 1663 high-fat responsive genes in each animal were used to perform regression analysis by random forest modeling [48], [49]. The genes with highest importance in the resulting model (n = 30) are shown in Figure 8B. The genes we identified as steatosis-associated by the network analysis (Figure 4, Figure 5), such as ACOX1, SCD, PPARγ, CFD and CIDEC were re-discovered by the regression analysis. Also, matrix metallopeptidase 13 (MMP13) implicated in liver fibrosis and ENTPD5, associated with hepatopathy and hepatocellular tumors were identified by the regression approach 50, 51. In addition to the known markers of hepatic pathology, random forest modeling uncovered putative novel candidate genes whose expression may predict the development of hepatic steatosis. This includes fatty acid metabolism associated genes (CRAT, ACAA1b and ACAT), matrix endopeptidase related proteins (ITIH5, MME) and several genes of unknown function. The gene identified as the best predictor of hepatic steatosis by the regression analysis, butyrylcholinesterase (BCHE) has been mostly studied in the context of its effect on the brain's cholinergic system [52]. However, there is evidence showing that serum butyrylcholinesterase is associated with adiposity, serum lipid profile and insulin resistance in humans [53]. Therefore, BCHE may be the potential marker of interest for further investigation of the hepatosteatosis development.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 8. Regression analysis of gene expression and hepatic triglyceride content in chow and high-fat-fed ApoE3L mice.

(A) Changes in hepatic triglyceride (TG) content induced by HFBT and HFP high-fat diets and a control (chow) diet during the 16-week time course. Plotted are average values of hepatic triglyceride levels per time point in each of the three diets. Statistically significant increase (p≤0.01) in hepatic TG content of HFBT and HFP fed mice compared to chow fed mice (marked with asterisk and hash symbol, respectively) was observed at the time-point 16 weeks, indicating development of hepatic steatosis. (B) The most important genes for the prediction of hepatic triglyceride levels, as assessed by Random Forests Regression analysis. The expression of 1663 high-fat responsive genes and the hepatic triglyceride levels in each animal were used as an input for the regression analysis. The plot shows the importance of the top 30 genes in the Random Forest regression model. The importance of a gene in the model is expressed as the increase of Mean Squared Error (MSE) when the gene is excluded from the analysis. Higher the percentage of increase of MSE, the more important the particular gene is for the prediction of the hepatic TG levels. The genes identified as steatosis-associated by the network analysis, such as ACOX1, SCD, PPARγ, CFD and CIDEC are also discovered among the top 30 genes resulting from the regression analysis of gene expression and hepatic TG levels.

https://doi.org/10.1371/journal.pone.0006646.g008

Together, inflammatory plasma proteome and hepatic triglyceride content support findings obtained by mRNA expression analysis, further substantiating the proposed shift between the hepatic inflammatory and steatotic transcriptional program during prolonged exposure to high-fat diets.

Hepatic pathophysiological changes induced by short- and long-term exposure to excess dietary fat

In this study, we focused on investigating the molecular mechanisms underlying the onset and progression of the metabolic syndrome during a 16-week high-fat feeding time-course in ApoE3L transgenic mice. The changes in hepatic transcriptome revealed that the adaptation to excess dietary fat proceeds in three phases: early, mid and late. The early (day 1 to week 1) and the late (week 8 to week 16) phases are characterized by the most prominent, and often reciprocal peaks in gene expression changes.

During the early phase, the initial sensing of the fat overload triggers the cellular stress response, characterized by activation of acute phase reactants, inflammatory and immune response (Figure 2, 3, 4). The deregulation of the cell cycle and certain apoptotic genes during the early phase of high fat feeding resemble hepatic regeneration response, further suggesting compromised liver integrity. Additionally, the immediate perturbation occurs in mitochondrial, microsomal and peroxisomal oxido-reductive processes, including the severe repression of specific members of cytochrome P450 family and glutathione transferase genes (Figure 2, 3, Table S2). The imbalance in oxido-reductive processes accompanied with a diminished protective function of glutathione transferases may create conditions of elevated liver sensitivity to oxidative damage [54]. The hypersensitivity to oxidative stress could synergize with the lipotoxic stress, manifested by the activation of unfolded protein response (UPR), endoplasmatic reticulum (ER) stress and DNA damage response (Figure 2, 3). Apart from the stress response, the early phase of the hepatic transcriptional response to high-fat diets is also characterized by metabolic adaptation. Of all the pathways involved in lipid metabolism, early activation is exclusively observed in lipid/lipoprotein binding and transport (Figure 2, 3). This suggests that the primary, short-term metabolic adjustment may involve hepatic elimination of fatty acid surplus by means of its excretion in the form of very low-density lipoprotein (VLDL) particles. Cholesterol efflux and biliary secretion are also immediately activated, coupled with repression of cholesterol biosynthesis (Figure 2, 3).

The major signature of the late phase of the hepatic adaptation to excess dietary fat is the transcriptional induction of nearly all aspects of lipid metabolism, including those regulated by PPARα (lipolysis, fatty acid beta-oxidation), PPARγ (adipogenesis and adipogenic transformation, lipogenesis, lipid accumulation, lipid uptake) and SREBP1 (lipogenesis, fatty acid synthesis, fatty acid desaturation, fatty acid elongation) (Figure 2, 3, 5). The need for long-term lipid management is, thus, likely resolved by an activation of both lipolysis and hepatic lipid storage. The pro-steatotic transcriptional program, manifested in the activation of lipogenic, adipogenic and lipid accumulation pathways and the activation of genes such as PPARγ gene itself, stearoyl-CoA desaturase 1 (SCD1) and CIDEC (Fsp27) (Figure 2, 3, 5, 6) suggest that the adipogenic transformation of hepatocytes and the development of liver steatosis may occur during the late phase of HF feeding [34]–[37], [55]. This finding is supported by the increase of total hepatic triglycerides in the late phase of the time-course and their significantly elevated levels at the week 16 in mice fed HF diets compared to chow fed mice (Figure 8A). The incidence of hepatic steatosis under the similar experimental conditions has also been reported in the literature [21], [35]. The regression analysis of association between hepatic gene expression and the lipid accumulation identified known and novel genes that may be valuable as predictors of hepatosteatosis development.

The observed hepatic adaptation to excess dietary fat occurs similarly in investigated beef tallow- (HFBT) and palm oil- fat (HFP) based diets. Nevertheless, specific differences in the gene expression response to these two diets do exist. This is particularly evident in the expression of pathways related to energy metabolism during the early and the long-term adaptation to high-fat feeding. In general, mice fed HFP diet show more pronounced transient induction of these pathways at the very beginning of the time course, but fail to activate them as efficiently as the mice fed HFBT diet at the late phase of the time-course (Figure 3). These differences in processes that are relevant for energy expenditure may account for a higher increase in body weight (Figure 1A) and significantly higher whole-body insulin resistance (data not shown) of HFP-fed mice compared to HFBT-fed mice at the end of the time-course. This observation may be relevant in the context of dietary recommendations, suggesting that the excess of palm oil-fat based diet may be at least as harmful, if not more so, than the cholesterol-containing beef tallow fat-based diet.

The hypothesized mechanism controlling the switch from an inflammatory to steatotic hepatic state during the high-fat feeding response

The majority of pathways affected by the high-fat diets exhibit opposite regulation during the early and the late phase of the high-fat feeding time-course. This synchronous swap of the major functional signatures between the early and the late phase and the fact that the key controllers of the reciprocally regulated processes (as revealed by the network analysis (Figure 4, Figure 5)) are implicated in the mutual repression suggest that the regulatory exchange may occur via tightly controlled reactions, limited to few master regulators.

The NF-κB and Akt regulators, the key controllers of the pathways showing early activation/late repression expression mode (i.e. inflammation, immune response, cell proliferation and cell differentiation), have been previously reported to act synergistically [56]–[59]. Similarly, PPARγ and SREBP1, identified in our study as the key regulators of pathways in the early repression/late activation transcriptional module are co-acting in regulating lipid metabolism and adipogenesis [60], [61]. In contrast to the synergistic activities within these transcriptional modules, there is emerging evidence of the antagonistic activity between them, particularly regarding NF-κB and PPAR regulators. There is limited evidence for the NF-κB mediated repression of PPARγ [62]. In turn, multiple mechanisms by which PPARs inhibit inflammatory gene expression through interference with NF-κB signaling have been reported (reviewed in [11], [15], [63]). Several studies investigated the inhibitory effect of PPARγ and its agonists on the activity of the inhibitor of nuclear factor κB kinase (IKK) complex [64]–[66]. The IKBKG (NF-κB essential modulator, NEMO) is the subunit of the IKK complex indispensable for the activation of NF-κB [42]. Considering its critical role, we propose the model in which down-regulation of NEMO in the mid and the late phases of HF response could be the most important regulatory event controlling the shut-down of NF-κB driven inflammatory response at the switch point between stressed and pathological hepatic state. In addition to NEMO, other two IKK related genes (IKBKE, IKBKAP) and the gene coding for NF-κB subunit RelB also show characteristic early induction and late repression transcription mode (Figure 6, Table S2). The coincident opposite transcriptional activity of PPARγ and its target pathways suggests that the two regulatory events may be interdependent (Figure 6). This hypothesis is supported by the recent studies showing that NEMO has an essential physiological role in preventing the spontaneous development of hepatic steatosis preceding hepatocellular carcinoma and that the hepatic ablation of NEMO in mice fed HF diet increases PPARγ mRNA levels and aggravates hepatic steatosis [39], [40]. The suggested mechanism of the mutual repression of PPAR and NEMO/NF-κB regulators during their coordination of the transition from early to the late phase of HF-response remains to be confirmed by biochemical studies.

The relevance of the proposed model of the transition from high-fat induced metabolic stress to metabolic syndrome in ApoE3L mice

The dynamic functional landscape of the hepatic transcriptional adaptation to excess dietary fat during the 16-week time-course suggests a model in which sequential physiological changes underlie the transition from metabolic stress to metabolic syndrome, summarized in Figure 9. These results provide novel insight into the delicate cross-talk between inflammation and lipid metabolism in controlling the progression of metabolic disease development. It is important to note that the late repression of the hepatic NF-κB driven inflammatory/immune response may not seem in line with the established model of the obesity-associated inflammation, well studied in the adipose tissue [10], [67]. The previously shown association of inflammatory signaling pathways with obesity and hepatic steatosis is the most prominent feature observed in our data. Nevertheless, the assumption that, similar to the situation in adipose tissue, hepatic inflammation is secondary to hepatic steatosis is not supported by our findings [68], [69]. In fact, our results show that the temporal order of events contradicts this assumption, at least during the examined time frame and in our mouse model. Extending the duration of treatment and/or increasing amount of excess dietary fat would likely provoke transitions to further grades of severity in hepatic pathology such as hepatosteatitis, fibrosis, cirrhosis and hepatocellular carcinoma. To fully understand the complex relationship between inflammation and metabolic syndrome, information originating from different organs and at various time points needs to be considered on the systems level [69].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 9. Model of the hepatic physiological response to high-fat diets during the 16-week time-course.

The summary of a proposed model for the hepatic physiological changes in response to high-fat diets during the 16-week time-course in ApoE3L mice. The initial perturbation of hepatic homeostasis by excess dietary fat triggers the stress response largely controlled by NF-κB and Akt regulators and manifested in activation of acute phase response, inflammation, immune response, hepatic regeneration-like response and lipotoxicity (day 1 to week 1). Upon prolonged high-fat feeding, liver fails to regain the basal state and consequently shifts to pathological state controlled by PPAR and SREBP regulators and characterized by hepatic lipid accumulation and adipogenic transformation, indicative of hepatic steatosis (late phase, week 8 to week 16). The flagship processes induced at the early and at the late phase are shown in boxes. The transition between the stressed and the pathological hepatic state may be controlled by trans-inhibitory interactions between NF-κB and PPARγ regulators, resulting in the tradeoff between inflammatory and steatotic transcription programs (mid phase, week 2 to week 4). On the systems level, the activation of steatotic program is followed by other metabolic syndrome associated pathologies such as obesity and whole-body insulin resistance.

https://doi.org/10.1371/journal.pone.0006646.g009

Finally, the identified central role of ppar and nemo/nf-κb regulators in coordinating the onset and progression of metabolic syndrome may have important implications in treatment of the disease. currently, pparγ ligands are used in clinics for their anti-inflammatory and insulin-sensitizing effects in diseases such as psoriasis, atherosclerosis, inflammatory bowel disease and type 2 diabetes [13], [14], [70]–[74]. our results suggest that pparγ activation could also have negative effects in disorders where the cross-talk between inflammation and lipid metabolism is essential to the development of pathogenesis, as it is in diet-induced metabolic syndrome. the activation of pparγ may inhibit nf-κb and therefore suppress inflammation, but in turn evoke transition to pathological state, in this case hepatic steatosis. despite the apparent harmful effects of inflammation, such as triggering insulin resistance, its protective physiological role in preventing transitions to even less preferable system states should not be neglected. thus, the tradeoff between the beneficial and harmful effects of altered pparγ activity should be carefully considered when using pparγ ligands and nf-κb-inhibiting agents to ameliorate metabolic syndrome associated pathologies. the presented findings demonstrate the use of high-throughput dataset analyses as a starting point for generating testable hypotheses that may open new avenues for dietary prevention strategies, clinical research and pharmaceutical therapies.

Ethics Statement

Animal experiments were approved by the Institutional Animal Care and Use Committee of the Netherlands Organization for Applied Scientific Research (TNO) and were in compliance with European Community specifications regarding the use of laboratory animals.

Animals and diets

The study involved 186 male ApolipoproteinE3-Leiden transgenic mice [24] at 14±2 weeks of age. Apolipoprotein E3-Leiden (ApoE3L) transgenic mice display lipoprotein profile that closely resembles that of humans and develop human-like dysbetalipoproteinemia and atherosclerotic lesions when fed Western-type diets. The age of 14 weeks was chosen as optimal because the animals are considered adult at that stage (neither juvenile at early time points nor aged at the late time points of the study). The first (control) group of ApoE3L mice was fed standard chow diet (RM3 (E) DU; Special Diet Services, Witham, Essex, UK), the second group of mice was fed a high-fat diet based on animal fats (HFBT, beef tallow fat high-fat diet +0.25% cholesterol, also referred to as Western-type diet, Arie Blok BV, The Netherlands) and the third group of mice was fed a high-fat diet based on plant fats (HFP, palm oil high-fat diet based on open source D12451 SDS diet in which lard is replaced with palm oil, Research Diet Services, Woerden, The Netherlands). The macronutrient content and the fatty acid composition of chow and high-fat diets are provided in Supporting .

Experimental design and sample preparation

From three weeks prior to diet intervention onwards, all animals were fed a standard chow diet. At the beginning of the study, mice were divided into three groups: (1) control group fed chow diet, (2) group fed HFBT diet and (3) group fed HFP diet. Because the interest of the study was to asses effects of high-fat diets under physiological conditions, animals were fed ad libitum. Series of control experiments employing a metabolic cage setup showed that C57Bl/6 mice, the genetic background of the APOE3L mice, have isocaloric food intake when fed low fat and HFBT and HFP diets. The light cycles were identical for all animals. For mRNA expression profiling, six mice from each diet group (total n = 150) were sacrificed at time points 0 days (chow only), 1 day, 3 days, 1, 2, 4, 8, 12 and 16 weeks, their livers were dissected after 4 hour fasting period (typically from 9 h AM to 1 h PM), snap frozen in liquid nitrogen and stored at −80°C until further processing.

RNA isolation, labeling and hybridization to microarrays

Total RNA was isolated using TRIzol reagent (Invitrogen, Breda, The Netherlands) according to the manufacturer's instructions. RNA was treated with DNAse and purified using the SV total RNA isolation system (Promega, Leiden, The Netherlands). Concentrations and purity of RNA samples were determined on a NanoDrop ND-1000 spectrophotometer (Isogen, Maarssen, The Netherlands). RNA integrity was checked on an Agilent 2100 bioanalyzer (Agilent Technologies, Amsterdam, The Netherlands) with 6000 Nano Chips according to the manufacturer's instructions. RNA was considered as suitable for array hybridization only if samples exhibited intact bands corresponding to the 18S and 28S ribosomal RNA subunits, displayed no chromosomal peaks or RNA degradation products, and had a RNA integrity number (RIN) >8. Applying this criterion, 142 RNA samples were used for hybridization to microarrays, including 5 to 6 biological replicate samples per diet, per time-point. RNA samples were hybridized to NuGO Affymetrix Mouse GeneChip arrays (NuGO_Mm1a520177) containing 23865 probesets including 73 control probesets [26]. Arrays were scanned on a GeneChip Scanner 3000 7G (Affymetrix). Detailed methods for labeling, hybridizations to the arrays and scanning are described in the eukaryotic section of the GeneChip Expression Analysis Technical Manual, Revision 3, from Affymetrix, and are available upon request. The gene expression data are made available via ArrayExpress repository [75].

Data preprocessing, differential expression analysis and Gene Set Enrichment Analysis

Quality control of microarray data, normalization, differential expression analysis and Gene Set Enrichment Analysis [30] were performed using packages from the R/Bioconductor project [76], [77] through the Management and Analysis Database for MicroArray eXperiments (MADMAX) analysis pipeline [78].

Quality control of the hybridized microarrays was performed using simpleaffy and affyplm packages. Upon rigorous examination of the resulting diagnostic plots, 116 microarrays of the supreme quality were taken for the further analysis. This resulted in analysis of 3 to 6 biological replicate samples per diet, per time-point. Gene expression estimates were calculated using the library GC-RMA, employing the empirical Bayes approach for background correction followed by quantile normalization. The custom MBNI CDF-file (MmNuGOMm1a520177 version 9.0.1), available at http://brainarray.mbni.med.umich.edu/Brainarray/Database/CustomCDF/CDF_download_v9.asp and http://nugo-r.bioinformatics.nl/NuGO_R.html [79], [80] was used to re-annotate the probes to new probesets, remove poor quality probes and derive unique signal values for different probesets representing the same gene. This resulted in gene expression values for 15105 genes with unique identifiers.

Differentially expressed genes between control and each of treatment groups per time point, as well as between each of time points and day 0, were identified using the limma package, applying linear models and moderated t-statistics that implement empirical Bayes regularization of standard errors [27]. False discovery rate of 10% (q-value<0.1) was used as a threshold for significance of differential expression. Significance of the overlap between differentially expressed genes in HFBT and HFP conditions (Figure 1B) was calculated by hypergeometric distribution.

The t-test values of differential expression between control and each of treatment groups per time point calculated using limma package were used as the input for the PreRanked scoring method within the Gene Set Enrichment Analysis (GSEA). Gene sets collection included 880 gene sets compiled from MSigDB C2, Biocarta, Kyoto Encyclopedia of Genes and Genomes (KEGG) and GenMAPP databases as well as the expert curated gene sets. Detailed information about gene sets used for GSEA analysis, including source websites is available upon request. Gene set size filter (min = 15, max = 500) resulted in filtering out 405 of 880 gene sets. The number of permutations for was set to 1000. Gene sets are considered significantly enriched at false discovery rate (FDR) smaller than 10% (q-value<0.1). In total, 314 gene sets are identified as significantly enriched in at least 1 of 16 comparisons (HFBT vs. chow and HFP vs. chow per time-point). Normalized enrichment scores (NES) of significantly enriched pathways and the corresponding FDR q-values across all experimental conditions are available upon request.

Hierarchical clustering and data visualization

Hierarchical clustering and visualization of gene expression changes in Figures 1C, 2 and 6 were performed in GeneSpring GX 7.3.1 software (Agilent), using average gene expression values under HFBT and HFP conditions divided by the median of chow gene expression, per time-point. Hierarchical clustering and visualization of NES values in Figure 3 were performed using HierarchicalClusteringViewer and HeatmapViewer modules within the GenePattern analysis suite [81]–[83].

Analysis of gene expression profiles of HF-responsive genes

Temporal analysis of the gene expression profiles where each of HF groups is compared to chow per time point (Figure S2, S3, S4 and S5) was performed using Smoothing Spline Clustering algorithm [28] through the R computing environment [76], applying the settings nchain = 5 and nclust = 24. Setting the number of clusters to 24 yielded Bayesian information criterion (BIC) values that were at the bottom of the U-shaped BIC curve (i.e. closely approaching the optimum), while resulting in images of reasonable complexity. As input values, average log₂ ratios of 1663 high-fat responsive genes were used (HFBT and HFP divided by median of chow per time point). To facilitate visualization of temporal trends, the starting time point (t = 0, only chow fed mice) was also included in the analysis.

Overrepresentation analysis of functional categories

Identification of overrepresented functional categories among 1663 high-fat responsive genes and their grouping into functionally related clusters was performed using DAVID Functional Annotation Clustering tool [29]. The analysis was performed using Gene Ontology, Protein domains, Pathways and Functional categories according to the default settings (Table S3, version date July 2007). Representative functional categories from the most statistically significant clusters are manually selected and listed in Figure 2. Equivalent analysis is performed using separately upregulated (922 genes) and downregulated (786 genes) high-fat responsive genes (Table S3, version date June 2009).

Network analysis

The network analysis was generated through the use of Ingenuity Pathways Analysis (version date March 2008) [31]. The data set containing gene identifiers and corresponding expression values for 1663 high-fat responsive genes was uploaded into the application. Of 1663 gene identifiers, 1660 were successfully mapped to its corresponding gene objects in the Ingenuity Pathways Knowledge Base and 1303 were identified as network eligible. These genes were overlaid onto a global molecular network developed from information contained in the Ingenuity Pathways Knowledge Base. Networks were then algorithmically generated and graphically represented based on connectivity of genes. The Functional Analysis of the network identified the biological functions and/or diseases that were most significant to the genes in the network (version date June 2009). Fischer's exact test was used to calculate a p-value determining the probability that each biological function and/or disease assigned to that network is due to chance alone. Genes or gene products are represented as nodes, and the biological relationship between two nodes is represented as an edge (line). The node color indicates up- (red) or down- (green) regulation. The direction of average expression changes in HFBT vs. chow comparisons (Table S2) at day 3 and week 12 was used for color coding in Figure 4 and Figure 5, respectively.

Accession numbers

For the microarray experiments described in this study, MIAME compliant protocols and datasets in Tab2MAGE are accessible from ArrayExpress microarray data repository [75] with the accession number E-TABM-553.

Plasma proteins and hepatic triglyceride content analysis

Plasma proteins were quantified by multiplex immunoassay measurements at Rules Based Medicine (Rules Based Medicine, Inc., Austin, Texas, USA). Plasma antigens immunoassay panel included in the Rodent Multi-Analyte Profile was used for measurement of expression levels of 58 proteins (RodentMAP version 2 antigen panel). Of these, 47 proteins had sufficient detectability of the expression signals and were included in the further analysis. Statistical significance of protein expression in HFBT and HFP fed mice compared to chow fed mice per time-point (marked with asterisk and hash symbol, respectively) was assessed by t-test. The p value of 0.05 was used as a threshold for significance.

Liver lipid content was defined as total triglyceride content (mmol) per mg of protein. Extraction was performed using a modified Bligh and Dyer extraction protocol, optimized for steatotic liver material. Triglyceride content was measured enzymatically using the Roche TG kit (Roche, cat. No. 11488872216). Protein content was measured by BCA analysis (BCA protein Assay Kit, Pierce, cat. No. 23225). Statistical significance of hepatic triglyceride content in HFBT and HFP fed mice compared to chow fed mice per time-point (marked with asterisk and hash symbol, respectively) was assessed by t-test. The p value of 0.01 was used as a threshold for significance.

Regression analysis of gene expression and hepatic triglyceride content

To determine the relation between gene expression data and triglyceride levels in the liver tissue, Random Forests regression was used through randomForest package of R statistical computing environment [48], [49], [76], [84]. The gene expression values of 1663 high-fat responsive genes and the hepatic triglyceride levels of the matching animals at each time-point were used as an input for the analysis. In Random Forests regression, a group of regression trees is used for the model. In every individual regression tree, the data are split at each node until all objects are in a single leaf. For each tree only a subset of the data is used and at each node only a subset of the variables is used to make the split. By intentionally incorporating only part of the data for each tree, the data which are not used to generate the tree can be used to assess its performance. The term used for this measure is the Out Of Box (OOB) error. Besides the OOB measure for the overall model, it is possible to asses the importance of individual variables in the model. Here we used the increase of Mean Squared Error (MSE). The Mean Squared Error is a general measure for the overall fit of a model and its increase measures the importance of individual variables (in this case genes) when they are excluded from the analyses. The higher the increase of the MSE when a gene is excluded from the analysis, the more important that particular gene is for the fit of the model (and for the prediction of triglyceride levels in liver). The MSE of the Random Forests regression model presented here is 49.9 percent. The top 30 genes that have the highest increase of MSE and therefore provide best prediction of triglyceride levels are shown in Figure 8B. The term prediction is used here to refer to the ability to predict the triglyceride levels as they were measured on the basis of gene expression, rather than in the sense of prognostic for the future progress of steatosis in time.

Random Forests is an improved Classification and Regression Trees (CART) method. It grows many classification trees or regression trees, hence the name ‘Forests’. For quantitative outcomes the forest is made of regression trees, where the tree predictor is the mean value of the training set observations in each terminal leaf. In our application of random forests the outcomes are quantitative, therefore the regression algorithm and not classification was used. Every tree is built using a deterministic algorithm and the trees are different owing to two factors. First, at each node, a best split is chosen from a random subset of the predictors rather than all of them. Second, every tree is built using a bootstrap sample of the observations. The out-of-bag (OOB) data, one-third of the observations, are then used to estimate the prediction accuracy. Unlike other tree algorithms, no pruning, trimming of the fully grown tree, is involved. Each observation is assigned to a leaf, the terminal node of a tree, according to the order and values of the predictor variables. For a particular tree, the predictions for observations are given only for the OOB data. The Random Forest predictor is computed by averaging the tree predictors over trees for which the given observation was OOB. Because the prediction for an observation is based on trees grown without the observation, an idea akin to cross-validation, the estimated errors are unbiased and the data were not divided in test and training sets.