Gene expression profiling reveals key genes and pathways related to the development of non-alcoholic fatty liver disease

Abstract

Background. This study aims to identify key genes and pathways involved in non-alcoholic fatty liver disease (NAFLD).

Material and methods. The dataset GSE48452 was downloaded from Gene Expression Omnibus, including 14 control liver samples, 27 healthy obese samples, 14 steatosis samples and 18 nonalcoholic steatohepatitis (NASH) samples. Differentially expressed genes (DEGs) between controls and other samples were screened through LIMMA package. Then pathway enrichment analysis for DEGs was performed by using DAVID, and alterations of enriched pathways were determined. Furthermore, protein-protein interaction (PPI) networks were constructed based on the PPI information from HPRD database, and then, networks were visualized through Cytoscape. Additionally, interactions between microRNAs (miRNAs) and pathways were analyzed via Fisher’s exact test.

Results. A total of 505, 814 and 783 DEGs were identified for healthy obese, steatosis and NASH samples in comparison with controls, respectively. DEGs were enriched in ribosome (RPL36A, RPL14, etc.), ubiquitin mediated proteolysis (UBE2A, UBA7, etc.), focal adhesion (PRKCA, EGFR, CDC42, VEGFA, etc.), FcγR-mediated phagocytosis (PRKCA, CDC42, etc.), and so on. The 27 enriched pathways gradually deviated from baseline (namely, controls) along with the changes of obese-steatosis-NASH. In PPI networks, PRKCA interacted with EGFR and CDC42. Besides, hsa-miR-330-3p and hsa-miR-126 modulated focal adhesion through targeting VEGFAand CDC42.

Conclusions. The identified DEGs (PRKCA, EGFR, CDC42, VEGFA), disturbed pathways (ribosome, ubi-quitin mediated proteolysis, focal adhesion, FcγR-mediated phagocytosis, etc.) and miRNAs (hsa-miR-330-3p, hsa-miR-126, etc.) might be closely related to NAFLD progression. These results might contribute to understanding NAFLD mechanism, conducting experimental researches, and designing clinical practices.

Introduction

Non-alcoholic fatty liver disease (NAFLD) is the most common cause of liver disease worldwide with a prevalence of 15–30% in Western populations.^1–4 In patients with NAFLD, nearly one quarter is nonalcoholic steatohepatitis (NASH),⁵ which is defined as the presence of hepatic steatosis and inflammation with hepatocyte injury (ballooning) with or without fibrosis.⁶ Many conditions, including obesity and location of fat storage, play a role in the pathogenesis of NAFLD, especially NASH.⁷ However, the molecular mechanism underlying NAFLD is still vague.

In the past years, NAFLD-related genes and microR-NAs (miRNAs) have been discovered. Previous studies have shown elevated levels of C-reactive protein (CRP), interleukin-6 (IL-6) and vascular endothelial growth factor (VEGF) in NASH patients compared with controls.⁸,⁹ The genetic variant I148M (rs738409) in Patatin-like phospholipase domain-containing 3 (PNPLA3) affects the secretion of hepatic very low density lipoproteins (VLDL) and promotes intracellular lipid accumulation in the liver by reducing the lipidation of VLDL, facilitating NAFLD development.¹⁰ Besides, miRNA-10b has been proven to regulate steatosis level by targeting peroxisome proliferator-activated receptors (PPAR)-α expression in NAFLD.¹¹ It’s reported that apoptosis-related p53 overexpression activated miR-34a/SIRT1/p53 pathway in NAFLD morbid obese patients.¹²

