Identification and integrated analysis of hepatocellular carcinoma-related circular RNA signature

Introduction

Hepatocellular carcinoma (HCC) is the fifth most common malignancy and the main cause of cancer-related death in the world (1). Owing to the lack of early diagnostic biomarkers with high specificity and sensitivity and the high frequency of tumor metastasis, the overall survival of patients with HCC is very poor (2). Therefore, identifying novel biomarkers and therapeutic targets is of great significance in improving the diagnosis and treatment of HCC and could lead to a better understanding of the pathology of HCC.

Accumulating evidence has shown that the dysregulation of circRNAs may contribute to tumorigenesis and cancer progression in humans (3). CircRNA is a newly discovered non-coding RNA with a special structure, which is characterized by a covalent closed-loop without 5’ caps or 3’ tails (4). In recent years, numerous pieces of evidence have revealed that circRNAs could be involved in gene regulation by acting as miRNA sponges, regulating splicing and transcription, or interacting with RNA-binding proteins (RBPs) (5). These findings suggest the potential regulatory role of circRNAs in biologic development and the pathogenesis and progression of diseases, especially types of cancer (6). Aberrant expression of circRNA also played an important role in the initiation and development of HCC (7). Huang et al., for example, found that the activation of the mTOR signaling pathway by circRNA-100338/miR-141-3p/RHEB axis is closely associated with poor prognosis of hepatitis B-related HCC (8). Another study reported that cicrRNA_101368 could function as a sponge of miR-200a to inhibit miR-200a expression, thus promoting the migration of HCC cells through HMGB1/RAGE signaling of miR-200a downstream (9).

In this study, we analyzed the expression profiles of circRNAs, miRNAs, and mRNAs of HCC downloaded from The Cancer Genome Atlas (TCGA) and the Gene Expression Omnibus (GEO). Differentially expressed that circRNAs (DEcircRNAs) were identified and validated. After predicting the sponging of miRNAs by circRNA and miRNA target genes, we established a circRNA-miRNA-mRNA network and a circRNA-miRNA-hub gene network for HCC. To predict the possible mechanisms and function of circRNAs in HCC, we also performed functional enrichment analyses. These results could improve our understanding of the regulatory mechanisms of circRNAs in the pathogenesis and progression of HCC and provide novel therapeutic targets or biomarkers for HCC.

Methods

Microarray data and RNA sequencing data

The microarray data of two circRNA expression profiles were collected from the GEO database (GSE94508 and GSE97332). RNA-sequencing (RNA-seq) data were downloaded from the TCGA database. miRNA sequencing data included 375 liver cancer tissues and 50 normal tissues, and the mRNA sequencing data included 374 liver cancer tissues and 160 normal tissues. We used the Limma R package to explore the significantly differentially expressed genes between tumorous and normal tissues. The cut-off value was |log₂FC| >1 and FDR <0.05 (FC, fold change; FDR, false discovery rate). Due to the public availability of data on the GEO and TCGA databases, ethical approval or informed consent was not required for this study.

Cell culture and validation of overlapping circRNAs with real-time Quantitative PCR

SMMC-7721 cells were cultured in Dulbecco’s Modified Eagle Medium (DMEM, Gibco) supplemented with 10% fetal bovine serum (FBS, Gibco) at 37 °C in a 5% CO₂ atmosphere in a humidified incubator.

Total RNAs from SMMC-7721 cells were extracted using TRIzol (Invitrogen) according to the manufacturer’s protocol. cDNA was synthesized from total RNA using M-MLV RT (Promega). Then real-time qPCR was applied with SYBR Master Mixture (Takara) on the LightCycler 480 II (Roche) according to standard protocols. The primer sequences for GAPDH acted as the control in these experiments. Primers used for detecting circRNAs were synthesized from GeneChem (Shanghai). Primer sequences were listed in Table 1. The quantity of each mRNA was calculated using the 2^−ΔΔCt method, and data were normalized using GAPDH as the loading control.

