Time-dependent network analysis reveals molecular targets underlying the development of diet-induced obesity and non-alcoholic steatohepatitis

Abstract

Prolonged high-fat diet leads to the development of obesity and multiple comorbidities including non-alcoholic steatohepatitis (NASH), but the underlying molecular basis is not fully understood. We combine molecular networks and time course gene expression profiles to reveal the dynamic changes in molecular networks underlying diet-induced obesity and NASH. We also identify hub genes associated with the development of NASH. Core diet-induced obesity networks were constructed using Ingenuity pathway analysis (IPA) based on 332 high-fat diet responsive genes identified in liver by time course microarray analysis (8 time points over 24 weeks) of high-fat diet-fed mice compared to normal diet-fed mice. IPA identified five core diet-induced obesity networks with time-dependent gene expression changes in liver. These networks were associated with cell-to-cell signaling and interaction (Network 1), lipid metabolism (Network 2), hepatic system disease (Network 3 and 5), and inflammatory response (Network 4). When we merged these core diet-induced obesity networks, Tlr2, Cd14, and Ccnd1 emerged as hub genes associated with both liver steatosis and inflammation and were altered in a time-dependent manner. Further, protein–protein interaction network analysis revealed Tlr2, Cd14, and Ccnd1 were interrelated through the ErbB/insulin signaling pathway. Dynamic changes occur in molecular networks underlying diet-induced obesity. Tlr2, Cd14, and Ccnd1 appear to be hub genes integrating molecular interactions associated with the development of NASH. Therapeutics targeting hub genes and core diet-induced obesity networks may help ameliorate diet-induced obesity and NASH.

Introduction

The global prevalence of overweight and obesity has increased every decade, with an estimated 1.46 billion adults overweight, 205 million men obese, and 297 million women obese (Finucane et al. 2011). Obesity is associated with increased risk of multiple chronic conditions included diabetes, cardiovascular diseases, stroke, fatty liver, non-alcoholic steatohepatitis (NASH), and other comorbidities (Guh et al. 2009; Strazzullo et al. 2010). Recent reports indicate current obesity trends will continue to cause substantial morbidity, mortality, and economic burden over the coming decades (Wang et al. 2011). Therefore, there is considerable interest in identifying molecular targets underlying the development of obesity and related comorbidities, which will facilitate the development of novel therapeutics.

Disease phenotypes are the result of complex molecular interactions between genes and proteins. The development of systems biology approaches has allowed candidate disease genes identified through genome-wide microarray studies to be considered within larger integrated molecular networks based on gene–gene interactions and protein–protein interactions (Ideker and Sharan 2008). Many studies have reported differences in lipid metabolism, inflammation, and insulin resistance-related genes based on cross-sectional microarray analysis of various murine diet-induced obesity models (Gregoire et al. 2002; Kim et al. 2004; Matsuzawa et al. 2007; Radonjic et al. 2009; Fraulob et al. 2010; Kleemann et al. 2010). Animal models of diet-induced obesity are characterized by hepatic steatosis, an early event in non-alcoholic fatty liver disease (NAFLD) and non-alcoholic steatohepatitis (NASH) development (Hebbard and George 2011). However, network studies of obesity and related comorbidities such as NASH have predominately been limited to static network analysis. In morbidly obese subjects with NAFLD undergoing bariatric surgery network analysis of genes differentially expressed compared to controls identified four primary gene networks associated with NAFLD including cell death, immunological disease, cellular movement, and lipid metabolism (Gawrieh et al. 2010). A more recent study reported gene expression differences in morbidly obese patients even with histologically normal liver and no overt sign of steatosis compared to controls (Bertola et al. 2010), indicating early transcriptional biomarkers may precede development of hepatic steatosis and non-alcoholic steatohepatitis (NASH). Unfortunately, these studies do not consider the dynamic changes in gene expression, which we know occur during the pathogenesis of diet-induced obesity and the development of obesity-related comorbidities (Radonjic et al. 2009; Do et al. 2011).

To date, only two microarray studies have considered the time course of global transcriptional changes in liver during the development of obesity using animal models (Radonjic et al. 2009; Do et al. 2011). In ApoE3Leden, mice fed a high-fat diet and sacrificed at multiple time points over 16 weeks, the transition from an inflammatory to a steatosis transcriptional program emerges as a principal signature of the hepatic adaptation to excess dietary fat; however, these changes may also be attributable to aging (Radonjic et al. 2009). In C57BL/6 J, mice fed a high-fat diet over 24 weeks, time course microarrays and metabolic analysis reveal evidence of compensatory transcriptional and metabolic changes in liver to resist the increase in energy balance, but adaptation appears to be disrupted later during the development of diet-induced obesity (Do et al. 2011). Interestingly, evidence from morbidly obese patients undergoing bariatric surgery indicates while some differentially expressed genes may plausibly cause NASH, some hepatic genes appear to play a hepatoprotective role and are part of a positive adaptation which may slow NASH development (Baranova et al. 2007). Taken together, these studies highlight the necessity of time-dependent diet-induced obesity models to improve our understanding of the dynamic changes occurring in molecular networks underlying the development of obesity and NASH.

