Comprehensive bioinformatics analysis of critical lncRNAs, mRNAs and miRNAs in non‑alcoholic fatty liver disease.

Non‑alcoholic fatty liver disease (NAFLD) is the most common fatty liver disease in developed countries, in which fat accumulation in the liver is induced by non‑alcoholic factors. The present study was conducted to identify NAFLD‑associated long non‑coding RNAs (lncRNAs), mRNAs and microRNAs (miRNAs). The microarray dataset GSE72756, which included 5 NAFLD liver tissues and 5 controls, was acquired from the Gene Expression Omnibus database. Differentially expressed lncRNAs (DE‑lncRNAs) and mRNAs (DE‑mRNAs) were detected using the pheatmap package. Using the clusterProfiler package and Cytoscape software, enrichment and protein‑protein interaction (PPI) network analyses were conducted to evaluate the DE‑mRNAs. Next, the miRNA‑lncRNA‑mRNA interaction network was visualized using Cytoscape software. Additionally, RP11‑279F6.1 and AC004540.4 expression levels were analyzed by reverse transcription quantitative polymerase chain reaction. There were 318 DE‑lncRNAs and 609 DE‑mRNAs identified in the NAFLD tissues compared with the normal tissues. Jun proto‑oncogene, AP‑1 transcription factor subunit (JUN), which is regulated by AC004540.4 and RP11‑279F6.1, exhibited higher degree compared with other nodes in the PPI network. Furthermore, miR‑409‑3p and miR‑139 (targeting JUN) were predicted as PPI network nodes. In the miRNA‑lncRNA‑mRNA network, miR‑20a and B‑cell lymphoma 2‑like 11 (BCL2L11) were among the top 10 nodes. Additionally, BCL2L11, AC004540.4 and RP11‑279F6.1 were targeted by miR‑20a, miR‑409‑3p and miR‑139 in the miRNA‑lncRNA‑mRNA network, respectively. RP11‑279F6.1 and AC004540.4 expression was markedly enhanced in NAFLD liver tissues. These key RNAs may be involved in the pathogenic mechanisms of NAFLD.

Introduction

As the most frequently diagnosed fatty liver disease in developed countries, non-alcoholic fatty liver disease (NAFLD) occurs when fat is enriched in the liver due to non-alcoholic factors (1,2). The risk factors for NAFLD include metabolic syndromes, for example combined hyperlipidemia, obesity, high blood pressure and type II diabetes mellitus and insulin resistance (3,4). Non-alcoholic steatohepatitis (NASH), which is the most severe type of NAFLD, is considered to be the primary cause of cirrhosis (4). The incidence of NAFLD is 9–36.9% globally (5,6), and ~20% of people in the United States of America (75–100 million people) are affected by the disease (7). Therefore, exploring the pathogenic mechanisms of NAFLD and developing novel treatment protocols are necessary.

Silencing of fatty acid transport protein 5 (FATP5) may reverse NAFLD; therefore, the activity of hepatic FATP5 is considered critical for maintaining fatty acid flux and caloric uptake during high-fat feeding (8). Patatin-like phospholipase domain containing 3 is an NAFLD-associated gene, which is closely associated with metabolic changes in hepatocytes and lipogenesis (9,10). Interleukin-17 is associated with the proinflammatory response and hepatic steatosis in NAFLD and contributes to the steatosis-steatohepatitis transition (11). Several microRNAs (miRNAs), including miR-21, miR-122, miR-451 and miR-34a, are overexpressed in patients with NAFLD; in particular, miR-122 levels are associated with the grades of fatty liver disease and are a potential marker of NAFLD (12). By inhibiting the expression of 3-hydroxy-3-methylglutaryl-co-enzyme A reductase, miR-21 was demonstrated to mediate cholesterol and triglyceride metabolism in an NAFLD model and may be a promising diagnostic and therapeutic marker for the disease (13). miR-34a/Sirtuin 1/tumor protein p53 (p53) signaling, which is correlated with liver cell apoptosis, is inhibited by ursodeoxycholic acid and activated upon the aggravation of illness in NAFLD (14). Despite these studies, the molecular mechanisms of NAFLD are not fully understood.

