FGFR Alterations in NMIBC and MIBC

The cohort in the present study included 124 non-muscle invasive bladder cancer (NMIBC: 5-year overall survival (OS) rate: 83%) and 265 muscle-invasive bladder cancer (MIBC: 5-year OS rate: 35%) with a median follow-up of 36 and 22 months, respectively (Supplementary Fig. 1A). We first assessed mRNA expression levels of FGFR family genes in the cohort. FGFR3 is actively transcribed among the FGFR families in the bladder, especially in BLCA (Supplementary Fig. 1B). FGFR3 was the most frequently mutated gene (81/389; 21%) among the FGFR families including FGFR1 (19/389; 5%), FGFR2 (8/389; 2%), and FGFR4 (12/389; 3%) (Supplementary Fig. 1C). A previous TCGA publication defined 58 significantly mutated genes (SMGs) in BLCA [9], which was largely recapitulated in our cohort including the clinicopathological characteristics (Fig. 1A, B).

Fig. 1

FGFR Alterations in NMIBC and MIBC. A Mutation landscape of 58 significantly mutated genes defined by the TCGA publication [9] in 389 bladder cancer (BLCA) samples from the OMPU-NCC cohort. The patients were classified into pTa (n = 59), pT1 (n = 65), and ≥ pT2 (n = 265, MIBC: muscle-invasive bladder cancer). B Recurrent mutation rate of 58 significantly mutated genes according to pathological T stages. C Schematic of the FGFR3 fusions identified in our cohort. FGFR3-TACC3 fusions were found in 14 of 289 patients, and the most frequent pattern (7 of 11) is shown. NSD2 and SPON2 are newly identified fusion partners. D FGFR3 mRNA expression levels according to the FGFR3 alterations. The difference was assessed by the Mann–Whitney U test; p < 0.05*, p < 0.001**, p < 0.0001***. E, F Kaplan–Meier curves demonstrating progression-free survival (PFS) in non-muscle-invasive bladder cancer (NMIBC) (E) and overall survival (OS) in MIBC (F). A log-rank test was used to assess the survival difference between the two groups; p < 0.05*

To delineate the allelic difference in FGFR3 among ethnicity, we sought to assess whether there is a specific allelic variant in the germline for the Asian population. GnomADv3.0, an integrative germline dataset of 71,702 individuals (mostly Western population), was utilized for the control [10]. We referred the Asian germline dataset (jMorp-14KJPN) [11] and identified five significantly enriched non-synonymous single nucleotide polymorphisms (SNPs) on the FGFR3 gene locus that are specific to the Asian population (Q29H, G65R, L164V, T450M, and A720S) (Supplementary Fig. 1D). However, these SNPs were not enriched in BLCA samples (Supplementary Table 1). Compared to BLCA with iFGFR3, FGFR3 mRNA expression level was significantly upregulated in patients with recurrent FGFR3 mutations but not in patients with the SNPs (Supplementary Fig. 1E). There seemed to be no survival difference based on the FGFR3 status (Supplementary Fig. 1F), indicating no clinical implication of these Asian-specific SNPs in FGFR3.

FGFR3 mutations was predominantly observed in cases with lower malignant properties such as NMIBC (pTa: 51%, pT1: 29%, more than pT2: 12%), low grade and negative lymph vascular invasions (Table 1). Regarding the mutational alleles, the TCGA publication, which only consists of MIBC samples, reports S249C and Y373C as the top two frequent FGFR3 mutations in BLCA (Supplementary Fig. 2A). We noted that recurrent K650E and T757P nonsynonymous mutations at the kinase domain (KD) were frequently observed in MIBC in our cohort compared to the TCGA cohort (Supplementary Fig. 2B). Although we examined the prognosis of five MIBC cases with mutation at KD (three in K650E and two in T757P), there was no difference in OS compared to that in cases with other FGFR3 alterations (Supplementary Fig. 2C). Interestingly, we found that FGFR3 mutations at KD were more prevalent in MIBC than in NMIBC cases (p = 0.021) (Supplementary Fig. 2D).

Table 1 Clinicopathological characteristics in 389 BLCA patients according to the FGFR3 status at the collection of biospecimens

