Bladder Cancer Exhibiting High Immune Infiltration Shows the Lowest Response Rate to Immune Checkpoint Inhibitors

Background: Bladder urothelial cancer (BLCA) treatment using immune checkpoint inhibitors (IMCIs) can result in long-lasting clinical benefits. However, only a fraction of patients respond to such treatment. In this study, we aimed to identify the relationships between immune cell infiltration levels (ICILs) and IMCIs and identify markers for ICILs.

Methods: ICILs were estimated based on single-sample gene set enrichment analysis. The response rates of different ICILs to IMCIs were calculated by combining the ICILs of molecular subtypes in BLCA with the response rates of different molecular subtypes of IMvigor 210 trials to a programmed cell death ligand-1 inhibitor. Weighted gene co-expression network analysis was used to identify modules of interest with ICILs. Functional enrichment analysis was performed to functionally annotate the modules. Screening of key genes and unsupervised clustering were used to identify candidate biomarkers. Tumor IMmune Estimation Resource was used to validate the relationships between the biomarkers and ICILs. Finally, we verified the expression of key genes in molecular subtypes of different response rates for IMCIs.

Findings: The basal squamous subtype and luminal infiltrated subtype, which showed low response rates for IMCIs, had the highest levels of immune infiltration. The neuronal subtypes, which showed the highest response rates to IMCIs, had low ICILs. The modules of interest and key genes were determined based on topological overlap measurement, clustering results, and inclusion criteria. Modules highly correlated with ICILs were mainly enriched in immune responses and epithelial–mesenchymal transition. After screening the key genes in the modules, five candidate biomarkers (CD48, SEPT1, ACAP1, PPP1R16B, and IL16) were selected by unsupervised clustering. The key genes were inversely associated with tumor purity and were mostly expressed in the basal squamous subtype and luminal infiltrated subtypes.

Interpretation: Patients with high ICILs may benefit the least from treatment with IMCIs. Five key genes could predict ICILs in BLCA, and their high expression suggested that the response rate to IMCIs may decrease.

Introduction

Bladder urothelial cancer (BLCA) is one of the 10 most common malignancies worldwide. Systemic platinum-based chemotherapy, introduced nearly 30 years ago, remains the standard of care for untreated patients with inoperable or advanced metastatic BLCA and is associated with a 5-year survival rate of less than or equal to 15% (1). The landscape for BLCA has recently shifted following the approval of therapies targeting the programmed cell death-1 (PD-1)/programmed cell death ligand 1 (PD-L1) axis (2). Anti-immune checkpoint therapy has tremendously expanded our knowledge of the immunobiology of cancer. However, response rates to immune checkpoint inhibitors (IMCIs) are approximately 20% in both the platinum refractory and previously untreated settings (3). Tumor immune cell therapy may be the key to combatting cancer; however, the relationship between immune cell infiltration levels (ICILs) and the therapeutic effects of IMCIs is not clear.

The therapeutic effects of IMCIs are associated with the expression of PD-L1 and tumor-infiltrating lymphocytes (TILs) (4), as well as tumor mutational burden (TMB) (5). Notably, even if TILs are similar, the tumor microenvironment (TME) may still vary, and differences in mixed infiltration vs. septal infiltration may explain these variations (6). Moreover, there are many types of TILs, including CD8⁺ T cells, which kill tumor cells, and regulatory T cells (Tregs), which inhibit CD8⁺ T cells. Subtle changes in the proportions of immune cells can have different effects on tumor progression (7). Overall, tumor escape from immune surveillance and tumor counterattacks against immune cells are extremely complex and multifactorial processes.

Molecular subtypes of muscle invasive bladder cancer (MIBC), classified using The Cancer Genome Atlas (TCGA) data, revealed different ICILs and suggested that the basal squamous subtype and luminal infiltrated subtype can be treated with IMCIs (8). Indeed, these subtypes have been reported to respond to immune checkpoint therapy (9, 10). However, a recent study found that the neuronal subtype, which is associated with a poor prognosis, had the highest response rate after using the PD-L1 inhibitor atezolizumab, and prognosis was significantly better than those of the other four subtypes (11). The complete response (CR) and partial response (PR) of the neuronal subtype reached 100%, whereas those of the other three subtypes, excluding the luminal subtype, did not exceed 20%.

IMCIs may have applications as novel treatments for patients with BLCA. However, the relationships between ICILs and the efficacy of IMCIs have not been elucidated. Accordingly, in this study, we examined that ICILs could be used to predict the therapeutic effects of IMCIs and evaluated candidate biomarkers for ICILs through network analysis and unsupervised clustering.

Materials and Methods

Data Sources and Preprocessing

