Comprehensive analysis of disulfidoptosis-related genes reveals molecular heterogeneity and key regulators in retinoblastoma progression

Introduction

Retinoblastoma is the most common intraocular malignancy in children, with a median age of diagnosis at approximately 2 years. It is the leading ocular disease causing blindness and disability in childhood (1). There was a global prevalence of 7.46 for retinoblastoma in 2021, and the global prevalence rate of retinoblastoma increased from 1990 to 2021, with average annual percentage changes of 0.68 (2). The primary treatment goal in retinoblastoma is to save the child’s life through early detection, treatment of the ocular tumour, and prevention of metastatic spread. With advances in medical technology, treatment strategies have shifted towards maximising visual potential and improving the quality of life while ensuring survival and preserving the eyeball (3,4). Survival and eye preservation rates vary significantly across regions due to disparities in economic conditions and healthcare infrastructure (5). There is a 3-year survival rate of 99.5% in high-income countries, while the 3-year survival rate is 80.3% and 57.3% in lower-middle-income countries and low-income countries, respectively (5), and such differences may be explained by the delay in diagnosis or lag time (6). Nevertheless, metastases still occur in some cases, either via scleral invasion into the orbit or through optic nerve extension into the central nervous system (7), which is life-threatening. Therefore, early identification of progressive retinoblastoma and timely initiation of treatment are critical.

Retinoblastoma typically results from biallelic mutations in the RB1 gene, located on chromosome 13q14.2 (8), and can be classified as inherited (approximately 45%) or non-inherited (approximately 55%) (9). Hereditary cases are usually bilateral, while non-hereditary cases tend to be unilateral. Several studies based on gene expression profiling have revealed genetic and molecular heterogeneity in retinoblastoma beyond RB1 mutations (10,11), suggesting the existence of distinct molecular subtypes—an insight that is critically important for determining optimal treatment strategies. Despite substantial progress in elucidating the mechanisms underlying retinoblastoma progression, the development of retinoblastoma-specific targeted therapies remains limited (12).

Cells can undergo diverse forms of programmed cell death (PCD), which are closely linked to regulation of the tumour immune microenvironment, metastasis, and progression (13,14). The induction of PCD—including necroptosis, autophagy, and ferroptosis—has emerged as a crucial strategy in cancer therapy (13). Disulfidoptosis, a newly recognised form of cell death (15), has been investigated in various cancers (16-18), but its role in retinoblastoma remains unexplored.

