Identification of the hsa_circ_0039466/miR-96-5p/FOXO1 regulatory network in hepatocellular carcinoma by whole-transcriptome analysis.

Background: Circular RNAs (circRNAs) are important for the process of cancer initiation and progression. However, the role of circRNAs in hepatocellular carcinoma (HCC) remains incompletely understood. Therefore, we further explored the expression network of circRNAs in HCC.
Methods: Whole-transcriptome microarrays of HCC and paired normal liver tissues were obtained from the Gene Expression Omnibus (GEO) database. The structures of tumor-associated circRNAs were acquired by the Cancer-Specific CircRNA Database (CSCD). StarBase, circBank, and R packages (miRNAtap and multiMiR) were used to predict miRNA targets of circRNAs and downstream molecules of miRNAs. Expression relationships between RNA-RNA interactions were evaluated by data from The Cancer Genome Atlas (TCGA) and GEO databases. ClusterProfiler and DOSE R packages were used for pathway enrichment to explore the biological functions of potential target genes. Finally, a possible circRNA-miRNA-mRNA regulatory network was established based on the competing endogenous RNA (ceRNA) hypothesis.
Results: The differentially expressed circRNAs (DECs) were matched with cancer-specific circRNAs in the CSCD database and a screening analysis was performed to obtain 5 cancer-specific circRNAs. A total of 329 possible target miRNAs for 5 cancer-specific circRNAs were predicted by the circBank database, and intersection analysis with differentially expressed miRNAs (DEmiRNAs) revealed that miR-6746-3p and miR-96-5p were the two most suitable miRNAs targets for our selected circRNAs. Further expression verification and prediction of base complementary paired binding sites demonstrated the hsa_circ_0039466/miR-96-5p axis as a crucial pathway in HCC. Next, we found that FOXO1 and LEPR were two potential downstream molecules of the hsa_circ_0039466/miR-96-5p axis through target gene prediction analysis, differential expression analysis, and intersection analysis. The pathway enrichment results suggested that the hsa_circ_0039466/miR-96-5p axis affects the progression and outcome of HCC through the insulin resistance pathway. Finally, through multi-data crossover analysis and data analysis of HCC samples further confirmed the existence of the hsa_circ_0039466/miR-96-5p/FOXO1 ceRNA regulatory network and that the axis was closely related to clinical stage. Conclusions: hsa_circ_0039466 facilitates the expression of FOXO1 by sponging miR-96-5p, and ultimately inhibits tumor progression. These results provide a theoretical basis for further understanding of the gene expression network of HCC. 2022 Annals of Translational Medicine. All rights reserved.

Introduction

Hepatocellular carcinoma (HCC) is the sixth most frequent malignancy and the third main cause of cancer-related death globally (1-3). Given the quiescent nature of HCC in the early stages, many patients may already be in advanced stages of the disease at the time of diagnosis, thus, limiting treatment options (4). Lack of early detection and prompt treatment is the main reason behind poor outcomes in patients with HCC (5,6). Thus, screening and identifying cancer biomarkers and understanding certain mechanisms of HCC are essential to facilitate early detection and guide treatment.

Although previous researchers demonstrated that several serum tumor markers including α-fetoprotein (AFP), des-carboxy prothrombin (DCP), and Lens culinaris agglutinin A-reactive fraction of alpha-fetoprotein (AFP-L3) might be used as prognostic variables for HCC, these biomarkers have low accuracy predicting HCC prognosis (7-10). It is therefore of great importance to screen and identify accurate biomarkers or clarify certain mechanisms of HCC to facilitate early detection and to guide clinical treatments for HCC.