FGFR3 mRNA expression levels were consistently upregulated in aFGFR3 compared to iFGFR3, regardless of the mutation sites (Supplementary Fig. 2E). The present study exhibited a frequency of 4% (17/389) for FGFR3 fusions (Table 1), including novel fusion partners (NSD2 and SPON2) (Fig. 1C). No histological variant was observed in cases with FGFR3 fusions (Table 1). The KD located at the C-terminus of FGFR3 has been retained in 13 of 17 (77%) fusions, and the mRNA expression level was significantly upregulated in cases with FGFR3 fusions compared to cases with iFGFR3 (Supplementary Fig. 2F). Upon stratifying our cohort into NMIBC and MIBC categories, FGFR3 mRNA expression levels were significantly higher in aFGFR3 cases than in iFGFR3 cases in both NMIBC and MIBC (Fig. 1D), with the highest median mRNA expression levels observed in MIBC patients with FGFR3 fusions. This finding underscores the clinical importance of detecting FGFR3 fusions, alongside mutations, in advanced MIBC patients and accentuates the importance of considering recently approved FGFR3 inhibitors [12]. The FGFR3 protein expression levels were increased in aFGFR3 compared to iFGFR3 cases (Supplementary Fig. 2G, H). We investigated the progression free survival (PFS) of 124 NMIBC patients (Fig. 1E). Patients with recurrent FGFR3 mutations showed a significantly better PFS compared to those with iFGFR3 (p = 0.037). However, this distinction was not evident in patients with FGFR3 fusions (p = 0.806). In the context of OS among 265 MIBC patients, no significant differences were observed based on FGFR3 status (Fig. 1F).

The Association between FGFR3 alteration and molecular Subtypes

We have adopted the established consensus MIBC subtype [6], the UROMOL subtype for NMIBC [13], and Baylor college subtype [14] (Fig. 2A and Supplementary Fig. 3A). As expected, distribution of subtypes significantly differed between NMIBC and MIBC across these three subtyping systems (Consensus MIBC: p < 0.0001, UROMOL: p < 0.0001, and Baylor: p < 0.0001) (Supplementary Fig. 3B-3D). In the overall cohort (n = 389), FGFR3 alterations were enriched in class_1 (54%) and class_3 (94%) for the UROMOL subtype and LumP (42%) for consensus MIBC subtype (Fig. 2B and Supplementary Fig. 3E). An elevated FGFR3 mRNA expression level was confirmed within these molecular subtypes (Fig. 2C, D and Supplementary Fig. 3F).

Fig. 2

Association between FGFR3 Alteration and Molecular Subtypes. A Summary of FGFR3 alterations, FGFR3 mRNA expression, histological variant, consensus MIBC subtypes [6], UROMOL NMIBC subtypes [13], and Baylor college [14]. B Pie chart of FGFR3 alterations including recurrent mutations and fusions in 389 BLCA cases. C FGFR3 mRNA expression levels according to UROMOL subtypes. D FGFR3 mRNA expression levels according to consensus MIBC subtypes. E Estimated proportion of UROMOL subtypes in 124 NMIBC cases. F Estimated proportion of consensus MIBC subtypes in 265 MIBC cases. G Estimated proportion of consensus MIBC subtypes in 265 MIBC cases categorized based on FGFR3 mutational status (KD: kinase domain; note that seven cases with FGFR3 fusions were classified into mut- group)

We further stratified the cohort into NMIBC (n = 124) and MIBC (n = 265) cohorts. In NMIBC, the UROMOL subtype showed a significant increase in class_1 and class_3 subtypes among aFGFR3 cases (p < 0.0001) (Fig. 2E). In MIBC, the consensus MIBC subtypes exhibited a significantly higher prevalence of the LumP subtype in aFGFR3 cases (49% vs 28%) (Fig. 2F). Interestingly, the proportion of Ba/Sq subtype in MIBC cases was similar between iFGFR3 and aFGFR3 (39% vs 38%). Other subtypes, such as LumU, appeared to decrease in aFGFR3 cases in place of an increase in LumP. Considering the enrichment of FGFR3 mutations at KD in MIBC cases (Supplementary Fig. 2D), we further assessed their association with molecular subtypes. Mutations at KD were notably inclined towards the Ba/Sq subtype (seven of eight: 88%), contrasting with other mutations (24 cases: 67% in LumP and 29% in Ba/Sq subtypes) (Fig. 2G). These findings illuminate the intricate interplay between aFGFR3 and BLCA subtypes.