Along with the development of high-throughput technology, microarray analysis and next-generation sequencing have been widely utilized to identify candidate disease-related genes and pathways. Based on RNA sequencing analysis of Shp(-/-) mice, novel gene signatures implicated in NASH have been identified and validated, including tetraspanin 4 (TSPAN4), peptidoglycan recognition protein 2 (PGLYRP2), thrombospondin 1 (THBS1), dual specific phosphatase-4 (DUSP4), etc.¹³ Through gene expression profiling and gene set enrichment analysis, Yoneda, et al. have found genes involved in NASH, including 27 up-regulated gene sets and 25 down-regulated gene sets.¹⁴ By analyzing gene expression profiles in high- and low-risk NAFLD patients, Moylan, et al. has found a 64-gene profile that can differentiate severe NAFLD from mild NAFLD and an independent 20-gene subset that correlates with NAFLD severity.¹⁵ Through array-based DNA methylation and mRNA expression profiling analyses of liver samples from healthy obese patients, steatosis patients, NASH patients, and controls, Ahrens, et al.¹⁶ have discovered NAFLD-specific mRNA expression and DNA methylation of nine genes coding for key enzymes in intermediate metabolism (including PC, ACLY and PLCG1) and insulin/insulin-like signaling (including IGF1, IGFBP2 and PRKCE). It’s also found that NAFLD-associated methylation changes were partially reversible. However, the pathways changed in NAFLD progression, and the miRNAs that regulate these altered pathways are still unclear.

In this study, the microarray dataset GSE48452 up-loaded by Ahrens, et al.¹⁶ was further analyzed to screen differentially expressed genes (DEGs) between healthy obese samples, steatosis samples, NASH samples, and controls. Then, pathway enrichment analysis was performed for DEGs, and variation of enriched pathways at different statuses was detected. Furthermore, protein-protein interaction (PPI) networks were constructed. Additionally, miRNAs targeting DEGs enriched in pathways were predicted. These results might be conducive to further experimental studies and provide better understanding of the nosogenesis of NAFLD.

Material and Methods

Gene expression data

The gene expression data set GSE48452 was downloaded from Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/), including 14 control liver samples (group C), 27 healthy obese samples (group H), 14 steatosis samples (group S) and 18 NASH samples (group N). It should be noted that only gene expression data were downloaded, and thus the sample size and information of enrolled individuals are different from the previous article of Ahrens, et al. In addition, information on sample number, age, body mass index, sex, fat area, inflammation, fibrosis, and non-alcoholic fatty liver disease activity score were collected from GEO (Table 1). Gene expression levels were measured by using [HuGene-1_1-st] Affymetrix Human Gene 1.1 ST Array [transcript (gene) version] (Affymetrix Inc., Santa Clara, California, USA). For the gene expression experiment, all patients provided written and informed consent, and its protocol was approved by the institutional review board (“Ethikkommission der Medizinischen Fakultät der Universität Kiel,” D425/07, A111/99).¹⁶

Table 1.

Overview for normal controls, healthy obese samples, steatosis samples, and NASH samples.

	Normal controls (C)	Healthy obese (H)	Steatosis (S)	NASH (N)
Number	14	27	14	18
Age	52 (23-77)	45 (33-60)	42 (24-65)	45 (32-58)
BMI	25 (17-31)	41 (29-56)	48 (40-60)	49 (24-70)
Sex (% male)	36	7	29	22
Fat (area in %)	1 (0-3)	2 (0-7)	36 (10-70)	72 (40-90)
Inflammation (0-3)	0 (0-1)	0 (0-1)	0 (0-1)	2 (1-3)
Fibrosis (0-4)	0 (0-2)	0 (0-1)	0 (0-1)	1 (0-4)
NAS (0-8)	0 (0-1)	0 (0-2)	2 (1-3)	5 (4-7)

The median and the interquartile range are provided for all numeric parameters. NASH: nonalcoholic steatohepatitis. BMI: body mass index. NAS: non-alcoholic fatty liver disease activity score.

DEGs screening

