Comprehensive analysis of long non‑coding RNA‑messenger RNA‑microRNA co‑expression network identifies cell cycle‑related lncRNA in hepatocellular carcinoma.

Long non‑coding RNAs (lncRNAs) have been shown to contribute to progression and prognosis of hepatocellular carcinoma (HCC). However, expression profiling and interaction of lncRNAs with messenger RNAs (mRNAs) and microRNAs (miRNAs) remain largely unknown in HCC. The expression profiling of lncRNAs, mRNA and miRNAs was obtained using microarray. The Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes pathway analysis were used to characterize potential functions of differentially expressed mRNAs. Cytoscape was applied to construct an lncRNA‑miRNA‑mRNA co‑expression network and candidate lncRNAs were validated via quantitative PCR in 30 pairs of HCC and adjacent tumor‑free tissues. In this study, 1,056 upregulated and 1,288 downregulated lncRNAs were identified, while 2,687 mRNAs and 6 miRNAs were aberrantly expressed in HCC compared with adjacent tumor‑free tissues. Potential functions of differentially expressed mRNAs were demonstrated to significantly participate in modulating critical genes in the cell cycle, such as cyclin E1 and cyclin B2. After screening, 95 lncRNAs, 5 miRNAs and 36 mRNAs were recruited for construction of lncRNA‑mRNA‑miRNA co‑expression network in the cell cycle pathway. Subsequently, the top 5 lncRNAs that potentially modulate critical genes in the cell cycle were selected as the candidates for further verification. Kaplan‑Meier curves using the Cancer Genome Atlas database showed that 13 targeted mRNAs were associated with overall survival of HCC patients. Finally, three lncRNAs, including ENST00000522221, lnc‑HACE1‑6:1 and lnc‑ICOSLG‑11:1, are significantly upregulated in HCC tissues compared with adjacent tumor‑free tissues. These findings suggest that lncRNAs play essential roles in the pathogenesis of HCC via regulating coding genes and miRNAs, and may be important targets for diagnosis and treatment of this disease.

Introduction

Hepatocellular carcinoma (HCC), the most prevalent subtype of liver cancer, has become the third leading cause of cancer-related death worldwide, accounting for nearly one-third of malignancies (1,2). Hepatitis B and C viral infection, aflatoxin B1 exposure, excessive alcohol consumption, obesity and some inherited metabolic disorders are generally considered to be major risk factors in the development of HCC (3). While current understanding of molecular mechanisms underlying HCC development has significantly improved, knowledge is limited in protein-coding genes that constitute the part of the genome. Recently, high-throughput studies of the human transcriptome have revealed that there is far more genomic transcription than previously anticipated, with the vast majority of the genome being transcribed into an astounding number of non-coding RNAs (ncRNAs) (4-6). A subset of these molecules displays distinct cancer-specific expression patterns, suggesting that they may be potential drivers of cancer biology and have increasing value as clinical biomarkers.

These ncRNAs are classified according to their sequence length into long ncRNA (lncRNA; >200 bp) and small ncRNA (<200 bp), such as microRNA (miRNA/miR), and they do not serve as a template for protein synthesis but nonetheless are potent regulators of gene expression at epigenetic (7), transcriptional (8) and post-transcriptional levels (9). For the past two decades, miRNAs have been recognized as controlling gene expression by binding to specific sites at the 3′ untranslated region of target mRNAs, causing translational repression or degradation (10), while the lesser-studied lncRNAs are implicated in almost every epigenetic regulation event (11). In addition, both miRNAs and lncRNAs have been implicated to modulate diverse cellular pathways that lead to oncogenesis, metastasis and progression (12,13). Despite growing evidence expanding the understanding of the complex role of individual lncRNAs or miRNAs in various steps of cancer development and progression, co-expression profiles, regulatory networks and interaction of lncRNA-mRNA-miRNA remain largely unexplored in cancer.