Time-dependent network approaches have been applied previously to elucidate the molecular interactions underlying the systemic inflammatory response in human leukocytes (Calvano et al. 2005). However, to date, the long-term dynamic molecular network underlying the development of diet-induced obesity in liver or diet-induced non-alcoholic steatohepatitis has not been previously reported. Accordingly, the aim of this study was to determine the core diet-induced obesity networks altered in a time-dependent manner in high-fat diet-fed C57BL/6 J mice compared to age-matched normal diet-fed controls. Furthermore, we identified hub genes in diet-induced obesity and an underlying protein–protein interaction network, which may cause the development of NASH.

Materials and methods

Experimental design, animals, and diets

A time course experimental design was used, based on a well-characterized long-term diet-induced obesity model, which included normal diet and age-dependent controls as previously reported (Do et al. 2011). Three hundred and sixty male inbred C57BL/6 J mice were obtained from the Jackson Laboratory (Bar Harbor, ME) at 4 weeks of age. All the mice were individually housed under a constant temperature (24 °C) and 12-h light/dark cycle, fed the AIN-76 semi-purified diet for 1 week acclimation period after arrival, and then randomly divided into an ND and an HFD group, with 18 mice per group. After the acclimation period, the HFD group was fed the AIN-76 semi-purified diet consisting of 20 % fat and 1 % cholesterol and ND group was fed the AIN-76 diet. (American Institute of Nutrition 1977). Mice were anaesthetized and killed at 0, 2, 4, 6, 8, 12, 16, 20, and 24 weeks. The liver tissue was then removed, rinsed, weighed, immediately frozen in liquid nitrogen, and stored at −70 °C until further analysis. The study protocol was approved by the Ethics Committee for animal studies at Kyungpook National University, Republic of Korea.

Histopathological analysis of hepatic fibrosis

Hepatic tissue was removed from each mouse and cleaned free of connective tissues. Samples were then fixed in 10 % (v/v) paraformaldehyde/PBS, embedded in paraffin, and sliced into 3–5-μm thick sections. Sections were stained by Masson’s trichrome staining to visualize the connective tissue with collagen fibers. The stained area was viewed using a microscope at a magnification of 200×.

Microarray analysis

Total RNA was amplified and purified using the Ambion Illumina RNA amplification kit (Ambion, USA) to yield biotinylated cRNA. A total of 750 ng biotinylated cRNA per sample was hybridized to Illumina Mouse WG-6 v2 Expression Beadchips (Illumina, USA) for 16 h at 58 °C according to the manufacturer’s instructions. Hybridized arrays were washed and stained with Amersham fluorolink streptavidinCy3 (GE Healthcare Bio-Science, UK) following the standard protocol in the bead array manual. The quality of hybridization and overall chip performance were determined by visual inspection of both internal quality controls and the raw scanned data in the Illumina BeadStudio software. Probe signal intensities significantly higher than background intensities were determined (detection p value <0.05), but probe signal data were not filtered to preserve probes with low expression at different time points. Probe signal intensities were quantile normalized and log-transformed. Microarray analysis was performed in the ArrayAssist® software (Stratagene, USA), Bioconductor and R programing language. LIMMA was used to determine significant differentially expressed genes (HFD responsive genes) between HFD- and ND-fed mice at each time point based on FDR < 5 %, Benjamin and Hochberg-adjusted p value <0.05, and log-fold change >1 (Smyth 2005).

Core diet-induced obesity network analysis

Three hundred and thirty-two high-fat diet responsive genes were used for core network analysis identified previously (Do et al. 2011). General study flow for the network analysis is illustrated in Fig. 1. To examine the interactions between the high-fat responsive genes, we used Ingenuity pathway analysis (Ingenuity® Systems, www.ingenuity.com, IPA version: 8.8 (2010), Content version: 3204 (2010)). IPA allows the identification of network interactions and pathway interactions between genes based on an extensive manually curated database of published gene interactions. We uploaded the ProbeID and the associated expression value from the microarray data into IPA. Each probeID was mapped to its corresponding gene object in the Ingenuity pathway knowledge base. These genes, called focus genes, were overlaid onto a global molecular network based on Ingenuity knowledge database. Networks of these focus genes were then algorithmically generated based on their types of interactions (direct and/or indirect). Scores were generated (based on Fisher’s test) to rank networks according to how relevant they are to the genes in the input dataset. The score takes into account the number of focus genes from our lists in the network, and the size of the network to approximate how relevant this network is to the original list of focus genes. The network is then presented as a graph indicating the molecular relationships between genes/gene products. Genes or gene products are represented as nodes, and the biological relationships between nodes are represented as an edge (line). The connectivity of genes (nodes) is based on the data in the IPA knowledge base, which is a large repository of gene–phenotype associations, molecular interactions, chemical knowledge, and regulatory events, manually curated from scientific publications. The node color indicates the degree of up-regulation (red) or down-regulation (green). Nodes are displayed using various shapes that represent the functional class of gene product. Edges are displayed with various labels that describe the nature of the relationship between the nodes. Time course microarray data were mapped onto the core networks to explore the dynamic changes during the development of diet-induced obesity.