Circular RNAs (circRNAs), produced by the alternative splicing of premature RNAs, are a class of conservative, stable, and abundant non-coding RNAs (ncRNAs) with tissue- or developmental-specific expression patterns (11-14). Numerous studies have found that a large number of circRNAs exert biological functions through sponge miRNAs. Specifically, they act as competitive endogenous RNAs (ceRNAs) by regulating the expression patterns of miRNA target genes (15). In other words, circRNAs with miRNA response elements can influence the stability or translation of target mRNAs by competitively binding with miRNAs (16). CircRNAs, as a kind of ceRNA, are more efficient than other ceRNAs in miRNA binding due to their high abundance and stability (17). More and more evidence indicates that circRNA/miRNA/mRNA expression networks widely participate in the initiation and progression of several malignancies (13). In recent years, a large number of studies have found that circRNAs can regulate the expression patterns of oncogenes/suppressor genes through ceRNA mechanisms (18), leading to development or progression of cancer. However, in this biological process, the stoichiometric relationship between the miRNA binding site of the circRNA and the mRNA target site of the miRNA must be considered (19,20). At the same time, with the application of high-throughput RNA sequencing technology and innovative bioinformatics algorithms, bioinformatics prediction is widely used to clarify the relationship between circRNAs and cancer by defining their function in the disease process (18). Numerous research projects have investigated the prediction of circRNAs on miRNA sponges and the establishment of circRNA-miRNA-gene regulatory networks, and demonstrated that circRNAs are an important class of transcripts in multicellular organisms (21,22). It was recently reported that the circRNA circUHRF1 functions as a sponge of miR-449c-5p in HCC (23). The circRNA cSMARCA5 inhibited the development of HCC by sponging miR-17-3p and miR-181b-5p to increase the expression of TIMP3 (24). Therefore, we believe that if a suitable circRNA-miRNA-mRNA regulatory network can be established, it may provide key information to accurately predict the prognosis of HCC and discover effective therapeutic targets.

