We first obtained ROS-related genes from the GeneCards database (https://www.genecards.org/) and Molecular Signature Database v7.1 (MSigDB; https://www.gsea-msigdb.org/gsea/msigdb). Then, we downloaded the level-three transcriptome RNA sequencing information and clinicopathological features of BLCA patients by accessing The Cancer Genome Atlas (TCGA) (https://gdc‐portal.nci.nih.gov/). Further, we utilized the GSE13507 acquired from the Gene Expression Omnibus database (GEO, https://www.ncbi.nlm.nih.gov/geo/) as the validation set.

Employing the R edge package (version R 4.0.5, https://bioconductor.org/packages/release/bioc/), the ROS‐related DEGs between BLCA and normal bladder samples were screened, setting the cutoff criteria as a False Discovery Rate (FDR) < 0.05 and a |log2-fold change (FC)| > 1.

Gene Ontology (GO) analysis that includes molecular function (MF), cell component (CC), and biological process (BP) analyses was implemented to explore the possible molecular mechanisms behind ROS‐related DEGs via the clusterProfiler package of R, utilizing the same approach for Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis (26–28) and considering P < 0.05 as significant enrichment.

ROS‐related DEGs were uploaded to the STRING database (http://www.string-db.org/) to obtain PPI information. The PPI network establishment and visualization were conducted via Cytoscape software, using the MCODE plug-in to screen the considerable PPI network modules.

The Connectivity Map (CMAP) database (http://www.broadinstitute.org) could be beneficial for researchers in the identification of probable molecular drugs closely associated with diseases, including cancer. The enrichment scores were -1–1, with a negative score showing that BLCA patients could benefit from this drug.

The prognosis-associated ROS-related genes were identified via Univariate Cox regression (UCR) analysis (survival package) and least absolute shrinkage and selection operator (LASSO) regression analysis (glmnet and survival package) with P < 0.05 in the TCGA dataset, followed by incorporating the results into the Multivariate Cox Regression (MCR) analysis. Finally, a ROS-correlated gene signature related, to OS was constructed based on MCR analysis results. The Risk Score (RS) was generated by this formula: RS = (Coef1*expression mRNA1) + (Coef2*expression mRNA2) + (Coef n * expression mRNA n), where Coef represents the MCR model coefficient of relevant mRNA. Based on the RS mean, patients were classified into high-risk group (HRG) and low-risk group (LRG), employing the Kaplan-Meier (K-M) method to compare the survival outcomes between different groups. Our study deployed time-related ROC analysis to determine the predictive prognostic value of the PS. Both T-distributed stochastic neighbor embedding (t-SNE) analysis alongside principal component analysis (PCA) were implemented to examine the risk signature classification capacity with the R packages “Rtsne” and “ggplot2”, employing the same approach to calculate the RS and then validated the ROS-related gene signature in the GSE13507 dataset.

We explored the relationships between the PS and clinical features (age, sex, T/N/TNM stages, and tumor grade) in the TCGA dataset via the chi-squared test. Then, stratified analysis was performed to further examine the PS reliability and stability of ROS in the prediction of BLCA patients’ OS. Additionally, we implemented UCR and MCR analyses to explore whether the RS was of independent prognostic value. Both RS and clinical features were incorporated to establish an OS-related nomogram, estimating the nomogram’s predictive capability by generating a calibration curve.

GSEA was implemented to investigate the latent mechanisms among different groups based on GSEA software (version 4.1.0). Then, we acquired mutation information for BLCA patients by accessing the TCGA database, calculated the total mutation number for each sample, and analyzed the top mutational genes among different risk groups using the maftools package. The TIMER, CIBERSORT, CIBERSORT-ABS, XCELL, QUANTISEQ, EPIC, and MCP-counter methods were utilized for the analysis of the ICI levels of 22 distinct leukocyte subsets in both groups. P < 0.05 indicated statistically significant.

The Genomics of Drug Sensitivity in Cancer (GDSC, http://www.cancerrxgene.org) database was accessed to predict BLCA patients’ response in both groups to chemotherapy drugs. Eventually, we assessed chemosensitivity by calculating the half-maximal inhibitory concentration (IC50) through the R package pRRophetic, with P < 0.05 indicating statistical significance.

Between 2022 and 2024, 20 BLCA tissue and their corresponding non-tumor tissue specimens were collected from Shaoxing People’s Hospital for immunohistochemical staining (IHC) and western blot. No patient in this study had received radiation therapy, adjuvant therapy, or preoperative chemotherapy. The samples from Shaoxing People’s Hospital were collected with informed consent, and the use of the stored cancer specimens and clinical data was granted clearance by the Academic Ethical Committee of Shaoxing People’s Hospital (ethical approval number: 2022-K-Y-054-01). The study was executed in a way that aligned with the Declaration of Helsinki.

IHC images of key genes in BLCA and normal tissue samples were acquired through the Human Protein Atlas (HPA) database while evaluating the staining intensity following the HPA database standard (https://www.proteinatlas.org/). Use anti-AKR1B1 (1:1000; Proteintech,15439-1-AP). Rabbit monoclonal antibody was used for immunohistochemistry of paraffin-embedded human and nude mouse BLCA specimens. In short, samples were processed using dewaxing, hydration, antigen extraction, IHC labeling, and pathology scores.

Human BLCA cells (T24 and 5637) were procured from Procell Life Science & Technology Company (Hubei, China). Herein, we grouped the logarithmic growth phase cells into the control, siAKR1B1-negative control (NC), and siAKR1B1 groups. The two cell lines were cultured in MCCOY’S 5A (Gibco, USA) and 1640 (Gibco, USA) medium supplemented with 10% fetal bovine serum (FBS, Gibco, USA) and 1% penicillin/streptomycin at 37°C and 5% CO₂. The cells went through treatment with 5 µl of siAKR1B1, using Lipofectamine 2000 to dilute the solution in Opti-MEM for 5 min. The solution was thereafter mixed and allowed to incubate at ambient temperature for a duration of 20 min, followed by introducing the composite into the cell culture plate. After a 48-h period of transfection, the cells were gathered for additional assessments.

Proteins were subjected to extraction using RIPA buffers and quantification by BCA kits. The protein eluate went through separation utilizing 10% SDS-PAGE and transferred to a PVDF membrane that was blocked and then incubated with primary and secondary antibodies. AKR1B1 (Proteintech, 15439-1-AP) and β-catenin (Abcam, ab32572) were detected by imaging with enhanced chemiluminescence reagents (Merck Millipore, Billerica, MA).

100 μL of suspension containing 5000 transfected cells was dispensed into each well of a 96-well plate along with 10 μL of CCK8 solution (MCE, HY-K0301). The plate was then placed in a cell culture incubator for 1 hour, following which absorbance readings were taken at 450 nm and recorded.

3600 BLCA cells were equally distributed into six-well plates and incubated at 37°C with 5% CO₂ for 14 days with regular medium changes. Post-incubation, the cells were fixed and stained using 4% paraformaldehyde and 0.1% crystal violet for 20 minutes each, after which images were captured and data documented.

5-ethynyl-2′-deoxyuridine (EdU) assay kit (MCE, China) was used as instructed by the manufacturer. In this experiment, BLCA cells were cultured in 96-well plates, with a seeding density of 4,000 cells per well, after incubation at 37°C for 72 hours. Next, BLCA cells were exposed to 10 μM EdU for 2 hours at 37°C. Subsequently, the cells were fixed using 4% paraformaldehyde and permeabilized with 0.5% Triton X-100 for 15 minutes at room temperature. After removing the fixatives, the cells were washed with PBS containing 1% BSA. Lastly, the cells were incubated in Click Additive Solution, protected from light, for 30 minutes and then stained with Hoechst to label the nucleus. Microscopic images were captured to observe the EdU detection samples. The proliferation of cells was further assessed by calculating the ratio of EdU-positive cells to the overall cell count.

The transfected HOS and 143B cell lines were cultured with the serum-free DMEM and serum-free 1640, respectively, in a Transwell upper chamber. Corresponding culture medium containing 10% FBS was added to the lower chamber. The cells were incubated at 37°C with 5% CO2 for 48 hours, fixed with formaldehyde, stained with crystal violet, and visualized under a microscope for analysis.

Statistical analysis was conducted using the R software (version 4.0.5). The significance of differential gene expression was ascertained using adjusted p-value to correct for the multiple testing phenomenon, with a significance threshold set at p-value < 0.05. Another statistical analysis was conducted with SPSS Statistics software version 20. The values were compared by one-way ANOVA or independent-samples Student’s t test. Statistical significance was determined at *p < 0.05, **p < 0.01, or ***p < 0.001. Values are presented as the mean ± SEM. Error bars indicate the SEM unless otherwise noted.

Supplementary Figure 1 illustrates the workflow diagram of this study. In total, we obtained 1749 ROS-related genes with relevance scores > 0.5 from the Gene Cards database and 70 ROS-related genes from the MSigDB database, acquiring 1767 genes after removing the overlapping genes. However, from the 1,767 genes, we eventually extracted the expression profiles of 1,719 ROS-associated genes identical to those in 412 and 19 BLCA and normal bladder tissue samples, respectively, in the TCGA dataset. By applying cutoff criteria of FDR < 0.05 and |log2 FC| > 1, 308 ROS-related Differentially Expressed Genes (DEGs) were identified; of them, 138 were downregulated, and 170 were upregulated. Moreover, GO, KEGG, and PPI analyses were deployed to explore the possible roles of ROS-associated genes. Both Univariate Cox regression (UCR) and LASSO analyses were performed to screen for prognostic ROS-related genes, and 71 genes were included in subsequent analyses (p<0.05). As a means to guarantee the clinical outcomes’ stability and reliability based on the 71 genes, we conducted LASSO analysis to further screen for prognostic ROS-related genes, identifying 31 genes related to OS. The MCR analysis identified 17 ROS-related genes (JUN, CALR, P4HB, ELN, MYC, FASN, REV3L, VHL, NID1, SLC38A1, TFRC, AKR1B1, ITGA3, CGB5, HLA-G, FADS1, and ORM1) that were utilized to construct a PS, which was subsequently validated in the GSE13507 cohort. Both K–M survival and time-dependent receiver operating characteristic (ROC) curves were employed to evaluate the prognostic value of the PS. A nomogram was constructed, aiming at predicting the outcomes of BLCA patients in combination with the PS and clinical factors. GSEA, ICI, and TMB analysis were implemented for the exploration of the molecular characteristics of the PS. The GDSC database was accessed for the prediction of chemotherapy response according to the PS. Candidate SMDs targeting ROS-related genes were screened against the CMAP database. To verify the authenticity of the data, the PS protein expression levels were detected through the HPA. AKR1B1 was selected for in vitro experimental validation, demonstrating the reliability of the ROS-related genes we constructed.

The GO analysis findings represented that ROS-associated genes were involved in multiple biological processes, including the response to toxic substances, aging, metal ions, oxidative stress, and ROS, besides cell cycle arrest and the cellular response to drugs ( Figure 1A ). The KEGG analysis findings showcased that these genes exhibited main involvement in multiple pathways, including the p53, platinum drug resistance, cell cycle, ErbB, PI3K-Akt, TNF, cellular senescence, IL-17, MAPK, HIF-1, and cGMP-PKG signaling pathways ( Figure 1B ) (26–28). For better comprehension of the involvements of ROS-associated genes in BLCA, a PPI network was established and visualized through the utilization of STRING database and Cytoscape software, which included 298 nodes and 2859 edges ( Figure 2A ). The MCODE plugin identified three crucial modules of target genes, and the critical modules consisted of 39 nodes and 321 edges, 29 nodes and 250 edges, and 31 nodes and 123 edges ( Figures 2B–D ).

Using the CMAP database, candidate SMDs for BLCA were identified based on ROS-related DEGs, identifying eleven SMDs (0297417-0002B, 5248896, puromycin, blebbistatin, anisomycin, STOCK1N-35215, methylergometrine, clofilium tosylate, verteporfin, withaferin A, and rottlerin) with anticancer functions in BLCA progression with enrichment scores < -0.8, p < 0.01, and n > 2 as the screening criteria ( Table 1 ).

Our study conducted UCR analysis to identify ROS-related DEGs notably correlated with OS, and 71 genes were included in subsequent analyses (p<0.05) ( Figure 3A ). Aiming to ensure the clinical outcomes stability and reliability based on the 71 genes, LASSO analysis was conducted to further screen for prognostic ROS-related genes, and we identified 31 genes related to OS ( Figures 3B, C ). Multivariate Cox Regression (MCR) analysis identified 17 ROS-related genes (JUN, CALR, P4HB, ELN, MYC, FASN, REV3L, VHL, NID1, SLC38A1, TFRC, AKR1B1, ITGA3, CGB5, HLA-G, FADS1, and ORM1) that were used to construct a PS ( Figure 3D ). A ROS-based RS was established depending on the coefficient of 17 genes according to this formula: risk score = (0.3078 × FASN expression) + (0.305 × CALR expression) + (0.3832 × P4HB expression) + (0.1599 × ELN expression) + (0.2941 × MYC expression) + (0.3702 × REV3L expression) + (-0.4548 × VHL expression) +(0.147 × NID1 expression) + (-0.2213 × SLC38A1 expression) + (0.1687 × TFRC expression) +(0.123 × AKR1B1 expression) + (-0.1371 × ITGA3 expression) +(0.1762 × CGB5 expression) + (-0.1483 × HLA-G expression) +(0.1368 × FADS1 expression) + (-0.252 × ORM1 expression) + (0.1274 × JUN expression). Subsequently, we classified patients into HRG and LRG in accordance with the median RS. The LRG patients had longer OS than those in the HRG (p < 0.05) ( Figures 4A, C ). Time‐dependent ROC analysis depicted that the signature AUC in the TCGA cohort was 0.78 at 5 years ( Figure 4B ). A heatmap was generated to show the differences in 17 ROS-related genes between the different groups ( Figure 4D ). PCA and t-SNE analyses indicated the signature’s good classification ability ( Figures 4E, F ). Additionally, the prognostic capacity of our constructed PS was validated in the GSE13507 dataset. The results in GSE13507 ( Supplementary Figure 2 ) were consistent with previous results, which demonstrated the good performance of the PS in predicting OS.

We first compared the ROS-based RSs among various subgroups classified by clinicopathological characteristics (TNM stage, sex, grade, age, T stage, and N stage). The RSs exhibited a significant correlation with clinicopathological factors and were markedly elevated in the following subgroups: >65 years of age, advanced T stage (T3/4), N stage (N1/2/3), pathological grade (High), and TNM stage (Stage III-IV) ( Figure 5 ). Subsequently, stratification analysis was conducted relying upon the clinical characteristics (age, sex, grade, and TNM/T/N stages). Male or female sex, age (>65 years) or (<=65 years), T stage (T3-T4), N stage (N0), pathological grade (High), and TNM stage (Stage III–IV) were associated with inferior OS in the high-risk subgroup (P < 0.05) ( Supplementary Figure 2 ), with no difference in OS in the T stage (T1/2), N stage (N1/2/3), or TNM stage (Stage I-II) subgroup ( Figure 6 ). Furthermore, to evaluate whether the RS was an autocephalous prognostic indicator for BLCA patients, univariate and multivariate Cox proportional hazard models were implemented. According to UCR analysis results, age, TNM/T/N stages, and RS were related to unfavorable OS ( Figure 7A ). According to the multivariate analysis, age, N stage, and RS were still associated with unfavorable OS ( Figure 7B ). Multiparameter ROC curve analyses also revealed that the AUC of the RS was 0.769 ( Figure 7C ), indicating that compared with traditional clinical prognostic indicators, the ROS-based RS exhibited remarkable performance in predicting prognosis. Collectively, the ROS-based RS was an autocephalous prognostic factor. A nomogram including RS, age, and N stage we established to predict the outcomes of BLCA patients ( Figure 8A ), with the calibration curve elucidating the nomogram’s good performance in predicting patient prognosis ( Figures 8B, C ).

The GSEA results demonstrated that carcinogenic signaling pathways, such as calcium, focal adhesion, ECM receptor interaction, MAPK, BLCA, GAP junction, Wnt, Hedgehog, cancer, and TGF-β signaling pathways, exhibited main enrichment in the HRG ( Figure 9 ). Several metabolism-associated signaling pathways, including autophagy regulation, peroxisomes, and oxidative phosphorylation, were highly enriched in the LRG.

A heatmap of the ICI data obtained via CIBERSORT, MCPcounter, XCELL, TIMER, CIBERSORT-ABS, QUANTISEQ, and EPIC analyses ( Figure 10 ) suggested that the RS was correlated with ICI in BLCA. Additionally, significant differences were observed in the fractions of distinct leukocyte subsets between both groups. The proportions of naive B cells and M0/M2 macrophages were lower in the HRG, whereas the proportions of CD8+/CD4+T cells/were greater in the LRG.

The mutation profile results among different risk groups in the TCGA dataset revealed somatic mutations in 93.53% (118) and 94.55% (191) of the BLCA patients in the HRG ( Figure 11A ) and LRG ( Figure 11B ), respectively. TP53, TTN, KMT2D, MUC16, ARID1A, KDM6A, PIK3CA, SYNE1, RB1, and KMT2C were the top 10 mutated genes in the HRG. TP53, TTN, KMT2D, MUC16, ARID1A, KDM6A, PIK3CA, SYNE1, RB1, and FGFR3 were the top 10 mutated genes in the LRG. Furthermore, the proportions of somatic mutations in KDM6A and FGFR3 significantly differed between both groups. Additionally, the LRG patients had more mutation events than the HRG ( Figure 11C ). Patients having a high TMB appeared to possess a better prognosis than those with a low TMB ( Figure 11D ). Further, we investigated the collaborative interaction effect of the ROS-based RS and TMB on prognosis ( Figure 11E ). We found that the HRG patients having a high TMB had shorter OS than those in the LRG with a high TMB, and the LRG patients having a low TMB had longer OS than those in the HRG with a low TMB. Interestingly, patients having a high TMB exhibited better OS than those having a low TMB in the HRG, and patients having a low TMB displayed worse OS than those having a high TMB in the LRG. Patients possessing high TMB in the LRG had a greater OS than patients in the other three patient groups, and patients having low TMB in the HRG tended to have a significantly worse OS than patients in the other three patient groups. Collectively, the ROS-based RS might be a probable biomarker for predicting OS in BLCA patients.

The GDSC database analysis findings depicted that the IC50 values of chemotherapy drugs, including GSK269962A, BMS.536924, JNJ.26854165, docetaxel, temsirolimus, cisplatin, thapsigargin, sunitinib, rapamycin, and paclitaxel, were greater in LRG patients than in those the HRG. In comparison, the IC50 values of BIBW2992 and gefitinib were greater in HRG patients than in those at LRG ( Figure 12 ).

IHC was utilized to additionally investigate the nine gene protein expression in the HPA database between BLCA and normal control tissues ( Figure 13 ). In line with the RNA sequencing data, P4BH, FASN, AKR1B1, and CBG5 proteins, which have a high prognostic risk, were upregulated in tumor tissues, and MYC proteins, which have a low prognostic risk, were downregulated in tumor tissues compared with normal controls.

While AKR1B1 has been documented in other types of cancer (29–31), its impact on BLCA remains unreported. Therefore, AKR1B1 was selected for further analysis. IHC and WB analyses elucidated AKR1B1 overexpression in BLCA tissues ( Figure 14A ). To understand the role of AKR1B1 in BLCA, we further investigated the effect of increased AKR1B1 levels in BLCA cell lines (T24 and 5637) via in vitro experiments. siRNA transfection successfully interfered with the mRNA expression of AKR1B1, which was confirmed by WB ( Figure 14B ). To further understand the effect of AKR1B1, colony formation analysis was also performed ( Figure 14C ), which showed that BLCA cell viability was significantly hindered after AKR1B1 was silenced by siRNA. The CCK-8 assay showcased that cell viability was impeded after the silencing of AKR1B1 expression ( Figure 14D ). In addition, an EdU proliferation assay showed that inhibiting AKR1B1 significantly lowered the percentage of EdU-positive BLCA cells ( Figure 14E ). To further test whether AKR1B1 affects BLCA cell metastasis, a Transwell assay was performed ( Figure 14F ), revealing that siRNA-mediated silencing of AKR1B1 inhibited BLCA cell migration and invasion.