Fig. 1

A study flow chart showing hepatic network analysis based on time course microarrays (2, 4, 6, 8, 12, 16, 20, 24 weeks) of liver from long-term high-fat diet-fed C57BL/6 J mice

The biological functions and/or diseases that were significantly associated with the genes in the core networks were identified by functional analysis based on the Ingenuity’s knowledge base. The right-tailed Fisher’s exact test was used to calculate a p value determining the probability that each biological function and/or disease assigned to a specific core diet-induced obesity network was due to chance alone.

Hepatic steatosis protein–protein interaction network analysis

Core networks were merged and overlayed with related “functions and diseases” to determine genes associated with two specific hepatic responses (immune/inflammatory and liver steatosis) according to the IPA knowledge base. Those genes associated with both liver steatosis and inflammation were used to build a protein–protein interaction (PPI) network using the STRING database which contains interactions between 2.5 million proteins from 630 organisms and is widely used to investigate protein regulatory networks in disease (Jensen et al. 2009). To enrich this PPI network, additional genes that are potentially associated with these overlapped genes were included through a broad mining of the String database and differentially expressed genes identified in our time-dependent diet-induced obesity model. The molecular interaction repository MiMI, which integrates data from multiple sources, was used to map protein to functional genes (Gao et al. 2009). Cytoscape was used to visual protein–protein interaction networks, which shows nodes and edges as proteins and interactions between two proteins, respectively (Shannon et al. 2003). A human protein–protein interaction network was also built with human orthologs of the mouse genes associated with both liver steatosis and inflammation in our diet-induced obesity model. Orthology between human and mouse was determined by searching ortholog genes via NCBI.

Results

Time-dependent changes in hepatic fibrosis and metabolic biomarkers in diet-induced obesity

To determine the time course of fibrosis during the development of diet-induced obesity, we performed Masson’s trichrome staining of livers from the HFD group compared with the ND group at 0, 2, 4, 6, 8, 12, 16, 20, and 24 weeks. The HFD-fed mice displayed signs of fibrosis from week 16 and accumulation of collagen surrounding the central venule (Fig. 2). However, fibrosis and collagen accumulation were more visible at week 24 in the HFD-fed mice. We also observed the C57BL/6 J mice used in this study developed signs of fibrosis as a result of aging even when fed a normal diet, although fibrosis was not morphologically detectable until week 16 in the ND group and was much less prominent than the HFD group at week 16 (Fig. 2). Time-dependent changes in blood glucose, plasma lipids, hepatic lipids, and body fat mass in this diet-induced obesity model have been previously reported (Do et al. 2011) and are shown in Table 1.

Fig. 2

Fibrotic changes of liver sections in time-dependent manner (Masson trichrome stain, original magnification ×200). a Normal diet groups, b High-fat diet groups. Red arrow fibrosis, black arrow lipid droplet (color figure online)

Table 1 Phenotype biomarkers of high-fat-fed mice

Core diet-induced obesity networks

Three hundred and thirty-two high-fat responsive genes were used for the core network analysis in IPA to identify significant networks and functions underlying diet-induced obesity. The top 5 networks exhibiting the highest scores and their associated network functions are described in Table 2 and Fig. 3. The top biological and toxicological functions associated with these core diet-induced obesity networks are shown in Table 3. Dynamic gene expression changes in the top 5 networks are shown in video format (Supplementary material). We also show the gene expression changes in the core diet-induced obesity networks at week 2 and week 24 in Fig. 3. Some specific network genes showed marked changes in a time-dependent manner based on time course microarrays as shown in Fig. 4.

Table 2 Top 5 core diet-induced obesity networks generated with IPA, based on high-fat diet responsive genes identified in long-term high-fat diet-fed C57BL/6 J mice

Fig. 3