Linear Models for Microarray Data (LIMMA) package¹⁷ of Bioconductor (available at http:// www.bioconductor.org/packages/release/bioc/html/) was used to identify genes that were differentially expressed between healthy obese and control samples (H vs. C), between steatosis and control samples (S vs. C), and between NASH and control samples (N vs. C). Because the DEGs with |log₂ fold change (FC)| > 1 were rare, |log₂ FC| ≠ 0 and p-value < 0.01 were chosen as the cut-off criteria.

Pathway enrichment analysis

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enriched by the three sets of DEGs were revealed by using the Database for Annotation, Visualization and Integrated Discovery (DAVID, available at http:// david.abcc.ncifcrf.gov/).¹⁸ The p-value < 0.1 was set as the cut-off criterion.

Scoring altered pathways

Three sets of significantly enriched pathways (H vs. C, S vs. C, and N vs. C) were combined, and a total of N pathways were obtained. The deviations of every pathway at status H, status S, and status N from status C were calculated based on gene expression levels.¹⁹ Firstly, DEGs from the gene pool of a specific pathway were selected out. Secondly, p-value was converted to Z-score (Z) by using inverse cumulative distribution function. Finally, an alteration score of a pathway, namely A (P), was determined for pathway P based on DEG Z-scores and following algorithm:

Several issues should be taken into account. First, not all DEGs involved in a pathway had significant influence on this pathway. For example, changes of downstream genes might be induced by abnormalities of upstream genes. In other words, changes of upstream genes involved in pathway could have more significant influence on pathways. Second, number of genes varies in different pathways. Thus, random perturbation treatment was applied to eliminate the influence of pathway size. During calculation, Z-score of DEGs was sorted in descending order. Bigger Z-score meant higher degree of differential expression. Assuming that pathway P contained k DEGs, when t genes (1 < t < k) had maximum average Z-score, these t genes were considered to have the most significant contribution to the pathway P, and then, A (P) was also calculated. In order to eliminate the influence of pathway size, normalization of each A (P) to A_corrected(P) score was accomplished by correcting background distributions based on the mean and standard deviation which were derived from A(P) scores of 10,000 random sets of k genes.

Construction of protein-protein interaction networks

PPIs were obtained from Human Protein Reference Database (HPRD, available at http://www.hprd.org/)²⁰ via Python program.²¹ If a protein normally expressed in human interacted with at least three DEGs, this protein would be selected to construct PPI networks together with DEGs. Then, PPI networks were visualized by using Cytoscape (available at http://cytoscape.org/).²² Moreover, in order to disclose changes of the networks at different statuses, topological properties of these networks were analyzed by using network analysis (a plugin in Cytoscape), such as degree, average shortest path length (ASPL), clustering coefficient (CC) and topological coefficient (TC).

Prediction of miRNAs associated with significant pathways

The common regulatory relationships between miRNAs and DEGs were abstracted from three databases, including miRecord,²³ miRTarBase²⁴ and Tarbase 6.0.²⁵ A total of 5,489 interactions were obtained, containing 482 miRNAs and 2,331 target genes. Then, Fisher’s exact test²⁶ was adopted to identify pairs of miRNA and pathway with significant relevance. Fisher’s exact test examines whether the proportion of target genes (namely, target genes of miRNA) in the gene group is same to the proportion of target genes in the entire genome. There are two hypotheses: first, whether the gene belongs to target gene group or not; second, whether the gene belongs to functional gene group (for example, group of pathway genes). Fisher’s score was defined as the probability of at least x genes can be identified in functional gene group in the K target genes.

N represents the number of total genes in the whole genome; M represents the number of genes in a functional genes group; and K represents the number of target genes. Finally, altered pathways at different statuses were identified, as well as the miRNAs regulating the corresponding pathways.

Results

Identification of DEGs

According to the threshold (p-value < 0.01 and |log₂ FC| ≠ 0), a total of 505 DEGs (including 205 up-regulated and 300 down-regulated DEGs), 814 DEGs (including 283 up-regulated and 531 down-regulated DEGs), and 783 DEGs (including 392 up-regulated and 391 down-regulated DEGs) were identified between healthy obese and control samples (H vs. C), between steatosis and control samples (S vs. C), and between NASH and control samples (N vs. C), respectively.

Pathway enrichment analysis for DEGs

In H vs. C contrast group, DEGs were significantly enriched in 6 KEGG pathways, such as FcyR-mediated phagocytosis (PRKCA, CDC42, etc.), ubiquitin mediated proteolysis (UBE2N, UBE2B, etc.), and so on. In S vs. C contrast group, DEGs were markedly enriched in ribosome (RPL36A, RPL14, etc.), ubiquitin mediated proteolysis (UBE2A, UBA7, etc.), focal adhesion (PRKCA, EGFR, CDC42, VEGFA, etc.), and so on. Meanwhile, in N vs. C contrast group, DEGs were significantly enriched in ribosome (RPL36A, RPL15, etc.), FcyR-mediated phagocytosis (PRKCE, PRKCD, etc.), and so on (Table 2).

Table 2.

Pathway enrichment analysis result for the DEGs between healthy obese, steatosis, NASH, and controls.

Group (term)	Count	p-value	DEGs
Obese
hsa00410: beta-alanine metabolism	5	0.001	CNDP1, EHHADH, ALDH2, HIBCH, DPYD
hsa00310: lysine degradation	5	0.017	EHHADH, WHSC1L1, ALDH2, ACAT1, BBOX1
hsa05120: epithelial cell signaling in pylori infection	6	0.010	EGFR, CDC42, IL8, ATP6V0D2, ATP6V0B, ATP6V1F
hsa00640: propanoate metabolism	4	0.036	EHHADH, ALDH2, HIBCH, ACAT1
hsa05110: Vibrio cholerae infection	5	0.039	PRKCA, ATP6V0D2, ATP6V0B, SEC61G, ATP6V1F
hsa05130: pathogenic Escherichia coli infection	5	0.041	PRKCA, CDC42, TUBA8, KRT18, ARPC4
hsa00380: tryptophan metabolism	4	0.063	EHHADH, ALDH2, CAT, ACAT1
hsa04666: FcyR-mediated phagocytosis	6	0.066	PRKCA, DNM3, CDC42, WASF3, ARPC4, CRK
hsa05219: Bladder cancer	4	0.071	EGFR, CCND1, IL8, MYC
hsa00280: valine, leucine and isoleucine degradation	4	0.079	EHHADH, ALDH2, HIBCH, ACAT1
hsa00510: N-glycan biosynthesis	4	0.088	MGAT5B, DPM2, DPAGT1, ALG14
hsa04120: ubiquitin mediated proteolysis	7	0.094	UBE2N, CBLC, UBE2I, FBXO4, ITCH,
			UBE2L3, UBE2B
Steatosis
hsa03010: ribosome	27	0.001	RPL17, RPL36A, RPL14, RPL15, RPL27A,
			RPS27L, RPS25, RPS27, RPL30, MRPL13.
hsa04120: ubiquitin mediated proteolysis	14	0.004	UBE2A, UBA7, SKP1, UBE2B, RBX1, UBE2N,
			CUL2, UBE2D2, CUL7, UBA2.
hsa03050: proteasome	7	0.012	PSMA1, PSMC6, PSMD12, PSMA6, PSMA4,
			POMP, PSMD6
hsa05211: renal cell carcinoma	8	0.026	CDC42, CUL2, PAK2, VEGFA, RAP1B,
			TCEB1, CRK, RBX1
hsa04510: focal adhesion	15	0.041	PRKCA, PARVG, EGFR, FLT4, MYL12B,
			MYL12A, PPP1CB, KDR, CDC42, VEGFA
hsa04115: p53 signaling pathway	7	0.063	CCNB3, ZMAT3, TSC2, MDM2, RRM2B,
			CCNG1, GADD45A
hsa04144: endocytosis	13	0.083	EGFR, PLD2, RET, TSG101, ADRBK1,
			KDR, CDC42.
hsa00563: (GPI)-anchor biosynthesis	4	0.084	PIGG, PIGX, PIGW, PIGA
hsa03018: RNA degradation	6	0.088	WDR61, EXOSC4, LSM6, LSM3, CNOT7, LSM1
NASH
hsa03010: ribosome	18	0.001	RPL36A, RPL15, RPS8, RPS7, RPL18A, RPS3A,
			RPL22, RPL9, FAU, RPS13.
hsa04610: complement and coagulation cascades	11	0.001	C8B, FGG, C9, MASP2, C6, C5, TFPI,
			SERPINC1, C4BPB, C4BPA...
hsa00140: steroid hormone biosynthesis	7	0.018	AKR1C2, HSD17B2, CYP7A1, SRD5A1,
			SRD5A2, UGT2A3, AKR1C1
hsa04666: FcyR-mediated phagocytosis	10	0.030	PLCG1, LIMK1, CFL2, PIP5K1C, ARPC5,
			INPP5D, VAV2, PRKCE, PRKCD, SYK
hsa04930: type II diabetes mellitus	6	0.064	SOCS2, PKLR, IKBKB, PRKCE, INSR, PRKCD
hsa04210: apoptosis	5	0.075	GTF2E2, TAF7, GTF2A2, TAF9B, GTF2B
hsa05020: prion diseases	5	0.075	NCAM2, C8B, C9, C6, C5
hsa00910: nitrogen metabolism	4	0.087	CTH, CA14, CA12, HAL
hsa03320: PPAR signaling pathway	7	0.097	ME1, LPL, PLIN1, CYP7A1, SCD, FADS2, FABP4

DEGs: differentially expressed genes. NASH: nonalcoholic steatohepatitis.

Alterations of pathways

After combining the above three sets of pathways enriched by DEGs, a total of 27 significant pathways were obtained. Taking normal controls as the baseline (Figure 1), alteration score was determined for each pathway. On the whole, all of the 27 pathways were gradually deviated from the baseline along with the change of normal-obesity-steatosis-NASH. Some pathways showed significant deviation from the baseline and obvious changes along with the change of normal-obesity-steatosis-NASH, such as ribosome, ubiquitin mediated proteolysis, focal adhesion and FcR-mediated phagocytosis.

Figure 1.

The line chart displaying the deviation of the 27 pathways at three statuses. Blue line: healthy obese; red line: steatosis; green line: NASH. The 27 biological pathways are in X-axis and alteration score of the pathway is in Y-axis. NASH: nonalcoholic steatohepatitis.

PPI networks analysis

The PPI network for DEGs in H vs. C contrast group consisted of 254 genes and 551 interactions (Figure 2), and some genes had high degrees, such as EGFR (degree = 49) and PRKCA (degree = 40) (Table 3). Besides, in the net-work for DEGs in S vs. C contrast group, there were 464 genes and 1,059 interactions (Figure 3), and some genes had high degrees, such as EGFR (degree = 62), PRKCA (degree = 47) (Table 3). Meanwhile, in these two networks, PRKCA interacted with numerous genes, such as EGFR and CDC42. Furthermore, the network for DEGs in N vs. C contrast group contained 471 genes and 995 interactions (Figure 4), and some genes had higher degrees, such as PRKCD (degree = 51) and PLCG1 (degree = 48) (Table 3). Additionally, along with the change of obesity-steatosis-NASH, the average degrees of PPI network decreased, ASPLs of PPI network increased, CCs and TCs of PPI network exhibited small declines (Table 4).

Figure 2.

The protein-protein interaction network for the differentially expressed genes between healthy obese samples and normal controls.

Table 3.

The top 10 genes in the PPI networks of obese, steatosis, and NASH

Group	Gene	Degree	Gene	Degree
Obese	EGFR	49	CDC42	22
	PRKCA	40	KRT18	19
	PTK2B	29	EIF2AK2	18
	CRK	26	MYC	18
	UBE2I	24	ATXN1	17
Steatosis	EGFR	62	PTK2B	29
	PRKCA	47	KDR	27
	CRK	36	CDC42	26
	MDM2	31	PPP2CA	23
	HDAC1	30	PAK2	23
NASH	PRKCD	51	HDAC1	23
	PLCG1	48	PPP2CA	23
	SYK	41	MET	20
	INSR	34	KRT18	20
	KDR	27	INPP5D	20

PPI: protein-protein interaction. NASH: nonalcoholic steatohepatitis.

Figure 3.

The protein-protein interaction network for the differentially expressed genes between steatosis samples and norma controls.

Figure 4.

The protein-protein interaction network for the differentially expressed genes between NASH samples and normal controls. NASH: nonalcoholic steatohepatitis.

Table 4.

Topological properties of the 3 networks.

Network features	Obese	Steatosis	NASH
Degree	4.338	4.267	3.869
CC	0.059	0.052	0.048
ASPL	3.881	4.251	4.419
TC	0.265	0.256	0.254

CC: clustering coefficient ASPL: average shortest path length. TC: topological coefficient NASH: nonalcoholic steatohepatitis. The value is a mean for all genes in each network.

miRNAs related to significant pathways

Some miRNAs were identified to regulate DEGs enriched in a set of significant pathways. For instance, hsa-miR-126 and hsa-miR-330-3p regulated focal adhesion via targeting VEGFA and CDC42 (Table 5).

Table 5.

The top 10 miRNA-target-pathway with the highest significant relevance.

miRNA	KEGG pathway	p-value	DEGs
hsa-miR-7	hsa04144: endocytosis	0.001555	EGFR
hsa-miR-34c-5p	hsa05219: bladder cancer	0.001578	MYC
hsa-miR-106b	hsa04120: ubiquitin mediated proteolysis	0.001630	ITCH
hsa-miR-29b	hsa05211: renal cell carcinoma	0.001702	CDC42, VEGFA
hsa-miR-7	hsa05219: bladder cancer	0.001719	EGFR
hsa-miR-126	hsa04510: focal adhesion	0.001878	VEGFA, CRK
hsa-miR-330-3p	hsa04510: focal adhesion	0.001941	VEGFA, CDC42
hsa-miR-15a	hsa05219: bladder cancer	0.002370	CCND1
hsa-miR-29a	hsa04144: endocytosis	0.002834	CDC42, RET
hsa-miR-34a	hsa05219: bladder cancer	0.003064	CCND1

miRNA: microRNA. KEGG: Kyoto Encyclopedia of Genes and Genomes. DEGs: differentially expressed genes. The genes enriched in pathways were regulated by the corresponding miRNAs.

Discussion

In the present study, gene expression profiles of healthy obese, steatosis and NASH samples were compared to healthy controls, and a total of 505, 814 and 784 DEGs were identified, respectively. According to the pathway deviation analysis, obvious changes of pathways enriched by DEGs along with the change of normal-obesity-steatosis-NASH were observed, such as ribosome (RPL36A, RPL14, etc.), ubiquitin mediated proteolysis (UBE2A, UBE2B, etc.), focal adhesion (PRKCA, EGFR, CDC42, VEGFA, etc.), and FcyR-mediated phagocytosis (PRKCA, CDC42, etc.). Reportedly, obesity is a risk factors of NAFLD,⁷ and NASH is always developed from steatosis. Therefore, it’s supposed that obesity develops into steatosis, steatosis develops into NASH, and these developments are associated with the deviation of the 27 pathways.

Among these pathways, ubiquitin mediated proteolysis and ribosome were both associated with protein metabolism. In the liver with non-alcoholic steatosis or steatohepatitis, inactivation of the ubiquitin-proteasome pathway mediates cell death caused by oxidative stress.²⁷ In this study, UBE2A (ubiquitin-conjugating enzyme E2A) and UBA7 (ubiquitin-like modifier activating enzyme 7) were demonstrated to take part in ubiquitin mediated proteolysis. Thus, UBE2A and UBA7 might be involved in NAFLD through ubiquitin mediated proteolysis.

There is no relevant evidence showing ribosome is associated with NASH or NAFLD, while ribosome is closely related to apoptosis,²⁸ which is a pathological feature of NAFLD and promote the development from simple steatosis to NASH.²⁹ In the present study, RPL36A (ribosomal protein L36a) and RPL14 (ribosomal protein L14) were found to participate in ribosome pathway. In addition, the expression levels of RPL36A were significantly down-regulated in healthy obese samples (log₂ FC = −0.186), steatosis samples (log₂ FC = −0.421), and NASH samples (log₂ FC = −0.340), in comparison with normal controls; the expression levels of RPL14 were significantly down-regulated in healthy obese samples (log₂ FC = −0.104), steatosis samples (log₂ FC = −0.221), and NASH samples (log₂ FC = −0.187), in comparison with normal controls. The down-regulation of RPL36A and RPL14 might further promote cell death. Therefore, these genes might play a role in NAFLD through ribosome pathway.

Focal adhesion is crucial in liver disease, because the stage of simple steatosis progressing over steatohepatitis to fibrosing steatohepatitis is characterised by adhesion of leucocytes to the sinusoidal endothelium.³⁰ In this study, PRKCA (protein kinase C α), EGFR (epidermal growth factor receptor), CDC42 (cell division cycle 42) and VEGFA (vascular endothelial growth factor A) were significantly enriched in focal adhesion. The expression levels of PRKCA were significantly up-regulated in healthy obese samples (log₂ FC = 0.210), steatosis samples (log₂ FC = 0.253), and NASH samples (log₂ FC = 0.063), in comparison with normal controls. In the PPI networks for DEGs of H vs. C contrast group and S vs. C contrast group, PRK-CA had high degrees and interacted with EGFR and CDC42. The study of Yang, et al. has shown that EGFs can stimulate proliferation of hepatic stellate cells,³¹ which is the primary effector cell, orchestrating the deposition of extracellular matrix (ECM) in fibrotic liver.³² In this study, the expression levels of EGFR were significantly up-regulated in healthy obese samples (log₂ FC = 0.224), steatosis samples (log₂ FC = 0.243), and NASH samples (log₂ FC = 0.059), in comparison with normal controls. Besides, CDC42 participates in the activation of saturated fatty acid-induced c-Jun N-terminal kinase (JNK) in hepatocytes,³³ which is evident in NASH patients.³⁴,³⁵ In addition, VEGFs play an important role in the development of liver fibrosis and hepatocarcinogenesis in NASH, coordinating with leptin-mediated neovascularization.³⁶ In this study, the expression levels of VEGFA were significantly up-regulated in healthy obese, steatosis, and NASH samples, in comparison with normal controls. This result was consistent with the study (at serum level) of Coulon, et al.³⁷ Hence, PRKCA, EGFR, CDC42, and VEGFA might participate in the development of NAFLD through focal adhesion.

In the progression of NAFLD, phagocytic dysfunction is aggravated.³⁸ In this study, DEGs PRKCA, CDC42, PRKCE and PRKCD were enriched in FcyR-mediated phagocytosis. It should be noted that PRKCA, PRKCE (protein kinase C ε), and PRKCD (protein kinase C δ) belong to the same protein family, and thus share similar bio-functions in liver.³⁹ As discussed above, PRKCA, CDC42, PRKCE, and PRKCD might participate in the development of NAFLD through FcyR-mediated phagocytosis.

In our study, it’s also discovered that CDC42 and VEG-FA were regulated by hsa-miR-330-3p, and VEGFA was regulated by hsa-miR-126, in the pathway of focal adhesion. Reportedly, the expression of hsa-miR-330-3p is altered in NAFLD,⁴⁰ and miR-126 is down-regulated in type 2 diabetes,⁴¹ which is commonly associated with NASH.⁴² Also, miR-126 can target insulin receptor substrate-1 (IRS-1),⁴³ which is correlated with the severity of liver damage in patients with NAFLD.⁴⁴ Thereby, hsamiR-330-3p and hsa-miR-126 might play pivotal roles in the development of NAFLD via regulating focal adhesion through VEGFA or CDC42.

In summary, a number of DEGs (including PRKCA, EGFR, CDC42, VEGFA, etc.), disturbed pathways (including ribosome, ubiquitin mediated proteolysis, FcγR-mediated phagocytosis, focal adhesion, etc.) and several miRNAs (including hsa-miR-330-3p, hsa-miR-126, etc.) were revealed in the development of NAFLD. The results of this study might be conducive to the further understanding of NAFLD mechanism, and provide a theoretical basis for further experimental researches and clinical practices. For example, the key genes and miRNAs identified in this study might serve as therapeutic targets for the treatment of NAFLD after generating pre-clinical evidences.

Disclosures of Conflict of Interest

The authors declare that they have no conflict of interests.

Authors’ Contributions

RFW proposed the study. RFW and XW performed research and wrote the first draft. SR collected and analyzed the data. All authors contributed to the design and interpretation of the study and to further drafts. LWZ is the guarantor.

Acknowledgements

This work was financially supported by the Natural Science Foundation of Heilongjiang Province of China (No. D201139).