Gene Expression Profiling and Network Analysis Reveals Lipid and Steroid Metabolism to Be the Most Favored by TNFα in HepG2 Cells

Background The proinflammatory cytokine, TNFα, is a crucial mediator of the pathogenesis of several diseases, more so in cases involving the liver wherein it is critical in maintaining liver homeostasis since it is a major determiner of hepatocyte life and death. Gene expression profiling serves as an appropriate strategy to unravel the underlying signatures to envisage such varied responses and considering this, gene transcription profiling was examined in control and TNFα treated HepG2 cells. Methods and Findings Microarray experiments between control and TNFα treated HepG2 cells indicated that TNFα could significantly alter the expression profiling of 140 genes; among those up-regulated, several GO (Gene Ontology) terms related to lipid and fat metabolism were significantly (p<0.01) overrepresented indicating a global preference of fat metabolism within the hepatocyte and those within the down-regulated dataset included genes involved in several aspects of the immune response like immunoglobulin receptor activity and IgE binding thereby indicating a compromise in the immune defense mechanism(s). Conserved transcription factor binding sites were identified in identically clustered genes within a common GO term and SREBP-1 and FOXJ2 depicted increased occupation of their respective binding elements in the presence of TNFα. The interacting network of “lipid metabolism, small molecule biochemistry” was derived to be significantly overrepresented that correlated well with the top canonical pathway of “biosynthesis of steroids”. Conclusions TNFα alters the transcriptome profiling within HepG2 cells with an interesting catalog of genes being affected and those involved in lipid and steroid metabolism to be the most favored. This study represents a composite analysis of the effects of TNFα in HepG2 cells that encompasses the altered transcriptome profiling, the functional analysis of the up- and down- regulated genes and the identification of conserved transcription factor binding sites. These could possibly determine TNFα mediated alterations mainly the phenotypes of hepatic steatosis and fatty liver associated with several hepatic pathological states.


Introduction
The multifunctional proinflamatory cytokine, TNFa (Tumor Necrosis Factor a) has effects on lipid metabolism, coagulation, insulin resistance, and endothelial function and therefore is important in cellular pathophysiology. It induces the activation of predominantly two signaling pathways within the cell -one involving the transcription factor NF-kB (Nuclear Factor-kB) that lies inactive in the cytosol and is activated on TNFa binding to its receptor on the membrane and the other involves the MAPK (Mitogen Activated Protein Kinase) pathway [1][2][3]. Both these pathways underlie the basic mechanism(s) of TNFa action in several pathological conditions. TNFa is widely implicated in various diseases. Although primarily known to be associated with inflammation, TNFa now occupies a central place in the crosstalk of many pathological states. It plays a major role in muscular abnormalities resulting in muscle wasting [4], idiopathic pulmonary fibrosis [5], rheumatoid arthritis [6], Crohn's disease [7] and psoriasis [8]. Several other diseases like cerebral malaria [9], asthma [10,11] and cystic fibrosis [12] also have a strong correlation to TNFa levels. Chronic inflammation is also frequently associated with the metabolic syndrome and TNFa is one of the significant mediators between the two states [13]. The metabolic syndrome is a cluster of metabolic abnormalities leading to increased risk for cardiovascular diseases and diabetes type 2 and herein visceral obesity and the resulting insulin resistance are the major determinants in the development of this syndrome [14].
The liver is central to metabolic control and actively participates in storing and releasing glucose, nitrogen metabolites and lipids as and when required. It is also majorly involved in detoxification mechanisms and is, therefore, no doubt significant for proper functioning of the body and any liver injury is extremely detrimental in preserving normal hepatic homeostasis. Liver injury and death are the major hallmarks in all liver diseases [15] and it is during these conditions that the liver is exposed to increased levels of TNFa apart from other inflammatory cytokines. Serum and hepatic TNFa levels are increased in patients with HBV and HCV infection [16][17][18] and also in subjects with detectable liver failure and other liver injuries [19,20]. In all these cases, TNFa along with other cytokines participate in initiating a wound healing response [15] or mediate hepatocyte damage [21,22]. Paradoxically, although harmful to the injured hepatocyte, TNFa induces healthy hepatocytes to proliferate rather than die [23]. This especially occurs as a conditioned response against any hepatectomy or liver regeneration. Thus to maintain normal liver integrity and homeostasis by either hepatocyte death as observed during liver injury or by inducing hepatocyte proliferation during liver regeneration depends on whether TNFa-induced apoptosis or cell survival pathways, respectively are predominant. All these imply that TNFa is very critical in liver function and its increased levels exert several pleiotropic effects thereby emphasizing the need for a detailed understanding of hepatic TNFa action [15].
In view of these varied effects of TNFa, we sought to study the TNFa-induced alteration in gene expression patterns within the hepatocyte, the common pathways they fall in and to find a correlation among co-regulated genes. Our results demonstrate a significant number of hepatic genes to be statistically altered by TNFa and the most significant pathway and network to be targeted predominantly involved those of lipid metabolism, cholesterol and steroid biosynthesis. Interestingly, conserved transcription factors' binding sites were found in these sets of genes that could determine the correlation between the altered gene set and the hepatic processes favored by TNFa. Although earlier studies [24,25] have shown some specific effects of TNFa on HepG2 cells, our results demonstrate a wholesome analysis of TNFa effect on HepG2 cells that correlates the altered transcriptome profile with identification of conserved transcription factor binding sites that could determine several hepatic pathological states.