Core networks and time-dependent changes in gene expression during the development of diet-induced obesity. Ingenuity pathway analysis was used to find core networks based on 332 high-fat diet responsive genes. Gene expression changes (log₂ fold change) are shown on as red (+) or green (−). Blue lines indicate genes, which were altered in a time-dependent manner from 2 weeks to 24 weeks in high-fat diet-fed mice. a Network 1; Cell-to-cell signaling and interaction, hematological system development and function, cellular movement. b Network 2; Lipid metabolism, small molecule biochemistry, endocrine system disorders. c Network 3; Hepatic system disease, lipid metabolism, molecular transport. d Network 4; Inflammatory response, cellular movement, immunological disease. e Network 5; Hepatic system disease, lipid metabolism, molecular transpor (color figure online)

Table 3 Top bio functions associated with high-fat diet responsive genes identified in long-term high-fat diet-fed C57BL/6 J mice

Fig. 4

Time-dependent gene expression changes in network 1–5 over 24 weeks in high-fat diet-fed mice compared to normal diet-fed mice. Based on probe intensity from Illumina Mouse WG-6 v2 beadchips

Network 1 was the most significant network composed of genes related to immune responses, cell cycle, and lipid homeostasis-related processes. Pro-inflammatory action of various chemokines accompanying toll-like receptor signaling pathway was central in Network 1. These inflammatory genes appear to interact with Cd14 and critical cell cycle regulators (Ccnd1 and Cdkn1a) as well as lipid biosynthetic and transport genes (Acacb and Apom) (Fig. 3a).

Genes in Network 2 (Fig. 3b) were related to lipid metabolism, the PPAR signaling pathway, and the adipokine signaling pathway, mediated by leptin (Lep). Figure 3b shows the expression of major genes in Network 2 where cell death-inducing DFFA-like effector a (Cidea), acyl-CoA oxidase 1 (Acox1), S100 calcium binding protein A11 (S100a11), lymphocyte antigen 6 complex D (Ly6d), lipoprotein lipase (Lpl), enoyl-Coenzyme A hydratase/3-hydroxyacyl Coenzyme A dehydrogenase (Ehhadh), and isocitrate dehydrogenase 1 (Idh1) were significantly up-regulated in a time-dependent manner in high-fat diet-fed mice compared to normal diet-fed controls.

Network 3 (Fig. 3c) was associated with hepatic system disease, lipid metabolism, and molecular transport. Network 3 was composed of hepatic-specific genes responsible for homeostasis of cholesterol and bile acids such as cytochrome P450 families and nuclear receptor subfamilies. Network 3 was also composed of genes related to “nuclear receptors in lipid metabolism and toxicity” such as Abcc2, Abcc3, Cyp7a1, Cyp8a1, Nr1l2, and Nr1l3. Among these genes, Abcb11, Abcc3, and Cyp7a1 were up-regulated, while Cyp8a1 and Nr1l3 were down-regulated when compared between week 2 and week 24. Tlr3 and Tlr7 were also slightly up-regulated genes among several toll-like receptors in Network 3.

Network 4 (Fig. 3d) was associated with the inflammatory response, cellular movement, and immunological disease and was composed of pro-inflammatory factors including tumor necrosis factor (Tnf), interferon gamma (Ifng), interleukins (Ils), and serum amyloid A2 (Saa2). In Network 4 (Fig. 3e), cytokines Il2, Il18, and Ccl2 were down-regulated, while liver fibrosis-associated genes, Csf1, Lgmm, and Gpnmb were up-regulated (Fig. 3).

Network 5, which was characterized by biological functions related to lipid and cholesterol metabolism, was composed of Ppara, Pparg, Acadm, Cyp7a1, Me1, Fdps, Fads1, Fads2, and Ptgs2. Network 5 included down-regulated genes related to the biosynthesis of cholesterol and fatty acids (Fdps and Me1) as well as genes playing a critical role in insulin and adipokine signaling pathways (Ppargc1a and Socs2). Furthermore, expression of Cyp7A1 responsible for bile acid synthesis, as well as Ptgs2, Acadm, and Gabpa responsible for inflammatory prostaglandin were down-regulated.

Hepatic steatosis hub genes connected via a protein–protein interaction (PPI) network