Different Pathways Modified by aFGFR3 between NMIBC and MIBC

Principal component analysis of the whole transcriptome exhibited the delineation of aFGFR3 (Fig. 3A). Thus, we examined gene set enrichment analysis (GSEA) for all human collections (H, C1-C8: 23,734 gene sets) (Supplementary Fig. 3G). There was no gene set showing a false discovery rate (FDR) of < 0.25 in C7 (immunologic_signature_gene_sets), whereas H (hallmark_gene_sets) represented epithelial-mesenchymal transition (EMT) pathway as the most down-regulated pathways in aFGFR3 (Supplementary Fig. 3H, Supplementary Table 2).

Fig. 3

Different Pathways Modified by aFGFR3 between NMIBC and MIBC. A Principal component analysis for whole transcriptome data in iFGFR3 and aFGFR3. B Gene set enrichment analysis plotting all human MSigDB collections (Hallmark, C1-8: 23,734 gene sets) by false discovery rate q-value (FDR-q) and normalized enrichment score (NES). The analysis was performed separately in NMIBC (left panel) and MIBC (right panel). C Top 10 gene sets upregulated in MIBC/aFGFR3 (n = 39) compared to MIBC/iFGFR3 (n = 226). D mRNA expression levels in putative basal markers (KRT5, 15, 6A, and 16) and luminal markers (FOXA1 and KRT19). The difference in the expression level was assessed by the Mann–Whitney U test

To further delve into the insight of affected pathways in aFGFR3, we re-ran GSEA after separating the cohort into NMIBC (n = 124) and MIBC (n = 265) groups. Interestingly, the pathways exhibiting significant difference between iFGFR3 and aFGFR3 cases (FDR_q value < 0.25) were notably distinct in NMIBC and MIBC (Fig. 3B). In NMIBC (69: iFGFR3 vs 55: aFGFR3), only one (1/23,734 gene sets) pathway (C2: LINDGREN_BLADDER_CANCER_CLUSTER_3_DN) was significantly upregulated in aFGFR3 cases, while 3,630 pathways (C1: 6, C2: 940, C3: 0, C4: 133, C5: 1372, C6: 120, C7: 596, C8: 457, H: 6) were significantly downregulated in aFGFR3 cases. In stark contrast, MIBC (226: iFGFR3 vs 39: aFGFR3) exhibited significant elevation in 140 pathways (C1: 2, C2: 45, C3: 0, C4: 46, C5: 3, C6: 0, C7: 0, C8: 44, H: 0) and downregulation in 95 pathways (C1: 0, C2: 0, C3: 0, C4: 44, C5: 0, C6: 5, C7: 0, C8: 46, H: 0) in aFGFR3 cases (Supplementary Table 3). Importantly, we discovered that the top 10 upregulated pathways in MIBC with aFGFR3 include pathways associated with “basal/squamous epithelial” characteristics (Fig. 3C). This observation was further substantiated by the significant increase in known basal markers such as KRT5, KRT15, KRT6A, and KRT16, as well as typical luminal markers like FOXA1 and KRT19 (Fig. 3D). These findings collectively suggest the substantial variation in pathways influenced by aFGFR3 between NMIBC and MIBC. Moreover, in addition to the previous findings indicating that aFGFR3 is associated with luminal subtypes [6], our results reveal that aFGFR3 in MIBC can influence both epithelial subtypes, including the basal/squamous type.

Immune checkpoint genes and immune-related cell compositions in aFGFR3

