Prognostic potential of m7G-associated lncRNA signature in predicting bladder cancer response to immunotherapy and chemotherapy

Introduction

Urothelial carcinoma is most commonly localized in bladder, both in China and on a global scale [1, 2]. Bladder cancer (BC) manifests with age-specific and gender-specific incidence disparities; in incidence, with notable gender differences as well [3]. Despite lower incidence rates among females, their prognostic outcomes are less favorable compared to males [3, 4]. For individuals diagnosed with muscle-invasive BC, the standard therapeutic intervention is radical cystectomy (RC). This procedure imposes considerable economic, social, and psychological burdens. Importantly, suboptimal five years survival rates persist even post-RC [5]. Given the range of existing therapeutic options, BC remains susceptible to recurrence, progression, and increased mortality rates, emphasizing the need for innovative therapeutic strategies [4]. In light of these clinical imperatives, efforts are actively searching for potent underway to identify effective pharmaceutical agents [6]. Concurrent research endeavors utilize various investigative strategies, including but not limited to single-cell RNA sequencing, artificial intelligence, and analyses of clinical and lifestyle parameters, lifestyle determinants, and beyond, to unravel the fundamental underlying mechanisms governing responsible for BC initiation and progression [6, 7]. This multidisciplinary approach has contributed to the development of an array of treatment modalities, including surgery, chemotherapy, immunotherapy, and radiotherapy. Of these, emerging immunotherapies represent a significant advancement. To optimize the utility of these treatments, there is an urgent requirement for the identification and validation of reliable prognostic biomarkers [8]. Such biomarkers would enable clinicians to formulate personalized treatment plans, thereby facilitating informed decision-making regarding aggressive interventions, such as RC.

RNA methylations play a critical role in the regulation of posttranscriptional gene expression [9]. The dysregulation of such methylations is notably implicated in disease onset and prognostic outcomes. Specifically, RNA methylations are integral to both the initiation and progression of various cancers [9]. While recent investigations have considerably enriched our comprehension of the link between mRNA N⁶-methyladenosine (m6A) and cancer, a significant knowledge gap persists concerning N⁷-methylguanosine (m7G) and cancer, particularly in long non-coding RNAs (lncRNAs) [10]. M7G is a key catalytic agent for methylation at the N7 position of ribo-guanosine, thereby augmenting the efficacy of mRNA translation efficiency [11, 12]. The role of m7G in cancer pathogenesis has elicited substantial scientific scrutiny, as numerous studies underscore its salient role in carcinoma development [13]. In the realm of prostate cancer, for example, the m7G-associated gene eIF4E has been identified to elevate the translation rates of oncogenic mRNAs, hence fostering tumorigenicity [14]. In the specific context of BC, our prior work introduced an m7G-related mRNA molecular subtype with predictive capacity for patient prognosis and therapeutic response [15]. Furthermore, METTL1, another m7G-associated gene, has been demonstrated to be overexpressed in tumor samples, thereby significantly accelerating BC progression [16]. However, the role of m7G-related lncRNAs in BC has been comparatively underexplored. LncRNAs critically modulate gene expression through a variety of mechanisms, including chromatin remodeling, transcriptional and post-transcriptional regulation, and protein metabolism [17, 18]. Preliminary studies have begun to examine the correlation between m7G and lncRNAs in other malignancies. For instance, in colon cancer, a signature based on m7G-related lncRNAs has been formulated to assess patient prognosis effectively [19]. Similar signatures have been constructed in endometrial cancer for predicting drug sensitivity, and in renal cancer for anticipating immunotherapeutic responses [20, 21].

Despite these advancements, a comprehensive understanding of the interplay between BC and m7G-associated lncRNAs remains to be achieved. Consequently, our research aim is to delineate the influence of m7G-related lncRNAs on BC through the establishment of a bioinformatically derived signature.

Materials and methods

Results

Constructing the risk score

The study workflow is outlined in Figure 1, and Supplementary Table 1 lists the 42 m7G-associated genes under investigation. Pearson correlation analysis yielded 271 m7G-related lncRNAs, as illustrated in Figure 2A. A Venn Diagram was subsequently utilized to filter 20 co-expressed, differentially regulated m7G-related lncRNAs. A lasso regression model was then applied to this subset, selecting seven lncRNAs (GATA3-AS1, LINC00930, LINC01341, MED14OS, MIR100HG, RUSC1-AS1, SNHG4) for the formulation of the m7G-related risk score, as detailed in Figure 2B. The prognostic value of these selected lncRNAs was evaluated via univariable Cox regression analysis. Figure 2C substantiates that GATA3-AS1, LINC00930, LINC01341, MED14OS, and MIR100HG serve as effective prognostic indicators of OS and BC patients within the TCGA dataset. Patients in both TCGA and IMvigor210 cohorts were stratified based on median risk score values into high and low-risk groups. Elevated risk scores were significantly correlated with poorer OS (Figure 2E, p<0.001) and CSS (Figure 2E, p<0.001) in the TCGA cohort. This trend was corroborated in the external IMvigor210 cohort, where a high-risk score was indicative of a worse prognosis (Figure 2F, p=0.022).