The RNA-sequencing (RNA-seq) results from 433 tissues and 408 cases of human bladder transitional cell carcinoma and papilloma samples, as well as data on the clinical characteristics of the patients, were obtained from TCGA database (portal.gdc.cancer.gov). Five MIBC molecular subtypes based on TCGA data were derived from a study by Robertson et al. (8). The response rates of different molecular subtypes to IMCIs were described by Kim et al. (11). In single-sample gene set enrichment analysis (ssGSEA), we normalized the expression data to Transcripts Per Million after calculating the total lengths of all exons for all transcripts. In weighted gene co-expression network analysis (WGCNA), we normalized the data using R package edgeR and filtered differentially expressed genes (DEGs) with an absolute fold change of >1 and a false discovery rate of <0.05.

Gene Signatures

A previously described procedure was used to determine the infiltration of immune cells in BLCA (12). We obtained a marker gene set for immune cell types from Bindea et al. (13). To calculate the one-sample gene set enrichment, we used the GSEA program to obtain the absolute enrichment scores from previously validated gene signatures.

Implementation of ssGSEA

The infiltration levels of immune cell types were quantified by ssGSEA in R package gsva. ssGSEA applies gene signatures expressed by immune cell populations (13) to individual cancer samples. The deconvolution approach used in our study included 24 immune cell types that were involved in innate immunity [natural killer (NK) cells, NK CD56^dim cells, NK CD56^bright cells, dendritic cells (DCs), plasmacytoid DCs, immature DCs, activated DCs (aDCs), neutrophils, mast cells, eosinophils, and macrophages] and adaptive immunity [B, T, CD8⁺ T, T helper (Th), Th1, Th2, Th17, T gamma delta, T central memory, T effector memory, T follicular helper (Tfh), Tregs, and cytotoxic T cells]. The obtained cytolytic activity (CYT) score obtained from the data set of Rooney et al. (14) consisted of cytolytic genes and was calculated as the geometrical means for perforin 1 and granzyme A. These two key cytolytic effectors are dramatically upregulated upon CD8⁺ T-cell activation and during productive clinical responses to anti-CTLA-4 and anti-PD-L1 immunotherapies; these two molecules are co-expressed in TCGA samples (14). We used the CYT scores to evaluate immune infiltration in different samples.

WGCNA and Module Preservation

WGCNA was performed using the WGCNA R package (15). Because some genes with no significant changes in expression between samples were highly correlated in WGCNA, genes with the most differential expression were used in subsequent WGCNA. Genes with the highest 25% of DEG variance were selected (16), guaranteeing the heterogeneity and accuracy of bioinformatics statistics for further co-expression network analysis. First, RNA-seq data were filtered to reduce outliers. The co-expression similarity matrix consisted of the absolute values of the correlation between transcript expression levels. A Pearson correlation matrix was constructed for paired genes. We constructed a weighted adjacency matrix using the power function amn = |cmn|β (cmn = Pearson correlation between gene m and gene n; amn = adjacency between gene m and gene n). The parameter β emphasized a strong correlation between genes and penalized a weak correlation. Next, an appropriate β value was selected to increase the similarity matrix and achieve a scale-free co-expression network. The adjacency matrix was then converted into a topological overlap matrix (TOM), which measured the network connectivity of genes defined as the sum of adjacent genes generated by all other networks. Average linkage hierarchical clustering was performed based on TOM-based dissimilarity measurements, and the minimum size (genome) of the gene dendrogram was 30. Through further analysis of modules, we calculated their dissimilarity and constructed module dendrograms.

Confirmation of Significant Modules

To determine the significance of each module, gene significance (GS) was calculated to measure the correlation between genes and sample traits. Module eigengenes were considered the main components in the principal component analysis of each gene module, and the expression patterns of all genes were summarized as a single feature expression profile within a given module. Next, GS was defined as the log10 conversion of the p-value in the linear regression between gene expression and clinical data (GS = lgP). Module significance (MS) was defined as the average GS within the module and calculated to measure the correlation between the module and sample traits. Statistical significance was determined using the relevant p-values. In order to increase the capacity of the modules, we selected a cutoff (<0.25) to merge some modules with similar heights. Next, we selected the ICILs previously calculated by ssGSEA and some clinical data for the clinical phenotype. The gene modules associated with the clinical phenotype were then analyzed.

Key Gene Identification

After selecting modules of interest, we calculated GS and module membership (MM, correlation between the module own genes and gene expression profiles) for each key gene and set their thresholds. In most similar studies, the thresholds for screening key genes in the module were defined as cor. gene MM >0.8 and cor. gene GS >0.2 (17, 18), and we adjusted these values appropriately to suit our study.

Functional Annotation—FunRich

