SLC7A11, a disulfidptosis-related gene, correlates with multi-omics prognostic analysis in hepatocellular carcinoma

Abstract

Background

This study sought to establish a risk score signature based on disulfidptosis-related genes (DRGs) to predict the prognosis of hepatocellular carcinoma (HCC) patients.

Methods

The expression data of DRGs from the Cancer Genome Atlas (TCGA) and the International Cancer Genome Consortium (ICGC) was analyzed to develop and validate a DRG prognostic signature (DRGPS). In vitro, experiments were conducted to explore DRG expressions and roles in HCC tissues and cell lines. HCC tissue microarrays were employed to analyze SLC7A11 expression and its association with clinicopathological characteristics.

Results

The DRGPS consisted of 5 DRGs (SLC7A11, MATN3, CLEC3B, CCNJL, and PON1). The survival rate of HCC patients in high-risk group was significantly lower than that in low-risk group. The DRGPS was also associated with the modulation of tumor microenvironment (TME), tumor mutation burden (TMB), stemness and chemosensitivity. Furthermore, pan-cancer analysis suggested that the DRGPS risk score was associated with immune infiltration and stemness in multiple cancers. Moreover, our DRGPS had potential for predicting treatment efficacy in HCC patients. Finally, we confirmed that downregulation of SLC7A11, a DRG, inhibited the proliferation and migration of HCC cells, while its high expression correlated with advanced TNM clinical stage and larger tumor size.

Conclusions

This study systematically describes a novel DRGPS constructed for predicting HCC prognosis, providing a new approach to risk stratification and treatment options. It also investigates the expression and function of SLC7A11, contributing to further exploration of the molecular mechanism underlying disulfidptosis in HCC, as well as its prognostic and therapeutic implications.

Graphical Abstract

Introduction

Hepatocellular carcinoma (HCC) represents the major histologic subtype of liver cancer, a third leading cause of cancer-related mortality in the world, accounting for approximately 80% of all primary liver cancer cases [1]. The heterogeneity of HCC renders conventional chemotherapy ineffective, resulting in a 5-year survival rate of only 18% [2, 3]. Hence, it is imperative that molecular mechanisms of HCC pathogenesis and development be explored and that novel early detection and treatment targets be identified for HCC patients.

There is a genetically regulated form of cell death called programmed cell death (PCD), which is mediated by the coordinated action of a variety of cellular signaling pathways. Recently, a growing number of PCD mechanisms, such as autophagy, ferroptosis, and cuproptosis, have been elucidated and substantiated as influential contributors to cancer initiation and progression [4, 5]. Likewise, exploiting the PCD mechanism of tumor cells has been proposed as a potentially beneficial anti-tumor strategy [6]. Disulfidptosis, a novel form of cell death triggered by glucose deprivation, results from the accumulation of disulfide stress and exhaustion of intracellular nicotinamide adenine dinucleotide phosphate (NADPH) [7, 8]. Through the formation of actin-rich protrusions, tumor cells migrate and invade, promoting tumor epithelial mesenchymal transformation (EMT) and metastasis [9]. Significantly, SLC7A11 has been documented to be a key player in disulfidptosis, and it is also involved in cystine uptake, antioxidative stress resistance, and a protective mechanism against ferroptosis, which is highly expressed in HCC and negatively impacts survival [10, 11]. A further finding is that glucose-6-phosphate dehydrogenase (G6PD) and transketolase (TKT), critical enzymes for NADPH production in disulfidptosis, are significantly up-regulated in HCC [12]. Nonetheless, little is known regarding the relationship between disulfidptosis and HCC.

As of now, therapeutic options for HCC are constrained. In individuals suffering from advanced stages, immunotherapy and chemotherapy remain the most appropriate treatment options [13]. It is imperative to note, however, that the tumor microenvironment (TME) exhibits significant immunosuppressive effects, posing challenges of chemoresistance and immune evasion in targeted and immunotherapeutic approaches for HCC [14]. Moreover, cancer stem cells (CSCs), which possess the characteristics of stem cells within tumors, play a critical role in the development of chemoresistance [15]. It is also believed that CSCs facilitate immune evasion through interactions with TME [16]. Recent investigations have endeavored to improve the therapeutic efficacy against HCC by inducing glucose deprivation to facilitate the demise of CSCs [17], highlighting the potential involvement of disulfidptosis in the stemness of HCC. Moreover, a novel prognostic model also indicates that a high hepatic arterial blood supply is significantly associated with poorer overall survival (OS) in HCC patients [18]. Nevertheless, there are still unmet needs for genetic markers to monitor prognostic outcomes and direct treatment options for HCC patients [19, 20].

In this study, we conducted a clustering analysis using the Cancer Genome Atlas–Liver Hepatocellular Carcinoma (TCGA–LIHC) cohort, based on disulfidptosis-related genes (DRGs). In addition, we analyzedthe expression profile of DRGs to idecipher their relationship between tumor stemness and therapeutic sensitivity. Aside from this, a prognostic signature risk score was devised as part of optimizing treatment and stratifying patients, with the potential to serve as a suitable biomarker for better molecular characterization of the TME in HCC. These findings may yield novel insights into personalized prognosis prediction and precise therapeutic strategies for HCC patients.

Materials and methods

Acquisition of data

