Integrated analysis of lncRNA-associated ceRNA network reveals potential biomarkers for the prognosis of hepatitis B virus-related hepatocellular carcinoma

Abstract

Background

There is evidence that abnormal expression of lncRNAs is associated with hepatitis B virus (HBV) infection-induced hepatocellular carcinoma (HCC). However, the mechanisms remain not fully elucidated. The study aimed to identify novel lncRNAs and explore their underlying mechanisms based on the ceRNA hypothesis.

Methods

The RNA and miRNA expression profiling in 20 tumor and matched adjacent tissues from HBV–HCC patients were retrieved from the Gene Expression Omnibus database under accession numbers GSE77509 and GSE76903, respectively. Differentially expressed lncRNAs (DELs), miRNAs (DEMs), and genes (DEGs) were identified using the EdgeR package. Protein–protein interaction (PPI) network was constructed for DEGs followed by module analysis. The ceRNA network was constructed based on interaction relationships between miRNAs and mRNAs/lncRNAs. The functions of DEGs were predicted using DAVID and BinGO databases. The prognosis values (overall survival [OS] and recurrence-free survival [RFS]) of ceRNA network genes were determined using The Cancer Genome Atlas (TCGA) data with Cox regression analysis and Kaplan–Meier method.

Results

The present study screened 643 DELs, 83 DEMs, and 1,187 DEGs. PPI network analysis demonstrated that CDK1 and CCNE1 were hub genes and extracted in functionally related modules. E2F2, CDK1, and CCNE1 were significantly enriched into cell cycle pathway. FAM182B-miR-125b-5p-E2F2 and LINC00346-miR-10a-5p-CDK1/CCNE1 ceRNA axes were obtained by constructing the ceRNA network. Patients with high expressions of DELs and DEGs in the above ceRNA axes had poor OS, while patients with the high expression of DEMs possessed excellent OS. CDK1 was also an RFS-related biomarker, with its high expression predicting poor RFS. The upregulation of LINC00346 and CDK1 but the downregulation of miR-10a-5p in HCC was validated in other microarray datasets and TCGA database.

Conclusion

The LINC00346-miR-10a-5p-CDK1 axis may be an important mechanism for HBV-related HCC, and genes in this ceRNA axis may be potential prognostic biomarkers and therapeutic targets.

Keywords:

• hepatocellular carcinoma
• hepatitis B virus
• ceRNA
• lncRNA
• miRNA
• prognosis
• bioinformatics analysis
• TCGA

Introduction

Hepatocellular carcinoma (HCC) is the fourth most prevalent human malignancy and the third cause of cancer-related deaths in China.¹ In 2015, it is estimated that there are 466,100 new cases and 422,100 deaths due to this disease.¹ Although patients with HCC can be managed with a series of therapeutic methods (including surgical resection, adjuvant chemotherapy, radiotherapy, and liver transplantation), the overall 5-year survival rate still remains poor (less than 20%).² Epidemiological studies have shown that chronic hepatitis B virus (HBV) infection is the predominant risk factor for the development, metastasis, and recurrence of HCC, accounting for about 80% of all HCC in China.^3-5 Thus, it is necessary to further investigate the molecular mechanisms of HBV-related HCC in order to screen novel prognostic biomarkers and develop effective therapeutic strategies.

Recently, there have studies to indicate that the abnormal expression of lncRNAs, a class of noncoding RNAs longer than 200 nt in length, is associated with the development of various cancers, including HBV-related HCC.^6,7 For example, Zuo et al found that lncRNA AX800134 was upregulated in HBV-positive HCC compared with HBV-negative HCC. Silencing AX800134 with siRNA interference significantly suppressed the growth and invasion but enhanced spontaneous apoptosis of HBx-expressing HepG2 cells.⁸ The study of Lv et al revealed that the expression of lncRNA DREH was frequently downregulated in HBV-associated HCC tissues in comparison with adjacent noncancerous hepatic tissues. Inhibition of DREH expression by HBx remarkably promoted the proliferation of HCC cells in vitro and in vivo.⁹ Yang et al identified that lncRNA-HEIH was highly expressed in liver samples from patients with HBV-related HCC. The expression level of lncRNA-HEIH in HBV-related HCC was significantly associated with recurrence and was an independent prognostic factor for survival.¹⁰ However, the mechanisms of lncRNAs in HBV-related HCC remain not fully elucidated.

Previously, emerging evidence has demonstrated that lncRNAs may function as molecular sponges for a miRNA through their miRNAs response elements (MREs) and thereby influence the translation inhibition or mRNA degradation of the transcript on the targets by the respective miRNAs, which is proposed as ceRNA hypothesis.¹¹ Accumulating data also indicated that this regulatory action plays important roles in HCC development.^12,13 For example, Lv et al¹⁴ showed that lnRNA Unigene56159 promoted the migration and invasion of HCC cells by acting as a ceRNA for miR-140-5p to de-repress the expression of Slug and induce the epithelial–mesenchymal transition (EMT). Mo et al¹⁵ observed the upregulated LINC01287 competitively bound to miR-298 and increased the expression of its target gene STAT3 to promote EMT and invasion of HCC cells. lncRNA n335586 was also reported to promote EMT of HCC cells and then migration as well as invasion through facilitating the expression of its host gene creatine kinase, mitochondrial 1A (CKMT1A) by competitively binding miR-924.¹⁶ lncRNA SNHG12 functioned as an oncogene to accelerate tumorigenesis and metastasis of HCC cells by sponging miR-199a/b-5p, which resulted in the high expression of MLK3 (mitogen-activated protein kinase kinase kinase 11) and activated the NF-κB pathway.¹⁷ HCAL directly interacted with and functioned as a sponge for miR-15a, miR-196a, and miR-196b to modulate lysosomal protein transmembrane 4 beta (LAPTM4B) expression in HCC.¹⁸ However, studies performed to investigate the ceRNA mechanisms of lncRNAs for HBV-related HCC were rare.^14,16

The goal of this study was to identify novel lncRNA– miRNA–mRNA interaction axes for explaining the development of HBV-associated HCC by constructing a ceRNA regulatory network using sequencing data collected from a public database. Also, the prognosis performance of related lncRNAs, miRNA, and mRNAs was also validated by utilizing The Cancer Genome Atlas (TCGA) datasets. We believe that our study may provide novel prognostic biomarkers and therapeutic targets for HBV-associated HCC.

Methods

Data collection and preprocessing