In 2015, Sun et al (15) investigated the role of long non-coding RNAs (lncRNAs) in NAFLD through microarray data analysis and identified that several differentially expressed lncRNAs (DE-lncRNAs) function in the pathogenesis of NAFLD. However, the molecular regulatory mechanisms in NAFLD have not been explored in detail. Based on the expression profiles deposited by Sun et al (15), the present study additionally identified the DE-lncRNAs and differentially expressed mRNAs (DE-mRNAs) between NAFLD and normal liver tissues. In addition, key mRNAs, lncRNAs and miRNAs involved in NAFLD were also identified through protein-protein interaction (PPI) network, enrichment and miRNA-lncRNA-mRNA interaction network analyses. The expression of key RNAs were detected by reverse transcription quantitative polymerase chain reaction (RT-qPCR).

Materials and methods

Data source

The normalized expression data and annotation data from the GSE72756 dataset were acquired from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/), which was generated using a GPL16956 Agilent-045997 Arraystar human lncRNA microarray V3 (Probe Name Version) platform. The sample set of GSE72756 included 5 NAFLD liver tissues (3 females and 2 males; mean age=38.8 years) and 5 normal liver tissues (3 females and 2 males; mean age=39.2 years). Samples used in this dataset were sourced from patients with NAFLD without other metabolic complications that were hospitalized in The Third Xiangya Hospital (Changsha, China) from March 2014 to November 2014, and NAFLD was confirmed independently by two senior pathologists by pathological examination. Liver tissues (50–100 mg) were isolated from the patients and then rapidly frozen in liquid nitrogen. Sun et al (15) deposited the microarray dataset GSE72756, and the study was approved by the Ethics Committee of the Third Xiangya Hospital of Central South University. Informed consent was obtained from all patients.

Differential expression analysis