Cell Culture and Treatment
Experiments were done in HepG2 (human hepatocellular carcinoma) cells obtained from the National Centre for Cell Science, Pune, India. HepG2 cells retain several functions of the normal human liver including synthesis of albumin, lipoprotein and several other liver specific functions. Among these functions notable are those related to cholesterol and triglyceride metabolism. In several hepatic physiological states, inflammation within the liver is a significant accompanying phenomenon along with the accumulation of lipid vesicles. Dysregulation of liver metabolism and hepatic inflammation run in parallel in several liver diseases and thereby inflammation within the liver predisposes it to the onset of a number of hepatic pathologies. HepG2 cells are recognized as an in vitro human model system and are identified as a possible bioartificial liver [26] thereby corroborating their use as a resource to study several hepatic functions and for metabolic studies [27]. They were maintained in DMEM supplemented with 10% fetal calf serum with 1% antibiotic-antimycotic (100 units/ml penicillin, 0.1 mg/ml streptomycin and 0.25 mg/ml amphoterecin B) at 37uC in a humified atmosphere of 5% CO 2 . On attaining confluence, cells were serum starved overnight and incubated in the absence (control) and presence of TNFa (Sigma Chemical Co., St. Louis, USA; 0.5 nM, 12 h). The dose and time were chosen becausein HepG2 cells, we had observed optimum inhibition of insulin-stimulated Akt activation that plateaued thereafter at further dose and time points. The total levels of the Akt protein, however, remained constant ( Figure S1).

RNA Isolation
On termination of incubation, total RNA was isolated from cells using Trizol (Invitrogen) according to the manufacturer's instructions and subjected to a cleanup using RNeasy protect Mini Kit (Qiagen, Germany). RNA was quantified in a NANODROP spectrophotometer (ND 1000) and the ratios of OD260/280 nm were between 1.8-2.0. The integrity of RNA samples was determined in the Bioanalyser 2100 System (Agilent Technologies, Palo Alto, CA) and all samples had RIN (RNA integrity number) values between 8.0-8.5.

Microarray and Data Analysis
Microarray experiments were performed in triplicate (three control and three TNFa-treated) at The Centre for Genomic Applications (TCGA), New Delhi, India according to the manufacturer's instructions (Affymetrix, Santa Clara, USA). Briefly, RNA (control and TNFa treated) was reverse transcribed to cDNA that was subjected to in vitro transcription to produce biotinylated cRNA that was hybridized to the human array chip (Human Genome U133 Plus 2.0, Affymetrix). Images were scanned using the GeneChip 3000 7G scanner (Affymetrix). Signal intensities and absolute call dataset were generated with Affymetrix Gene Chip Operating Software (GCOS) using the MAS5.0 algorithm. This resulted in normalized values of each chip to a target intensity of 500 as detailed in the Affymetrix statistical algorithm. The efficiency of the labeling reaction and hybridization performance was assessed by the optimal 39/59 hybridization ratios for the housekeeping genes, poly-A-spike in controls and the prokaryotic controls. Detection calls identified transcripts to be either reliably detected (Present) or below the background levels (Absent) if the detection p-value was above or below 0.05 respectively. Genes with present calls in at least 2 out of the three replicates were included in the analysis that was done with Microsoft Excel 2007 using the steps as described in http:// jura.wi.mit.edu/bio/education/bioinfo2007/arrays/. The data were normalized such that the average expression of the present genes in all the samples was identical, log transformed and statistical analysis was done with student's t-test. Genes were filtered on the basis of call, fold change and p value. All genes with a signal log ratio .60.5 (1.41 fold) and p value ,0.05 were chosen to be statistically altered by TNFa. The 1.41-fold criterion as a minimum cutoff for change is chosen based on the literature reported earlier in HepG2 cells [28]. The study was in accordance to the MIAMI guidelines and the microarray data have been submitted to the Gene Expression Omnibus (GEO) database [29]. Accession number for the data is GSE12161.