In this study, we aimed to examine the role of disulfidoptosis-related genes (DRGs) in retinoblastoma, and their potential as therapeutic targets. Using data from the Gene Expression Omnibus (GEO) database, we systematically analysed the dysregulation of DRGs in retinoblastoma. Potential key DRGs were identified through machine learning algorithms. Additionally, we characterized distinct molecular subtypes based on differentially expressed DRGs and investigated their specific signalling pathways and immune profiles. This study offers novel insights into the molecular heterogeneity of retinoblastoma and highlights potential targets for its diagnosis and treatment. We present this article in accordance with the TRIPOD reporting checklist (available at https://tp.amegroups.com/article/view/10.21037/tp-2024-596/rc).

Methods

The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Data acquisition and preprocessing

Three microarray datasets relevant to retinoblastoma—GSE208143, GSE97508, and GSE24673—were obtained from the GEO (https://www.ncbi.nlm.nih.gov/geo/) database. Detailed information regarding sample sizes and platforms for these three datasets is presented in Table 1. The data from GSE208143 and GSE97508 were merged into a larger dataset after batch effect correction using the combat function from the sva package (version 3.50.0), which was then used as the training set in this analysis. The GSE24673 dataset was retained for validation purposes.

Screening of differentially expressed DRGs (DE-DRGs)

Differentially expressed genes (DEGs) between retinoblastoma and control samples were identified using the limma package (version 3.58.1), followed by adjustment for multiple testing using the Benjamini-Hochberg method. A threshold of adjusted P<0.05 and |log₂fold change (FC)| >1 was applied to screen for DEGs. Sixteen DRGs were retrieved from a previous study (19). The intersection of the DEGs and DRGs yielded the final list of DE-DRGs.

Consensus clustering

To explore potential molecular subtypes of retinoblastoma, unsupervised hierarchical clustering was performed using the ConsensusClusterPlus package (version 1.66.0) based on the expression profiles of the DE-DRGs. Consensus Clustering is a resampling-based approach that involves subsampling 80% of the samples at each iteration, and partitioning each subsample into up to k (maximum, k=5) groups using the k-means algorithm. This procedure was repeated 1,000 times. The resulting consensus score matrix and cumulative distribution function (CDF) curve were visualised and used to determine the optimal number of clusters.

Machine learning algorithm

Key genes were further identified from the DE-DRGs using two machine learning algorithms: support vector machine-recursive feature elimination (SVM-RFE) and least absolute shrinkage and selection operator (LASSO) regression analysis. SVM-RFE (20) is a feature selection method based on SVMs. It evaluates the importance of each feature by training an SVM classifier and retaining the most important features. This process is then recursively applied to the remaining feature sets. RFE is a backward feature elimination technique that begins with all features and progressively removes those with the smallest weights. SVM-RFE is particularly effective for feature selection in linearly inseparable data, and has been widely applied in this context (21). In this study, SVM-RFE was performed using the e1071 package (version 1.7-14) to calculate the importance and rank of each DE-DRG, and the optimal gene signature was selected based on the highest accuracy and lowest error rate. LASSO is a generalized linear model that incorporates an L1-norm penalty as a regularization term to the ordinary least squares method. This approach reduces model complexity and overfitting by shrinking coefficients of insignificant variables to zero (22,23). In this analysis, LASSO regression was performed with 10-fold cross-validation using the glmnet package (version 4.1-8) and repeated for 1,000 iterations to select feature variables based on the minimum criteria. Genes selected by both algorithms were considered key genes.

Enrichment analyses

The biological functions associated with the DE-DRGs were examined using enrichment analysis via the clusterProfiler package (version 4.4.4). Gene set enrichment analysis (GSEA) was conducted to investigate the mechanisms associated with the key genes. Briefly, genes were ranked based on fold changes in expression, and GSEA was performed using adjusted P<0.05 and | normalized enrichment score (NES) | >1 as threshold criteria.

Gene set variation analysis (GSVA)

GSVA was performed to calculate the disulfidoptosis score and to explore functional pathway differences between the two molecular subtypes. The predefined gene sets from the “KEGG” and “Hallmark” collections in the MSigDB database were used as the enrichment background. To calculate the disulfidoptosis score, a gene set defined by the DRGs was employed. Enrichment scores for each pathway were computed using the GSVA package (version 1.50.1), and comparisons were made between the two subtypes or between retinoblastoma and control samples.

Evaluation of immune infiltration

The infiltration levels of 28 types of immune cells were quantified using the single-sample GSEA (ssGSEA) algorithm. In addition, stromal and immune scores for each retinoblastoma tissue sample were estimated based on gene expression profiles using the ESTIMATE algorithm. The sum of these two scores served as an indirect indicator of tumour purity. Differences in immune cell infiltration levels and each score between the two subtypes or between retinoblastoma and controls were assessed using the Wilcoxon test. Spearman correlation analysis was conducted to evaluate associations between key genes and immune cell infiltration levels.

Cell culture and transfection

The human retinoblastoma cell lines Weri-Rb1 (#CL-0465, Wuhan Pricella Biotechnology Co., Ltd., Wuhan, China), Y79 (Wuhan SUNNCELL Biotechnology Co., Ltd.), and HXO-Rb44 (Central Laboratory of Hunan Medical University, Changsha, China) were cultured in RPMI-1640 medium (#11875119, Gibco, California, USA). The normal retinal epithelial cell line ARPE-19 (#CL-0026, Pricella) was maintained in DMEM/F12 medium (#PM150312, Pricella) at 37 ℃ in a humidified atmosphere with 5% CO₂. All media were supplemented with 10% foetal bovine serum, penicillin-streptomycin.

Lentiviral transduction was employed to overexpress EPAS1 in HXO-Rb44 cells to explore its functions. The coding sequence of EPAS1 was obtained from the NCBI database. HXO-Rb44 cells were seeded into 6-well plates and infected with lentivirus (1×10⁸ TU/mL) carrying either oe-EPAS1 or a negative control sequence when cell confluence reached 70–90%. Stably transfected cells were selected 72 h after infection.

RNA extraction and quantitative reverse transcription polymerase chain reaction (qRT-PCR)

Total RNA was isolated from cells using TRIzol reagent and reverse-transcribed into complementary DNA. PCR amplification was then carried out using the SYBR Green PCR Master Mix under reaction conditions: 95 ℃ for 10 min and 40 cycles of 95 ℃ for 12 s and 60 ℃ for 40 s. The sequences of the used primers were as follows: EPAS1-F: 5'-TGACAGCTGACAAGGAGAAGAA-3', EPAS1-R: 5'-TGTGTTCGCAGGAAGCTGAT-3', SLC7A11-F: 5'-TGGTCAGAAAGCCTGTTGTGT-3', SLC7A11-R: 5'-CATGGAGCCAAAGCAGGAGA-3', GAPDH-F: 5'-GAAGGTCGGAGTCAACGGAT-3', GAPDH-R: 5'-CTTCCCGTTCTCAGCCATGT-3'.

Western blotting

Total protein was isolated following cell lysis and quantified using the bicinchoninic acid assay. The proteins were electrophoretically separated, transferred onto polyvinylidene fluoride (PVDF) membranes, and blocked. Afterwards, primary antibodies: EPAS1 (#DF2928, Affinity Biosciences, Jiangsu, China) and SLC7A11 (#DF12509, Affinity Biosciences) and GAPDH (ab181602, Abcam, Cambridge, UK), co-incubation with the membranes were carried out overnight at 4 ℃, followed by incubation with secondary antibodies for 60 min at dark. Subsequently, the blotting bands were visualised using enhanced chemiluminescence (ECL) reagent, followed by film exposure.

Cell viability and apoptosis assays

Cell viability was determined using a cell counting kit-8 (CCK-8) assay. The cells were seeded in 96-well plate for 24 h of culture. Subsequently, 10 µL of CCK-8 solution was added for each well for another 2 h of culture, and absorbance at 450 nm was measured to evaluate the cells viability. An annexin V/PI assay was used to detect apoptosis. Briefly, the cells were collected and resuspended in binding buffer. Apoptosis was determined by flow cytometry after staining with Annexin V-fluorescein isothiocyanate (FITC) and propidium iodide (PI).

Cell migration and invasion assays

A wound-healing assay was performed to assess cell migration. Cells were seeded in 6-well plate (5×10⁵ cells/well) overnight and wounded using a pipette tip. Detached cells were washed by PBS, and serum-free medium was added, followed by incubation in an incubator (37 ℃ and 5% CO₂). The plates were photographed at 0 and 24 h, respectively. A transwell assay was used to assess cell invasion. Briefly, cells were inoculated into the upper chamber of a Matrigel gel-pre-coated Transwell, and RPMI-1640 complete medium containing 10% fetal bovine serum (FBS) was added to the lower chamber. Following 24 h of incubation, cells that migrated to the lower chamber were fixed with absolute methanol and stained with crystal violet. After washing with PBS, the cells were photographed and counted under a microscope using the ImageJ software.

Nicotinamide adenine dinucleotide phosphate (NADP)+/reduced NADP (NADPH) measurement

The NADP+/NADPH ratio in the cells was determined using an NADP+/NADPH assay kit (#S0179, Beyotime) with a colour reaction based on WST-8, according to the manufacturer’s instructions.

Statistical analysis

All experimental data were presented as mean ± standard deviation and were statistically analysed and graphed using GranphPad 8.0 software. Differences among groups were detected using one-way analysis of variance and Tukey’s multiple comparison test. Statistical significance was set at P<0.05.

Results

Dysregulation of DRGs in retinoblastoma

After batch effect correction, the GSE208143 and GSE97508 datasets were merged into a larger dataset. Principal component analysis (PCA) confirmed a uniform distribution of samples, indicating that batch effects had been substantially mitigated (Figure 1A). This combined dataset was used for subsequent analyses.

Figure 1 Dysregulation of DRGs in retinoblastoma. (A) PCA illustrating the distribution of samples in the GSE208143 and GSE97508 datasets before and after batch effect correction; (B) disulfidoptosis scores of retinoblastoma and control samples. ***P<0.001. Heatmap (C) and volcano plot (D) showing the expression patterns of DEGs in retinoblastoma and control samples. (E) Venn diagram identifying the shared genes between DEGs and DRGs; (F-G) Enriched Gene Ontology terms (F) and KEGG pathways (G) for the DE-DRGs. ATP, adenosine triphosphate; BP, biological process; CC, cellular component; DE-DRG, differentially expressed-DRG; DEG, differentially expressed genes; DRG, disulfidoptosis-related genes; FC, fold change; GEO, Gene Expression Omnibus; GSVA, gene set variation analysis; KEGG, Kyoto encyclopedia of genes and genomes; MF, molecular function; NADH, nicotinamide adenine dinucleotide phosphate; NADPH, reduced nicotinamide adenine dinucleotide phosphate; PCA, principal component analysis; PC, principal component; RB, retinoblastoma.

Sixteen DRGs were retrieved from a previous study (19). The score was significantly different between retinoblastoma and control samples (Figure 1B, P<0.001), suggesting the potential involvement of disulfidoptosis in retinoblastoma. We subsequently examined the expression patterns of DRGs. Initially, 4,487 DEGs were identified, including 2,013 up-regulated and 2,474 down-regulated genes in retinoblastoma compared to controls (Figure 1C,1D). Intersection of the 16 DRGs with the DEGs yielded 8 DE-DRGs (Figure 1E). Enrichment analysis revealed that these DE-DRGs were primarily associated with fatty acid biosynthesis and multiple oxidative phosphorylation-related functions, including mitochondrial respiratory chain complex assembly, reduced nicotinamide adenine dinucleotide (NADH) dehydrogenase complex assembly, mitochondrial electron transport from NADH to ubiquinone, and mitochondrial adenosine triphosphate (ATP) synthesis coupled electron transport (Figure 1F,1G).

Disulfidoptosis-based molecular subtypes of retinoblastoma

Based on the expression profiles of the 8 DE-DRGs, 33 retinoblastoma samples were classified into two clusters (C1 and C2) using consensus clustering (Figure 2A,2B). PCA demonstrated clear separation between the two clusters, confirming the successful identification of disulfidoptosis-based molecular subtypes (Figure 2C). The C2 subtype exhibited activation of several tumourigenesis-associated hallmark pathways, including MYC targets, glycolysis, the p53 pathway, and mTORC1 signalling, while the C1 subtype showed activation of pathways related to inflammatory response, myogenesis, and coagulation (Figure 2D). In terms of KEGG pathways, 109 pathways were significantly different between the two clusters (Figure 2E). The C1 subtype showed activation of 49 pathways, such as neuroactive ligand-receptor interaction, calcium signalling, and cytokine-cytokine receptor interaction. In contrast, the C2 subtype showed activation of 60 pathways, including mismatch repair, DNA replication, and cell cycle.

Figure 2 Disulfidoptosis-based molecular subtypes of retinoblastoma. (A) Consensus CDF curve used to determine optimal k values; (B) consensus matrix heatmap for k=2; (C) PCA plot illustrating the separation of samples into two clusters; (D) significantly enriched hallmark pathways differing between the two clusters, identified by GSVA; (E) significantly enriched KEGG pathways differing between the two clusters, identified by GSVA. CDF, cumulative distribution function; GSVA, gene set variation analysis; KEGG, Kyoto encyclopedia of genes and genomes; PC, principal component; PCA, principal component analysis.

Immune infiltration characteristics of the two subtypes were subsequently evaluated. Among the 28 immune cell types analysed, 20 showed significant differences in infiltration levels between the C1 and C2 subtypes (P<0.05), with most immune cells exhibiting markedly lower infiltration in the C2 subtype (Figure 3A). This suggests a potentially immunosuppressive microenvironment in the C2 subtype. Consistently, ESTIMATE analysis revealed that the C2 subtype exhibited higher tumour purity, whereas the C1 subtype displayed higher immune and stromal scores (Figure 3B, P<0.001). Moreover, the expression levels of multiple immune checkpoint genes were significantly lower in the C2 subtype compared with the C1 subtype (Figure 3C, P<0.05).

Figure 3 Immune characteristics of the two molecular subtypes. (A) Heatmap showing the infiltration levels of 28 immune cell types across the two clusters; (B) box plots showing differences in ESTIMATE score, immune score, stromal score, and tumour purity between the two clusters; (C) box plots showing differences in the expression of immune checkpoint genes between the two clusters. *, P<0.05; **, P<0.01; ***, P<0.001; ns, non-significance.

Machine learning algorithms identify key genes in retinoblastoma

From the 8 DE-DRGs, key genes were further identified using machine learning algorithms. The SVM-RFE algorithm selected three feature genes—SLC7A11, EPAS1, and OXSM—which yielded the highest classification accuracy (0.997) and the lowest error (0.00314) during 5-fold cross-validation (Figure 4A). In LASSO regression, the minimum residual sum of squares was observed at a λ (lambda.min) value of 3e−4, identifying EPAS1 and SLC7A11 as feature genes (Figure 4B). Two overlapping genes—EPAS1 and SLC7A11—were consequently designated as key genes in retinoblastoma (Figure 4C). In both the training and validation sets, the expression levels of EPAS1 and SLC7A11 were significantly down-regulated in retinoblastoma samples compared with controls (Figure 4D, P<0.001 in the training set and P<0.05 in the validation set).

Figure 4 Machine learning algorithms identify key genes in retinoblastoma. (A) Accuracy (left) and error rate (right) of the SVM-RFE model for identifying feature genes; (B) coefficient distribution and binomial deviance profile of the LASSO regression model; (C) venn diagram showing overlapping genes identified by both SVM-RFE and LASSO; (D) box plots showing the expression levels of EPAS1 and SLC7A11 in retinoblastoma and control samples. *, P<0.05; ****, P<0.0001. CV, cross validation; LASSO, least absolute shrinkage and selection operator; RB, retinoblastoma; SVM-RFE, support vector machine-recursive feature elimination.

Functional roles of the key genes in retinoblastoma

GSEA was performed to investigate the biological functions associated with EPAS1 and SLC7A11 in retinoblastoma. The top 5 pathways were displayed in Figure 5A. Expression of both EPAS1 and SLC7A11 was positively correlated with oxidative phosphorylation (NES =3.44 for EPAS1; NES =2.72 for SLC7A11), while negatively correlated with olfactory transduction (NES =−3.61 for EPAS1; NES =−3.00 for SLC7A11). Besides, expression of SLC7A11 was also linked to cell cycle (NES =−3.03) and DNA replication (NES =−2.80) (Figure 5A). GeneMANIA analysis identified the top 20 genes interacting with EPAS1 and SLC7A11, which were predominantly involved in cellular responses to hypoxia or oxygen levels, protein hydroxylation, and amino acid transmembrane transporter activity (Figure 5B).

Figure 5 Functional analysis of EPAS1 and SLC7A11. (A) Top five KEGG pathways associated with the expression of EPAS1 and SLC7A11; (B) PPI network for EPAS1 and SLC7A11, analysed using GeneMANIA, showing the top 20 interacting genes and their associated functions. KEGG, Kyoto encyclopedia of genes and genomes; NES, normalized enrichment score; PPI, protein-protein interaction.

Correlations of key genes with immune cells infiltration

Immune cell infiltration in retinoblastoma was analysed, and significant differences were observed in 19 immune cell types between retinoblastoma and control samples (Figure 6A, P<0.05). Notably, several immune cell types exhibited higher infiltration levels in retinoblastoma samples (Figure 6A). The expression levels of both EPAS1 and SLC7A11 were significantly correlated with multiple immune cell populations (Figure 6B). For example, EPAS1 expression was positively correlated with plasmacytoid and immature dendritic cells, and negatively correlated with activated CD4⁺ T cells and B cells. Similarly, SLC7A11 expression was positively correlated with CD56^bright natural killer cells and negatively correlated with memory B cells.

Figure 6 Correlation of immune cell infiltration with EPAS1 and SLC7A11 expression. (A) Infiltration levels of 28 immune cell types in retinoblastoma and control samples, analysed using ssGSEA; (B) correlations between immune cell infiltration levels and the expression of EPAS1 and SLC7A11. *, P<0.05; **, P<0.01; ***, P<0.001; ns, not significant. ssGSEA, single-sample gene set enrichment analysis.

In vitro validation of EPAS1 and SLC7A11 expression

EPAS1 and SLC7A11 expression were validated in three retinoblastoma cell lines: Weri-Rb1, Y79, and HXO-Rb44. Compared to the normal retinal epithelial cell line ARPE-19, both the mRNA and protein levels of EPAS1 and SLC7A11 were reduced in all three retinoblastoma cell lines, particularly in HXO-Rb44 cells (Figure 7A,7B). This is consistent with the findings of bioinformatics analysis.

Figure 7 Expression of EPAS1 and SLC7A11 in retinoblastoma cells. (A) mRNA expression of EPAS1 and SLC7A11 in retinoblastoma cell lines (Weri-Rb1, Y79, and HXO-Rb44) and the normal retinal epithelial cell line ARPE-19, determined by qRT-PCR; (B) protein expression of EPAS1 and SLC7A11 in the same cell lines, determined by Western blot analysis. *, P<0.05; **, P<0.01. qRT-PCR, quantitative reverse transcription polymerase chain reaction.

EPAS1 overexpression inhibited the malignant cytological behaviors of retinoblastoma cells

Given the previously reported role of SLC7A11 in retinoblastoma in a previous study (24), we investigated the role of EPAS1 in retinoblastoma. EPAS1 was successfully overexpressed in HXO-Rb44 cells (Figure 8A,8B). EPAS1 overexpression markedly suppressed the viability (Figure 8C), migration (Figure 8D), and invasion (Figure 8E) of HXO-Rb44 cells, while enhanced apoptosis of HXO-Rb44 cells (Figure 8F). Emerging evidence has demonstrated a close link between disulfidoptosis and the depletion of NADPH, which serves as an important reducing agent that mitigates disulfide stress and promotes cell survival (25-27). EPAS1 overexpression elevated the relative NADP+/NADPH ratio in the cells (Figure 8G), indicating an increased risk of disulfidoptosis.

Figure 8 Functional role of EPAS1 in retinoblastoma cells. (A,B) Successful overexpression of EPAS1 in HXO-Rb44 cells confirmed by qRT-PCR (A) and Western blot (B). (C) Cell viability of HXO-Rb44 cells following EPAS1 overexpression, assessed by CCK-8 assay. (D) Migration capacity of HXO-Rb44 cells after EPAS1 overexpression, assessed by wound healing assay. (E) Invasion capacity of HXO-Rb44 cells following EPAS1 overexpression, assessed using Transwell assay. The invasive cells were stained by crystal violet. Scale bar: 100 μm. (F) Apoptosis of HXO-Rb44 cells after EPAS1 overexpression, determined by flow cytometry. (G) Relative NADP+/NADPH ratio in HXO-Rb44 cells following EPAS1 overexpression. **, P<0.01. CCK-8, cell counting kit-8; FITC-A, fluorescein isothiocyanate-area; LL, lower left; LR, lower right; NADH, nicotinamide adenine dinucleotide phosphate; NADPH, reduced nicotinamide adenine dinucleotide phosphate; NC, negative control; oe, overexpression; PE-A, phycoerythrin-area; Q, quadrant; qRT-PCR, quantitative real-time reverse transcription polymerase chain reaction; UL, upper left; UR, upper right.

Discussion

Retinoblastoma is a malignant paediatric tumour and a leading ocular disease that causes blindness and disability in childhood (12). Early diagnosis is critical, as timely treatment can cure the tumour and preserve both the eye and vision (28). To facilitate early diagnosis, timely treatment, and reduced mortality in retinoblastoma, the identification of key genes is particularly important. Disulfidoptosis, a newly recognised form of cell death (29), has not yet been explored in retinoblastoma, and the role of DRGs remains largely unclear. Therefore, this study aimed to elucidate the functions of DRGs in retinoblastoma using gene expression profiling.

Currently, the availability of high-throughput data resources is rapidly expanding, and numerous researchers now employ a range of bioinformatics methods to mine large-scale datasets in search of key genes implicated in the pathogenesis of various diseases. These efforts aim to identify novel biomarkers relevant to disease diagnosis and prognosis (30-32). Accordingly, bioinformatics approaches were utilised in the present study to identify key genes associated with retinoblastoma. Given the low incidence of retinoblastoma and the associated challenges in obtaining tissue samples, the sample size of individual datasets is typically small, which may compromise the accuracy of analytical results. To address this limitation, we merged two independent datasets to increase the overall sample size. Based on this combined dataset, we first calculated a disulfidoptosis score for each sample and found that the score was significantly lower in retinoblastoma samples than in control samples. Further differential expression analysis revealed that eight of the sixteen DRGs were dysregulated in retinoblastoma. These findings suggest a potential involvement of disulfidoptosis in the pathogenesis of retinoblastoma.

Tumour subtyping of heterogeneous retinoblastoma is considered critically important for developing optimal treatment strategies and predicting disease progression risk (33). Previous studies have identified two main molecular subtypes of retinoblastoma (10,34). One subtype exhibits few genetic alterations, a high expression score for mature cone markers, and maintains a cone-differentiated phenotype. In contrast, the other subtype harbours frequent recurrent genetic alterations, elevated expression of ganglion cell markers, and displays a more aggressive behaviour (10,34). Consistently, Zeng et al. (35) also identified these two primary subtypes based on gene expression profiles of retinoblastoma patients and further delineated several subgroups within each, differing in photoreceptor gene expression, retinal markers, and immune response. These findings underscore the distinct cells of origin in retinoblastoma, which may derive from either retinal progenitor cells or cone photoreceptor precursors (36). In the present study, two disulfidoptosis-related molecular subtypes of retinoblastoma were identified, which differed significantly in terms of signalling pathways and immune characteristics. Notably, the C2 subtype demonstrated activation of multiple hallmark oncogenic pathways, including MYC targets, glycolysis, the p53 pathway, and mTORC1 signalling, all of which are well-established tumourigenesis-associated pathways (37-40). In addition, the C2 subtype exhibited significantly lower levels of immune cell infiltration, as well as reduced immune and stromal scores and decreased expression of various immune checkpoint genes, indicating higher tumour purity. Both the pathway activation and immune landscape suggest that the C2 subtype may represent a more aggressive and immunosuppressive tumour phenotype. However, due to the lack of clinical and survival data, we were unable to correlate our findings with patient pathological features or prognostic outcomes.

Machine learning is a rapidly advancing field in medicine that offers superior predictive performance compared with conventional statistical models by capturing non-linear relationships between predictive factors and outcomes, as well as complex interactions among predictive variables (41,42). Owing to these advantages, machine learning approaches are being increasingly applied in the medical field (43). In the present study, two machine learning algorithms were employed, which identified two key genes in retinoblastoma: EPAS1 and SLC7A11. SLC7A11 is a key regulator of disulfidoptosis. Specifically, under glucose starvation, excessive accumulation of disulfides in SLC7A11-high-expressing tumour cells induces aberrant disulfide bond formation among actin cytoskeleton proteins, resulting in F-actin network collapse and subsequent cell death (44). Overexpression of SLC7A11 has been frequently reported in various human cancers (45), and disulfidoptosis can be triggered in these cells by glucose transporter inhibitors, thereby suppressing tumour growth (44). These findings suggest novel therapeutic avenues that could exploit such metabolic vulnerabilities. However, the expression of SLC7A11 decreased in retinoblastomas in this study. We speculate that SLC7A11 may function differently in retinoblastoma due to tumour-specific metabolic regulation. For example, a recent study has shown that SLC7A11 is involved in mediating the effects of arginine deprivation during the growth of retinoblastoma cells (24). Arginine deprivation effectively inhibits RB cell growth in vitro. During arginine deprivation, the activation of SLC7A11 is regulated by GCN2, and knockdown of SLC7A11 makes retinoblastoma cells partially resistant to arginine deprivation (24). EPAS1, also known as hypoxia-inducible factor 2 alpha (HIF-2α), is primarily expressed in endothelial cells (46). Hypoxia is a fundamental hallmark of tumours. Intra-tumoural hypoxia induces the expression of various genes through transcriptional activation mediated by HIF-1 and HIF-2, thereby contributing to the establishment of an immunosuppressive tumour microenvironment (47). Unlike HIF-1α, the role of EPAS1/HIF-2α in human cancers remains less well characterised. Although HIF-2α activation has been proposed as a common feature of tumours (48), some studies report a tumour-suppressive role in neuroblastoma (49), lung cancer (50), thyroid carcinoma (51) and glioma (52). In this study, we found that EPAS1 was down-regulated in retinoblastoma, consistent with the reduced expression of EPAS1/HIF-2α in papillary thyroid carcinoma (51). Therefore, we hypothesised that EPAS1/HIF-2α may play a tumour-protective role in retinoblastoma, probably due to the tissue-specific hypoxia response. Our in vitro experiments confirmed that EPAS1 functions as a tumour suppressor in retinoblastoma and that its overexpression was associated with the inhibition of retinoblastoma cell growth.

We found that the expression of both EPAS1 and SLC7A11 was associated with oxidative phosphorylation and immune cell infiltration in retinoblastoma. The oxidative phosphorylation metabolic pathway consumes oxygen to generate ATP through the transfer of electrons across a series of transmembrane protein complexes—collectively known as the electron transport chain—located in the inner mitochondrial membrane (53). Cancer cells frequently exhibit altered energy metabolism to support their high proliferative rates, typically characterised by enhanced glycolysis rather than reliance on mitochondrial respiration and oxidative phosphorylation (54). Consequently, oxidative phosphorylation has emerged as a promising target for cancer therapy (55). Immune cell infiltration is also implicated in the progression of retinoblastoma. For instance, Sarver et al. (56) reported a robust immune gene expression signature and increased infiltration of monocytes, macrophages, dendritic cells, and T lymphocytes in early-onset retinoblastoma, whereas these immune features were absent in normal retinae and late-onset tumours. Myeloid-derived suppressor cells have also been shown to play key roles in the progression and invasion of retinoblastoma (57). The significant associations of EPAS1 and SLC7A11 with oxidative phosphorylation and immune cell infiltration underscore their potential involvement in retinoblastoma pathogenesis. As the roles of SLC7A11 and EPAS1 in retinoblastoma remain largely unexplored, further investigations are warranted to elucidate their precise functions in this malignancy.

Although we demonstrated for the first time that DRGs are involved in the tumourigenesis and tumour microenvironment of retinoblastoma, we acknowledge some limitations of the current study. While we combined the two datasets to generate a larger dataset, the overall sample size remained small, particularly the number of control retinal samples, likely due to the inherent difficulty of collecting normal retinal tissues. Although the ComBat algorithm was applied to mitigate batch effects the merged dataset, the batch effect may not have been limited completely and other artefacts may have been introduced during this process. In addition, detailed clinical information about the patients included in the datasets was lacking, which limited tour ability to present a more comprehensive understanding of the clinical relevance of the identified subtypes and key genes. Disulfidoptosis is a relatively newly recognised form of cell death. Understanding of the genes that are directly and indirectly involved is still being elucidated. Therefore, several relevant genes may have been missed from this study. This seems to be an inherent limitation of bioinformatics analyses that rely on currently available data.

Conclusions

In summary, this study identified two disulfidoptosis-based molecular subtypes of retinoblastoma, each characterised by distinct oncogenic signalling patterns and immune infiltration profiles. Additionally, we identified EPAS1 and SLC7A11 as critical genes implicated in retinoblastoma progression. To our knowledge, this is the first study to report the involvement of disulfidoptosis in retinoblastoma. These findings may contribute to a more comprehensive understanding of the molecular heterogeneity and pathogenesis of retinoblastoma.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tp.amegroups.com/article/view/10.21037/tp-2024-596/rc

Data Sharing Statement: Available at https://tp.amegroups.com/article/view/10.21037/tp-2024-596/dss

Peer Review File: Available at https://tp.amegroups.com/article/view/10.21037/tp-2024-596/prf

Funding: This study was supported by the Project for Improving Evidence-Based Competence of the National Administration of Traditional Chinese Medicine Specialty (No. 2019XZZX-YK005).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tp.amegroups.com/article/view/10.21037/tp-2024-596/coif). Y.C. is a current employee of AIER Eye Hospital Group. The other authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Nag A, Khetan V. Retinoblastoma - A comprehensive review, update and recent advances. Indian J Ophthalmol 2024;72:778-88. [Crossref] [PubMed]
2. Chen J, Cao X, Xu S, et al. Global, regional, and national burden of retinoblastoma in infants and young children: findings from the global burden of disease study 1990-2021. EClinicalMedicine 2024;76:102860. [Crossref] [PubMed]
3. Schaiquevich P, Francis JH, Cancela MB, et al. Treatment of Retinoblastoma: What Is the Latest and What Is the Future. Front Oncol 2022;12:822330. [Crossref] [PubMed]
4. Pareek A, Kumar D, Pareek A, et al. Retinoblastoma: An update on genetic origin, classification, conventional to next-generation treatment strategies. Heliyon 2024;10:e32844. [Crossref] [PubMed]
5. The Global Retinoblastoma Outcome Study: a prospective, cluster-based analysis of 4064 patients from 149 countries. Lancet Glob Health 2022;10:e1128-40. [Crossref] [PubMed]
6. Utomo PT, Respatika D, Ardianto B, et al. Lag time, high-risk histopathological features, metastasis, and survival interrelation in retinoblastoma: a perspective from lower-middle income country. Int J Ophthalmol 2022;15:1994-2000. [Crossref] [PubMed]
7. Apumayta ED, Buitrago M, Rioja M, et al. Liver metastasis of retinoblastoma. Ecancermedicalscience 2025;19:1824. [Crossref] [PubMed]
8. Nag A, Khetan V. Genetics of Retinoblastoma - An Update. Semin Ophthalmol 2025; Epub ahead of print. [Crossref]
9. Roohollahi K, de Jong Y, van Mil SE, et al. High-Level MYCN-Amplified RB1-Proficient Retinoblastoma Tumors Retain Distinct Molecular Signatures. Ophthalmol Sci 2022;2:100188. [Crossref] [PubMed]
10. Liu J, Ottaviani D, Sefta M, et al. A high-risk retinoblastoma subtype with stemness features, dedifferentiated cone states and neuronal/ganglion cell gene expression. Nat Commun 2021;12:5578. [Crossref] [PubMed]
11. Yang J, Li Y, Han Y, et al. Single-cell transcriptome profiling reveals intratumoural heterogeneity and malignant progression in retinoblastoma. Cell Death Dis 2021;12:1100. [Crossref] [PubMed]
12. Rathore S, Verma A, Ratna R, et al. Retinoblastoma: A review of the molecular basis of tumor development and its clinical correlation in shaping future targeted treatment strategies. Indian J Ophthalmol 2023;71:2662-76. [Crossref] [PubMed]
13. Tong X, Tang R, Xiao M, et al. Targeting cell death pathways for cancer therapy: recent developments in necroptosis, pyroptosis, ferroptosis, and cuproptosis research. J Hematol Oncol 2022;15:174. [Crossref] [PubMed]
14. Peng F, Liao M, Qin R, et al. Regulated cell death (RCD) in cancer: key pathways and targeted therapies. Signal Transduct Target Ther 2022;7:286. [Crossref] [PubMed]
15. Zheng T, Liu Q, Xing F, et al. Disulfidptosis: a new form of programmed cell death. J Exp Clin Cancer Res 2023;42:137. [Crossref] [PubMed]
16. Chen H, Yang W, Li Y, et al. Leveraging a disulfidptosis-based signature to improve the survival and drug sensitivity of bladder cancer patients. Front Immunol 2023;14:1198878. [Crossref] [PubMed]
17. Wang T, Guo K, Zhang D, et al. Disulfidptosis classification of hepatocellular carcinoma reveals correlation with clinical prognosis and immune profile. Int Immunopharmacol 2023;120:110368. [Crossref] [PubMed]
18. Zhao S, Wang L, Ding W, et al. Crosstalk of disulfidptosis-related subtypes, establishment of a prognostic signature and immune infiltration characteristics in bladder cancer based on a machine learning survival framework. Front Endocrinol (Lausanne) 2023;14:1180404. [Crossref] [PubMed]
19. Tang J, Peng X, Xiao D, et al. Disulfidptosis-related signature predicts prognosis and characterizes the immune microenvironment in hepatocellular carcinoma. Cancer Cell Int 2024;24:19. [Crossref] [PubMed]
20. Sanz H, Valim C, Vegas E, et al. SVM-RFE: selection and visualization of the most relevant features through non-linear kernels. BMC Bioinformatics 2018;19:432. [Crossref] [PubMed]
21. Chen Y, Zhang Z, Hu X, et al. Epigenetic characterization of sarcopenia-associated genes based on machine learning and network screening. Eur J Med Res 2024;29:54. [Crossref] [PubMed]
22. Han X, Song D. Using a Machine Learning Approach to Identify Key Biomarkers for Renal Clear Cell Carcinoma. Int J Gen Med 2022;15:3541-58. [Crossref] [PubMed]
23. Daneshvar A, Mousa G. Regression shrinkage and selection via least quantile shrinkage and selection operator. PLoS One 2023;18:e0266267. [Crossref] [PubMed]
24. Wang D, Chu WK, Yam JCS, et al. GCN2-SLC7A11 axis coordinates autophagy, cell cycle and apoptosis and regulates cell growth in retinoblastoma upon arginine deprivation. Cancer Metab 2024;12:31. [Crossref] [PubMed]
25. Wang J, Wang M, Wu S, et al. Tumor suppressor BAP1 suppresses disulfidptosis through the regulation of SLC7A11 and NADPH levels. Oncogenesis 2024;13:31. [Crossref] [PubMed]
26. Wang K, Li L, Liang G, et al. Sonodynamic activated nanoparticles with Glut1 inhibitor and cystine-containing polymer stimulate disulfidptosis for improved immunotherapy in bladder cancer. Biomaterials 2025;319:123178. [Crossref] [PubMed]
27. Yao HF, Ge J, Chen J, et al. CASC8 activates the pentose phosphate pathway to inhibit disulfidptosis in pancreatic ductal adenocarcinoma though the c-Myc-GLUT1 axis. J Exp Clin Cancer Res 2025;44:26. [Crossref] [PubMed]
28. Sheehan AP. Retinoblastoma: Early Diagnosis is Crucial. J Pediatr Health Care 2020;34:601-5. [Crossref] [PubMed]
29. Zheng P, Zhou C, Ding Y, et al. Disulfidptosis: a new target for metabolic cancer therapy. J Exp Clin Cancer Res 2023;42:103. [Crossref] [PubMed]
30. Zhang W, Li K, Li S, et al. High-throughput sequencing reveals hub genes for human early embryonic development arrest in vitro fertilization: a pilot study. Front Physiol 2023;14:1279559. [Crossref] [PubMed]
31. Li X, Li M, Zhao T, et al. The discovery of promising candidate biomarkers in kidney renal clear cell carcinoma: evidence from the in-depth analysis of high-throughput data. Am J Cancer Res 2023;13:4288-304.
32. Soni U, Singh A, Soni R, et al. Identification of candidate target genes of oral squamous cell carcinoma using high-throughput RNA-Seq data and in silico studies of their interaction with naturally occurring bioactive compounds. J Biomol Struct Dyn 2024;42:8024-44. [Crossref] [PubMed]
33. Kooi IE, Mol BM, Moll AC, et al. Loss of photoreceptorness and gain of genomic alterations in retinoblastoma reveal tumor progression. EBioMedicine 2015;2:660-70. [Crossref] [PubMed]
34. Kapatai G, Brundler MA, Jenkinson H, et al. Gene expression profiling identifies different sub-types of retinoblastoma. Br J Cancer 2013;109:512-25. [Crossref] [PubMed]
35. Zeng Q, Wang S, Chen L, et al. Transcriptome analysis reveals molecularly distinct subtypes in retinoblastoma. Sci Rep 2023;13:16475. [Crossref] [PubMed]
36. Xu XL, Singh HP, Wang L, et al. Rb suppresses human cone-precursor-derived retinoblastoma tumours. Nature 2014;514:385-8. [Crossref] [PubMed]
37. Dhanasekaran R, Deutzmann A, Mahauad-Fernandez WD, et al. The MYC oncogene - the grand orchestrator of cancer growth and immune evasion. Nat Rev Clin Oncol 2022;19:23-36. [Crossref] [PubMed]
38. Zhao M, Wang T, Gleber-Netto FO, et al. Mutant p53 gains oncogenic functions through a chromosomal instability-induced cytosolic DNA response. Nat Commun 2024;15:180. [Crossref] [PubMed]
39. Marcucci F, Rumio C. On the Role of Glycolysis in Early Tumorigenesis-Permissive and Executioner Effects. Cells 2023;12:1124. [Crossref] [PubMed]
40. Armijo ME, Campos T, Fuentes-Villalobos F, et al. Rheb signaling and tumorigenesis: mTORC1 and new horizons. Int J Cancer 2016;138:1815-23. [Crossref] [PubMed]
41. Greener JG, Kandathil SM, Moffat L, et al. A guide to machine learning for biologists. Nat Rev Mol Cell Biol 2022;23:40-55. [Crossref] [PubMed]
42. Lo Vercio L, Amador K, Bannister JJ, et al. Supervised machine learning tools: a tutorial for clinicians. J Neural Eng 2020;
43. Heo J, Yoon JG, Park H, et al. Machine Learning-Based Model for Prediction of Outcomes in Acute Stroke. Stroke 2019;50:1263-5. [Crossref] [PubMed]
44. Liu X, Nie L, Zhang Y, et al. Actin cytoskeleton vulnerability to disulfide stress mediates disulfidptosis. Nat Cell Biol 2023;25:404-14. [Crossref] [PubMed]
45. Koppula P, Zhuang L, Gan B. Cystine transporter SLC7A11/xCT in cancer: ferroptosis, nutrient dependency, and cancer therapy. Protein Cell 2021;12:599-620. [Crossref] [PubMed]
46. Wang N, Hua J, Fu Y, et al. Updated perspective of EPAS1 and the role in pulmonary hypertension. Front Cell Dev Biol 2023;11:1125723. [Crossref] [PubMed]
47. Semenza GL. Targeting intratumoral hypoxia to enhance anti-tumor immunity. Semin Cancer Biol 2023;96:5-10. [Crossref] [PubMed]
48. Toledo RA, Jimenez C, Armaiz-Pena G, et al. Hypoxia-Inducible Factor 2 Alpha (HIF2α) Inhibitors: Targeting Genetically Driven Tumor Hypoxia. Endocr Rev 2023;44:312-22. [Crossref] [PubMed]
49. Westerlund I, Shi Y, Toskas K, et al. Combined epigenetic and differentiation-based treatment inhibits neuroblastoma tumor growth and links HIF2α to tumor suppression. Proc Natl Acad Sci U S A 2017;114:E6137-46. [Crossref] [PubMed]
50. Mazumdar J, Hickey MM, Pant DK, et al. HIF-2alpha deletion promotes Kras-driven lung tumor development. Proc Natl Acad Sci U S A 2010;107:14182-7. [Crossref] [PubMed]
51. Zhang R, Zhao J, Zhao L. EPAS1/HIF-2α Acts as an Unanticipated Tumor-Suppressive Role in Papillary Thyroid Carcinoma. Int J Gen Med 2023;16:2165-74. [Crossref] [PubMed]
52. Acker T, Diez-Juan A, Aragones J, et al. Genetic evidence for a tumor suppressor role of HIF-2alpha. Cancer Cell 2005;8:131-41. [Crossref] [PubMed]
53. Nolfi-Donegan D, Braganza A, Shiva S. Mitochondrial electron transport chain: Oxidative phosphorylation, oxidant production, and methods of measurement. Redox Biol 2020;37:101674. [Crossref] [PubMed]
54. Chelakkot C, Chelakkot VS, Shin Y, et al. Modulating Glycolysis to Improve Cancer Therapy. Int J Mol Sci 2023;24:2606. [Crossref] [PubMed]
55. Ashton TM, McKenna WG, Kunz-Schughart LA, et al. Oxidative Phosphorylation as an Emerging Target in Cancer Therapy. Clin Cancer Res 2018;24:2482-90. [Crossref] [PubMed]
56. Sarver AL, Xie C, Riddle MJ, et al. Retinoblastoma tumor cell proliferation is negatively associated with an immune gene expression signature and increased immune cells. Lab Invest 2021;101:701-18. [Crossref] [PubMed]
57. Li Q, Han J, Yang Y, et al. PD-1/PD-L1 checkpoint inhibitors in advanced hepatocellular carcinoma immunotherapy. Front Immunol 2022;13:1070961. [Crossref] [PubMed]