FunRich (www.funrich.org) is an independent software tool for functional enrichment and interaction network analysis of genes and proteins (19). To explore the biological functions of genes in modules, we used the FunRich enrichment analysis function for analysis.

Tumor IMmune Estimation Resource (TIMER)

TIMER (cistrome.shinyapps.io/timer) provides a user-friendly web interface for dynamic analysis and visualization of the associations between immune infiltrates and gene expression (20). We used the Gene module to verify the association between genes and immune infiltration. The scatterplots were generated and displayed, showing Spearman's correlation and statistical significance.

Validation of Key Genes in the Different Molecular Subtypes

Few expression profiling, sequencing, or array data have been published for different prognostic outcomes after treatment with IMCIs. Therefore, according to the different responses of various molecular subtypes to IMCIs, we identified data with molecular subtypes to verify the key genes.

Results

Immune Phenotype Landscape in BLCA

Diverse immune cell populations infiltrate the TME and activate or suppress the antitumor response. To assess the spectrum of immune cell infiltration, the ssGSEA (21) approach was utilized to deconvolve the relative abundance of each cell type based on expression profiling data retrieved from the TCGA database. In this analysis, 408 patients with BLCA, for whom transcriptome profiling data and clinical characteristics were available, were included in this study. We found significant heterogeneity in terms of infiltration of numerous immune cell types among the cohort (Figure 1A). To facilitate further characterization, unsupervised clustering was applied to categorize the cohort into three infiltration subgroups, i.e., termed high (n = 55), median (n = 195), and low (n = 183) infiltration.

FIGURE 1

Figure 1. Immune infiltration landscape of BLCA. (A) Unsupervised clustering of 408 patients from The Cancer Genome Atlas cohort using single-sample gene set enrichment analysis scores from 24 immune cell types. Molecular subtype, sample type, stage, histological, survival, gender, anatomic location, and race are shown in the lower panel. Hierarchical clustering was performed with Euclidean distance and Ward linkage analyses. Three distinct immune infiltration clusters, termed high infiltration, moderate infiltration, and low infiltration, were defined. (B) Relative cytolytic scores (CYT scores) for tumors with high, moderate, and low immune infiltration were clustered according to overall immune cell infiltration. Kruskal–Wallis tests were used for analyses, and results with P-values of < 0.0001 were considered significant. Error bars represent the medians, interquartile ranges, and minimum to maximum values. (C) The evolution of immune infiltration levels in 19 cases with paracancerous tissues.

Considering the concomitant infiltration of activating and suppressive immune cell types, we investigated whether higher ICILs were correlated with elevation of cytotoxic function. To this end, CYT, which can serve as a surrogate index for quantifying the magnitude of the antitumor response, was examined. Patients with high infiltration status showed the highest CYT scores, indicating that cytotoxic function was efficiently elicited in those patients (Figure 1B).

A comparative analysis of 19 cases with normal tissues revealed that most normal tissues showed moderate infiltration. The corresponding cancer tissues showed a partial decrease in the ICILs (Figure 1C). Although the difference in ICILs was minor, we found that this result was related to lack of sufficient post-cluster ICILs grading. After increased grading, the ICILs in all cases were different. Therefore, in the next WGCNA, we analyzed DEGs.

ICILs of Different Molecular Subtypes

Although there were few cases of the neuronal subtype, the high response rate to IMCIs was interesting. There were only 20 samples of neuronal subtypes; 15 had low ICILs, and 5 had moderate ICILs. Thus, although 75% of cases of the neuronal subtype had low ICILs, we could not assume that low ICILs were responsive to IMCIs because most low ICILs were observed for the luminal-papillary subtype, which had a low response rate to IMCIs (Figure 2A). In addition to the neuronal subtype, the luminal subtype showed a better response rate than the other three subtypes, and 68 and 32% of cases of this subtype had low and moderate ICILs, respectively (Figure 2B). There were only three molecular subtypes showing high ICILs, and their response rates to IMCIs were all low; the basal squamous subtype and the luminal infiltrated subtype accounted for 98% of cases. Thus, these findings showed that high ICILs were associated with low response rate to IMCIs.

FIGURE 2

Figure 2. Molecular subtypes and immune infiltration. (A) Different ICILs showed variations in the proportions of molecular subtypes. (B) ICILs of different molecular subtypes were altered. (C) TCGA was similar to the IMvigor 210 trial in terms of total sample size and molecular subtype ratio. (D) The low ICIL group had the highest efficacy after receiving IMCIs. High ICILs had the lowest efficiency. ICILs, immune cell infiltration levels; IMCIs, immune checkpoint inhibitors.

CR and PR of Different ICILs