Real Time PCR
Real Time quantitative PCR was done for validation of gene expression microarray data. Total RNA (Control and TNFatreated) was isolated as described above and 2 mg RNA from each set was reverse transcribed using MMuLV reverse transcriptase and random hexamers. 2 ml of this was amplified in a final volume of 25 ml using the SYBR Green Real Time PCR Master Mix (Applied Biosystems) and gene specific primers. After an initial activation at 50uC for 2 min and a preincubation step at 95uC for 10 min, reaction mixtures were subjected to 40 cycles of denaturation at 95uC for 15 sec, annealing at 56uC for 30 sec and extension at 72uC for 35 sec (ABI 7500). This was followed by a dissociation curve analysis at 95uC for 15 sec, 60uC for 20 sec and 95uC for 15 sec. Data analysis was done as described by Pfaffl, 2001 [30] and 18S rRNA was taken as the internal control. No RT control (without the reverse transcriptase) and a no template control (without the RNA template) were taken as the negative controls for the experiments. All experiments were done in triplicate.

Gene Ontology (GO) Classification and Cluster Analysis
Genes that were found to be significantly altered by TNFa were classified into categories of biological processes and molecular functions using Gene Ontology Tree Machine (GOTM, [31]) and DAVID [32]. Genes over-represented in our dataset in the context of the reference set (Human Genome U133 Plus 2.0 chip, Affymetrix) were determined with the GO Tool Box [33] using Bonferroni correction. Genes were then clustered on the basis of their normalized expression intensities using the Avadis software (version 3.1, Stratagene).

Promoter and Pathway Analysis
To identify the likely identical regulatory controls among the altered genes, genes common to a single cluster and Gene Ontology term were taken for promoter analysis. 5.0 kb upstream regions from the transcription start sites of such genes were retrieved using Ensembl [34] and the transcription factor(s) binding within this region was determined using Over-represented Transcription Factor Binding Site Prediction (OTFBS, Release 6.0.) tool [35] that detects over-represented motifs of known transcription factors within groups of related sequences. This was done considering the fact that genes with similar function as well as expression profiles are more likely to be under the same regulatory control than genes within either of the categories [36]. Also to categorically determine which specific pathways the altered genes map to and to determine the interacting network among them, if any, the dataset containing the gene identifiers alongwith the corresponding fold changes were uploaded into the Ingenuity Pathways Analysis program. Each identifier was mapped to its corresponding gene object in the Ingenuity knowledge base. A 1.41 fold change (log ratio 60.5) was set to identify genes whose expression was significantly differentially regulated. These genes, called focus genes, were overlaid onto a global molecular network developed from information contained in the Ingenuity knowledge base. Networks of these focus genes were then algorithmically generated based on their connectivity. A graphical representation of the molecular relationships between these genes is shown where genes or gene products are represented as nodes, and the biological relationship between two nodes is represented as an edge (line). All edges are supported by at least 1 reference from the literature, from a textbook, or from canonical information stored in the Ingenuity knowledge base. Human, mouse, and rat orthologs of a gene are stored as separate objects in the Ingenuity knowledge base, but are represented as a single node in the network. The intensity of the node color indicates the degree of up-(red) or down-(green) regulation. Canonical Pathways Analysis identified the pathways from the Ingenuity Pathways Analysis library of canonical pathways that were most significant to the dataset. Genes from the dataset with a 1.41 fold change (log 2 ratio 60.5) were associated with a canonical pathway in the Ingenuity knowledge base and considered for the analysis. The significance of the association between the dataset and the canonical pathway was measured in 2 ways: 1) A ratio of the number of genes from the dataset that met the expression value cutoff that map to the pathway divided by the total number of molecules that exist in the canonical pathway is displayed. 2) Fischer's exact test was used to calculate a p-value determining the probability that the association between the genes in the dataset and the canonical pathway is explained by chance alone.

