Cuproptosis-Related lncRNAs Modulate the Prognosis of MIBC by Regulating the Expression Pattern of Immunosuppressive Molecules Within the Tumor Microenvironment

Introduction

Bladder carcinoma (BLCA) ranks among the top ten most prevalent and costly cancers globally, accounting for approximately 550,000 new cases annually and significant morbidity and mortality.^1,2 Muscle-invasive bladder cancer (MIBC) is initially diagnosed in 25% of patients, while 10–20% of non-muscle invasive bladder cancer (NMIBC) cases progress to MIBC.³ Despite advancements in surgical and treatment modalities, MIBC prognosis remains poor, with a 5-year post-therapy survival rate of 40%-60%.^4,5 Investigating MIBC development mechanisms and identifying prognostic biomarkers are vital for enhancing patient outcomes and tailoring postoperative treatment strategies.

Cuproptosis, a recently discovered form of programmed cell death, involves copper, an essential cofactor in enzymatic activity and physiological homeostasis.^6–9 Previous studies have shown that c opper induces tumor cell death, inhibits tumor progression, contributes to proteotoxic stress, impedes proteasomes, and promotes anti-angiogenesis.^10–12 However, few studies have explored the BLCA microenvironment or its biological functions, warranting further investigation into cuproptosis’ relationship with BLCA development and progression. In the meantime, long non-coding RNA (lncRNA) dysregulation is associated with various tumor types,^13,14 but the complete role of cuproptosis-related gene regulatory factors in lncRNA regulation remains unclear. Uncovering the mechanisms by which cuproptosis-related gene modification influences MIBC development may reveal potential therapeutic targets.

Using The Cancer Genome Atlas (TCGA) and IMvigor210 cohorts, we analyzed MIBC patient data to determine the prognostic significance of cuproptosis-related lncRNA. Our findings indicated that six cuproptosis-related lncRNAs possess prognostic value in MIBC patients, and we constructed a prognostic signature based on their ability to predict overall survival. Crucially, our study suggests that cuproptosis-related lncRNA modulates immunosuppressive molecule expression, impacting MIBC prognosis.

Materials and Methods

Tissue Sample Collection and RT-qPCR Analysis

Ten pairs of fresh tumor tissue samples from MIBC patients and paired normal tissue samples were collected at the Xinhua Hospital Affiliated Shanghai Jiao Tong University School of Medicine. This research was approved by the Medical Ethics Committee of Xinhua Hospital of Affiliated Shanghai Jiao Tong University School of Medicine and informed consent was obtained from all patients. The study complies with the Declaration of Helsinki.

Total RNA was isolated from tissues using TRIzol reagent (Invitrogen, CA, USA) following the manufacturer’s protocol. Then, cDNA was used as a qPCR template to quantify genes using TB Green^® Premix Ex TaqTM II (Takara, Shanghai, China). The relative gene expression was normalized to β-actin according to the cycle threshold (2−ΔCt) method. The experiment was performed in triplicate, and the primer sequences employed for analysis are listed in Table 1.

Table 1 Primer Sequences of lncRNA

MIBC Datasets and Samples

