IGFBP-2 as a biomarker in NAFLD improves hepatic steatosis: an integrated bioinformatics and experimental study.

BACKGROUND AND AIMS: Non-alcoholic fatty liver disease (NAFLD) has become a common chronic liver disease in the world. Simple steatosis (SS) is the early phase of NAFLD. However, the molecular mechanisms underlying the development of steatosis have not yet been fully elucidated.
METHODS: Two public datasets (GSE48452 and GSE89632) through the Gene Expression Omnibus (GEO) database were used to identify differentially expressed genes (DEGs) in the development of steatosis. A total of 72 participants including 38 normal histological controls and 34 SS patients were included in this study. Gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and protein-protein interaction (PPI) network analysis were performed to explore the function of DEGs. The results were further confirmed in high-fat diet (HFD)-fed mice and oleate-treated HepG2 cells.
RESULTS: Total 57 DEGs including 31 up- and 26 down-regulated genes between SS patients and healthy controls were determined. GO and KEGG analysis showed that most of the DEGs were enriched in the ligand-receptor signaling pathways. PPI network construction was used to identify the hub genes of the DEGs. MYC, ANXA2, GDF15, AGTR1, NAMPT, LEPR, IGFBP-2, IL1RN, MMP7, and APLNR were identified as hub genes, and IGFBP-2 expression was found to be reversely associated with hepatic steatosis, fasting insulin, HOMA-IR index, and ALT levels. In HFD-fed mice, hepatic IGFBP-2 was also downregulated and negatively associated with hepatic triglyceride (TG) levels. Moreover, overexpression of IGFBP-2 ameliorated the oleate induced accumulation of TGs in hepatocytes.
CONCLUSIONS: This study identified novel gene signatures in the hepatic steatosis and will provide new understanding and molecular clues of hepatic steatosis.

Introduction

Non-alcoholic fatty liver disease (NAFLD), which is considered to be the hepatic manifestation of the metabolic syndrome is defined as the abnormal hepatic triglyceride (TG) accumulation without excessive alcohol use. NAFLD disease spectrum ranges from simple steatosis (SS) to non-alcoholic steatohepatitis (NASH) with or without fibrosis that can further progress to cirrhosis, and hepatocellular carcinoma (HCC). It is estimated that about 25% of the population in the world has some stage of NAFLD (1). Moreover, NAFLD is not only a liver disease but also increases the risk of other metabolic diseases such as diabetes, cardiovascular disease, and chronic kidney disease (2). Thus, it is important to investigate the mechanism of NAFLD occurrence for NAFLD prevention and treatment.

Although the pathogenesis of NAFLD is still unclear, the 'two-hit hypothesis' has been proposed for pathophysiology of NAFLD (3). The first hit is the elevated hepatic lipid accumulation caused by insulin resistance. On the basis of the first hit, liver is more sensitive to the second hit such as oxidative stress and inflammation, and the second hit can further lead to NAFLD progression from SS to NASH. Abnormal lipid accumulation plays the central role in the pathogenesis of NAFLD (4). Excessive lipid accumulation in the liver is caused by dysregulation of lipid metabolism (5). However, the mechanisms of lipid metabolism dysregulation in the NAFLD are still needed to be further clarified.

In the past few years, liver transcriptome technology including microarray and RNA-seq has been widely used for gene expression profiling in the liver tissues from NAFLD patients. However, inconsistencies exist in these studies and some datasets can be further mined to identify hub genes that can be targeted to treat hepatic steatosis. In this study, we combined two public datasets which have full histological data to identify differentially expressed genes (DEGs) between SS and healthy controls (HC). Based on the result of DEGs, gene enrichment and pathway annotation analysis were conducted. Next, protein–protein interaction (PPI) networks and hub gene identification were performed, and IGFBP-2, one of the hub genes was found to be negatively correlated with steatosis grade, fasting insulin levels, HOMA-IR index, and ALT levels. In high-fat diet (HFD)-induced SS mice models, we also found that hepatic IGFBP-2 was negatively associated with hepatic TG levels. Moreover, overexpression of IGFBP-2 in hepatocytes ameliorated the oleate induced accumulation of TGs.

Materials and methods

Datasets of NAFLD