In the authors' previous study, the microarray profiles of miRNAs in HCC were reported (14). In the present study, high-throughput technologies were employed to simultaneously establish expression profiles of lncRNAs and mRNAs in the same 7 pairs of HCC and adjacent tumor-free specimens. An integrative method was then applied to comprehensively analyze the relationships between aberrantly expressed lncRNAs, mRNAs and miRNAs. The significant involvement of one novel lncRNA in the regulation of cell cycle was highlighted. The comprehensive analysis by the present study provided novel insights into the lncRNA-mRNA-miRNA co-expression network at the transcriptional and post-transcriptional levels in modulating HCC growth. These findings may serve as the foundation for future studies of lncRNA-mRNA-miRNA interactions in HCC and contribute to further systematic studies of tumorigenesis in the future.

Materials and methods

Sample collection

A total of 37 paired HCC specimens and matched adjacent tumor-free tissue counterparts were collected from patients (average age: 56.92 years; male/female =6.4) admitted to the Zhongshan Hospital of Fudan University between January 2008 and December 2009. The diagnosis was based on imaging presentation with a typical pattern of HCC and clinicopathological examination after surgery (Fig. S1). Among these specimens, seven paired HCC specimens and matched adjacent tumor-free tissue counterparts were randomly selected for microarray analysis. The clinical information of the seven patients were provided in Table SI.

The present study was approved by the Institutional Ethics Committee of Zhongshan Hospital of Fudan University and conformed to the provisions of the Declaration of Helsinki. Written informed consent was obtained from all patients who participated in this study.

RNA preparation

Total RNAs were extracted and purified from 7 pairs of specimens for lncRNA, mRNA and miRNA, respectively, using an RNeasy mini kit for total RNA (Qiagen, GmBH) and mirVana™ miRNA Isolation kit without phenol (Ambion; Thermo Fisher Scientific, Inc.). The whole process was conducted according to the manufacturer's protocols. NanoDrop ND-2000 (Thermo Fisher Scientific, Inc.) was applied for quantification of total RNAs and RNA Integrity Number was detected using Agilent Bioanalyzer 2100 (Agilent technologies, Inc.).

Microarray analysis

Two microarrays, including human lncRNA plus mRNA microarray v5.0 (4x180 K) and human miRNA microarray v21.0 (Agilent Technologies, Inc.), were used for detection of differential expression of lncRNAs plus mRNAs, and miRNAs, respectively. The threshold set for differentially expressed genes (DEGs) was P<0.05, fold-change ≥2.0 and false discovery rate (FDR) value <0.05. Vioplot was conducted to explore signal intensity of DEGs. Hierarchical clustering was performed to display the distinguishable gene expression pattern among samples. The data discussed in this publication have been deposited in NCBI's Gene Expression Omnibus and are accessible through GEO Series accession number GSE101728 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE101728) and GSE108724 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE108724).

Gene ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes pathway (KEGG) pathway analysis

