The TCGA-LIHC database (https://portal.gdc.cancer.gov/) was used to download the associated clinical data along with the mRNA expression profiling data (including 374 HCC and 50 paracancer tissues samples) and miRNA expression profiling data (including 375 HCC and 50 paracancer tissue samples). The gene expression profiling data of the LIRI-JP (Liver cancer-RIKEN, Japan) cohort, comprising 232 tumor samples, was acquired from the ICGC database (https://dcc.icgc.org/) together with its corresponding clinical information. Demographic and clinicopathological characteristics of HCC patients in the TCGA-LIHC and LIRI-JP cohorts are summarized in Supplementary Table S1. Additionally, gene expression profiles for the GSE54236 and GSE101685 datasets were retrieved from the GEO database.

Firstly, genes with expression level of 0 in more than half of the samples were excluded. The differential expression analysis of mRNA and miRNA between HCC and paracancer tissues was then calculated using DESeq2 R package. The following criteria were used to identify differential genes: adjusted p < 0.05 and |log2 (fold change)| >1.

The downstream target miRNAs of circELMOD3 were initially predicted using the circMine database (http://www.biomedical-web.com/circmine/home) with default parameters. Subsequently, the predicted miRNAs were intersected with upregulated miRNAs identified in the TCGA-LIHC dataset. The downstream target mRNAs of these miRNAs were predicted using TargetScan 8.0 (https://www.targetscan.org/vert_80/) with default parameters and miRWalk 2.0 databases (http://mirwalk.umm.uni-heidelberg.de/), with binding sites restricted to the 3′UTR region. Only mRNAs predicted in both databases were retained. These selected mRNAs overlapped with downregulated mRNAs in TCGA-LIHC. Finally, the circELMOD3-miRNA and miRNA-mRNA interactions were used to construct the ceRNA network via Cytoscape v3.9.1 software.

To evaluate the potential biological functions and pathways of mRNAs in the ceRNA network, we performed Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichment studies using the “clusterProfiler” package. p-value less than 0.05 was defined as cutoff criterion.

We initially collected clinical data, excluded samples with missing survival time or event status, and removed patients with an overall survival time of less than 1 month prior to performing survival analysis. The Kaplan-Meier curve was employed to perform survival analysis on miRNA and mRNA through the utilization of the “survival” and “survminer” packages. p-value less than 0.05 was defined as cutoff criterion.

LASSO is a penalized regression technique that can be employed to screen variables from high-dimensional data for constructing prognostic models. Firstly, prognostic genes in the ceRNA network were identified using univariate Cox regression, with p < 0.05 being defined as cutoff criterion. Subsequently, the relevant genes that influence prognosis were identified via the glmnet LASSO regression method. The following formula is utilized to calculate the risk score: risk score = (EXP_gene1 × coefficient_gene1) +(EXP_gene2 × coefficient_gene2) + (EXP_genex × coefficient_genex). Patients were categorized into high-risk and low-risk groups according to their median risk score, and the survival difference between the two groups was then assessed using a log-rank test. Furthermore, model performance was evaluated by employing the ‘timeROC’ package to draw the receiver operating characteristic (ROC) curve. We validated the prognostic risk model from the LIRI-JP cohort by calculating individual patient risk scores according to the aforementioned formula. Furthermore, we evaluated predictive value using time-dependent ROC curve analysis and Kaplan-Meier analysis.

The nomogram is widely used for prognostic prediction in tumor patients, univariate and multivariate Cox regression were performed to identify independent prognostic predictors in HCC patients, including gender, stage, age and risk score. A nomogram incorporating these prognostic predictors was constructed using the “rms” package. The fitting and prediction power of our prognostic model were evaluated by plotting calibration curves.

The TCGA-LIHC dataset was subjected to GSEA. Based on the median gene expression levels, HCC patients were stratified into high- and low-expression groups to conduct GSEA. Pathways exhibiting significant enrichment in these two groups were subsequently analyzed. The criteria for statistical significance were FDR < 0.25 and p < 0.05.

HCC and adjacent tissues were collected as previously described (Lai, et al., 2023). The informed consent was obtained from each patient and the study was approved by the Ethics Committee of Guangxi Medical University (20220146). Human liver cancer cell line MHCC97H was purchased in CellCook Biotech Co., Ltd. and Hep3B cell line from Cell bank of the Chinese academy of sciences. The cell culture and cell transfection procedure were performed following earlier depicted protocols (Lai, et al., 2023). Trizol reagent was used to extract total RNA, and the absorbance was measured by a microplate reader. Using PrimeScript RT Master Mix (Takara, RR036A), total RNA was reverse-transcribed into cDNA, followed by qRT-PCR analysis on a Stepone Plus system (Applied Biosystem, United States) with a qPCR detection kit (GenStar, A303). GAPDH served as the internal reference for mRNA quantification. We used the 2^-△△CT method to perform a relative gene expression study. Supplementary Table S2 lists the primer sequences that were used in this research.

The Western blotting analysis was performed following earlier depicted protocols (Lai, et al., 2023). The primary antibodies used were N-cadherin (1:1000, 22018-1-AP), E-cadherin (1:2500, 20874-1-AP), CDK4 (1:2000, 11026-1-AP), CDK6 (1:2500, 14052-1-AP), CyclinD1 (1:10000, 60186-1-lg), GAPDH (1:10000, 60004-1-lg), Tubulin (1:6000, 11224-1-AP) and all antibodies are sourced Proteintech. Subsequently, the membrane was incubated with the secondary antibody for 1 hour at room temperature on a shaker. Protein bands were visualized by adding an enhanced chemiluminescence (ECL) substrate and imaging with an iBright FL1000 system.

R software (version 4.3.2) and GraphPad Prism software (version 8.3) were used to perform data analyses. The relationship between two gene variables was investigated using Spearman correlation analysis, and the difference between the two groups was assessed using the Wilcoxon rank-sum test. Univariate Cox analysis was applied to identify potential risk factors affecting prognosis. Using multivariate Cox analysis, independent prognostic factors were identified. The predictive accuracy of the HCC prognostic model was validated through receiver operating characteristic (ROC) curve analysis. All in vitro experiments included three independent biological replicates, and data are presented as mean ± standard deviation. p-value less than 0.05 was defined as cutoff criterion.

The detailed flow chart of study process is shown in Figure 1. Differentially expressed mRNAs and miRNAs in HCC were performed using the TCGA-LIHC dataset. The results showed a total of 3982 significantly dysregulated mRNAs, with 2853 showing upregulated and 1130 displaying downregulated (Supplementary Figures S1A, B). Furthermore, 225 differentially expressed miRNAs were found, of which 196 were upregulated and 29 were downregulated (Supplementary Figures S1C, D).

To identify potential miRNAs that might bind to circELMOD3, the circMine database was utilized. The ceRNA hypothesis suggests that the expression of miRNA and circRNA should be negatively related. By intersecting upregulated miRNAs from the TCGA-LIHC dataset with circMine-predicted miRNAs, we identified 17 candidate miRNAs (Figure 2A). Considering circELMOD3 was previously found to affect various biological functions of HCC cells, we further explored the relationship between circELMOD3-binding miRNAs and the survival of HCC patients. Log-rank tests and univariate Cox regression analysis were carried out on these miRNAs (Supplementary Tables S3, 4). Prognostic miRNAs were selected based on significance thresholds of p < 0.05 in both analyses, yielding five miRNAs (including hsa-miR-106b-5p, hsa-miR-301a-5p, hsa-miR-760, hsa-miR-3127-5p, hsa-miR-3677-3p) with potential circELMOD3 interactions. The relative expression levels of five prognosis-related miRNAs in HCC were significantly upregulated compared to paracancer tissues in the TCGA-LIHC database (Figures 2B–F). Kaplan-Meier analysis displayed that high expression levels of five identified miRNAs were associated with poor overall survival in HCC patients (Figures 2G–K), The findings suggest that these five miRNAs may have oncogenic properties in HCC.

To comprehensively comprehend the role of circELMOD3 through ceRNA network in HCC, we constructed a circELMOD3-centered miRNA-mRNA interaction network. Candidate mRNAs potentially binding to the five miRNAs were predicted using miRWalk and TargetScan databases. Subsequently these predicted mRNAs were further intersected with downregulated mRNA identified in the TCGA-LIHC database, resulting in a set of candidate mRNAs (Supplementary Figures S2A–E). Finally, by integrating the predictive relationships between circELMOD3-miRNA and miRNA-mRNA interactions, we constructed and visualized a ceRNA regulatory network comprising circELMOD3, 5 miRNAs, and 274 mRNAs using Cytoscape software (Figure 3A).

To explore the regulatory role of circELMOD3 in ceRNA network, we carried out functional enrichment studies using KEGG and GO. The results of GO analysis exhibited that biological functions (BP) primarily encompassed organic acid catabolic process, carboxylic acid catabolic process, and epithelial cell proliferation. In terms of cellular component (CC), it was mainly associated with membrane raft, membrane microdomain, and phosphatidylinositol 3-kinase complex. Regarding molecular function (MF), the circELMOD3 regulatory gene was predominantly linked to metal ion transmembrane transporter activity, phosphatidylinositol 3-kinase regulator activity, and transaminase activity (Figure 3B). Furthermore, the KEGG pathway enrichment analysis demonstrated that circELMOD3 regulatory genes are mainly involved in cell adhesion molecules, JAK-STAT signaling pathway, as well as cell adhesion molecules (Figure 3C). These findings suggest a pivotal role for circELMOD3 in HCC development.

From the TCGA-LIHC cohort, 342 HCC patients with complete survival data and a minimum survival time of 1 month were included based on predefined screening criteria. We employed univariate Cox regression analysis with p < 0.05 and HR < 1 as the thresholds to screen for prognostic related genes in the circELMOD3-regulated ceRNA network, yielding 46 candidate genes (Supplementary Tables S5, 6). Subsequently, we further employed LASSO regression analysis to screen genes, selected the optimal penalty coefficient λ (λ = 0.072) by 10-fold cross-validation, and obtained 4 genes with non-zero coefficients (ALDH2, DNASE1L3, STARD5, ACOT12), and the following formula was used to determine the risk score: risk score = (EXP _ALDH2 *-0.02341441) + (EXP _DNASE1L3 *-0.09186373) + (EXP _STARD5 *-0.15196182) + (EXP _ACOT12 *-0.02724502) (Figures 4A, B). HCC Patients were categorized as low-risk or high-risk depending on their median risk scores. Patients in the low-risk group exhibit higher survival probabilities compared to those in the high-risk group (Figure 4C). Furthermore, expression levels of all four prognostic genes were considerably lower in the high-risk group than in the low-risk group; patient mortality increased significantly with increasing risk values (Figure 4D). The time-dependent receiver operating characteristic (time-ROC) curve was used to assess the predictive power of this model; the area under the 1-year, 3-year, and 5-years ROC curves, respectively, were 0.734, 0.718, and 0.707 (Figure 4E).

To validate the robustness of the model, patients in the LIRI-JP cohort were stratified into low- and high-risk groups based on the same risk score calculation method used for the TCGA-LIHC cohort. The results showed that compared to those with higher risk scores patients, with lower risk scores patients exhibited higher survival rates and decreased mortality (Figure 4F). Furthermore, prognostic genes were considerably downregulated in the high-risk group compared to the low-risk group, patient mortality increased significantly with increasing risk values (Figure 4G). Additionally, time-ROC curve analysis displayed strong predictive ability for the LIRI-JP cohort, with AUC values of 0.742, 0.740, and 0.730 for predicting one-, two-, and three-year survival respectively (Figure 4H).

To assess prognostic stability of the model across diverse clinical subgroups, we performed stratified analyses based on age (≥65 years vs. < 65 years), sex (female vs. male), and stage (III-IV vs. I-II). The findings consistently showed a significantly lower risk of death in the low-risk group compared with the high-risk group. Although the survival difference in the female subgroup was not statistically significant, the low-risk group still showed a lower mortality risk when follow-up duration was limited to less than 70 months. (Supplementary Figures S3A–C). Similarly consistent outcomes were observed when performing stratified analyses within LIRI-JP cohort (Supplementary Figures S3D–F).

We developed a nomogram integrating the risk score (calculated as the sum of the expression levels of four LASSO-selected genes multiplied by their respective coefficients) and stage to enhance the clinical utility of the prognostic model in the TCGA-LIHC cohort. Firstly, the study employed univariate Cox regression analysis to investigate the potential prognostic predictive value of the risk score and specific clinical factors. The results demonstrated significant association between both the risk score (HR = 9.489, 95%CI: 4.098–21.970) and stage (HR = 1.784, 95%CI: 1.445–2.202) with overall survival in HCC patients (Figure 5A). Furthermore, multivariate cox regression analysis demonstrated that both the risk score (HR = 5.898, 95%CI: 2.344–14.837) and stage (HR = 1.608, 95%CI: 1.289–2.006) independently influenced OS among HCC patients (Figure 5B). Finally, we constructed nomograms incorporating risk scores and stage to accurately predict one-, three-, and five-year survival for HCC patients (Figure 5C). Higher total points on the nomogram correlated with poorer survival outcomes; moreover, the calibration plots for one-, three-, and five-year survival probabilities showed great concordance with actual observations (Figure 5D). Moreover, time-ROC analysis shown that, in comparison to the risk score and stage, the combined nomogram performed best in predicting the overall survival of HCC patients (AUC = 0.761 at 1 year, 0.768 at 3 years, and 0.725 at 5 years) (Figure 5E).

To further characterize the four prognostic genes, we observed their downregulation in HCC using data from the TCGA-LIHC database (Figure 6A), as well as in paired samples (Figure 6B). These findings were validated in two independent GEO datasets (Supplementary Figures S4A, B). Correlation analysis revealed a positive association among the four prognostic genes (Figure 6C). Additionally, Kaplan-Meier survival analysis demonstrated that high expression of ALDH2, DNASE1L3, STARD5, and ACOT12 was associated with prolonged survival in HCC patients (Figures 6D–G). GSEA results showed that all four prognostic genes were negatively correlated with the cell cycle pathway, indicating that relatively high expression levels of the four prognostic genes could inhibit the cell cycle pathway compared with low expression levels. (Figures 6H–K).

Based on the GSEA results indicating the association of the four prognostic genes with the cell cycle pathway, and KEGG enrichment analysis showed that they were related to cell adhesion molecule pathway. We hypothesized that circELMOD3 might be related to cell cycle and EMT pathways. Subsequently, Western blot was performed to verify the expression levels of cell cycle pathway-related proteins CDK4, CDK6, CyclinD1 and EMT-related markers N-cadherin and E-cadherin. The results exhibited that silencing circELMOD3 in Hep3B cells led to an increase in the relative expression levels of CDK4, CDK6, CyclinD1 and N-cadherin proteins while decreasing the relative expression levels of E-cadherin proteins (Figures 7A, C). Conversely, overexpression of circELMOD3 in MHCC97H cells resulted in a decreased in the relative expression levels of CDK4, CDK6, CyclinD1 and N-Cadherin proteins with an increase observed for E-cadherin protein (Figures 7B, D). To further validate the expression levels of the four prognostic related genes in HCC tissues, we performed qPCR validation using 23 pairs of HCC tissues. The results showed that ALDH2, ACOT12, STARD5 and DNASE1L3 were significantly downregulated in HCC tissues compared with paracancer tissues. (Figures 7E–H).