Characterization of the tumor-infiltrating immune repertoire in muscle invasive bladder cancer

Introduction: Muscle-invasive bladder cancer (MIBC) is a heterogeneous disease with several taxonomic molecular subtypes showing different genetic, clinical, and epidemiological profiles. It has been suggested that MIBC-subtypes follow different tumorigenesis pathways playing decisive roles at different stages of tumor development, resulting in distinct tumor microenvironment containing both innate and adaptive immune cells (T and B lymphocytes). We aim to characterize the MIBC tumor microenvironment by analyzing the tumor-infiltrating B and T cell repertoire according to the taxonomic molecular subtypes.

Methods: RNAseq data from 396 MIBC samples included in TCGA were considered. The subtype information was collected from the international consensus taxonomic classification describing six subtypes: Basal/Squamous-like (Ba/Sq), Luminal papillary (LumP), Luminal non-Specify (LumNS), Luminal unstable (LumU), Stroma-rich, and Neuroendocrine-like (NE-like). Using MiXCR, we mapped the RNA read sequences to their respective B-cell receptor (BCR) and T-cell receptor (TCR) clonotypes. To evaluate the BCR and TCR differences among subtypes, we compared diversity measures (richness and diversity) using a Wilcoxon test and we performed a network analysis to characterize the clonal expansion. For the survival analysis stratified by subtypes, Cox regression models adjusted for age, region, and pathological stage were performed.

Results: Overall, we found different patterns of tumor-infiltrating immune repertoire among the different MIBC subtypes. Stroma-rich and Ba/Sq tumors showed the highest BCR and TCR infiltration while LumP showed the lowest. In addition, we observed that the Ba/Sq and Stroma-rich tumors were more clonally expanded than the Luminal subtypes. Moreover, higher TCR richness and diversity were significantly associated with better survival in the Stroma-rich and Ba/Sq subtypes.

Discussion: This study provides evidence that MIBC subtypes present differences in the tumor microenvironment, in particular, the Ba/Sq and the Stroma-rich are related with a higher tumoral-infiltrating immune repertoire, which seems to be translated into better survival. Determining the causes of the different tumoral-infiltrating immune repertoire according to the MIBC molecular subtypes will help to improve our understanding of the disease and the distinct responses to immunotherapy of MIBC.

Introduction

Bladder cancer (BC) is the fourth most common cancer in Northern America and Europe among men and its incidence is still rising (1, 2)⁠. Urothelial bladder cancer (UBC) morphology represents 95% of BC. Overall, UBC is considered an immunogenic tumor due to its relatively high tumor mutational burden (3)⁠ and its responsiveness to Bacillus Calmette–Guerin (BCG) bladder instillations and checkpoint inhibitors (4)⁠. However, not all patients benefit from these therapies, possibly, because BC is not a single disease.

Most of UBC (80%) are diagnosed as non-muscle invasive tumors. While this is a milder form of UBC, a large proportion (40%) of patients suffer of multiple recurrences with some of them invading the detrusor muscle (MIBC), this being a life-threat event requiring a more aggressive treatment (5, 6)⁠. MIBC has further been classified according to somatic DNA-based and RNA-based features. Regarding the latter, several taxonomic classifications with different numbers and names of MIBC subtypes have been proposed (7–17). Recently, Kamoun et al. published an international consensus paper based on a network-based analysis done with 1750 MIBC transcriptomic profiles from six independent MIBC studies. The authors reported up to six subtypes: Luminal papillary (LumP), Luminal unstable (LumU), Luminal Non-Specified (LumNS), Stroma-rich, Basal/Squamous-like (Ba/Sq), and Neuroendocrine-like) (17). Interestingly, it has been observed that each of the subtypes has distinct differentiation patterns, oncogenic mechanisms, tumor microenvironments, as well as histological and clinical associations.

The network of immunoregulatory pathways in bladder cancer is quite complex, with several immune mechanisms arresting the effective antitumor T-cell response. T cells and dendritic cells (DCs) expressed inhibitory receptors in their membranes, repressing tumor growth. Type 1 T helper (TH1) cells favor the generation of an anti-tumor immune response. However, type 2 T helper (TH2) cells favor pro-tumor immune responses. Other immune cell types that promote tumor development are myeloid-derived suppressor cells (MDSCs), M2 macrophages and regulatory T (Treg) cells. In addition, Mast cells have been implicated in an indirect pro-tumor role, although the mechanism remains unclear (18).

This immune infiltration varies according to bladder cancer stages. Thus, non-muscle invasive bladder cancer (NMIBC) and MIBC show significant differences in the infiltration of immune cells. Furthermore, Kamoun et al. characterized the tumor microenvironment for the different MIBC subtypes using cell deconvolution tools and they observed that the Ba/Sq and the Stroma-rich subtypes had higher immune and stromal infiltration as well as distinct immune cell populations than the rest of subtypes. Even though the immune infiltration was mainly found within these two subtypes, it showed distinct immune cell populations. Ba/Sq tumors were enriched in cytotoxic lymphocytes and natural killer cells, whereas stroma-rich tumors overexpressed T- and B-cell markers. LumNS tumors were the only luminal type associated with immune infiltration signals; these were mainly for B and T lymphocytes (17). These differences could help in the selection of the patients for immunotherapies. (19)

