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

BACKGROUND: Circular RNAs (circRNAs), a novel type of non-coding RNA, play a vital role in the pathogenesis and development of cancer. CircRNAs signatures may be useful as prognostic and predictive factors as well as clinical tools for evaluating disease status and prognosis. This study was carried out to explore novel circRNA signatures in hepatocellular carcinoma (HCC).
METHODS: The expression profiles of circRNAs were retrieved from the Gene Expression Omnibus (GEO). The expression profiles of miRNAs and mRNAs were obtained from The Cancer Genome Atlas (TCGA) and the Genotype-Tissue Expression (GTEx) database. The results of the microarray were validated by quantitative real-time RCR (qPCR). Based on circRNA-miRNA pairs and miRNA-mRNA pairs, a competitive endogenous RNA (ceRNA) network was constructed. Functional analysis was performed via Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), gene set enrichment analysis (GSEA), and gene set variation analysis (GSVA). Furthermore, survival analysis was carried out using the Kaplan-Meier curve and the log-rank test.
RESULTS: Differentially expressed circRNAs in HCC from GEO databases (GSE94508 and GSE97332) were screened and analyzed using the bioinformatics method. We detected a total of 26 differentially expressed circRNAs by qPCR and then selected 6 circRNAs to construct the circRNA-miRNA-mRNA networks. Through prognostic analysis, 3 target hub genes (AURKA, KIF5B, and RHOA) of circRNAs were discovered. Moreover, GSEA and GSVA were used to reveal the functions of AURKA, KIF5B, and RHOA in HCC.
CONCLUSIONS: We identified three hub genes, and our results suggest that the circHMGCS1/miR-581/AURKA, circHMGCS1/miR-892a/KIF5B, and circTMCO3/miR-577/RHOA axes may play a vital role in HCC progression. 2020 Annals of Translational Medicine. All rights reserved.

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 |log2FC| >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% CO2 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 . 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

Gene	Forward primer (5’-3’)	Reverse primer (5’-3’)
circRNA(hsa_circ_0004913)	TCACTCTGACGGAACTTGACA	GCTCCTTCTGCTTTGGCTGT
circRNA(hsa_circ_0006248)	AAAGACCACATACGCCAACA	TCAGCCCAAATACTCCAAGA
circRNA(hsa_circ_0038929)	CATGCCTTGGGACCTGTTAA	AGCAAAGTTGAGTGCGGAGA
circRNA(hsa_circ_0005428)	AACAGCTCTGGCAAGAAACT	TTTGATGGAGGAAAGTGGAA
circRNA(hsa_circ_0006913)	AACCCAAGAAACTGCTGTCG	GAAATCCCTGATGTGCTCACT
circRNA(hsa_circ_0007248)	TTCTGCGGGCATTGGTGATA	CTCCCTGAGCATTGCATTGGAC
circRNA(hsa_circ_0059342)	GGGACCTCAAATCAAAGAAA	TCCTGAAGTTGTAATGTGGG
circRNA(hsa_circ_0007762)	GGGACCTCAAATCAAAGAAA	GCCCACATAGAGCCACTTAC
circRNA(hsa_circ_0059859)	AGGTGCGGACTCTTGTTCTT	GAGTGCAGATGATGTGGAGG
circRNA(hsa_circ_0056548)	AGGAGGAGAAGTTGGAGGCT	GAGGCAATGGAGTTAAAGGACA
circRNA(hsa_circ_0043302)	TTGACGCCTTAGATAACAGC	ATTCAATCCAGGATGCAAAA
circRNA(hsa_circ_0055033)	AAGAAGTCCTGCATCGAAAT	CACATGATCCTTACTGTCCG
circRNA(hsa_circ_0001359)	CCTGGAATCACGAAGCACAG	TGGGTAATACTGCCGCTGGT
circRNA(hsa_circ_0067323)	ATCCACAACTGATGCACCTA	TGCTTCAACAGTGACTACGC
circRNA(hsa_circ_0001360)	GTTGTCTCGTCGTCATCGTC	GCTTATTAACTGTGCAGGTGT
circRNA(hsa_circ_0000204)	CCTCCTAGAAGTAATGCCCACA	CAGCACCTCACTGTCAAAGAAT
circRNA(hsa_circ_0009581)	CGACTGTGACCTCCTTATGT	ACTGTTCAGCCTCCTTGTCT
circRNA(hsa_circ_0072389)	TGGCTGCCACTCTGTACTCT	TTCCAGCATCTACACCATCA
circRNA(hsa_circ_0007540)	CAGCAGTGGAGGCAGAAGTT	AAACCAAGGGATGGCATAGA
circRNA(hsa_circ_0027774)	TGTGCTGGAATTGATGTTCG	CTGAGGGAGGGTCTGTTTGA
circRNA(hsa_circ_0000517)	TGAGCTTCGGGGAGGTGAGTT	GGGAGAGCCCTGTTAGGGC
circRNA(hsa_circ_0030130)	CTGCCTGTAAGTGCCAAGT	GAGGCCAGGGTAGTCCATT
circRNA(hsa_circ_0084429)	GGCAACTGTATGAGCCACT	GGACTAACTCCCTCGAAAT
circRNA(hsa_circ_0031027)	GGACCCTATTATCGGAAGC	CACCTATGGAGGCGAAGAA
circRNA(hsa_circ_0008160)	GTGCTCTTGTTTGGGATTA	TTTGAGGCCAATAGGTTAA
circRNA(hsa_circ_0011536)	GACAATGCTCAAGTGGTGGTA	TGCTGTCAAGGGTAGGAGTTA
GAPDH	TGACTTCAACAGCGACACCCA	CACCCTGTTGCTGTAGCCAAA

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 (). A total of 869 DEcircRNAs (421 up-regulated and 448 down-regulated) were identified in the GSE97332 dataset (). 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 (), indicating that they were appropriate for further research. Moreover, the basic characteristics of the six DEcircRNAs are listed in . Their basic structural patterns are shown in 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 (|log2[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

circRNA ID	Position	Genomic length	Strand	Best transcript	Gene symbol
hsa_circ_0004913	chr17:62248459-62265775	17,316	−	NM_018469	TEX2
hsa_circ_0072389	chr5:43294157-43299077	4,920	−	NM_001098272	HMGCS1
hsa_circ_0000517	chr14:20811404-20811492	88	−	NR_002312	RPPH1
hsa_circ_0031027	chr13:114164552-114193822	29,270	+	NM_017905	TMCO3
hsa_circ_0008160	chr21:38439561-38441924	2,363	−	NR_028352	PIGP
hsa_circ_0011536	chr1:35824525-35827390	2,865	+	NM_005095	ZMYM4

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 |log2FC| >1 and FDR <0.05 ().

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

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 (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 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 |log2FC| >1 and FDR <0.05, we obtained a total of 2,135 differentially expressed mRNAs (DEGs) from the TCGA-LIHC and GTEx databases (). 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 (). Then, Cytoscape was used to construct the circRNA-miRNA-hub gene networks (). 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 (, 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 (, P<0.001). In addition, elevated expressions of KIF5B and RHOA correlated with the clinical tumor stage (, 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 (). 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 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.