Comparison of hepatic gene expression profiles between three mouse models of Nonalcoholic Fatty Liver Disease.

Non-alcoholic fatty liver disease (NAFLD) has become the most common chronic liver disorder worldwide. Murine models of NAFLD have been widely used to explore its pathogenesis. In this study, we performed a systematic evaluation of hepatic genome-wide mRNA expression by RNA-Sequencing using three mouse models of NAFLD: leptin receptor deficient db/db mice, high-fat high-sugar diet (HSHF)-induced obese mice, and dexamethasone (DEX)-induced NAFLD mice. As a result, we found both distinct and common pathways in the regulation of lipid metabolism from transcriptomes of three mouse models. Moreover, only a total of 12 differentially expressed genes (DEGs) were commonly detected among all three mouse groups, indicating very little overlap among all three models. Therefore, our results suggest that NAFLD is a heterogeneous disease with highly variable molecular mechanisms.

Introduction

Nonalcoholic fatty liver disease (NAFLD) has become the leading cause of liver damage and dysfunction worldwide and is strongly associated with metabolic diseases. The spectrum of NAFLD ranges from simple steatosis to nonalcoholic steatohepatitis (NASH) and end-stage liver diseases. The pathogenesis of NAFLD is multi-factorial and triggered by various environmental factors in the context of a genetic predisposition.,

Murine models are important platforms to investigate a deep mechanistic understanding and test therapeutic interventions for human diseases. Most common mouse models of NAFLD widely used are diet-induced models, including high fat diet, high fructose diet. Besides, leptin receptor deficient db/db mice, is also a frequently used obese model, which develop hepatic steatosis., Moreover, stress induced increased glucocorticoid levels also play a critical role in the development of NAFLD. Indeed, dexamethasone (DEX), a class of frequently prescribed drugs for anti-inflammation in patients, is associated with metabolic side-effects, including dyslipidemia and hepatic steatosis. Overall, the cause of NAFLD could be heterogeneous, while current published research is not enough to elucidate its exact mechanisms. Therefore, classifying the cause of NAFLD and its underlying molecular mechanism could help to gain better understanding into its pathogenesis and identify potential therapeutic targets for its treatment.

Transcriptomic analysis provide a technically and theoretically well-developed global method for the comparison of different murine pathway activation, which may provide a reasonable proxy for such analysis., Over the past decade, several microarray profiling studies have been performed to investigate specific aspects of NAFLD biology.,, However, very few studies have ever simultaneously investigated the full spectrum of fatty liver disease, or performed a comparison between different mouse models. In the present study, we performed a systematic transcriptome-wide analysis of three murine models of NAFLD, including a typical high sugar high fat diet induced NAFLD (HSHF), a dexamethasone-induced NAFLD (DEX), as well as leptin-receptor deficient mice (db/db). We analyzed the biological function of dysregulated genes (DEGs) and associated pathways in each mouse model of NAFLD. Based on these results, we aim to provide novel insights into the molecular characteristics and pathways of NAFLD pathogenesis, and potential therapeutic targets for its treatment.

Materials and methods

Animal experiments

All of animal experiments were approved by the Institutional Care and Use Committee at Shanghai Jiao Tong University School of Medicine with the “Guide for the Care and Use of Laboratory Animals” published by the National Institutes of Health. Wild-type male C57BL/6J mice at the age of 8 weeks were purchased from Shanghai Laboratory Animal Company (SLAC). Male db/db mice aged 8 weeks were purchased from the Nanjing Biomedical Research Institute of Nanjing University (Nanjing, China). One day after arrival, they were habituated for 1 week for regular chow. All mice were housed at 21 °C ± 1 °C with 55% ± 10% humidity and a 12-h light/12-h dark cycle. HSHF group of 8 week old male C57BL/6J mice received a high-sucrose/high fat diet (HSHF diet; 45 kcal% fat, 17 kcal% sucrose) for 12wks as described. For DEX treatment, 8-week old male C57BL/6J mice received 2 weeks of dexamethasone (1 mg/kg; Sigma Aldrich, St. Louis, MO) or vehicle (0.25% PEG400, 0.25% Tween20 in ddH2O) treatment via subcutaneous injection.

Sample collection

Mice were anesthetized with sodium-pentobarbital (Nembutal, 80 mg/kg, ip). The abdominal cavity was opened and blood was withdrawn from the inferior vena cava. The liver was then perfused through the portal vein with ice-chilled saline. Liver from the left lobe was harvested and snap frozen in liquid nitrogen.

Liver triacylglycerol assay