In UBC, there are extensive evidences for an overall suppression of immunosurveillance responses within the tumor. However, little is known about the antigen specific responses (20). Among all the immune cell populations, B and T cells are key components of the adaptive immune response. T cells are involved in cell-mediated immunity, whereas B cells are primarily responsible for antibody responses against the specific antigens recognition through the B cell receptors (BCR) or immunoglobulins (Ig) Both receptors can recognize a large number of molecules. The BCR and TCR loci are form by recombining a set of variable (V), diversity (D) and joining (J) gene segments and its diversity is mainly concentrated in the complementary-determining region 3 (CDR3).

The BCR are made up of two heavy chains (IGH) and two light chains, the kappa (κ) chains (IGK) and the lambda (λ) chains (IGL). The receptor diversification arises from two different processes: somatic recombination and somatic hypermutation (21)⁠. By contrast, T cell receptors (TCR) are either TCRαβ or TCRγδ. Approximately 95% of T cells express a TCRαβ receptor, consisting of a TCRα (TRA) and a TCRβ (TRB) chain. The remaining 5% are made by a TCRγ (TRD) and a TCRδ (TRG) chain. These TCR chains are highly diverse in their variable domains (22)⁠.

The immune cell infiltration harboring these receptors may play decisive roles at different stages of tumor development resulting in a tumor microenvironment containing different balances of T and B cell receptors, in addition to the cancer cells and surrounding stroma. Their impact on tumor progression and treatment response has been suggested. In fact, T-cell infiltration play a central role in modern immunotherapy response in bladder cancer, among other cancers (4, 23, 24)⁠, whereas the role of B-cell infiltration has yet to be defined (25)⁠. Furthermore, this infiltration could be used as immunological biomarkers which may drive towards patient stratification (26).

Thus, given the fact the MIBC subtypes presented different immune microenvironment and clinical behavior, we hypothesized that the adaptive immune infiltration of these receptors is also distinct across the MIBC subtypes. Therefore, our aim was to further characterize the MIBC immune microenvironment by analyzing the tumor-infiltrating B- and T- cell repertoire according to the tumor taxonomic molecular subtypes towards a better understanding of MIBC progression pathways.

Material and methods

TCGA data