The present study aims to provide a theoretical basis for in-depth exploration of the occurrence and development mechanisms of HCC, and potential diagnostic markers and therapeutic targets. We present the following article in accordance with the MDAR reporting checklist (available at https://atm.amegroups.com/article/view/10.21037/atm-22-3147/rc).

Methods

Tissue samples

A total number of 124 pairs of HCC tissues and adjacent normal tissue specimens used in this study were obtained from patients undergoing surgical resection at The Third Affiliated Hospital of Sun Yat-sen University from December 2012 to September 2018. All samples were placed and stored in liquid nitrogen for long periods.

Research on gene chip data

We collected the RNA-seq data of HCC patients from The Cancer Genome Atlas (TCGA, http://gdc.cancer.gov/) and Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo) databases. Then, we analyzed the data to verify the differentially expressed genes (DEGs) in HCC patients. We used the GSE128274 dataset of the GEO chip for further analysis, which includes expression data of circRNAs, miRNAs, and mRNAs. Survival analysis images [Kaplan-Meier (KM) curves] were generated using the website tool Gene Expression Profiling Interactive Analysis (GEPIA, http://gepia.cancer-pku.cn/). Immunohistochemical (IHC) images were acquired from The Human Protein Atlas database (https://www.proteinatlas.org/).

Differential analysis of circRNAs, miRNAs, and mRNAs

Differential analysis of GSE128274 was performed using the Limma software package. The |log2fold change (FC)|>1 were used to identify differentially expressed circRNAs (DECs), differentially expressed miRNAs (DEmiRNAs), and DEGs in HCC tissues and paired normal tissues. The parental genes of potential circRNAs were identified using circBase. Cancer-Specific CircRNA Database (CSCD, http://gb.whu.edu.cn/CSCD/), an online tool for the study of cancer-specific circRNAs, was used to acquire the structures of potential circRNAs.

Prediction of target miRNAs for circRNAs

The target miRNAs of the candidate circRNAs were predicted by using the circBank database (http://www.circbank.cn/). CircBank is a comprehensive human circRNA database containing more than 140,000 annotated human circRNAs from different source (25). The information on potential miRNAs of circRNAs can easily be obtained by entering the circBase ID of the circRNAs into the search bar of the circBank database.

Identification and validation of miRNA expression

We considered miRNAs included in the target miRNA set and the target DEmiRNA set as potential miRNA targets for specific circRNAs in HCC. Furthermore, TCGA was used to verify the expression of these miRNAs.

Prediction and validation of target mRNAs for miRNAs

The target mRNAs of the candidate miRNAs were predicted using R packages miRNAtap and multiMiR. After merging the target mRNAs predicted by these two tools, target mRNAs potentially involved in regulating HCC progression by miRNAs were obtained by intersection of the DEGs and target mRNAs. In other words, the genes contained in both target genomes and target DEGs were considered potential target genes for specific miRNAs in HCC. Similarly, the TCGA dataset was used to confirm the expression of these target genes.

Pathway enrichment analysis

In order to explore the biological pathways of potential target genes, pathway enrichment analysis of target genes was carried out through the ClusterProfiler and DOSE R language packages. The biological pathways associated with the potential target genes were identified using the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. In addition, the biological processes, molecular functions, and cellular components that changed in potential target genes were analyzed by Gene Ontology (GO) enrichment analysis. Subsequently, the top 10 enrichment pathways after enrichment analysis were selected for correlation analysis. Through correlation analysis, we found the significant target genes (FOXO1 and LEPR) of the hsa_circ_0039466/miR-96-5p axis in HCC. After comparing the different characteristics of FOXO1 and LEPR in terms of prognostic value and gene expression level, we finally chose FOXO1 as the core target of the hsa_circ_0039466/miR-96-5p axis in HCC.

Quantitative real‑time polymerase chain reaction (qRT-PCR)

In this study, total RNA was isolated by homogenizing the obtained clinical HCC samples using Trizol reagent (Invitrogen, USA). Then, the extracted RNA was reverse transcribed into cDNA by the HisScript® III RT SuperMix for qPCR (+gDNA wiper) kit (Vazyme, China). The qPCR experiments were performed according to the instructions of the ChamQ Universal SYBR qPCR Master Mix kit (Vazyme, China). MiRNA reverse transcription and qPCR experiments were performed using the Bulge-Loop miRNA qRT-PCR Starter kit (RiboBio, China). Roche Light-Cycler®480 was used to perform qPCR. GAPDH and U6 were used as endogenous controls for mRNA and miRNA, respectively. Primers used in this experiment were designed and purchased from GENEray (Shanghai, China) (Table S1).

Ethical aspects

The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). This study was approved by the Clinical Research Ethics Committee of The Third Affiliated Hospital of Sun Yat-sen University (No. 2022-02-135), and informed consent was obtained from all participants.

Statistical analysis

R software version 4.10 (http://www.r-project.org) and GraphPad Prism 8.0 were used for statistical analysis. Student’s t-test or Fisher’s exact test was used to investigate differences between groups. The Limma R package was used to analyze RNA-seq data. Pathway enrichment analysis was performed by using the ClusterProfiler and DOSE packages. Data were presented as the mean values ± the standard deviation (SD). A t-test was used to compare means. Spearman correlation analysis was performed among hsa_circ_0039466, miR-96-5p, and FOXO1. P<0.05 was considered significant for all analyses.

Results

Selection of cancer-specific circRNAs in HCC

The circRNA dataset of HCC patients from the GEO database (GSE128274) was chosen for further genetic analysis. A total of 82 significant DECs, including 47 downregulated DECs and 35 upregulated DECs, were identified ( and Table S2). displays the 10 molecules with the greatest differences in expression. Subsequently, the DECs were matched with cancer-specific circRNAs in the CSCD database to obtain 6 genomic positions containing 8 cancer-specific circRNAs ( and ). The names of circBase and parental genes of the 8 circRNAs are shown in .

Figure 1

Identification of potential circRNAs in HCC. (A) Volcano plot of DECs in HCC from the GSE128274 dataset. The dark orange dots and dark green dots represent upregulated DECs and downregulated DECs with significance (adjusted P<0.05 and |log2FC| >1), respectively. The grey dots represent DECs without significance. The maroon dots represent 6 cancer-specific circRNAs with significance. (B) The heatmap of DECs in HCC from the GSE128274 dataset. (C) The intersection analysis of DECs and cancer-specific circRNAs in the CSCD. HCC, hepatocellular carcinoma; DECs, differentially expressed circRNAs; CSCD, Cancer-Specific CircRNA Database.

Table 1

The top 10 DECs between HCC tissues and adjacent normal tissues

Genomic position	Log2FC	AveExpr	t	P	Adj.P	B
chr7:99306657-99358604	−5.165	−0.48	−4.0128	<0.001	0.015	−0.056
chr12:21015689-21175930	−4.276	0.359	−4.27	<0.001	0.009	0.783
chr5:38523520-38530768	−3.649	−1.033	−5.587	<0.001	<0.001	4.514
chr19:10183599-10184111	−3.522	0.469	−4.719	<0.001	0.005	2.102
chr4:100231921-100257933	−3.375	2.522	−4.35	<0.001	0.008	1.094
chr6:3155858-3177876	4.091	0.671	5.969	<0.001	<0.001	5.577
chr6:3155858-3226045	4.234	0.647	6.301	<0.001	<0.001	6.498
chr7:134215394-134250298	5.56	2.843	5.405	<0.001	0.002	3.563
chr7:134222331-134260679	8.163	2.322	7.566	<0.001	<0.001	8.593
chr7:134217755-134254281	8.166	2.25	9.247	<0.001	<0.001	13.164

DECs, differentially expressed circRNAs; HCC, hepatocellular carcinoma.

Table 2

Information of circRNAs in

Genomic position	circBase ID	Log2FC	P	Adj.P
chr16:56660791|56678088	hsa_circ_0039466	−2.76	<0.001	0.008
chr17:40765866|40817580	hsa_circ_0043863/hsa_circ_0043864	3.22	<0.001	0.004
chr17:40766513|40818257	hsa_circ_0000767/hsa_circ_0000768	2.64	<0.001	0.009
chr19:10183599|10184111	novel	−3.52	<0.001	0.004
chr7:134215394|134250298	hsa_circ_0082476	5.56	<0.001	0.001
chr7:134217755|134254281	hsa_circ_0082477	8.17	<0.001	<0.001

circRNAs, circular RNAs.

Prediction and analysis of miRNA targets of potential circRNAs in HCC

We used the circBank database to predict the potential target miRNAs of the remaining 7 circRNAs because the novel circRNA derived from chr19:10183599|10184111 has not been annotated and named in circBase. We found that in the circBank database, 5 of the 7 circRNAs had 329 miRNA targets, and the other 2 circRNAs (hsa_circ_0000767 and hsa_circ_0000768) had no predicted miRNA targets (available online: https://cdn.amegroups.cn/static/public/atm-22-3147-1.pdf). Cytoscape software was adopted to establish a circRNA-miRNA connecting network (). In order to identify the potential target miRNAs of the 5 circRNAs, we screened out DEmiRNAs between HCC tissues and adjacent normal tissues using the dataset GSE128274 from HCC patients. A total of 55 significant DEmiRNAs (Table S3), containing 15 upregulated DEmiRNAs and 40 downregulated DEmiRNAs, were identified. The 55 significant DEmiRNAs are presented in Table S3. We intersected the predicted miRNAs targets of circRNAs and significant DEmiRNAs to find the core target miRNAs of the DECs. As previous studies have shown that ceRNAs are negatively related to miRNAs, we chose the principle of “low expression of circular RNA, high expression of miRNA, and high expression of circular RNA, low miRNA expression” to determine the target miRNAs (13,16). According to this principle, the results showed that miR-6746-3p was included in the target miRNAs of both the upregulated DEC set and the downregulated DEmiRNA set (). Additionally, miR-96-5p was included in the target miRNAs of both the downregulated DEC set and the upregulated DEmiRNA set (). The base complementary paired binding sites between hsa_circ_0039466 and miR-96-5p and between hsa_circ_0043863/ hsa_circ_0043864 and miR-6746-3p are shown in and , respectively. Furthermore, the expression levels of miR-96-5p and miR-6746-3p were confirmed using the expression profile of HCC from TCGA, which demonstrated that only miR-96-5p was significantly upregulated in HCC, in contrast, miR-6746-3p expression was not significantly different in HCC tissues and adjacent normal tissues (). The high expression trend of miR-96-5p in HCC shown in is consistent with that shown in . Based on our analysis, hsa_circ_0039466/miR-96-5p might be an important axis in the carcinogenesis of HCC.

Figure 2

The circRNA-miRNA network of 5 cancer-specific circRNAs. The network consisting of 5 circRNAs (hsa_circ_0039466, hsa_circ_0082476, hsa_circ_0082477, hsa_circ_0043863, and hsa_circ_0043864) and 329 miRNAs was generated by Cytoscape 3.7.1 software.

Figure 3

Identification of potential miRNA targets of selected circRNAs in HCC. (A) The intersection analysis of target miRNAs of circRNAs and DEmiRNAs; (B) the intersection analysis of miRNA targets of upregulated circRNAs and downregulated DEmiRNAs; (C) the intersection analysis of target miRNAs of downregulated circRNAs and upregulated DEmiRNAs; (D) the dual coordinate plot of DECs and DEmiRNAs in HCC from the GSE128274 dataset. The dark orange dots and forest green dots represent DECs and DEmiRNAs with significance, respectively. The x-axis represents 5 potential miRNA targets with significance. HCC, hepatocellular carcinoma; DECs, differentially expressed circRNAs.

Figure 4

The detailed potential circRNA-miRNA interaction sites. (A,B) The detailed potential circRNA-miRNA interaction sites of targeted miRNAs with the highest context score percentile were obtained using TargetScan and miRanda data; (C) the expression level of miR-96 in TCGA_LIHC; (D) the expression level of miR-6746 in TCGA_LIHC. TCGA, The Cancer Genome Atlas; LIHC, liver hepatocellular carcinoma.

Identification of downstream target genes of miRNAs in HCC

It has been generally confirmed that miRNAs achieve their biological functions by inhibiting the expression of downstream molecules. To ascertain the downstream target mRNAs of the hsa_circ_0039466/miR-96-5p axis, the gene expression data from GSE128274 was analyzed. We obtained 1353 DEGs (available online: https://cdn.amegroups.cn/static/public/atm-22-3147-2.pdf) between HCC tissues and adjacent normal tissues, including 661 upregulated DEGs and 692 downregulated DEGs (). Subsequently, by intersecting the target mRNAs of hsa-miR-96-5p with the DEGs, 44 genes were screened out for further analysis (), and displays the 10 mRNAs with the greatest difference in expression.

Figure 5

Identification of potential target mRNAs of miR-96-5p in HCC. (A) Volcano plot of DEGs in HCC from the GSE128274 dataset. The red dots and blue dots represent upregulated DEGs and downregulated DEGs with significance (adjust P<0.05 and |log2FC| >1), respectively. The black dots represent DEGs without significance. (B) The intersection analysis of target genes of downregulated miR-96-5p and DEGs. (C) The network consisting of target miRNAs (miR-96-5p) and 44 mRNAs was generated by Cytoscape 3.7.1 software. HCC, hepatocellular carcinoma; DEGs, differentially expressed genes.

Table 3

The top 10 DEGs between HCC tissues and adjacent normal tissues

miRNA	Gene	Log2FC	AveExpr	t	P	Adj.P
miR-96-5p	LEPR	−3.431	3.044	−6.138	<0.001	0.006
miR-96-5p	EPHA3	−2.678	3.431	−4.151	0.002	0.022
miR-96-5p	CSRNP1	−2.667	5.999	−6.188	<0.001	0.006
miR-96-5p	SLC25A25	−2.457	4.475	−4.383	0.001	0.019
miR-96-5p	RGPD1	−2.317	4.895	−4.528	0.001	0.017
miR-96-5p	VSIG10	2.204	5.415	6.57	<0.001	0.005
miR-96-5p	TRIB3	2.348	6.185	3.706	0.005	0.032
miR-96-5p	SLC50A1	2.555	2.587	4.683	0.001	0.015
miR-96-5p	FN1	2.743	5.483	5.755	<0.001	0.008
miR-96-5p	ANLN	4.403	3.076	7.601	<0.001	0.003

DEGs, differentially expressed mRNAs; HCC, hepatocellular carcinoma.

Establishment of a possible circRNA-miRNA-mRNA regulatory network in HCC

In order to demonstrate the biological mechanism of the hsa_circ_0039466/miR-96-5p axis in HCC, we conducted pathway enrichment analysis for these 44 target genes. The results of GO analysis showed that positive regulation of cell projection organization, cortical actin cytoskeleton, and carbohydrate derivative transmembrane transporter activity were activated (). The results of KEGG analysis indicated that the insulin signaling pathway, insulin resistance, human papillomavirus infection, and the AMPK signaling pathway were enriched, which suggests that the hsa_circ_0039466/miR-96-5p axis is closely related to insulin metabolism (). In addition, the insulin signaling pathway and insulin resistance had the highest weight in the enrichment map (), which further suggests that insulin metabolism is closely related to the hsa_circ_0039466/miR-96-5p axis in HCC.

Figure 6

Biological pathways regulated by 44 target mRNAs of miR-96-5p in HCC. (A) The GO enrichment analysis plot displayed the BPs, MFs, and CCs regulated by 44 target mRNAs of miR-96-5p in HCC. (B) The KEGG plot exhibited enriched signaling and cellular pathways regulated by 44 target mRNAs of miR-96-5p in HCC. (C) Enrichment map of the top 10 enriched terms. The enrichment map organized the 10 enriched terms into a network with edges connecting overlapping gene sets. (D) Network plot of the top 10 enriched terms. The dots represent genes which are involved in these significant terms. (E) Correlation matrix plot showed the correlation of 8 overlapping genes and miR-96-5p. HCC, hepatocellular carcinoma; GO, Gene Ontology; BPs, biological processes; MFs, molecular functions; CCs, cellular components; KEGG, Kyoto Encyclopedia of Genes and Genomes.

We conducted correlation analysis on the genes related to the top 10 enrichment pathways to construct the circRNA-miRNA-mRNA regulatory subnetwork in HCC (). Only 8 genes were enriched in the top 10 enrichment pathways, among which 8 pathways were enriched in mTOR and 7 pathways were enriched in FOXO1 (). Subsequently, we conducted correlation analysis between miR-96-5p and the gene expression values of these 8 targets (TCGA dataset, ). The results showed that FOXO1 and LEPR negatively correlated with miR-96-5p, a result which was significant. According to the regulation principles of ceRNA, the hsa_circ_0039466/miR-96-5p/FOXO1 axis and hsa_circ_0039466/miR-96-5p/LEPR axis might play an important role in HCC.

Finally, we verified the feasibility of the hsa_circ_0039466/miR-96-5p/FOXO1 axis and hsa_circ_0039466/miR-96-5p/LEPR axis in HCC. The results demonstrated that FOXO1 has a low expression in HCC tissues in the TCGA database (P<0.001). Additionally, FOXO1 had weaker staining in HCC tissues than in adjacent normal tissues, with fewer positive areas in IHC, and the OS rates of patients with high FOXO1 expression levels were higher than those of low expression levels (P=0.016) (). We also found that LEPR, which is another mRNA target of miR-96-5p, has a low expression in HCC tissues in the TCGA database (P<0.001) (), however, there was no statistical difference in LEPR protein level between HCC tissues and adjacent normal tissues as determined by IHC (). Furthermore, no statistical difference was found between OS rates for low and high LEPR (P=0.83) ().

Figure 7

Validation of target mRNA (FOXO1 and LEPR) gene expression and survival analysis in HCC using TCGA data. (A) Validation of the expression of FOXO1 using TCGA data, where P<0.05 represented a significant difference. (B) Immunohistochemistry images of FOXO1 in normal (left) and tumor (right) tissues. The protein expression level of FOXO1 was significantly higher in HCC tissues than in normal tissues. (C) Kaplan-Meier survival curves displayed OS for HCC patients with high and low levels of FOXO1 mRNA expression. (D) Validation of the expression of LEPR using TCGA data, where P<0.05 represented a significant difference. (E) Immunohistochemistry images of LEPR in normal (left) and tumor (right) tissues. (F) Kaplan-Meier survival curves demonstrated OS for HCC patients with high and low levels of LEPR mRNA expression. HCC, hepatocellular carcinoma; TCGA, The Cancer Genome Atlas; OS, overall survival.

In summary, the hsa_circ_0039466/miR-96-5p/FOXO1 axis might play a key role in HCC. Hsa_circ_0039466 might promote the expression of FOXO1 by suppressing the inhibitory effect of miR-96-5p on FOXO1, ultimately halting tumor progression.

Validation of the existence of the hsa_circ_0039466/miR-96-5p/FOXO1 ceRNA regulatory network in HCC tissues

To further verify the hsa_circ_0039466/miR-96-5p/FOXO1 ceRNA regulatory network in HCC, we detected the expression levels of hsa_circ_0039466, miR-96-5p, and FOXO1 in 124 pairs of HCC tissues and adjacent normal tissues. Correlation analysis among hsa_circ_0039466, miR-96-5p, and FOXO1 was performed, and the relationship between hsa_circ_0039466/miR-96-5p/FOXO1 and the clinical stage of HCC patients was also analyzed. The qRT-PCR results showed that hsa_circ_0039466 and FOXO1 had a low expression in HCC tissues, while miR-96-5p was highly expressed in HCC tissues compared with adjacent normal tissues (, n=124). In addition, the results of Spearman correlation analysis showed that the expression level of hsa_circ_0039466 correlated negatively with the expression level of miR-96-5p (Spearman correlation =−0.574; P<0.001), but correlated positively with the expression level of FOXO1 (Spearman correlation =0.647; P<0.001), while the miR-96-5p expression correlated negatively with FOXO1 expression (Spearman correlation =−0.576; P<0.001) (). These results further demonstrated the possibility of the existence of the hsa_circ_0039466/miR-96-5p/FOXO1 ceRNA regulatory axis. Next, we further analyzed the relationship between the expression level of hsa_circ_0039466/miR-96-5p/FOXO1 and clinical data such as gender, age, cancer stage, Barcelona clinic liver cancer (BCLC) stage and T stage. The results showed that the expression levels of hsa_circ_0039466, miR-96-5p, and FOXO1 were not related to age and gender (). On the other hand, hsa_circ_0039466, miR-96-5p, and FOXO1 expression were significantly associated with cancer stage, BCLC and T stage (). In particular, the expression level of hsa_circ_0039466 or FOXO1 in advanced HCC stage III + IV, BCLC stage C, and stage T (T3 + T4) was significantly lower than that in stage I + II, BCLC stage A + B, and stage T (T1 + T2), while miR-96-5p showed the opposite trend (). These clinical data further support the existence of the hsa_circ_0039466/miR-96-5p/FOXO1 ceRNA regulatory network in HCC and that it may play an important role in prognosis of HCC patients.

Figure 8

Validation of the existence of the hsa_circ_0039466/miR-96-5p/FOXO1 ceRNA regulatory network in HCC tissues. (A) The relative expression levels of hsa_circ_0039466, miR-96-5p, and FOXO1 in HCC tissues and adjacent normal tissues measured by qRT-PCR (n=124). (B) Spearman correlation analysis among hsa_circ_0039466, miR-96-5p, and FOXO1 expression detected by qRT-PCR (n=124). (C) The correlation between the expression level of hsa_circ_0039466, miR-96-5p, or FOXO1 and the gender and age of HCC patients. (D) The correlation between the expression level of hsa_circ_0039466, miR-96-5p, or FOXO1 and the cancer stage, BCLC stage, and T stage of HCC patients. Data were presented as mean ± SD. ns, no significance; **P<0.01; ***P<0.001. ceRNA, competing endogenous RNA; HCC, hepatocellular carcinoma; qRT-PCR, quantitative real-time polymerase chain reaction; BCLC, Barcelona Clinic Liver Cancer.

Discussion

The current study presented a potential circRNA-miRNA-mRNA ceRNA regulatory network in the pathogenesis of HCC through bioinformatics analysis. In general, hsa_circ_0039466 promotes the expression of FOXO1 by sponging miR-96-5p, and ultimately inhibits tumor progression.

Growing evidence suggests that circRNAs play a critical role in the initiation and development of human cancers including HCC (26). Although the circRNA-miRNA-mRNA regulatory network has been reported by some studies in HCC (27,28), the mechanism and pathways of the circRNA-related ceRNA network in HCC is yet not fully understood and requires further research.

After using differential expression analysis and retrieval of matched CSCD data, 8 cancer-specific circRNAs were selected for further research. A series of studies have verified that circRNAs, acting as miRNA sponges, promote the expression of downstream genes (29-31). Therefore, in our study, potential miRNA targets of the 8 circRNAs were predicted. We found that the hsa_circ_0039466/miR-96-5p axis may play a key role in HCC through step-by-step analysis, including miRNA expression analysis and circRNA-miRNA correlation analysis.

miRNAs play important roles in genetic regulatory networks and are closely related to many biological processes including cell proliferation, apoptosis, tumor genesis, development, and immune response (32,33). The principle of the biological function of miRNAs can be summarized as degrading the target mRNA by binding to the 3'-UTR, so any ceRNA that can “sponge” miRNA binding may “protect” the mRNA in this way (16). As a kind of ceRNA, circRNAs, due to their high abundance and stability, are more effective than other ceRNAs in miRNA binding (17). In this biological process, the stoichiometric relationship between the miRNA binding site of the circRNA and the mRNA target site of the miRNA must be considered (19,20). There is no relevant data in the literature discussing hsa_circ_0039466, so it is urgent to explore its potential role in HCC.

miR-96-5p may play a role as a tumor suppressor molecule in pancreatic cancer and nasopharyngeal carcinoma (34,35). miR-96-5p may also play a carcinogenic role in breast cancer, gastic cancer, and head and neck squamous cell carcinoma (HNSCC) (36,37). Downregulated miR-96-5p has been shown to induce cell apoptosis in gastric adenocarcinoma (38). miR-96-5p inhibits hepatic stellate cell (HSC) activation by inhibiting autophagy, thereby participating in the regulation of liver fibrosis (39). miR-96-5p may play an important role in HCC by decreasing caspase-9 expression (40). However, the exact role and mechanism of miR-96-5p in HCC is still unclear and needs further investigation.

To further explore the role of the hsa_circ_0039466/miR-96-5p axis in HCC, we analyzed the specific downstream molecular mechanisms of this regulatory axis. Through cross-analysis of the target gene set of miR-96-5p and the DEG dataset, we obtained 44 possible targets. Subsequently, we conducted pathway enrichment for these 44 target genes. The hsa_circ_0039466/miR-96-5p axis affects the progression of HCC by regulating insulin metabolism. The KEGG enrichment pathways found in this study, namely the insulin signaling pathway, insulin resistance, and AMPK signaling pathway, were found to be strongly associated with cancer progression, including that of HCC (41-44).

In order to find the core target genes of the hsa_circ_0039466/miR-96-5p axis in HCC, we conducted correlation analysis on the genes related to the top 10 enrichment pathways. FOXO1 and LEPR were confirmed to significantly correlate negatively with miR-96-5p. Combining these data, a circRNA-miRNA-mRNA (hsa_circ_0039466/miR-96-5p/FOXO1 and hsa_circ_0039466/miR-96-5p/LEPR) triple ceRNA subnetwork regulated by hsa_circ_0039466 was successfully established in HCC. Besides, the expression correlation analysis and the expression levels of RNAs in the network were further verified using TCGA data, confirming the analytical accuracy of bioinformatics analysis. Finally, we found that the hsa_circ_0039466/miR-96-5p/FOXO1 axis may be the most important one related to progression of HCC. FOXO1 is a transcriptional regulator of the G1/S checkpoint and apoptosis and acts as a tumor suppressor in many tumors, such as cervical cancer, prostate cancer, and gastric cancer (45-48). A large number of studies have reported that FOXO1 plays an important role in the occurrence and development of HCC and sorafenib resistance (45,49-57). However, the related role of FOXO1 in HCC is still not fully understood and needs further investigation.

We successfully established a hsa_circ_0039466/miR-96-5p/FOXO1 axis triple ceRNA subnetwork regulated by hsa_circ_0039466 in HCC. As shown in , FOXO1 is the key gene of the insulin signaling pathway, insulin resistance, the AMPK signaling pathway, among other important pathways. FOXO1 is under-expressed in HCC tissues () and low expression of FOXO1 in tumor tissues suggests worse prognosis in HCC patients (). Furthermore, we detected the expression of hsa_circ_0039466, miR-96-5p, and FOXO1 in HCC tissue samples. We found that hsa_circ_0039466 and FOXO1 have a low expression in HCC tissues, while miR-96-5p was highly expressed as compared with adjacent normal tissues. Moreover, the expression level of miR-96-5p correlated negatively with the expression levels of hsa_circ_0039466 or FOXO1, while there was a positive correlation between hsa_circ_0039466 and FOXO1. HCC patients in earlier clinical stages had higher expression levels of hsa_circ_0039466 or FOXO1 than those in advanced stages, whereas miR-96-5p showed the opposite. We have every reason to believe that hsa_circ_0039466 promotes the expression of FOXO1 by suppressing the inhibitory effect of miR-96-5p on FOXO1, therefore, inhibiting tumor progression. We analyzed the results of each component, which mutually verified the above conclusions to be reliable.

Our study provides a new idea for the in-depth understanding of the specific molecular mechanism of the circRNA-miRNA-mRNA network in hepatocarcinogenesis. The research approach may also be used to probe the ceRNA molecular interaction mechanism of this network in the occurrence and development of other human cancer types. Certainly, more in vitro and in vivo experiments are needed to further elucidate the role of the ceRNA regulatory network established by our study in HCC. Meanwhile, extensive future multicenter HCC clinical data is needed to evaluate the diagnostic and prognostic value of each component of this established regulatory network, which will help uncover promising biomarkers for the diagnosis and prognosis of HCC patients.

Conclusions

In summary, our study revealed hsa_circ_0039466, miR-96-5p, and FOXO1 ceRNA have a potential regulatory network in HCC by whole transcriptome analysis. In general, hsa_circ_0039466 might promote the expression of FOXO1 by suppressing the inhibitory effect of miR-96-5p on FOXO1, ultimately inhibiting HCC tumor progression. Targeting components in this ceRNA regulatory network may be a promising approach to assist in providing diagnostic and prognostic parameters as well as uncover potential targets for treatment of HCC.

The article’s supplementary files as