Figure 1:

The workflow of this study.

Figure 2:

Construction of the m7G-related risk score: the co-expressed differentially lncRNAs (A), the result of the LASSO regression model (B), the result of the univariable Cox regression model in TCGA (C). Validation of the risk score: the Kaplan–Meier method results of overall survival (D) and cancer-specific survival (E) in TCGA, and the Kaplan–Meier method result of overall survival in IMvigor210 (F). The correlation between the risk score and clinical parameters: TCGA: WHO grade (G), T stage (H), lymph node metastasis (I), age (J), sex (K); IMvigor210: sex (L).

Basic clinical parameters for both TCGA and IMvigor210 cohorts are tabulated in Table 1 and Supplementary Table 2, respectively. Within the TCGA dataset, Table 1 demonstrates a statistically significant relationship between the risk score and multiple clinical variables, including WHO grade, AJCC stage, T stage, lymph node metastasis, OS, and CSS. In the IMvigor210 cohort, the risk score was found to significantly discriminate differences in OS.

Table 1:

The clinicopathologic characteristics of the TCGA included patients.

Characteristic	High risk-score	Low risk-score	p
n	201	197
Age, mean±SD	68.44±9.92	67.16±11.08	0.227
Follow-up time (median, month)	48.67	27.73	<0.001
Survival time (median, month)	22.83	88.03	<0.001
BMI, mean±SD	27.32±5.71	26.61±5.61	0.244
Gender, n, %			0.914
Female	54 (13.6 %)	51 (12.8 %)
Male	147 (36.9 %)	146 (36.7 %)
Smoking history, n, %			0.333
No	49 (12.7 %)	56 (14.5 %)
Yes	148 (38.4 %)	132 (34.3 %)
Lymphatic vessel invasion, n, %			0.490
No	60 (22.2 %)	63 (23.3 %)
YES	79 (29.3 %)	68 (25.2 %)
WHO grade, n, %			<0.001
High grade	199 (50.4 %)	178 (45.1 %)
Low grade	1 (0.3 %)	17 (4.3 %)
AJCC stage, n, %			<0.001
AJCC stage I–II	38 (9.6 %)	88 (22.2 %)
AJCC stage III–IV	163 (41.2 %)	107 (27 %)
T stage, n, %			<0.001
T stage 3_4	154 (42 %)	95 (25.9 %)
T stage a_2	39 (10.6 %)	79 (21.5 %)
Lymph node metastasis, n, %			0.015
YES	78 (21.8 %)	49 (13.7 %)
NO	109 (30.5 %)	121 (33.9 %)
Distant metastasis, n, %			0.747
YES	5 (2.5 %)	5 (2.5 %)
NO	81 (40.1 %)	111 (55 %)
Overall survival, n, %			<0.001
Alive	88 (22.1 %)	135 (33.9 %)
Dead	113 (28.4 %)	62 (15.6 %)
Cancer-specific survival, n, %			<0.001
Alive	111 (28.9 %)	153 (39.8 %)
Dead	79 (20.6 %)	41 (10.7 %)

1. AJCC, American Joint Committe on Cancer; BMI, body mass index; SD, standard deviation; WHO, World Health Organization; n, number.

Validating the risk score

In the analysis of the TCGA dataset, we identified a statistically significant positive correlation between the risk score and pivotal clinical parameters, specifically WHO grade (Figure 2G, p=2.4e-06), T stage (Figure 2H, p=7.8e-07), and lymph node metastasis (Figure 2I, p=2.1e-04). Conversely, the data revealed no significant associations between the risk score and age (Figure 2J, p=0.15) or sex (Figure 2K, p=0.58). These observations were corroborated by the IMvigor210 dataset, which also exhibited an absence of significant correlation between the risk score and sex (Figure 2L, p=0.64).