We conducted a systematic search within the GEO database (http://www.ncbi.nlm.nih.gov/geo/) (6) by using the terms: 'Homosapiens', 'steatosis', and 'NAFLD'. Datasets meeting the following eligibility criteria were included in this study: (1) mRNA expression profiling was performed using liver tissues; (2) datasets contained SS subject and matched HCs; and (3) datasets contained steatosis grade.

Dataset analyses

Two datasets: GSE48452 (7) and GSE89632 (8) were included in this study from 28 studies initially identified by screening GEO. In GSE48452 dataset, 14 normal controls and 14 SS patients were included in this study. In GSE89632, 24 normal controls and 20 SS patients were included. The gene expression matrix was downloaded from GEO database. The data were then annotated by using corresponding annotation document to map the probe to the gene symbols. If multiple probes mapped to the one gene, the mean value was used. The DEGs were identified between SS tissues and normal control tissues by using 'limma' R package (9). Genes with |log2 fold change (FC)| > 0.5 and P-value < 0.05 were considered as DEGs. NASH patients and obese patients in these two datasets were excluded in the analysis. Detailes of included datasets information are shown in Supplementary Table 1 (see section on supplementary materials given at the end of this article).

Functional and pathway enrichment analysis

To understand the biological function and pathways of the DEGs, gene ontology (GO) functional enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed using 'ClusterProfiler' R package (10). P-value <0.05 was regarded as the cut-off criteria.

The protein–protein interaction network construction and hub genes identification

PPI networks was constructed by STRING online software to analyze interaction. Confidence >0.15 was used as the selection threshold. Cytoscape software was used to visualize PPI network (11). Maximal clique centrality (MCC) method was used by cytohubba, a plug-in of cytoscape, to determine the centrality of each gene in the PPI network (12). Genes with top MCC values were regarded as hub genes related with steatosis.

Animal models

All mice were housed in a specific pathogen-free facility at the Sun Yat-sen University under a 12 h darkness:12 h light cycle with access to food and water ad libitum. All studies were approved by the Animal Care and Protection Committee of Public Health School, Sun Yat-sen University. Eight-week-old male C57BL/6J mice were fed with HFD (60% fat kcal%) (MD12032, medicience, Jiangsu, China) for 0, 4 or 8 weeks (n =5). Liver samples were obtained, fixed in 4% formalin for 24 h and dehydrated by a series of ethanol solution. After that, liver tissues were embedded in paraffin and sectioned at 4 μm thickness. The sections were used for hematoxylin and eosin (H & E) staining.

Overexpression of IGFBP-2

The cDNA of IGFBP-2 was amplified from HepG2 cells by PCR and then cloned into the pCDH-CMV vector. Lentiviruses were produced by co-transfecting specified pCDH-CMV-IGFBP-2, pMD2.G, and psPAX2 into 293T cells. Viruses supplemented with 10 μg/mL polybrene was transferred into a well seeded with 1 × 105 HepG2 cells.

Quantitative RT-PCR analysis

Trizol reagent (Invitrogen) was used to extract total RNA from liver tissues and cells. And RT was performed with a PrimeScript RT Reagent Kit (TaKaRa) according to the manufacturer’s instructions. The mRNA expression was quantified in a 96-well plate using a 7500 Real-Time PCR System (Applied Biosystems, Inc.) with SYBR Green Master Mix (TaKaRa). The relative mRNA expression levels were determined using the 2-ΔΔCt method with GAPDH as an internal control for normalization. The primers sequences used in this study are shown in Supplementary Table 2.

Results

Identification of DEGs in nonalcoholic simple steatosis patients

As shown in Fig. 1A, we observed 143 up- and 81 down-regulated mRNA expression from the GSE48452 dataset in SS patients compared with HC. And we also found 1755 up- and 2055 down-regulated mRNA expression from the GSE89632 dataset in SS patients compared to HC (Fig. 1B). Moreover, venn diagrams display the overlap of 31 up- and 26 down-regulated genes between SS and HC using GSE48452 and GSE89632 (Fig. 1C and D).

Figure 1

Identification of differentially expressed genes (DEGs) in the development of steatosis. (A) Volcano plot of GSE48452. (B) Volcano plot of GSE89632. (C) Venn diagram displaying 31 up-regulated genes between healthy controls (HCs) and simple steatosis (SS) subjects using GSE48452 and GSE89632. (D) Venn diagram displaying 26 down-regulated genes between HC and SS subjects using GSE48452 and GSE89632.

Functional annotation

GO analysis including biological process, molecular function (MF) and cellular component (CC) were performed to elucidate the biological function of the 57 DEGs between SS and HC (Fig. 2A, B and C). In terms of biological process (BP) ontology, type I interferon signaling pathway was the most significantly enriched (GO: 0060337, P = 7.54 × 10−6). The most significantly enriched GO term of MF was receptor ligand activity (GO: 0048018, P = 6.69 × 10−4). And lipid droplet was the most significantly enriched GO term of CC (GO: 0005811, P = 1.46 × 10−4). Furthermore, KEGG pathways analysis was undertaken to investigate the signaling pathways of the 57 DEGs between SS and HC. As shown in Fig. 2D, the top 3 enriched pathways were PPAR signaling pathway (hsa03320, P = 2.68 × 10−3), Cytokine–cytokine receptor interaction (hsa04060, P = 4.17 × 10−3) and Apelin signaling pathway (hsa04371, P = 1.38 × 10−2).

Figure 2

Functional annotation of DEGs in the development of steatosis. (A) The enriched biological process (BP) terms in gene ontology (GO) analysis. (B) The enriched molecular function (MF) terms in GO analysis. (C) The enriched cellular component (CC) terms in GO analysis. (D) Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis.

PPI network analysis and hub genes identification

STRING database was used to construct PPI network involving the 57 DEGs and Cytoscape software was used to visualize the network as shown in Fig. 3A. The genes in the central of the network were regarded as the hub genes which may have an important role in the pathophysiological process. In the analysis of PPI network by the MCC method, the genes with top 10 MCC value were MYC, ANXA2, GDF15, AGTR1, NAMPT, LEPR, IGFBP-2, IL1RN, MMP7, and APLNR. The hub genes relationship in the STRING database is shown in Fig. 3B.

Figure 3

Construction of protein–protein interaction (PPI) networks for DEGs in the development of steatosis. (A) The PPI network of DEGs in steatosis. (B) Visualization for hub genes.

Correlation of hub genes expression with hepatic steatosis

Hepatic steatosis grade and percentage data were extracted from GSE48452 and GSE89632. Linear regression was performed to examine the association between expression of hub genes and hepatic steatosis. In terms of steatosis grade, all hub genes were significantly correlated with steatosis grade in GSE89632 dataset (Table 1). MYC, ANXA2, GDF15, AGTR1, NAMPT, LEPR, IGFBP-2, IL1RN, and MMP7 expression were negatively correlated with steatosis grade. APLNR expression was positively associated with steatosis grade in the GSE48452 dataset. As regards to steatosis percentage, the hepatic expression of MYC, GDF15, AGTR1, NAMPT, LEPR, and IGFBP-2 were significantly negatively associated with hepatic steatosis percentage (Table 2). Furthermore, IGFBP-2 had the highest absolute value of correlation coefficient between steatosis percentage and the hub genes expression. The correlation between IGFBP-2 expression and hepatic steatosis is shown in Fig. 4.

Table 1

The correlation of hub genes expression with steatosis grade.

Hub genes	Steatosis grade (GSE48452)		Steatosis grade (GSE89632)
	r	P	r	P
MYC	−0.405	0.032	−0.812	<0.001
ANXA2	−0.405	0.033	−0.467	0.001
GDF15	−0.244	0.210	−0.475	0.001
AGTR1	−0.459	0.014	−0.444	0.003
NAMPT	−0.484	0.009	−0.702	<0.001
LEPR	−0.468	0.012	−0.587	<0.001
IGFBP-2	−0.459	0.014	−0.753	<0.001
IL1RN	−0.375	0.049	−0.659	<0.001
MMP7	−0.245	0.208	−0.426	0.004
APLNR	0.374	0.050	0.423	0.004

Table 2

The correlation of hub genes expression with steatosis percentage.

Hub genes	Steatosis percentage (GSE48452)		Steatosis percentage (GSE89632)
	r	P	r	P
MYC	−0.390	0.041	−0.478	0.001
ANXA2	−0.421	0.026	−0.252	0.099
GDF15	−0.264	0.175	−0.401	0.001
AGTR1	−0.382	0.045	−0.186	0.227
NAMPT	−0.441	0.019	−0.394	0.008
LEPR	−0.401	0.035	−0.357	0.017
IGFBP-2	−0.433	0.021	−0.479	0.001
IL1RN	−0.393	0.039	−0.346	0.021
MMP7	−0.182	0.355	−0.204	0.185
APLNR	0.308	0.110	0.200	0.193

Figure 4

Correlation of hepatic IGFBP-2 expression levels with steatosis. (A) The correlation of hepatic IGFBP-2 expression levels with steatosis grade in GSE48452. (B) The correlation of hepatic IGFBP-2 expression levels with steatosis grade in GSE89632. (C) The correlation of hepatic IGFBP-2 expression levels with steatosis percentage in GSE48452. (D) The correlation of hepatic IGFBP-2 expression levels with steatosis percentage in GSE89632.

Correlation of hepatic IGFBP-2 expression with anthropometric and biochemical variables

Next, the clinical data were derived from GSE89632 dataset. Spearman correlation analysis was conducted to analyze the association between hepatic IGFBP-2 expression with anthropometric and biochemical variables. Hepatic IGFBP-2 expression was significantly associated with waist, serum aspartate aminotransferase (AST), serum alanine aminotransferase (ALT), fasting insulin levels, and HOMA-IR index (r = −0.515, r = −0.445, r = −0.626, r = −0.514, and r = −0.518, respectively) (Table 3). After adjusted for age, gender, and fasting insulin levels, hepatic IGFBP-2 expression was still significantly negatively associated with serum ALT and hepatic steatosis percentage (r = −0.453 and r = −0.446, respectively). The correlation of hepatic IGFBP-2 expression with serum AST, serum ALT, fasting insulin levels, and HOMA-IR index is shown in Fig. 5.

Table 3

The correlation of hepatic IGFBP-2 expression with biochemical parameters (GSE89632).

Variables	IGFBP-2 expression		IGFBP-2 expression (age and gender adjusted)		IGFBP-2 expression (age, gender, and insulin adjusted)
	r	P	r	P	r	P
Gender	−0.049	0.753	–	–	–	–
Age (years)	0.019	0.911	–	–	–	–
Fasting insulin (pmol/L)	−0.514	0.001	−0.623	0.001	–	–
BMI (kg/m2)	−0.287	0.062	−0.306	0.146	0.221	0.311
Waist (cm)	−0.515	0.001	−0.332	0.112	0.043	0.846
AST (U/L)	−0.445	0.003	−0.378	0.069	−0.242	0.265
ALT (U/L)	−0.626	<0.001	−0.620	0.001	−0.453	0.030
ALP (U/L)	0.007	0.963	0.105	0.626	−0.051	0.816
Triglycerides (mmol/L)	−0.205	0.217	−0.389	0.060	−0.072	0.744
Total cholesterol (mmol/L)	−0.183	0.270	−0.316	0.133	−0.225	0.302
LDL cholesterol (mmol/L)	−0.193	0.273	−0.361	0.083	−0.267	0.219
HDL cholesterol (mmol/L)	0.322	0.059	0.289	0.171	0.050	0.821
Fasting glucose (mmol/L)	−0.241	0.139	−0.173	0.419	−0.227	0.297
HOMA-IR	−0.518	0.001	−0.655	0.001	−0.306	0.156
Hemoglobin a1c	−0.253	0.120	−0.045	0.833	0.066	0.763
Steatosis percentage	−0.491	0.001	−0.691	<0.001	−0.446	0.033

Figure 5

Correlation of hepatic IGFBP-2 expression levels with biochemical variables. (A) The correlation of hepatic IGFBP-2 expression levels with AST levels. (B) The correlation of hepatic IGFBP-2 expression levels with ALT levels. (C) The correlation of hepatic IGFBP-2 expression levels with fasting insulin levels. (D) The correlation of hepatic IGFBP-2 expression levels with HOMA-IR index.

Expression of hub genes in HFD mice

To further confirm the expression of hub genes in the mice model, C57BL/6J mice were fed with HFD for 0, 4, and 8 weeks. H & E staining showed that HFD induced hepatic SS in mice without obvious inflammation (Fig. 6A). Moreover, hepatic TG levels and serum fasting glucose were significantly increased in mice after 8 weeks of HFD feeding (Fig. 6B and C). And compared to the controls, HFD-fed mice exhibited significantly higher expression of APLNR and lower expression of NAMPT, LEPR, AGTR1, and IGFBP-2 in the liver (Fig. 6D). Correlation analysis showed that the IGFBP-2 expression was significantly negatively associated with hepatic TG levels (Fig. 6E).

Figure 6

Hub genes expression in high-fat diet (HFD)-fed mice. Eight-week-old male C57BL/6J mice were fed with HFD for 0, 4, or 8 weeks (n = 5). (A) Hematoxylin and eosin (H & E) staining of liver slices. (B) Hepatic triglyceride (TG) levels. (C) Fasting glucose levels. (D) Hepatic hub genes mRNA expression. (E) The correlation of IGFBP-2 expression levels with hepatic TG levels. Values are shown as the mean ± s.d. *P < 0.05 vs HFD 0 week.

Overexpression of IGFBP-2 ameliorated OA-induced accumulation of triglycerides in HepG2 cells

Lentivirus particles containing pCDH-CMV-IGFBP-2 (Lv-IGFBP-2) transfection induced significantly higher expression of IGFBP-2 in HepG2 cells (Fig. 7A). Furthermore, Oil red O staining showed that IGFBP-2 overexpression by Lv-IGFBP-2 decreased the lipid accumulation in OA-treated cells (Fig. 7B). Cellular TG levels were also significantly decreased by IGFBP-2 overexpression (Fig. 7C). Furthermore, IGFBP-2 overexpression decreased the expression of lipogenesis related genes (SREBP1, FAS, and ACC) and increased the expression of lipid oxidation related genes (PPARα and CPT1) (Fig. 7D).

Figure 7

Overexpression of IGFBP-2 ameliorated steatosis in HepG2 cells. HepG2 cells preloaded with oleate (200 μM) were treated with Lv-IGFBP-2. (A) The mRNA expression of IGFBP-2. (B) Oil red O staining of cells. (C) Cellular TG levels. (D) The mRNA expression of lipid metabolism related genes. Values are shown as the mean ± s.d. *P < 0.05 vs OA-treated cells.

Discussion

Despite the increasing prevalence of NAFLD all over the world, the molecular mechanism regulating the steatosis and initiation of NAFLD is still unclear. In this study, we combined two different GSE datasets which have histological data of NAFLD and identified a series of DEGs between controls and SS patients. Functional annotation and PPI network were performed to investigate the roles of DEGs. Hub genes were identified by the PPI network analysis and were significantly correlated with steatosis. And IGFBP-2, one of the hub genes was found to be negatively correlated with serum fasting insulin, HOMA-IR index and ALT levels. Moreover, overexpression of IGFBP-2 ameliorated accumulation of TGs in the HepG2 cells.

IGFBP-2 is a 36 kDa protein which belongs to one of six similar proteins that bind to insulin-like growth factors I and II (IGF-I and IGF-II). IGFBP-2 has an RGD integrin-recognition motif and a heparin-binding domain (13). Studies have demonstrated that IGFBP-2 may be involved in the metabolic disease. A prospective study showed that IGFBP-2 levels were strongly and inversely associated with the risk of type 2 diabetes mellitus (T2DM) (14). Moreover, serum IGFBP-2 concentrations were inversely correlated with BMI and insulin resistance (15, 16). Overexpression of IGFBP-2 in mice showed resistance to obesity and improved insulin resistance (17, 18), and IGFBP-2 was able to prevent adipogenesis and stimulate glucose uptake in adipocytes (17, 19). In accordance with the above findings, we found that hepatic IGFBP-2 expression was reversely associated with fasting insulin levels and HOMA-IR index. However, whether the protective effect of IGFBP-2 against hepatic steatosis is dependant on the its beneficial effects on glucose metabolism is still unclear. Future studies are needed to better clarify this issue. The role of IGFBP-2 in the NAFLD development is also partly reported. Adenoviral overexpression of IGFBP-2 has been shown to reverse hepatic steatosis in obese mice and hepatic expression of IGFBP-2 was epigenetically inhibited in patients with hepatic steatosis (7, 18). Recently, one cohort study demonstrated that circulating IGFBP-2 levels were negatively associated with NAFLD incidence (20). Pia et al. found that IGFBP-2 is related to hepatic de novo lipogenesis and may serve as a biomarker of NAFLD progression (21). Here, we combined two datasets from GSE database and found that hepatic IGFBP-2 expression as a hub gene was also reversely associated with steatosis and serum ALT/AST levels. Moreover, downregulation of IGFBP-2 was confirmed in the HFD-induced steatosis mice model. Overexpression of IGFBP-2 attenuated OA-induced TGs accumulation in HepG2 cells. Taken together, these results suggested that IGFBP-2 might be a promising target for SS in NAFLD progression.

In all, 57 DEGs were found between SS and HC with 31 up- and 26 down-regulated genes. In the functional annotation of DEGs, the GO and KEGG analysis showed that most of DEGs were enriched in the ligand–receptor signaling pathways including type I interferon, cytokine, PPAR, hormone activity, and apelin signaling. Nowadays, researches have delineated numerous ligand–receptor signaling among the cells in the liver which play a critical role in the liver homeostasis and injury response including NAFLD (22, 23, 24). However, the detailed construction and deep understanding of ligand–receptor signaling especially in the SS state of NAFLD need to be further investigated. The genes with top 10 MCC value in the PPI analysis were treated as hub genes including MYC, ANXA2, GDF15, AGTR1, NAMPT, LEPR, IGFBP-2, IL1RN, MMP7, and APLNR. Among the hub genes, MYC, ANXA2, GDF15, AGTR1, NAMPT, LEPR, IGFBP-2, IL1RN, and MMP7 were downregulated in SS. And APLNR expression was upregulated in the SS patients.

MYC is a proto-oncogene which plays an important role in cell cycle, apoptosis, and cellular transformation. Studies showed that MYC expression and activation promote NASH and hepatocarcinogenesis (25, 26). However, we found that hepatic expression of MYC was lower in SS and one study demonstrated that MYC deficiency can induce abnormal lipid accumulation in the hepatocytes (27). It seems that MYC may play multiple roles in the development of NAFLD.

All, ANXA2, GDF15, NAMPT, and MMP7 can function as secreted proteins in the biological process. ANXA2 as an autocrine factor was found to stimulate osteoclasts formation (28), and the role of ANXA2 in the NAFLD is still unknown. GDF15 which is the ligand of TGF-beta receptors acts as a pleiotropic cytokine involving in many diseases (29). Overexpression of GDF15 ameliorated the liver injury and NASH progression (30, 31), and we also found hepatic GDF15 was lower in the SS patients. NAMPT have two major functions. The enzymatic function is to participate one step in the biosynthesis of NAD(32). On the other hand, extracellular NAMPT has cytokine-like activity (32). Lower NAMPT expression has been found in steatosis and was able to aggravate steatosis which is consistent with our findings (33, 34). MMP7 which is involved in the breakdown of extracellular matrix was found to be positively associated with fibrosis in the NAFLD patients (35) but the role of MMP7 in the steatosis is unidentified.

AGTR1, LEPR, IL1RN, and APLNR are all receptors which are expressed by liver cells. AGTR1 is the type 1 receptor of angiotensin II. Studies have demonstrated AGTR1 variant and polymorphisms were correlated with the risk of NAFLD (36, 37). Lower hepatic AGTR1 expression was found in SS patients in our analysis. As for LEPR, the leptin receptor, lower expression of LEPR and LEPR mutation will lead to leptin resistance which promotes the metabolic disease progression including NAFLD (38, 39). Here, we found hepatic LEPR expression was decreased in the SS patients. IL1RN as the inhibitor of interleukin 1, was found to be lower in SS patients. However, the role of IL1RN in the steatosis needs to be further verified. Apelin–APLNR signaling is involved in many metabolic disorders (40). Plasma apelin levels were higher in the NAFLD and apelin protected against steatosis (41, 42). In this study, we found higher hepatic expression of APLNR in SS as compared to HC. Future studies are needed to elucidate whether the higher expression of APLNR is compensatory or dentrimental in the progression of NAFLD.

In summary, we have investigated the molecular changes in hepatic SS by combining two public datasets from GSE. We also successfully identified several candidate therapeutic targets including MYC, ANXA2, GDF15, AGTR1, NAMPT, LEPR, IGFBP-2, IL1RN, MMP7, and APLNR as hub genes. However, more experiments are needed to elucidate the functions of these hub genes in the development of steatosis, and more large and independent liver-biopsy NAFLD populations are needed to further confirm these results in the future.

Declaration of interest

The authors declare that there is no conflict of interest that could be perceived as prejudicing the impartiality of the research reported.

Funding

This study was supported and funded by the National Key R&D Program of China (SQ2019YFC1604603), National Natural Science Foundation of China (Grant No. 82073583).

Institutional review board statement

All studies were approved by the Animal Care and Protection Committee of Sun Yat-sen University. All animal procedures were conducted in accordance with the National Institutes of Health guide for the care and use of Laboratory animals (NIH Publications No. 8023, revised 1978).

Author contribution statement

Xu Chen: Methodology, Investigation, Funding acquisition, and Writing – Original Draft; Yi Tang: Formal analysis, Data Curation, and Investigation; Shen Chen: Investigation; Wenhua Ling: Project administration; Qing Wang: Conceptualization, Supervision, Funding acquisition, and Writing – Review and Editing.