We next examined TME according to the FGFR3 status. First, PD-L1 expression as evaluated by the combined positive score (CPS) was positively correlated with CD274 mRNA expression (Supplementary Fig. 4A-B), and the CPS seemed to be lower in aFGFR3 than in iFGFR3 (p = 0.07) (Supplementary Fig. 4C-D). Similar to the previous studies, CD8 + T-cell counts in the specimens seemed lower in aFGFR3 (p = 0.088) (Supplementary Fig. 4E-F). We explored the correlation between FGFR3 status and expression levels of putative immune checkpoint genes in 389 BLCA (Supplementary Fig. 4G, Supplementary Table 4). There was a positive correlation among immune checkpoint genes, whereas the correlation of these genes with FGFR3 was modest (Fig. 4A). We analyzed the potential candidates that are differentially expressed between iFGFR3 and aFGFR3 by analyzing the expression of the immune checkpoint genes, and identified that T-cell exhaustion markers, including TIM3, were most upregulated in iFGFR3, whereas HVEM and dendritic cell marker such as CD40 were increased in aFGFR3 (Fig. 4B, Supplementary Fig. 4H).

Fig. 4

Immune Checkpoint Genes and Immune-Related Cell Compositions in aFGFR3. A Spearman rank correlation among FGFR3 status, tumor mutation burden (TMB), and immune-checkpoint genes in 389 BLCA. B Bubble plots of difference in the immune checkpoint genes expression levels according to FGFR3 status. C CIBERSORTx analysis estimating composition of immune-related cells in 389 BLCA samples [15]. D The estimated proportion of each immune cell type from CIBERSORTx in iFGFR3 and aFGFR3 (Mann–Whitney U test; p < 0.05*, p < 0.001**, p < 0.0001***, n.s: non-significant). E The difference in correlation coefficient among estimated proportions of immune-related cells between NMIBC/iFGFR3 (n = 69) and NMIBC/aFGFR3 (n = 55). F The difference in correlation coefficient among estimated proportions of immune-related cells between MIBC/iFGFR3 (n = 226) and MIBC/aFGFR3 (n = 39)

We next estimated immune cell composition by using CIBERSORTx, a digital cytometry from bulk tissues [15] (Supplementary Table 5), and revealed distinct infiltration patterns of various immune cell types in BLCA according to FGFR3 status, including B naïve cells, B memory cells, T-CD4 + memory resting cells, T-follicular helper cells, M0 macrophages, dendritic cell proportions (Fig. 4C-D). We also confirmed that the actual cell count (CD8 and FOXP3) in the specimens and the estimated proportion was significantly correlated (Supplementary Fig. 5A).

Since the comprehensive GSEA showed differential pathways influenced by aFGFR3 between NMIBC and MIBC (Fig. 3B), we compared immune cell composition based on FGFR3 status within NMIBC (n = 124) and MIBC (n = 265) (Supplementary Fig. 5B-C). The correlation between immune cell proportions varied between iFGFR3 and aFGFR3 in NMIBC and MIBC cases (Supplementary Fig. 5D-E). Specifically, in NMIBC, we observed a negative correlation between the estimated populations of M0 macrophages and activated dendritic cells in iFGFR3 (r = -0.2844, 95% CI: -0.4883 to -0.05121), while a positive association between these cells was evident in aFGFR3 (r = 0.3525, 95% CI: 0.09615 to 0.5649) (Fig. 4E). In MIBC, a negative correlation between the estimated populations of CD8 T cells and memory B cells was found in iFGFR3 (r = -0.1180, 95% CI: -0.2447 to -0.01272), whereas a positive association between these cells was observed in aFGFR3 (r = 0.3805, 95% CI: 0.07380 to 0.6214) (Fig. 4F). These data collectively indicate differential immune cell compositions between iFGFR3 and aFGFR3, highlighting the distinct biological impact of aFGFR3 on the TME, which further varies between NMIBC and MIBC.

TME Heterogeneity in aFGFR3

To further explore the TME differences according to FGFR3 status, we employed the EcoTyper (Fig. 5A) [16]. EcoTyper is a machine learning pipeline for identifying cell states from bulk expression data, which covers 12 major cell lineages, including immune-related cells, fibroblasts, endothelial cells, and epithelial cells (Supplementary Fig. 5F). This methodology delineated 69 transcriptionally distinct cell states, unveiling ten clinically relevant multicellular communities known as Cellular Ecotypes (CE1-10). These ecotypes showed significant correlations with OS (shorter in CE1 and longer in CE10) and response to CPIs across various cancer datasets.