The study population included 404 MIBCs patients from TCGA with available consensus taxonomic subtype data by Kamoun et al. (17)⁠ based on RNAseq tumor data. Information on tumor gene expression (RNA-seq), demographic, and clinicopathological characteristics were retrieved through TCGA data portal (https://tcga-data.nci.nih.gov/tcga/). Mutational rate data was retrieved from Thorsson et al. (27)⁠. All of the patients provided informed consent to TCGA. Subtype information was directly extracted from the original published paper (17)⁠. The final number of patients analyzed was 396 (LumP=126, LumNS=20, LumU=53, Stroma-rich=45, and Ba/Sq=152). We excluded 6 neuroendocrine-like tumors because of the insufficient number for the posterior subtype analysis. Two patients were also discarded based on the out-ranged and very low values presented in the TCR reads (Table S1).

BCR and TCR data extraction

B-cell receptors (BCR) and T-cell receptors (TCR) were extracted from the RNAseq FASTQ files using the bioinformatic software MiXCR (28)⁠. We applied the pipeline described in https://mixcr.readthedocs.io/en/master/ for alignments using paired-end RNA-seq. MiXCR captures all CDRs and framework regions of immune genes and permits the assembly of full-length clonotypes. In this paper, we extracted BCR and TCR and defined the clonotypes according to their CDR3 sequences retrieved from the bulk RNA-seq data. The median number of reads and clones for the four datasets are displayed in Table S2. TRD and TRG reads and clones were very few, therefore we decided to filter out these receptors for the analysis.

Richness and diversity analyses

To evaluate the number of BCR/TCR clones and their frequency we assessed the richness and the diversity. They were calculated through Expression and Entropy measurements, respectively (29)⁠. The number of BCR/TCR reads can be highly dependent on the sequencing depth. Therefore, we accounted for this by calculating the expression dividing the number of reads by the total number of sequenced reads in the RNA-seq FASTQ files. Expression was estimated with the following formula:

$$I{g}_i/T{R}_i=\frac{M_i}{N_i+M_i};i=1,\text{...},N$$

Where M_i is the number of reads that map to a specific VDJ recombination and N_i is the number of reads that map to a anything else in the genome in n samples.

In addition, to take into consideration not only the number of clones but its frequency, Shannon entropy (H index) was estimated as:

$$H=-{\displaystyle \sum _{i=1}^N{p}_i{\log }_2(p_i);i=1,\text{...},N}$$

Where N is the number of unique clones and p_i is the frequency of clone i. We defined a clone as those reads that had the same V and J gene, same CDR3 length, and 90% of nucleotide identity for BCR, and 95% for TCR. We restricted this analysis to those reads that estimated the CDR3 region.

Network analysis

The network used to assess the overall clonal nature and the dominance of a clone was generated applying an algorithm similar to that already described in previous publications (29, 30)⁠. Briefly, each vertex represents a B-cell or T-cell sequence where the size indicates the number of identical chains. An edge (defined by the clone definition: same V and J segments, same CDR3 length, and 90%/95% nucleotide identity between CDR3s for BCR/TCR, respectively) between two vertex indicates that the sequence belongs to the same clone and clusters define each clone in the repertoire. The analysis was done using igraph package in R using the layout with_graphopt option to generate the plot.

Then, the network was quantified calculating the Gini Index for vertex size (clonal expansion) and cluster size (clonal diversification). Gini Index is a measure of unevenness extensively used to measure wealth distribution. A Gini coefficient of zero expresses perfect equality and a Gini coefficient of 1 expressed maximal inequality. It measures the inequality among values of frequency distribution. We used the Gini function from ineq package in R to calculate the Gini coefficient for vertex size and cluster size distribution. When applied to vertex size, the overall clonal nature is represented. If it was closer to 1, vertices were unequal, showing expansion of some of them, and closer to 0, otherwise. When applied to cluster size, clonal dominance was represented. If closer to 1, clusters were unequal and therefore represented dominant clones; if closer to 0, all clusters were of equal size. Finally, all these information was considered together to compare the clonal expansion and diversification trends by subtypes.

Statistical analysis

To evaluate the BCR and TCR differences among subtypes, we compared diversity measures (expression and entropy) using a Wilcoxon rank test. We also checked the correlation between the clinic-pathological variables available in TCGA and diversity measures for all receptors (IGH, IGK, IGL, TRA and TRB) stratifying by MIBC subtype and using Wilcoxon rank test when the variable was categorical and Spearman correlation test when continuous. In order to assess the correlation by MIBC subtypes between the mutational rates (silent and non-silent mutation rates, SMR and NSMR, respectively) and the diversity measures, a Spearman correlation test was applied.

The inflammatory score was calculated based on Thorsson et al. (27)⁠ scores calculation approach. The authors applied CIBERSORT to the TCGA RNA-seq data to estimate the relative fraction of 22 immune cell types within the leukocyte compartment. More specifically, they aggregate the cell types of interest to generate the different scores. According to this, we calculated the inflammatory score by adding the relative fraction of inflammation related immune cells (Inflammatory score = Monocytes + Macrophages.M0 + Macrophages.M1 + Macrophages.M2 + Dendritic.cells.resting + Dendritic.cells.activated + Mast.cells.resting + Mast.cells.activated + Neutrophils + Eosinophils + B.cells.naive + B.cells.memory + T.cells.CD4.naive + T.cells.CD4.memory.resting + T.cells.CD4.memory.activated + T.cells.follicular.helper + T.cells.regulatory + T.cells.gamma.delta + T.cells.CD8 + T.helper + NK.cells.resting + NK.cells.activated). To check the correlation between the inflammatory score and both measures, richness and diversity for all receptors (IGH, IGK, IGL, TRA and TRB) stratified by subtype, a Spearman correlation test was used.

Survival analysis were performed to assess the association between the diversity measures (expression and entropy) and overall survival (OS). Hazard ratios and 95% confidence intervals were estimated with Cox regression models adjusted for age, region, and pathological disease stage. The different BCR and TCR measurements were included in the models to evaluate their prognostic value. Cox model considering the potential interaction between the BCR and TCR measurements and the subtypes were also calculated.

Results

Sociodemographic features were similar across subtypes, the LumNS and LumU being slightly older than the rest of the subtypes (Table 1). As expected, the Stroma-rich and Ba/Sq subtypes had less papillary features. The LumNS subtype was the most advanced at diagnosis (65% of the tumors in Stage IV). The summary of the sequenced reads by MIBC subtypes is showed in Table S2. For all of MIBC subtypes, the number of reads was higher for the BCR compared to the TCR clonotypes. Stroma-rich subtype (mean(sd) = 41954 [2290 – 1312355]) and LumNS (mean(sd) = 24998 [2156 - 304394]) were the subtypes with the highest number of BCR reads, whereas the highest amount of TCR reads was shown by the Ba/Sq (mean(sd) = 404.5 [13 - 3975]) and Stoma-rich (mean(sd) = 308 [65 - 5753]) subtypes (Table S2).

TABLE 1

Table 1 Sociodemographic and clinical features for MIBC cases according to subtypes.

BCR/TCR infiltration is significantly different across MIBC subtypes

Overall, we found significant differences of BCR and TCR richness and diversity across MIBC subtypes (Figures 1, 2). The highest BCR infiltration was observed for the Stroma-rich subtype that, jointly with Ba/Sq tumors, showed the highest TCR richness and diversity, too. On the other hand, LumP tumors presented the lowest BCR and TCR infiltration pattern. Wilcoxon rank test results comparing the TCR and BCR diversity measures between the different subtypes are reported in Tables S3, S4.

FIGURE 1

Figure 1 BCR and TCR richness among the different MIBC subtypes. BCR related results are plotted in purple and TCR in yellow. In the Y axis the logarithm of the expression is represented. Each box of the boxplots is a subtype (see legend). The differences between richness across the subtypes are displayed only when significant. * → 0.05 > p > 0.01; ** → 0.01 > p > 0.001; *** → p < 0.001.

FIGURE 2

Figure 2 BCR and TCR diversity among the different MIBC subtypes. BCR related results are plotted in purple and TCR in yellow. In the Y axis entropy is represented. Each box of the boxplots is a subtype (see legend). The differences between diversity across the subtypes are displayed only when significant. Wilcoxon test: * → 0.05 > p > 0.01; ** → 0.01 > p > 0.001; *** → p < 0.001.

Significant differences in BCR clonal expansion and clonal diversification between the different subtypes were also found through the network analyses (Figure 3; Tables S5, S6). The Stroma-rich subtype showed the highest BCR clonal expansion levels while the LumP showed the lowest. The rest of the subtypes behaved similar in terms of BCR clonal expansion. The Stroma-rich subtype showed the highest levels of dominant clones for IGL and IGK and LumP showed the lowest. Higher levels of IGH dominant clones were observed for the LumNS subtype and the lowest levels were observed for the LumP, again.

FIGURE 3

Figure 3 BCR clonal expansion and diversification by subtypes. (A) On the Y axis, clonal expansion is represented by using the Gini Index and on the X axis, clonal diversification is displayed by Gini Vertex. The plot is colored by subtype, each dot is a MIBC sample and the boxplot represented the distribution of each clonal measurement by subtype. (B) B-cell repertoire networks from two samples representing one Ba/Sq (blue) and one LumP (green). Each vertex represents a unique BCR being the vertex size defined by the number of identical BCRs considering the nucleotide sequences. An edge exists between vertices when they belong to the same clone as defined before, so clusters are groups of interconnected vertices forming a clone.

Association between mutational rates and BCR/TCR infiltration varies across the different MIBC subtypes

While no correlation was observed between the mutational rates and the infiltration measures overall, the correlation patterns were highly heterogeneous across the different MIBC subtypes (Figures 4A, B, S1A, B; Tables 2, 3). In the Stroma-rich subtype, mutational rates showed a negative correlation trend with both BCR and TCR richness (NSMR rho: IGH=-0.24, IGK=-0.28, IGL=-0.27, TRA=-0.26, TRB=-0.24; SMR rho: IGH=-0.25, IGK=-0.29, IGL=-0.29, TRA=-0.30, TRB=-0.28) and TCR diversity (NSMR rho: TRA=-0.12, TRB=-0.26; SMR rho: TRA=-0.17, TRB=-0.31). A slightly positive correlation between the mutational rates and the TCR richness (NSMR rho: TRA=0.18, p.value=0.03; TRB=0.19, p.value=0.02; SMR rho: TRA=0.16, p.value=0.05; TRB=0.16, p.value=0.05) and diversity (NSMR rho: TRA=0.18, p.value=0.03; TRB=0.21 p.value=0.01; SMR rho: TRA=0.17, p.value=0.04; TRB=0.18 p.value=0.03) was found in the Ba/Sq. In addition, LumP showed a positive correlation with both BCR and TCR richness (NSMR rho: IGH=0.21, p.value=0.02; IGK=0.22, p.value=0.01; IGL=0.21, p.value=0.02; TRA=0.23, p.value=9.7e-03; SMR rho: IGH=0.19, p.value=0.03; IGK=0.20, p.value=0.02; IGL=0.20, p.value=0.03; TRA=0.22, p.value=0.01) and with TRA diversity (NSMR rho: 0.30, p.value=0.02; SMR rho: 0.28, p.value=4.1e-03) (Tables 2, 3).

FIGURE 4

Figure 4 Correlation between mutational rates and inflammatory score with richness by subtypes. BCR related results are plotted in purple and TCR in yellow. In the Y axis the logarithm of the expression is displayed. On the X axes, the logarithm 10 of the (A) non-silent (B) silent mutational rates, (C) inflammation score are plotted. Each line, is the regression line assessed in the correlation test performed by subtypes and they are colored by them.

TABLE 2

Table 2 Correlation between the two mutational rates and inflammatory score and richness for all receptors by subtypes.

TABLE 3

Table 3 Correlation between the two mutational rates and inflammatory score and diversity for all receptors by subtypes.

Ba/Sq BCR/TCR infiltration is significantly associated with inflammatory score

Ba/Sq subtype showed a significant positive correlation with both BCR and TCR richness (IGH: rho=0.20, p.value=1.6e-02; IGK: rho=0.22, p.value=8.2e-03; IGL: rho=0.21, p.value=9.7e-03, TRA: rho=0.47, p.value=2.4e-09; TRB: rho=0.36, p.value=7.4e-06) and TCR diversity (TRA: rho=0.43, p.value=9.6e-08; TRB: rho=0.40, p.value=1.4e-06). In addition, a significant positive correlation between the inflammatory score and TCR richness (TRA: rho=0.29, p.value=1.2e-03; TRB: rho=0.28, p.value=1.6e-03) and diversity (TRA: rho=0.33, p.value=6.1e-06; TRB: rho=0.37, p.value=9.1e-05) was observed for the LumP subtype (Figures 4C, S1C; Tables 2, 3). Interestingly, BCR and TCR richness and diversity did not show any pattern of correlation with the inflammatory score in the Stroma-rich subtype.

There is no association between the clinical data and BCR/TCR infiltration

We explored the association between the clinic-pathological variables available in TCGA and both measures, richness and diversity, for all receptors (IGH, IGK, IGL, TRA and TRB) by stratifying by subtype. No associations for any of the clinical variables was observed (Table S7).

BCR and TCR infiltration is associated with overall survival in Ba/Sq and Stroma-rich subtypes

Richness and diversity were associated with OS when we stratified by subtypes (Table S8 and Figure 5) showing an interaction effect between the BCR and TCR measurements with the different subtypes (Table S9). Ba/Sq subtype showed significantly better OS for higher TCR richness (HR[95% CI]: TRA=0.57 [0.37-0.86], TRB=0.53 [0.34-0.81]) and diversity (HR[95% CI]: TRA=0.79 [0.69-0.9], TRB=0.81 [0.71-0.92]) and BCR richness (HR[95% CI]: IGH=0.81 [0.63-1.03], IGK=0.77 [0.6-0.98], IGL=0.77 [0.59-1.00]), and Stroma-rich subtype showed significantly better OS for higher TCR richness (HR[95% CI]: TRA=0.23 [0.08-0.7], TRB=0.2 [0.06-0.62]) and diversity (HR[95% CI]: TRA=0.59 [0.43-0.82], TRB=0.62 [0.46-0.84]). Interestingly, OS was not associated with diversity measures if stratification by subtypes was not performed. BCR and TCR infiltration were not significantly associated with OS in either of Luminal subtypes.

FIGURE 5

Figure 5 Survival analyses considering richness and diversity for all receptors results by subtype. BCR related results are plotted in purple and TCR in yellow. Each square represents the Hazard Ratio and its corresponding lines the 95%CI. The different colors indicates the different subtypes (see legend) and all cases together are display in gray.

Discussion

Cancer classifications, such as the pathological, the molecular, and the taxonomic, aim to improve the patient management. Molecular subtyping studies have allowed the allocation of cancer into homogeneous groups that are considered to harbor similar molecular and clinical characteristics. Furthermore, this has helped researchers to identify both actionable targets for drug design as well as biomarkers for response prediction. In deep, molecular subtyping studies have allowed to better correlate cancer cases with clinical outcomes than the traditional classifications of cancer (31–33)⁠. MIBC is a heterogeneous disease with several taxonomic molecular subtypes showing different genetic, clinical, and epidemiological profiles (17, 34)⁠. It has also been suggested that MIBC-subtypes follow different tumorigenesis pathways playing decisive roles at different stages of tumor development and resulting in a tumor microenvironment containing different balances of adaptive immune cells (T and B lymphocytes) (17)⁠. In addition, the different MIBC subtypes have been associated with different therapeutics options (16). However, despite the growing evidences of subtyping clinical implications, MIBC subtypes have yet to enter into routine clinical practice (35).

While a very recent study detailed the UBC immune-profile, it did not considered the different molecular subtypes (19). To our knowledge, this is the first study characterizing the MIBC immune microenvironment by analyzing the tumor-infiltrating B and T repertoire according to the taxonomic molecular subtypes using RNAseq data from 396 MIBC samples included in TCGA. Furthermore, we report on the association of tumor-infiltrating immune repertoire with mutational rates, inflammatory score, overall survival, and clinico-pathological features, across the different MIBC subtypes.

We observed large differences on BCR and TCR richness and diversity, as well as on BCR clonal expansion and diversification among the different MIBC subtypes. Stroma-rich and Ba/Sq tumors showed the highest BCR and TCR infiltration while LumP subtype showed the lowest. In addition, we observed that the Ba/Sq and Stroma-rich tumors were more clonally expanded than the Luminal subtypes. Moreover, the correlation between mutational rates and inflammatory score with BCR/TCR measurements highly varied across subtypes.

A high infiltration in Ba/Sq tumors is further supported by the observation that Schistosoma-associated bladder cancer arises from a chronic granulomatous inflammation and irritation leading to squamous cell carcinoma subtype of the bladder (36)⁠ Urinary tract infections (UTIs) have been controversially established as risk factor for bladder cancer (37)⁠. Since it has been shown that during UTIs an adaptive immune response is generated (38)⁠, we could think that subtyping bladder cancer could help to better establish their relationship with bladder cancer.

The high infiltration of BCR and TCR in the Stroma- rich subtypes relies on the fact that these tumors displayed overexpression of smooth muscle, myofibroblast, fibroblast and endothelial gene signatures, intermediate urothelial differentiation and overexpression of B-cell markers (17, 39). The fact that molecular subtyping is performed on biopsy specimens representing only a fraction of the tumor mass, warrants particular caution when considering this subtype. While some tumors are actually stroma-rich, some biopsy specimens are stroma-rich due to chance or sampling at the tumor margin (40). Hence, it is necessary to be careful in drawing conclusions about this subtype.

Nevertheless, the fact that immune infiltration is associated with mutational rate in Stroma-rich tumors and with the inflammatory score in Ba/Sq subtype suggest that the source and type of immune infiltration may be different in these two subtypes. Our observation is also supported by Kamoun et al. that reported distinct immune cell populations in these two MIBC subtypes (17)⁠. Specifically, we found that the tumor-infiltrating immune repertoire inversely correlates with the mutational rates in the Stroma-rich subtype. This inverse correlation could be explained because this subtype is characterized by a high stromal infiltration, mainly of smooth cells, fibroblasts, and myofribroblasts that could cover the cancer cells where the mutations mainly originate (17)⁠.

In general, high mutational rates lead to the formation of tumor neoantigens, making tumors more immunogenic and more sensitive to immunotherapy (41)⁠. The above observations, jointly with the fact that higher TCR richness and diversity was associated with better survival only in the Stroma-rich and the Ba/Sq subtypes, reinforce the need to further explore the joint role of the tumor subtype and the immune infiltration in the anti-tumor response (42)⁠.

The three luminal subtypes showed the lowest BCR and TCR infiltration. This group of subtypes is characterized by a less aggressive presentation of the disease and a better prognosis (11, 43). Moreover, the immune infiltration patterns positively correlated with the mutational rates. Intriguingly, the correlation was only observed for the BCR richness and diversity, but not for the TCR measurements. Whether the type of mutations in the luminal subtypes are more immunogenic than those in the other subtypes requires further exploration.

The impact of the immune infiltration pattern on prognosis varies across the different subtypes. While BASQ-like subtype has been characterized by a more aggressive presentation of the disease and worse prognosis (11, 44)⁠, we were further able to differentiate a Ba/Sq tumors subset with better prognosis associated with TCR infiltration. Similar results were observed for Stroma-rich tumors, characterized by a better survival. The fact that BCR and TCR richness and diversity were not associated with OS in the luminal subtypes may explain why these tumors are poor responders to immunotherapy (45, 46)⁠, a fact that could be used for patient stratification in the clinics.

There are some limitations that should be considered when interpreting these results. First, we applied MiXCR tool to map the read sequences using RNAseq data to their respective BCR and TCR clonotypes as done elsewhere (29, 42)⁠. In this line, we have previously explored the tumor infiltrating B cell repertoires across tumor types (42)⁠ showing a large variability on BCR infiltration across tumor types and an increase clonality in primary tumors compared to adjacent non-tumor tissues. Despite these tools provide accurate annotations, further studies with targeted sequencing are necessary to validate the extracted features and associations. Another limitation of this study was the limited clinical data available in TCGA since the consortium’s main objective was the detail molecular characterization of the tumors. This impaired dawning clear conclusions from the lack of association between the tumor-infiltrating immune repertoire and clinic-pathological conditions. Hence, additional analyses in independent studies with more completed data are needed. We are also aware that some of the subtypes extracted from the MIBC samples had limited sample size (LumNS=20, NE-like=6). Indeed, NE-like tumors were excluded from the analysis due to this reason. Therefore, these findings will need to be further validated in a larger sample sized study. In addition, the characterization of the risk factors associated with the different subtypes would strongly increase the clinical and treatment potential significance of the findings.

This study provides sound evidence that MIBC subtypes present differences in the tumor B- and T-cell immune repertoire. In particular, the Stroma-rich and Ba/Sq tumors are related with a higher tumoral-infiltrating immune repertoire, however the origin of the immune infiltration may be different in these two subtypes. Interestingly, the Ba/Sq subtype immune infiltration correlated with inflammation-related cells infiltrating the tumor while the Stroma-rich immune infiltration correlated with the mutational rates. Importantly, BCR and TCR infiltration was associated with a better overall survival in both Ba/Sq and Stroma-rich subtypes. A better definition of the different immune-related mechanism leading to several MIBC taxonomic subtypes will improve our understanding of the disease and the identification of novel therapeutic strategies.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

Ethics statement

The studies involving human participants were reviewed and approved by TCGA data. The patients/participants provided their written informed consent to participate in this study.

Author contributions

SP, NM, and RB conceived the study design and analysis plan. RB performed the data analysis. KY, SP, and RB extracted the data with MiXCR and helped in the data analysis. SP, NM, and MS supervised the work. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the 2019 AACR-AstraZeneca Immuno-oncology Research Fellowship (19-40-12-PINE). The funder was not involved in the study design, collection, analysis, interpretation of data, the writing of this article or the decision to submit it for publication.

Acknowledgments

We thank to Lola Alonso for regarding the technical support and for useful discussion.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2023.986598/full#supplementary-material

References

1. Ferlay J, Soerjomataram I, Dikshit R, Eser S, Mathers C, Rebelo M, et al. Cancer incidence and mortality worldwide: Sources, methods and major patterns in GLOBOCAN 2012. Int J Cancer (2015) 136(5):E359-E386. doi: 10.1002/ijc.29210

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Ferlay J, Colombet M, Soerjomataram I, Mathers C, Parkin DM, Piñeros M, et al. Estimating the global cancer incidence and mortality in 2020: GLOBOCAN sources and methods. Int J Cancer (2019) 144(8):1941–53. doi: 10.1002/ijc.31937

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Schrock AB, Fabrizio D, He Y, Chung J, Resnick M, Stephens PJ, et al. Analysis of POLE mutation and tumor mutational burden (TMB) across 80,853 tumors: Implications for immune checkpoint inhibitors (ICPIs). Ann Oncol (2017) 28:517-529. doi: 10.1093/annonc/mdx376.035

CrossRef Full Text | Google Scholar

4. van Wilpe S, Gerretsen ECF, van der Heijden AG, de Vries IJM, Gerritsen WR, Mehra N. Prognostic and predictive value of tumor-infiltrating immune cells in urothelial cancer of the bladder. Cancers (2020) 12(9). doi: 10.3390/cancers12092692

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Rouprêt M, Babjuk M, Compérat E, Zigeuner R, Sylvester R, Burger M, et al. European Guidelines on upper tract urothelial carcinomas: 2013 update. Eur Urol (2013) 63(6):1059-1071. doi: 10.1016/j.eururo.2013.03.032

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Smith AB, Deal AM, Woods ME, Wallen EM, Pruthi RS, Chen RC, et al. Muscle-invasive bladder cancer: Evaluating treatment and survival in the national cancer data base. BJU Int (2014) 114(5):719–26. doi: 10.1111/bju.12601

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Sjodahl G, Lauss M, Lovgren K, Chebil G, Gudjonsson S, Veerla S, et al. A molecular taxonomy for urothelial carcinoma. Clin Cancer Res (2012) 18(12):3377–86. doi: 10.1158/1078-0432.CCR-12-0077-T

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Damrauer JS, Hoadley KA, Chism DD, Fan C, Tiganelli CJ, Wobker SE, et al. Intrinsic subtypes of high-grade bladder cancer reflect the hallmarks of breast cancer biology. Proc Natl Acad Sci (2014) 111(8):3110–5. doi: 10.1073/pnas.1318376111

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Rebouissou S, Bernard-Pierrot I, de Reynies A, Lepage M-L, Krucker C, Chapeaublanc E, et al. EGFR as a potential therapeutic target for a subset of muscle-invasive bladder cancers presenting a basal-like phenotype. Sci Trans Med (2014) 6(244):244ra91. doi: 10.1126/scitranslmed.3008970

CrossRef Full Text | Google Scholar

10. Weinstein JN, Akbani R, Broom BM, Wang W, Verhaak RGW, McConkey D, et al. Comprehensive molecular characterization of urothelial bladder carcinoma. Nature (2014) 507(7492):315–322. doi: 10.1038/nature12965

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Choi W, Porten S, Kim S, Willis D, Plimack ER, Roth B, et al. Identification of distinct basal and luminal subtypes of muscle- invasive bladder cancer with different sensitivities to frontline chemotherapy. Cancer Cell (2015) 25(2):152–65. doi: 10.1016/j.ccr.2014.01.009.Identification

CrossRef Full Text | Google Scholar

12. Dyrskjøt L, Reinert T, Algaba F, Christensen E, Nieboer D, Hermann GG, et al. Prognostic impact of a 12-gene progression score in non–muscle-invasive bladder cancer: A prospective multicentre validation study. Eur Urol (2017) 72(3):461-469. doi: 10.1016/j.eururo.2017.05.040

PubMed Abstract | CrossRef Full Text | Google Scholar

13. March-Vila E, Pinzi L, Sturm N, Tinivella A, Engkvist O, Chen H, et al. On the integration of in silico drug design methods for drug repurposing. Front Pharmacol (2017) 8:298(MAY). doi: 10.3389/fphar.2017.00298

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Marzouka NAD, Eriksson P, Rovira C, Liedberg F, Sjödahl G, Höglund M. A validation and extended description of the Lund taxonomy for urothelial carcinoma using the TCGA cohort. Sci Rep (2018) 8(1). doi: 10.1038/s41598-018-22126-x

CrossRef Full Text | Google Scholar

15. Mo Q, Nikolos F, Chen F, Tramel Z, Lee YC, Hayashi K, et al. Prognostic power of a tumor differentiation gene signature for bladder urothelial carcinomas. J Natl Cancer Institute (2018) 110(5):448–459. doi: 10.1093/JNCI/DJX243

CrossRef Full Text | Google Scholar

16. Robertson AG, Kim J, Al-Ahmadie H, Bellmunt J, Guo G, Cherniack AD, et al. Comprehensive molecular characterization of muscle-invasive bladder cancer. Cell (2018) 171(3):540–56. doi: 10.1016/j.cell.2017.09.007

CrossRef Full Text | Google Scholar

17. Kamoun A, de Reyniès A, Allory Y, Sjödahl G, Robertson AG, Seiler R, et al. A consensus molecular classification of muscle-invasive bladder Cancer[Formula presented]. Eur Urol (2020) 77(4):420-433. doi: 10.1016/j.eururo.2019.09.006

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Schneider AK, Chevalier MF, Derré L. The multifaceted immune regulation of bladder cancer. Nat Rev Urol (2019) 16(10):613–630. doi: 10.1038/s41585-019-0226-y

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Viveiros N, Flores BC, Lobo J, Martins-Lima C, Cantante M, Lopes P, et al. Detailed bladder cancer immunoprofiling reveals new clues for immunotherapeutic strategies. Clin Transl Immunol (2022) 11(9):e1402. doi: 10.1002/cti2.1402

CrossRef Full Text | Google Scholar

20. Joseph M, Enting D. Immune responses in bladder cancer-role of immune cell populations, prognostic factors and therapeutic implications. Front Oncol (2019) 9:1270. doi: 10.3389/fonc.2019.01270

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Georgiou G, Ippolito GC, Beausang J, Busse CE, Wardemann H, Quake SR. The promise and challenge of high-throughput sequencing of the antibody repertoire. Nat Biotechnol (2014) 32(2):158–168. doi: 10.1038/nbt.2782

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Garcia KC, Teyton L, Wilson IA. Structural basis of T cell recognition. Annu Rev Immunol (1999) 17:369-397. doi: 10.1146/annurev.immunol.17.1.369

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Zou W. Regulatory T cells, tumour immunity and immunotherapy. Nat Rev Immunol (2006) 6(4):295–307. doi: 10.1038/nri1806

PubMed Abstract | CrossRef Full Text | Google Scholar

24. O’Donnell JS, Teng MWL, Smyth MJ. Cancer immunoediting and resistance to T cell-based immunotherapy. Nat Rev Clin Oncol (2019) 16(3):151–167. doi: 10.1038/s41571-018-0142-8

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Largeot A, Pagano G, Gonder S, Moussay E, Paggetti J. The b-side of cancer immunity: The underrated tune. Cells (2019) 8(5):449. doi: 10.3390/cells8050449

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Batista R, Vinagre N, Meireles S, Vinagre J, Prazeres H, Leão R, et al. Biomarkers for bladder cancer diagnosis and surveillance: A comprehensive review. Diagn (Basel Switzerland) (2020) 10(1):39. doi: 10.3390/diagnostics10010039

CrossRef Full Text | Google Scholar

27. Thorsson V, Gibbs DL, Brown SD, Wolf D, Bortone DS, Ou Yang TH, et al. The immune landscape of cancer. Immunity (2018) 48(4):812-830. doi: 10.1016/j.immuni.2018.03.023

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Bolotin DA, Poslavsky S, Mitrophanov I, Shugay M, Mamedov IZ, Putintseva EV, et al. MiXCR: Software for comprehensive adaptive immunity profiling. Nat Methods (2015) 12(5):380–381. doi: 10.1038/nmeth.3364

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Pineda S, López de Maturana E, Yu K, Ravoor A, Wood I, Malats N, et al. Tumor-infiltrating b- and T-cell repertoire in pancreatic cancer associated with host and tumor features. Front Immunol (2021) 12:730746. doi: 10.3389/fimmu.2021.730746

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Bashford-Rogers RJM, Palser AL, Huntly BJ, Rance R, Vassiliou GS, Follows GA, et al. Network properties derived from deep sequencing of human b-cell receptor repertoires delineate b-cell populations. Genome Res (2013) 23(11):1874-1884. doi: 10.1101/gr.154815.113

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Nutt CL, Mani DR, Betensky RA, Tamayo P, Cairncross JG, Ladd C, et al. Gene expression-based classification of malignant gliomas correlates better with survival than histological classification. Cancer Res (2003) 63(7):1602–7. Available at: https://aacrjournals.org/cancerres/article/63/7/1602/511083/Gene-Expression-based-Classification-of-Malignant.

PubMed Abstract | Google Scholar

32. Fallahpour S, Navaneelan T, De P, Borgo A. Breast cancer survival by molecular subtype: a population-based analysis of cancer registry data. CMAJ Open (2017) 5(3):E734–9. doi: 10.9778/cmajo.20170030

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Lotan Y, de Jong JJ, Liu VYT, Bismar TA, Boorjian SA, Huang HC, et al. Patients with muscle-invasive bladder cancer with nonluminal subtype derive greatest benefit from platinum based neoadjuvant chemotherapy. J Urol (2022) 207(3):541–50. doi: 10.1097/JU.0000000000002261

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Sun X, Hoadley KA, Kim WY, Furberg H, Olshan AF, Troester MA. Age at diagnosis, obesity, smoking, and molecular subtypes in muscle-invasive bladder cancer. Cancer Causes Control (2017) 28(6):539–544. doi: 10.1007/s10552-017-0885-z

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Sjödahl G, Jackson CL, Bartlett JM, Siemens DR, Berman DM. Molecular profiling in muscle-invasive bladder cancer: more than the sum of its parts. J Pathol (2019) 247(5):563–73. doi: 10.1002/path.5230

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Zaghloul MS, Gouda I. Schistosomiasis and bladder cancer: Similarities and differences from urothelial cancer. Expert Rev Anticancer Ther (2012) 12(6):753-763. doi: 10.1586/era.12.49

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Bayne CE, Farah D, Herbst KW, Hsieh MH. Role of urinary tract infection in bladder cancer: a systematic review and meta-analysis. World J Urol (2018) 36(8):1181–1190. doi: 10.1007/s00345-018-2257-z

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Mora-Bau G, Platt AM, van Rooijen N, Randolph GJ, Albert ML, Ingersoll MA. Macrophages subvert adaptive immunity to urinary tract infection. PloS Pathog (2015) 11(7). doi: 10.1371/journal.ppat.1005044

PubMed Abstract | CrossRef Full Text | Google Scholar

39. McConkey DJ, Choi W. Molecular subtypes of bladder cancer. Curr Oncol Rep (2018) 20(10):77. doi: 10.1007/s11912-018-0727-5

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Sjödahl G, Liedberg F, Höglund M, Eriksson P. When the molecular subtype is hidden behind a veil of stroma. Eur Urol (2021) 80(2):160–1. doi: 10.1016/j.eururo.2021.04.014

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Wu Z, Zhu K, Liu Q, Liu Y, Chen L, Cui J, et al. Profiles of immune infiltration in bladder cancer and its clinical significance: an integrative genomic analysis. Int J Med Sci (2020) 17(6):762-772. doi: 10.7150/ijms.42151

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Yu K, Ravoor A, Malats N, Pineda S, Sirota M. A pan-cancer analysis of tumor-infiltrating b cell repertoires. Front Immunol (2022) 12:790119. doi: 10.3389/fimmu.2021.790119

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Choi W, Czerniak B, Ochoa A, Su X, Siefker-Radtke A, Dinney C, et al. Intrinsic basal and luminal subtypes of muscle-invasive bladder cancer. Nat Rev Urol (2014) 11(7):400–410. doi: 10.1038/nrurol.2014.129

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Yuk HD, Jeong CW, Kwak C, Kim HH, Moon KC, Ku JH. Clinical outcomes of muscle invasive bladder cancer according to the BASQ classification. BMC Cancer (2019) 19(1). doi: 10.1186/s12885-019-6042-1

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Rosenberg JE, Hoffman-Censits J, Powles T, van der Heijden MS, Balar AV, Necchi A, et al. Atezolizumab in patients with locally advanced and metastatic urothelial carcinoma who have progressed following treatment with platinum-based chemotherapy: A single-arm, multicentre, phase 2 trial. Lancet (2016) 387(10031):1909-1920. doi: 10.1016/S0140-6736(16)00561-4

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, et al. TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature (2018) 554(7693):544–548. doi: 10.1038/nature25501

PubMed Abstract | CrossRef Full Text | Google Scholar