Based on the CR and PR ratios of different molecular subtypes, we calculated the response rates of different ICILs to IMCIs. The ratio of CR to PR was derived from the IMvigor 210 cohort study (11), which included 348 cases, and the proportion of each molecular subtype was similar to that of the 408 cases from the TCGA database (Figure 2C). Finally, we calculated the ratio of CR to PR for different ICILs in TCGA in combination with the ratio of each molecular subtype (Figure 2D). Although the response rates of different ICILs did not exceed 30%, the response rate to IMCIs increased with the decrease in ICILs. The percentage of low ICILs was 7% higher than that of high ICILs.

WGCNA: Identification of the Most Significant Modules and Genes

Wgcna was performed to construct a gene co-expression network for identification of biologically significant gene modules and elucidation of genes associated with ICILs. Eliminating outlier samples (Supplementary Figure 1A), the 8510 DEGs with the highest 25% of variance by cluster analysis were placed in a module. In this study, we selected β = 3 (scale-free R² = 0.957) as a soft threshold to ensure a scale-free network (Supplementary Figure 1B) and obtained 11 modules without merging for subsequent analysis (Figure 3A).

FIGURE 3

Figure 3. Weighted gene co-expression network of bladder urothelial cancer. (A) Identification of a co-expression module in bladder urothelial cancer. The branches of the cluster dendrogram correspond to the 11 different gene modules. Each piece of the leaves on the cluster dendrogram corresponds to a gene. (B) Correlations between the gene module and clinical traits. The correlation coefficient in each cell represents the correlation between the gene module and the clinical traits, which decreased in size from red to blue. The corresponding P-values are also shown. (C) Scatter plot of module eigengenes in pink, brown, and yellow modules. The red line indicates the screening threshold. The burgundy box identifies the key genes for each module.

To analyze the relationships between the modules and ICILs of the samples, we used MS as the overall gene expression level of the corresponding module to calculate the correlations with clinical phenotypes. The pink module was most significantly associated with the ICILs, with a correlation close to 0.65. In addition, the brown and yellow modules exhibited relatively high positive correlations with immune infiltration (Figure 3B). Thus, we chose these modules as the modules of interest and used them for subsequent analyses.

To determine the importance of genes in the module and the correlations between gene expression and ICILs, we appropriately adjusted the screening thresholds for MM and GS to select key genes. We identified 9 key genes (ACAP1, CD48, CXCR4, GYPC, IL16, NCKAP1L, PPP1R16B, SEPT1, and WIPF1, cor.MM >0.8 and cor.GS > 0.5) in the pink module, 4 key genes (PARP12, TAP1, TAP2, and TYMP, cor.MM >0.7 and cor.GS >0.4) in the brown module, and 10 key genes (COL6A2, CTHRC1, EMILIN1, FBN1, FSTL1, GLT8D2, LRRC32, MMP2, TGFB, and TIMP2, cor.MM >0.8 and cor.GS >0.4) in the yellow module (Figure 3C).

Functional Annotation of Modules

Although the three modules had high correlations with immune infiltration, the functions were not the same (Figure 4). Both the pink module and the brown module were related to immune responses. The pink module was also related to chemokine activity, and the brown module was related to interferon signaling, both of which are immune-related functions. However, the yellow module was different, showing associations with the extracellular matrix and the epithelial–mesenchymal transition (EMT). Although this result was surprising, it was consistent with previous reports demonstrating that EMT-related genes are related to resistance to PD-1 blockade in bladder cancer (21).

FIGURE 4

Figure 4. Enriched functions of the modules highly correlated with immune infiltration. Biological processes (BPs), cellular components (CCs), molecular functions (MFs), and biological pathways are listed. For detailed enrichment results, refer to the Supplementary Material.

Cluster Analysis

To assess the relationships among immune cells, immune-related genes, and key genes, unsupervised clustering was performed using ssGSEA (Figure 5). We noted that the genes from different modules were clustered in different zones that were close to each other. However, the key genes in the pink module were distributed in two cluster zones, suggesting that differences existed among these key genes, even though their biological functions were similar. CD48, SEPT1, and ACAP1 in the pink module were strongly associated with the toxicity of T cells, whereas PPP1R16B and IL16 also showed similar associations but were more closely related to B cells and Tfh cells. Among the key genes in the brown module, TAP1 and TAP2 encode components of MHC I. Together with PARP12 and TYMP, these key genes were closely related to aDCs and NK CD56^dim cells. The key genes in the yellow module, which was related to the mechanisms of the EMT and extracellular matrix, were clustered even more closely together and had the closest relationship with NK cells. Based on the results of these cluster analyses, we considered CD48, SEPT1, ACAP1, PPP1R16B, and IL16 in the pink module as important candidate genes because their expression levels best reflected the ICILs of the samples and cytotoxic function.