Next, we sought to identify hub genes that may drive the development of obesity-related comorbidities focusing specifically on hepatic steatosis, an early event in the development of NASH. When we overlaid “diseases and functions” on the merged core networks in IPA, forty-one genes were associated with the “inflammatory response” (p value 9.98E−30) and fourteen genes were associated with “hepatic steatosis” (p value 2.90E−12). We found nine genes associated with both the “inflammatory response” and “hepatic steatosis,” these were Ccnd1, Tnf, Tlr2, Il10, Lep, Cd14, Fabp4, Pparg, and Ppara. Ccnd1, Tlr2, and Cd14 were selected as potential hepatic steatosis hub genes as these genes showed distinct time-dependent expression in the high-fat diet-fed mice. To explore further potential molecular interactions between these hub genes, we first built three separate PPI networks (one depth degree) using the Mus musculus protein–protein interaction (PPI) database (Jensen et al. 2009) to generate further seed genes (Fig. 5a). Of the seed genes identified using this approach, eighteen seed genes were differentially expressed in the high-fat diet-fed mice, including Ccnd1, Cd14, Cd44, Csf1r, Cyp8b1, Egfr, H2-Aa, Junb, Lgals3, Lpl, Nsdhl, S100a8, and Tlr2 (Fig. 5a), of which thirteen genes were mapped to proteins via MiMI in Cytoscape. We then built an expanded hepatic steatosis PPI network based on these thirteen genes using the Cytoscape plug-in MiMI (Fig. 5b), with genes without a linkage to the central PPI network removed. Accordingly, as shown in Fig. 5c, Tlr2, Cd14, and Ccnd1 were connected via a PPI network associated with the “ErbB and insulin signal pathway (yellow dotted line).” In addition, these hub genes interacted with proteins forming a cluster associated with “lipid transport system” and “DNA damage responses or repair system (green dotted line).”

Fig. 5

Hepatic steatosis protein–protein interaction (PPI) network for (I) mice and (II) humans. IPA revealed Tlr2, Cd14, and Ccnd1 were high-fat diet responsive hub genes associated with both the inflammatory response and hepatic steatosis. Using the PPI database (I), eighteen high-fat responsive genes were selected from PPI network built with Tlr2, Cd14, and Ccnd1 using String database (a). Final PPI network was constructed with 13 matched genes among 18 genes in Cytoscape plug-in for MiMI and rearranged for easy interpretation (b). Significantly differentially expressed genes linking the inflammatory response and hepatic steatosis were Tlr, Cd14, and Ccnd1 (red nodes). The protein–protein interaction network revealed two gene clusters associated with “lipid transport and DNA damage repair systems” (green nodes) and “ErbB/insulin signaling pathways” (yellow nodes). Lpl, Egfr, Lgals3, and Csf1r may play a role in the progression of hepatic insulin resistance, macrophage activation to severe non-alcoholic fatty liver disease in diet-induced obesity (c). (II) A human PPI network was built with 9 genes linking between “inflammation response” and “hepatic steatosis” using Cytoscape plug-in for MiMI. Tlr2 and Cd14 were associated with Ccnd1 by Egfr/insulin signaling molecules such as Egfr, Pik3r1, Prkaca, Src, and Stat3 (color figure online)

To check whether a similar PPI network exists in humans, we mapped the nine mouse genes associated with both the “inflammatory response” and “hepatic steatosis” to their human orthologs and built a human hepatic steatosis PPI network (Fig. 5). Similarly, Tlr2, Cd14, and Ccnd1 were associated through Egfr-mediated signaling molecules, such as Src, Pi3 k, ser/thr protein kinase, and Stat3 in the human hepatic steatosis PPI network (Fig. 5).

Figure 6 shows the dynamic changes related to hepatic steatosis PPI network based on time course microarrays of HFD-fed mice compared to ND-fed mice. During progression of diet-induced obesity, 6 genes (Ccnd, Cd14, Csf1r, Lgals3, Lpl, and Tlr2) among 7 were up-regulated in the hepatic steatosis PPI network.

Fig. 6

Time-dependent changes in central molecules of the hepatic steatosis PPI network which are involved in protein–protein interactions linking steatosis and inflammation in mice fed high-fat diet over 24 weeks

Discussion

We reveal dynamic changes in the core molecular networks underlying the development of diet-induced obesity. Importantly, we used a progressive diet-induced obesity model involving C57BL/6 J mice fed a physiological relevant high-fat diet (20 %kcal fat) over 24 weeks compared to age-matched normal diet-fed controls. Morphological analysis revealed hepatic fibrosis was evident after 16 weeks, while we previously reported hepatic lipid accumulation was evident after 2 weeks in this diet-induced obesity model (Do et al. 2011). Hepatic lipid accumulation is an early event in the development of non-alcoholic hepatic steatosis (NASH), which can later lead to hepatic fibrosis and eventually liver cirrhosis (Cohen et al. 2011). Inflammation is strongly associated with NASH development in animals and humans. In the present study, we identified Tlr2, Cd14, and Ccnd1 as hub genes associated with both hepatic steatosis and inflammation, and built a protein–protein interaction network to determine further interrelated proteins and biological functions.

Core molecular networks underlying the development of diet-induced obesity