Two datasets under accession numbers GSE77509¹⁹ and GSE76903¹⁹ were retrieved from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). These two datasets examined the RNA expression profiling and noncoding (miRNA) expression profiling in 20 primary tumors and 20 matched adjacent normal tissues from patients with HBV-induced HCC by high-throughput sequencing via HiSeq 2500 System (Illumina, San Diego, CA, USA). The samples were the same for the two datasets. The fragment per kilobase per million mapped reads (FPKM) expression data in TXT files were downloaded and preprocessed by removing low abundance genes with an FPKM of <1. The lncRNA and mRNA genes in RNA expression profiling were annotated based on the Ensembl Gene ID and HUGO Gene Nomenclature Committee (HGNC; http://www.genenames.org/).²⁰

Differentially expressed RNA analysis

The differentially expressed genes (DEGs), differentially expressed lncRNAs (DELs), and differentially expressed miRNAs (DEMs) between primary tumors and adjacent normal tissues were identified using the EdgeR package of R software (Version 3.22.3; http://www.bioconductor.org/packages/release/bioc/html/edgeR.html).²¹ P-value was adjusted to false discovery rate (FDR) with multitest package (Version 2.36.0; http://bioconductor.org/packages/release/bioc/html/multtest.html).²² The FDR of <0.05 and |logFC(fold change)| >1 were set as the statistical threshold value. Hierarchical cluster heatmap representing the expression intensity and direction of DEGs, DELs, and DEMs was generated using the pheatmap R package (Version: 1.0.8; https://cran.r-project. org/web/packages/pheatmap) based on Euclidean distance.

Protein–protein interaction (PPI) network

The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING; Version 10.0; http://stringdb.org/) database²³ was used to assess the direct and indirect correlations between DEGs. The screened interaction pairs among DEGs were used to construct the PPI network with the Cytoscape software (Version 3.6.1; www.cytoscape.org/).²⁴ The topological features of the PPI network, consisting of degree (the number of edges [interactions] of a node protein]), betweenness centrality (BC; the number of shortest paths that run through a node), closeness centrality (CC; the average length of the shortest paths to access all other proteins in the network), and average path length (APL; the average of distances between all pairs of nodes), were then calculated using the CytoNCA plugin in cytoscape software (http://apps.cytoscape.org/apps/cytonca)²⁵ to determine which genes were hub nodes. Functionally related clusters with well-interconnected genes were further identified from the PPI network using the Molecular Complex Detection (MCODE; Version: 1.4.2, http://apps.cytoscape.org/apps/mcode) algorithm²⁶ with default scoring options. Modules with score >4 and node >6 were considered to be significant.

Function enrichment analysis

Gene Ontology (GO) Biological Process term and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of DEGs were conducted using the Database for Annotation, Visualization and Integrated Discovery (DAVID) online tool (Version 6.8; http://david.abcc.ncifcrf. gov)²⁷ and BinGO²⁸ plugin in Cytoscape to predict their underlying functions. Statistical significance was defined as FDR <0.05.

lncRNA–miRNA–mRNA ceRNA regulatory network construction

The miRcode database (Version 11; http://www.mircode.org/)²⁹ was used to screen the interaction relationships between DELs and DEMs, and then, the DELs–DEMs interaction network was constructed using the Cytoscape software (Version 3.6.1; www.cytoscape.org/).²⁴ The target genes of DEMs in the DELs–DEMs interaction network were predicted using the miRwalk database (Version 2.0; http://www.zmf.umm. uni-heidelberg.de/apps/zmf/mirwalk2),³⁰ which were then overlapped with the DEGs to obtain the DELs–DEMs–DEGs regulatory relationships. The negative interaction pairs between DEMs and DEGs/DELs were integrated to construct the DELs–DEMs–DEGs ceRNA network using the Cytoscape software (Version 3.6.1; www.cytoscape.org/).²⁴

Furthermore, all known HCC-related pathways were downloaded from Comparative Toxicogenomics Database (CTD; http://ctd.mdibl.org/),³¹ which was then overlapped with the pathways enriched by the genes in the ceRNA network to obtain potentially HCC-related ceRNA network.

Prognosis values of DELs, DEMs, and DEGs in ceRNA network

The miRNAs and mRNAs expression profile data of HCC were also collected from TCGA (https://gdc-portal.nci.nih. gov/) database, with only the HBV samples having survival information, included. Univariate Cox regression analysis was performed to screen prognosis-related DELs, DEMs, and DEGs using the survival package (Version 2.40.1; https://cran.r-project. org/package=survival). The samples were then classified into a low-expression group (< median) and a high-expression group (> median) based on the expression of each prognosis-related DEL, DEM, and DEG. The Kaplan–Meier (KM) method with the log-rank test was employed to compare the overall survival (OS) and recurrence-free survival (RFS) between the high- and low-expression groups through the GraphPad Prism software (Version 5; GraphPad Software, Inc., La Jolla, CA, USA). P<0.05 was considered to be statistically significant.

Validation of expressions of crucial DELs, DEMs, and DEGs in ceRNA network

The expressions of crucial DELs, DEMs, and DEGs were also validated in TCGA dataset and other microarray datasets that detected the mRNA (GSE121248: 70 vs 37; GSE94660³²: 21 vs 21; GSE25599³³: 10 vs 21 normal), miRNA (GSE69580: 5 vs 5), and lncRNA (GSE27462³⁴: 5 vs 5) expression profile between tumor and matched adjacent tissues from HBV– HCC patients. All microarray datasets were also collected from the GEO database. The expression difference was tested by t-test. P<0.05 was set as statistical significance.

Results

Differential expression analysis

A total of 133 lncRNAs, 18,628 protein-coding mRNAs, and 2,578 miRNAs were annotated in mRNA-seq and miRNA-seq data. After removing the low abundance genes with an FPKM of <1, 80 lncRNAs, 16,169 protein-coding mRNAs, and 874 miRNAs were left for differential expression analysis. Based on the given threshold (FDR <0.05 and |logFC| >1), a total of 43 DELs (14 downregulated and 29 upregulated), 83 DEMs (10 downregulated and 73 upregulated), and 1,187 DEGs (650 downregulated and 537 upregulated) were identified between adjacent normal tissues and primary tumors (Table 1). The heat map analysis showed that the samples with similar features tended to be clustered according to the expressions of DELs (Figure 1A), DEMs (Figure 1B), and DEGs (Figure 1C).

Figure 1 Hierarchical clustering and heat map analysis of differentially expressed lncRNAs (A), miRNAs (B), and genes (C).

Notes: Each row represents a sample, and each column represents an lncRNA, miRNA, or gene. High- or low-relative expression is displayed as a red or blue strip, respectively. Each group contained 20 different samples.

Abbreviations: T, tumor; N, normal.

Table 1 Differentially expressed genes, miRNAs, and lncRNAs

Symbol	LogFC	FDR	Symbol	LogFC	FDR	Symbol	LogFC	FDR
CDC6	1.05	2.01E–03	hsa-miR-215-3p	29.41	3.08E–06	LINC01662	4.62	1.21E–07
CCNE1	1.02	1.22E–02	hsa-miR-301b	3.80	8.46E–11	DSCR8	4.49	1.73E–17
CDK1	1.04	2.02E–03	hsa-miR-483-3p	3.42	3.40E–09	LINC01976	4.30	7.43E–06
E2F2	1.03	7.83E–03	hsa-miR-410-3p	3.25	1.76E–08	LINC00632	3.70	8.51E–13
BUB1	1.11	6.22E–04	hsa-miR-7974	3.25	3.05E–08	DSCR4	3.63	2.03E–12
CCNB1	1.03	2.17E–03	hsa-miR-483-5p	3.22	2.62E–08	LINC02089	3.50	9.85E–05
SFN	1.54	4.46E–06	hsa-miR-200c-3p	3.17	4.31E–08	MIR2052HG	3.34	4.32E–11
GNG4	1.48	6.54E–05	hsa-miR-183-5p	3.11	6.10E–08	PRNT	2.54	1.17E–03
UBE2C	1.23	2.15E–04	hsa-miR-1910-5p	3.09	4.44E–04	LINC00346	2.02	3.33E–05
GNGT1	2.04	1.49E–05	hsa-miR-493-5p	3.08	8.66E–08	FAM182B	1.30	1.04E–02
KIF4A	1.15	5.96E–04	hsa-miR-139-5p	−2.22	1.49E–04	PACRG-AS3	−1.24	2.72E–02
GNG13	1.40	1.87E–02	hsa-miR-30c-2-3p	−1.59	1.04E–02	LINC01558	−1.59	1.17E–03
PF4	−1.36	2.03E–02	hsa-miR-378i	−1.54	1.32E–02	LINC02312	−1.60	1.47E–03
IL-6	−1.04	2.94E–02	hsa-miR-199a-5p	−1.45	2.11E–02	LINC02453	−1.74	3.36E–04
CXCR1	−1.056	1.64E–02	hsa-miR-30a-5p	−1.45	2.16E–02	LILRB1-AS1	−1.87	0.013104
INS	−3.6	9.75E–11	hsa-miR-125b-5p	−1.43	2.29E–02	LINC01561	−2.05	7.28E–05
RD3L	−3.16	2.70E–10	hsa-miR-378d	−1.43	2.43E–02	LINC01550	−2.09	1.89E–05
VGLL1	−2.96	4.88E–06	hsa-miR-10a-5p	−1.41	2.62E–02	LINC01620	−2.16	2.77E–05
TH	−2.85	7.02E–15	hsa-miR-133a-3p	−1.38	3.15E–02	LINC01554	−2.19	5.61E–06
CLDN8	−2.83	3.13E–06	hsa-miR-101-3p	−1.37	3.12E–02	B3GALT5-AS1	−2.31	3.83E–06

Abbreviations: FC, fold change; FDR, false discovery rate.

Function enrichment for DEGs

The upregulated and downregulated DEGs were subjected to the DAVID to predict their functions. The results indicated that 16 GO biological process terms were obtained for the upregulated genes, mainly involving cell cycle (ie, E2F2, CDK1, CCNE1, BUB1, UBE2C, and CCNB1), while 27 GO biological process terms were enriched for the downregulated genes, mainly involving inflammatory response (PF4 and CXCR1) (Table 2 and Figure 2A). Furthermore, the KEGG pathway enrichment analysis was performed. In line with the GO enrichment results, the cell cycle (E2F2, CDK1, CCNE1, BUB1, and CCNB1) and p53 signaling pathway (CDK1, CCNE1, and CCNB1) were also obtained for the upregulated genes, while cytokine–cytokine receptor interaction was enriched for the downregulated genes (PF4 and CXCR1) (Table 3 and Figure 2B).

Figure 2 Function enrichment analyses for the differentially expressed genes.

Note: (A) GO enrichment and (B) KEGG pathways enrichment.

Abbreviations: GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Table 2 GO enrichment for differentially expressed genes using the DAVID database

Category	Term	FDR	Genes
Up	GO:0007049~cell cycle	4.24E–10	E2F1, KIF23, E2F2, KIFC1, XRCC2, E2F7, E2F8, MAEL, PKMYT1, TTK, PTTG1, AURKB, GTSE1, KIF2C, CCNE1, CDCA8, CDC45, CDCA2, PIWIL3, CDCA5, CDC6, CDK1, EGFL6, MND1, PBK, HMGA2, UBE2C, PRDM9, BUB1B, ERN2, NEK2, ANLN, CEP55, SPC24, SPC25, NCAPH, DUSP13, HJURP, NCAPG, CENPA, BUB1, FBXO43, SKA3, SKA1, TRIP13, EXO1, DLGAP5, NUF2, KIF18A, BIRC5, NDC80, CDC20, CDKN3, CDC25C, CDC25A, GSG2, CCNB1, CCNB2, TEX11
	GO:0022402~cell cycle process	1.13E–11	KIF23, E2F1, KIFC1, XRCC2, MAEL, PKMYT1, TTK, AURKB, PTTG1, GTSE1, CCNE1, KIF2C, CDCA8, CDCA2, PIWIL3, CDCA5, CDC6, CDK1, MND1, PBK, UBE2C, HMGA2, PRDM9, BUB1B, ERN2, NEK2, ANLN, CEP55, SPC24, SPC25, NCAPH, DUSP13, CENPA, NCAPG, FBXO43, BUB1, SKA3, SKA1, TRIP13, EXO1, DLGAP5, KIF18A, NUF2, NDC80, BIRC5, CDC20, CDC25C, CDKN3, CDC25A, CCNB1, CCNB2, TEX11
	GO:0022403~cell cycle phase	1.74E–15	E2F1, KIF23, KIFC1, XRCC2, NEK2, MAEL, TTK, PKMYT1, ANLN, PTTG1, CEP55, AURKB, GTSE1, SPC24, CCNE1, KIF2C, SPC25, NCAPH, CDCA8, DUSP13, NCAPG, BUB1, FBXO43, CDCA2, SKA3, PIWIL3, SKA1, CDCA5, TRIP13, EXO1, CDK1, CDC6, DLGAP5, NUF2, KIF18A, MND1, CDC20, BIRC5, NDC80, PBK, HMGA2, CDKN3, CDC25C, UBE2C, CDC25A, CCNB1, PRDM9, CCNB2, BUB1B, TEX11
	GO:000279~M phase	3.28E–16	KIF23, KIFC1, XRCC2, NEK2, MAEL, TTK, PKMYT1, ANLN, PTTG1, CEP55, AURKB, SPC24, KIF2C, SPC25, NCAPH, CDCA8, DUSP13, NCAPG, BUB1, FBXO43, CDCA2, SKA3, PIWIL3, SKA1, CDCA5, TRIP13, EXO1, CDK1, CDC6, DLGAP5, NUF2, KIF18A, MND1, CDC20, BIRC5, NDC80, PBK, HMGA2, CDC25C, UBE2C, CDC25A, CCNB1, PRDM9, CCNB2, BUB1B, TEX11
	GO:000278~mitotic cell cycle	2.13E–11	E2F1, KIF23, KIFC1, NEK2, TTK, PKMYT1, ANLN, PTTG1, CEP55, AURKB, GTSE1, SPC24, CCNE1, KIF2C, SPC25, CDCA8, NCAPH, NCAPG, CENPA, BUB1, CDCA2, SKA3, SKA1, CDCA5, CDK1, CDC6, DLGAP5, NUF2, KIF18A, CDC20, BIRC5, NDC80, PBK, HMGA2, CDKN3, CDC25C, UBE2C, CDC25A, CCNB1, CCNB2, BUB1B
	GO:0007067~mitosis	6.55E–14	KIF23, KIFC1, NEK2, PKMYT1, ANLN, CEP55, AURKB, PTTG1, SPC24, KIF2C, SPC25, CDCA8, NCAPH, NCAPG, BUB1, CDCA2, SKA3, SKA1, CDCA5, CDK1, CDC6, DLGAP5, NUF2, KIF18A, BIRC5, NDC80, CDC20, PBK, HMGA2, CDC25C, UBE2C, CDC25A, CCNB1, CCNB2, BUB1B
	GO:000280~nuclear division	6.55E–14	KIF23, KIFC1, NEK2, PKMYT1, ANLN, CEP55, AURKB, PTTG1, SPC24, KIF2C, SPC25, CDCA8, NCAPH, NCAPG, BUB1, CDCA2, SKA3, SKA1, CDCA5, CDK1, CDC6, DLGAP5, NUF2, KIF18A, BIRC5, NDC80, CDC20, PBK, HMGA2, CDC25C, UBE2C, CDC25A, CCNB1, CCNB2, BUB1B
	GO:0000087~M phase of mitotic cell cycle	9.83E–14	KIF23, KIFC1, NEK2, PKMYT1, ANLN, CEP55, AURKB, PTTG1, SPC24, KIF2C, SPC25, CDCA8, NCAPH, NCAPG, BUB1, CDCA2, SKA3, SKA1, CDCA5, CDK1, CDC6, DLGAP5, NUF2, KIF18A, BIRC5, NDC80, CDC20, PBK, HMGA2, CDC25C, UBE2C, CDC25A, CCNB1, CCNB2, BUB1B
	GO:0048285~organelle fission	1.97E–13	KIF23, KIFC1, NEK2, PKMYT1, ANLN, CEP55, AURKB, PTTG1, SPC24, KIF2C, SPC25, CDCA8, NCAPH, NCAPG, BUB1, CDCA2, SKA3, SKA1, CDCA5, CDK1, CDC6, DLGAP5, NUF2, KIF18A, BIRC5, NDC80, CDC20, PBK, HMGA2, CDC25C, UBE2C, CDC25A, CCNB1, CCNB2, BUB1B
	GO:0051301~cell division	2.47E–08	KIF23, KIFC1, NEK2, ANLN, CEP55, AURKB, PTTG1, SPC24, CCNE1, SPC25, CDCA8, NCAPH, NCAPG, CDCA2, BUB1, SKA3, POU3F2, SKA1, CDCA5, CDK1, CDC6, NUF2, BIRC5, NDC80, CDC20, HMGA2, UBE2C, CDC25C, CDC25A, CCNB1, CCNB2, BUB1B
Down	GO:0007267~cell–cell signaling	1.82E–07	EDN3, GABRB3, FCRL2, VIPR1, VIPR2, GDNF, WNT2, KCNQ5, WISP2, SLC1A2, GRIN2B, CHRNA4, EFNB3, NPBWR1, IL26, NRXN1, NTSR1, IL22, SIGLEC6, GRM7, WNT9A, DRD1, GRAP, OXT, DRD5, TH, MME, RIMS1, CCL24, INS, CCL21, PRIMA1, BMP3, IL6, PLP1, NOS1, NTF3, DLGAP2, GABRA5, NPY5R, KCNK3, CCL17, WNT7B, CXCL14, PNOC, GRIA1, NTRK2, ADRA1A, SLC5A7, WNT7A, IL2, HTR2A
	GO:0006952~defense response	2.83E–07	KLRC4, KLRC2, KLRC3, CXCR1, CFP, GRIN2B, HAMP, RNASE7, IFNG, XCR1, NLRP7, CAMP, PRG2, PSG3, IL22, NCR1, NCR3, CCR9, PROK2, PSG8, PPBP, CD40LG, GRM7, DEFA3, PLA2G2D, CTSG, NGF, KIR3DL2, CLEC1B, C7, DRD1, CCK, CCL24, AZU1, FCN3, INS, CCL21, FCN2, CNR2, SFTPD, IL1RAPL2, SELP, IL6, IL5, IL1RL1, GABRA5, CCL19, CD5L, STAB2, S100A12, CCL17, MPO, SELE
	GO:0044057~regulation of system process	3.67E–05	BMP10, EDN3, DRD1, CCK, ERBB4, MYL3, EDN2, DRD5, OXT, TH, GDNF, DES, GRIN2B, INS, IFNG, LGI1, ARC, GNAO1, NOS1, NTF3, NPY1R, NPY5R, PROK2, TNNT3, CHRM2, NTRK2, TBXA2R, AVPR1A, IL2, HTR2A, NGF
	GO:0007268~synaptic transmission	1.17E–04	DRD1, GABRB3, DRD5, OXT, TH, VIPR1, RIMS1, WNT2, KCNQ5, SLC1A2, GRIN2B, CHRNA4, PRIMA1, PLP1, NOS1, NTF3, DLGAP2, NPBWR1, GABRA5, NRXN1, NTSR1, NPY5R, KCNK3, PNOC, GRIA1, GRM7, SLC5A7, WNT7A, HTR2A
	GO:0051046~regulation of secretion	1.79E–03	EDN3, IL6, EDN2, OXT, FGF23, NPY1R, GDNF, NPY5R, GCK, GRIN2B, CD40LG, INS, GRM7, IFNG, NTRK2, AVPR1A, CHRNA4, TRPV6, IL2, HTR2A, NGF
	GO:0006935ĉhemotaxis	3.12E–03	EDN3, IL6, EDN2, CXCR1, CCL19, PF4, CCL17, CCR9, AZU1, CCL24, PROK2, PPBP, CXCL14, CCL21, IFNG, SFTPD, XCR1, LECT2
	GO:0042742~defense response to bacterium	8.64E–03	SELP, IL6, CAMP, PRG2, STAB2, S100A12, AZU1, CFP, PPBP, HAMP, RNASE7, IFNG, DEFA3, CTSG
	GO:0050900~leukocyte migration	9.20E–03	AZU1, EDN3, SELP, IL6, EDN2, IFNG, ELANE, SFTPD, PF4, SELE
	GO:0022610~biological adhesion	1.39E–02	CLDN8, CLSTN2, OPCML, DSCAML1, CLDN10, L1CAM, MEGF10, WISP2, SRPX, TNR, DPT, DSCAM, TECTA, SELP, MAG, HAPLN4, RET, CDHR1, PCDH11X, CDHR2, IGFALS, SIGLEC11, AJAP1, STAB2, NRXN1, PCDH19, CTNNA3, CLEC4M, DSG4, LYVE1, FREM3, HEPACAM, SIGLEC6, FREM2, CD40LG, DSG1, CDH19, ITGAD, CNTN3, SELE, COL20A1, IL2
	GO:0006954înflammatory response	2.80E–02	SELP, C7, IL6, IL5, CXCR1, CCL19, IL22, S100A12, CCL17, NCR3, AZU1, CFP, CCL24, PROK2, FCN3, CD40LG, INS, CCL21, FCN2, CNR2, XCR1, SELE, PLA2G2D, NGF

Note: Top 10 terms were listed.

Abbreviations: DAVID, Database for Annotation, Visualization and Integrated Discovery; FDR, false discovery rate; GO, Gene Ontology.

Table 3 KEGG pathway enrichment for differentially expressed genes using the DAVID database

Category	Term	FDR	Genes
Up	hsa04110:cell cycle	2.16E–08	E2F1, E2F2, CDK1, CDC6, PKMYT1, TTK, CDC20, PTTG1, SFN, CDC25C, CDC25A, CCNB1, CCNE1, CDC45, CCNB2, BUB1, BUB1B
	hsa04080:neuroactive ligand–receptor interaction	2.66E–04	GABRD, CGA, GABRG2, GABRA2, GABRA3, GLRA2, LEP, GRM4, SSTR5, KISS1R, HTR1B, GABRR1, GRID2, NPFFR2, TAAR1, HTR1D, GLP1R
	hsa04115:p53 signaling pathway	2.08E–02	CCNB1, CDK1, CCNE1, CCNB2, SFN, GTSE1
	hsa04060:cytokine–cytokine receptor interaction	2.69E–02	LEP, CCR8, CCL20, GDF5, EGF, BMP7, CCL7, IL11, CCL26
	hsa04062:chemokine signaling pathway	2.73E–02	CCR8, GNGT1, CCL20, GNG13, GNG4, CCL7, CCL26
	hsa00604:glycosphingolipid biosynthesis	2.99E–02	ST6GALNAC5, B4GALNT1
	hsa04512:ECM–receptor interaction	3.14E–02	IBSP, COMP, COL2A1, COL11A1
	hsa04350:TGF-beta signaling pathway	3.33E–02	COMP, GDF5, BMP7, PITX2
	hsa03320:PPAR signaling pathway	4.83E–02	LPL, MMP1, FABP6
Down	hsa00830:retinol metabolism	8.15E–10	CYP3A4, CYP1A1, CYP2B6, CYP2C8, ADH1B, CYP26A1, ADH7, CYP1A2, CYP3A43, CYP2A13, CYP4A11, UGT2B17, CYP4A22, LRAT, ADH4, CYP2A6, CYP2A7, RDH16
	hsa00982:drug metabolism	2.89E–05	CYP3A4, CYP2B6, CYP2C8, ADH1B, ADH7, CYP1A2, CYP2E1, GSTM5, CYP3A43, CYP2A13, UGT2B17, ADH4, CYP2A6, CYP2A7
	hsa04080:neuroactive ligand–receptor interaction	1.79E–04	GPR83, DRD1, GABRB3, DRD5, PTH1R, PRSS1, VIPR1, VIPR2, GCGR, GRIN2B, CNR2, GLP2R, GABRP, NPBWR1, GABRA5, NPY1R, NTSR1, NPY5R, GRIA1, CHRM2, GRM7, PTGDR, TBXA2R, AVPR1A, ADRA1A, CTSG, HTR2A
	hsa04060:cytokine–cytokine receptor interaction	2.07E–04	CXCR1, CNTFR, PF4, PF4V1, CCL24, CXCR5, CCL21, IFNG, XCR1, AMHR2, IL6, IL5, FLT3, TNFRSF13B, TNFRSF13C, IL26, CCL19, IL22, CCL17, CCR9, TSLP, CXCL14, PPBP, CD40LG, IL5RA, NGFR, IL2
	hsa00980:metabolism of xenobiotics by CYP/CYP450	3.64E–04	CYP3A43, CYP3A4, UGT2B17, CYP1A1, CYP2B6, ADH4, CYP2C8, ADH1B, ADH7, CYP2E1, CYP1A2, GSTM5
	hsa00232:caffeine metabolism	1.50E–03	CYP2A13, NAT2, CYP2A6, CYP2A7, CYP1A2
	hsa00983:drug metabolism	2.10E–02	CYP3A43, CYP3A4, CYP2A13, UGT2B17, NAT2, UPP2, CYP2A6, CYP2A7
	hsa04340:Hedgehog signaling pathway	2.00E–02	WNT2, WNT1, WNT10A, WNT7B, WNT9A, HHIP, GAS1, WNT7A, BMP5

Abbreviations: DAVID, Database for Annotation, Visualization and Integrated Discovery; FDR, false discovery rate; KEGG, Kyoto Encyclopedia of Genes and Genomes.

PPI network construction

Using the STRING database, 2,065 interaction relationship pairs (eg, BUB1–CDK1) were obtained between the 357 DEGs (182 downregulated and 175 upregulated), which were used for constructing a PPI network (Figure 3). After calculating the topological features for each protein in PPI network, CCNB1, CDK1, G protein subunit gamma 4 (GNG4), UBE2C, G protein subunit gamma transducin 1 (GNGT1), kinesin family member 4A (KIF4A), PF4, and G protein subunit gamma 13 (GNG13) were found to be shared in four topological characteristics and ranked in the top 30, suggesting that they may be hub genes (Table 4).

Figure 3 Protein and protein interaction network for the differentially expressed genes.

Notes: Red, upregulated; green, downregulated. The larger size of node (protein) indicates the higher degree (interaction relationships) of it.

Abbreviation: FC, fold change.

Four significant functionally related modules (Figure 4) were subsequently screened using the MCODE, among which the module 1 was the most significant with score =14.474 and node =38, followed by module 2 with score =13.152 and node =46. Also, GO analysis of genes in modules 1 and 2 with BinGO plugin of Cytoscape indicated that they were involved in mitotic cell cycle (CCNB1, CDK1, BUB1, and UBE2C) and cell surface receptor-linked signaling pathway (PF4 and CXCR1) (Table 5).

Figure 4 Modules extracted from the protein and protein interaction network.

Notes: (A) module 1; (B) module 2; (C) module 3; and (D) module 4. Red, upregulated; green, downregulated. The larger size of node (protein) indicates the higher degree (interaction relationships) of it.

Table 4 Topological features for each protein in PPI network

Gene	Degree	Gene	CC	Gene	BC	Gene	APL	Overlapped	Expression
GNGT1	87	IL6	0.3179	IL6	0.2515	IL6	3.1461	CCNB1	Up
GNG13	85	CHRM2	0.2989	UBE2C	0.1845	CHRM2	3.3455	CDK1	Up
GNG4	64	SH3GL2	0.2959	NGF	0.1695	SH3GL2	3.3792	GNG4	Up
CDK1	56	EGF	0.2947	TH	0.1586	EGF	3.3933	UBE2C	Up
CCNB1	52	ESR1	0.2890	GNGT1	0.1413	ESR1	3.4607	GNGT1	Up
CDC20	48	CCL20	0.2848	FOSB	0.1406	CCL20	3.5112	KIF4A	Up
BUB1	48	GNGT1	0.2846	MMP1	0.1388	GNGT1	3.5140	PF4	Down
CCNB2	47	GNG13	0.2823	CHRM2	0.1186	GNG13	3.5421	GNG13	Up
KIF2C	47	INS	0.2821	SH3GL2	0.1164	INS	3.5449
PMCH	45	UBE2C	0.2808	EGF	0.1106	UBE2C	3.5618
AURKB	44	PF4	0.2803	E2F1	0.1043	PF4	3.5674
CENPA	43	PTHLH	0.2786	ESR1	0.1013	PTHLH	3.5899
KIF20A	39	GNG4	0.2777	GNG13	0.0956	GNG4	3.6011
BUB1B	39	NGF	0.2766	INS	0.0904	NGF	3.6152
NCAPG	38	GNAO1	0.2728	KIF4A	0.0670	GNAO1	3.6657
CDCA8	38	KIF4A	0.2722	PTHLH	0.0669	KIF4A	3.6742
GCGR	37	PPBP	0.2713	CDK1	0.0637	PPBP	3.6854
NPSR1	37	GNA14	0.2709	ACAN	0.0632	GNA14	3.6910
BIRC5	37	LEP	0.2707	CCL20	0.0620	LEP	3.6938
TTK	37	TERT	0.2707	GNAO1	0.0593	TERT	3.6938
NDC80	37	GNAZ	0.2669	CAMK2B	0.0547	GNAZ	3.7472
DLGAP5	36	NGFR	0.2661	LPL	0.0444	NGFR	3.7584
UBE2C	35	SYT9	0.2655	CYP19A1	0.0430	SYT9	3.7669
CDC45	35	E2F1	0.2651	CCNB1	0.0417	E2F1	3.7725
PF4	34	HTR1B	0.2645	GNG4	0.0402	HTR1B	3.7809
KIF4A	34	CDK1	0.2639	CHRNA1	0.0394	CDK1	3.7893
PBK	34	PMCH	0.2639	PF4	0.0393	PMCH	3.7893
KIF23	34	HTR1D	0.2639	RHO	0.0370	HTR1D	3.7893
CXCR1	33	CCNB1	0.2625	NOS1	0.0353	CCNB1	3.8090
XCR1	33	CXCR1	0.2597	NGFR	0.0348	CXCR1	3.8511

Note: Top 30 genes were listed.

Abbreviations: APL, average path length; BC, betweenness centrality; CC, closeness centrality; PPI, protein–protein interaction.

Table 5 BinGO enrichment for genes in modules

Module	Description	FDR	Genes in test set
M1	Mitotic cell cycle	9.74E–37	CDCA5|BUB1B|CDCA8|NCAPG|TTK|CENPA|SKA1|AUR KB|CDC20|CCNB2|CCNB1|PTTG1|NUF2|PBK|NEK2|BU B1|CEP55|DLGAP5|UBE2C|KIF23|NDC80|ANLN|KIF18 A|CDK1|BIRC5|KIF2C|SPC24|CDKN3|SPC25
	Nuclear division	9.74E–37	CDCA5|BUB1B|CDCA8|NCAPG|SKA1|AURKB|CDC20|C CNB2|CCNB1|PTTG1|NUF2|PBK|NEK2|BUB1|CEP55|D LGAP5|UBE2C|KIF23|NDC80|ANLN|KIF18A|CDK1|BIRC 5|KIF2C|SPC24|SPC25
	Mitosis	9.74E–37	CDCA5|BUB1B|CDCA8|NCAPG|SKA1|AURKB|CDC20|C CNB2|CCNB1|PTTG1|NUF2|PBK|NEK2|BUB1|CEP55|D LGAP5|UBE2C|KIF23|NDC80|ANLN|KIF18A|CDK1|BIRC 5|KIF2C|SPC24|SPC25
	M phase of mitotic cell cycle	1.65E–36	CDCA5|BUB1B|CDCA8|NCAPG|SKA1|AURKB|CDC20|C CNB2|CCNB1|PTTG1|NUF2|PBK|NEK2|BUB1|CEP55|D LGAP5|UBE2C|KIF23|NDC80|ANLN|KIF18A|CDK1|BIRC 5|KIF2C|SPC24|SPC25
	Organelle fission	1.65E–36	CDCA5|BUB1B|CDCA8|NCAPG|SKA1|AURKB|CDC20|C CNB2|CCNB1|PTTG1|NUF2|PBK|NEK2|BUB1|CEP55|D LGAP5|UBE2C|KIF23|NDC80|ANLN|KIF18A|CDK1|BIRC 5|KIF2C|SPC24|SPC25
	M phase	2.75E–36	CDCA5|BUB1B|CDCA8|NCAPG|TTK|SKA1|AURKB|CDC 20|CCNB2|CCNB1|PTTG1|NUF2|PBK|NEK2|BUB1|CEP 55|DLGAP5|UBE2C|KIF23|NDC80|ANLN|KIF18A|CDK1 |BIRC5|KIF2C|TRIP13|SPC24|SPC25
	Cell cycle phase	1.03E–35	CDCA5|BUB1B|CDCA8|NCAPG|TTK|SKA1|AURKB|CDC 20|CCNB2|CCNB1|PTTG1|NUF2|PBK|NEK2|BUB1|CEP 55|DLGAP5|UBE2C|KIF23|NDC80|ANLN|KIF18A|CDK1 |BIRC5|KIF2C|TRIP13|SPC24|CDKN3|SPC25
	Cell cycle process	5.27E–34	CDCA5|BUB1B|CDCA8|NCAPG|TTK|CENPA|SK A1|AURKB|CDC20|CCNB2|CCNB1|PTTG1|NUF 2|PBK|NEK2|BUB1|CEP55|DLGAP5|UBE2C|KIF 23|NDC80|ANLN|KIF18A|CDK1|BIRC5|KIF2C| TRIP13|SPC24|CDKN3|SPC25
	Cell cycle	7.78E–34	CDCA5|HJURP|BUB1B|CDCA8|NCAPG|TTK|CENPA|SK A1|AURKB|CDC20|CCNB2|CCNB1|CDC45|PTTG1|NU F2|PBK|NEK2|BUB1|CEP55|DLGAP5|UBE2C|KIF23|ND C80|ANLN|KIF18A|CDK1|BIRC5|KIF2C|TRIP13|SPC24| CDKN3|SPC25
	Cell division	8.53E–28	UBE2C|CDCA5|BUB1B|CDCA8|NCAPG|KIF23|SKA1|A URKB|NDC80|CDC20|CCNB2|ANLN|CCNB1|PTTG1| NUF2|CDK1|BIRC5|NEK2|KIF2C|BUB1|CEP55|SPC24 |SPC25
M2	G-protein-coupled receptor protein signaling pathway	4.81E–36	NPFFR2|CHRM2|PMCH|CXCR5|HTR2A|ADRA1A|GNGT 1|GNA14|GRM4|CXCR1|CNR2|TBXA2R|NPBWR1|KISS 1R|CCR9|CCR8|PROK2|TAC3|NTSR1|P2RY12|XCR1|ED N2|NPY5R|GCGR|HTR1D|NPY1R|HTR1B|CCK|AVPR1A| SSTR5|GNG13|PPYR1|PNOC
	Cell surface receptor-linked signaling pathway	8.03E–26	NPFFR2|CHRM2|PMCH|CXCR5|HTR2A|ADRA1A|GNGT 1|GNA14|GRM4|CXCR1|CNR2|TBXA2R|NPBWR1|KISS 1R|CCR9|CCR8|PROK2|TAC3|NTSR1|P2RY12|XCR1|ED N2|NPY5R|EDN3|GCGR|HTR1D|NPY1R|HTR1B|CCK|A VPR1A|SSTR5|GNG13|PPYR1|PNOC|PF4
	Signaling	6.64E–24	NPFFR2|CHRM2|PMCH|CXCR5|OXT|HTR2A|ADRA1A |GNGT1|GNA14|GRM4|CXCR1|CNR2|GRM7|TBXA2R| GNG4|NPBWR1|KISS1R|CCR9|CCR8|PROK2|TAC3|CCL 19|NTSR1|P2RY12|XCR1|EDN2|NPY5R|CCL21|EDN3| CCL20|GCGR|HTR1D|NPY1R|HTR1B|CCK|AVPR1A|SST R5|GNG13|APLN|PPYR1|PNOC|GNRH2|PF4
	Behavior	1.44E–23	PMCH|OXT|HTR2A|GRM4|CXCR1|CNR2|KISS1R|CCR 9|CCR8|PROK2|CCL19|NTSR1|XCR1|EDN2|NPY5R|C CL21|EDN3|CCL20|NPY1R|HTR1B|CCK|AVPR1A|PPBP |PPYR1|PF4
	Signaling pathway	3.33E–20	NPFFR2|CHRM2|PMCH|CXCR5|HTR2A|ADRA1A|GNGT 1|GNA14|GRM4|CXCR1|CNR2|TBXA2R|GNG4|NPBW R1|KISS1R|CCR9|CCR8|PROK2|TAC3|NTSR1|P2RY12|X CR1|EDN2|NPY5R|EDN3|GCGR|HTR1D|NPY1R|HTR1B |CCK|AVPR1A|SSTR5|GNG13|PPYR1|PNOC|PF4
	Signaling process	2.04E–17	CHRM2|PMCH|OXT|HTR2A|ADRA1A|GNGT1|GNA14| GRM4|CNR2|GRM7|NPBWR1|KISS1R|PROK2|CCL19|N TSR1|P2RY12|XCR1|EDN2|NPY5R|CCL21|EDN3|CCL2 0|GCGR|HTR1D|NPY1R|HTR1B|CCK|AVPR1A|SSTR5|G NG13|APLN|PNOC|GNRH2|PF4
	Signal transmission	2.04E–17	CHRM2|PMCH|OXT|HTR2A|ADRA1A|GNGT1|GNA14| GRM4|CNR2|GRM7|NPBWR1|KISS1R|PROK2|CCL19|N TSR1|P2RY12|XCR1|EDN2|NPY5R|CCL21|EDN3|CCL2 0|GCGR|HTR1D|NPY1R|HTR1B|CCK|AVPR1A|SSTR5|G NG13|APLN|PNOC|GNRH2|PF4
	Signal transduction	1.43E–14	CHRM2|PMCH|OXT|HTR2A|ADRA1A|GNGT1|GNA14| GRM4|CNR2|KISS1R|PROK2|CCL19|P2RY12|XCR1|ED N2|CCL21|EDN3|CCL20|GCGR|HTR1D|NPY1R|HTR1B| CCK|AVPR1A|SSTR5|GNG13|APLN|PNOC|GNRH2|PF4
	Response to stimulus	2.05E–12	NPFFR2|CHRM2|PMCH|OXT|HTR2A|ADRA1A|GNGT1| GRM4|CXCR1|CNR2|GRM7|GNG4|KISS1R|CCR9|CCR8 |PROK2|CCL19|NTSR1|P2RY12|XCR1|EDN2|NPY5R|C CL21|EDN3|CCL20|GCGR|HTR1D|NPY1R|HTR1B|CCK| AVPR1A|PPBP|GNG13|APLN|PPYR1|PF4
	Response to chemical stimulus	2.96E–11	XCR1|EDN2|CCL21|EDN3|CCL20|GCGR|HTR1D|NP Y1R|HTR1B|OXT|CCK|HTR2A|AVPR1A|PPBP|ADRA1A |GNG13|CXCR1|CNR2|GNG4|CCR9|CCR8|PROK2|CC L19|PF4
M3	Cyclic-nucleotide-mediated signaling	6.43E–09	GLP1R|VIPR1|GLP2R|PTH1R|DRD1|PTHLH|DRD5
	G-protein-coupled receptor protein signaling pathway	2.93E–08	VIPR1|GPR15|GLP2R|PTH1R|DRD1|TAAR1|PTHLH|DR D5|PTGDR
	Second messenger-mediated signaling	9.78E–08	GLP1R|VIPR1|GLP2R|PTH1R|DRD1|PTHLH|DRD5
	G-protein signaling, coupled to cyclic nucleotide second messenger	9.78E–08	VIPR1|GLP2R|PTH1R|DRD1|PTHLH|DRD5
	Cell surface receptor-linked signaling pathway	2.27E–05	VIPR1|GPR15|GLP2R|PTH1R|DRD1|TAAR1|PTHLH|DR D5|PTGDR
	Signaling	3.01E–05	GLP1R|VIPR1|GPR15|VIPR2|GLP2R|PTH1R|DRD1|CGA|T AAR1|PTHLH|DRD5|PTGDR
	Signaling pathway	5.37E–05	GLP1R|VIPR1|GPR15|GLP2R|PTH1R|DRD1|TAAR1|PTH LH|DRD5|PTGDR
	Intracellular signal transduction	8.25E–05	GLP1R|VIPR1|GLP2R|PTH1R|DRD1|PTHLH|DRD5
	Cell–cell signaling	1.41E–04	VIPR1|VIPR2|DRD1|CGA|PTHLH|DRD5
	Signal transduction	1.59E–04	GLP1R|VIPR1|VIPR2|GLP2R|PTH1R|DRD1|CGA|PTHLH |DRD5
	Signaling process	4.51E–04	GLP1R|VIPR1|VIPR2|GLP2R|PTH1R|DRD1|CGA|PTHLH |DRD5
	Cell communication	7.10E–04	VIPR1|VIPR2|DRD1|CGA|PTHLH|DRD5
M4	Drug metabolic process	3.44E–16	CYP2A6|CYP2C8|CYP2B6|CYP1A2|CYP1A1|CYP2E1|C YP3A4
	Secondary metabolic process	4.89E–11	CYP26A1|CYP2A6|CYP1A2|CYP1A1|CYP2E1|CYP3A4
	Oxidation reduction	2.91E–10	CYP26A1|CYP2A6|CYP2C8|CYP2B6|CYP4A11|CYP1A2|C YP1A1|CYP2E1|CYP3A4
	Steroid metabolic process	4.84E–10	CYP2A6|CYP2B6|CYP1A2|CYP1A1|UGT2B17|CYP2E1| CYP3A4
	Lipid metabolic process	1.79E–09	CYP26A1|CYP2A6|CYP2B6|CYP4A11|CYP1A2|CYP1A1| UGT2B17|CYP2E1|CYP3A4
	Response to drug	1.12E–07	CYP2A6|CYP2C8|CYP2B6|CYP1A2|CYP1A1|CYP3A4
	Cellular catabolic process	2.17E–06	CYP26A1|CYP2A6|CYP2C8|CYP2B6|CYP1A2|CYP1A1| CYP3A4
	Small molecule metabolic process	3.78E–06	CYP26A1|CYP2A6|CYP2C8|CYP2B6|CYP4A11|CYP1A2| CYP1A1|CYP3A4
	Cellular lipid metabolic process	6.87E–06	CYP26A1|CYP4A11|CYP1A2|CYP1A1|CYP2E1|CYP3A4
	Catabolic process	9.30E–06	CYP26A1|CYP2A6|CYP2C8|CYP2B6|CYP1A2|CYP1A1| CYP3A4

Abbreviation: FDR, false discovery rate.

CeRNA network construction

By searching the miRcode database, 14 DELs–DEMs interaction relationship pairs (including five DELs, all upregulated, and 10 DEMs, five upregulated and five down-regulated) were predicted, which were used for constructing the DELs–DEMs network. Subsequently, the target genes of these 10 DEMs were predicted with the miRwalk database. After removal of the positive–negative relationships between DEMs and DEGs, 113 DELs–DEMs interaction relationship pairs (including eight DEMs, three upregulated and five downregulated, and 82 DEGs, 67 upregulated and 15 downregulated) were left for constructing the DEMs–DEGs network. By integrating the DELs–DEMs network and DEMs–DEGs network, a DELs–DEMs–DEGs ceRNA network was established containing 95 nodes (five DELs, eight DEMs, and 82 DEGs) and 239 edges (14 DELs–DEMs, 113 DELs–DEGs, and 112 DEGs–DEGs) (Figure 5).

Figure 5 ceRNAs interaction network of lncRNA–miRNA–mRNA.

Notes: Square nodes represent lncRNAs; triangle nodes represent miRNAs; circular nodes represent mRNAs. Red, upregulated; green, downregulated.

Abbreviation: FC, fold change.

Function enrichment analysis with DAVID showed the genes in the ceRNA network participated in four significant KEGG pathways, including cell cycle, p53 signaling pathway, neuroactive ligand–receptor interaction, and pathways in cancer (Table 6). By searching the CTD database with “Hepatocellular Carcinoma” as the keyword, 244 KEGG pathways were found to be associated with HCC. Among them, three were common with the enrichment results of the genes in the ceRNA network, including cell cycle (CCNE1, E2F2, and CDK1), p53 signaling pathway (CCNE1 and CDK1), and pathways in cancer (E2F2). Thus, the DELs–DEMs–DEGs interaction relationship pairs associated with these three pathways were extracted to form the HCC-related ceRNA network (Figure 6), in which four DELs (DSCR4, FAM182B, PRNT, and LINC00346), five DEMs (hsa-miR-199a-5p, hsa-miR-30a-5p, hsa-miR-125b-5, hsa-miR-10a-5p, and hsa-miR-133a-3p), and seven DEGs (CDC6, CCNE1, CDK1, E2F2, BUB1, CCNB1, and SFN) were involved.

Figure 6 HCC-related ceRNAs interaction network of lncRNA–miRNA–mRNA.

Notes: Square nodes represent circRNAs, triangle nodes represent miRNAs, circular nodes represent mRNAs, and rhombus nodes represent HCC pathways. Red, upregulated; green, downregulated.

Abbreviations: FC, fold change; HCC, hepatocellular carcinoma.

Table 6 KEGG pathways for genes in ceRNA network

Term	P-value	Genes
hsa04110:cell cycle	2.85E–06	CCNB1, E2F2, CDK1, CCNE1, CDC6, BUB1, SFN
hsa04115:p53 signaling pathway	0.001617	CCNB1, CDK1, CCNE1, SFN
hsa04080:neuroactive ligand–receptor interaction	0.04606	GPR83, GABRB3
hsa05200:pathways in cancer	0.046995	E2F2, CCNE1

Abbreviation: KEGG, Kyoto Encyclopedia of Genes and Genomes.

Prognosis prediction for DELs, DEMs, and DEGs

Ninety-eight HBV-related HCC samples, which have been used for mRNA and miRNA sequencing, were collected from TCGA database. Univariate Cox regression analysis was then used to screen OS- and RFS-related DELs, DEMs, and DEGs from HCC-related ceRNA network in these samples. The results showed that two DELs, four DEMs, and seven DEGs were significantly associated with OS, but only five DEGs were significantly associated with RFS (Table 7). KM curve was subsequently drawn according the expression level of each DEL, DEM, and DEG in the sequencing data. In line with the Cox regression analysis results, KM curve analysis also (Figure 7) showed that two DELs (FAM182B and LINC00346), four DEMs (hsa-miR-30a-5p, hsa-miR-125b-5p, hsa-miR-10a-5p, and hsa-miR-133a-3p), and seven DEGs (CDC6, CCNE1, CDK1, E2F2, BUB1, CCNB1, and SFN) were significantly associated with OS but not PRNT (expression value =0 in TCGA data), DSCR4 (P=0.493), and hsa-miR-199a-5p (P=0.101). Also, all the relationships between their expressions and the prognosis results were in line with our expectation, that is, patients with the high expression of the DELs and DEGs (all were upregulated genes in HBV-related HCC) had the poor survival, while patients with the high expression of DEMs (all were downregulated genes in HBV-related HCC) possessed excellent survival. As shown in Figure 8, KM curve analysis also showed that the highly expressed five DEGs (CDC6, CDK1, BUB1, CCNB1, and SFN) were significantly associated with RFS.

Figure 7 Kaplan–Meier analysis to display the correlation of differentially expressed lncRNAs (A), miRNAs (B), and genes (C) with overall survival outcomes for patients with HBV-related HCC.

Abbreviations: HBV, hepatitis B virus; HCC, hepatocellular carcinoma.

Figure 8 Kaplan–Meier analysis to display the correlation of differentially expressed genes with recurrence-free survival outcomes for patients with HBV-related HCC.

Abbreviations: HBV, hepatitis B virus; HCC, hepatocellular carcinoma.

Further combination with their interaction relationships in the ceRNA network suggested that FAM182B-miR-125b-5p-E2F2 and LINC00346-miR-10a-5p-CDK1/CCNE1 ceRNA axes were especially important for the development and prognosis of HBV-related HCC.

Table 7 Cox regression analysis to screen survival-related genes

ID	Overall survival		Recurrence-free survival
	HR	P-value	HR	P-value
E2F2	1.24	0.049	1.19	0.16
CDC6	1.28	0.0372	1.41	0.0099
CCNE1	1.05	0.0456	1.08	0.32
CDK1	1.32	0.022	1.48	0.0036
BUB1	1.31	0.013	1.37	0.0054
CCNB1	1.4	0.014	1.36	0.023
SFN	1.22	0.0023	1.19	0.0055
hsa-miR-10a-5p	0.882	0.04	0.854	0.27
hsa-miR-125b-5p	0.835	0.019	0.952	0.74
hsa-miR-133a-3p	0.11	0.045	0.879	0.39
hsa-miR-199a-5p	0.986	0.88	0.907	0.28
hsa-miR-30a-5p	0.911	0.0452	0.866	0.32
LINC00346	1.67	0.0051	0.99	0.85
DSCR4	1.03	0.68	0.994	0.95
FAM182B	1.152	0.047	0.84	0.19

Validation of expressions of crucial DELs, DEMs, and DEGs in ceRNA network

The upregulation of LINC00346, CDK1, and CCNE1 but the downregulation of miR-10a-5p and miR-125b-5p was also validated in other microarray datasets and TCGA data. FAM182B was not found to be differentially expressed in GSE27462 and TCGA data. E2F2 was demonstrated to be differentially expressed in GSE94660, GSE25599, and TCGA data but not in GSE121248 (Table 8). These findings indicated that LINC00346-miR-10a-5p-CDK1 ceRNA axis may be a potentially verifiable mechanism for HBV-related HCC.

Table 8 Confirmation of expressions of crucial lncRNAs, miRNAs, and mRNAs using other datasets

RNA type	Dataset	Symbol	Tumor (mean ± SD)	Control (mean ± SD)	P-value
lncRNA	GSE27462	FAM182B	75.834±27.913	71.793±23.127	0.8096
		LINC00346	124.665±28.955	72.848±14.226	0.01208
	TCGA	FAM182B	3.403±1.255	2.539±0.944	0.081
		LINC00346	14.834±3.436	12.679±1.403	0.0397
miRNA	GSE69580	miR-125b-5p	39.384±23.416	355.428±82.423	0.000606
		miR-10a-5p	4.229±2.103	10.111±2.647	0.00507
	TCGA	miR-125b-5p	8.962±1.092	10.203±0.352	4.09E–08
		miR-10a-5p	13.607±1.198	14.779±0.472	4.89E–06
mRNA	GSE121248	E2F2	4.252±0.181	4.239±0.157	0.696
		CDK1	7.224±1.105	5.379±0.799	2.68793E–16
		CCNE1	7.056±1.032	6.415±0.216	0.00000344
	GSE94660	E2F2	0.83±0.443	0.0611±0.028	1.234E–07
		CDK1	2.141±0.769	0.251±0.121	2.9E–10
		CCNE1	1.412±0.669	0.242±0.129	9.01E–08
	GSE25599	E2F2	0.769±0.496	0.302±0.255	0.01958
		CDK1	5.116±2.925	1.152±0.923	0.005978
		CCNE1	2.015±1.820	0.430±0.242	0.02254
	TCGA	E2F2	5.612±1.296	2.959±1.076	0.000926
		CDK1	8.018±1.428	4.579±1.266	0.00374
		CCNE1	6.322±1.993	3.162±0.979	0.00193

Abbreviation: TCGA, The Cancer Genome Atlas.

Discussion

In the present study, we identified FAM182B-miR-125b-5p-E2F2 and LINC00346-miR-10a-5p-CDK1/CCNE1 ceRNA axes as important mechanisms for the development of HBV-related HCC. They were involved in HBV-related HCC by influencing cell cycle. Also, the genes in these two axes were significantly associated with the OS of patients. LINC00346-miR-10a-5p-CDK1 may be especially crucial because CDK1 was considered as a hub gene in the PPI network and was also associated with RFS as well as the expressions of all of them confirmed in other datasets.

Numerous studies have shown that HBV infection of hepatocytes promotes cell cycle progression by accelerating G1/S and G2/M transition and thus increases cell proliferation ability, ultimately inducing the development of HCC.^35,36 It is well accepted that CCNE1 is a positive regulator of G1/S phase transition³⁷ and CCNB1 is required for G2/M transition and mitosis resumption by forming a maturation promoting factor with CKD1.³⁸ Transcriptional factor E2F2 can be activated by Cyclin-CDK enzymatic complex after phosphorylating the protein retinoblastoma (Rb), which promotes the transcription of E2F2 target genes to regulate the G1/S-phase transition.³⁹ Thus, CCNE1, CCNB1, CKD1, and E2F2 genes are suggested to be upregulated in HBV-related HCC. These hypotheses have been demonstrated by previous studies. For example, Sung et al⁴⁰ used the RNA sequencing (RNA-seq) and Sanger sequencing to confirm that CCNE1 gene was highly expressed in HBV integrated tumors compared with adjacent normal tissue. Chin et al⁴¹ observed that delivery of a replication competent HBV genome into hepatocyte lines Huh7 and PMH induced the expression of CCNB1. Cheng et al⁴² also used in vitro experiments to prove HBV persistently activated the CCNB1-CDK1 kinase in HCC cells. In line with these studies, our study also found CCNE1, CCNB1, and CKD1 were upregulated in tumor samples of patients with HBV-related HCC and the high expression of them predicted poor prognosis. The CKD1 may be especially important because it was associated with both OS and RFS. Although there was a study to indicate E2F2 upregulation in HCC,⁴³ its relationship with HBV has not been investigated. Our study may be the first to reveal that HBV infection may trigger E2F2 upregulation and lead to the development of HCC and poor prognosis for patients. miRNAs are the class of small RNAs (18–25 nucleotides) that downregulate target gene expressions via binding to the 3′-untranslated region (UTR). Thus, the upregulation of cell cycle-related genes may be attributed to the downregulation of miRNAs. In this study, we also investigated the DEMs between tumor and normal samples and predicted their interaction with target genes by the miRwalk database. Our results indicated that miR-10a-5p and miR-125b-5p could regulate CDK1/CCNE1 and E2F2, respectively. There have studies to explore the miRNAs to regulate these target genes in HCC, such as miR-7/497/195-CCNE1,^44,45 miR-582-5p-CDK1,⁴⁶ and miR-214/490-E2F2,^47,48 but not focused on the relationships of our prediction. However, the studies on the expressions of miR-10a-5p and miR-125b-5p in HCC may indirectly illuminate their underlying negative relations. Zhu et al⁴⁹ identified the DEMs in seven paired specimens of HCC using the microarray technique and found that miR-10a-5p and miR-125b-5p were significantly downregulated. Over-expression of miR-10a⁵⁰ and miR-125b⁵¹ was reported to suppress the metastasis of HCC cells in vivo. In line with these studies, our study also showed that miR-10a-5p and miR-125b-5p were downregulated in HBV-related HCC and high expression of them predicted excellent prognosis. lncRNAs are proposed to act as a ceRNA to involve in the regulation effects of miRNAs on the expression of target genes. Thus, the upregulation of cell cycle-related genes may also be attributed to the upregulation of lncRNAs that sponged the miRNAs. In this study, we also investigated the DELs between tumor and normal samples and predicted their interaction with miRNAs by the miRcode database. Our results indicated that upregulated FAM182B and LINC00346 may regulate cell cycle-related genes by interacting with miR-125b-5p and miR-10a-5p, resulting in poor prognosis. In line with our study, there has been a study to demonstrate that LINC00346 was upregulated in bladder cancer tissues compared to normal tissues. Knockdown of LINC00346 inhibited bladder cancer cell proliferation and migration and induced cell cycle arrest and cell apoptosis.⁵² The high expression of LINC00346 was also found to be significantly associated with poor OS in HCC⁵³ and breast cancer samples.⁵⁴ Nevertheless, no studies were performed to investigate the interaction of LINC00346 with miRNAs. Also, any investigation related to FAM182B has not been found until now. These implied that our identified ceRNA axes (FAM182B-miR-125b-5p-E2F2 and LINC00346-miR-10a-5p-CDK1/CCNE1) may be novel mechanisms for HBV-related HCC.

There were some limitations in this study. First, our study has preliminarily predicted that these ceRNA axes may be associated with the development of HBV-related HCC and some of them were confirmed in some other microarray datasets. Thus, further clinical, in vitro (dual luciferase reporter assay), and in vivo (loss of function) experiments are necessary to validate the expressions of controversial genes (such as FAM182B and E2F2) and regulatory relationships between DELs and DEMs as well as between DEMs and DEGs, and their roles for the proliferation, metastasis, and invasion of HBV infection hepatocytes. Second, there were no clinical data in our used datasets (GSE77509 and GSE76903) and, thus, we only preliminarily predicted the associations between prognosis (OS and RFS) and each of our identified DEL/DEM/DEG using TCGA data via univariate cox regression analysis. Whether these genes were independent biomarkers needed further clinical trials with multivariate Cox’s model that integrated all the clinical information (such as HBV DNA level, liver function parameters, pathologic stage, pathologic node, and pathologic metastasis, grade, therapeutic strategies including hepatectomy, radiofrequency ablation, and solafenib)^55–57 and all DELs/DEMs/DEGs expression levels.

Conclusion

The present study preliminarily indicates that FAM182B and LINC00346 may be novel prognostic biomarkers and therapeutic targets for HBV-associated HCC. They function as a ceRNA to sponge miR-125b-5p and miR-10a-5p to de-repress cell cycle-related genes (E2F2, CDK1, and CCNE1) and promote the cell growth of HCC cells.

Ethics approval and informed consent

As the data used in this study were downloaded from GEO or TCGA database, and no human experiment was involved in this study, there were no ethical approval and informed consent. Thus, the principles of the Declaration of Helsinki were also not followed.

Availability of data and materials

All the microarray data were downloaded from the GEO database in NCBI (http://www.ncbi.nlm.nih.gov/geo/). The mRNA and miRNA Seq-data were obtained from TCGA (https://tcga-data.nci.nih.gov/).

Author contributions

HL and ZB participated in the conception and design of this study. HL and XZ performed the acquisition of data. CL and CS were involved in the analysis and interpretation of data. HL and CL performed the statistical and bioinformatics analyses. HL drafted the manuscript. ZB revised the manuscript for important intellectual content. All authors read and approved the final manuscript. All authors contributed toward data analysis, drafting and critically revising the paper and agree to be accountable for all aspects of the work.

Disclosure

The authors report no conflicts of interest in this work.