FIGURE 5

Figure 5. Five candidate biomarkers were found by unsupervised clustering, and the sources of genes identifying different immune function modules can be found in the Discussion section.

TIMER Validation and Co-expression

TIMER utilizes data from the TCGA database to facilitate the visualization of associations among genes and immune infiltration. CD48, SEPT1, ACAP1, PPP1R16B, and IL16 were all negatively correlated with tumor purity, supporting the results for highly correlating with ICILs (Figure 6A). Additionally, the co-expression of these five genes was analyzed using the R corplot package, and no correlations were <0.78; the correlation between ACAP1 and SEPT1 reached 0.97 (Figure 6B).

FIGURE 6

Figure 6. Verification of five candidate markers. (A) The correlation of five biomarkers with tumor purity and fraction of immune cells. (B) The five candidate biomarkers were highly correlated with each other at the transcript level. Correlation analysis using Spearman correlation; results with P-values of < 0.0001 were considered significant. (C) Candidate biomarkers were highly expressed in the basal squamous subtype and the luminal infiltrated subtype, and the ICILs of these two subtypes were the highest. Gene expression differed among subtypes (P < 0.01). Differences were analyzed by one-way ANOVA. (D) In GSE87304, the expression of the candidate genes in the previously described subtypes was still the highest.

Key Genes in Molecular Subtypes

There was a significant difference in the expression of the key genes in different molecular subtypes based on TCGA data, using one-way analysis of variance (ANOVA). Among the five molecular subtypes, five key genes were highly expressed in the basal and luminal infiltrated subtypes, both of which had low responses to IMCIs (Figure 6C). In Gene Expression Omnibus, GSE87304 with molecular subtype was screened out for key genes. The expression levels of the five key genes were higher in the basal and luminal infiltrated subtypes than in the luminal subtype (Figure 6D). The basal and luminal infiltrated subtypes had the highest ICILs and the worst IMCI response rate.

Discussion

The use of IMCIs as therapeutic agents for BLCA has resulted in favorable outcomes in some patients; however, most patients do not benefit from treatment with these agents. Tumor immunity is an extremely complex biological process, and factors affecting the efficacy of IMCIs include PD-L1 expression levels and TMB. In our study, we found that ICILs affected the response rate of IMCIs, and higher ICILs led to lower response rates. In order to easily identify the ICILs of BLCA, the modules obtained by WGCNA were related to the ICILs calculated by ssGSEA, and key genes were obtained through further screening and unsupervised clustering. High expression of the key genes reduced the effectiveness of IMCI treatment.

Our results showed that increased expression of key genes indicated decreased tumor purity, which was not solely caused by ICILs because there are many other cells, such as cancer-associated fibroblasts (CAFs), also found in the tumor (22). CAFs can promote the EMT and are also associated with tumor progression and treatment resistance (23, 24). Although CAFs are likely to be associated with ICILs, the yellow module, which was enriched in EMT and extracellular matrix functions, was highly correlated with ICILs. Therefore, CAFs may be responsible for the high resistance of IMCIs in high ICILs.

In this study, we identified a number of representative genes related to immunity and found that these genes were involved in the immune system through different mechanisms. During the development of tumor immunity, tumors show many alterations, including suppression of identity characteristics and stress, induction of immune cell apoptosis, attraction of Tregs, and neutralization of complement. These processes may block the immune system from clearing tumor cells, and the genes involved in these processes are initially aggregated for subsequent cluster analysis. The process of hiding the tumor identity involves major histocompatibility complex (MHC) I-related genes (e.g., HLA-A/B/C, TAP1/2, TAPBP, and B2M). The genes that hide stress to prevent NK cells from attacking include MICA, MICB, and NKG2D (25). Additionally, genes involved in counterattacking immune cells and inducing apoptosis include PD-1 (PDCD1), PD-L1 (CD274), and CTLA4 (26). FAS also binds to FASL to activate the non-inherent apoptotic pathway (27, 28), and tumor cells release transforming growth factor (TGF)-β and interleukin (IL)-10 to counterattack immune cells (29, 30). Tregs block lymphocytes from attacking tumor cells, which release C–C motif chemokine ligand 22 to bind to Tregs expressing C–C motif chemokine receptor 4, thereby attracting Tregs to tumors (31). CD46, CD55, and CD59 are membrane-bound complement regulatory proteins that prevent complement-mediated cytolysis (31, 32).