GO enrichment analysis was evaluated by Database for Annotation, Visualization and Integrated Discovery bioinformatics tool (https://david.ncifcrf.gov) to display the potential biological functions of DEGs. Enrichment map by Cytoscape Desktop program (http://baderlab.org/Software/EnrichmentMap) was applied to visualize the results of gene-set enrichment as a network in the aspects of biological process, cellular component and molecular function. KEGG pathway was used to identify significantly enriched signaling transduction pathways. The significance of the pathway was measured by P-value. Furthermore, Gene Set Enrichment Analysis (GSEA) for all detected mRNAs was performed in order to validate the consistence of the findings from GO and KEGG pathway enrichment analysis.

Construction of the lncRNA-mRNA-miRNA interaction network

To establish a co-expression network of lncRNAs, mRNAs and miRNAs, differential expressed mRNAs involved in the most significant pathway from the analysis above were recruited. During this procedure, the interactions of miRNA-mRNA were predicted by miRanda database tool (http://www.microrna.org). Furthermore, LncTar (http://www.cuilab.cn/lnctar) was employed to predict the candidate miRNAs and mRNAs that may have underlying relationships with lncRNAs (free energy cutoff <-20). Pearson correlation coefficient was calculated to measure the intensity between lncRNA, mRNA and miRNA. Cytoscape was applied to construct and visualize the potential networks between the three types of RNAs.

Overall survival data collection

To estimate the potential roles of the targeted mRNA in the prognosis of HCC patients, the RNA sequencing data of 365 HCC samples and related clinical data were retrieved from the Cancer Genome Atlas (TCGA) database for later Kaplan-Meier tests.

Quantitative PCR (qPCR)

In the validation step, total RNAs were extracted from 30 pairs of HCC and matched tumor-free specimens by TRIzol® reagent (Thermo Fisher Scientific, Inc.). The reverse transcription and the detection of qualification were conducted by the assistance of the PrimeScript™ RT Master Mix kit and the SYBR® Premix Ex Taq™ kit (Takara Biotechnology Co., Ltd.), respectively. Reverse transcription was conducted at 37°C for 15 min, and then at 85°C for 5 sec for heat inactivation. The conditions for qPCR were as follows: Stage 1: Initial denaturation, 1 cycle, 95°C for 30 sec; stage 2: PCR, 40 cycles, 95°C for 5 sec, 60°C for 30 sec. Data were analyzed using the 2−ΔΔCq method (15). Primers for the lncRNAs were designed and specialized by Primer Premier version 5.0 (Premier Biosoft International; Table SII). The expression level of the candidate lncRNAs was measured by the ABI 7900HT Fast Real-Time PCR system, with GAPDH as the internal control. qPCR was repeated three times.

Statistical analysis

All statistical analysis of microarray data was conducted using the R package, including vioplot (CRAN, v0.2), pcaMethods (Bioconductor, v1.76.0), pheatmap (CRAN, v1.0.8) and limma (Bioconductor, v3.40.0). Fold-change ≥2.0, P<0.05 and FDR <0.05 were set up for further analysis. Expression of the candidate lncRNAs was presented as the scatter plot with mean ± standard deviation and the range of the log-transformed expression levels was estimated by Wilcoxon rank-sum test. Kaplan-Meier estimator was applied to analyze the overall survival data from TCGA database. P<0.05 was considered to indicate a statistically significant difference.

Results

Strategy

To gain insights into possible interactions between lncRNAs, mRNAs and miRNAs in HCC, the present study first sought to create a high-throughput map of such an expression profile for use in subsequent analyses. A five-step approach was applied to identify targets of HCC-specific lncRNAs, mRNAs and miRNAs (Fig. 1). Firstly, expression profiles of lncRNAs, miRNAs and mRNAs in HCC and matched tumor-free specimens were obtained using high-throughput microarray assays. Secondly, GO enrichment, KEGG pathway analysis and GSEA were used to enrich and analyze differentially expressed gene sets. Thirdly, the significant correlation networks between lncRNAs, mRNAs and miRNAs, were explored using online analysis tools (miRanda and LncTar). Fourthly, the differentially expressed mRNAs involved in the most significant pathways were retained to construct the lncRNA-mRNA-miRNA interaction networks. Finally, some of the predicted direct target-related lncRNAs were selected to be validated in thirty paired HCC specimens and matched adjacent tumor-free tissues.

Figure 1

Overview of the study design for lncRNA-mRNA-miRNA co-expression network in HCC. HCC, hepatocellular carcinoma; miRNA, microRNA; lnc, long noncoding RNA; KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, gene ontology; q, quantitative.

Evaluation of microarray raw data

Vioplot was applied for visualization of the distribution of the signal intensity of lncRNA-mRNA microarray from seven paired HCC specimens and matched adjacent tumor-free tissue counterparts. Specifically, lncRNAs from the microarray data were divided into six types, including bidirectional, exonic antisense, exonic sense, intergenic, intronic antisense and intronic sense lncRNAs. The mean expression level of these six types of lncRNAs was slightly lower than that of mRNAs, while the standard error expression between lncRNA and mRNA was almost the same (Fig. S2).

In the lncRNA-mRNA microarray, a total of 91,007 lncRNA and 29,859 mRNA species were detected in the HCC and matched tumor-free specimens. Meanwhile, 2,570 miRNA species were identified in the miRNA microarray. After filtering out the low signal intensity, 32,793 lncRNA, 22,224 mRNA and 369 miRNA species were retained (Table SIII). By the application of PCA, top 1,000 lncRNA and mRNA species that showed the greatest variation and 369 miRNAs were visualized in a heatmap. The findings revealed that the distribution of lncRNAs and mRNAs in HCC group was separated remarkably from the controls, implicating that these abnormally expressed lncRNAs and mRNAs may play an important role in HCC (Fig. S3A-C). However, for the miRNA microarray, the distribution did not present well (Fig. S3D).

Differentially expressed lncRNAs, mRNAs and miRNAs in HCC specimens

With the inclusion criteria of fold-change ≥2.0 and FDR <0.05, 2,344 lncRNA species (1,056 being upregulated and 1,288 being downregulated), 2,687 mRNA (1,171 being upregulated and 1,516 being downregulated) and 6 miRNA (5 being upregulated and 1 being downregulated) were confirmed to be differentially expressed in the present study (Table SIV).

GO enrichment and KEGG pathway analysis

To uncover the involvement of target genes for lncRNAs and miRNAs in various biological processes, GO enrichment and KEGG pathway analysis were used for differentially expressed mRNAs. During this procedure, GO enrichment covers three domains, including molecular functions (MF), biological processes (BP) and cellular components (CC). As shown in Fig. 2A-C and Table SV, predicted target genes of upregulated mRNAs fall into purine nucleotide binding in MF, cell cycle in BP and intracellular non-membrane-bound organelle in CC, such as cell cycle-related genes CDKN2A (cyclin-dependent kinase inhibitor 2A), cyclin-dependent kinase inhibitor 2B and Aurora kinase A in BP. Furthermore, predicted target genes of down-regulated mRNAs involved in carbohydrate binding in MF, immune response in BP and extracellular region in CC, such as immune response-related gene Mannan-binding lectin serine protease 1 and Mannan-binding lectin serine protease 2 in BP.

Figure 2

GO and KEGG pathway analysis for differentially expressed mRNAs from microarray test. GO analysis consisted of three aspects, including (A) molecular functions, (B) cellular components, (C) biological processes and (D) KEGG pathway. The green binding network enriched the upregulated DEGs while the blue binding network gathered the downregulated ones. Every red spot in the separated networks represented each of the involved GO analysis pathways. The number of red spots a network possessed reflects associations with the number of DEGs related to the binding circle. Meanwhile, the thickness of binding lines indicated the binding intensity of the DEGs from the two connected red spots. The major pathways are provided in text as examples of possible molecular interactions for which non-coding RNAs may have regulatory effects. For full spectrum of GO and KEGG pathway analysis information, please refer to Table SV and SVI, respectively. DEG, differential expressed genes; KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, gene ontology.

In the KEGG pathway analysis, upregulated mRNAs were predicted to modulate the cell cycle, such as cyclin-dependent kinase 1 (CDK1), CCNE1, cyclin E2 (CCNE2), proliferating cell nuclear antigen (PCNA) and ubiquitin-conjugating enzyme E2C; while the downregulated mRNAs are associated with metabolic pathways, such as formimidoyltransferase cyclodeaminase, glutamic-pyruvic transaminase, glutamic-oxaloacetic transaminase 1, glutamic-oxaloacetic transaminase 2 and spermidine N1-acetyltransferase 1 (Fig. 2D and Table SVI). All significant pathways for mRNAs that may exert modulating functions are shown in Tables SVII and SVIII. Notably, among these related pathways, target genes in controlling cell cycle were validated as the most significant data set (P=8.75x10-21).

To validate the results of GO enrichment and KEGG pathway analysis, all of the mRNAs detected by microarray were recruited into the GSEA based on fold-change value (16). The findings demonstrated that mRNAs involved in the cell cycle were markedly upregulated, which is consistent with the results of GO enrichment and KEGG pathway analysis. Meanwhile, when the R language package was applied to visualize the GSEA of top three significant pathways in modulating cell cycle, G2-M checkpoint and cell cycle mitotic, are prominent, which further confirms that the genes controlling cell cycle pathways are very active in HCC development and progression (Fig. S4).

Construction of lncRNA-mRNA-miRNA co-expression network

Based on the above network, it was shown that differentially expressed mRNAs involved in the cell cycle pathway played a crucial role in the pathogenesis of HCC and they were then selected as basic elements for further exploration of the coordination and orchestra of lncRNAs, mRNAs and miRNAs in specific cell functions, such as cell cycle control. With the assistance of miRanda and Lnctar, the significantly expressed miRNAs and lncRNAs binding to the cell cycle-related genes were screened out. After the calculation of Pearson correlation coefficients, the screening results showed that a co-expression network consists of 193 nodes and 655 connections between 95 lncRNAs, 36 mRNAs and 5 miRNAs. Within this co-expression network, 187 pairs presented as positive and 468 pairs presented as negative. This co-expression network indicated that one lncRNA could target 10 mRNA at most in a given process and that one mRNA could be under control by three lncRNAs at most; while one miRNA could target 31 coding genes at most, and one mRNA could correlate with three lncRNAs at most in a specific biological activity. This provides a snapshot for gene expression control for a particular gene in normal or malignant processes. This multiple-layer modulation at epigenetic, transcriptional (mRNA), post-translational levels (miRNA or lncRNAs) constitutes a very comprehensive signaling network governing the cell cycle in a fine-tuned fashion. Taken together, a total of 95 lncRNAs, 36 mRNAs and 5 miRNAs are involved in the predicted lncRNA-mRNA-miRNA interaction network after combining the proposed binding targets (Fig. 3). Next, these lncRNAs, miRNAs and mRNAs were categorized into different phases of the cell cycle based on the predicted target molecules (Fig. 4). Among these 95 lncRNAs, the top five lncRNAs associated with cell cycle modulation, including ENST00000522221, ENST00000577319, lnc-GADD45B-4:1, lnc-HACE1-6:1 and lnc-ICOSLG-11:1, were selected for validation analysis, and all of targeted mRNAs, miRNAs and lncRNAs for these five lncRNAs are also shown (Table I). Kaplan-Meier curves of overall survival of targeted mRNAs, including cyclin B1 (CCNB1), CCNB2, CCNE1, CCNE2, cell division cycle 7 (CDC7), CDK1, mitotic arrest deficient 2 like 2, minichromosome maintenance complex component 2 (MCM2), minichromosome maintenance complex component 3 (MCM3), minichromosome maintenance complex component 4 (MCM4), minichromosome maintenance complex component 7 (MCM7), PCNA, polo like kinase 1 (PLK1), are also shown using TCGA database (Fig. 5).

Figure 3

Construction of the lncRNA-mRNA-miRNA co-expression network. All RNAs in this plot were screened out from the genes in association with the cell cycle after validation. The red spots stood for upregulated genes and blue ones for downregulated. The bigger circle (left) is composed of 95 lncRNAs and the smaller one (right) contains 36 mRNAs. One downregulated miRNA (blue) and four upregulated miRNA (red) are shown in the upper part. The red lines are applied for the binding relationship between miRNAs and lncRNAs and mRNA. The green lines show the associations between mRNAs and lncRNAs. miRNA, microRNA; lnc, long noncoding RNA.

Figure 4

Distribution diagram of lncRNAs, mRNAs and miRNAs in the four phases of the cell cycle. At first, 36 mRNAs of cell cycle-related genes were distributed into four visual stages of cell cycle. They were displayed in the relevant cycle parts, indicating possible target labeled genes. A total of five miRNAs are presented in the middle of the circle for that they were predicted to be associated with all phases of cell cycle. The categorization of lncRNAs was based on their predicted relationship with 36 cell cycle-related mRNAs. Similarly, they were labeled as lncRNA including in their own cycle phase. The lncRNA, mRNA and miRNAs that had been cited in previous studies are underlined and labelled with the references. miRNA, microRNA; lnc, long noncoding RNA.

Table I

The network of the interactions between lncRNAs and miRNAs and mRNAs.

Gene ID (lncRNA)	Regulation	Target lncRNA	Target miRNA	Target mRNA
ENST00000522221	Up	-	hsa-miR-106b-5p, hsa-miR-93-5p	CDC7, CCNE1, MCM3, PCNA, MCM2, CCNB2, PLK1, CCNB1, CCNE2
ENST00000577319	Down	-	-	MAD2L2, CDC7, CDK1, MCM3, PCNA, MCM2, MCM4
lnc-GADD45B-4:1	Up	ENST00000502102, lnc-POLR1E-3:2, lnc- WDSUB1-2:2, NR_103821	-	MAD2L2, CCNE1, MCM7, MCM3, CCNB2, MCM4
lnc-HACE1-6:1	Up	ENST00000548595, NR_038131, NR_103821,	hsa-miR-106b-5p, hsa-miR-93-5p	CCNE1, MCM7, MCM2, CCNB2, PLK1, CCNB1
lnc-ICOSLG-11:1	Up	ENST00000502102, lnc-POLR1E3:2, NR_038131	hsa-miR-106b-5p, hsa-miR-93-5p	CDC7, CCNE1, MCM7, MCM3, MCM2, CCNE2

lnc, long non-coding; hsa, homo sapiens; miRNA/miR, microRNA.

Figure 5

Kaplan-Meier survival analysis of targeted mRNAs based on the data from TCGA database. Overall survival rate of targeted mRNAs, including (A) CCNB1, (B) CCNB2, (C) CCNE1, (D) CCNE2, (E) CDC7, (F) CDK1, (G) MAD2L2, (H) MCM2, (I) MCM3, (J) MCM4, (K) MCM7, (L) PCNA and (M) PLK1. Adapted from UALCAN: http://ualcan.path.uab.edu/index.html.

Validation of candidate lncRNAs

To validate the reliability of the results obtained from the microarray, top five differentially expressed lncRNAs was selected for qPCR validation in 30 pairs of HCC and matched tumor-free specimens. The results found that ENST00000522221, lnc-HACE1-6:1 and lnc-ICOSLG-11:1 are significantly upregulated in HCC tissues compared with that in adjacent tumor-free tissues (P=0.0080, P=0.0034 and P=0.0003, respectively; Fig. 6).

Figure 6

Validation of five candidate lncRNAs in 30 pairs of HCC and adjacent tumor-free tissues. Expression levels of top five candidate lncRNAs, including ENST00000522221 (A), ENST00000577319 (B), lnc-GADD45B-4:1 (C), lnc-HACE1-6:1 (D) and lnc-ICOSLG-11:1 (E). Data are presented as scatter plot of mean values with standard deviation. **P<0.01, ***P<0.001. HCC, hepatocellular carcinoma; lnc, long non-coding.

Discussion

Aberrant expression of ncRNAs, including lncRNAs and miRNAs, as well as protein-coding genes have been widely used to uncover underlying molecular mechanisms contributing to tumorigenesis, progression, metastasis and drug resistance or relapse. A comprehensive framework of orchestrated expression and interaction of lncRNAs, mRNAs, miRNAs in HCC, however, remains to be elucidated. In this study, an integrated lncRNA-mRNA-miRNA signature that consists of co-expressed lncRNA, mRNA and miRNA datasets was examined. Furthermore, GO enrichment and KEGG pathway analyses of differentially expressed mRNAs point to potential functions and pathways critical for HCC cell cycle, and a co-expression network that was composed of 193 network nodes and 655 connections between 95 lncRNAs, 36 mRNAs and 5 miRNAs was constructed to reflect possible interactions of these molecules.

Involved mRNAs and ncRNAs in the interaction network, which potentially modulate critical genes in four phases of the cell cycle, including G1, S, G2 and M were listed. The mRNAs that engaged in the four parts of the cycle belong to prominent and classic cell cycle-related genes, such as CDKN2A/B (17), PCNA (18) and CDC7 (19) in G1/S phase and cell division cycle 25C (20), CDK1 (21) and PLK1 (22) in G2/M phase. The miRNAs and lncRNAs were categorized based on their predicted targeting mRNA encoding genes involved in the cell cycle. It turned out that all the five enriched miRNAs (including miR-4270, miR-3651, miR-93-5p, hsa-miR-106b-5p and miR-221-3p) have potential target molecules critical for the entire cell cycle. Notably, no prior studies have proved that these miRNAs are associated with the cell cycle in HCC specimens. In 2012, a study demonstrated that downregulated miR-106b induces G1/S arrest in HCC cell lines (23). Another study showed that overexpression of miR-93-5p promotes G1 or S arrest in ovarian carcinoma (24). Additionally, upregulated miR-221-3p promotes accumulation of human fibroblasts in the G1/S phase (25). It is noteworthy that none of the 95 lncRNAs have been reported to participate in the modulation of the cell cycle in HCC. Only ENST00000557691 and ENST00000609697 have been reported in studies of microarray expression profiles in glioblastoma multiforme (26) and multiple sclerosis (27), respectively. Finally, ENST00000522221, lnc-HACE1-6:1 and lnc-ICOSLG-11:1 were validated as cell cycle-related regulators in an extended 30 pairs of HCC specimens. Therefore, the present study provides a snapshot of how mRNAs, miRNAs and lncRNAs govern a given process, such as the cell cycle, at epigenetic, transcriptional and post-translational levels in an orchestrated fashion in HCC.

Although the study of lncRNAs is still in its infancy, it has become increasing apparent that lncRNAs are critical regulators for cellular and physiological processes by modulating gene and other ncRNA expression. Previous studies have separately demonstrated the involvement of miRNAs (28) and lncRNAs (29) in the development of cancer. Since the first cancer-targeted MRX34, a liposome-based miR-34 mimic, entered Phase I clinical trials in patients with advanced HCC in 2013, understanding the biology of lncRNAs and its connections with mRNAs and miRNAs may lead to promising clinical applications in HCC. Moreover, several studies are available that focus on individual lncRNA expression profiles and regulatory cascades in the tumorigenesis of HCC, such as the HCC upregulated long non-coding RNA (30), H19 (31), metastasis associated lung adenocarcinoma transcript 1 (32) and HCC upregulated EZH2-associated long non-coding RNA (33). However, no studies provide a snapshot of multilayer modulation and interaction networks of lncRNAs, mRNAs and miRNAs in hepatic carcinogenesis, although dissecting each network or interactions involved is far beyond the scope of this manuscript, and requires more advanced approaches and extensive resources.

In the present study lnc-HACE1-6:1, lnc-ICOSLG-11:1 and ENST00000522221 in 30 pairs of HCC specimens were validated. Located on chr6, lnc-HACE1-6:1 (865 bp) was identified as a transcript of nucleophosmin 1 pseudogene 10 (ENST00000398310) from the Ensemble database (http://asia.ensembl.org/index.html). Evidence exists that NPM1 is associated with various cancers, such as prostate cancer and acute myeloid leukemia (34-36). Especially, it was found that overex-pression of NPM1 suppresses p53 via blocking cell cycle-related proteins, such as ARF/MDM2 in colorectal carcinoma (37). Similarly, lnc-ICOSLG-11:1 (395 bp), located in chr21, was identified as a transcript of H2A histone family member Z pseudogene 1 (ENST00000416034) from the Ensemble database. The reports revealed that H2AFZ was overexpressed in HCC patients and promoted the cell growth by affecting cell cycle-related proteins (38). In addition, there is a growing belief that lncRNAs can function as endogenous miRNA sponges as a part of competing endogenous RNA (ceRNA) network. For example, the PTEN pseudogene 1, a lncRNA sharing a high degree of sequence homology with tumor suppressor gene PTEN, acts as a decoy for PTEN-targeting miRNAs (39). In the present study, both lnc-HACE1-6:1 and lnc-ICOSLG-11:1 were predicted to bind with miR-106b-5p and miR-93-5p, and to interact with cell cycle-related genes, such as CCNE1 (40), MCM7 (41) and MCM2 (42). It was postulated that high levels of lnc-HACE1-6:1 and lnc-ICOSLG-11:1 affect cell cycle-related genes by acting as the ceRNA via mediating NPM1 and H2AFZ expression, respectively. ENST00000522221, known as a transcript of THUMP domain containing 3 antisense RNA 1, has never been reported to associate with HCC. In the present study, it was predicted that ENST00000522221 was associated with miRNAs, including miR-106b-5p and miR-93-5p, and cell cycle-related genes. Nevertheless, more studies are warranted to gain better understanding of the role of ENST00000522221, lnc-HACE1-6:1 and lnc-ICOSLG-11:1 in hepatic carcinogenesis.

In conclusion, the present study has laid a framework to understand miRNA-lncRNA-mRNA co-expression profiles, interaction and orchestrated modulation of cellular processes, such as the cell cycle in hepatic carcinogenesis. Probing complicated interplays between lncRNAs and mRNAs/miRNAs holds promise to identify new diagnostic and prognostic markers for HCC and, at least in part, novel targets for antitumor therapies.