To elucidate the core molecular networks in liver underlying the development of diet-induced obesity, high-fat diet responsive genes were analyzed using the Ingenuity knowledge base, which consists of expert curated molecular interactions. The high-fat diet responsive genes formed five core networks. Network 1 consisted of genes associated with movement, trafficking, and differentiation of immune cells, as well as fatty acid oxidation, lipid transport, and cell cycle. Tlr2, Ccl5, Cxcl9, Itch, and Cd14 were transcriptionally up-regulated, while Acacb and Apom were down-regulated consistently in Network 1 from week 2 to week 24 in HFD-fed mice. Acacb and Apom play a role in fatty acid oxidation and lipid transport. Acacb catalyzes the carboxylation of acetyl-CoA to malonyl-CoA, the rate-limiting step in fatty acid synthesis. Acacb is also reported to regulate fatty acid oxidation via malonyl-CoA inhibition of carnitine palmitoyltransferase I, the rate-limiting step in fatty acid uptake, and oxidation by mitochondria (Abu-Elheiga et al. 2001). Genes involved in the innate immune system and the trafficking of immune cells, including Tlr2, Cd14, Icam1, and various chemokines were consistently up-regulated in Network 1. Toll-like receptors (Tlr) are associated with liver diseases including alcoholic liver injury, ischemia/reperfusion liver injury, liver fibrosis, and liver cancer (Seki et al. 2007; Imaeda et al. 2009). Tlr2 recognizes conserved pathogen-associated molecular patterns and plays a critical role in the activation of the innate immune system in response to invading microorganisms (Akira 2003). In addition, blockade of Tlr2 signaling is reported to prevent insulin resistance in HFD-fed mice (Caricilli et al. 2008; Himes and Smith 2010). The role of Tlr4 and Tlr9 in NASH has not been well studied; nevertheless, Cd14 mediates the signaling of Tlr2 and Tlr4; therefore, Cd14-mediated Tlr activation may contribute to the metabolic complications associated with obesity via increased inflammation (Roncon-Albuquerque et al. 2008). Two critical cell cycle regulators, CyclinD1 (Ccnd1) and cyclin-dependent kinase inhibitor 1A (Cdkn1a) in Network 1, are also associated with innate immune system activation and were increased in a time-dependent manner in high-fat diet-fed mice.

Genes in Network 2 were related to lipolysis, lipid uptake and lipid transport, and interestingly involved in liver fibrosis, steatotic liver, and NAFLD (Shechter et al. 2003; Lydatakis et al. 2006; Guillén et al. 2009; Jinno et al. 2010; Yang et al. 2010). Tnf, Cidea, Acox1, S100a11, Ly6d, Lpl, and Idh1 were up-regulated in a time-dependent manner compared to normal diet—control. The genes in Network 2 appear to be regulated by peroxisome proliferator-activated receptor-gamma (Pparg), a critical metabolic regulator implicated in pathology of numerous diseases including obesity, diabetes, atherosclerosis, and cancer (Chawla et al. 2001; Libby 2001; Boitier et al. 2003). Pparg activates the formation of lipid droplets from triglycerides in hepatic tissues. The surface of the lipid droplet is coated with proteins of the perilipin family including Cide-a, Cide-b, and Cide-c, which can regulate lipid metabolism and mediate excess accumulation of intracellular lipids associated with obesity via hydrolysis of hepatic lipid droplets (Zhou et al. 2003; Li et al. 2007). Cidea-null mice are lean and resistant to diet-induced obesity and diabetes, partly due to a significantly higher metabolic rate compared with wild-type mice (Zhou et al. 2003).

Idh1 gene expression is regulated by sterols and activated by SREBP-1a and SREBP-2 in human hepatoma HepG2 cells and appears to regulate lipogenesis in hepatic cells (Shechter et al. 2003). Idh1 activity can be coordinately regulated with the cholesterol and fatty acid biosynthetic pathways. Idh1 gene expression was up-regulated from week 16, which may have accelerated hepatic lipogenesis, together with Lpl in our diet-induced obesity model. Fabp4 plays a role in fatty acid binding, uptake, transport, and was down-regulated in high-fat diet-fed mice after 20 weeks.