In our cluster analysis, we included genes that have been shown to be associated with ICILs. Among these genes, the well-known T-cell surface receptors PD-1 and CTLA4 (33) were clustered in the same group as the five important candidate genes: CD48 (encoding CD48), SEPT1 (encoding Septin 1), ACAP1 (encoding ArfGAP with coiled-coil, Ankyrin repeat, and PH domains 1), PPP1R16B (encoding protein phosphatase 1 regulatory subunit 16B), and IL16 (encoding IL-16). These genes were also closely related to CD8 T cells and cytotoxic cells. Moreover, in the final molecular subtype validation, their expression was highest in the basal and infiltrated subtypes of high ICILs. High expression of these molecules could indicate that the effects of IMCIs may decrease.

To elucidate the associations between these five genes and TILs, further investigations showed that except for PPP1R16B, all other genes had been previously reported to be associated with immune cells. CD48, a member of the signaling lymphocyte activation molecule family, participates in the adhesion and activation of immune cells, thereby contributing to T-cell activation and proliferation (34). SEPT1, a member of the septin family of GTPases, contributes to cancer progression and proliferation in oral squamous cell carcinoma (35). Septins are cytoskeletal proteins that provide compression and rigidity and support efficient motion of motile T cells (36). PPP1R16B (TIMAP), an endothelium-enriched TGF-β downstream protein that structurally belongs to the targeting subunit of myosin phosphatase, regulates macrophage M2 phenotypic phagocytosis (37). IL-16, a cytokine known for its chemotactic and inflammatory properties, induces proliferation in cutaneous T-cell lymphoma T cells and plasma cells in multiple myeloma and recruits CD4⁺ protumor macrophages in breast cancer (38, 39).

The brown module, which included the MHC components TAP1 and TAP2, and the yellow module, which was closely related to the mechanism of the EMT, were both associated with ICILs. Cytotoxic T cells are known to destroy cancer cells by recognizing the peptides presented via MHC-I on the tumor cell surface (40). In addition, the EMT not only is involved in tumor metastasis but also promotes immune escape and enhances immunosuppressive signals (41, 42). Therefore, ICILs in tumors are a manifestation of a very complex biological process. Our findings showed that high ICILs may have no advantage in the treatment of IMCIs.

In the calculation of the response rates of different ICILs for IMCIs, we showed that the response rates of different molecular subtypes in various ICILs were the same. This may have caused some bias in the results. However, in BLCA, few studies have evaluated IMCIs using sequencing or array profiling. Further studies are needed to verify our conclusions.

In summary, our results suggested that patients with high ICILs may benefit the least from treatment with IMCIs in BLCA. CD48, SEPT1, ACAP1, PPP1R16B, and IL16 were candidate genes for determining ICILs in BLCA, and these genes may have application as biomarkers to guide treatment with IMCIs.

Data Availability Statement

Publicly available datasets were analyzed in this study. These data can be found here: https://portal.gdc.cancer.gov.

Author Contributions

BL proposed the study concept, design, and drafted the manuscript. BL and SP collected, analyzed, and interpreted the data. YZ, XC, and BW participated in revising the manuscript. All authors have read and approved the final manuscript.

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.

Supplementary Material

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

Supplementary Figure 1. (A) Clustering of samples and removal of outliers. (B) Analysis of network topology for various soft-thresholding powers in scale independence and mean connectivity.

Supplementary Table 1. Enrichment analysis of brown module.

Supplementary Table 2. Enrichment analysis of pink module.

Supplementary Table 3. Enrichment analysis of yellow module.

References

1. Massard C, Gordon MS, Sharma S, Rafii S, Wainberg ZA, Luke J, et al. Safety and efficacy of durvalumab (medi4736), an anti-programmed cell death ligand-1 immune checkpoint inhibitor, in patients with advanced urothelial bladder cancer. J Clin Oncol. (2016) 34:3119–25. doi: 10.1200/JCO.2016.67.9761

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Chism DD. Urothelial carcinoma of the bladder and the rise of immunotherapy. J Natl Comprehen Cancer Netw. (2017) 15:1277–84. doi: 10.6004/jnccn.2017.7036

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Szabados B, van Dijk N, Tang YZ, van der Heijden MS, Wimalasingham A, Gomez de Liano A, et al. Response rate to chemotherapy after immune checkpoint inhibition in metastatic urothelial cancer. Eur Urol. (2018) 73:149–52. doi: 10.1016/j.eururo.2017.08.022

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Chen YP, Zhang Y, Lv JW, Li YQ, Wang YQ, He QM, et al. Genomic analysis of tumor microenvironment immune types across 14 solid cancer types: immunotherapeutic implications. Theranostics. (2017) 7:3585–94. doi: 10.7150/thno.21471

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Samstein RM, Lee CH, Shoushtari AN, Hellmann MD, Shen R, Janjigian YY, et al. Tumor mutational load predicts survival after immunotherapy across multiple cancer types. Nat Genet. (2019) 51:202–6. doi: 10.1038/s41588-018-0312-8

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Keren L, Bosse M, Marquez D, Angoshtari R, Jain S, Varma S, et al. A Structured tumor-immune microenvironment in triple negative breast cancer revealed by multiplexed ion beam imaging. Cell. (2018) 174:1373–87.e19. doi: 10.1016/j.cell.2018.08.039

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Buchan SL, Dou L, Remer M, Booth SG, Dunn SN, Lai C, et al. Antibodies to costimulatory receptor 4-1BB enhance anti-tumor immunity via T regulatory cell depletion and promotion of CD8 T cell effector function. Immunity. (2018) 49:958–70.e7. doi: 10.1016/j.immuni.2018.09.014

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Sharma P, Callahan MK, Bono P, Kim J, Spiliopoulou P, Calvo E, et al. Nivolumab monotherapy in recurrent metastatic urothelial carcinoma (CheckMate 032): a multicentre, open-label, two-stage, multi-arm, phase 1/2 trial. Lancet Oncol. (2016) 17:1590–8. doi: 10.1016/S1470-2045(16)30496-X