Fig. 5

TME Heterogeneity in aFGFR3. A EcoTyper analysis identifying cell states from bulk expression data, which covers 12 major cell lineages, including immune-related cells, fibroblasts, endothelial cells, and epithelial cells [16]. A total of 69 transcriptional cell states were identified, and an analysis on the cell-state co-occurrence patterns offered ten clinically distinct multicellular communities known as Cellular Ecotypes (CE1-10). CEs had been shown to correlate with OS (shorter in CE1 and longer in CE10) and treatment response to checkpoint inhibitors (CPIs) in various types of cancer data sets. Note that In CE2 and CE8, different cell states were depicted in specific cell types, including fibroblasts, epithelial cells, and mono/macrophages. B Pie charts showing the proportion of patients with aFGFR3 in each CE. C Kaplan–Meier curves of OS according to the cell ecotypes (CE1-10). D Estimated proportion of NMIBC/aFGFR3 (n = 55) and MIBC/aFGFR3 (n = 39) among CE2, CE7, CE8, and other CEs. E Distribution of CE2, CE7, CE8, and other CEs in NMIBC/aFGFR3 (n = 55) and MIBC/aFGFR3 (n = 39). F Cell state for the fibroblasts was defined by EcoTyper in 264/389 (68%) of the present cohort. Among them, we compared the mRNA expression of 806 genes determining the state of fibroblasts in CE2 and CE8 for MIBC/aFGFR3 patients. G Cell state for the monocytes/macrophages was defined by EcoTyper in 232/389 (60%) of the present cohort. Among them, we compared the mRNA expression of 452 genes determining the state of monocytes/macrophages in CE2 and CE8 for MIBC/aFGFR3 patients

In our BLCA cohort, we noticed an enrichment of aFGFR3 (mutations/fusions) in specific CEs (CE2, CE7, and CE8) (Fig. 5B). Kaplan–Meier curves for each CE exhibited distinct OS favoring CE7 and CE8, whereas CE2 showed the worst OS with a median of 13 months (Fig. 5C). Upon further analysis of the NMIBC (n = 124) and MIBC (n = 265) cases, we found that CE2 predominantly comprised MIBC patients, whereas CE7 and CE8 were associated with NMIBC patients in aFGFR3 (n = 94) (Fig. 5D).

To explore the association between the TME in aFGFR3 and the response to CPIs, we focused on MIBC patients within CE2 and CE8. Of 39 MIBC/aFGFR3 cases, 11 (28%) and 13 (33%) were defined as CE2 and CE8, respectively (Fig. 5E). In these two CEs, different cell states were depicted in specific cell types, including fibroblasts, epithelial cells, and mono/macrophages (Fig. 5A). By extracting the expression matrix from the EcoTyper pipeline, we identified unique expression patterns in fibroblasts 806 genes (Fig. 5F) and mono/macrophages 452 genes (Fig. 5G) between CE2/aFGFR3 and CE8/aFGFR3 cases. These data collectively suggest the remarkable TME heterogeneity within the aFGFR3 subgroup, potentially influencing the efficacy of CPIs.

Response of CPIs according to molecular subtypes and FGFR3 status

In the present cohort, 72 of 389 patients were treated with CPIs (pembrolizumab: 60 patients and avelumab: 12 patients) (Fig. 6A). The objective response rate (ORR) was 22% (pembrolizumab: 20% and avelumab 33%) (Supplementary Fig. 6A). We first assessed the ORR according to CEs from EcoTyper (Supplementary Fig. 6B). In line with the original report [16], a favorable ORR in CE10 (40% in 10 cases) and a poor ORR in CE6 (0% in 8 cases) were confirmed. The ORR in the top two allocated ecotypes were 21% in CE2 (n = 14) and 29% in CE8 (n = 14). Patients achieving CR/PR exhibited significantly higher PD-L1 CPS (p = 0.013) and TIM3 positive cells (p = 0.004) than those with SD/PD/unknown response (Fig. 6B). Tumor mutation burden (TMB) was also positively correlated with PD-L1 CPS (p = 0.021) and TIM3 positive cells (p = 0.009) (Supplementary Fig. 6C).

Fig. 6

