Data acquisition involved three steps. First, HCC gene expression profiles and matched clinical data were acquired from The Cancer Genome Atlas (TCGA; https://www.cancer.gov), including 424 tissue samples (374 tumor, 50 normal) and corresponding clinical data for 377 patients. Second, survival-specific expression data for 118 HCC cases were independently obtained from the Gene Expression Omnibus (GEO) dataset (GSE10186; https://www.ncbi.nlm.nih.gov). Third, a list of 246 SRGs was generated by screening the GeneCards database (https://www.genecards.org), retaining only genes with a relevance score $>$15 [37]. The matrix data pertaining to TCGA, GEO, and SRGs are comprehensively detailed in the Appendix.

Bioinformatics analysis was conducted using Perl (v5.30.0) and R (v4.3.2). Gene expression data were formatted into a matrix, featuring genes as rows and samples as columns. For genes with multiple probes, expression values were averaged to produce a single value per gene. Genes with zero mean expression across all samples were removed. Finally, the analysis was limited to the intersection of this filtered expression matrix and the curated SRG set. Differential expression analysis between normal and tumor tissues was conducted using the Wilcoxon rank-sum test. Differentially expressed genes (DEGs) were defined as those with a false discovery rate (FDR) 0.05 and logFC $>$1. These DEGs were subsequently subjected to downstream functional and prognostic evaluations. Univariate Cox proportional hazards regression analysis was used to evaluate the impact of gene expression on overall survival (OS). Genes with a Cox p-value 0.05 were considered significantly associated with prognosis. Hazard ratios (HRs) and 95% confidence intervals (CIs) were calculated for each. Results were visualized using forest plots, with HR $>$1 indicating risk factors and HR 1 indicating protective factors.

Tumor transcriptomic data from TCGA (https://www.cancer.gov) and GEO (https://www.ncbi.nlm.nih.gov), encompassing survival-associated gene expression profiles, were subjected to molecular subtyping using an unsupervised consensus clustering approach. This was implemented via the ConsensusClusterPlus R package with the following parameters: clusterAlg = “km”, 80% pItem = 0.8, pFeature = 1, and reps = 500 to ensure robustness. The consensus matrix and cumulative distribution function (CDF) were evaluated across a predefined range of cluster numbers (k = 2 to 9) to determine the optimal k value. The optimal number of stable clusters (k = 2) was selected based on the criteria of maximal intercluster divergence and intracluster concordance, as reflected by the consensus matrix heatmap and relative change in area under the CDF curve. The resulting subtypes were designated Cluster A and Cluster B. Subsequently, to evaluate the prognostic value of the identified subtypes, we integrated OS data. Survival differences between clusters were visualized using Kaplan–Meier analysis and statistically compared via the log-rank (Mantel-Cox) test. Furthermore, HRs and their CIs were also determined using a univariate Cox proportional hazards model, with all analyses conducted using the survival and survminer R packages.

To develop a prognostic risk signature based on SRGs, a combination of least absolute shrinkage and selection operator (LASSO) regression and multivariate Cox proportional hazards regression was employed. The data were randomly divided into training (70%) and testing (30%) sets using the caret package in R, which was used for all statistical analyses. To identify genes linked to survival, a LASSO-penalized Cox proportional hazards model was applied to the training cohort using the glmnet package. The optimal regularization parameter ($\lambda$) was selected based on the minimum criteria within a 10-fold cross-validation framework. Genes with non-zero coefficients at optimal $\lambda$ (lambda.min) were selected to refine the list. Multivariate Cox proportional hazards regression analysis was performed on these retained genes to construct the definitive prognostic model. To refine the model, bidirectional stepwise regression (in both directions) was performed, and the resulting final multivariate Cox model was used to calculate a prognostic risk score for each patient (Mathematical Components 1).

(1)$$Riskscore={\sum }_{i=1}^{n}Coe{f}_i\times x_i$$

$Coe{f}_i$ represents the correlation coefficient of prognostic genes; $x_i$ represents the expression of prognostic genes. Ultimately, multivariate Cox regression analysis was conducted to determine significant clinical prognostic genes.

(2)$$Surv(futime,fustat)=Age+Sex+Stage+RiskScore$$

(3)$$GAPD{H}_{\text{activity}}\left(U/\text{mg prot}\right)=\frac{\mathrm{\Delta }A}{\epsilon \times d}\times {10}^9\times \frac{V_{\text{total}}}{V_{\text{sample}}\times C_{\text{pr}}}\times \frac{1}{T}$$

The absorbance at 340 nm was recorded immediately ($A_1$), and again after a 5-min incubation at 37 °C ($A_2$). $\mathrm{\Delta }A=A_1-A_2$; $\epsilon$ = molar extinction coefficient of NADH at 340 nm (6.22  $\times$  10³ L/mol/cm); $d$ = light-path length (0.6 cm for a 96-well UV plate); $V_{\text{total}}$ = total reaction volume (0.2 mL); $V_{\text{sample}}$ = volume of sample added to the reaction (0.006 mL); $C_{\text{pr}}$ = protein concentration of the sample (mg/mL); $T$ = reaction time (5 min); 10⁹ = 1 mol = 10⁹ nmol.

(4)$$F/F_0=\frac{F_t}{F_0}$$

Here, F₀ represents the baseline fluorescence intensity prior to stimulation, and F_t represents the fluorescence intensity at time t. The resulting F/F₀ versus time curves were used to characterize dynamic intracellular Ca²⁺ responses.

Using the median training set risk score as a cutoff, patients in both the training and testing cohorts were stratified into high- and low-risk groups. Survival differences between these groups were assessed using the log-rank test. Time-dependent receiver operating characteristic (ROC) analysis (timeROC package) was conducted to assess the 3-year prognostic accuracy of the model. The model was validated across both cohorts based on specific performance metrics: training set (log-rank p 0.01, AUC $>$0.68) and testing set (p 0.05, AUC $>$0.65).

The prognostic capacity of the risk score was assessed separately in the training, testing, and combined cohorts. For each cohort, OS probabilities between the predefined high- and low-risk groups were estimated via the Kaplan–Meier method. To evaluate the statistical significance of survival differences, the log-rank (Mantel-Cox) test was employed, with survival curves visualized with the survminer package. To quantify the time-dependent predictive accuracy of the risk score, a time-dependent ROC analysis was performed using the timeROC package. This analysis was conducted at three clinically relevant time points: 1, 3, and 5 years post-diagnosis. To evaluate the model’s discriminative power, the area under the curve (AUC) and its 95% CIs were calculated for each time point. A custom function was used to plot all three ROC curves on a single graph.

To assess whether the SRG-derived risk score serves as an independent predictor of OS, multivariate Cox proportional hazards regression analysis was performed. Clinical covariates (including age, sex and tumor stage) were incorporated alongside the continuous risk score. Prior to analysis, samples with any missing clinical information (coded as “unknown”) were excluded. To ensure a consistent, matched cohort, only samples present in both the risk score and complete clinical datasets were retained. The Cox proportional hazards model, defined by Mathematical Components 2, was then fitted using a combined dataframe that included survival time (futime), status (fustat), all clinical variables, and the riskScore (Mathematical Components 2). HRs, 95% CIs, and Wald test p were extracted for each variable from the fitted model, and then summarized and exported. Finally, a forest plot was generated using the ggforest function from the survminer package to visually present the HRs and 95% CIs of all covariates, with the risk score highlighted.

To visually depict the expression profiles of the signature genes across the patient cohort, a heatmap was constructed using the pheatmap package. Patients were sorted by ascending risk order, after which the model gene expression matrix was extracted and transposed to align genes as rows and ordered patients as columns. To visualize relative expression patterns, row-wise z-score normalization (scale = “row”) was applied and mapped to a continuous blue-white-red color scale.

A prognostic nomogram was developed using the regplot package to provide quantitative estimates of 1-, 3-, and 5-year OS. This tool was based on a final multivariate Cox proportional hazards model that incorporated all significant independent predictors, including clinical variables and the continuous risk score. To assess the accuracy of the nomogram’s predictions, calibration curves were generated for the same three time points using the rms::calibrate function with 1000 bootstrap iterations. This method compares the predicted survival probabilities against the observed Kaplan–Meier estimates to determine performance. The model’s goodness-of-fit was visually verified by the proximity of the calibration curves to the 45-degree line. Cumulative hazard curves (via survminer::ggsurvplot with fun = “cumhaz”) were used to compare risk over time between groups divided by their median nomogram score. Prior to all analyses, data integrity was ensured by excluding samples with any missing clinical information and retaining only those patients present in both the risk score and clinical datasets.

To explore the immune microenvironment associated with high-risk genes, immune cell enrichment analysis was performed leveraging the Tumor Immune Single Cell Hub 2 (TISCH2) database (https://tisch.compbio.cn/gallery/). Separately, to elucidate the functional interactions among these genes, a protein–protein interaction network was constructed. The network data were retrieved from the STRING database (version 12.0; https://cn.string-db.org/), applying a minimum required interaction confidence score threshold to filter for high-confidence interactions. The resulting network was subsequently imported into Cytoscape software (version 3.9.0) for advanced visualization and topological analysis, enabling clearer representation of molecular interaction hubs and modules [38].

The human HCC cell lines Huh-7 (STCC10102; Servicebio, Wuhan, China) and HepG2 (STCC10114P; Servicebio) were cultured in Dulbecco’s Modified Eagle Medium (DMEM, G4511-500ML; Servicebio) supplemented with 10% fetal bovine serum (FBS, FSP500; ExCell Bio, Shanghai, China) and 1% penicillin–streptomycin (G4003; Servicebio). Cells were maintained under standard conditions at 37 ℃ in a 5% CO₂ humidified atmosphere, with daily monitoring of confluence and morphology, ensuring optimal growth before subsequent experimental use. Both Huh-7 and HepG2 cell lines have been authenticated by STR profiling and confirmed to be mycoplasma-negative.

To investigate the functional impact of ANXA5, Huh-7 and HepG2 cells were transduced with a recombinant lentiviral vector targeting ANXA5 (single hairpin RNA [shRNA]-ANXA5; Genebio, Shanghai, China). A non-targeting, scrambled shRNA lentivirus (shRNA-negative control [NC]) was used as the NC.

The specific shRNA sequences are listed in Table 1. Briefly, Huh-7 and HepG2 cells were seeded in 6-well plates; at approximately 80% confluence, they were transduced with lentiviruses in a mixture of DMEM and polybrene (40804ES76; Yeasen, Shanghai, China). After 12 h, the transduction mixture was replaced with fresh complete culture medium. Following an additional 3–4 days of culture, stable cell pools were selected using puromycin (HY-K1057; MedChemExpress, Shanghai, China) over a 3- to 4-day period. Transduction efficiency was initially assessed by observing green fluorescent protein expression under a fluorescence microscope. The knockdown (KD) efficiency of ANXA5 was subsequently confirmed at the mRNA and protein levels using quantitative PCR (qPCR) and Western blotting, respectively. All procedures were performed according to the manufacturer’s protocol (Genebio) and established laboratory methods [39, 40].

Total protein was isolated from control and transfected cells (Huh-7 and HepG2) with RIPA lysis buffer (G2002-100ML; Servicebio). Following extraction, lysates were clarified by centrifugation at 12,000 rpm for 10 min at 4 ℃. The protein concentration in the supernatants was quantified using a BCA assay kit (G2026-200T; Servicebio). Equal protein concentrations were resolved by 10% sodium dodecyl sulfate polyacrylamide gel electrophoresis (SDS-PAGE) and electrotransferred onto PVDF membranes. After blocking in 5% non-fat milk for 30 min at room temperature, the membranes were probed overnight at 4 °C with primary antibodies against ANXA5 (1:1000, 11060-1-AP; Sanying, Wuhan, China), GAPDH (1:1000, GB15004; Servicebio), and $\beta$-actin (1:5000, GB15003; Servicebio). Post-incubation, membranes were washed three times with Tris-buffered saline with Tween 20 and then incubated with horseradish-peroxidase (HRP)-conjugated goat anti-rabbit secondary antibody (1:5000, GB23303; Servicebio). Signal detection was performed using an enhanced chemiluminescence substrate, and bands were imaged on X-ray film.

Co-immunoprecipitation (Co-IP) assays were performed to assess the interaction between ANXA5 and GAPDH under different ANXA5 expression conditions. Huh-7 and HepG2 cells from the wild-type (WT), NC, and ANXA5-KD groups were lysed using ice-cold IP lysis buffer supplemented with protease inhibitors. After centrifugation at 12,000 $\times$g for 15 min at 4 ℃, the clarified supernatants were collected, and protein concentrations were quantified by the BCA assay. Equal amounts of total protein from each group were incubated overnight at 4 ℃ with an anti-ANXA5 antibody (IP-grade) or control IgG under gentle rotation. Protein A/G agarose beads were then added and incubated for an additional 2–4 h to capture immune complexes. The beads were washed repeatedly with lysis buffer to reduce nonspecific binding. IP’d proteins were eluted by boiling in SDS sample buffer and subjected to SDS-PAGE, followed by electrotransfer to PVDF membranes. Western blotting was performed using an anti-GAPDH antibody (Western blot-grade) to detect co-precipitated proteins. A fraction of whole cell lysate (input) was included to verify protein expression levels across different groups.

Total RNA extraction was performed using TRIzol Reagent (G3013; Servicebio) on control and transfected cells samples (Huh-7 and HepG2). RNA quality was determined by spectrophotometry. Subsequently, cDNA was synthesized from the RNA templates with the SweScript All-in-One RT SuperMix Kit (G3337; Servicebio). Gene expression was quantified by qPCR using the 2$\times$ Universal Blue SYBR Green qPCR Master Mix (G3326; Servicebio) and the 2^{-Δ⁢Δ⁢CT} calculation method. The specific primers utilized are listed in Table 2.

Following trypsinization and resuspension in complete medium, control and transfected cells (Huh-7 and HepG2) were plated in 96-well plates at a density ranging from 5 $\times$ 10³ to 5 $\times$ 10⁴ cells per well. After a 24-h incubation under standard conditions (37 ℃, 5% CO₂), 10 µL of Cell Counting Kit-8 (CCK-8, G4103-5; Servicebio) reagent was added to each well. Then the plates were incubated for an additional 4 h, after which the absorbance at 450 nm was recorded using a microplate reader to assess cell viability.

The enzymatic activity of GAPDH in cultured control and transfected cells (Huh-7 and HepG2) were measured using a commercial GAPDH activity assay kit (BC2215; Solarbio) according to the manufacturer’s instructions [41]. GAPDH activity was expressed as units per mg of protein, where one unit was defined as the amount of enzyme that oxidizes 1 nmol of NADH per minute under the assay conditions. For measurements performed in a 96-well ultraviolet (UV) plate (light path = 0.6 cm), activity was calculated using Mathematical Components 3.

For IP₃ quantification, control and transfected cells (Huh-7 and HepG2) lysates were prepared, and a standard curve was established using serial dilutions of the supplied standard. The biotinylated detection antibody and HRP conjugate were diluted to 1$\times$ working solutions according to the manufacturer’s instructions. The enzyme-linked immunosorbent assay (ELISA) assay was conducted using a commercial kit (E-EL-0059; Elabscience, Wuhan, China) following the recommended protocol. Briefly, 50 µL of standards or cell lysates were added to a 48-well plate, followed by 50 µL of 1$\times$ biotinylated antibody. Following a 45-min incubation, the wells were washed three times with phosphate-buffered saline (PBS). Next, 100 µL of diluted HRP conjugate was added and incubated for 30 min, after which the plate was washed five times. Following the addition of 90 µL of substrate solution, the reaction was incubated for 30 min and then terminated with the stop solution. The absorbance was subsequently measured at 450 nm using a microplate reader [42].

For rescue and inhibition assays, Huh-7-WT and HepG2-WT cells were assigned to three groups: vehicle control (dimethyl sulfoxide [DMSO] (GC203006-10 mL; Servicebio), PLC inhibitor-treated (U73122, HY-13419; MedChemExpress) and inactive analog-treated (U73343, HY-108630; MedChemExpress), However ANXA5-KD (Huh-7-ANXA5 and HepG2-ANXA5) was another group. For pharmacological modulation, above three groups WT cells were pretreated with U73122 or U73443 (5 µM) for 30 min prior to subsequent analyses, whereas control cells received an equivalent volume of DMSO (0.1%). After treatment, intracellular IP₃ levels of four were quantified by ELISA according to the manufacturer’s instructions.

Control and transfected cells (Huh-7 and HepG2) were seeded into black 96-well clear-bottom plates at a density of 1 $\times$ 10⁵ cells per well and cultured for 24 h. After two to three washes with PBS, cells were loaded with 5 µM Rhod-4 AM (indicator used to detect intracellular Ca²⁺, MX4508-50UG; Maokangbio, Shanghai, China) premixed 1:1 (v/v) with 20% (w/v) Pluronic F-127 (MS4301-1G; Maokangbio), resulting in a final Pluronic concentration of 0.02%. Cells were incubated with this solution at 37 ℃ for 30 min, followed by three washes with Ca²⁺-free Ringer solution.

Fluorescence measurements were recorded using the SpectraMax iD3 multi-mode microplate reader (MD SpectraMax iD3; Molecular Devices, Shanghai, China) with a 540 $\pm$ 10 nm excitation filter and a 590 nm emission filter, acquiring signals every 2 s. Baseline fluorescence (F₀) was recorded for 0–60 s. At 60 s, ATP (D7378-2ml; Beyotime, Shanghai, China) (100 µM) in Ca²⁺-free Ringer solution (155 mM NaCl, 4.5 mM KCl, 3 mM MgCl₂, 10 mM D-glucose, 5 mM Na-HEPES) was applied to induce ER Ca²⁺ store depletion, and fluorescence was monitored from 60 to 154 s. Subsequently, 2 mM CaCl₂ was added (final [Ca²⁺] = 1 mM) to trigger SOCE, which was recorded from 154 to 430 s. The total recording duration was 0–600 s.

For each well, time-dependent fluorescence intensity (F_t) was collected, and intracellular Ca²⁺ changes were expressed as the ratio (Mathematical Components 4).

Cell migratory capacity was evaluated using a wound healing assay. WT, NC, and ANXA5-KD cells (Huh-7 and HepG2) were plated in 6-well plates at a density of 5 $\times$ 10⁵ cells per well and grown to complete confluence. A straight scratch was generated across the cell monolayer using a sterile 200 µL pipette tip. After gentle washing to remove detached cells, cultures were maintained in medium containing a reduced serum concentration. Images of identical wound regions were acquired immediately after scratching and again at 48 h. Cell migration was quantified by calculating the change in wound width with ImageJ software (National Institutes of Health, Bethesda, MD, USA).

Apoptosis was evaluated by flow cytometry using an Annexin V-FITC/PI apoptosis detection kit (BL107B; Biosharp, Guangzhou, China). To ensure inclusion of detached apoptotic cells, both culture supernatants and adherent cells were collected. Control and transfected cells (Huh-7 and HepG2) were detached with EDTA-free trypsin (G4002-100ML; Servicebio) to avoid nonspecific staining, followed by centrifugation at 2000 $\times$g for 5 min. After removal of the supernatant, cell pellets were washed twice with ice-cold PBS and resuspended in 500 µL of 1$\times$ binding buffer at a final concentration of 1 $\times$ 10⁶ cells/mL. Then cells were incubated with Annexin V-FITC and propidium iodide for 15 min at room temperature in the dark in accordance with the manufacturer’s instructions. Flow cytometry was conducted immediately after staining.

All animal procedures were reviewed and approved by the Animal Care and Use Committee of Bengbu Medical University (Approval No. [2025]514), and conducted in compliance with relevant institutional and national regulations. Fifteen male BALB/c nude mice (3 weeks old) obtained from Jiangsu Cavens Biotechnology (Shanghai, China) were randomly assigned to three groups (n = 5 per group): Huh-7-WT, Huh-7-NC, and Huh-7-ANXA5. Due to constraints in time and resources, the study was conducted exclusively with Huh-7 cells, representing a potential experimental limitation. For xenograft establishment, mice were anesthetized with isoflurane (induction at 3–4% and maintenance at 1.5–2% in oxygen, administered via inhalation; Hengfeng Qiang Biotechnology Co., Ltd., Jiangsu, China, batch No. 2401) and 5 $\times$ 10⁶ corresponding cells suspended in 100 µL of DMEM were injected subcutaneously into the right flank. Tumor growth and body weight were monitored daily. After a 3-week experimental period, animals were deeply anesthetized with isoflurane (4–5%, inhalation) and euthanized by cervical dislocation. Death was confirmed by cessation of respiration and heartbeat. Tumors were subsequently excised and weighed for further analyses.

Tumor and adjacent tissue specimens were fixed in 4% polyoxymethylene (G1101; Servicebio) at 37 ℃ with gentle agitation overnight. Following fixation, tissues were dehydrated, paraffin-embedded, and cut into sections with a thickness of 4–5 µm. Histological evaluation was carried out by hematoxylin and eosin (H&E) staining using a commercial H&E staining kit (G1076; Servicebio) in accordance with the manufacturer’s instructions.

Immunohistochemical staining was performed on paraffin sections following deparaffinization and antigen retrieval. To reduce nonspecific binding, sections were blocked in 10% bovine serum albumin at 37 ℃ for 1 h. Slides were subsequently incubated with primary antibodies at 4 ℃ overnight and then treated with the corresponding HRP-conjugated secondary antibodies for 30 min at 37 ℃. Signal development was carried out using the GTvision IHC detection system (Gene Tech Company Ltd., Shanghai, China) in accordance with the manufacturer’s protocol, with PBS washes performed between steps.

Statistical analyses were conducted using Origin 2022, ImageJ, and R Studio. Quantitative data obtained from qPCR and Western blot analyses are expressed as the mean $\pm$ standard deviation. Differences between the two groups were evaluated using a two-tailed Student’s t-test, whereas one-way analysis of variance was applied for comparisons involving more than two groups. Statistical significance was defined as p 0.05.

The bioinformatics analysis workflow of this study is summarized in Fig. 2A. To explore the involvement of SRGs in HCC, we compared their expression profiles between tumor and adjacent normal tissues. Among the 106 differentially expressed SRGs (FDR 0.05, $|$logFC$|$ $>$1), 92 were upregulated in tumor samples, and 14 were downregulated (Fig. 2B). To assess the prognostic relevance of these genes, univariate Cox regression analysis was conducted. The resulting forest plot (Fig. 2C) highlighted 32 genes significantly correlated with OS in patients with HCC, which were subsequently selected for model construction. Of these, 28 genes functioned as risk factors (HR $>$1). Then, consensus clustering analysis was conducted to further validate the prognostic impact of these genes. The CDF curve, evaluated for K = 2 to K = 9, indicated that K = 2 was the optimal number to define two distinct subgroups (Fig. 2D). The consensus matrix at K = 2 confirmed clear separation between the two clusters with high consensus, which demonstrated clear segregation between the two clusters (Fig. 2E). The clinical relevance of these two gene clusters was assessed using Kaplan–Meier survival analysis, which revealed a significant difference in OS between the clusters (p 0.001). Patients in Cluster B exhibited a more favorable prognosis than those in Cluster A (Fig. 2F).

Fig. 2.

DEG heatmap to illustrate prognostic genes. (A) Bioinformatics analysis workflow of this study. (B) Water blue represents normal sample tissues; pink represents tumor tissues; and red and blue in the DEG heatmap correspond to high and low gene expression, respectively. (C) Prognostic forest plot with color HRs: blue indicates genes with protective roles (HR 1); red indicates genes with risk roles (HR $>$1). (D) The consensus CDF plot showed optimal cluster stability at K = 2. (E) Cluster separation was most distinct at K = 2. (F) Kaplan–Meier curve showing that group B had significantly better OS than group A, with group B shown in orange and group A in blue. DEG, differentially expressed gene; HRs, hazard ratios; CDF, cumulative distribution function; OS, overall survival.

In summary, differential expression and univariate Cox regression analyses identified a subset of SRGs associated with HCC prognosis, including 28 risk genes. Consensus clustering based on these prognostic SRGs divided patients into two distinct subgroups with significantly different OS, indicating that these gene clusters may serve as robust indicators of clinical outcomes in HCC.

To optimize the gene set, minimize overfitting, and select the most informative prognostic markers, we applied the LASSO regression. Based on the optimal penalty parameter ($\lambda$) determined by the minimum criteria, a 17-gene prognostic signature was established (Fig. 3A,B). The individual risk scores for patients were calculated using a formula combining the LASSO coefficients with the expression levels of these 17 genes.

Fig. 3.

Identification of risk signature. (A) LASSO regression of the OS-related genes. (B) Cross-validation for tuning the parameter selection in the LASSO regression. LASSO, least absolute shrinkage and selection operator.

In summary, LASSO regression enabled the identification of a concise 17-gene signature, and the derived risk scores provide a quantitative measure of patient prognosis in HCC.

To evaluate the prognostic value of the refined SRG signature, a risk score was calculated for each patient in the TCGA and GEO cohorts. Patients were stratified into low- and high-risk groups based on the median risk score. Kaplan–Meier analysis revealed a significant difference in OS, with patients in the low-risk group exhibiting markedly longer survival times than those in the high-risk group (p 0.001; Fig. 4A). To assess the robustness of the prognostic model, the cohort was randomly split into a training set and testing set. In the training cohort, survival analysis revealed a significantly better outcome for patients classified as low risk (p 0.001; Fig. 4B). The predictive performance of the risk score was further evaluated using time-dependent ROC analysis. In the whole cohort, the AUC for predicting OS was 0.755 at 1 year, 0.692 at 3 years, and 0.723 at 5 years (Fig. 4C). Similarly, in the training set, the AUC values were 0.821, 0.720, and 0.749 for 1-, 3-, and 5-year survival, respectively (Fig. 4D).

Fig. 4.

Validation of risk signature in TCGA and GEO cohort. (A) Kaplan–Meier curves for the entire cohort show blue (low-risk) and red (high-risk) groups (p 0.001). (B) Kaplan–Meier curves for the training group use the same label color as in (A) (p 0.001). (C) ROC curve of the full cohort depicts 1-year survival in green, 3-year survival in blue, and 5-year survival in red. (D) In the training cohort ROC curve, the label color is the same as in (C) (p 0.001). TCGA, The Cancer Genome Atlas; GEO, Gene Expression Omnibus; ROC, receiver operating characteristic.

Overall, the SOCE-related prognostic signature successfully categorized patients into distinct risk groups that exhibited significant differences in clinical outcomes.

Multivariate Cox regression analysis was conducted to identify which genes within the 17-gene LASSO-derived prognostic signature had independent prognostic significance. This analysis identified nine genes as factors with independent prognostic value. A heatmap illustrating the expression patterns of these nine genes was created to visualize their relationship with patient outcomes (Fig. 5A).

Fig. 5.

Refinement of the prognostic signature and construction of a nomogram. (A) The risk heatmap shows nine genes. (B) Forest plot illustrating independent prognostic factors. (C) The nomogram incorporates four risk factors: age, T stage, N stage, and risk score. A specific point value is assigned to each factor on its corresponding scale. The total points are the sum of all four factors. A vertical line drawn from the total points axis to the outcome axes provides approximate predictions of the patient’s 1-, 3-, and 5-year OS probabilities, as illustrated by the red vertical line in the figure. (D) The calibration curve of the nomogram shows that the green, blue, and red lines represent 1-, 3-, and 5-year survival probability calibrations, respectively.

Subsequently, the relationship between these nine genes and key clinicopathological factors was evaluated using multivariate analysis, and the results were visualized as a forest plot (Fig. 5B). The analysis revealed significant correlations between these genes and tumor T stage and N stage (p 0.001 and p 0.05, respectively), suggesting a potential role in local tumor invasion and progression. No statistically significant associations were observed with age or sex (p $>$ 0.05). Age was retained in the subsequent model to reflect its established clinical relevance to cancer incidence, whereas sex was excluded due to lack of statistical significance, resulting in a more concise predictive tool.

By integrating the remaining variables, an individualized prognostic nomogram was generated to predict 1-, 3-, and 5-year OS probabilities (Fig. 5C). Curve analysis demonstrated strong concordance between the model’s predicted and actual survival outcomes (Fig. 5D). Closer examination of the heatmap focused on risk-associated genes, and found that six (rearranged during transfection [RET], secreted phosphoprotein 1/osteopontin [SPP1], vitamin D receptor [VDR], annexin A5 [ANXA5], menin 1 [MEN1], and GNAS complex locus [GNAS]) were upregulated in the high-risk group, marking them as candidates for further analysis.

In summary, multivariate Cox regression identified nine genes with independent prognostic value, and their expression patterns were associated with key clinicopathological features such as tumor T and N stages. An individualized nomogram incorporating these variables reliably predicted 1-, 3-, and 5-year OS, with calibration curves confirming strong agreement between predicted and observed outcomes. Among the risk-associated genes, six (RET, SPP1, VDR, ANXA5, MEN1, and GNAS) were notably upregulated in the high-risk group, highlighting their potential relevance in HCC progression and patient prognosis.

To systematically investigate the enrichment patterns of the aforementioned six high-risk genes within the tumor immune microenvironment, we conducted a detailed analysis of each gene using the TISCH2 database, with representative results illustrated using the HCC dataset LIHC-GSE140228 (Fig. 6A). The analysis revealed that the key gene ANXA5 is highly expressed across multiple immune cell types, with notable enrichment particularly in dendritic cells and macrophages (Fig. 6B,C), suggesting its potential involvement in antigen presentation and innate immune regulation. In contrast, the other five genes showed generally low expression levels across various immune cells without significant enrichment. Although GNAS exhibited relatively high expression in certain cell populations, its expression pattern lacked cell type specificity and appeared diffuse, and thus was not pursued further. Based on these findings and established evidence linking ANXA5 to tumor progression, metastasis, and poor clinical outcomes in diverse cancers—including HCC—we selected it as a key gene for subsequent mechanistic investigation. Studies have shown that overexpressed ANXA5 in HCC is linked to malignant phenotypes and poor prognoses in in vitro and in vivo models [14, 43, 44, 45]. Furthermore, high ANXA5 expression is significantly correlated with aggressive clinicopathological features and poor survival in various solid tumors [46, 47, 48, 49]. To further elucidate the biological functions and regulatory networks of ANXA5, we constructed a protein–protein interaction network. Our analysis demonstrated a high-confidence interaction between ANXA5 and GAPDH (Fig. 6D), indicating their potential synergy in cellular metabolic regulation and immune responses. This finding provides a critical molecular clue for future functional studies.

Fig. 6.

Immune cell distribution map in LIHC-GSE140228. (A) Example of an immune cell distribution map; different colors represent different immune cells. (B) The map shows the distribution of the key gene ANXA5, with color intensity representing the degree of enrichment. (C) Immune cell expression of the key gene ANXA5, with column height representing enrichment degree. (D) The protein–protein interaction network demonstrated a strong interaction between ANXA5 and GAPDH. Node size represents the number of interactions, with larger nodes indicating greater connectivity. Since most proteins (except F3) had similar interaction counts, size differences were minimal. The most connected nodes, ANXA5 and GAPDH, are marked in red. Panels (A–C) were generated using the TISCH2 database (https://tisch.compbio.cn/gallery/). Panel (D) was constructed based on data obtained from the STRING database (version 12.0; https://cn.string-db.org/) and visualized using Cytoscape (v3.9.0).

In summary, while ANXA5 expression was detected in immune cell subsets within the tumor microenvironment—a topic for future research—this study primarily focused on investigating its functional role in key HCC malignant phenotypes, including proliferation, invasion, and metastatic potential.

To assess transfection efficiency, ANXA5 expression was evaluated in WT, NC, and three ANXA5-targeting shRNA groups (shRNA-ANXA5(h)-1, shRNA-ANXA5(h)-2, shRNA-ANXA5(h)-3) using Western blotting and qPCR. The results demonstrated that shRNA-ANXA5(h)-1 was the most effective sequence for KD (Fig. 7A–C), whereas no significant difference was observed between the WT and NC groups (p $>$ 0.05). Consequently, we renamed shRNA-NC and shRNA-ANXA5(h)-1 to NC and ANXA5-KD, respectively, for simplicity.

Fig. 7.

Western blotting and qPCR analysis of ANXA5 and GAPDH expression. (A) Western blot bands of ANXA5. (B) Quantitative analysis of ANXA5 protein expression based on Western blot band intensity (p 0.05). (C) Relative mRNA expression level of ANXA5 measured by qPCR (p 0.05). (D) Western blot bands of GAPDH. (E) Quantitative analysis showing no significant difference in GAPDH protein expression among groups (p $>$ 0.05). (F) Relative mRNA expression level of GAPDH measured by qPCR, with no statistically significant differences between groups (p $>$ 0.05). (i) Huh-7, (ii) HepG2. Statistical significance is indicated as follows: *, p 0.05; **, p 0.01; ***, p 0.001; ns, not significant. Western blot bands were imaged using a chemiluminescence imaging system (SCG-W3000, Servicebio, Wuhan, China). Quantitative bar graphs were generated using GraphPad Prism 10.

Western blotting and qPCR analyses across WT, NC, and ANXA5-KD groups consistently demonstrated that ANXA5 KD did not significantly alter GAPDH expression levels (p $>$ 0.05; Fig. 7D–F).

To determine whether the interaction between ANXA5 and GAPDH depends on ANXA5 expression levels, co-IP assays were conducted in WT, NC, and ANXA5-KD cells. Input analysis confirmed comparable GAPDH expression levels across all groups, whereas ANXA5-KD cells exhibited a significant reduction in ANXA5 levels. Co-IP using an anti-ANXA5 antibody followed by immunoblotting for GAPDH revealed a notable ANXA5–GAPDH interaction in WT and NC cells, which was substantially diminished in the ANXA5-KD group. No GAPDH signal was detected in the IgG control IP (Fig. 8A). Quantitative analysis of input samples showed no significant differences in GAPDH expression among the groups (p $>$ 0.05), whereas ANXA5 levels were significantly reduced in the ANXA5-KD group (p 0.05; Fig. 8B). Quantitative analysis of IP’d samples revealed a significant difference in GAPDH levels between the ANXA5-KD group and both the WT and NC groups (p 0.05), with no detectable GAPDH signal in the IgG control group (Fig. 8C).

Fig. 8.

Interaction of ANXA5 with GAPDH by Co-IP. (A) Input bands showed comparable GAPDH expression across all groups, whereas ANXA5 expression was markedly reduced in ANXA5-KD cells. IP with an anti-ANXA5 antibody followed by immunoblotting for GAPDH revealed a strong ANXA5–GAPDH interaction in WT and NC cells, which was substantially diminished upon ANXA5-KD, whereas no GAPDH signal was detected in the IgG IP control. (B) Quantitative analysis showed no significant difference in GAPDH protein expression among input sample groups (p $>$ 0.05), whereas there was a significant difference in ANXA5 protein expression in the ANXA5-KD group (p 0.05). (C) Quantitative analysis of the ANXA5-KD group showed significant differences in GAPDH signals among IP-ANXA5 groups (p 0.05), with no signal detected in the IP-IgG control. (i) Huh-7, (ii) HepG2. Co-IP, Co-immunoprecipitation; KD, knockdown; WT, wild-type. Statistical significance is indicated as follows: ***, p 0.001; ns, not significant. Co-IP bands were imaged using a chemiluminescence imaging system (SCG-W3000, Servicebio, Wuhan, China). Quantitative bar graphs were generated using GraphPad Prism 10.

In summary, these findings indicate that ANXA5 forms a specific complex with GAPDH in HCC cells in an ANXA5-dependent manner. Importantly, ANXA5 depletion markedly weakened the ANXA5–GAPDH interaction without altering total GAPDH protein levels, suggesting that ANXA5 regulates GAPDH primarily via protein–protein interaction rather than expression control. This observation provides a mechanistic basis for subsequent analyses exploring how ANXA5 modulates GAPDH enzymatic activity and downstream metabolic and Ca²⁺ signaling pathways.

To systematically investigate the role of ANXA5 in the malignant behavior of HCC cells, functional assays including wound healing, CCK-8 proliferation, and flow cytometry-based apoptosis analysis were performed.

In the wound healing assay, ANXA5-KD cells exhibited significantly slower wound closure compared with WT and NC cells (Fig. 9A). Quantitative analysis confirmed a marked reduction in cell migration rate in the KD group (p 0.05; Fig. 9B), indicating that ANXA5 depletion effectively inhibited HCC cell motility.

Fig. 9.

ANXA5 regulates HCC cells in vitro. (A) Micrographs of wound healing assays: (i) Huh-7, (ii) HepG2. (B) Quantification of cell migration in the wound healing assay: (i) Huh-7, (ii) HepG2. (C) CCK-8 assay quantifying cell proliferation: (i) Huh-7, (ii) HepG2. (D) Flow cytometry scatter plots of apoptosis detection. The quadrants represent: lower left (viable cells), upper left (mechanically damaged cells), upper right (late apoptotic cells), and lower right (early apoptotic cells); (i) Quantitative analysis of Huh-7-WT, (ii) Huh-7-NC, (iii) Huh-7-ANXA5, (iv) Huh-7 cells. (E) Flow cytometry scatter plots of apoptosis detection. (i) Quantitative analysis of HepG2-WT, (ii) HepG2-NC, (iii) HepG2-ANXA5, (iv) HepG2 cells. Statistical significance is indicated as follows: **, p 0.01; ***, p 0.001; ns, not significant. Scale bar = 50 µm. Wound healing images were captured using a microscope (model: XSP-C204), and apoptosis data were analyzed using FlowJo (v10.8.1). Quantitative bar graphs were generated using GraphPad Prism 10.

Cell proliferation was assessed using the CCK-8 assay, which showed that absorbance values in the ANXA5-KD group showed significantly decreased levels relative to the control groups (p 0.05; Fig. 9C), suggesting substantially suppressed proliferative activity and highlighting the role of ANXA5 in maintaining HCC cell proliferation.

Flow cytometry with Annexin V-FITC/PI staining was used to measure the apoptosis level. Scatter plots revealed a notable increase in the proportions of apoptotic cells in the ANXA5-KD groups (Fig. 9D). Quantitative analysis confirmed that the proportion of apoptotic cells was markedly increased in the KD group relative to the control groups (p 0.05; Fig. 9E), indicating a pronounced pro-apoptotic effect of ANXA5 KD.

Together, these findings establish that ANXA5 plays a critical role in regulating migration, proliferation, and apoptosis in HCC cells. ANXA5 KD led to reduced cell migration and growth while significantly increasing apoptosis, highlighting its potential as a therapeutic target for HCC.

Based on prior evidence indicating a strong interaction between ANXA5 and GAPDH, and having confirmed that ANXA5 does not alter GAPDH expression, we investigated whether ANXA5 influences GAPDH enzymatic activity. An ELISA-based assay revealed a positive correlation between ANXA5 levels and GAPDH activity (p 0.05; Fig. 10A).

Fig. 10.

ANXA5 modulates Ca²⁺ signaling through the GAPDH/IP₃ pathway. (A) GAPDH enzyme activity was determined by spectrophotometrically monitoring the rate of NADH consumption during the reaction. (B) Measurement of IP₃ level. (C) IP₃ rescue experiment. There were no significant differences between the vehicle and U73443 groups (p $>$ 0.05), unlike the comparison between the ANXA5-KD and U73122 groups (p 0.05). The rescue assay shows that U73122 (a PLC inhibitor) and ANXA5-KD similarly impair IP₃ production. (D) In the Ca²⁺ oscillation curve, the first peak represents ER Ca²⁺ release, whereas the second peak corresponds to the SOCE. The first dashed line indicates the time point of ATP stimulation, and the second dashed line marks the addition of extracellular Ca²⁺. (E) AUC of ER Ca²⁺ release and SOCE peak. (i) Huh-7, (ii) HepG2. IP₃, inositol 1,4,5-trisphosphate; PLC, phospholipase C; AUC, area under the curve. Statistical significance is indicated as follows: *, p 0.05; **, p 0.01; ***, p 0.001; ns, not significant. This panel was generated using GraphPad Prism 10.

To investigate the potential mechanism, we examined whether ANXA5 may influence PIP₂ metabolism through GAPDH enzymatic activity. Since direct quantification of PIP₂ was technically challenging, we measured the levels of its downstream metabolite, IP₃, as an indirect indicator of PIP₂ hydrolysis. ELISA analysis showed a significant decrease in IP₃ concentration following ANXA5 knockdown (p 0.05; Fig. 10B).

To investigate whether ANXA5 regulates intracellular IP₃ levels via the PLC pathway, we performed pharmacological inhibition and genetic KD experiments. Cells were divided into four groups: vehicle control, ANXA5-KD, U73122, and U73343 (Fig. 10C). Measurement of intracellular IP₃ revealed that U73122 treatment and ANXA5-KD both significantly reduced IP₃ levels compared with vehicle control (p 0.05), whereas U73343 had no significant effect (p $>$ 0.05). Notably, U73122 exhibited the strongest inhibition of IP₃ production, consistent with its role as a PLC inhibitor. These results indicate that ANXA5 contributes to IP₃ generation in HCC cells, likely via the PLC pathway, and that pharmacological PLC inhibition mimics the effect of ANXA5 depletion, supporting the proposed ANXA5/GAPDH/IP₃/SOCE working model.

To further examine the functional impact, cytosolic Ca²⁺ dynamics were assessed using a fluorescent Ca²⁺ indicator in Huh-7 and HepG2 cells with ANXA5-KD. Ca²⁺ oscillation curves showed that ANXA5-KD attenuated both ER Ca²⁺ release and the SOCE peak compared with WT and NC groups (Fig. 10D). The AUC for ER release and SOCE peak was quantified from the Ca²⁺ oscillation curves (Fig. 10E), confirming that decreased ANXA5 expression compromised ER Ca²⁺ release and consequently diminished SOCE (p 0.05).

In summary, our findings suggest that ANXA5 may positively influence GAPDH enzymatic activity without affecting its expression levels. Moreover, ANXA5 knockdown was associated with a marked reduction in the downstream metabolite IP₃, which may reflect indirect alterations in PIP₂ metabolism. Pharmacological inhibition of PLC with U73122 mimicked the effects of ANXA5 depletion, supporting the involvement of the PLC pathway in ANXA5-mediated IP₃ generation. Functionally, ANXA5 deficiency impaired intracellular Ca²⁺ signaling, as evidenced by diminished ER Ca²⁺ release and reduced SOCE in Huh-7 and HepG2 cells. Collectively, these results suggest that ANXA5 regulates HCC cell Ca²⁺ dynamics through GAPDH-dependent modulation of the PLC/IP₃ axis, consistent with the proposed ANXA5/GAPDH/IP₃/SOCE working model.

To evaluate the in vivo function of ANXA5, nude mice were subcutaneously inoculated with WT, NC, or ANXA5-KD cells (Huh-7) in the right flank. Measurement of resected tumors showed that knocking down ANXA5 significantly reduced tumor volumes compared to the WT and NC groups (Fig. 11A,B). Quantitative analysis further demonstrated that ANXA5 KD markedly reduced the rate of tumor progression (Fig. 11C) and tumor burden at sacrifice (Fig. 11D), and was associated with a slower decline in mouse body weight (Fig. 11E), suggesting a less aggressive tumor burden.

Fig. 11.

ANXA5 promotes tumor growth and malignancy in vivo. (A) Images of tumors in situ. (B) Excised tumor specimens from each group. (C) Tumor growth curves over time. (D) Bar graph showing final tumor mass. (E) Body weight changes of mice during the study. (F) Histological, cross-sectional analysis (top to bottom respectively) of Huh-7-WT, Huh-7-NC, and Huh-7-ANXA5 cells; (left to right respectively) H&E staining of Ki-67, vimentin, and CD206. Scale bar = 50 µm. (G) Quantification of Ki-67-positive cells. (H) Quantification of vimentin expression. (I) Quantification of CD206-positive cells. H&E, hematoxylin and eosin. Statistical significance is indicated as follows: **, p 0.01; ***, p 0.001; ns, not significant. Line graphs and bar charts in this panel were generated using GraphPad Prism 10, and HE and IHC images were acquired using CaseViewer.

Histological examination using H&E staining showed that tumors from the control groups exhibited more compact and dense architecture compared with the more loosened structure observed in ANXA5-KD tumors (Fig. 11F). Immunohistochemical analysis of Ki-67, a proliferation marker, and vimentin, a marker of epithelial–mesenchymal transition (EMT), revealed substantially higher expression levels in the WT and NC groups than in the ANXA5-KD group. Additionally, staining for the M2 macrophage marker CD206 showed denser positive signals in control tumors. These findings suggest a potential association between ANXA5 expression and M2 macrophage–related signals within the tumor microenvironment (Fig. 11F). Quantitative analysis of staining intensity confirmed significant reductions in proliferation, malignancy, and M2 marker expression upon ANXA5 KD (p 0.05; Fig. 11G–I). Collectively, these findings indicate that ANXA5 contributes to tumor growth, malignancy, and the M2-polarized tumor microenvironment in vivo.

In summary, the in vivo experiments demonstrated that ANXA5 depletion significantly suppressed tumor growth and progression. Tumors derived from ANXA5-KD cells exhibited smaller volumes, slower growth kinetics, and reduced tumor burden at the endpoint compared with WT and NC controls. Histological assessment revealed a less compact tumor architecture in ANXA5-KD tumors, accompanied by decreased Ki-67 and vimentin expression, indicating diminished proliferation and EMT activity. Furthermore, reduced CD206 staining suggested a decrease in M2-associated macrophage signals within the tumor tissue. Taken together, these findings indicate that ANXA5 expression is associated with tumor growth, malignant phenotypes, and changes in immune-related markers in vivo.