Table 1 Primer sequences of circRNAs used for qPCR
Full table

Predictions of target miRNA and miRNA target genes

The Circular RNA Interactome (CircInteractome) (10) and Cancer-Specific CircRNA (CSCD) (11) were used to predict target miRNAs and their miRNA binding sites (MREs). These target miRNAs were further selected by DEmiRNA based on TCGA, and the screened miRNAs were identified as potential target miRNAs of the DEcircRNAs.

Interactions between miRNA and mRNA were predicted based on the TargetScan (12), miRTarBase (13), and miRDB (14) databases. Only the target mRNAs recognized by the three databases were considered as candidate mRNAs. Moreover, intersections with the DEmRNAs were identified to screen out the hub genes targeted by the DEmiRNAs.

Construction of the CeRNA network

The circRNA-miRNA-mRNA regulatory network was constructed by using a combination of circRNA-miRNA pairs and miRNA-mRNA pairs and visualized by Cytoscape 3.6.1.

GO and KEGG pathway enrichment analysis

mRNAs in the ceRNA network were evaluated by GO (15) annotation and KEGG (16) pathway analyses with the clusterProfiler package in R. P value <0.05 served as the cut-off point to assess the functional pathways.

Gene functional analysis

GSEA (17) is a computational method for assessing whether a set of genes defined by a priori show statistically significant, consistent differences between two biological states. GSEA was performed to analyze the enrichment of datasets between high- and low-expression groups of hub genes. A false discovery rate (FDR) <25% and nominal P value <5% were set as the cut-off criteria. In addition, GSVA (18) was used to find the most relevant pathways of hub genes, according to the gene sets files from the KEGG databases.

Results

Identification of differentially expressed circRNAs and miRNAs

Two circRNA microarray datasets (GSE94508 and GSE97332) were downloaded from the GEO database, and the DEcircRNAs were generated by using the limma R package. A total of 270 DEcircRNAs (25 up-regulated and 245 down-regulated) were identified in the GSE94508 dataset (Figure 1A,B). A total of 869 DEcircRNAs (421 up-regulated and 448 down-regulated) were identified in the GSE97332 dataset (Figure 1C,D). We found that 26 overlapping DEcircRNAs were obtained from two datasets, and these were chosen for further validation. The expressions of these DEcircRNAs in SMMC-7721 cells were detected using qPCR. Among them, the mRNA expression abundance of six DEcircRNAs was validated to be up-regulated in SMMC-7721 cells (Figure 2), indicating that they were appropriate for further research. Moreover, the basic characteristics of the six DEcircRNAs are listed in Table 2. Their basic structural patterns are shown in Figure 3 based on the data from CSCD.

Figure 1 Heat maps and volcano plots for differentially expressed circRNAs in LIHC based on GSE94508 and GSE97332 from the GEO database. Red represents up-regulated circRNAs and green represents down-regulated circRNAs. Black represents circRNAs without statistically significant change (|log₂[fold change]| >1 and FDR <0.05).

Figure 2 Differentially expressed circRNAs in the two microarray datasets. (A) Heat map of overlapping differentially expressed circRNAs in LIHC based on the GSE94508 and GSE97332 from the GEO database; (B) histogram of mRNA expression abundance of the six target circRNAs in SMMC-7721 cell lines by real-time quantitative PCR (Y-axis is ΔCt).

Table 2 Basic characters of the six hub circRNAs
Full table

Figure 3 Structural graphs of the six circRNAs. Red spots represent miRNA response elements, blue spots represent RNA binding protein, and green curves represent the open reading frame.

Based on the TCGA-LIHC database, we obtained 251 differentially expressed miRNAs (DEmiRNAs, 229 up-regulated and 22 down-regulated) using edgeR package analysis with |log₂FC| >1 and FDR <0.05 (Figure S1).

Figure S1 Heat maps and volcano plots of the differentially expressed miRNAs between LIHC tissues and normal controls in the TCGA dataset.