Response of CPIs According to Molecular Subtypes and FGFR3 Status. A Of 389 BLCA patients, 72 were treated with CPIs including PD-1 inhibitor pembrolizumab (n = 60) and PD-L1 inhibitor avelumab (n = 12). Oncoprint sorted by the treatment response in those 72 patients is shown. B PD-L1 combined positive score (CPS), Cell count of TIM3 + cells in high power field, and FGFR3 mRNA expression were compared according to the response to CPIs. The difference was assessed by the Mann–Whitney U test. C The Objective response rate (ORR) in 72 patients treated with CPIs according to consensus subtypes (left panel) and FGFR3 alterations (right panel). D The estimated proportion of consensus subtypes in the IMvigor210 trial (PD-L1 inhibitor atezolizumab in patients with metastatic urothelial carcinoma) [17]. Pie charts show the proportion of FGFR3 mutations (not included for FGFR3 fusions) among the subtypes. E The ORR in 274 patients treated with atezolizumab according to consensus subtypes (left panel) and FGFR3 mutations (right panel). F The ORR in the IMvigor210 trial (n = 274) and the present cohort (OMPU: n = 72) treated with CPIs in Ba/Sq subtype (left panel) and LumP subtype (right panel). Fisher’s exact test was performed to assess the difference of the ORR according to FGFR3 status. Note that the data from IMvigor210 does not include the information on FGFR3 fusions. G Differentially expressed gene (DEG) analysis between MIBC/LumP/iFGFR3 (n = 64) and MIBC/LumP/aFGFR3 (n = 19) in the present cohort

Regarding molecular subtypes, the ORR was 23%, 14%, 42%, and 20% in Ba/Sq, LumP, LumU, and other subtypes, respectively (p = 0.331) (Fig. 6C). The ORR to CPIs was comparable in aFGFR3 compared to iFGFR3 (31% vs 20%; p = 0.467). To validate our result, we analyzed the data from the IMvigor 210 trial [17], exploring the PD-L1 inhibitor atezolizumab in patients with metastatic urothelial carcinoma (UC) using RNA-seq and hybrid capture-based next-generation sequencing for 274 patients. Despite lacking information on FGFR3 fusions, we assessed the distribution of consensus MIBC molecular subtypes and FGFR3 mutations (Fig. 6D). Of 49 patients with FGFR3 mutations, 11 (23%), 32 (65%), 1 (2%), and 5 (10%) were assigned to Ba/Sq, LumP, LumU, and other subtypes, respectively, exhibiting a higher LumP proportion in aFGFR3 than in iFGFR3 (p < 0.0001) (Supplementary Fig. 6D). LumU was rarely observed in aFGFR3, which was consistent with our cohort (Fig. 2F and Fig. 6D). The ORR in the IMvigor210 trial was 20%, 20%, 35%, and 19% in Ba/Sq, LumP, LumU, and other subtypes, respectively (p = 0.20) (Fig. 6E). Notably, the ORR in patients with FGFR3 mutations (n = 49) was 25%, which was comparable to 21% in iFGFR3 cases (some of whom presumably harbored FGFR3 fusions).

Our study revealed the remarkable heterogeneity within the TME even among the aFGFR3 cases. Thus, we stratified the cohort of 72 patients treated with CPIs based on molecular subtypes and FGFR3 status (Fig. 6F and Supplementary Fig. 6E). In the Ba/Sq subtype, the ORR ranged from 17 to 27% regardless of FGFR3 alterations, with no significant difference. However, the LumP subtype presented a striking contrast: a significantly higher ORR in aFGFR3 (50%) cases compared to a 5% in iFGFR3 cases in the present cohort (p = 0.022). This trend was also observed in the IMvigor 210 cohort with an ORR of 25% in FGFR3 mutations and 12% in iFGFR3. We performed gene expression analysis between LumP/iFGFR3 (n = 64) and LumP/aFGFR3 (n = 19) in the present cohort (Fig. 6G). Strikingly, several immune-related genes were significantly upregulated in LumP/iFGFR3 cases, including IDO1, CCL24, IL1RL1, LGALS4, and NCAM (CD56). These findings underscore the potential of these genes as promising targets for immunotherapy.