PubMed Abstract | CrossRef Full Text | Google Scholar

10. 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:1909–20. doi: 10.1016/S0140-6736(16)00561-4

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Kim J, Kwiatkowski D, McConkey DJ, Meeks JJ, Freeman SS, Bellmunt J, et al. The cancer genome atlas expression subtypes stratify response to checkpoint inhibition in advanced urothelial cancer and identify a subset of patients with high survival probability. Eur Urol. (2019) 75:961–4. doi: 10.1016/j.eururo.2019.02.017

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Zhang L, Zhao Y, Dai Y, Cheng JN, Gong Z, Feng Y, et al. Immune landscape of colorectal cancer tumor microenvironment from different primary tumor location. Front Immunol. (2018) 9:1578. doi: 10.3389/fimmu.2018.01578

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Bindea G, Mlecnik B, Tosolini M, Kirilovsky A, Waldner M, Obenauf AC, et al. Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity. (2013) 39:782–95. doi: 10.1016/j.immuni.2013.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Rooney MS, Shukla SA, Wu CJ, Getz G, Hacohen N. Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell. (2015) 160:48–61. doi: 10.1016/j.cell.2014.12.033

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinform. (2008) 9:559. doi: 10.1186/1471-2105-9-559

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Pan S, Zhan Y, Chen X, Wu B, Liu B. Identification of biomarkers for controlling cancer stem cell characteristics in bladder cancer by network analysis of transcriptome data stemness indices. Front Oncol. (2019) 9:613. doi: 10.3389/fonc.2019.00613

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Wang Y, Chen L, Wang G, Cheng S, Qian K, Liu X, et al. Fifteen hub genes associated with progression and prognosis of clear cell renal cell carcinoma identified by coexpression analysis. J Cell Physiol. (2018) 234:10225–37. doi: 10.1002/jcp.27692

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Chen L, Yuan L, Qian K, Qian G, Zhu Y, Wu CL, et al. Identification of biomarkers associated with pathological stage and prognosis of clear cell renal cell carcinoma by co-expression network analysis. Front Physiol. (2018) 9:399. doi: 10.3389/fphys.2018.00399

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Pathan M, Keerthikumar S, Ang CS, Gangoda L, Quek CY, Williamson NA, et al. FunRich: an open access standalone functional enrichment and interaction network analysis tool. Proteomics. (2015) 15:2597–601. doi: 10.1002/pmic.201400515

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, et al. TIMER: a web server for comprehensive analysis of tumor-infiltrating immune cells. Cancer Res. (2017) 77:e108–10. doi: 10.1158/0008-5472.CAN-17-0307

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Wang L, Saci A, Szabo PM, Chasalow SD, Castillo-Martin M, Domingo-Domenech J, et al. EMT- and stroma-related gene expression and resistance to PD-1 blockade in urothelial cancer. Nat Commun. (2018) 9:3503. doi: 10.1038/s41467-018-05992-x

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Gascard P, Tlsty TD. Carcinoma-associated fibroblasts: orchestrating the composition of malignancy. Genes Dev. (2016) 30:1002–19. doi: 10.1101/gad.279737.116

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Erdogan B, Webb DJ. Cancer-associated fibroblasts modulate growth factor signaling and extracellular matrix remodeling to regulate tumor metastasis. Biochem Soc Trans. (2017) 45:229–36. doi: 10.1042/BST20160387

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Kalluri R. The biology and function of fibroblasts in cancer. Nat Rev Cancer. (2016) 16:582–98. doi: 10.1038/nrc.2016.73

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Raulet DH. Roles of the NKG2D immunoreceptor and its ligands. Nat Rev Immunol. (2003) 3:781–90. doi: 10.1038/nri1199

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Havel JJ, Chowell D, Chan TA. The evolving landscape of biomarkers for checkpoint inhibitor immunotherapy. Nat Rev Cancer. (2019) 19:133–50. doi: 10.1038/s41568-019-0116-x