The construction of the ceRNA network and functional assessment

To better comprehend the role of circRNA in miRNA mediated mRNA, we constructed a circRNA-miRNA-mRNA (ceRNA) network. Firstly, using the CircInteractome database, we predicted the miRNAs targeted by the 6 DEcircRNAs and identified 100 circRNAs-miRNAs pairs. After intersecting with the DEmiRNAs, only 13 circRNA-miRNA pairs remained, including 6 circRNAs and 11 DEmiRNAs. We then identified mRNAs targeted by the 11 DEmiRNAs in 3 databases (miRDB, miRTarBase, and TargetScan). Lastly, we established a ceRNA network based on 6 circRNA nodes, 11 miRNA nodes, and 114 mRNA nodes in HCC (Figure 4).

Figure 4 The ceRNA network of circRNA-miRNA-mRNA in LIHC. Red represents circRNAs, green represents differentially expressed miRNA, and sky blue represents mRNA.

To achieve a better understanding of the underlying biological processes (BP) and pathways associated with the circRNAs-related target genes in the ceRNA network, we applied the GO and KEGG pathway enrichment analyses. The results of BP, cellular component (CC), molecular function (MF), and KEGG pathways are indicated in Figure 5A,B,C,D (P<0.05). In addition, GSVA further confirmed that a high level of circRNAs-related target genes was significantly enriched in 9 KEGG terms (thyroid cancer, pathogenic Escherichia coli infection, hypertrophic cardiomyopathy hcm, lysosome, hedgehog signaling pathway, oocyte meiosis, natural killer cell-mediated cytotoxicity, base excision repair, and cell cycle) (Figure 5E).

Figure 5 Pathway enrichment analyses of circRNA-related target genes. Functional and signaling pathway analysis of these target genes in LIHC according to biological processes (A), cellular components (B), molecular function (C), and KEGG pathway (D). (E) GSVA was used to analyze these target genes in TCGA.

The construction of the CircRNA-miRNA-hub gene network and functional enrichment analysis of hub genes

Based on the thresholds of |log₂FC| >1 and FDR <0.05, we obtained a total of 2,135 differentially expressed mRNAs (DEGs) from the TCGA-LIHC and GTEx databases (Figure S2). The above results showed that 7 DEGs were involved in the ceRNA network. Then, to further explore the functions of these hub genes in the network in carcinogenesis and the development of HCC, we also established a CircRNA-miRNA-hub gene regulatory network. Firstly, we presented 106 circRNAs-related target genes and 2,135 DEGs in TCGA and GTEx and generated 7 common target genes DEGs (CITED2, ACSL4, MARCKS, KIF5B, AURKA, SMO, and RHOA) using Venn diagram analysis (Figure 6A). Then, Cytoscape was used to construct the circRNA-miRNA-hub gene networks (Figure 6B). Moreover, the Kaplan-Meier curves, which were plotted using the GEPIA database, indicated that higher expression of AURKA, KIF5B, and RHOA were significantly correlated with poorer overall survival (Figure 6C,D,E, P<0.05), however, no significant prognostic significance was shown in CITED2, ACSL4, MARCKS, or SMO. We further verified the expression levels of 7 hub genes using EdgeR software in TCGA and GTEx data and found that an up-regulation of AURKA, KIF5B, and RHOA in tumor tissues compared with normal tissues (Figure 6F,G,H, P<0.001). In addition, elevated expressions of KIF5B and RHOA correlated with the clinical tumor stage (Figure S3, P<0.05), suggesting LIHCs with high KIF5B and RHOA expression are prone to progress to a more advanced stage.

Figure S2 Heat maps and volcano plots of the DEGs between LIHC tissues and normal controls in the TCGA and GTEx datasets.