To substantiate the prognostic utility of the risk score, Kaplan–Meier survival analyses were executed across multiple subgroups within the TCGA and IMvigor210 datasets. In the TCGA dataset, it was unequivocally observed that patients with elevated risk scores manifested significantly inferior OS in various demographic and clinical subgroups: females (Figure 3A, p=0.035), smokers (Figure 3B, p<0.001), individuals aged 70 years or below (Figure 3C, p=0.01), males (Figure 3D, p<0.001), patients with WHO high-grade tumors (Figure 3E, p<0.001), those devoid of lymphatic vessel invasion (Figure 3F, p<0.001), individuals at Ta_2 stage (Figure 3G, p=0.021), patients without lymph node metastasis (Figure 3H, p=0.002), those lacking distant metastasis (Figure 3I, p<0.001), those classified under AJCC stage III–IV (Figure 3J, p=0.026), and individuals with a body mass index of 25 or below (Figure 3K, p<0.001). Concomitant analysis of the IMvigor210 cohort corroborated the significance of the risk score, revealing a positive correlation with smoking status (Figure 3L, p=0.002), AJCC stage I–II (Figure 3M, p=0.015), and male gender (Figure 3N, p=0.032).

Figure 3:

The prognostic ability of the risk score: Kaplan–Meier analysis results of subgroups in TCGA: Female (A), smoker (B), age ≤70 (C), male (D), WHO high grade (E), no lymphatic vessel invasion (F), T a_2 (G), no lymph node metastasis (H), no distant metastasis (I), AJCC stage III–IV (J), body mass index ≤25 (K); IMvigor210: smoker (L), AJCC stage I–II (M), male (N).

Constructing and validating a nomogram

In the TCGA dataset, univariable Cox regression analysis results (Figure 4A) guided the development of a multivariable Cox regression model incorporating both the risk score and pertinent clinical variables such as age, smoking history, AJCC stage, T stage, lymph node metastasis stage, and distant metastasis stage. The multivariable analysis underscored the independent prognostic utility of the risk score in OS (Figure 4B, p<0.001). Similarly, in the IMvigor210 cohort, the risk score demonstrated significant prognostic implications in both univariable (Figure 4C, p=0.008) and multivariable (Figure 4D, p=0.015) Cox regression models.

Figure 4:

Construction and validation of a nomogram: the results of univariable (A) and multivariable (B) Cox regression model in TCGA, the results of univariable (C) and multivariable (D) Cox regression model in IMvigor210, the TCGA nomogram (E), the results of multiparameter ROC analysis without (F) and with (G) the risk score in TCGA nomogram, the calibration curves of TCGA nomogram (H), the DCA curve of TCGA nomogram (I).

A nomogram was constructed based on the multivariable Cox regression model from the TCGA dataset (Figure 4E). The C-index was calculated as 0.72 (0.689–0.754). For further nomogram validation, a multiparameter ROC analysis was conducted, revealing that the inclusion of the risk score elevated the area under the curve value from 0.73 (Figure 4F) to 0.8 (Figure 4G). Calibration curves for 1, 3, and 5 years exhibited strong alignment between predicted and observed outcomes (Figure 4H). Additionally, decision curve analysis (DCA) confirmed the favorable performance of the risk score (Figure 4I). Supplementary Figure 1 provides the IMvigor210 nomogram and corresponding validation data.

Enrichment analyses for the risk score

Utilizing TCGA data, we conducted GO, KEGG, and GSEA to elucidate functional insights. GO analysis indicated significant enrichment in biological processes (BP) encompassing neutrophil-mediated immunity, granulocyte chemotaxis, and neutrophil chemotaxis, as well as in cellular components including the external side of the plasma membrane and the cornified envelope. Molecular function (MF) was enriched in categories such as chemokine activity, immunoglobulin binding, and IgG binding (Figure 5A). Similarly, KEGG analysis highlighted immune-related pathways, including cytokine–cytokine receptor interaction, the PI3K-Akt signaling pathway, the B cell receptor signaling pathway, the chemokine signaling pathway, and the IL-17 signaling pathway (Figure 5B). To enhance the credibility of our GSEA results, we enriched pathways from both REACTOME and KEGG databases. As depicted in Figure 5C, increased REACTOME pathways were primarily related to immune processes, cytokine signaling in the immune system, adaptive immune system, signaling by interleukins, and neutrophil degranulation pathways. In contrast, downregulated REACTOME pathways were primarily involved in metabolic processes, encompassing cytochrome P450 – organized by substrate type, fatty acid metabolism, regulation of lipid metabolism via PPAR alpha, and lipid metabolism (Figure 5D). Similarly, enhanced KEGG pathways were predominantly immune-related (Figure 5E), while downregulated KEGG pathways were primarily associated with metabolism, notably drug metabolism (Figure 5F).