Ly6d, which marks the earliest stage of B cell specification and identifies the branch point between B cell and T cell development, was found to increase in a time-dependent manner during diet-induced obesity. Interestingly, Acox1 (acyl-CoA oxidase 1) and Lamb3 gene expressions were only up-regulated at week 24. Acox1, a rate-limiting enzyme in peroxisomal fatty acid beta-oxidation, was consistently down-regulated up to week 20 but was up-regulated at week 24. Interestingly, mice nullizygous for both Pparα and Acox1 lack spontaneous hepatic peroxisome proliferation and show a marked decrease in hepatic steatosis (Hashimoto et al. 1999). The time course of Lamb3 expression was very similar to that of Acox1, Lamb3 was consistently down-regulated up to week 20 but appeared to be increased at week 24. Lamb3 appeared to be one of the early genes closely involved in the development of liver fibrosis. Lamb3 expression in NASH patients with fibrosis is reported to be significantly higher than in NASH patients without fibrosis (Lydatakis et al. 2006). In the current study, morphological analysis of liver tissue indicated central venule fibrosis developed in HFD-fed mice after 16 weeks compared to ND-fed mice. S100a11 expression showed a time-dependent increase over 24 weeks, which was more prominent later during the development of diet-induced obesity. S100A11 is reportedly involved in tumor promoter activity as elevated protein level of S100A11 was found in hepatocellular carcinoma (HCC) cell lines with high metastatic potential. Thus, it is plausible that S100A11 may be an early biomarker of progression from liver steatosis to metastasis. Lpl gene expression was also consistently high over 24 weeks. Lpl plays roles in the clearance of plasma triglycerides from VLDL and hence may cause accumulation of hepatic triglycerides and FFA, which can be toxic to hepatocytes.

Genes in Network 3 were mainly hepatic-specific genes responsible for the homeostasis of cholesterols and bile acids such as cytochrome P450 families and nuclear receptor subfamilies (Fig. 3). Cyp7A1 gene expression was consistently up-regulated, and Cyp7A1 encodes cholesterol-7α-hydroxylase, a rate-limiting enzyme in the bile acid biosynthesis pathway that converts cholesterol into bile acids. Nr5a2 (Lrh-1, liver receptor homolog-1) is a bile acid transporter, which mediates bile acid homeostasis and steroidogenesis together with the Cyp7A1. Since Lrh-1 plays a key role in hepatic sterol metabolism as a bile acid transporter, several of its target genes are crucial to the regulation of reverse cholesterol transport, bile acid synthesis, and enterohepatic circulation (Fayard et al. 2004). Lrh-1 down-regulation is probably due to the decreased cholesterol biosynthesis and reverse cholesterol transport caused by high-fat feeding. Cholesterol feeding in rodents is reported to decrease sterol 12α-hydroxylase (Cyp8B1) and increase Cyp7A1, mediated by SREBP regulation (Le Martelot et al. 2009), indicating these are regulatory enzymes of bile acid synthesis. The dynamic transcriptional changes in Network 3 genes indicate hepatic sterol biosynthesis was decreased and bile acid biosynthesis was increased progressively in the HFD-fed mice, most likely attributable to excess cholesterol influx into liver from excess dietary intake.

Genes forming Network 4 were associated with the inflammatory response, cellular movement, and immunological disease. Network 4 was composed of pro-inflammatory factors such as Tnf, Ifng, Il1A, and Saa2. Network 4 also contained a few specific genes that were altered in a time-dependent manner. Gpnmb, Lgnm, Lyz, and Saa2 gene expression was up-regulated. Although Gpnmb is a novel osteoclastic protein that plays a key role in osteoclast differentiation and activity (Sheng et al. 2008), Gpnmb is also found in inflammatory liver macrophages (Haralanova-Ilieva et al. 2005). Lgnm encodes the legumain protein, which is involved in the processing of bacterial peptides and endogenous proteins via major histocompatibility complex (MHC) class II proteins in the lysosomal/endosomal degradation pathway. Accordingly, dynamic changes in genes in Network 4 may underlie the development of chronic inflammation, which occurs as a consequence of a long-term high-fat diet.

Genes in Network 5 were collectively associated with hepatic system disease, lipid metabolism, and molecular transport. Network 5 was composed of genes related to various metabolic processes of macromolecules such as gluconeogenesis (Ppargc1a), oxidation (Ppara), as well as biosynthesis (Me1) of fatty acids, cholesterols (Fdps), bile acids (Cyp7A1), and prostaglandin (Ptgs2). Expression of Fdps was decreased in a time-dependent manner along with its neighbors (Ppara, Ppargc1a, Scd, and Socs2; Fig. 3). The transcriptional profile of Ppara, Ppargc1a, Socs2, and Me1 were all decreased in time-dependent manner, which may be attributable to negative feedback due to elevated plasma and hepatic cholesterol, as well as hepatic triglycerides due to the high-fat feeding. The Ppara and Ppargc1a transcriptions were related to Fads, which regulates unsaturation of fatty acids, Ppara was also associated with Fads1 via Socs2. In rat islets, Socs2 is reported to suppress glucose-stimulated insulin secretion and Socs2 is related to increased susceptibility to T2D in humans (Kato et al. 2006). We found Nsdh, a non-specific aldehyde dehydrogenase, was down-regulated via Ppara in Network 5. Endogenous aldehydes are formed during the metabolism of alcohols, amino acids, biogenic amines, vitamins, steroids, and lipids. Network 5 also included Acadm (acyl-CoA dehydrogenase for medium fatty acids), an enzyme essential for converting specific fatty acids to energy during fasting particularly in the liver, Ptgs2 (prostaglandin-endoperoxide synthase 2) which plays a role in inflammation as indicated by severe nephropathy developed in Ptrs2-null mice (Norwood et al. 2000) and Gabpa(GA-binding protein alpha chain) which is involved in the activation of cytochrome oxidase expression and nuclear control of mitochondrial function (Virbasius et al. 1993; Virbasius and Scarpulla 1994).