Figure 6 The identification of hub genes in LIHC. (A) VennDiagram analysis of circRNA-related target genes and differentially expressed genes in TCGA and GTEx; (B) the CircRNA-miRNA-hubgene network. The network consists of 5 circRNAs, 6 miRNAs, and 7 hub genes. Red represents circRNAs, blue represents miRNAs, and yellow represents differentially expressed mRNAs; (C-E) Kaplan-Meier curves for survival analyses of three hub-genes in patients with TCGA-LIHC. (C) AURKA (P=0.00022); (D) KIF5B (P=0.047); (E) RHOA (P=0.046); (F-H) the box plot showing the mRNA expression levels of hub-genes detected in the TCGA-LIHC database (P<0.05). (F) AURKA; (G) KIF5B; (H) RHOA.

Figure S3 The association between the expression of hub genes and clinical stage (P<0.05). (A) AURKA; (B) KIF5B; (C) RHOA.

To further study the potential roles of AURKA, KIF5B, and RHOA in HCC, we performed an analysis of GSEA and GSVA on the TCGA-LIHC database. The results of GSEA showed that the hub genes with the highest enrichment score were closely associated with signaling pathways related to cancer (Figure 7A,B,C). Furthermore, GSVA confirmed that many cell cycle-related KEGG pathways, including DNA replication and mismatch repair, were enriched in the groups with high-expression of AURKA, KIF5B, and RHOA, further suggesting that the activation of these hallmark genes might participate in the process of tumor proliferation (Figure 7D,E,F).

Figure 7 The functional analysis of hub genes in the TCGA-LIHC dataset. (A-C) GSEA applied to validate the hub gene signatures. (A) AURKA; (B) KIF5B; (C) RHOA; (D-F) GSVA-derived clustering heat maps of differentially expressed pathways for single hub genes. (D) AURKA; (E) KIF5B; (F) RHOA.

Discussion

Over recent decades, studies have emerged, showing the abnormal expression profiles of non-coding RNAs in HCC (19,20). Some investigations have demonstrated circRNA to be closely associated with the occurrence of HCC and its dysfunction to inhibit tumor growth and metastasis (21,22). However, the molecular mechanisms underlying the participation of circRNAs in HCC progression remain elusive.

To explore the molecular regulatory mechanisms of HCC, we studied circRNA microarray profiles and identified 26 DEcircRNAs in HCC tissues and paired normal tissues. Then, qPCR was applied to detect the mRNA expression abundance of the 26 DEcircRNAs in SMMC-7721 cell lines. Finally, hsa_circ_0004913, hsa_circ_0072389, hsa_circ_0000517, hsa_circ_0031027), hsa_circ_0008160, and hsa_circ_0011536 were confirmed by qPCR to be significantly dysregulated in HCC cells. In addition, studies have indicated that circRNAs could act as promising biomarkers for cancer diagnosis (23). Therefore, our findings indicated that circTEX2 (hsa_circ_0004913), circHMGCS1 (hsa_circ_0072389), circRPPH1 (hsa_circ_0000517), circTMCO3 (hsa_circ_0031027), circPIGP (hsa_circ_0008160), and circZMYM4 (hsa_circ_0011536) could represent potentially valuable diagnostic biomarkers for HCC treatment. The previous study demonstrated that circRNAs had been deemed to possess the function of miRNA sponges, which can suppress the binding of miRNAs to target genes, regulating their expression (24). Therefore, we proposed the hypothesis that the circRNA-miRNA-mRNA axis might be the possible molecular regulatory mechanism underlying HCC. In this study, we constructed a ceRNA network based on the above 6 circRNAs and identified 11 DEmiRNAs and 114 mRNAs in this network. Through GO and KEGG pathway analysis, these mRNAs were suggested to participate in many crucial cancer-related biological functions and pathways.