Figure 5:

Function analysis: Gene ontology results (A), Kyoto Encyclopedia of Genes and Genomes results (B), Gene Set Enrichment Analysis results in REACTOME (C and D) and KEGG (E and F).

The risk score in immune checkpoints and immune infiltration analyses

Regarding immune checkpoints, samples with higher risk scores exhibited positive correlations with the expression levels of CD274, CD40, CTLA4, PDCD1, and PDCD1LG2, among others. Inversely, the risk score maintained a positive association with the expression of CD96, SIGLEC15, TNFRSF14, TNFRSF25, and TNFSF15 (Figure 6A). Regarding immune cell infiltration, high-risk score samples exhibited positive correlations with CD4+ T memory activated cells, T gamma delta cells, Macrophage M0 cells, Macrophage M1 cells, Macrophage M2 cells, and neutrophil cells. Conversely, these high-risk samples displayed negative associations with B plasma cells, NK resting cells, and monocytes (Figure 6B).

Figure 6:

Immune checkpoints analysis results (A), immune infiltration (B). The TCIA score result (C). Chemotherapy response (D). IC₅₀: The half-maximal inhibitory concentration; ns, p≥0.05; *, p<0.05; **, p<0.01; ***, p<0.001.

Discussion

The m7G cap serves as an indispensable element for post-transcriptional modifications and is a key factor for efficient RNA processing [17]. Numerous m7G-associated lncRNAs have been implicated as prognostic indicators across a range of malignancies [36]. Moreover, m7G-related mRNA has proven to be of prognostic relevance in breast cancer patients [20]. Building on this foundation, the current study endeavors to delineate the relationship between m7G-related lncRNAs and BC. We have formulated a BC-specific signature, and our empirical data corroborate the prognostic utility of this m7G-related lncRNA signature in BC cohorts, a validation that extends to both internal and external samples. Additionally, we discerned a meaningful association between this signature and the likely effectiveness of both immunotherapeutic and chemotherapeutic interventions, thereby endorsing the prospective utility of these lncRNAs and the signature as valuable biomarkers in BC.

GATA3-AS1 is observed to be markedly upregulated in breast cancer and is implicated in both tumor progression and mechanisms of immune evasion [37]. It also has a contributory role in the carcinogenesis of endometrial carcinoma [38]. LINC00930 is upregulated and exhibits a relationship with tumorigenesis in nasopharyngeal carcinoma [39]. LINC01341 has been characterized as a 5-methylcytosine-associated lncRNA in lung squamous cell carcinoma, thereby offering insights into its plausible involvement in RNA methylation processes [40]. MIR100HG is elevated in BC and plays a role in disease progression [41], while concurrently influencing the formation of cancer stem cells and affecting the distant metastasis of lung cancer [42]. RUSC1-AS1 accelerates the pathogenesis of both osteosarcoma [43] and breast cancer [44, 45]. Finally, the expression of SNHG4 has been found to impact multiple cancer types [46], [47], [48], [49].

A formidable obstacle in harnessing lncRNA-related signatures for BC research lies in the paucity of external validation, largely due to limited lncRNA expression data in publicly available datasets. To address this, our study employed external validation using the IMvigor210 dataset, which incorporates lncRNA expression metrics, thereby fortifying the credibility of our conclusions. We found that the risk score not only demonstrated statistically significant variations across an array of clinical parameters but also sustained positive correlations with multiple clinical variables. Notably, no significant associations were discerned between the risk score and either age or smoking history across both the TCGA and IMvigor210 datasets. Within the TCGA dataset, the risk score was relatively consistent between female and male cohorts. Therefore, the applicability of the risk score appears to be expansive, unhindered by variables such as age, smoking history, or gender. In prognostic assessments, the risk score independently evinced substantial prognostic merit for BC patients in the TCGA dataset, a finding corroborated in the IMvigor210 cohort. Age is a salient risk factor for BC [50] and influences the effectiveness of BCG treatment, thereby informing therapeutic decisions [51, 52]. Fortunately, the risk score exhibited no significant differences between patients aged 70 years or less and those aged 70 years or older, implying that the risk score can be applied without age-related restrictions. Importantly, we observed no marked differences in the risk score between patients aged 70 years or younger and those older than 70 years, suggesting the score’s applicability across age groups. Furthermore, the risk score demonstrated prognostic significance across heterogeneous subgroups in both the TCGA and IMvigor210 datasets. This corroborates earlier studies highlighting the prognostic utility of m7G-related lncRNAs in a variety of malignancies, including lung adenocarcinoma [53], colon carcinoma [54], hepatocellular carcinoma [36], and endometrial cancer [20]. In urological malignant tumors, the m7G-related lncRNA risk score has also independently prognosticated patient outcomes in renal cancer [21]. Based on these findings, we postulate that our risk score could function as a viable biomarker for BC.