Electrophoretic Mobility Shift Assay (EMSA)
The electrophoretic mobility shift assay analysis was done essentially as described by Shiraga et al. [37]. Briefly nuclear extracts were prepared from control and TNFa-treated (0. 5  . Either the wild type or mutated labeled double stranded oligonucleotide (40,000 cpm) was then added to the reaction mixture and incubated for 45 min at room temperature. On termination of incubation, samples were loaded onto a non-reducing 5% polyacrylamide gel and electrophoresed in 0.5XTBE. Gels were then dried and subjected to phosphorimager analyses (FLA 2000, Fujifilm, Japan). The densitometric analyses were done using the Alpha DigiDoc 1201 software (Alpha Innotech Corporation, CA, USA). The same size rectangle box was drawn surrounding each band and the intensity of each was analyzed by the program after subtraction of the background intensity.

Cholesterol Synthesis
On attaining confluence, HepG2 cells were incubated in the absence or presence of TNFa (0.5 nM, 12 h). Additionally to determine which specific step of the steroid biosynthetic pathway is the most significant contributor for TNFa action herein, gene specific siRNAs were used to categorically knock down specific genes. HepG2 cells were transfected with 50 nM of either control or specific siRNAs (Santa Cruz, CA, USA) according to the manufacturer's instructions. After allowing the cells to grow in fresh DMEM for 48 h, they were incubated with TNFa (0.5 nM, 12 h) and on termination of incubation, cells were lysed and lipids were extracted by addition of 1.0 ml of chloroform: methanol (2:1). These lipid extracts were vacuum dried and the residues were dissolved in 2-propanol containing 10% Triton-X-100. Cholesterol was estimated in an analyser (COBAS INTEGRA 400 plus, Roche, Germany) according to the manufacturer's instructions. Results are expressed after normalization to the total protein content (DC Protein Assay Kit, BIORAD, CA).

Data Analysis
Microarray data from three replicate experiments generated by GCOS were exported to Microsoft Excel 2007 for further analysis and is described above. A value of p,0.05 was taken to be statistically significant. Only those genes with present calls in at least 2 out of 3 replicates were analyzed and a cutoff value for differential gene expression was defined as transcripts that showed a.1.41 fold change and a t-test p value of ,0.05. RT-PCR was done in triplicates and data presented are means6 S.E.M of three independent experiments. Statistical significance between the control and treated groups was determined using Student's t-test where p#0.05 was considered as statistically significant.

TNFa Induced Alteration in Gene Expression in HepG2 Cells
To identify genes altered by TNFa in HepG2 cells, mRNA expression in these cells incubated in the absence and presence of TNFa (0.5 nM, 12 h) was analysed using whole genome oligonucleotide expression arrays (Human Genome U133 Plus 2.0, Affymetrix). TNFa treatment significantly (p,0.05) altered the expression of 140 genes in HepG2 cells. Figure 1 shows the volcano plot of the microarray hybridization with the significantly up-regulated genes being represented with red colored boxes and those down regulated with green ones. Those that were not significantly altered are depicted in black dots. These altered genes were classified into functional categories of biological processes using Gene Ontology Tree Machine (GOTM) and DAVID. The GO database has a hierarchical nature such that genes annotated with a specific node (GO term) are also annotated with every ancestor of that node with the result that the genes map onto more than one GO term. A functional GO class thereby included genes assigned in the listed GO term by the GO consortium or a more specific break-up of the listed term that reflects in the gene being classified in more than one GO term. From our dataset, genes that clearly mapped onto a single GO term were assigned as such and for others that mapped onto more than one GO term, we put genes in the term where there was experimental validation available and when not, we put them in a particular GO class considering their functional closeness to other members of that class. Table S1 gives the Gene Ontology terms, the fold change and the p value of the significantly altered genes. Of these, 67 genes were up-regulated and 73 were down-regulated with a minimal fold change of 61.41 (log 2 ratio $60.5). Genes involved in cell adhesion like Laminin (gamma 1, 21.4 fold), Tao kinase2 (21.6 fold) and Fibroblast Growth Factor Receptor-Like 1 (21.5 fold) were significantly down regulated by TNFa. Almost all the genes involved in response to stimulus (HSPA1B, ORM2, MST1, H3/O and C4B) and regulation of biological processes (PUM1, CEBPa and BRPF1) were down regulated. While some genes of signal transduction like FCER1G, RASSF7, NR1H4 and RGL3 were significantly (p,0.05) down regulated, others of the same category like MKNK2, RYK, NFKB2, IL22RA1, THRAP6 and KISS1R were upregulated. Further upregulated genes were involved in transport such as GOSR2 (+1.4 fold), HCN3 (+1.7 fold), SYTL1 (+1.4 fold), SLC43A2 (+1.9 fold) and SLC27A5 (+2 fold) and lipid and carbohydrate metabolism like FDPS (+1.6 fold), PAFAH2 (+2 fold), FADS (+1.7 fold), SQLE (+2.5 fold), HSD17B7 (+1.7 fold) and IDH2 (+1.5 fold). While genes involved in carbohydrate metabolism, for e.g. AKR1B1, GNPDA2, and PDK3 were significantly downregulated, some genes of nucleic acid and amino acid metabolism like PAH, NDUFS1, DHFR and TNRC9 were upregulated. Others like STOM, UBE2J1, SCRN3 (protein metabolism), SIRT6, ELF5, DNMT3B (nucleic acid metabolism), AP3B1, SLC6A12, STXBP5, AP1S2, AQP3 (transport) and TIMP1 were significantly downregulated. Interestingly genes involved in cellular homeostasis i.e. calreticulin (21.6 fold) and protein disulfide isomerase (21.9 fold) were also downregulated.

Real Time PCR Validation of Microarray Data
Quantitative real time PCR was done to validate a group of selected up-and down-regulated genes. The primer sets used are shown in Table 1. The experiment was done for FADS1, SQLE, PDIA4, HSP90B, CEBPa, CALR, FDPS, MST1, PDK3, AQP3 and NFKB2. As in the results of the microarray analysis, all these genes depicted an identical pattern of alteration as compared to the control. In TNFa treated cells, PDIA4, HSP90B1, CEBPa, CALR, MST1, PDK3 and AQP3 showed an inhibition of the transcript levels while FADS1, SQLE, FDPS and NFKB2 showed increased levels ( Figure 2).

Functional Classification of Genes Altered by TNFa in HepG2 Cells
All the 67 up-regulated and 73 down-regulated genes were separately grouped into molecular functions using the GO Tool Box. At Level 3, molecular functions that emerged specifically as being altered among the up-regulated genes were those of lipid binding and ligase activity while those of transcription repressor and enzyme inhibitor activities were specific to the down-regulated genes. Figures 3A and B depict respectively the fraction of TNFa induced up-and down-regulated genes in HepG2 cells among the various molecular functions as identified by the GO Tool Box. Significantly (p,0.01) over represented Gene Ontology terms (molecular function) in down-and up-regulated genes are shown in Tables 2A and B respectively. The tables give the number of genes with that gene ontology term on the Human Genome U133 Plus 2.0 array, the number of genes in the same gene ontology term in our altered dataset and the p value of the occurrence of that number of genes in the dataset. For genes up-regulated by TNFa, several gene ontology terms related to lipid and fat metabolism like keto-steroid reductase activity, cholate-CoA ligase activity, squalene monooxygenase activity, oxidoreductase activity were significantly (p,0.01) over represented indicating an overall favor of genes involved in fat metabolism within the hepatocyte by TNFa. GO terms involved with immune responses like IgE receptor activity, immunoglobulin receptor activity and IgE binding were also highly over represented in the down-regulated set indicating a compromise in immune defense mechanism(s). An identical analysis using the GO class: biological processes identified several genes belonging to the steroid and lipid metabolic processes to be significantly (p,0.01) over-represented among the up-regulated set and those very specifically involved in diverse aspects of the overall immune response to be significantly (p,0.01) over represented in the down-regulated set of genes (Table S2).

Clustering and Promoter Analyses
TNFa induced up-and down-regulated genes were clustered hierarchically on the basis of their normalized expression intensities ( Figure S2) and genes with the same GO functional term and within the same cluster were selected for further promoter analysis. Example results of genes from three GO terms i.e. ''Metabolism'', ''Signal Transduction'' and ''Regulation of Biological Process'' that clustered together were taken for promoter analysis. A 5kb upstream sequence from the transcription start site of each of these genes was retrieved and the binding sites of transcription factors within this region were determined as described above. Table 3 is a list of the transcription factor binding sites in a set of clustered genes within the indicated GO term. Significant transcription factors consistently occurring in the promoters of genes altered by TNFa and clustering together with common GO terms included SREBP-1, CEBPa, MEF2 and AREB6. All these are in several ways implicated in cellular metabolic processes and immune response, two phenomena that are modulated by TNFa. Surprisingly, binding sites of several members of the forkhead family of transcription factors, that are however not well characterised were also significantly found to be overrepresented in our dataset. These included FOXJ2, FOXD3 and the fork head homolog-3 (HFH-3). The biological relevance of these relatively less characterized transcription factors remain to be investigated.

EMSA for Validation of Candidate Transcription Factors
To validate the transcription factors that were found to possess binding sites in the promoters of genes that were altered by TNFa as described above, electrophretic mobility shift assay was performed for two of the candidate transcription factors namely, SREBP-1 and FOXJ2. Wild type and mutated oligonucleotide sequences that were used are given in Table 4. As compared to the incubation in the presence of the free labeled probe alone, the addition of the nuclear extract from the control cells caused significant formation of the DNA-protein complex both with SREBP-1 ( Figure 4A, B) and FOXJ2 ( Figure 4C, D) oligonucleotides. This complex formation was significantly enhanced in the presence of nuclear extracts from TNFa-treated cells indicating that TNFa increased the binding of these transcription factors to their binding elements that corroborates with our predicted findings detailed above. No detectable DNA-protein complexes were identified with the mutated oligonucleotides of either SREBP-1 or FOXJ2 that additionally suggests the specificity of the DNA-protein complexes observed.

Network and Pathway Analyses
We then wanted to identify the interacting network among the altered genes in the context of other large biological pathways. We began by analyzing TNFa-induced altered genes using the Ingenuity Pathway Analysis software that generates networks, associated biological functions and canonical pathways by overlaying them onto a global molecular network developed from information contained in the Ingenuity knowledge base. Five highly significant networks with scores above 20 were obtained from the set of TNFa induced altered genes. These scores, derived from p values, indicated the likelihood of the focus genes belonging to a network versus those obtained by chance alone thereby eliminating the probability of their occurrence in a network to be due to noise. The network with the highest score comprised that of hematological system development and function, immune and lymphatic system development and function and tissue morphology. Other subsequent networks with decreasing scores included those of (i) immunological disease, skeletal and muscular disorders and inflammatory disease (ii) lipid metabolism, small molecule biochemistry, nervous system development and function, (iii) cell signaling, post translational modification, dermatological diseases and conditions and (iv) cellular growth and proliferation, cancer and hematological disease. This implies that the altered genes were distributed among diverse networks that could be expected considering the varied cellular actions of TNFa. Interestingly, while ''biosynthesis of steroids'' was the top canonical pathway that was significantly (with a p value of 1.27E-03) altered by TNFa in HepG2 cells, ''lipid metabolism'' was the most significant (1.70E-04-4.39E-02) molecular and cellular function to be affected. So we further analyzed the network 3 of the list of significantly affected networks i.e of lipid metabolism, small molecule biochemistry, nervous system development and function since it could also be correlated well with the most significantly affected canonical pathway and molecular and cellular function. Figure 5 shows the associated network functions of this interaction. One of the key central points in this network is that of insulin (INS1). Apart from regulating FADS1, INS1 indirectly also regulates HRAS that inturn regulates SQLE, these are the two significantly upregulated genes in our dataset. Insulin was postulated some time ago by Satoh et al., 1990 [38] to regulate SQLE and insulin regulation of hepatic desaturases has also been described [39]. INS1 is directly regulated by HNF1A that in turn interacts with AQP3 and PAFAH2 that were down and up regulated, respectively, in our experiment. HRAS shows interactions with SQLE, PCYT1A and 1B, genes that were differentially regulated by TNFa. Another central point in this interaction network was occupied by RHOA (ras homolog gene family, member A). Two of the significantly altered genes from our list, HSPA1A and RGNEF, were found to be binding strongly with RHOA. Interestingly, RHOA is also regulated by HRAS on one end and INS1 on the other, although indirectly. All these indicate a cross talk among molecules that otherwise seem to be distantly placed in biological pathways and therefore many positions of this  Table 2. Values of each gene were normalized to those of 18S rRNA and are plotted with respect to the control that has arbitrarily been taken as 1. Each value is the mean 6 SEM of three independent experiments. *p,0.01 as compared to control. **p,0.05 as compared to control. doi:10.1371/journal.pone.0009063.g002 network being occupied by genes altered by TNFa possibly explains to some extent the diverse effects of TNFa on cellular physiology. Figure 6 shows the ''biosynthesis of steroids'' pathway obtained using the Ingenuity Pathway Analysis tool that was derived as the most significant pathway being affected by our set of TNFa induced altered genes. Most of the affected genes clustered at significant points of the cholesterol biosynthetic pathway. These include three isoforms of farnesyl diphosphate synthase (FDPS, EC nos. 2.5.1.1; 2.5.1.10; 2.5.1.29), squalene expoxidase (SQLE, EC no. 1.14.99.7) and emopamil binding protein (EBP, sterol isomerase, EC. no. 5.3.3.5). All these were upregulated and have been marked with red diamond boxes in the figure with the intensity of the color being an indicator of the fold of upregulation. Two of these namely SQLE and FDPS were also validated by Real-Time PCR (Figure 2). Phenotypically, as a  consequence of this, we observed enhanced cholesterol accumulation in HepG2 cells on treatment with TNFa (0.5 nM) for 12 h ( Figure 7A). In an attempt to understand which of the above TNFa-altered genes of this pathway is the most significant contributor for the observed increase in the cholesterol levels in the presence of TNFa, siRNAs were used to knock down the specific genes namely SQLE, FDPS and EBP. Use of the siRNAs could knock down the levels of the respective proteins by 75% (data not shown). Figure 7B shows the effect of TNFa on the cellular cholesterol levels in the presence of the siRNAs. In cells transfected with SQLE or FDPS siRNAs, the TNFa-mediated increase in the cellular cholesterol content was significantly (p,0.01) prevented with SQLE siRNA demonstrating a higher fold of prevention. EBP siRNA however depicted a very mild preventive effect. All these indicate that of these three enzymes of the steroid biosynthetic pathway that are up-regulated by TNFa, SQLE is the most significant contributor of the elevated synthesis of cholesterol by TNFa in HepG2 cells. This pathway emerging as the top canonical pathway is therefore indicative of a correlation between cellular lipid homeostasis and the inflammatory molecule that may possibly underlie the pathogenesis of several diseases.

Discussion
TNFa is a major proinflammatory cytokine that activates varied cellular signaling pathways resulting in diverse biological responses. Within the liver, TNFa is known to regulate contrasting processes. On the one hand, during liver regeneration for e.g.after partial hepatectomy TNFa induces hepatocyte proliferation thereby regulating normal cell function while on the other hand, in the presence of any hepatocyte injury, TNFa production can also trigger cell death rather than proliferation. TNFa therefore is a crucial deciding factor between hepatocyte life and death and therefore of global liver homeostasis. In the context of these paradoxically contrasting hepatic responses towards TNFa, we employed high density oligonucleotide arrays (Human Genome U133 plus 2.0, Affymetrix) to identify hepatic genes altered by TNFa and elucidate the common transcriptional regulatory controls among the altered genes.
Microarray analysis of HepG2 genes altered by TNFa identified 67 genes to be up-and 73 genes to be down-regulated. Very significantly and consistently altered were genes of lipid metabolism mainly squalene epoxidase, fatty acid desaturase, farnesyl diphosphate synthase and hydroxysteroid (17-b) dehydrogenase 7 thus indicating lipid metabolism to be crucially targeted by TNFa. This suitably agrees with the fact that TNFa is very significant in the pathogenesis of fatty liver diseases that accompanies a host of metabolic disorders. TNFa has been shown to induce fat accumulation in the liver resulting in hepatic steatosis [40]. In fact, Ruan et al., [41] also reported that squalene epoxidase and farnesyl diphosphate synthase were increased by 2.4 and 1.7 fold respectively in the rat liver on infusion with TNFa. The transcription factor, CEBPa is implicated in the regulation of numerous liver-specific genes and it is a critical regulator of hepatic metabolism [42][43][44][45][46] and our results are in agreement with earlier reports [47,48] that show a significant inhibition of CEBPa by TNFa thereby implying a major alteration of cellular metabolism.
The transcription factor, ATF5 binds to the cAMP response element (CRE) on cAMP inducible promoters and regulates Table 3. Over representation of conserved transcription factor binding sites within putatively co-regulated genes. GO   several hepatic specific enzymes [49] and its loss accompanies hepatic carcinomas [50]. Not reported earlier, our demonstration of ATF5 being regulated by TNFa in HepG2 cells, therefore, is indicative of its participation in the cellular role of TNFa in hepatic physiology. Our results also demonstrated that genes involved in cell death and cell proliferation were significantly altered by TNFa. Earlier reports also demonstrated the same with the transcription factor NF-kB and c-jun playing a significant mediatory role therein [51,52]. TNFa acts primarily through the NFkB and MAPK pathways and their activation is observed immediately as initiators of downstream signaling cascades. In our study, we observed activation of two downstream molecules of these pathways namely, NFKB2 (Nuclear Factor of Kappa light Polypeptide Gene Enhancer in B-cells 2) and MKNK2 (MAPK interacting serine/threonine kinase 2). NFKB2 is a subunit of NFkB that has been shown to be increased at the mRNA and protein level by TNFa [53] and its inappropriate activation is linked to varied pathophysiological states. MKNK2 contains a MAPK binding and activating motif and it participates in the regulation of protein biosynthesis and translation, in the protein kinase cascade and in response to stress. The involvement of the MAPK pathway in TNFa-induced toxic liver injury was reported in an in vivo mouse model [54] and our study has demonstrated an up-regulation of MKNK2 by 1.8-fold in HepG2 cells by TNFa.
Also, MKNK2 is known to interact with TRAF2 (TNF receptorassociated factor 2) and thereby play a significant role in TNFa receptor activity associated cellular phenomena. Genes involved in ''Cell Adhesion'' were significantly down regulated implying a compromise in these phenomena by TNFa within the hepatocyte. In our study we found laminin gamma 1 (LAMC1) and TIMP 1 to be significantly down regulated (21.4 fold and 21.5 fold respectively) by TNFa indicating a modulation of the extracellular matrix. The regulation of laminin, TIMP1 and MMP9 by TNFa has been described [55][56][57][58][59] and here our results show decreased expression of TIMP1 that together with the decreased laminin gamma 1 expression indicate towards a modulation of the extracellular matrix and thereby of cell adhesion within the hepatocyte by TNFa.
Interestingly, in the present study, a catalog of genes not earlier reported to be regulated by TNFa were found to be altered in HepG2 cells. Significant among these were those from the GO term: nucleic acid metabolism mainly DHFR, MTHFS, SIRT6 and DNMT3B; almost all of these are involved in the synthesis, degradation or function of nucleic acids. Nucleic acids and their metabolites are critical in the elucidation of an immune response [60] and the alteration of these genes by TNFa as identified in our study could at least in part possibly explain the role of nucleic acids in the effects of TNFa in HepG2 cells. Also several altered genes that categorized within the GO term ''Signal transduction'' like RASSF7, RGL3, KISSIR and RYK are being reported here for the first time to be regulated by TNFa. RASSF7 and RGL3 are significant signaling intermediates of the RAS pathway and their alterations indicate towards modulation of this important cellular signaling cascade by TNFa in HepG2 cells. Another important molecule affected by TNFa is the gene RYK receptor like tyrosine kinase (RYK) that like RTKs are deregulated in cancers [61] and our revelation of the regulation of its gene expression by TNFa in HepG2 cells could provide information surrounding these phenomena. Our study, therefore, apart from demonstrating the effects of TNFa on genes already known to be regulated as such, has also unearthed several genes that were hitherto not known to be regulated by TNFa.
Several well characterized and other not so well characterized transcription factors were identified with binding elements in sets of co-regulated genes. Myocyte enhancer factor 2 or MEF-2 family of transcription factors that bind as homo-and heterodimers to the MEF2 sites in the promoter regions of target genes were found to be the most enriched in genes that belonged to the ''Regulation of Biological Processes'' GO term (Table 3) suggesting that this may be a major determiner of these processes within the hepatocyte.
Sterol Regulatory Element Binding Proteins (SREBPs) are a basichelix-loop-helix leucine zipper class of transcription factors that bind to the sterol regulatory element DNA sequence TCACNC-CAC. Within the liver, this transcription factor is known to regulate several hepatic processes mainly gluconeogenesis [62], lipid and cholesterol biosynthesis [63] and thereby attains significance in diverse metabolic disturbances. TNFa has been shown to stimulate the maturation of SREBP-1 in hepatocytes [64] and our results demonstrated herein imply that SREBP-1 could be important in mediating the effects of TNFa on these genes.
Interferon regulatory factor 1 or IRF-1 is a member of the interferon regulatory transcription factor (IRF) family and activates interferon transcription as well as of genes induced by interferons [65]. IRF1 exhibits regulatory roles in diverse cellular processess [66,67] and it has been associated with some forms of cancer mainly lung and gastric cancer. Our demonstration of IRF-1 depicting a high probability of over-representation in genes that belonged to the ''Metabolism'' and ''Signal Transduction'' GO termsindicates that this factor might modulate genes involved with these GO terms and thereby underlie the crosstalk between inflammation and cell cycle phenomena on the one end and general cellular metabolism on the other.  Other than these, the less well characterized but emerging to be possessing significant binding elements in our dataset include FOXD3, FOXJ2 and HFH-3. These belong to the forkhead box proteins (FOX proteins) that were initially known as regulators of embryonic development and have now been recognized to play important roles in regulating the expression of genes involved in cell growth, proliferation and differentiation. FOXD3 has recently been shown to be involved in the maintenance of neural crest and pluripotent cells [68,69]. The expression of FOXJ2 starts early in life and continues till adulthood. It binds to two different DNA binding sites [70] and our depiction of its enhanced binding to its element in the presence of TNFa suggests that it might be critical in the TNFa action in HepG2 cells. The hepatocyte nuclear factor-3 (HNF-3) fork head homolog (HFH) proteins share homology to the winged helix DNA binding domain and have been implicated in cell development and organogenesis [71]. Study of the role of these factors in the context of various cellular processes, alterations of which, underlie the pathogenesis of several diseases is an area worthy of investigation.
One of the most significant networks that emerged from the TNFa induced altered genes was that of lipid metabolism and this subsequently correlated well with the top canonical pathway, the biosynthesis of steroids. TNFa increases hepatic cholesterol synthesis [72] and interferes with lipid homeostasis [73]. Modulation of SQLE and FDPS has been shown to alter intracellular hepatic cholesterol synthesis [74][75][76] implying their significance in this pathway. Also, SREBP-1 has been demonstrated to regulate the expression of these genes in HepG2 cells and the liver [77,78] that also additionally validates our transcription factor predictions and analyses. In our study, TNFa significantly elevated intracellular cholesterol levels in HepG2 cells which validates our claim of this pathway being the most favored and in these series of events, SQLE is the most significant mediator of the effects of TNFa on cholesterol synthesis. TNFa induces hepatic steatosis in mice [40] with increased hepatic lipid synthesis [72,79]. Such an effect of TNFa on hepatic lipid homeostasis was also reported recently by Qin et al. [80] where they demonstrated that in vivo infusion of TNFa significantly increased hepatic production of total circulatory apolipoprotein B100 and VLDL apolipoprotein B100 in fasting and post prandial states in a hamster model.
To conclude, our results demonstrate that TNFa induced alteration in the hepatocyte transcriptome majorly includes genes of lipid metabolism that acutely targets biosynthesis of steroids and cholesterol. While several characterized transcription factors were consistently found to be having binding elements in these sets of genes, others mainly those of the fork head family that are not too well characterized also appeared to be substantially over represented. Our results presented here put forth an overall detailed analysis of the action of TNFa in HepG2 cells by encompassing the altered gene expression profile, the different GO functional terms these altered genes categorise in and the consequent predicted signatures of conserved transcription factor binding sites in these sets of genes. Given that TNFa is now invariably being associated with several diseases, our results provide information towards identification of the underlying genes and transcription factors that determine among other significant liver phenotypes, hepatic steatosis and fatty liver that are critically associated with several pathological states.   A.HepG2 cells were incubated in the absence (C) or presence of TNFa (0.5 nM, 12 h) and cells were lysed and lipids were extracted by addition of 1.0 ml of chloroform: methanol (2:1). B. HepG2 cells were transfected with 50 nM of either control (Con siRNA) or SQLE or FDPS or EBP siRNAs. After 48 h, they were incubated either alone or in the presence of TNFa as in ''A''. Cholesterol was estimated in the extracted lipid fraction as described in the Methods section. Results are expressed after normalization to the total protein content. Each value is the mean6 SEM of three experiments. *p,0.05 as compared to control (C), determined using Student's t test (A). *p,0.001 as compared to ConsiRNA; **p,0.01 and ***p,0.05 as compared to the incubation of TNFa alone in the presence of ConsiRNA (B). doi:10.1371/journal.pone.0009063.g007