To further identify the crucial circRNAs involving in the network, we screened 7 hub genes, and thereby established the circRNA-miRNA-hub gene network. Among this network, several miRNAs have been identified to play vital roles in the development and progression of HCC, such as miR-577 (25), miR-326 (26), miR-182 (27), miR-892a (28), and miR-581 (29). Moreover, using the TCGA and GTEx databases, we detected the expression of 7 hub genes and found that the expression levels of AURKA, KIF5B, and RHOA were elevated in tumor tissues compared with normal tissues and that higher expression of these genes was dramatically associated with poorer overall survival. Previous studies have shown that circRNAs could regulate host gene expression or work as a miRNA sponge, regulating miRNA-mediated gene expression and thus leading to tumor progression and development (30). In this study, we investigated DEcircRNAs in HCC and then selected the DEcircRNAs that had host genes that were differentially expressed in TCGA-LIHC data and correlated with the prognosis as a candidate circRNAs. Following this, we speculated that these candidate DEcircRNAs (circHMGCS1 and circTMCO3) might participate in the pathogenesis of HCC. In addition, GSEA and GSVA analyses revealed that the functions of AURKA, KIF5B, and RHOA were also closely associated with cancer progression, which indicated that the circRNAs were meaningfully associated with HCC metastasis.

Conclusions

In conclusion, our research revealed circHMGCS1/miR-581/AURKA, circHMGCS1/miR-892a/KIF5B, and circTMCO3/miR-577/RHOA regulatory axes, which could serve as novel biomarkers in the detection of HCC. Finally, bioinformatics analysis indicated that the interaction pairs of the circRNAs mentioned above might be involved in the initiation and development of HCC.

Acknowledgments

Funding: This study was supported by the clinical medical science and technology of Jiangsu Technical Department (BL2014060) and the “Medical Leading Talent and Innovation Team” project (LJ201134) of Jiangsu Province.

Footnote