GO analysis revealed enrichment in neutrophil-mediated immunity, neutrophil chemotaxis, and chemokine activity. Concurrently, KEGG and GSEA identified enriched pathways pertinent to immune functions and metabolic processes, including cytokine–cytokine receptor interaction, B cell receptor signaling pathway, adaptive immune system, and neutrophil degranulation, among others. In the ensuing checkpoint analysis, a positive correlation was observed between the risk score and the expression levels of key immune markers such as CD274, CD40, CTLA4, PDCD1, and PDCD1LG2. In contrast, samples characterized by low-risk scores demonstrated a positive correlation with the expression levels of CD96, SIGLEC15, TNFRSF14, TNFRSF25, and TNFSF15. The current literature presents no consensus on the relationship between PD-L1 expression and response to immunotherapy [55]. However, findings from the phase 3 IMvigor130 trial indicated that patients with elevated PD-L1 expression experienced favorable outcomes with atezolizumab (anti-PD1/PD-L1) monotherapy [56]. These data corroborate the notion that higher risk scores may predict increased efficacy from anti-PD1/PD-L1 immunotherapy. M6A RNA modification has been implicated as a pivotal factor influencing the tumor immune microenvironment in BC [57, 58]. Our analysis of immune cell infiltration in TCGA samples, approached from the perspective of m7G modifications, revealed a positive association between risk scores and increased infiltration of CD4+ T memory activated cells, M2 macrophages, and neutrophils. Elevated infiltration of these specific immune cell types generally indicates a higher likelihood of positive response to immunotherapy [58], [59], [60]. In alignment with our prognostic outcomes, the presence of Macrophage M0 cells was associated with an unfavorable prognosis [61]. This result lends further credence to the hypothesis that individuals with higher risk scores are more likely to derive pronounced benefits from anti-PD1/PD-L1 immunotherapy. Finally, we employed data from the Tumor Cancer Immunome Atlas (TCIA) to project immunotherapy responses. The analysis further substantiated that high-risk score patients are more likely to demonstrate favorable responses to anti-PD1 immunotherapy.

The empirical data suggest that patients categorized with higher risk scores may exhibit an enhanced responsiveness to anti-PD1 immunotherapy. This supposition is congruent with the upregulated expression of PD-1 and PD-L1 observed in this high-risk cohort. Conjoint analyses of gene expression patterns, immune cell infiltration metrics, and data from the TCIA collectively fortify the proposition that patients with elevated risk scores are more likely to benefit from anti-PD1/PD-L1 immunotherapeutic interventions as opposed to those with lower scores. Moreover, extant research lends credence to the idea that patients with lower risk scores could potentially derive therapeutic advantages from anti-CTLA4 immunotherapy. Investigations led by Wang et al. demonstrate that the lncRNA GATA3-AS1 modulates immune evasion mechanisms in breast cancer via the regulation of GATA3 and PD-L1 expression [62]. Similarly, aberrant expression of MIR100HG has been implicated in immune escape phenomena in the context of gastric cancer [63]. The regulatory role of SNHG4 in immune evasion in colorectal cancer is enacted through its modulation of the miR-144-3p/MET signaling axis [64]. These molecular mechanisms serve to elucidate the utility of the risk score as a predictive metric for immunotherapy responsiveness. In pharmacological responsiveness, patients with elevated risk scores displayed heightened sensitivity to a panel of chemotherapeutic agents, including cisplatin, docetaxel, doxorubicin, gemcitabine, and vinblastine. Conversely, methotrexate was identified as a chemotherapeutic agent to which patients with lower risk scores exhibited sensitivity. Notably, GATA3-AS1 has been ascertained to possess predictive value in gauging the response of breast cancer to neoadjuvant chemotherapy [37]. MIR100HG, which is derived from cancer stem cells, has been substantially implicated in chemotherapeutic resistance mechanisms [42, 65].

Although the expression and prognostic value of the risk score were exhibited in TCGA and external cohorts. We still need to demonstrate these findings using basic experiments in our future work.

Conclusions

Our results validated the risk score had significant clinical value in BC patients. Additionally, a notable correlation was observed between the risk score and responses to treatments. The included m7G-related lncRNAs and the signature show promising potential for managing BC.