Based on the TCGA–LIHC (https://portal.gdc.cancer.gov), this study retrieved the expression of genes in patients, somatic mutation profile, copy number variation (CNV) data and their matching clinical traits, which consists of 371 tumor samples and 50 controls. Furthermore, raw count data of a Japanese HCC population, the ICGC–LIRI–JP cohort (including 231 tumor samples), were acquired from the International Cancer Genome Consortium (ICGC) (https://dcc.icgc.org/projects/LIRI-JP). Further analysis did not include patients without survival information.

Evaluation of DRG clusters

Genes identified in the literature that are involved in the disulfidptosis pathway are considered DRGs [7]. The DRGs were clustered using the “ConsensusClusterplus” R package, from which the number and consistency of clusters were built. Meanwhile, differentially expressed genes (DEGs) within DRG clusters A and B were identified using the R package “limma” based on the following criteria of |logFC|> 1, and adjusted P < 0.01.

Analysis of functional traits

Analyses of the Kyoto Encyclopedia of Genes and Genomes (KEGG) and gene ontology (GO) were conducted using “clusterProfiler” and “ggplot2” R package. From MSigDB, the reference gene set “c2.cp.kegg.v7.4.symbols.gmt” in conjunction with the R package “GSVA” were selected to conduct Gene Set Variation Analysis (GSVA) (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp). Adjusted P value < 0.05 was set as significant screen criteria.

Development of the DRG prognostic signature (DRGPS)

DRGPS was established to assess the prognosis of HCC patients. Univariate cox regression analysis was performed to determine the impact of DEGs associated with DRG clusters on HCC survival (P < 0.05). After incorporating the gene risk data, the Least Absolute Shrinkage and Selector Operator (LASSO) analysis was conducted to eliminate the possibility of overfitting. Finally, a DRGPS was then established through a multivariate cox regression analysis. For the calculation of the DRGPS risk score, a formula was requisite, as delineated below: ${\text{risk score}} = \mathop \sum \limits_{{{\text{i}} = 1}}^{{\text{n}}} \beta {\text{i*}}\left( {\text{expression of DRGi}} \right)$, where n denotes the amount of DRGs associated with the construction of the DRGPS, and β means the regression coefficient of DRGs. The median risk scores of DRGPS were used to distinguish low-risk and high-risk patients with HCC by the “survminer” package. The constructed DRGPS was validated by the ICGC cohort, in addition to the TCGA–LIHC cohort. The survival outcomes of the risk subgroups were compared using Kaplan–Meier survival analysis, and the prognostic performance of DRGPS was evaluated using the Area Under Curves (AUCs), Receiver Operating Characteristic Curves (ROCs) and decision curve analysis.

Clinical correlation and stratification analysis of DRGPS

Cox regression analysis was performed from univariate and multivariate perspectives to determine if the DRGPS risk score can serve as an independent prognostic indicator of HCC patients. Based on clinicopathological strata, Kaplan–Meier analysis was carried out to evaluate the DRGPS risk score. Through the use of the R packages “survival” and “rms”, a variety of clinicopathological factors were examined to develop a predictive nomogram for survival outcomes at 1 year, 3 years, and 5 years. An analysis of calibration plots illustrated how predicted results differed from actual results.

Infiltration of immune cells

To estimate the specific types of immune cells in DRG A and B clusters, the CIBERSORT algorithm was employed for our analysis. The ESTIMATE algorithm was also applied to compute immune scores, stromal scores, and ESTIMATE scores. Furthermore, the single-sample gene set enrichment analysis (ssGSEA) algorithm was utilized to compare immune functions within the HCC TME between the two risk groups. For the purpose of comparing the immunological landscape between the two risk groups, seven algorithms were employed, including EPIC, CIBERSORT, QUANTISEQ, CIBERSORT–ABS, ESTATE, TIMER and MCPcounter. In addition, different immune checkpoint genes (ICGs) were selected to determine if their expression levels differed between high- and low-risk individuals.

An assessment of the tumor mutation burden (TMB) and the tumor stemness

Using the “maftools” R package, somatic mutations in HCC patients from high- and low-risk groups were analyzed with mutation annotation formats (MAF) compiled from the TCGA database. We scored each HCC sample according to the number of nonsynonymous mutations/exon length (35 million) and performed a survival analysis. In addition, stemness indices mRNAsi and mDNAsi were calculated for each HCC sample.

Differences in chemotherapeutic efficacy

“pRRophetic” and “oncoPredict” R package from the Cancer Therapeutics Response Portal (CTRP) and Genomics of Drug Sensitivity in Cancer (GDSC; http://www.cancerrxgene.org/) were used to predict chemotherapy response for HCC patients. The half-maximum inhibitory concentrations (IC50s) for each sample were determined using a ridge regression algorithm. The Wilcoxon signed-rank test was conducted to determine whether the high- and low-risk groups were sensitive to chemotherapeutic drugs differently.

Tissue samples, cell lines and informed consent

At Xiangya Hospital, Central South University, we collected specimens of HCC and normal paracancerous tissues to further verify the expression of DRGs in our signature. Prior to surgery, none of the patients had received anti-tumor treatment. It was objectively impossible to obtain informed consent considering that the experiment was performed using paraffin-embedded blocks from the Pathology Department. Therefore, an informed waiver was requested and approved by Central South University’s Ethics Committee (No.202303036). Human normal liver cells (L02) and HCC cell lines (HepG2 and Hep3B were generously provided by the Cancer Research Institute of Central South University (Hunan, China). The cells were cultured in Dulbecco’s modified Eagle’s medium (Hyclone, Logan, UT, USA) supplemented with 10% fetal bovine serum (FBS) and maintained at 37 ℃ in a 5% CO₂ atmosphere. The culture medium was added with 100 mg/ml streptomycin and 100 U/ml penicillin for optimal growth conditions.

Quantitative reverse transcription polymerase chain reaction (qRT-PCR)

Cell and tissue RNA was isolated by TRIzol Reagent (Invitrogen, Carlsbad, USA). The initial cDNA synthesis using 1 μg of total RNA was achieved through the Reverse Transcription Kit (TransGen Biotech, Beijing, China). A qPCR analysis was then performed with the SYBR-Green Master Mix (TransGen Biotech, Beijing, China) in the following conditions: 30 s at 94 ℃, followed by 40 cycles of 5 s at 94 ℃, 15 s at 60 ℃, and 10 s at 72 ℃. Data were normalized using GAPDH as a reference. Please refer to the Supplementary Table 1 for details on these primer sequences.

Immunohistochemistry (IHC)

IHC was conducted on a tissue microarray (OD-CT-DgLivT10-024, Shanghai Outdo Biotech CO., LTD., China) consisting of 80 pairs of HCC and normal paracancerous tissues. The tissues were incubated overnight at 4 ℃ with an anti-SLC7A11 antibody (1:500, HA601071, Huaan Biotechnology Co., Ltd., China). After applying the appropriate secondary antibody, the labeled antigen was visualized using diaminobenzidine. The expression of SLC7A11 was scored semiquantitatively for positive cell ratio (0, 0–5%; 1, 6–25%; 2, 26–50%; 3, 51% ~ 75%; or 4, > 75%) and positive intensity (0, no positive; 1, weak; 2, moderate; or 3, strong). The final score was defined as positive cell ratio multiplied by positive intensity.

Cell transfection

The cells were seeded in 6-well culture plates, and cell transfection was performed when the concentration of the cells reached 40–60%. For SLC7A11 downexpression, the small interfering RNA (siRNA) targeting SLC7A11 and corresponding negative control were synthesized from Ribobio (Guangzhou, China). The cells were transfected by Lipo2000™ (Invitrogen,Carlsbad, CA) following the manufacturer’s protocol. Furthermore, 24 h after transfection, the cells were collected for cell proliferation assays and wound healing assay. The sequence of siSLC7A11 would be found in the Supplementary Table 1.

Cell proliferation assay

A total of 1 × 10³ transfected cells were plated in 96-well plates and incubated for 12, 24, 48, and 72 h. To assess cell viability, the Cell Counting Kit 8 (CCK-8, APExBIO, Houston, USA) was utilized, with 10 μl of CCK-8 added to each well. Following a 1-h incubation at 37 ℃, the absorbance was measured at a wavelength of 450 nm.

Migration assay

A total of 2 × 10⁵ cells were seeded in the upper chamber using 200 µl of serum-free medium, while 600 µl of medium containing 10% FBS was added to the lower chamber. The cells were incubated at 37 ℃ with 5% CO₂ for 48 h. After incubation, non-penetrating cells were gently removed from the membrane using a cotton swab, while migrated or invaded cells were fixed with 0.1% crystal violet.

Wound healing assay

Cells were plated at a density of 1 × 10⁶ cells per well in 6-well plates and allowed to grow until they reached 90% confluence. The cell layers were gently scratched using a sterile 100 µm pipette tip to create wounds. The plates were washed gently with PBS and then cultured for 48 h. The wound gaps were photographed at specified timepoints.

Statistical analysis

The statistical analysis and visualization were carried out primarily using R (version 4.2.2), and some statistical analyses were carried out using by GraphPad Prism 9 software (GraphPad Software, Inc.). Enumeration data was analyzed by Chi-square test or Fisher’s precision probability test. The Wilcox rank-sum test function in R software was employed to investigate differences between samples of HCC versus controls. Using the Chi-square test, clinicopathological parameter distributions were compared between high- and low-risk groups in LIHC. The Kaplan–Meier method was used to generate survival curves, and the log-rank test was used to determine the prognostic significance of the DRGPS risk scores. The correlation between immune scores and DRGPS prognostic characteristics was determined by Spearman’s correlation analysis. In addition, paired t tests were performed to analyze the results of the RT-qPCR. P < 0.05 indicated statistical significance.

Results

Expression profiles of DRGs in HCC patients

To begin with, we examined the TCGA–LIHC data set to determine if any DRGs were expressed differently in HCC samples compared to healthy individuals. Except for NUBPL and NDUFS1, other DRGs were differentially expressed between the two groups (Fig. 1a). As illustrated in Fig. 1b, c, a network loop diagram among these DRGs was constructed to further investigate the interaction relationship and function among them. Our next step was to summarize the CNV and somatic mutation frequencies for each DRG in HCC at the genetic level. Consequently, 19 of 371 HCC samples (5.12%) displayed gene mutations (Fig. 1d), missense mutations were the most common, and the base mutation pattern was predominantly T > A for all genes except LRPPRC. Moreover, 7 DRGs, i.e., NUBPL, NDUFS1, NCKAP1, LRPPRC, OXSM, SLC7A11, and SLC3A2, were detected with copy number amplification (Fig. 1e). In addition, Fig. 1f shows the location of these DRGs on 23 chromosome pairs. To elucidate the relationship between these 10 DRGs and HCC patients’ prognosis, Kaplan–Meier analysis showed that high expression of the DRGs was associated with poorer clinical outcomes and shorter OS (Fig. 1g). According to our findings, there is a significant difference in the expression of DRGs between HCC and normal samples, which may contribute to tumorigenesis and heterogeneity.

Fig. 1

Expression profiles of DRGs in HCC. a Differential expression of 10 DRGs in HCC and normal tissues. Blue: normal tissue; and red: tumor tissue; b network relationship loop diagram of the DRGs. The thickness of the line indicates the strength of association between DRGs. Red: positive correlation; and blue: negative correlation; c function and expression co-localization of DRGs; d mutation frequency and classification of 10 DRGs in HCC; e frequency of CNVs of the DRGs; f location of CNV alterations on chromosomes of the DRGs; g Kaplan–Meier survival analysis of 10 DRGs on the prognosis of HCC patients. DRGs, disulfidptosis-related genes; HCC, hepatocellular carcinoma; CNV, copy number variant. The P values were shown as ***P < 0.001

Identification of DRG clusters

Further exploration of the relationship between DRG expression and HCC subtype classification was conducted using unsupervised consensus clustering analysis. According to empirical CDF plots and Delta area plots, k = 2 provided the highest clustering stability for clusters A (n = 116) and B (n = 255) (Fig. S1a–c). Moreover, Kaplan–Meier curves predicted better OS for cluster B than cluster A (Fig. S1d). Significant differences in immune cell infiltration also existed between the two clusters (Fig. S1e), indicating a relationship between TME and DRG clusters. The CD4 T cells in cluster A were more likely to be activated, while the CD8 T cells were likely to be more prevalent in cluster B, as were dendritic cells, natural killer cells, T follicular helper cells, and type 2 T helper cells. Subsequent analysis of pathway mechanisms using GSVA demonstrated that cluster A was significantly enriched in “mammalian target of rapamycin (mTOR) signaling pathway”, “endocytosis”, “non-homologous end-joining”, “aminoacyl tRNA biosynthesis”, and “ubiquitin-mediated proteolysis”. In contrast, cluster B was significantly enriched in “α-linolenic acid metabolism”, “arachidonic acid metabolism”, “complement and coagulation cascade”, multiple amino acid metabolism, and other metabolism-related pathways (Fig. S1f). The heatmap in Fig. S1g suggested that DRGs were differentially expressed in different clusters. Furthermore, it demonstrated the relationship between the two DRG clusters and patients’ clinical characteristics, such as gender, age, T-stage, pathological grade, clinical stage, and survival outcome. Together, these results indicate that clustering of DRGs is of vital importance in patients with HCC.

Comprehensive analysis of HCC gene clusters based on DRG cluster-associated DEGs

To further explore the genes in DRGs clusters, we characterized 858 DRG cluster-associated DEGs using the “limma” R package and carried out secondary clustering analysis to identify gene clusters [21, 22], namely, gene clusters A and B (Fig. 2a–c). Kaplan–Meier survival curves indicated that patients in cluster B exhibited improved OS compared to those in cluster A (Fig. 2d). To examine the expression of DRGs in gene clusters A and B, it was revealed that there were significant differences in the expression of eight DRGs between the gene clusters, besides NUBPL and NDUFA11 (Fig. 2e). Analysis of DEGs combined with clinical features showed that gene cluster A was associated with advanced HCC stages and adverse clinical outcomes (Fig. 2f). To gain insight into the biological functions and pathways involved in gene clusters A and B based on DRG cluster-associated DEGs, we then proceeded with enrichment analyses using GO and KEGG in Fig. S2. GO enrichment analysis showed that GO terms were mainly engaged in α-amino acid metabolism, carboxylic acid catabolism, and cellular amino acid catabolic processes (Fig. S2a). Meanwhile, the KEGG pathway analysis demonstrated the predominant involvement of DEGs in pathways including retinol metabolism, ECM–receptor interaction, cholesterol metabolism, glycolysis, gluconeogenesis, amino acid biosynthesis, and the PPAR signaling pathway (Fig. S2b).

Fig. 2

Development and validation of the DRGPS. a K-means clustering analysis based on DRG-cluster associated DEGs identified two gene clusters; b empirical CDF plots; c delta area plots; d survival analysis; e differential expression of the 10 DRGs in two gene clusters; f correlation between two gene clusters and clinical features; g, h LASSO regression analysis; k qRT-PCR validated the expression of DRGPS risk genes in human HCC cell lines and normal liver cell lines, as well as in HCC and its adjacent tissues; i, j correlation between geneCluster, DRGcluster, and risk score. DRGs, disulfidptosis-related genes; DRGPS, DRG prognostic signature; LASSO, least absolute shrinkage and selection operator; qRT-PCR, quantitative reverse transcription polymerase chain reaction; HCC, hepatocellular carcinoma; the P values were shown as *P < 0.05; **P < 0.01; and ***P < 0.001

Development and validation of DRGPS

Kaplan–Meier analysis was applied to gene cluster-related DEGs, resulting in the identification of 96 prognosis-associated genes (P < 0.05). LASSO regression was also employed to prevent overfitting and identified 8 risk genes after univariate Cox regression (Fig. 2g, h). The DRGPS based on five risk genes (SLC7A11, MATN3, CLEC3B, CCNJL, and PON1) was derived by multivariate Cox regression analysis, followed by RT-qPCR validation of risk gene expression (Fig. 2k). Based on the coefficients and expressions of five risk genes in DRGPS, each HCC patient’s risk score was derived as follows: the risk score = (0.229) * SLC7A11 + (0.200) * MATN3 + (− 0.434) * CLEC3B + (0.391) * CCNJL + (− 0.104) * PON1. Then, patients were classified into high- and low-risk groups based on their calculated scores. In addition, risk scores were analyzed in relation to gene clusters and DRG clusters (Fig. 2i, j). A Sankey diagram shows patient survival outcomes according to DRG cluster, gene cluster, and risk group (Fig. S3a). Furthermore, we examined the expression levels of 10 DRGs, all of which were elevated in the high-risk group (Fig. S3b). After that, we divided the TCGA–LIHC data set into two groups at a ratio of 6:4, one for training and the other for testing. Risk distribution maps for the training group (Fig. S4a) and the testing group (Fig. S4b) demonstrated an ascending association between risk score and the likelihood of patient mortality, accompanied by a corresponding reduction in survival time. In the training group (Fig. S4e) and the testing group (Fig. S4f), patients with high-risk scores had a lower survival rate than those with low scores, with Kaplan–Meier survival curves showing a significant difference (P < 0.001) (Fig. S4q, r). The risk heat map showed that SLC7A11, MATN3, and CCNJL were high-risk genes, while CLEC3B and PON1 were low-risk genes (Fig. S4i, j). AUC values were further calculated evaluate the predictive capacity of ROC curves for OS in HCC patients within the training and testing cohorts (Fig. S4m, n). We enrolled all TCGA–LIHC patients again in the total test group. Subsequently, the association between DRGPS risk score and OS was validated using an external ICGC data set. As the disulfidptosis risk score escalated, there was a concurrent rise in the likelihood of patient mortality, coupled with a reduction in survival time (Fig. S4c, g, k for total TCGA, and d, h, t for ICGC). The sensitivity and specificity of DRGPS were also assessed using ROC analysis. AUC values of 0.795, 0.734, and 0.713 were calculated at 1, 3, and 5 years in the total test group, respectively (Fig. S4o). Finally, our risk score signature may be significantly associated with HCC survival and show promising predictive ability.

The clinical utility of DRGPS

Univariate and multivariate Cox regression analyses were conducted to evaluate the prognostic significance of DRGPS in patients with HCC (Figs. 3a and 3b). DRGPS was identified as an independent prognostic factor for survival (P < 0.001, hazard ratio = 1.0406, 95% confidence interval [1.018–1.063]). Furthermore, we incorporated the risk score into a ROC curve analysis alongside other clinicopathological factors (Fig. 3c) and found that the risk score achieved an AUC of 0.795, outperforming both tumor grade and stage. To improve patient prognosis prediction accuracy, a nomogram was constructed, which considered age, gender, stage, grade, and risk score as variables (Fig. 3d). The validity of this nomogram was confirmed by comparing different parameters in Fig. 3e–h with the actual survival rates. These results suggested that the DRGPS had a high degree of reliability and sensitivity. In addition, we conducted a clinicopathological stratification analysis of DRGPS. The AUC of our DRGPS was higher than other published articles [23,24,25] (Fig. S5a–c). Besides, HCC patients were separated into different groups by Kaplan–Meier analysis and that all high-risk patients exhibited significantly worse survival rates than those of low-risk patients. The DRGPS may be a reliable tool for predicting the survival of HCC patients on the basis of subclinical stratification by age, gender, TNM stage, pathological grading, and clinical stage (Fig. S6).

Fig. 3

Clinical application value of DRGPS risk score. a, b Univariate and multivariate Cox regression analyses determined the prognostic value of DRGPS risk score; c multi-ROC curves of the DRGPS risk score and other clinical factors; d a predictive nomogram predicted 1-, 3-, and 5-year OS in HCC patients; e–h calibration curve (e), clinical impact curve (f), decision curve analysis (g), concordance index (h) validated the accuracy of the nomogram. DRGPS, DRG prognostic signature; ROC, receiver operating characteristic; OS, overall survival; HCC, hepatocellular carcinoma

Comparison of TMB and stemness analysis of DRGPS

Oncological maps of tumor somatic mutations were used to identify whether mutation patterns differed between high- and low-risk groups (Fig. 4a, and Fig. 4b). The top five most commonly mutated genes in the high-risk group were TP53 (40%), TTN (24%), CTNNB1 (30%), MUC16 (20%), and ALB (10%). While the top five most commonly mutated genes in the low-risk group were TTN (23%), CTNNB1 (22%), TP53 (17%), ALB (14%), and MUC (12%). In addition, TMB was more prevalent among those at high risk (P = 0.002) (Fig. 4c). A significant difference was also observed between the high and low TMB groups as far as survival time was concerned (P = 0.021) (Fig. 4d). As can be found in Fig. 4e, when the DRGPS was combined with patients’ TMB, the prognosis of the high-risk group with high TMB was significantly worse than that of the low-risk group with low TMB. Given the critical mechanism of tumor stemness in HCC progression and drug resistance [26], we analyzed the stemness of the TCGA–LIHC cohort based on DRGPS risk scores. It was found that our risk score was positively correlated with the stemness index RNAss (R = 0.29) and EREG.EXPss (R = 0.12) (Fig. 4f, g). We further extended the risk score to pan-cancer stemness analysis. Our risk score was associated with stemness in a variety of tumors including bladder uroepithelial carcinoma (BLCA), low-grade glioma (LGG), pheochromocytoma and paraganglioma (PCPG), skin cutaneous melanoma (SKCM), stomach adenocarcinoma (STAD), and thyroid cancer (THCA).

Fig. 4

Comparison of TMB and stemness analysis based on DRGPS risk score. a, b Waterfall maps of the somatic mutations in the high-risk group and the low-risk group; c difference in TMB between high and low risk groups; d difference in survival probability between high and low TMB groups; e Kaplan–Meier analysis of the combination of TMB group and DRGPS risk score. f, g Correlation analysis of risk scores with the stemness index RNAss and EREG.EXPss; h stemness analysis of risk score in pan-cancer. TMB, tumor mutation burden; DRGPS, DRG prognostic signature. The P values were shown as *P < 0.05; **P < 0.01; and ***P < 0.001

Correlation between immune infiltration and DRGPS risk score

The complicated TME of HCC affects the efficiency of immune checkpoint inhibitors (ICIs) [27]. It was, therefore, necessary to analyze the relationship between our risk score and its immune status. The CIBERSORT analysis revealed an association between certain immune cells and our 5 risk genes (Fig. 5a). According to TIMER analysis, our risk score showed a positive correlation with M2 macrophages, while displaying a negative correlation with M1 macrophages and naive B cells (Fig. 5b–d). TME scores were determined using ESTIMATE, which revealed that our risk score was associated with both stromal scores and ESTIMATE scores (Fig. 5e). After creating a heat map of immune cell infiltration via seven algorithms, we visually depicted the correlation between our risk score and the TME (Fig. 5f). A comparison of the expression of common ICGs between the two risk groups was conducted to investigate the capacity of our risk score to predict treatment options and the efficacy of ICIs. A notable disparity in the expression of ICGs including CD274 (PDL1), IDO1, and CD44 was observed between the high- and low-risk groups (P < 0.001) (Fig. 5g), lending support to the use of ICIs in patients with HCC at different risk. Furthermore, ssGSEA revealed distinct immune function disparities between the two risk groups in four categories that included cytolytic activity, major histocompatibility complex (MHC) molecules of class I, and type I and type II interferon responses (Fig. 5h). Last but not least, immune cells in pan-cancers were correlated with DRGPS risk scores (Fig. 5i, j). As a result of these findings, DRGPS is believed to be reliable for predicting HCC patients' immune status.

Fig. 5

Correlation between DRGPS risk score and immune infiltration. a CIBERSORT algorithm analyzed the correlation of immune cells with the risk genes [2]. TIMER algorithm analyzed the correlation between risk scores and macrophages M1, M2, and B cells naïve immune cell subpopulations; (e) ESTIMATE algorithm analyzed the correlation between risk scores and stromal score, immune score, and ESTIMATE score; (f) immune landscape heatmap between high- and low-risk groups based on DRGPS; (g) differences in expression levels of different immune checkpoint genes between high- and low-risk groups; (h) ssGSEA analysis of differences in immune function between the two risk groups; Correlation of high risk scores (i) and low risk scores (j) with immune cells in pan-cancer. DRGPS, DRG prognostic signature; ssGSEA, single sample gene set enrichment analysis. The P values were shown as *P < 0.05; **P < 0.01; and ***P < 0.001

Chemosensitivity determined in HCC patients using DRGPS risk score

The development of clinically effective chemotherapeutic agents for HCC patients has been slow due to the lack of guidance from validated prognostic biomarkers [28]. Disulfidptosis is different from the traditional PCD mode and represents a new therapeutic opportunity in the era of apoptosis resistance. Therefore, our objective was to determine the chemosensitivity of each risk group based on GDSC results. In comparison with the low-risk group, the high-risk group was more sensitive to Alisertib, Bortezomib, Crizotinib, Daporinad, Docetaxel, Ruxolitinib and Sorafenib. However, Selumetinib exhibited a low sensitivity in the high-risk group, while the low-risk group would probably benefit from therapy with these drugs (Fig. S7). Collectively, these results suggest that DRGPS may be used to guide chemosensitivity in patients with HCC.

SLC7A11 induced HCC cell proliferation and migration in vitro

SLC7A11 expression was analyzed in an HCC tissue microarray consisting of 80 pairs of HCC and normal paracancerous tissues. Significant upregulation of SLC7A11 was observed in human HCC tissues compared with paired paracancerous tissues (P = 0.0073) (Fig. 6a). We further analyzed the relationship between the expression level of SLC7A11 and the clinicopathological factors of HCC patients. We defined SLC7A11 scores greater than 3 as high expression and vice versa as low expression. The results showed that high expression of SLC7A11 was significantly correlated with more advanced TNM stage and larger tumor size (Table 1, both P < 0.05).

Fig. 6

SLC7A11 promotes the proliferation and migration of HCC cells in vitro. a Immunohistochemistry of SLC7A11 in HCC and paracancerous tissues. Hep3B and HepG2 cell were selected and transfected with siSLC7A11. Evaluation of migration and proliferation capacity by wound healing assay (b, c), transwell assay (d, e) and CCK-8 assay (f, g). HCC, hepatocellular carcinoma; CCK-8, Cell Counting Kit 8

Table 1 Relationship between SCL7A11 expression and clinicopathological factors of HCC patients

To further investigate the biological role of SLC7A11, we initially generated HCC cell lines with reduced SLC7A11 expression and subsequently conducted a series of experiments to examine whether SLC7A11 could regulate tumor cell proliferation and migration. Consequently, the closure of the wound area was significantly delayed in siSLC7A11 cells compared to the control cells (Fig. 6b, c). As anticipated, this outcome was further confirmed by the Transwell assay (Fig. 6d, e). Furthermore, the CCK-8 assay revealed that reduced expression of SLC7A11 suppressed the proliferation activity of HepG2 and Hep3B cells (Fig. 6f, g). Taken together, these results indicate that SLC7A11 may play a role in promoting the proliferation, migration, and invasion of HCC cells.

Discussion

Due to the heterogeneous nature of HCC, poor prognosis have been reported, necessitating a more detailed understanding of its genetic characteristics and prognosis [29]. Disulfidptosis is a novel form of fast cellular death resulting from the collapse of the actin cytoskeleton under disulfide stress, which holds significant promise for early diagnosis, treatment selection, and prognosis prediction [7].

This study initiated with the analysis of the differential expression of 10 DRGs in the TCGA–LIHC cohort, of which 8 DRGs were prominently associated with OS in HCC patients. Through consensus clustering analysis, we categorized HCC patients into two subgroups in accordance with these 10 DRGs. As shown in Kaplan–Meier analysis, the prognosis of DRGs-A cluster was significantly less favorable than that of DRGs-B cluster, with significant differences in clinical indicators and immune infiltration. Specifically, the proportion of activated B cells and activated CD8 T cells in infiltration was higher in the DRGs-B cluster. As documented previously, the cross-talk between B cells and T cells could enhance local immune activation in HCC, and activated CD8 T cells were regarded as the mainstay of anti-tumor immunity, which would contribute to a better prognosis for HCC patients [30].

To develop a DRGPS, we identified 5 important DRGs (SLC7A11, MATN3, CLEC3B, CCNJL, and PON1) based on secondary clustering of gene cluster-associated DEGs sequentially with Kaplan–Meier analysis, LASSO, as well as univariate and multivariate Cox regression analysis. These 5 genes have been identified as tightly implicated in cancer deterioration. SLC7A11, as a cystine transporter protein, protects tumor cells against oxidative stress by mediating the synthesis of glutathione and has a crucial function in ferroptosis and disulfidptosis [8]. In our in vitro experiments, SLC7A11 promoted HCC cell proliferation and migration. Moreover, the IHC results of our HCC microarray indicated that high expression of SLC7A11 correlated with advanced TNM clinical stage and larger tumor size, which contributed to the maintenance of stemness in CSCs, TME and chemoresistance [31, 32]. PON1 is a serum enzyme related to high-density lipoprotein. It has been determined that AFP-negative HCC is associated with N-glycosylation of PON1, which could be utilized as a potential biomarker [33]. There is evidence that downregulation of CLEC3B in HCC patients has a negative impact on prognosis, and CLEC3B may also serve as a potential biomarker of HCC [34]. CLEC3B exerts its influence in tumor aggression and metastasis; moreover, exosomal CLEC3B attenuates angiogenesis and inhibits tumor cell migration, reshaping the TME [35]. Besides, Hu B et al. employed CCNJL into a prognostic model to evaluate the OS of HCC patients [36]. In addition, MATN3 is a structural domain member of von Willebrand factor A, which participates in the filamentous network formation of the extracellular matrix [37]. Dai et al. utilized MATN3 as one of the EMT-related genes to develop a prognostic model of gastric cancer patients [38]. All these 5 genes are associated with survival, hence a disulfidptosis risk signature containing these 5 genes would be a promising prognostic tool for HCC.

Early intervention in high-risk populations for HCC can effectively reduce mortality and prolong survival time [39]. Therefore, we further categorized HCC patients into two DRGPS risk score groups, i.e., high- and low-risk, which showed significant differences in prognosis. We also examined the performance of our DRGPS risk score in the TCGA–LIHC and ICGC data sets. With the incorporation of the clinical characteristics, the overall AUC of our risk score was 0.795, even better than TNM stage. To further investigate the influence of disulfidptosis on HCC heterogeneity and consequent clinical outcomes, the risk score was accurately quantified using a nomogram model [40, 41]. Corresponding 1-, 3-, 5-year AUC of our DRGPS were higher than other published articles [23,24,25]. Therefore, our DRGPS could serve as an independent prognostic biomarker and exhibited superior diagnostic capacity for HCC patients.

Tumor evolution is motivated by intrinsic factors, namely, genomic instability caused by DNA deleterious mutations in tumor cells, and extrinsic factors, such as immune editing induced by immune cells within the TME [42]. As a highly heterogeneous tumor, HCC develops through the acquisition of multiple genetic and molecular variations during the whole oncogenic process [43]. High TMB can induce tumor chemoresistance, T-cell impairment, and chromosomal destabilization, often suggesting a worse prognosis [44]. However, in patients receiving ICIs, high TMB is associated with potent anti-tumor immunity and a better immunotherapeutic response. Also, high TMB may become a new predictive biomarker for patients receiving ICIs (such as PD1 or CTLA4 inhibitors) [45]. To examine the potential relevance between risk score and TMB, survival analysis in this study indicated that HCC patients with high risk and high TMB would have less survival probability and a worse prognosis. Mutations in TP53 and CTNNB1 were more prominent in the high-risk group compared to the low-risk group. Mutations in TP53 and CTNNB1 are also the ubiquitous variants of HCC, and exert pivotal effects in the immunosuppressive microenvironment [46]. Given that CSCs are crucial in tumorigenesis, treatment resistance, and tumor recurrence in HCC [26], we then conducted a correlation analysis between risk score and HCC stemness. Our risk score displayed a positive correlation with tumor stemness index (TSI), and a comprehensive pan-cancer analysis revealed a robust and statistically significant association between our risk score and TSI across a spectrum of malignancies. This observation substantiates that tumor cells in the high-risk group exhibited augmented aggressiveness and heightened resistance to therapeutic interventions, in contrast to their low-risk counterparts. These findings collectively underscore the predictive capacity of our DRGPS in delineating both the TMB and stem cell characteristics specific to HCC patients with a high degree of precision.

The TME has an essential function in the genesis and advancement of HCC, showing a complicated relationship with tumor cells, stromal cells, and immune cells [47]. Approximately 50% of HCC patients have been reported to receive systemic therapy, and ICIs have revolutionized the administration [48]. Our analysis of immune infiltration revealed a negative correlation between our risk score and M1 macrophages, as well as a positive correlation with M2 macrophages. Immune function enrichment unveiled elevated MHC class I activity in the high-risk group, whereas low-risk group exhibited a propensity for heightened interferon signaling and cytolytic activity. Previous scholarly inquiries have delineated that M2 macrophages exert impact on tumor invasion, advancement, and immune evasion mechanisms by upregulating MHC class I expression and releasing anti-inflammatory cytokines [49]. M2 macrophages could collaborate with HCC stem cells to facilitate resistance to Sorafenib [50]. Instead, the shift of M2 to M1 macrophage polarization in HCC could activate an anti-tumor immune response [51]. A recent study classified HCC patients via immunophenotype [46]. Patients with higher interferon signaling and cytolytic activity and a higher proportion of M1 macrophages were identified as the inflammatory group of HCC, consistent with our low-risk group. Conversely, the heightened frequency of TP53 and CTNNB1 gene mutations, coupled with the absence of interferon signaling and a greater prevalence of immunosuppressive M2 macrophages, collectively suggested a strong association of our high-risk group with the non-inflammatory subtype of HCC. These findings support that our risk score signature demonstrates exceptional predictive efficacy in characterizing the immune classification of HCC.

HCC therapy is marked by intricate challenges and heightened chemoresistance. Presently, the comprehension of the heterogeneous mechanisms underlying HCC is constrained, and there is a dearth of validated prognostic biomarkers to inform clinical decisions. The advancement of clinically validated pharmaceuticals for HCC is progressing at a sluggish pace [52]. Disulfidptosis stands apart from traditional PCD patterns and represents a novel therapeutic opportunity in the landscape of chemoresistance. Therefore, this study utilized three drug databases to evaluate the chemosensitivity between the high- and low-risk groups. In addition to Sorafenib, the standard treatment for advanced HCC, drug sensitivity analysis revealed that the high-risk group exhibited increased susceptibility to compounds, such as Alisertib, Bortezomib, Crizotinib, Daporinad, Docetaxel and Ruxolitinib. Alisertib, an Aurora A kinase inhibitor, induces cell cycle arrest in the G2/M phase and promotes autophagy via the PI3K/Akt/mTOR pathway, showing significant antiproliferative effects on HCC cell lines Hep3B and HepG2 [53, 54]. Bortezomib, an FDA-approved proteasome inhibitor, can block the cell cycle and inhibit stemness-induced immune responses. However, the efficacy of Bortezomib, either as monotherapy or in combination with Sorafenib, remains limited and is dependent on more precise patient stratification and the identification of sensitivity biomarkers [55]. Patients in the high-risk group identified by this study may benefit from these therapeutic agents. Crizotinib has also shown promising effects in the treatment of HCC. It inhibits HCC cell proliferation by suppressing the phosphorylation of ALK, AKT, and ERK while promoting apoptosis [56]. When combined with doxorubicin, Crizotinib synergistically enhances HCC cell death by reducing multidrug resistance protein 1 and activating the autophagic cell death pathway [57]. Daporinad, a Notch signaling inhibitor, significantly suppresses HCC proliferation and migration by regulating the extracellular matrix and inhibiting the Hes1/PTEN/AKT/mTOR signaling pathway [58]. Docetaxel induces G2/M phase arrest and promotes reactive oxygen species generation and glutathione depletion, leading to apoptosis in HCC cells [59]. However, a phase II clinical trial demonstrated limited efficacy and significant toxicity of Docetaxel in patients with advanced liver cancer [60]. The heterogeneity of HCC restricts the effectiveness of Docetaxel, necessitating treatment selection based on the molecular and genetic characteristics of the patient's tumor [61]. Moreover, in vitro study has shown that Ruxolitinib inhibits JAK/STAT signaling, reducing proliferation and colony formation in HCC cells [62]. It also targets JAK1/2 to regulate hepatic stellate cell function, inhibiting liver fibrosis progression and accelerating its reversal [63]. Experiments in vivo have confirmed that Ruxolitinib can slow disease progression and prolong survival in mice with metastatic HCC [64]. The low-risk group revealed heightened sensitivity to Selumetinib. Selumetinib, a MEK1/2 inhibitor, has been found with limited single-agent activity in advanced HCC, despite evidence of target inhibition [65]. However, in a Phase Ib study for advanced HCC, the combination of Selumetinib with Sorafenib showed an acceptable safety profile and encouraging anti-tumor activity [66]. Hence, we conducted a chemosensitivity analysis based on our risk scores of HCC patients, facilitating the optimization of chemotherapeutic selection between the high- and low-risk groups, thereby augmenting the efficacy of chemotherapy and guiding the individualized treatment of HCC.

Conclusion

This study systematically depicts the DRGPS developed for HCC and investigates its potential clinical implications. It demonstrates independent prognostic value, along with strong sensitivity and reliability in predicting HCC patients’ survival. This DRGPS sheds light on the potential pathogenesis of TMB, TSI, and immune infiltration, enabling the identification of HCC patients who may benefit from enhanced responses to immunotherapy and chemotherapy. Meanwhile, SLC7A11 may played a critical role in promoting the development of HCC. Therefore, it is of paramount significance for gaining an in-depth comprehension of the intricate mechanisms inherent in our risk signature, thereby fostering individualized therapeutic strategies for HCC patients.

Limitation

This study still has several limitations. First, this was a bioinformatic analysis based primarily on data sourced from TCGA and ICGC databases, with the lack of larger sample size. Therefore, further validation of the constructed model with prospective real-world data is necessary to assess its clinical utility. Second, there is a paucity of relevant studies on disulfidptosis, and the DRGs used in our analysis may not fully represent the entire profile of HCC-specific disulfidptosis. Finally, this study conducted a preliminary exploration into the expression and functionality of SLC7A11 in HCC. Further efforts are warranted to elucidate comprehensively the molecular mechanisms underlying disulfidptosis in HCC.