Conflicts of Interest: The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Bray F, Ferlay J, Soerjomataram I, et al. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2018;68:394-424. [Crossref] [PubMed]
2. Siegel RL, Miller KD, Jemal A. Cancer statistics, 2016. CA Cancer J Clin 2016;66:7-30. [Crossref] [PubMed]
3. Memczak S, Jens M, Elefsinioti A, et al. Circular RNAs are a large class of animal RNAs with regulatory potency. Nature 2013;495:333-8. [Crossref] [PubMed]
4. Lasda E, Parker R. Circular RNAs: diversity of form and function. RNA 2014;20:1829-42. [Crossref] [PubMed]
5. Li X, Yang L, Chen LL. The Biogenesis, Functions, and Challenges of Circular RNAs. Mol Cell 2018;71:428-42. [Crossref] [PubMed]
6. Yin Y, Long J, He Q, et al. Emerging roles of circRNA in formation and progression of cancer. J Cancer 2019;10:5015-21. [Crossref] [PubMed]
7. Yu J, Xu QG, Wang ZG, et al. Circular RNA cSMARCA5 inhibits growth and metastasis in hepatocellular carcinoma. J Hepatol 2018;68:1214-27. [Crossref] [PubMed]
8. Huang XY, Huang ZL, Zhang PB, et al. CircRNA-100338 Is Associated With mTOR Signaling Pathway and Poor Prognosis in Hepatocellular Carcinoma. Front Oncol 2019;9:392. [Crossref] [PubMed]
9. Li S, Gu H, Huang Y, et al. Circular RNA 101368/miR-200a axis modulates the migration of hepatocellular carcinoma through HMGB1/RAGE signaling. Cell Cycle 2018;17:2349-59. [Crossref] [PubMed]
10. Dudekula DB, Panda AC, Grammatikakis I, et al. CircInteractome: A web tool for exploring circular RNAs and their interacting proteins and microRNAs. RNA Biol 2016;13:34-42. [Crossref] [PubMed]
11. Xia S, Feng J, Chen K, et al. CSCD: a database for cancer-specific circular RNAs. Nucleic Acids Res 2018;46:D925-9. [Crossref] [PubMed]
12. Tan LP, Seinen E, Duns G, et al. A high throughput experimental approach to identify miRNA targets in human cells. Nucleic Acids Res 2009;37:e137. [Crossref] [PubMed]
13. Chou CH, Shrestha S, Yang CD, et al. miRTarBase update 2018: a resource for experimentally validated microRNA-target interactions. Nucleic Acids Res 2018;46:D296-302. [Crossref] [PubMed]
14. Wong N, Wang X. miRDB: an online resource for microRNA target prediction and functional annotations. Nucleic Acids Res 2015;43:D146-52. [Crossref] [PubMed]
15. Ashburner M, Ball CA, Blake JA, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet 2000;25:25-9. [Crossref] [PubMed]
16. Kanehisa M, Furumichi M, Tanabe M, et al. KEGG: new perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res 2017;45:D353-61. [Crossref] [PubMed]
17. Reimand J, Isserlin R, Voisin V, et al. Pathway enrichment analysis and visualization of omics data using g:Profiler, GSEA, Cytoscape and EnrichmentMap. Nat Protoc 2019;14:482-517. [Crossref] [PubMed]
18. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 2013;14:7. [Crossref] [PubMed]
19. Klingenberg M, Matsuda A, Diederichs S, et al. Non-coding RNA in hepatocellular carcinoma: Mechanisms, biomarkers and therapeutic targets. J Hepatol 2017;67:603-18. [Crossref] [PubMed]
20. O'Brien A, Zhou T, Tan C, Alpini G, Glaser S. Role of Non-Coding RNAs in the Progression of Liver Cancer: Evidence from Experimental Models. Cancers (Basel) 2019. [Crossref] [PubMed]
21. Han D, Li J, Wang H, et al. Circular RNA circMTO1 acts as the sponge of microRNA-9 to suppress hepatocellular carcinoma progression. Hepatology 2017;66:1151-64. [Crossref] [PubMed]
22. Chen Q, Chen Z, Cao S, et al. Role of CircRNAs_100395 in Proliferation and Metastases of Liver Cancer. Med Sci Monit 2019;25:6181-92. [Crossref] [PubMed]
23. Bach DH, Lee SK, Sood AK. Circular RNAs in Cancer. Mol Ther Nucleic Acids 2019;16:118-29. [Crossref] [PubMed]
24. Xiong DD, Dang YW, Lin P, et al. A circRNA-miRNA-mRNA network identification for exploring underlying pathogenesis and therapy strategy of hepatocellular carcinoma. J Transl Med 2018;16:220. [Crossref] [PubMed]
25. Wang LY, Li B, Jiang HH, et al. Inhibition effect of miR-577 on hepatocellular carcinoma cell growth via targeting β-catenin. Asian Pac J Trop Med 2015;8:923-9. [Crossref] [PubMed]
26. Hu S, Ran Y, Chen W, et al. MicroRNA-326 inhibits cell proliferation and invasion, activating apoptosis in hepatocellular carcinoma by directly targeting LIM and SH3 protein 1. Oncol Rep 2017;38:1569-78. [Crossref] [PubMed]
27. Cao MQ, You AB, Zhu XD, et al. miR-182-5p promotes hepatocellular carcinoma progression by repressing FOXO3a. J Hematol Oncol 2018;11:12. [Crossref] [PubMed]
28. Jia B, Tan L, Jin Z, et al. MiR-892a Promotes Hepatocellular Carcinoma Cells Proliferation and Invasion Through Targeting CD226. J Cell Biochem 2017;118:1489-96. [Crossref] [PubMed]
29. Wang YQ, Ren YF, Song YJ, et al. MicroRNA-581 promotes hepatitis B virus surface antigen expression by targeting Dicer and EDEM1. Carcinogenesis 2014;35:2127-33. [Crossref] [PubMed]
30. Kristensen LS, Hansen TB, Venø MT, et al. Circular RNAs in cancer: opportunities and challenges in the field. Oncogene 2018;37:555-65. [Crossref] [PubMed]