Based on the normalized probe expression data, the R package in Linear Models for Microarray Data (16) (http://www.bioconductor.org/packages/2.9/bioc/html/limma.html) was utilized to identify and annotate the differentially expressed probes between NAFLD and normal liver tissues. A false discovery rate (FDR; adjusted P-value) <0.01 and |log fold-change (FC)|>1 were set as thresholds. Using the R package pheatmap (17) (http://cran.r-project.org/web/packages/pheatmap/index.html), a clustering heatmap was drawn for differentially expressed probes. According to the annotation information, the differentially expressed probes were divided into DE-lncRNAs and DE-mRNAs.

Enrichment analysis and PPI network analysis of DE-mRNAs

Using the R package clusterProfiler (18) (http://bioconductor.org/packages/release/bioc/html/clusterProfiler.html), Gene Ontology (GO; http://www.geneontology.org) (19) functional and Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.genome.ad.jp/kegg) (20) pathway enrichment analyses were performed for the DE-mRNAs. The terms with FDR <0.05 were considered to be significant. Based on the Search Tool for the Retrieval of Interacting Genes (http://www.string-db.org/) (21) database, the interaction pairs among the DE-mRNAs were predicted with the threshold of required confidence >0.4. The PPI network was constructed using Cytoscape software (version 3.0.1, http://www.cytoscape.org) (22). Using the CytoNCA plugin (http://apps.cytoscape.org/apps/cytonca) (23) in Cytoscape, Closeness centrality (CC) (24), Degree Centrality (DC) (25), and Betweenness centrality (BC) (26) scores were calculated. The nodes with increased CC, DC and BC scores compared with other nodes were identified as hub nodes (27) in the PPI network.

Construction of the lncRNA-mRNA regulatory network and functional prediction of DE-lncRNAs

Based on the Pearson product-moment correlation coefficient (28), the target genes of upregulated and downregulated lncRNAs were predicted using a threshold of FDR <0.05 and correlation coefficient >0.995. Next, Cytoscape software (22) was utilized to draw lncRNA-mRNA regulatory networks. Additionally, enrichment analysis was conducted for the target genes of each DE-lncRNA using FDR <0.05 as a threshold. Using the R package clusterProfiler (18), the enriched pathways were compared to identify the significant pathways of target genes for each DE-lncRNA. An FDR <0.05 was used as the cut-off criterion.

Construction of miRNA-lncRNA-mRNA interaction network

Combined with the WEB-based gene set analysis toolkit (Webgestalt; http://www.webgestalt.org) (29), the miRNAs targeting DE-mRNAs involved in the PPI network were predicted. A count (number of target genes) ≥4 and FDR <0.05 were set as the thresholds. According to the annotation information of the DE-lncRNAs, their corresponding sequences were extracted from human reference genome hg19 (30). The mature sequences of the predicted miRNAs were downloaded from the miRbase database (http://www.mirbase.org/) (31). Using MiRanda (http://www.microrna.org) (32) and RNAhybrid (http://bibiserv.techfak.uni-bielefeld.de/rnahybrid/) (33) software, lncRNAs containing significant binding sites for the aforementioned miRNAs were predicted. The predicted results of the miRanda and RNAhybrid analyses were compared to obtain the intersecting miRNA-lncRNA pairs. Based on the obtained miRNA-lncRNA and miRNA-mRNA pairs, and the associated lncRNA-mRNA and mRNA-mRNA pairs, an miRNA-lncRNA-mRNA interaction network was constructed and subjected to topological property analysis using Cytoscape software.

Sample information and RT-qPCR

A total of 2 normal liver tissues (2 females; age range: 47–60 years; mean age=53.5 years; obtained from March to June 2017 via surgical resection) and 2 NAFLD liver tissues (2 males; age range: 41–47 years; mean age=44 years; obtained from April to May 2017 via surgical resection) were provided by Nanjing First Hospital (Nanjing, China). The present study was approved by the Ethics Committee of Nanjing First Hospital. Informed consent was obtained from all patients.

Using TRIzol® reagent (Thermo Fisher Scientific, Inc., Waltham, MA, USA), total RNA was extracted from the samples. The PrimeScript RT Master MIX kit (Takara Bio, Inc., Otsu, Japan) was used to synthesize first-strand cDNA. RT-qPCR was performed using Power SYBR Green PCR Master Mix (cat. no., 4367659; Thermo Fisher Scientific, Inc.) on an ABI 7500 FAST real-time PCR system (Applied Biosystems; Thermo Fisher Scientific, Inc.). The thermocycler conditions for qPCR were as follows: Initial denaturation (50°C, 3 min); 40 cycles of denaturation (95°C, 3 min), annealing (95°C, 10 sec), and extension (60°C, 30 sec). The specificity of the primer amplicons was examined by melting curve analysis. The comparative Cq method (34) was employed to quantify target mRNA and miRNA expression. mRNA expression was normalized to that of GAPDH. The primers used in the present study were as follows: GAPDH-forward (F)-5′-TGACAACTTTGGTATCGTGGAAGG-3′; GAPDH-reverse (R)-5′-AGGCAGGGATGATGTTCTGGAGAG-3′; RP11-279F6.1-h-F-5′-CGGACATAGCCAACGCACCT-3′; RP11-279F6.1-R-5′-TTCATACTTCTGCTGCGTCCA-3′; AC004540.4-F-5′-TTCACAACACACTCAAAGCCT-3′; AC004540.4-R-5′-CAACTGCACTCCAAATGGCTA-3′.

Statistical analysis

All data are expressed as the mean ± standard error of the mean, and the differences between the two groups were compared by Student's t-test. Statistical analyses were performed using SPSS 22.0 software (IBM Corp., Armonk, NY, USA), and GraphPad Prism 5 (GraphPad Software, Inc., La Jolla, CA, USA) was used to visualize the results. P<0.05 was considered to indicate a statistically significant difference.

Results

The clustering heatmap for differentially expressed probes is presented in the Fig. 1. There were 318 DE-lncRNAs, including 105 upregulated and 213 downregulated, and 609 DE-mRNAs, including 353 upregulated and 256 downregulated, in the NAFLD liver tissues compared with the normal liver tissues.

Figure 1.

Clustering heatmap for the differentially expressed probes. GSM1860540, GSM1860541, GSM1860543, GSM1860545 and GSM1860548 represent normal liver tissue samples. GSM1860542, GSM1860544, GSM1860546, GSM1860547 and GSM1860549 represent non-alcoholic fatty liver disease liver tissue samples. The blue-yellow-red variation indicates alterations in expression from low to high. The y- and upper x-axes represent the cluster of genes and samples with similar alterations in expression, respectively.

The top 5 GO ‘biological process’ (BP), GO ‘cellular component’ (CC), GO ‘molecular function’ (MF) and KEGG terms enriched in the DE-mRNAs analysis are summarized in Table I. For the upregulated mRNAs, the enriched GO and KEGG terms primarily included ‘carboxylic acid metabolic process’ (GO_BP; FDR=1.24×10−20), ‘cytoplasmic part’ (GO_CC; FDR=8.66×10−12), ‘molecular function’ (GO_MF; FDR=1.73×10−14), and ‘metabolic pathways’ (KEGG pathway; FDR=2.65×10−09). Downregulated mRNAs were primarily enriched in ‘single-organism process’ (GO_BP; FDR=1.52×10−12), ‘extracellular space’ (GO_CC; FDR=5.95×10−11), ‘molecular function’ (GO_MF; FDR=6.89×10−09), and ‘adrenergic signaling in cardiomyocytes’ (KEGG pathway; FDR=4.34×10−2). The PPI network for the DE-mRNAs involved 442 nodes and 1,409 edges, and it was identified that Jun proto-oncogene, AP-1 transcription factor subunit (JUN) interacted with B-cell lymphoma 2 (Bcl-2)-like 11 (BCL2L11), as demonstrated in Fig. 2. The top 15 nodes, including JUN, with the highest BC, CC and DC scores are summarized in Table II.

Table I.

Top 5 functions and pathways enriched in the analysis of differentially expressed mRNAs.

Category	Term	ID	Description	P-value	FDR	Count
Upregulated	GO_BP	GO:0019752	Carboxylic acid metabolic process	2.90×10−24	1.24×10−20	69
	GO_BP	GO:0032787	Monocarboxylic acid metabolic process	7.56×10−24	1.61×10−20	53
	GO_BP	GO:0006082	Organic acid metabolic process	1.14×10−23	1.62×10−20	73
	GO_BP	GO:0043436	Oxoacid metabolic process	2.13×10−23	2.27×10−20	72
	GO_BP	GO:0044281	Small molecule metabolic process	7.36×10−22	6.27×10−19	110
	GO_CC	GO:0044444	Cytoplasmic part	1.98×10−14	8.66×10−12	206
	GO_CC	GO:0044432	Endoplasmic reticulum part	1.50×10−12	3.27×10−10	55
	GO_CC	GO:0005783	Endoplasmic reticulum	3.69×10−12	5.38×10−10	68
	GO_CC	GO:0044421	Extracellular region part	5.42×10−12	5.93×10−10	120
	GO_CC	GO:0005789	Endoplasmic reticulum membrane	1.20×10−11	1.05×10−9	48
	GO_MF	GO:0003674	Molecular function	2.54×10−17	1.73×10−14	324
	GO_MF	GO:0003824	Catalytic activity	1.18×10−15	4.01×10−13	163
	GO_MF	GO:0048037	Cofactor binding	4.10×10−12	9.30×10−10	25
	GO_MF	GO:0016491	Oxidoreductase activity	8.89×10−12	1.51×10−9	41
	GO_MF	GO:0050662	Coenzyme binding	2.34×10−8	3.18×10−6	17
	KEGG Pathway	hsa01100	Metabolic pathways	1.02×10−11	2.65×10−9	72
	KEGG Pathway	hsa00040	Pentose and glucuronate interconversions	2.24×10−5	2.59×10−3	7
	KEGG Pathway	hsa01200	Carbon metabolism	3.94×10−5	2.59×10−3	12
	KEGG Pathway	hsa00830	Retinol metabolism	4.32×10−5	2.59×10−3	9
	KEGG Pathway	hsa05204	Chemical carcinogenesis	4.99×10−5	2.59×10−3	10
Downregulated	GO_BP	GO:0044699	Single-organism process	4.23×10−16	1.52×10−12	197
	GO_BP	GO:0032501	Multicellular organismal process	1.48×10−12	1.85×10−9	128
	GO_BP	GO:0044707	Single-multicellular organism process	1.54×10−12	1.85×10−9	125
	GO_BP	GO:0008150	Biological process	6.60×10−12	5.93×10−9	210
	GO_BP	GO:0051128	Regulation of cellular component organization	1.24×10−11	8.92×10−9	58
	GO_CC	GO:0005615	Extracellular space	1.54×10−13	5.95×10−11	52
	GO_CC	GO:0044449	Contractile fiber part	1.04×10−11	2.02×10−9	19
	GO_CC	GO:0005576	Extracellular region	2.82×10−11	3.64×10−9	103
	GO_CC	GO:0043292	Contractile fiber	4.73×10−11	4.58×10−9	19
	GO_CC	GO:0030017	Sarcomere	1.73×10−10	1.34×10−8	17
	GO_MF	GO:0003674	Molecular function	1.33×10−11	6.89×10−9	213
	GO_MF	GO:0005515	Protein binding	8.48×10−9	2.19×10−6	159
	GO_MF	GO:0008307	Structural constituent of muscle	1.73×10−8	2.99×10−6	8
	GO_MF	GO:0008092	Cytoskeletal protein binding	3.61×10−8	4.67×10−6	29
	GO_MF	GO:0005488	Binding	1.29×10−6	1.34×10−6	187
	Pathway	hsa04261	Adrenergic signaling in cardiomyocytes	2.73×10−4	4.34×10°	8
	Pathway	hsa00512	Mucin type O-glycan biosynthesis	4.43×10−4	4.34×10°	4

GO, Gene Ontology; BP, biological process; CC, cell component; MF, molecular function; KEGG, Kyoto Encyclopedia of Genes and Genomes; FDR, false discovery rate.

Figure 2.

Protein-protein interaction network of differentially expressed mRNAs. Red and green dots represent upregulated and downregulated mRNAs, respectively. Node size represents the degree of the node.

Table II.

Top 15 nodes with increased DC, BC and CC scores in the protein-protein interaction network.

DC		BC		CC
mRNA	Score	mRNA	Score	mRNA	Score
JUN	59	JUN	33383.63	JUN	0.069956
ACACA	52	ACACA	23963.01	FOS	0.069231
FOS	51	CTNNB1	18141.12	CTNNB1	0.069133
CDK1	45	FOS	16604.96	DECR1	0.068756
DECR1	42	DECR1	15936.45	MAPK3	0.068542
CTNNB1	41	CDK1	12244.42	ACACA	0.068531
MAPK3	36	IL8	11466.93	ACTA1	0.068425
ACTA1	32	MAPK3	10861.69	ACTA2	0.068383
ENO3	31	ACTA1	9482.17	CREM	0.068351
CCNB1	31	GART	8161.233	CDK1	0.068277
IL8	29	MBP	7587.735	AR	0.068087
ACTA2	29	AR	6690.792	GART	0.068024
GART	28	ENO3	6230.197	IL8	0.068014
EHHADH	24	DCN	6073.046	CCNB1	0.067804
DCN	23	RTN4IP1	5927.582	ENO3	0.067784

BC, Betweenness centrality; CC, Closeness centrality; DC, Degree Centrality.

Following prediction of the target genes of the upregulated and downregulated lncRNAs, lncRNA-mRNA regulatory networks were constructed. For the upregulated lncRNAs, the lncRNA-mRNA regulatory network contained 182 nodes (including 37 lncRNAs and 145 target genes) and 672 interactions. For the downregulated lncRNAs, the lncRNA-mRNA regulatory network contained 140 nodes (including 47 lncRNAs and 93 target genes) and 450 interactions, among which AC004540.4 and RP11-279F6.1 were identified to target JUN. The top 10 nodes (including AC004540.4 and RP11-279F6.1) in the lncRNA-mRNA regulatory networks are summarized in Table III. Subsequent to enrichment analysis, the enriched pathways for the target genes of each upregulated and downregulated lncRNA were compared to identify significant pathways.

Table III.

Top 10 nodes in the lncRNA-mRNA regulatory networks for the upregulated and downregulated lncRNAs.

Upregulated		Downregulated
Symbol	Degree	Symbol	Degree
AC004540.4	60	XLOC_014103	47
RP11-279F6.1	54	C17orf76-AS1	45
AQP7P1	53	RP11-279F6.3	41
XLOC_007896	52	AK027145	41
AK025288	51	RP11-13L2.4	34
RP1-60O19.1	47	RP11-120K9.2	30
TRHDE-AS1	46	LOC100505806	29
RP11-345M22.2	42	BC073897	26
LOC100422737	38	RP11-471J12.1	21
AK055386	29	RMST	14

lncRNA, long non-coding RNA.

The miRNAs (including miR-409-3p and miR-139, which targeted JUN) targeting the DE-mRNAs involved in the PPI network are summarized in Table IV. Subsequent to prediction of the binding sites between DE-lncRNAs and the miRNAs associated with the miRNA-mRNA pairs, an miRNA-lncRNA interaction network was constructed (involving 26 miRNAs, 111 lncRNAs and 224 interactions). In the miRNA-lncRNA interaction network, miR-409-3p and miR-139 interacted with AC004540.4 and RP11-279F6.1, respectively. The top 10 nodes exhibiting the highest degrees are summarized in Table V. Finally, an miRNA-lncRNA-mRNA interaction network was constructed, which contained 249 nodes, including 36 miRNAs, 95 lncRNAs, 118 mRNAs and 845 interactions (Fig. 3). The top 10 nodes [including miR-20a; solute carrier family 30, member 10, (SLC30A10); and BCL2L11] with the highest degrees in the miRNA-lncRNA-mRNA interaction network are summarized in Table VI. In particular, BCL2L11 was targeted by miR-20a in the interaction network.

Table IV.

miRNAs targeting the differentially expressed mRNAs involved in the protein-protein interaction network.

miRNA	Gene count	P-value
hsa_CAGTATT, miR-200B, miR-200C, miR-429	19	adjP=0.00004
hsa_ACCAAAG, miR-9	16	adjP=0.00220
hsa_AACATTC, miR-409-3P	8	adjP=0.00220
hsa_GGCAGCT, miR-22	10	adjP=0.00220
hsa_ACTGTAG, miR-139	7	adjP=0.00240
hsa_ACATTCC, miR-1, miR-206	11	adjP=0.00240
hsa_AAAGGAT, miR-501	7	adjP=0.00240
hsa_TGTTTAC, miR-30A-5P, miR-30C, miR-30D, miR-30B, miR-30E-5P	16	adjP=0.00240
hsa_GTGCCTT, miR-506	19	adjP=0.00240
hsa_TACTTGA, miR-26A, miR-26B	11	adjP=0.00240
hsa_ACTTTAT, miR-142-5P	11	adjP=0.00240
hsa_GTGCCAT, miR-183	8	adjP=0.00360
hsa_TGGTGCT, miR-29A, miR-29B, miR-29C	14	adjP=0.00670
hsa_ACTGTGA, miR-27A, miR-27B	13	adjP=0.00680
hsa_GCACTTT, miR-17-5P, miR-20A, miR-106A, miR-106B, miR-20B, miR-519D	15	adjP=0.00700
hsa_GTGCAAT, miR-25, miR-32, miR-92, miR-363, miR-367	10	adjP=0.00760
hsa_TCTGATC, miR-383	4	adjP=0.00820
hsa_CCTGTGA, miR-513	6	adjP=0.00920

hsa, Homo sapiens; miR, microRNA; adjP, adjusted P-value.

Table V.

Top 10 nodes with highest degrees in the miRNA-lncRNA interaction network.

miRNA		lncRNA
Symbol	Degree	Symbol	Degree
hsa-miR-92a	32	C17orf76-AS1	9
hsa-miR-367	21	AC079776.2	8
hsa-miR-20b	21	RP4-669L17.4	6
hsa-miR-17-5p	21	XLOC_001223	5
hsa-miR-145	18	XLOC_000638	5
hsa-miR-25	15	RP11-365O16.3	5
hsa-miR-20a	15	RP11-315I20.1	5
hsa-miR-106a	11	RP11-261C10.3	5
hsa-miR-501	10	RP11-249C24.3	5
hsa-miR-106b	9	BC005927	5

hsa, Homo sapiens; miRNA/miR, microRNA; lncRNA, long non-coding RNA.

Figure 3.

miRNA-lncRNA-mRNA interaction network. Orange diamonds, pink triangles and green circles represent lncRNAs, miRNAs and mRNAs, respectively. Node size represents the degree of the node. lncRNA, long non-coding RNA; miRNA, microRNA.

Table VI.

Top 10 nodes with highest degrees in the miRNA-lncRNA-mRNA interaction network.

miRNA		lncRNA		mRNA
Symbol	Degree	Symbol	Degree	Symbol	Degree
hsa-miR-92	42	AQP7P1	15	FOS	25
hsa-miR-20b	36	C17orf76-AS1	14	JUN	24
hsa-miR-17-5p	36	AK055386	8	BCL2L11	20
hsa-miR-367	31	RP11-365O16.3	7	ACACA	19
hsa-miR-20a	30	AC079776.2	7	STK39	17
hsa-miR-106a	26	SLC30A10	6	SLC38A2	15
hsa-miR-25	25	RP4-669L17.4	6	CDK16	15
hsa-miR-106b	24	XLOC_014103	5	GFPT2	14
hsa-miR-506	21	XLOC_000638	5	PTHLH	14
hsa-miR-200c	20	XLOC_001223	5	RPS6KA2	13

Hsa, Homo sapiens; miRNA/miR, microRNA; lncRNA, long non-coding RNA.

RP11-279F6.1 and AC004540.4 expression

The RT-qPCR results revealed that RP11-279F6.1 and AC004540.4 expression levels were markedly enhanced in the liver tissues from patients with NAFLD compared with the control liver samples (Fig. 4; P<0.01).

Figure 4.

Expression levels of (A) RP11-279F6.1 and (B) AC004540.4 in the control and NAFLD liver tissues by reverse transcription quantitative polymerase chain reaction. ***P<0.01 vs. control. NAFLD, non-alcoholic fatty liver disease.

Discussion

In the present study, 318 DE-lncRNAs, including 105 upregulated and 213 downregulated lncRNAs, and 609 DE-mRNAs, including 353 upregulated and 256 downregulated mRNAs, were screened in the NAFLD liver tissues compared with normal liver tissues. In the PPI network for DE-mRNAs, JUN, the targeting gene of BCL2L11, was among the top 15 nodes. JUN was targeted by the lncRNAs RP11-279F6.1 and AC004540.4 and miRNAs miR-409-3p and miR-139. Additionally, miR-409-3p and miR-139 were predicted as the DE-mRNAs involved in the PPI network. Additionally, miR-409-3p and miR-139 were regulated by AC004540.4 and RP11-279F6.1, respectively. In the miRNA-lncRNA-mRNA interaction network, miR-20a, SLC30A10 and BCL2L11 were among the top 10 nodes. RP11-279F6.1 and AC004540.4 expression was markedly increased in the NAFLD patient liver tissues compared with the control liver samples.

JUN, also known as AP1, was identified to be increased in a previous study of NAFLD (35). Phosphorylation of the transcriptional activation domain of AP1 is conducted by JNKs to enhance its activity, thereby accelerating the progression and development of NASH (36,37). JUN is considered to be an oncogene in the liver, and its expression is enhanced in response to inflammatory stimuli, promoting liver tumorigenesis (38,39). JUN serves an important role in hepatitis B virus-associated tumorigenesis by promoting the proliferation of hepatocytes and dysplasia progression, indicating that JUN is a useful treatment target for preventing hepatitis-associated tumorigenesis (40). The results from the present study indicated that JUN is involved in the pathogenesis of NAFLD.

In the present study, JUN interacted with BCL2L11 and was targeted by miR-20a in the miRNA-lncRNA-mRNA interaction network. The serum/plasma level of miR-20a has potential value for detecting hepatitis C virus (HCV) infection, and therefore circulating miR-20a may be useful as a predictive marker in liver fibrosis mediated by HCV (41). Apoptosis is a major cause of hepatocyte elimination in NAFLD, and inhibition of the anti-apoptotic protein Bcl-2 and activation of the pro-apoptotic protein p53 promotes inflammation in NAFLD (42). Overexpression of Bcl-2 results in resistance to reperfusion injury in the ischemic liver by suppressing apoptosis and is associated with increased caspase 3 and cytoplasmic cytochrome c and a deficiency of Bcl-extra large (43). By targeting the anti-apoptotic gene Bcl-2, miR-15b and miR-16 regulate tumor necrosis factor-mediated hepatic apoptosis in the process of acute liver failure (44). The data from the present study suggest that miR-20a serves a role in NAFLD by targeting BCL2L11 and affecting the expression of JUN.

JUN was also regulated by the miRNAs of miR-409-3p and miR-139. In hepatoma HuH-7 cells, miR-409-3p decreased the production of fibrinogen by downregulating fibrinogen beta chain precursor expression (45). A previous study suggested that miR-409-3p may be utilized to detect the progression of NAFLD (46). miR-409-3p was also identified as a biomarker for the therapeutics and diagnosis of a number of heart failure-associated diseases and a risk factor of NAFLD (47). Overexpression of miR-139, which was downregulated in hepatocellular carcinoma (HCC) samples, suppresses the progression and metastasis of HCC by downregulating rho-kinase 2, indicating that miR-139 may be used to predict the outcome of HCC (48). In the present study, miR-409-3p and miR-139 interacted with the lncRNAs of AC004540.4 and RP11-279F6.1, respectively. Furthermore, the expression levels of AC004540.4 and RP11-279F6.1 were markedly enhanced in liver tissues from patients with NAFLD compared with the control liver samples. Therefore, we hypothesized that miR-409-3p, miR-139, AC004540.4 and RP11-279F6.1 co-regulated JUN expression in patients with NAFLD.

Although the present study succeeded in identifying specific key miRNAs and lncRNAs in the development of NAFLD, there were also certain limitations. For example, the analyses were based on the microarray dataset GSE72756 downloaded from the GEO database. However, the stages of NAFLD (fatty liver, steatohepatitis or fibrosis/cirrhosis) were not clearly described in the GEO database. Therefore, the degree of NAFLD was not clear. In addition, the sample size was too small, and the predicted regulatory associations were not validated. Future studies will aim to confirm the predicted regulatory associations using cell line experiments.

In conclusion, 318 DE-lncRNAs and 609 DE-mRNAs were identified in NAFLD liver tissues by bioinformatics analysis. Additionally, specific mRNAs (including JUN and BCL2L11) and miRNAs (including miR-20a, miR-409-3p and miR-139) may serve essential roles in the pathogenesis of NAFLD. The lncRNAs AC004540.4 and RP11-279F6.1 were also implicated in the mechanisms of NAFLD.