PubMed Abstract | CrossRef Full Text | Google Scholar

27. O'Connell J, Houston A, Bennett MW, O'Sullivan GC, Shanahan F. Immune privilege or inflammation? Insights into the Fas ligand enigma. Nat Med. (2001) 7:271–4. doi: 10.1038/85395

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Ryan AE, Shanahan F, O'Connell J, Houston AM. Addressing the “Fas counterattack” controversy: blocking fas ligand expression suppresses tumor immune evasion of colon cancer in vivo. Cancer Res. (2005) 65:9817–23. doi: 10.1158/0008-5472.CAN-05-1462

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Gorelik L, Flavell RA. Immune-mediated eradication of tumors through the blockade of transforming growth factor-beta signaling in T cells. Nat Med. (2001) 7:1118–22. doi: 10.1038/nm1001-1118

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Rojas JM, Avia M, Martin V, Sevilla N. IL-10: a multifunctional cytokine in viral infections. J Immunol Res. (2017) 2017:6104054. doi: 10.1155/2017/6104054

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Curiel TJ, Coukos G, Zou L, Alvarez X, Cheng P, Mottram P, et al. Specific recruitment of regulatory T cells in ovarian carcinoma fosters immune privilege and predicts reduced survival. Nat Med. (2004) 10:942–9. doi: 10.1038/nm1093

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Chen S, Caragine T, Cheung NK, Tomlinson S. CD59 expressed on a tumor cell surface modulates decay-accelerating factor expression and enhances tumor growth in a rat model of human neuroblastoma. Cancer Res. (2000) 60:3013–8.

PubMed Abstract | Google Scholar

33. Brunner-Weinzierl MC, Rudd CE. CTLA-4 and PD-1 control of T-cell motility and migration: implications for tumor immunotherapy. Front Immunol. (2018) 9:2737. doi: 10.3389/fimmu.2018.02737

PubMed Abstract | CrossRef Full Text | Google Scholar

34. McArdel SL, Terhorst C, Sharpe AH. Roles of CD48 in regulating immunity and tolerance. Clin Immunol. (2016) 164:10–20. doi: 10.1016/j.clim.2016.01.008

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Kato Y, Uzawa K, Yamamoto N, Kouzu Y, Koike H, Shiiba M, et al. Overexpression of Septin1: possible contribution to the development of oral cancer. Int J Oncol. (2007) 31:1021–8. doi: 10.3892/ijo.31.5.1021

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Tooley AJ, Gilden J, Jacobelli J, Beemiller P, Trimble WS, Kinoshita M, et al. Amoeboid T lymphocytes require the septin cytoskeleton for cortical integrity and persistent motility. Nat Cell Biol. (2009) 11:17–26. doi: 10.1038/ncb1808

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Yang J, Yin S, Bi F, Liu L, Qin T, Wang H, et al. TIMAP repression by TGFbeta and HDAC3-associated Smad signaling regulates macrophage M2 phenotypic phagocytosis. J Mol Med. (2017) 95:273–85. doi: 10.1007/s00109-016-1479-z

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Ghigo E, Barry AO, Pretat L, Al Moussawi K, Desnues B, Capo C, et al. IL-16 promotes T. whipplei replication by inhibiting phagosome conversion and modulating macrophage activation. PLoS ONE. (2010) 5:e13561. doi: 10.1371/journal.pone.0013561

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Richmond J, Tuzova M, Cruikshank W, Center D. Regulation of cellular processes by interleukin-16 in homeostasis and cancer. J Cell Physiol. (2014) 229:139–47. doi: 10.1002/jcp.24441

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Garrido F, Ruiz-Cabello F, Aptsiauri N. Rejection versus escape: the tumor MHC dilemma. Cancer Immunol Immunother. (2017) 66:259–71. doi: 10.1007/s00262-016-1947-x

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Terry S, Savagner P, Ortiz-Cuaran S, Mahjoubi L, Saintigny P, Thiery JP, et al. New insights into the role of EMT in tumor immune escape. Mol Oncol. (2017) 11:824–46. doi: 10.1002/1878-0261.12093

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Datar I, Schalper KA. Epithelial-mesenchymal transition and immune evasion during lung cancer progression: the chicken or the egg? Clin Cancer Res. (2016) 22:3422–4. doi: 10.1158/1078-0432.CCR-16-0336

CrossRef Full Text | Google Scholar