Hepatic steatosis hub genes interact via a protein–protein interaction network

Diet-induced obese animals are characterized by early evidence of multiple obesity comorbidities, dependent on genetic background, diet composition, and duration. Hepatic steatosis and hepatic fibrosis are involved in the pathogenesis of NASH and were evident in our model of diet-induced obesity, despite modest high-fat diet intake. By merging the core diet-induced obesity networks, we identified Tlr2, Cd14, and Ccnd1 as hub genes related to both liver steatosis and the inflammatory response according to the IPA Knowledge base. We found these liver steatosis hub genes were interconnected at the molecular level via a protein–protein interaction network. Functional analysis of the protein–protein interaction network revealed Tlr2, Cd14, and Ccnd1 interact with protein clusters associated with the ErbB and insulin signaling pathway, lipid transport, and DNA damage response.

ErbB is a downstream effector of stress-induced insulin resistance, and stress-dependent transcription of ErbB2/ErbB3 receptors triggers a PI3 K cascade that induces serine phosphorylation of IRS proteins culminating in insulin resistance (Hemi et al. 2002). Insulin resistance is a pathological state in which target cells fail to respond to changes in circulating insulin and is associated with obesity, type 2 diabetes, as well as conditions of acute and chronic stress caused by high-fat feeding. The PPI network revealed the hepatic steatosis hub genes were interconnected via Egfr, Lgals3, and Csf1r, which form a protein cluster involved in ErbB and insulin signaling, as well as Lpl which belongs to a protein cluster associated with lipid transport. Lpl gene expression was consistently up-regulated in a time-dependent manner over 24 weeks, while Csf1r and Lgals gene expressions were up-regulated later after 16 week in high-fat diet-fed mice, conversely Egfr gene expression was down-regulated. Epidermal growth factor receptor (Egfr) is one of the proinflammatory factors which is reported to play a role in the activation of the carcinogenesis-associated signal transducer and transcription-3 (Stat3) (Starska et al. 2009). Cyclin D1 plays a critical role in transition from the G1 to S Phase of the cell cycle. Amplification, overexpression, and variations of the cyclin D1 gene have been detected in several types of cancers, including liver cancers (Zhang et al. 1993). Lgals3 is associated with cirrhosis in hepatocytes and has been reported not only to activate the expression of cyclin D1 but also to regulate cell proliferation (Matsuda et al. 2008). Lgals3 and Csf1r appeared to contribute to ErbB/insulin signaling pathways via maintaining the negative association with Egfr (Merlin et al. 2011), as well as recruiting insulin signaling molecules (Hamilton 1997).

We acknowledge the hepatic steatosis PPI network we present here remains to be verified in diet-induced obesity in vivo. Nevertheless, our findings highlight, hub genes related to steatosis, a common obesity comorbidity, appear to be highly connected via underlying protein–protein interaction networks. Therefore, perturbation of these hub genes may have far wider functional effects on insulin signaling, lipid transport, and DNA damage. Conversely, small molecule inhibitors of nodes in PPI networks may lead to dynamic changes in hub gene expression. A major strength of the current study was the use of a long-term diet-induced obesity model which recapitulated morphological changes in liver observed in human NASH, such as hepatic lipid accumulation, inflammation, and fibrosis. We used a physiologically relevant high-fat (20 %kcal) diet to induce obesity and NASH over 24 weeks, rather than using excessively high-fat (80 %kcal) diets suggested by recent murine NASH models (Hebbard and George 2011; Ogasawara et al. 2011). In future studies, it may be worthwhile to identify hub genes and protein–protein interaction networks related to the development of other obesity comorbidities including insulin resistance and type 2 diabetes.

Conclusion

High-fat diet responsive genes form several core networks underlying the development of diet-induced obesity. We identified Tlr2, Cd14, and Ccnd1 as major hub genes associated with both hepatic steatosis and inflammation. Tlr2, Cd14, and Ccnd1 are interconnected via a protein–protein interaction network, which form clusters associated with the ErbB/insulin signaling pathway, lipid transport, and DNA damage. Therefore, Tlr2, Cd14, and Ccnd1 represent possible targets to prevent the development of NASH in diet-induced obesity. However, further study is needed to elucidate the exact roles of these proteins in the development of NASH and progression to liver cirrhosis.