For liver triacylglycerol assay, hepatic tissues (~100 mg) were homogenized in 50 mM NaCl. The homogenate (500 μl) was mixed with chloroform/methanol (2:1, 4 ml) and incubated overnight at room temperature with gentle shaking. Homogenates were vortexed and centrifuged for 5 min at 3000 g. The lower lipid phase was collected and concentrated. The lipid pellets were dissolved in 1% Triton X100 in phosphate-buffered saline, and hepatic triacylglycerol content was determined via enzymatic colorimetric methods.

RNA extraction and transcriptome sequencing

Approximately 5 mg of frozen liver powder was completely dissolved in 500 μl of TRIzol reagent at room temperature. Total RNA was extracted with TRIzol (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's protocol. The quality of total RNA was assessed by the Agilent Bioanalyzer 2100 (Santa Clara, CA, USA). RNA libraries were constructed by using rRNA-depleted RNAs with TruSeq Stranded Total RNA Library Prep Kit (Illumina, San Diego, CA, USA) according to the manufacturer's instructions. Libraries were controlled for quality and quantified using the BioAnalyzer 2100 system (Agilent Technologies, Inc., USA). 10 pM libraries were denatured as single-stranded DNA molecules, captured on Illumina flow cells, amplified in situ as clusters and finally sequenced for 150 cycles on Illumina HiSeq Sequencer according to the manufacturer's instructions.

Data processing and analysis

Paired-end reads were harvested from Illumina HiSeq 4000 sequencer, and were quality controlled by Q30. After 3’ adaptor-trimming and low quality reads removing by cutadapt software (v1.9.3), the high quality clean reads were aligned to the reference genome (UCSC MM10) with hisat2 software (v2.0.4). Then, guided by the Ensembl GTF gene annotation file, cuffdiff software (part of cufflinks) was used to get the gene level FPKM as the expression profiles of mRNA, and fold change and P-value were calculated based on FPKM, differentially expressed mRNA were identified. GO (Gene Ontology) and KEGG Pathway (Kyoto Encyclopeida of Genes and Genomes) enrichment analysis were performed based on the differentially expressed mRNAs. FPKM values were calculated to normalize read counts. Based on FPKM of each gene of each sample, the squares of Pearson coefficient r values were calculated to show correlations between samples and reproducibility. FPKM of each gene was averaged in each group, and log2 (FPKM+1) values were used to generate a heatmap, in which groups and genes were clustered, respectively. The Euclidean distances between groups were also calculated based on log2 (FPKM+1) of all genes. Differentially expressed genes were annotated and classified with GO and KEGG terms by using the R clusterprofile package.

Statistical analysis

Data were expressed as mean ± S.E.M. Student's t test was performed to evaluate significant differences between the compared groups, which were defined as P < 0.05.

Results

RNA-Seq data validation

As shown in Table 1, in this study, we included three different mouse models, including dexamethasone treated group (DEX) and high sugar high fat induced NAFLD group (HSHF), as well as leptin receptor-deficient obese group (db/db), with two to four replicates for each group. As expected, we observed a significant elevation of hepatic triglyceride contents in all three mouse models of NAFLD (Fig. S1). The reproducibility and reliability of these data sets are demonstrated in Fig. S2, showing that gene expression reads as fragments per kilobase of transcript per million mapped reads (FPKM) 3 of all samples in three groups were well-correlated. Pearson correlation r values between all samples in three groups are listed in Table S1, all well correlated. The correlation coefficient matrix of all samples shown in Fig. S3 demonstrates clear patterns, confirming that all of the data within each group were well-correlated. Principal correlation analysis (PCA) of all three groups also showed samples within each group shared similarities (Fig. S3). These data demonstrated that all of the RNA-Seq results obtained were reproducible and reliable.

Table 1

Group information of three mouse models of NAFLD used in the study.

Group	Strain	Treatment	Number of samples (n)
			CON vs. Treatment
1	C57BL6/J WT	DEX for 2wks	4 vs. 4
2	C57BL6/J WT	HSHF for 12wks	2 vs. 2
3	C57BL6/J db/db	/	4 vs. 4

Genome-wide gene expression and clustering analysis

We then carried out hierarchical clustering analysis of gene expressions in the livers of all three different groups. Normalized gene counts were detected after batch-effect correction by IQR (Inter-Quartile Range) boxplot and PCA analysis to confirm it's well normalized with no obvious bias in each group (Fig. 1A, B). Normalized gene counts were further calculated and plotted in a heatmap (Fig. 1C), which shows different gene expression patterns among three different groups. To quantify the differentially expressed genes, we used the Cutdiff to compare the expression level of each gene in each group of different treatment vs. WT C57BL/6J mice and plotted log10 (p.adj) versus log2 (-fold change).

Figure 1

Gene expression patterns in the liver of three different mouse models of NAFLD by RNA-Seq analysis. (A) IQR boxplot and (B) PCA analysis shows normalized read counts in each mouse sample after batch-effect correction. (C) Heatmap shows expression patterns between different models after batch-effect correction.

Differential expression analysis of RNA transcripts from the liver of NAFLD models

To identify DEGs in each group, we used a fold change of |logFC| > 1 and a FDR of q = 0.05 as the threshold. Numbers of DEGs detected were summarized in Fig. 2A and Table 2. Among DEX group, 704 genes were significantly differentially expressed (P.adj <0.05). Of those, 93 genes (13.2%) were significantly upregulated compared with genes in control group; 611 genes (86.8%) were downregulated. In the HSHF group, 403 genes were significantly differentially expressed. Of those, 341 genes (84.6%) were significantly upregulated compared with genes in control group; 62 genes (15.4%) were downregulated. In db/db group, 948 DEGs were detected. Of those, 589 genes (62.13%) were significantly upregulated; and up to 359 genes (37.9%) were downregulated (Fig. 2B).

Figure 2

Analysis of differentially expressed genes (DEGs) in the liver of three different mouse models of NAFLD. (A) Distribution of up-regulated and down-regulated DEGs detected in each model. (B) Quantification of DEGs in three models. Bar graph shows the number of up-regulated (green color) and down-regulated genes (blue color) in each group compared with control group.

Table 2

Global similarity of three NAFLD models using adjusted significance thresholds.

Overview of the pair-wise comparison of the transcriptomic analysis of NAFLD models on a single gene (lower triangle of the table, in grey) and pathway-level (upper triangle, in white). The column “regulated genes” provides the total number of differentially expressed genes for murine phenotype at a genome-wide.

We further analyzed the common DEGs between each of the two groups (Fig. 3 and Table 4). 91 DEGs were detected in both db/db mice and HSHF group, while 64 DEGs in both db/db mice with DEX group. However, only 24 DEGs were detected in both DEX group and HSHF group. The average ratio of overlapping DEGs between each mouse models ranged between 6% and 13%, suggesting the heterogeneity of these three NAFLD models (Fig. 3A–C). Of 91 common DEGs detected between HSHF and db/db mice, some of the genes are related to fatty acid transporter, lipoprotein transport and other lipid metabolic process (Pparg, CD36, Apoa4, Elovl5, Abdcd1) (Fig. 3A). Top genes related to lipid metabolism in each group were summarized in Table 3. Only 12 DEGs were detected in all three groups, with five up-regulated genes, which has been shown to be participated in lipid metabolism (Apoa4, Thrsp, Sqle) (Fig. 3D).

Figure 3

Common key DEGs detected in three different mouse models of NAFLD. (A) Common up- and down-regulated genes detected between db/db and HSHF group. (B) Common up- and down-regulated genes detected between db/db and DEX group. (C) Common up- and down-regulated genes detected between DEX and HSHF group. (D) Common DEGs detected in three NAFLD groups.

Table 4

Common key DEGs in each moue model of NAFLD in comparison with each other.

Types of DEGs	Gene Names
HSHF vs. db/db
DEGs upregulated in both HSHF induced NAFLD and db/db mice (n = 91)	Ly6d, Cyp2a22, Themis, Ttc39a, Acnat2, Acss3, Hamp2, Cd36, Rcan2, Lcn2, Cidec, Cyp17a1, Extl1, Osbpl3, Robo1, Tmem28, Smpd3, Eda, Fam131c, Scd1, Tubb2a, Ubd, Gal3st1, Slc35f2, Elovl5, Tceal8, Tmem98, Abcb1a, Dntt, Pparg, Saa2, 2010003K11Rik, Sulf2, Mest, Acot11, Vnn1, Fam83f, S100a10, Lgalsl, Chchd6, Limk1, Mmd2, Thrsp, Gdf15, Anxa2, Aqp8, Cyp2a4, Them6, Tenm3, Slc17a4, Scd2, Apoa4, Rxrg, Fam124a, Mms22l, Bhlhb9, Acyp2, Acaa1b, Tnfrsf19, Mup19, Ifi27l2b, Saa1, Enc1, Myo1d, Ptp4a3, Anxa5, Copz2, Mgst3, S100a11, Sorbs1, Cldn2, Synj2, Lgals1, Adam11, Gsta1,Abcd1, Ppm1h, Mcm5, Wfdc2, Gale, Srxn1, Ermp1, Abcd2, Samd9l, Fabp4, Serpina7, Unc119, Pgd, Tnip1, Maged2, Gas6
DEGs downregulated in both HSHF induced NAFLD and db/db mice (n = 30)	Slc13a2, Cyp1a2, Gadd45g, Ppp1r3g, Mup13, BC089597, Igfbp1, Fam222a, Cyp51, Sqle, Mup15, C8b, Hsd17b6, Mup14, Idi1, Lama3, Cxcl13, Lifr, Susd4, E130012A19Rik, Cabyr, Mup11, Cyp7b1, Mup12, Cadm4, Selenbp2, Gpr110, Mup7, Enho, Serpina1e
DEGs upregulated in HSHF induced NAFLD but downregulated in db/db mice (n = 15)	Pld4, Lsp1, H2-Eb1, Dsg1c, Smoc2, Nrep, Arhgdib, Tnfsf12, Selplg, Coro1a, Ccl6, H2-Ab1, H2-Aa, H2-DMb1, Tlr12
DEGs upregulated in db/db but downregulated in HSHF induced NAFLD mice (n = 8)	Acot6, Raet1d, Raet1e, Cyp2c40, Etnppl, Dct, G6pc, Cyp26a1
HSHF vs. DEX
DEGs upregulated in both HSHF induced NAFLD and DEX-treated mice (n = 9)	Apoa4, Lcn2, Saa1, Gramd1b, Gdf15, Thrsp, Lyve1, Rgs16, Plin3
DEGs downregulated in both HSHF induced NAFLD and DEX-treated mice (n = 15)	Cyp51, Fbf1, Sqle, Sik1, Cabyr, Hsd17b6, Nsdhl, Sc4mol, Ccbl2, Lifr, Ugt2b37, Idi1, Ugt2b38, Cyp2c54, BC089597
DEGs upregulated in HSHF induced NAFLD but downregulated in DEX-treated mice (n = 24)	Aqp8, Srebf1, Serpina7, H2-Ab1, Hamp2, H2-Aa, Cd74, Tmem106a, Nrep, Mpeg1, Krt18, Cldn2, Abcg5, Paqr7, Trp53i11, Abcd1, Acaa1b, Tenm3, Ddx58, Fads2, Abcg8, Sorbs1, Klhdc7a, Cybb
DEGs upregulated in DEX-treated but downregulated in HSHF induced NAFLD mice (n = 5)	Tat, Arrdc2, Igfbp1, Gadd45g, Scara5
DEX vs. db/db
DEGs upregulated in both DEX-treated and db/db mice (n = 24)	Acacb, Acot1, Treh, Mt2, Obp2a, Agpat9, Cyp2b10, Lcn2, Mt1, Ppp1r3c, Gdf15, Apoa4, Saa1 Slc45a3, Cyp4a14, Prob1, Insig2, Chd1l, Mfsd2a, Fmo5, Thrsp, Gldc, 1810011O10Rik, Bhmt
DEGs downregulated in both DEX-treated and db/db mice (n = 40)	Zfp707, Marco, Ccnd1, Dna2, Irgm2, Rapgef4, Cyp51, Lifr, Dbp, Ttc39c, Apol9b, H2-Aa, Mup3, Irf6, Elovl2, Mug2, BC089597, Fbxo21, Cdh1, Slco1a1, Hsd17b6, Rnase4, Sel1l3, Adh6-ps1, Ces1c, Phlda1, Keg1, Ces1b, Nrp1, Nudt7, Sqle, H2-Ab1, Idi1, Mapk15, Sort1, Cabyr, Aox3, Marveld2, Cyp2u1, Nrep
DEGs upregulated in DEX-treated but downregulated in db/db mice (n = 6)	Porcn, Serpine2, Lipin1, Gadd45g, Igfbp1, AA986860
DEGs downregulated in DEX-treated but upregulated in db/db mice (n = 42)	Ccng1, Hnmt, Gtdc1, 5033411D12Rik, Caln1, Tenm3, Gm13152, Acad10, Hacl1, Abcd1, Rtn4, Hamp2, Cd59a, Btbd7, Adra1b, Pls3, Entpd5, Tlr5, Rwdd3, E2f3, Sorbs1, Pla2g15, Fasn, Nat8, Slc25a10, Aldh1a1, Bche, Paqr9, Ebpl, Rdh16, Prkdc, Slc16a7, Aqp8, Nfe2l2, Acly, Lrit1, Cyp2c38, Acaa1b, Serpina7, Mme, Prlr, Cldn2
DEX vs. HSHF vs. db/db
DEGs upregulated in all three groups (n = 5)	Apoa4, Lcn2, Saa1, Gdf15, Thrsp
DEGs downregulated in all three groups (n = 7)	Cyp51, Sqle, Cabyr, Hsd17b6, Lifr, Idi1, BC089597
DEGs upregulated in DEX-treated but downregulated in both db/db and HSHF induced NAFLD mice (n = 2)	Gadd45g, Igfbp1
DEGs downregulated in DEX-treated but upregulated in both db/db and HSHF induced NAFLD mice (n = 8)	Tenm3, Abcd1, Hamp2, Sorbs1, Aqp8, Acaa1b, Serpina7, Cldn2
DEGs upregulated in HSHF-induced NAFLD but upregulated in both db/db and DEX-treated mice (n = 3)	H2-Ab1, H2-Aa, Nrep

Table 3

Key DEGs significantly changed involved in lipid metabolic process in the liver of each mouse model.

Metabolic process	Changes of top genes in each group
	DEX induced NAFLD			HSHF induced NAFLD			db/db
	Gene name	Log2(Fold Change)	FDR	Gene name	Log2(Fold Change)	FDR	Gene name	Log2(Fold Change)	FDR
Lipid	Srebf1	−2.28955	0.00116535	Apoa4	3.95879	0.00256724	Pnpla3	4.94022	0.00057796
	Elovl2	−1.01537	0.0011653	Scd1	2.75418	0.00256724	Scd1	3.04407	0.0005776
	Acnat1	−1.01528	0.019002	Srebf1	2.15722	0.00256724	Elovl6	3.11782	0.00057796
	Sqle	−1.97855	0.00116535	Elovl5	1.89249	0.00256724	Elovl5	2.46632	0.00057796
	Plin5	2.37117	0.00116535	Lpl	1.94893	0.00256724	Cd36	4.45396	0.0005779
	Acly	−1.65391	0.00116535	Smpd3	2.88631	0.0203045	Pparg	2.28283	0.00057796
	Insig1	−2.18166	0.00116535	Cd36	1.79329	0.00256724	Ppargc1a	2.11752	0.03424480
	Insig2	1.96634	0.0011653	Pparg	1.65057	0.00256724	Acly	1.45867	0.00057796
	Fasn	−1.68574	0.00298341	Apoc2	1.41882	0.0104423	Acnat2	4.83398	0.0005779
	Apoa4	3.01287	0.00116535	Cyp17a1	1.42665	0.00256724	Insig2	2.36451	0.00057796
	Fabp1	−1.0047	0.00116535				Fasn	2.07406	0.000577962
	Acot1	1.32154	0.0170247				Apoa4	1.6458	0.00057796
	Lpin1	1.99705	0.00116535				Scd2	1.66453	0.00057796
	Dhcr7	−1.36218	0.00116535				Pck2	−1.86523	0.00057796
	Hmgcr	−1.51934	0.00116535				Hacl1	1.3138	0.00155599
	Pcsk9	−2.11944	0.00116535				Pdk4	3.26674	0.00057796
	Acsl3	−2.48382	0.0011653				Gpam	2.4615	0.00057796
	Acox1	−1.09071	0.0052477				Pdk1	1.47899	0.00057796
	Fads2	−1.20849	0.00116535				Gpat2	−1.23393	0.02466420
	Fads3	−1.71553	0.00116535				Fabp5	−1.53929	0.000577962
	Soat2	−1.00581	0.00593422				Fabp4	1.10654	0.00057796
							Acot1	1.70572	0.00057796
							Lpin1	−2.43886	0.00057796
							Elovl2	−1.84094	0.00057796
							Acsl5	1.03864	0.00057796

GO analysis

To further understand common and distinct pathways involved, GO tool was used to compare DEGs in all three groups. DEGs detected in each mouse model in GO analysis were classified into three groups, named as, cellular component, biological process and molecular function. Fig. 4A–C and Table S2 showed the top 20 pathways involved in each model. As shown in Fig. 4A–C, DEGs in three mouse models were mainly enriched in biological process. In db/db mice model, up-regulated DEGs were mainly enriched in GO terms fatty acid metabolic process, xenobiotic metabolic process, sulfur compound metabolic process, and lipid catabolic process (Fig. 4A). In HSHF group, up-regulated DEGs were mainly enriched in GO terms leukocyte migration, positive regulation of lipid metabolic process, collagen metabolic process, and antigen processing and presentation of exogenous peptide antigen (Fig. 4B). In DEX group, however, most of the DEGs were down-regulated, which were mainly enriched in GO terms fatty acid metabolic process, steroid metabolic process, cellular ketone metabolic process and regulation of small molecule metabolic process (Fig. 4C).

Figure 4

Over-represented pathways of Gene Ontology (GO) analysis and significantly enriched GO terms of DEGs in different models of NAFLD. (A) Over represented pathways of Gene Ontology (GO) analysis and significantly enriched GO terms of DEGs in db/db mice; (B) Over represented pathways of Gene Ontology (GO) analysis and significantly enriched GO terms of DEGs in HSHF induced NAFLD mice; (C) Over represented pathways of Gene Ontology (GO) analysis and significantly enriched GO terms of DEGs in DEX treated mice.

KEGG signaling enrichment analysis

KEGG signaling enrichment analysis was performed to investigate the pathways based on DEGs identified as well. Table S3 show the top 20 pathways involved in each cohort. In DEX group, down-regulated DEGs were mainly mapped to retinol metabolism, bile secretion, chemical carcinogenesis, and steroid hormone biosynthesis. In HSHF group, up-regulated DEGs were mainly overrepresented in PPAR signaling pathway, cholesterol metabolism, and biosynthesis of unsaturated fatty acid; while down-regulated DEGs were mainly mapped to steroid hormone biosynthesis, retinol metabolism, steroid biosynthesis, and chemical carcinogenesis. In db/db group, up-regulated DEGs were mainly mapped to fatty acid metabolic process, sulfur compound metabolic process, nucleoside bisphosphate metabolic process and ribonucleoside bisphosphate metabolic process.

Comprehensive analysis of distinct and similar pathways in three different mouse models of NAFLD

Since comparing DEGs between different murine groups might overlook common biological pathways, we further clustered GO and KEGG pathways detected in three models. Clustering heatmap showed that three models displayed distinct pathway pattern (Fig. 5D and Table 5). Metabolic pathways, fatty acid degradation, fatty acid elongation, histidine metabolism, glycerolipid metabolism, glutathione metabolism, pyruvate metabolism and other pathways related to metabolism were distinctively enriched in db/db group. DEX group displayed distinct metabolic pathways as well, including cellular ketone metabolic process, lipid catabolic process, triglyceride biosynthesis process, neutral lipid biosynthesis process, acylglycerol biosynthesis process, as well as cofactor biosynthesis process. In HSHF group, distinct pathways detected mainly focus on cholesterol biosynthetic process, cholesterol mechanism, fat digestion and absorption, and steroid biosynthesis process. Also, HSHF group showed specific changes in immune related metabolic pathways, such as cytokine–cytokine receptor interaction, positive regulation of adaptive immune response, and leukocyte activation involved in inflammatory response.

Figure 5

Comparison of over-presented functional enrichment analysis between three models performed on shared pathways. (A) Comparison of over-presented functional enrichment analysis between HSHF and db/db group. (B) Comparison of over-presented functional enrichment analysis between HSHF and DEX group. (C) Comparison of over-presented functional enrichment analysis between DEX and db/db group. (D) Selected functional enrichment analysis show distinct pattern in each group. Over-represented functions are shown in the heat map with adjusted P value < 0.05. Red color representing significantly changed pathway, blue color representing insignificant result.

Table 5

Comparison of over-presented functional enrichment analysis performed on the shared pathways between three mouse models of NAFLD.

Pathways	Different Models
GO and KEGG pathway terms	DEX		HSHF		db/db
	Regulation	BH.p value	Regulation	BH.p value	Regulation	BH.p value
GO:0008202 Steroid Metabolic Process	down	1.10783E-11	down	3.9539E-05	up	0.000592725
GO:0010876 Lipid Localization	down	1.20163E-05	up	7.07082E-05	up	5.54001E-05
GO:0006805 Xenobiotic Metabolic Process	down	6.87516E-07	down	0.015039461	up	4.63612E-09
GO:0062012 Regulation Of Small Molecule Metabolic Process	down	1.97349E-06	up	1.42514E-05	up	0.005334651
GO:0044282 Small Molecule Catabolic Process	up	0.003188734	up	0.006204059	up	4.63612E-09
GO:0035634 Response To Stilbenoid	down	0.00084057	up	1.76742E-06	up	2.68169E-06
GO:0042737 Drug Catabolic Process	up	0.003188734	down	0.02905059	up	9.41721E-08
GO:0062013 Positive Regulation Of Small Molecule Metabolic Process	down	0.006111985	up	4.26301E-06	up	0.003422671
GO:1901615 Organic Hydroxy Compound Metabolic Process	down	1.14533E-05	down	0.002631838	up	1.68246E-06
Mmu00980 Metabolism Of Xenobiotics By Cytochrome P450	down	0.00366925	down	0.04096318	up	0.007078388
mmu02010 ABC Transporters	down	0.005205567	up	0.00718615	up	0.013547692
mmu03320 PPAR signaling pathway	/	/	up	4.45039E-05	up	1.96799E-08
mmu01212 Fatty acid metabolism	down	0.020176687	/	/	up	0.000577213
mmu01040 Biosynthesis of unsaturated fatty acids	/	/	up	0.001703037	up	3.42971E-06
mmu00830 Retinol metabolism	down	1.44555E-11	down	9.7534E-08	up	2.35428E-10
mmu00591 Linoleic acid metabolism	/	/	down	0.019186006	up	0.049841292
mmu00100 Steroid biosynthesis	down	0.000107417	down	2.97813E-05	/	/
GO:1990845 Adaptive thermogenesis	/	/	up	9.50564E-05	up	1.467E-05
GO:1901568 Fatty acid derivative metabolic process	down	0.000479109	/	/	up	9.22256E-06
GO:0120162 Positive regulation of cold-induced thermogenesis	/	/	up	0.001272215	up	0.003804817
GO:0009410 Response to xenobiotic stimulus	down	0.000143497	/	/	up	9.50984E-07
GO:0006869 Lipid transport	down	1.37567E-05	/	/	up	0.000688304
GO:0006790 Sulfur compound metabolic process	down	5.15667E-05	/	/	up	6.22691E-15
GO:0006732 Coenzyme metabolic process	down	9.23777E-05	/	/	up	4.63612E-09
GO:0001659 Temperature homeostasis	/	/	up	7.54435E-05	up	0.003811025

GO, Gene Ontology; NAFLD, nonalcoholic fatty liver disease; HSHF, high sucrose high fat; BH, Benjamini and Hochberg.

Venn diagram showed common pathways detected between each group (Fig. 5A–C). Among all three groups, HSHF group and db/db group displayed most shared pathways based on DEGs detected (Fig. 5A). Interestingly, up to 34 common pathways were detected between db/db group and DEX group even with very few common DEGs detected. Total numbers of GO and KEGG pathways detected in each model were summarized in Table 2. To further understand common pathways detected among all three groups, we clustered resulted from GO and KEGG analysis (Fig. 6). With red dot representing down-regulated pathways and blue representing up-regulated pathways, we could easily recognize different patterns of pathways enriched in three groups. Pathway analysis revealed that steroid metabolic process and regulation of small molecule metabolic pathway were common pathways detected in all three groups. Even though common pathways were detected between DEX group with either db/db or HSHF group, most pathways show opposite regulation status.

Figure 6

Functional enrichment pathway analysis show distinct regulation status of common pathways in three different mouse models of NAFLD.

Discussion

NAFLD is a continuous spectrum of diseases characterized by excessive lipid accumulation in hepatocytes. NAFLD progresses from simple liver steatosis to non-alcoholic steatohepatitis and, in more severe cases, to liver fibrosis, cirrhosis, and hepatocellular carcinoma (HCC). Since mouse models are the cornerstone for mechanistic and interventional studies of disease pathogenesis, many animal models of NAFLD have already been developed., In this study, we analyzed three different murine models to evaluate differences in gene expression and pathways between them. We chose commonly used HSHF diet as diet-induced NALFD model, genetic obese db/db mice and DEX treated mouse model as drug-induced model. Through comprehensive analysis from transcriptome data, we mapped DEGs detected in all three groups to GO and KEGG database, and further compared DEGs and pathways in each group with each other. Only 12 DEGs were detected in all three groups, five of them were up-regulated genes including Apoa4, Thrsp, which are known to participate in lipid metabolism. Pathway analysis revealed that metabolic pathways, steroid metabolic process, regulation of small molecule metabolic pathway were common biological functions in all three groups. We further compared results from pathway analysis in a more specific way, and detected a more similar regulation pattern of metabolic pathways between HSHF and db/db group compared to DEX group. Overall, our clustering analysis showed a more similar gene expression pattern between db/db mice with HSHF group.

As known to all already, most commonly used methods to generate NAFLD models are western diet models. Western-diet is a diet high in saturated fat, trans-fat, and sugar. In humans, this type of diet has been shown to induce obesity, metabolic syndrome, NAFLD, and potentially, NASH. In the liver of mice fed with a high-fat diet (HFD), exaggerated free fatty acid levels were detected, inducing hepatic insulin resistance, decreased fatty acid oxidation, and increased de novo lipogenesis in hepatocytes. HSHF diet, however, is a combination of both HFD and high fructose diet. Findings obtained with C57BL/6 mice fed HSHF diet suggested that fructose consumption is necessary for the progression of fat deposition to fibrogenesis in liver., Though weight gain, body fat and liver steatosis were both detected in HFD and HSHF mouse models, hepatic oxidative stress, liver macrophage numbers, fibrogenesis and collagen deposition were only detected to be increased in mice fed the HSHF diet, suggesting altered immune signaling pathway involved apart from lipid metabolism change in the liver of HSHF diet induced NAFLD model., Genetic modified db/db mice are homozygous for the autosomal recessive diabetic gene (db). The db gene encodes a point mutation of the leptin receptor (Ob-Rb), which leads to defective leptin signalling. These mice have persistent hyperphagia, are obese and diabetic, and develop hepatic steatosis. db/db mice do not spontaneously develop inflammation when fed a normal control diet. So, when fed a standard diet without additional hits, these mice are useful models of NAFLD. Glucocorticoids, one of the most used classes of corticosteroids used in everyday clinic, are commonly used at high doses and for prolonged periods in the treatment of a variety of diseases. Among the many side effects of glucocorticoids, most known related to metabolism are increased insulin resistance with disturbances in glucose homeostasis and increased deposition of lipids in the liver, mostly triglycerides.,, Since different murine models of NAFLD focus on different aspect of the disease pathology, the pathway analysis may provide in itself additional insights.

In defining molecular mechanisms of obesity-associated NAFLD, it has been well-established that steatosis occurs when the rate of fatty acid synthesis is greater than the rate of fatty acid oxidation and secretion. For instance, our recent studies showed that upregulation of lipogenic genes or defects in the fatty acid oxidation could induce liver steatosis and contribute to fatty liver., However, lipogenic genes were repressed by GC-GR signaling in the liver, suggesting that GCs-induced hepatosteatosis could be attributed to alternative mechanisms. Indeed, recent studies showed that excess glucocorticoids could promote hepatic lipid deposition through enhancing fatty acid uptake in the liver and repressing fatty acid oxidation. Together, these reports suggest that the molecular mechanisms of different NAFLD models might be distinct, which needs further intensive studies. In the present study, we provided a comprehensive transcriptome analysis of liver tissue between three totally different mouse models of NAFLD. We detected very few common genes related to metabolism in all three models, including Apoa4, Thrsp and Sqle. Apoa4, apolipoprotein A-IV (ApoA-IV), has been well known for its role in lipid metabolism and metabolic regulation in vitro and in vivo. Previous studies have already indicated that high fat diet, steatosis or dexamehasone treatment could increase Apoa4 expression level in the liver. Besides, hepatic ApoA-IV expression could serve as an early diagnostic biomarker in liver fibrosis. Thrsp, known as a nuclear protein in the regulation of lipid metabolism, has been reported to be highly expressed in the liver and adipose tissue. Sqle, also known as squelene monooxygenase, a key rate-enzyme in cholesterol biosynthesis, has been revealed as a top outlier metabolic gene overexpressed in NAFLD-induced hepatocellular carcinoma.

Results from GO and KEGG pathway analysis among three groups showed very few common pathways. We used clustering heatmap and dot map to show that three mouse models displayed distinct pathway pattern. Apart from the expected enrichment of gene sets associated with lipid metabolism confirming changes in fatty acid biology, an enrichment of several gene sets associated with inflammation was observed in HSHF group mainly. As previously reported in other study, inflammation response pathways in NAFLD might be considered as one of the major relevance and may already be active on a subcellular level, even without yet visible change on a cellular or morphological level., This subcellular change might explain the underlying mechanism of prolonged HSHF diet induced NAFLD development and possibly progression to NASH., As for DEX induced NAFLD, previous studies already showed that pathways of lipid deposition stimulated by excess glucocorticoids include stimulation of de novo lipogenesis, increased release of free fatty acids from adipose stores and stimulation of their uptake by the liver, etc., In our pathway analysis, we found DEX up-regulates triglyceride metabolic process, neutral lipid biosynthesis, as well as other lipid metabolic process, while down-regulates fatty acid metabolic process and steroid metabolic process. Thus, our study provided information of gene transcription changes in the liver of different NAFLD models, indicating very few common pathways cross different models.

Several limitations should be pointed out in our study. Firstly, the number size is relatively small. Secondly, we used only pathway analysis, which may possibly overlook certain similarities between different models. Thirdly, our present study is a descriptive study, functional experiments should be performed to further explore the role and mechanisms of target genes in the development of fatty liver disease.

Conclusions

The present study reveals both transcriptional similarities and differences between three different NAFLD models. Our results suggest that NAFLD is a heterogeneous disease with highly variable molecular mechanisms.

Conflict of interests

The authors have no conflict of interests to declare.