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 [33]GSE77509 and [34]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.[35]^1
   In 2015, it is estimated that there are 466,100 new cases and 422,100
   deaths due to this disease.[36]^1 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%).[37]^2 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.[38]^3^-[39]^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.[40]^6^,[41]^7 For example, Zuo et al found that lncRNA
   [42]AX800134 was upregulated in HBV-positive HCC compared with
   HBV-negative HCC. Silencing [43]AX800134 with siRNA interference
   significantly suppressed the growth and invasion but enhanced
   spontaneous apoptosis of HBx-expressing HepG2 cells.[44]^8 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.[45]^9 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.[46]^10 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.[47]^11 Accumulating data
   also indicated that this regulatory action plays important roles in HCC
   development.[48]^12^,[49]^13 For example, Lv et al[50]^14 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[51]^15
   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.[52]^16
   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.[53]^17 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.[54]^18 However, studies performed to investigate the
   ceRNA mechanisms of lncRNAs for HBV-related HCC were
   rare.[55]^14^,[56]^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 [57]GSE77509[58]^19 and
   [59]GSE76903[60]^19 were retrieved from the Gene Expression Omnibus
   (GEO) database ([61]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; [62]http://www.genenames.org/).[63]^20

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;
   [64]http://www.bioconductor.org/packages/release/bioc/html/edgeR.html).
   [65]^21 P-value was adjusted to false discovery rate (FDR) with
   multitest package (Version 2.36.0;
   [66]http://bioconductor.org/packages/release/bioc/html/multtest.html).[
   67]^22 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;
   [68]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; [69]http://stringdb.org/) database[70]^23 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;
   [71]www.cytoscape.org/).[72]^24 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
   ([73]http://apps.cytoscape.org/apps/cytonca)[74]^25 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,
   [75]http://apps.cytoscape.org/apps/mcode) algorithm[76]^26 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;
   [77]http://david.abcc.ncifcrf. gov)[78]^27 and BinGO[79]^28 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; [80]http://www.mircode.org/)[81]^29
   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; [82]www.cytoscape.org/).[83]^24 The
   target genes of DEMs in the DELs–DEMs interaction network were
   predicted using the miRwalk database (Version 2.0;
   [84]http://www.zmf.umm. uni-heidelberg.de/apps/zmf/mirwalk2),[85]^30
   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;
   [86]www.cytoscape.org/).[87]^24

   Furthermore, all known HCC-related pathways were downloaded from
   Comparative Toxicogenomics Database (CTD;
   [88]http://ctd.mdibl.org/),[89]^31 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 ([90]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;
   [91]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
   ([92]GSE121248: 70 vs 37; [93]GSE94660[94]^32: 21 vs 21;
   [95]GSE25599[96]^33: 10 vs 21 normal), miRNA ([97]GSE69580: 5 vs 5),
   and lncRNA ([98]GSE27462[99]^34: 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 ([100]Table 1). The heat map analysis showed that the
   samples with similar features tended to be clustered according to the
   expressions of DELs ([101]Figure 1A), DEMs ([102]Figure 1B), and DEGs
   ([103]Figure 1C).

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
   [104]Open in a new tab

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

Figure 1.

   [105]Figure 1
   [106]Open in a new tab

   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.

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) ([107]Table 2
   and [108]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) ([109]Table 3 and [110]Figure 2B).

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:0000279~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:0000278~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:0000280~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
   [111]Open in a new tab

   Note: Top 10 terms were listed.

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

Figure 2.

   [112]Figure 2
   [113]Open in a new tab

   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 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
   [114]Open in a new tab

   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
   ([115]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 ([116]Table 4).

Figure 3.

   [117]Figure 3
   [118]Open in a new tab

   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.

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
   [119]Open in a new tab

   Note: Top 30 genes were listed.

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

   Four significant functionally related modules ([120]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) ([121]Table
   5).

Figure 4.

   Figure 4
   [122]Open in a new tab

   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 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
   [123]Open in a new tab

   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) ([124]Figure 5).

Figure 5.

   [125]Figure 5
   [126]Open in a new tab

   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 ([127]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 ([128]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.

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
   [129]Open in a new tab

   Abbreviation: KEGG, Kyoto Encyclopedia of Genes and Genomes.

Figure 6.

   [130]Figure 6
   [131]Open in a new tab

   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.

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 ([132]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 ([133]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 [134]Figure 8, KM curve analysis also showed that the
   highly expressed five DEGs (CDC6, CDK1, BUB1, CCNB1, and SFN) were
   significantly associated with RFS.

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
   [135]Open in a new tab

Figure 7.

   [136]Figure 7
   [137]Open in a new tab

   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.

   [138]Figure 8
   [139]Open in a new tab

   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.

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 [140]GSE27462 and TCGA data. E2F2 was demonstrated to be
   differentially expressed in [141]GSE94660, [142]GSE25599, and TCGA data
   but not in [143]GSE121248 ([144]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 [145]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 [146]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 [147]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
   [148]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
   [149]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
   [150]Open in a new tab

   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.[151]^35^,[152]^36 It is well accepted that CCNE1 is
   a positive regulator of G1/S phase transition[153]^37 and CCNB1 is
   required for G2/M transition and mitosis resumption by forming a
   maturation promoting factor with CKD1.[154]^38 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.[155]^39 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[156]^40
   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[157]^41 observed that delivery of a
   replication competent HBV genome into hepatocyte lines Huh7 and PMH
   induced the expression of CCNB1. Cheng et al[158]^42 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,[159]^43 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,[160]^44^,[161]^45 miR-582-5p-CDK1,[162]^46 and
   miR-214/490-E2F2,[163]^47^,[164]^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[165]^49
   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[166]^50 and
   miR-125b[167]^51 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.[168]^52 The high expression of LINC00346 was also
   found to be significantly associated with poor OS in HCC[169]^53 and
   breast cancer samples.[170]^54 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 ([171]GSE77509
   and [172]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)[173]^55^–[174]^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
   ([175]http://www.ncbi.nlm.nih.gov/geo/). The mRNA and miRNA Seq-data
   were obtained from TCGA ([176]https://tcga-data.nci.nih.gov/).

Footnotes

   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.

References