The RNA-seq transcriptome data of patients with BLCA and the corresponding clinical data were downloaded from the TCGA data portal (https://portal.gdc.cancer.gov/). Finally, 396 MIBC patients were included for further analysis based on the simultaneous availability of mRNA expression profile data and complete clinical information. The IMvigor210 dataset was obtained from http://research-pub.gene.com/IMvigor210CoreBiologies/. The detailed clinical patient characteristics are summarized in Table 2. Patients with locally advanced or metastatic uroepithelial carcinoma were included in the study and further analyzed.

Table 2 Clinical Characteristics of BLCA Patients

Cuproptosis-Related Genes Detection

Expression matrixes of 19 cuproptosis-related genes were extracted from the TCGA databases based on previous publications, including NFE2L2, NLRP3, ATP7B, ATP7A, SLC31A1, FDX1, LIAS, LIPT1, LIPT2, DLD, DALAT, PDHA1, PDHB, MTF1, GLS, CDKN2A, DBT, GCSH, DLST.

Identification of Cuproptosis-Related lncRNAs

The long non-coding RNA annotation file of Genome Reference Consortium Human Build 38 (GRCh38) was obtained from the GENCODE website for the annotation of lncRNAs in both the TCGA dataset and the Imvigor210 dataset. Pearson correlation analysis was conducted in the TCGA and IMvigor210 datasets to identify cuproptosis-related lncRNAs. Following this, 71 survival-associated cuproptosis-related lncRNAs were selected for further investigation.

Establishment of the Prognostic Model of Cuproptosis-Related lncRNAs

The least absolute shrinkage and selection operator (LASSO) Cox regression algorithm was used to establish a prognostic model based on the normalized expression matrix of candidate prognostic 6 cuproptosis-related lncRNAs. The patient’s risk score was calculated according to the expression level of each gene and its corresponding regression coefficient. The formula was established as follows: . MIBC samples were divided into high-risk and low-risk groups based on the median value of the risk score. Ttryyyyhhe Kaplan-Meier (KM) curve, the receiver operating characteristic (ROC) and the principal component (PCA) were applied to evaluate the difference in the survival status of patients between the two groups. Finally, a nomogram model was established to enhance the predictive efficacy of our model in the TCGA dataset.

Screening and Functional Enrichment of Differential Genes Between the Two Groups

The differential expression genes (DEGs) were filtrated by “Limma” packages between‐group (high‐risk vs low‐risk). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses based on the DEGs (|log2FC| ≥ 1, FDR < 0.05) between the high-risk and low-risk groups.

Estimation of Infiltrating Cells in TME in MIBC Samples

Single-sample gene set enrichment analysis (ssGSEA) was used to calculate the infiltration score of 16 immune cells and the activity of 13 immune-related pathways.

Calculation of Tumor Mutation Burden Scores

Tumor mutational burden (TMB) data from BLCA samples downloaded from TCGA were analyzed using the “maftools” R package. The waterfall plot displayed the relationship between the risk signature model and TMB for BLCA patients.

Statistical Analysis

Kaplan-Meier analysis and Log rank test were used to compare OS between groups. Univariate and multivariate Cox regression analyses were used to determine independent predictors of OS. The chi-square test was used for the comparison of proportional differences. All statistical analysis is carried out using R software.

Results

Identification of Cuproptosis-Related lncRNA in MIBC Patients

To ascertain the biological function of cuproptosis-related lncRNAs in the development of MIBC, we identified these lncRNAs from the TCGA dataset and the IMvigor210 dataset. The study flowchart is depicted in Figure 1A. We identified 71 common cuproptosis-related lncRNAs in both the TCGA and IMvigor210 datasets. Subsequently, univariate Cox regression analysis revealed 13 cuproptosis-related lncRNAs that were significantly associated with survival prognosis (Figure 1B). We systematically examined the expression profiles of these 13 prognostic-related cuproptosis-related lncRNAs between MIBC and adjacent normal tissue pairs, finding that their expression patterns varied between tumor and normal tissues. Some cuproptosis-related lncRNAs were highly expressed in tumor tissues, while others were more highly expressed in normal tissues (Figure 1C). Our results support the notion that cuproptosis-related lncRNAs may constitute one of the components influencing the initiation and progression of MIBC.

Figure 1 Predictive modeling of Cupropsisi-related lncRNAs. (A) Experimental flow chart. (B) Cupropsisi-related lncRNAs associated with survival. (C) Expression of Cupropsisi-related lncRNAs associated with survival in tumor tissue versus normal tissue. (D-E) The predictive modeling of Cupropsisi-related lncRNAs was constructed via lasso regression. (F) Association of Cupropsisi-related lncRNAs with Cupropsisi-related genes. *P < 0.05, **P < 0.01, ***P <0.001.

Construction and Validation of Prognostic Signatures for Cuproptosis-Related lncRNAs

To explore the contribution of cuproptosis-related lncRNAs in the prognosis of MIBC patients, 396 MIBC patients were randomly assigned to TCGA training cohort and validation cohort at a ratio of 5:5. Meanwhile, we performed LASSO regression analysis based on the expression values of 13 cuproptosis-related lncRNAs previously screened in the TCGA training cohort, and finally selected 6 cuproptosis-related lncRNAs, including FAM13A-AS1, GHRLOS, LINC00456, OPA1-AS1, RAP2C-AS1, and UBE2Q1-AS1 (Figure 1D–F).

Subsequently, patients were divided into either the high-risk group or the low-risk group according to the median value of the risk score of prognostic cuproptosis-related lncRNAs. In both the training and validation cohorts, overall survival (OS) was longer in the low-risk group than in the high-risk group (Figure 2A–C). Moreover, the survival rate was almost double in the low-risk group compared with that in the high-risk group (Figure 2D). The PCA analyses indicated that patients in different risk groups were distributed in three directions (Figure 2E). In addition, the predictive value of the signature for OS and clinical characteristics was evaluated using time-dependent ROC curves, showing that the predictive value of the risk signature model was higher than that of clinical traits (Figure 2F and G). For further evaluating the relationship between the level of lncRNAs expression and clinical features, the heatmap was used to display the relationship between the expression profiles of 6 cuproptosis-related lncRNAs and their clinical characteristics in all MIBC patients (Figure 2H).

Figure 2 Verification of the prognostic model. (A-B) Survival curve between low-risk and high-risk groups in the training cohort and validation cohort. (C) Distribution of the risk score, OS status, and heatmap of the expression of Cuproptosis-related lncRNAs. (D) Comparison of survival between the high- and low-risk groups. (E) Validation of the model using PCA (F and G) The AUC of ROC curves showed the predictive abilities of the prognostic model. (H) The heatmap displayed the relationship between the expression of cuproptosis-related lncRNAs and clinicopathologic features in high- and low-risk groups.

Construction of the Cuproptosis-Related lncRNAs Nomogram

Univariate and multivariate Cox regression analysis revealed that the risk score of cuproptosis-related lncRNAs was a substantial independent prognostic factor for overall survival in the TCGA dataset (Figure 3A and B). To improve the accuracy of the prediction model, we constructed a nomogram for MIBC based on the risk signature and clinical characteristics. The total score was the sum of the individual scores calculated based on the use of the nomogram, and the OS of the patient gradually diminish as the score increased. Figure 3C illustrates an example of employing the nomogram to predict the survival probability for a randomly selected patient. ROC analysis demonstrated that the nomogram surpassed other factors in prognosticating MIBC outcomes (Figure 3D). Additionally, the nomogram exhibited high concordance between the predicted and observed survival probabilities (Figure 3E).

Figure 3 The predictive of the model was verified in the nomogram model. Univariate (A) and multivariate (B) analyses confirmed that risk scores were independent prognostic factors. (C) Nomogram combining the risk score model. (D) The ROC curves of the nomogram. (E) Calibration plots of 1-, 3-, and 5-year survival probabilities. *P < 0.05, **P < 0.01, ***P <0.001.

Functional Annotation of Cuproptosis-Related lncRNAs Risk Score

To investigate the function of cuproptosis-related lncRNAs, the “Limma” package was applied to analyze DEGs between the high‐risk and low‐risk groups, and 445 DEGs were filtrated. Interestingly, GO enrichment showed that DEGs were significantly enriched in many areas related to the regulation of immune responses, including antigen binding, antibody production, humoral immune response, and lymphocyte-mediated immunity (Figure 4A). KEGG pathway analyses also indicated that the cytokine-cytokine receptor interaction, IL-17 signaling pathway, and TNF signaling pathway were enriched in these DEGs (Figure 4B).

Figure 4 The exploration of the function of cuproptosis-related lncRNAs. GO (A) enrichment and KEGG (B) pathways analysis for DEGs. Expression of tumor-infiltrating immune cells (C) and enrichment of immune signaling pathways (D). (E) The expression of immune checkpoint inhibitory molecules between the two groups. *P < 0.05, **P < 0.01, ***P <0.001. Abbreviation: ns, no significance.

GO and KEGG analysis displayed that DEGs in the risk score were enriched in immune interaction pathways. We further studied the relationship between the risk score of cuproptosis-related lncRNAs and tumor-infiltrating immune cells. Intriguingly, the results showed that both pro-tumor-related immune factors and anti-tumor-related immune parameters were upregulated in the high-risk group (Figure 4C and D). Subsequently, we analyzed the expression of tumor immunosuppressive checkpoints between the two groups and the majority of tumor immunosuppressive checkpoints were detected to be highly expressed in the high-risk group (Figure 4E). This observation may explain the limited clinical efficacy of single immune checkpoint inhibitors.

The Relationship Between Cuproptosis-Related lncRNAs Risk Scores and the Tumor Mutation Burden

Previously reported that the TMB can influence the infiltration of tumor immune cells.¹⁵ We explored the expression pattern of TMB in the TME between the high- and low-risk groups. The results depicted that the mutation rates of the top 20 most substantially variable genes were significantly distinct between the high- and low-risk groups, as shown in Figure 5A and B. Upon further analysis, we observed that a substantial number of patients in the low-risk group exhibited low TMB, and their overall survival (OS) was significantly compromised, indicating that low TMB adversely affected patient OS. This observation may also provide insight into why some patients in the previous low-risk group did not exhibit prolonged OS (Figure 5C and D). Finally, we integrated the risk signature, patient survival status, tumor stage, molecular subtype, and TMB, identifying a robust correlation among these factors (Figure 5E).

Figure 5 The association between risk models and TMB. Expression patterns of TMB between low- (A) and high-risk groups (B). (C–E) Survival analysis of TMB combined risk score.

Validation of Cuproptosis-Related lncRNAs and Association with Immunotherapy Response in IMvigor210 Cohort

Our prior investigations established a robust association between risk scores and immune checkpoints. To validate the constructed model, BLCA patients from the IMvigor210 cohort were selected. Employing the same formula, a corresponding risk score was calculated for each patient within the IMvigor210 cohort. The complete response (CP) and stable disease (SD) rates were higher in the low-risk score group compared to the high-risk score group (Figure 6A). Moreover, overall survival (OS) was prolonged in the low-risk score group relative to the high-risk score group (Figure 6B). Regrettably, the scores of TCGA subtypes did not exhibit significant differences, but the risk scores of patients with CR and SD were notably lower (Figure 6C and D). Utilizing the Lund subtype, the risk scores were highest in the genomically unstable subtype and lowest in the infiltrated subtype (Figure 6E and F).

Figure 6 Validation of fractional risk score model using IMvigor210 dataset. (A) Comparison of the anti-PD-L1 response rates between the high- and low-risk groups. (B) Survival curves between high - and low-risk groups. Risk scores among TCGA subtypes (C), immune response groups (D), and Lund subtypes (E and F).

Relative mRNA Expression of Cuproptosis-Related lncRNAs

The RT-qPCR validation results from MIBC patient tissues revealed differential expression of FAM13A-AS1, GHRLOS, LINC00456, OPA1-AS1, RAP2C-AS1, and UBE2Q1-AS1 between normal and tumor samples (Figure 7A–F). Corresponding P-values for FAM13A-AS1, GHRLOS, LINC00456, and OPA1-AS1 were 0.0098, 0.0273, 0.0137, and 0.0273, respectively. Although UBE2Q1-AS1 expression did not display significant differences, a trend was observed between tumor and normal tissues, potentially warranting sample size expansion for further investigation (Figure 7F). Collectively, these validation results suggest that most lncRNAs in the model exhibit distinct expression patterns, and additional samples may be required for comprehensive validation.

Figure 7 RT-qPCR results for relative mRNA expression in cuproptosis-related lncRNAs. Relative mRNA expression of (A) FAM13A-AS1, (B) GHRLOS, (C) LINC00456, (D) OPA1-AS1, (E) RAP2C-AS1, and (F) UBE2Q1-AS1. *P < 0.05, **P < 0.01.

Discussion

The high prevalence of BLCA poses a significant threat to human health and imposes a substantial economic burden. As the tumor progresses to MIBC, the specific survival rate declines to 65–78%.^5,16 Necessitating further exploration of the biological mechanisms underlying MIBC development, particularly at the genetic level.

Cuproptosis, a novel copper-dependent cell death modality,⁸ is regulated by copper ions binding to lipid-acylated components of the tricarboxylic acid cycle (TCA) in mitochondrial respiration, subsequently causing proteotoxic stress and cell death.^17,18 Numerous cuproptosis-related genes have been implicated in cancer progression.^19,20 LncRNAs have emerged as crucial regulators in tumorigenesis by modulating nuclear structure, transcription, mRNA stability, translation, and post-translational modifications.^21,22 It has been reported that cuproptosis-related genes can regulate the stability of lncRNA to promote tumor progression. Exploring the expression patterns of cuproptosis-related lncRNAs may provide us with new ideas and targets for diagnosis, treatment, and prognosis assessment in MIBC.

In this study, we investigated the expression and functions of cuproptosis-related lncRNAs in MIBC and BLCA patients using TCGA and IMvigor210 databases. LASSO Cox regression analysis identified six lncRNAs (FAM13A-AS1, GHRLOS, LINC00456, OPA1-AS1, RAP2C-AS1, and UBE2Q1-AS1) with independent prognostic values for MIBC. A risk model based on these lncRNAs classified MIBC patients into high- and low-risk groups, with the latter exhibiting significantly better survival outcomes. The risk score, an independent prognostic factor, demonstrated good predictive ability for MIBC patients’ clinical characteristics.

Cuproptosis-related lncRNAs may act as oncogenes or tumor suppressors. Our study identified 6 cuproptosis-related lncRNAs which had prognostic values. FAM13A-AS1 is an important tumor suppressor factor.^23,24 Previous research found that the upregulation of FAM13A-AS1 inhibits the proliferation, migration, and invasion of breast cancer cell lines by targeting the miRNA-205-3p/DDI2 axis.²⁵ GHRLOS expression levels in colorectal cancer tissues were significantly lower than in matched normal tissues, and decreased GHRLOS expression was significantly associated with the occurrence of lymph node metastasis and distant metastasis in colorectal cancer.²⁶ There are relatively few reports on UBE2Q1-AS1, and studies have shown that UBE2Q1-AS1 may contribute to GC development.²⁷ RAP2C-AS1 is widely expressed in tissues and organs. Studies have indicated that RAP2C-AS1 regulates the occurrence and development of esophageal cancer through the lncRNA-miRNA-mRNA axis.²⁸ Our RT-qPCR results confirmed the low expression of RAP2C-AS1 in tumor tissues, which may imply that it is a tumor suppressor gene.

To elucidate the role of cuproptosis-associated lncRNAs in MIBC, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses based on differentially expressed genes (DEGs) between distinct risk groups. Unexpectedly, we discovered that numerous cytokine activities, lymphocyte-mediated immune responses, and cytokine-cytokine receptor interactions were enriched. It is plausible to hypothesize that cuproptosis-related lncRNAs are intimately linked to tumor immune regulation, which aligns with prior findings.^29,30 Delving further into the impact of cuproptosis-associated lncRNAs on immune cell infiltration, we observed that various tumor-associated immune cells, including Treg cells, mast cells, and follicular helper T cells, were highly infiltrated in the high-risk group, alongside an enrichment of numerous anti-tumor immune-related signaling molecules. Intriguingly, we identified a substantial presence of immunosuppressive molecules, such as PD-1, IDO1, CTLA4, and PD-L1, in high-risk patients. Previous research has substantiated that elevated PD-L1 expression in MIBC patient tumor tissues enhances the efficacy of anti-PD-L1 therapy.³¹ I Our findings indicate that high-risk patients in our cohort may exhibit improved responses to anti-PD-L1 immunotherapy. We subsequently validated our hypothesis utilizing the IMvigor210 dataset of BLCA patients undergoing anti-PD-L1 treatment. The constructed predictive model also demonstrated strong prognostic accuracy for anti-PD-L1 treatment response; however, we were surprised to find that patient survival in the high-risk group did not improve following anti-PD-L1 therapy. This may be attributable to the elevated expression of immunosuppressive molecules within the high-risk group, rendering nanotherapeutic anti-immune treatments ineffective.³² These observations suggested that future therapeutic approaches may be directed toward combination therapies for the clinical benefit of more patients.^33,34 For patients with high risk factors, we may need to manage their treatment modalities more closely. Targeting six cuproptosis-related lncRNAs may become alternative treatment options for them in the future.

Nonetheless, this study is not without limitations, including the restricted sample size and the absence of methylation level analysis for the candidate lncRNAs. To comprehensively elucidate the regulatory mechanisms of candidate lncRNAs in MIBC, future investigations should encompass an expanded sample size and the acquisition of solid tumor specimens. Naturally, the regulator mechanism of the six Cuproptosis-related lncRNAs screened for MIBC are theoretical, and numerous of experiments are needed to verify their roles in future. Another major limitation is that genetic testing is more expensive and complex than such several readily available, easy measurable blood-based biomarkers, and the expression patterns of the target genes may vary widely due to differences in the assay platforms.

Conclusion

In summary, we identified six prognostic lncRNA markers associated with cuproptosis-related genes and validated their predictive significance in the context of MIBC patient prognosis by leveraging publicly available databases. Our findings suggest that candidate cuproptosis-related lncRNAs may modulate the expression of immunosuppressive molecules, thereby influencing tumor initiation and progression. This investigation has enhanced our comprehension of cuproptosis-related lncRNAs and their mechanistic roles within the tumor microenvironment, potentially unveiling novel therapeutic targets and prognostic biomarkers for further exploration.