Gemcitabine-Resistant Biomarkers in Bladder Cancer are Associated with Tumor-Immune Microenvironment

Introduction

Bladder cancer (BCa) is the ninth most prevalent cancer globally (Antoni et al., 2017; Bray et al., 2018). Approximately 540,000 new cases have been diagnosed and 195,000 patients die of BCa each year (Antoni et al., 2017; Bray et al., 2018). The BCa incidence is elevated worldwide and the tumor burden increases due to population aging and environmental pollution during the past 2 decades (Ebrahimi et al., 2019; Yang et al., 2019).

Although surgical operation has been utilized in BCa treatment, the prognosis is still poor (Kirkali et al., 2005; Goebell et al., 2008; Humphrey et al., 2016). The 5-year relapse rate after initial treatment is from 55% to 85% in non-muscle-invasive BCa (NMIBC) and the 5-year survival rate is from 30% to 45% in muscle-invasive BCa (MIBC) (Choueiri and Raghavan, 2008; Lightfoot et al., 2014; Antoni et al., 2017). Chemotherapy is a promising treatment for reducing the recurrence rate and improving the survival rate of BCa patients (Coen et al., 2019). Gemcitabine (GEM) is a kind of cytosine analogue that inhibits DNA synthesis. Combination therapy of GEM and other chemotherapeutic drugs has been widely utilized in the treatment of MIBC (Oh et al., 2009).

However, GEM-resistance causes a severe challenge in the treatment of MIBC. It is reported that the response rate of advanced MIBC with GEM treatment is less than 40%, which indicates a limited efficacy of GEM treatment (Kaufman et al., 2009; Sternberg et al., 2013). The long-term curative effects of GEM declined sharply with the extension of treatment time (Cao et al., 2018). In addition, inherent or acquired drug resistance is usually observed in clinical practice (Bergman et al., 2002; Wang and Lippard, 2005). As a consequence, it is necessary and vital to explore potential mechanisms of GEM-resistance. On the one hand, MIBC has the nature of high somatic-mutation frequency and molecular heterogeneity, which exerts a critical role in drug resistance (Giannopoulou et al., 2019). On the other hand, dysfunction of immune system exerts a crucial role in tumor resistance (Yu et al., 2018). Co-delivery of GEM and small interfering RNA targeting IDO1 could relieve the immune brakes and further alleviate the immune inhibition of M2 macrophages, which indicates that these immune cells are associated with regulation of immune response to GEM (Chen et al., 2020). In addition, bioinformatics analyses constructs a microRNA (miRNA)-gene regulatory network associated with alteration of memory CD4⁺ T cells in GEM-resistant pancreatic cancer cells, which suggests the immune system is implicated in the microenvironment of GEM-resistance (Gu et al., 2020).

The present study focused on the key genes, miRNA-gene regulatory network, and their immune microenvironment based on GEM-resistant BCa in order to explore reliable prognostic indicators and provide treatment targets for GEM-resistant BCa.

Materials and methods

The Cancer Genome Atlas (TCGA) Bladder Urothelial Carcinoma (BLCA) dataset

The gene expression profiling dataset and clinical data of TCGA BLCA (accessed September 1, 2021) were downloaded from the TCGA website (http://portal.gdc.cancer.gov/). TCGA BLCA comprised BCa tissues (n = 404) and adjacent normal tissues (n = 18). Among all samples, there were 18 pairs of BCa tissues (n = 18) and matched adjacent normal tissues (n = 18).

GSE77883 dataset from Gene Expression Omnibus (GEO) database

Gene expression profiles of GSE77883 (accessed September 1st, 2021) were downloaded from GEO (http://www.ncbi.nlm.nih.gov/geo/). Six cells containing untreated T24 cells (n = 3) and GEM-resistant T24 cells (n = 3) were enrolled. RNA was extracted and measured through microarray (Platform: GPL17077 Agilent-039494 SurePrint G3 Human GE v2 8 × 60 K Microarray 039381).

Gene Set Variation Analysis (GSVA)

We downloaded the gene set of TOOKER_GEMCITABINE_RESISTANCE_UP (M19654) (Tooker et al., 2007) from the Molecular Signatures Database (MSigDB: http://www.gsea-msigdb.org/gsea/index.jsp) (Liberzon et al., 2015), and M19654 is the key GEM-related gene set from chemical and genetic perturbations of MSigDB. The normalized GEM-resistance GSVA score of the gene set was measured for each BCa tissue from TCGA BLCA using the GSVA algorithm with the GSVA R package (Hänzelmann et al., 2013). The median value of the GEM-resistance score was used to divide all TCGA BCa tissues into high score of the GEM-resistance group (n = 202) and low score of the GEM-resistance group (n = 202).

Immune cells analysis in tumor-immune microenvironment

Tumor-infiltrating immune cells were measured and analyzed with the MCPcounter (Microenvironment Cell Populations-counter) R package (Becht et al., 2016). In order to explore whether GEM-resistance is immune-related, we compared the immune cells between the high score of the GEM-resistance group (n = 202) and low score of the GEM-resistance group (n = 202) in TCGA BCa tissues. The Pearson method was adopted to assess the correlation between the normalized GEM-resistance GSVA score and immune cells. To further clarify the correlation between immune system and BCa development, we compared the immune cells between 18 pairs of BCa tissues (n = 18) and matched adjacent normal tissues (n = 18) in TCGA BLCA.

Identification of differentially expressed genes (DEGs) in GSE77883 and TCGA BLCA datasets

Limma R package was adopted to detect DEGs between GEM-resistant T24 cells and untreated T24 cells in the GSE77883 dataset (Ritchie et al., 2015). In addition, we also identified DEGs between 18 pairs of BCa tissues and adjacent tissues in TCGA BLCA. Benjamini-Hochberg method (Benjamini and Hochberg, 1995; Benjamini and Hochberg, 2000) was used to adjust the p-values for multiplicity and control false discovery rate. The thresholds were |logFC| ≥ 1 and adjusted p-value (adj. p-value) < 0.05. Venn diagram was adopted to find overlapped DEGs in the GSE77883 and TCGA BLCA datasets. These overlapped DEGs were considered as GEM-resistant genes in BCa.

Functional enrichment analysis

ClusterProfiler R package (Yu et al., 2012) was applied to identify the biological functions of overlapped DEGs through Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway collections. Metascape (http://metascape.org) was utilized to identify the most closely enriched clusters (Zhou et al., 2019).

Protein–protein interaction (PPI) network and selection of hub genes

We applied String (version 11.0: http://string-db.org/) (Szklarczyk et al., 2017) to construct interactions among proteins on the basis of the overlapped DEGs with the interaction score of 0.400 was set as threshold. In addition, cytoHubba in Cytoscape software screened 10 genes with highest connection degrees. Univariate Cox regression analysis further detected prognosis-related genes among the top 10 genes based on BCa patients (n = 404) from the TCGA BLCA dataset.

Immunohistochemistry and validation by TCGA BLCA and Genotype-Tissue Expression (GTEx) datasets

Immunohistochemistry was extracted and analyzed from The Human Protein Atlas (THPA) database (http://www.proteinatlas.org/) (Uhlen et al., 2010). We evaluated expression levels of the identified six prognosis-related genes (CAV1, CSPG4, FBLN1, COL6A2, FABP4, and PCOLCE) between tumor and normal tissues at protein level.

To confirm the differential expression of the six prognosis-related genes in a larger sample size, box plots were adopted to compare the expression of CAV1, CSPG4, FBLN1, COL6A2, FABP4, and PCOLCE between BCa tissues and normal tissues from the TCGA BLCA dataset and the GTEx dataset through Gene Expression Profiling Interactive Analysis (GEPIA) (http://gepia.cancer-pku.cn/) (Tang et al., 2017).

Survival analysis

Based on the TCGA BLCA dataset, we analyzed the six selected genes (CAV1, CSPG4, FBLN1, COL6A2, FABP4, and PCOLCE) with overall survival (OS) and disease-free survival (DFS) through the Kaplan–Meier (KM) Plotter (http://kmplot.com/analysis/) (Nagy et al., 2018).

Predictive value of the six hub genes in immunotherapy

CAMOIP (Comprehensive Analysis on Multi-Omics of Immunotherapy in Pan-cancer) is a tool for analyzing the expression data and mutation data from the immunotherapy-treated projects, using a standard processing pipeline (Lin et al., 2021). The IMvigor210 cohort to investigate the clinical activity of immunotherapy with atezolizumab in metastatic BCa was used for an integrated biomarker evaluation (Mariathasan et al., 2018). We used gene expression profiling from the IMvigor210 cohort to evaluate the predictive value of the six key genes (CAV1, COL6A2, FABP4, FBLN1, PCOLCE, and CSPG4) in OS after immunotherapy through CAMOIP.

Pearson correlation analysis explored the six hub genes with tumor-infiltrating immune cells

TIMER (Tumor Immune Estimation Resource) (http://timer.cistrome.org/) (Li et al., 2016; Li et al., 2017) was adopted to measure the impacts of immune cells on BCa OS through separating all BCa samples (n = 404) into high and low abundance groups based on median value of each immune cell abundance. In addition, Pearson method measured correlations between immune cells and expression levels of the six genes (CAV1, CSPG4, FBLN1, COL6A2, FABP4, and PCOLCE) in BCa patients (n = 404) from the TCGA BLCA dataset through Pearson correlation analysis.

Genetic mutation analysis

The cBioPortal database (http://www.cbioportal.org/) was utilized to investigate mutations of the six hub genes (CAV1, CSPG4, FBLN1, COL6A2, FABP4, and PCOLCE) in BCa patients (n = 404) from the TCGA BLCA dataset. Survival analysis and immune cells were also explored based on genetic mutations.

miRNA-Gene regulatory network

Two databases including miRTarbase (http://mirtarbase.mbc.nctu.edu.tw/php/) (Chou et al., 2018) and Targetscan (http://www.targetscan.org/vert_72/) (Agarwal et al., 2015) were experimentally validated miRNA-target gene interaction databases and they were adopted to predict upstream miRNAs and to build the miRNA-gene regulatory network. Venn diagram was used to identify overlapped miRNAs as key miRNAs and KM Plotter was utilized to evaluate the effects of key miRNAs on BCa OS based on BCa patients from the TCGA BLCA dataset as mentioned before. Pearson method was performed to evaluate the pairwise gene correlation in BCa samples from the TCGA BLCA dataset.

Gene set enrichment analysis (GSEA)

To clarify the roles of the six hub genes (CAV1, CSPG4, FBLN1, COL6A2, FABP4, and PCOLCE), we applied GSEA to analyze the enrichment of BCa samples (n = 404) in the TCGA BLCA dataset by assessing the normalized enrichment score (NES) (Subramanian et al., 2007).

Screening of potential targeted drugs in GEM-resistant BCa

The Connectivity Map (CMAP) database (http://portals.broadinstitute.org/cmap) (Lamb et al., 2006; Lamb, 2007) was utilized to explore the potential drugs with antagonism or synergism to GEM-resistance. This database used details of DEGs to identify potential targeted drugs. Eighty-two key DEGs in GEM-resistant BCa were uploaded to the CMAP database. Next, the 82 key DEGs were compared to expression profiles stored in CMAP in order to select potential drugs. Drugs with p-value <0.05 were considered as significant targeted drugs to GEM-resistance.

Statistical analysis

Statistical analyses were performed by R software (v3.6.1: http://www.r-project.org), GraphPad Prism 7.0, Metascape, GEPIA, KM Plotter, TIMER, and cBioPortal. Univariate Cox regression, KM method, and log-rank test were adopted for survival analysis by calculating hazard ratio (HR) and 95% confidence interval (CI). Student’s t test was applied to evaluate quantitative variables. The correlation coefficient, R-value, was used to estimate the strength of Pearson correlation analysis. The two-sided p-value <0.05 was set as the threshold. Cytoscape software (v.3.6.1) was adopted for visualization of networks.

Results

Figure 1 showed the workflow.

FIGURE 1

FIGURE 1. Workflow of this study. TCGA: The Cancer Genome Atlas; BLCA: Bladder urothelial carcinoma; GEO: Gene Expression Omnibus; BCa: Bladder cancer; GSVA: Gene set variation analysis; GTEx: Genotype-tissue expression; DEGs: Differentially expressed genes; PPI: Protein–protein interaction; miRNA: microRNA; GSEA: Gene set enrichment analysis.

Tumor immune microenvironment in GEM-resistance

Based on 404 BCa tissues from the TCGA BLCA dataset, GSVA analysis indicated that BCa patients with high score of GEM-resistance had worse OS compared with those with low score (HR = 1.57, 95%CI = 1.14–2.16) (p = 0.006) (Figure 2A). Furthermore, correlation analysis revealed GEM-resistance score was positively associated with cytotoxicity scores (R = 0.376, p < 0.001), macrophages/monocytes (R = 0.317, p < 0.001), NK cells, and cancer-associated fibroblasts. GEM-resistance score was negatively associated with myeloid dendritic cells (Figures 2B,C). In addition, BCa patients with high score of GEM-resistance had higher abundance of the cytotoxicity scores, macrophages/monocytes, NK cells, and cancer-associated fibroblasts as well as lower abundance of myeloid dendritic cells in comparison with those with low score, which confirmed the results of correlation analysis (Figures 2D,E). From the above results, we identified that GEM-resistance score in BCa is closely related to immune microenvironment. In addition, GSVA analysis based on the GSE77883 dataset indicated that GEM-resistant T24 cells had higher score of GEM-resistance than untreated T24 cells (p = 0.03) (Supplementary Figure S1).

FIGURE 2

FIGURE 2. Gene set variation analysis identified that gemcitabine (GEM)-resistance was associated with prognosis and immune microenvironment in 404 bladder cancer (BCa) patients from the TCGA BLCA dataset. (A) Kaplan–Meier survival indicated BCa patients with high score of GEM-resistance had poor overall survival; (B,C) Correlations between GEM-resistance score and immune cells; (D,E) Differences in abundance of immune cells between high score and low score of GEM-resistance. *p < 0.05.

Tumor-immune microenvironment in BCa development

Based on 18 pairs of BCa tissues and matched adjacent normal tissues, immune cells infiltration analysis suggested that the abundance of B cells, myeloid dendritic cells, endothelial cells, and cancer associated fibroblast cells was obviously down-regulated in BCa tissues compared with matched adjacent normal tissues (p < 0.05). However, for other immune cells, no change was observed (p > 0.05) (Supplementary Figure S2). From the above results, we identified that BCa development is closely related to an immune microenvironment.

Identification of key genes in GSE77883 and TCGA BLCA datasets

We firstly compared GEM-resistant T24 cells with untreated T24 cells in the GSE77883 dataset; 2,176 DEGs containing 888 up-regulated and 1,289 down-regulated genes were obtained in GEM-resistant BCa cells (Figure 3A and Figure 3C).

FIGURE 3

FIGURE 3. Identification of key genes in the GSE77883 and TCGA BLCA datasets. (A) Heat map of differentially expressed genes (DEGs) between gemcitabine (GEM)-resistant T24 cells and untreated T24 cells based on the GSE77883 dataset; (B) Heat map of DEGs between bladder cancer (BCa) tissues and matched adjacent normal tissues based on the TCGA BLCA dataset; (C) Volcano plot of DEGs between GEM-resistant T24 cells and untreated T24 cells based on the GSE77883 dataset; (D) Volcano plot of DEGs between BCa tissues and matched adjacent normal tissues based on the TCGA BLCA dataset; (E) Venn diagram identified overlapped DEGs in both GSE77883 and TCGA BLCA datasets; (F) Heat map of 82 overlapped DEGs and they were key genes in both GEM-resistance and BCa development. adj.P.-value was adjusted p-value for Benjamini-Hochberg (BH) method.

In addition, we compared 18 pairs of BCa tissues and matched adjacent normal tissues in the TCGA BLCA dataset. We obtained 1,398 DEGs containing 468 up-regulated and 930 down-regulated genes in BCa tissues (Figure 3B and Figure 3D).

In order to investigate which genes were associated with both GEM-resistance and BCa development, Venn diagram combined the DEGs from two datasets and identified the overlapped DEGs obtained in both datasets. Ultimately, 82 overlapped DEGs (11 overlapped up-regulated and 71 overlapped down-regulated genes) were obtained and were considered as GEM-resistant genes in BCa (Figures 3E,F) (Supplementary Table S1).

Functional enrichment analysis

Biological process analysis indicated that the 82 overlapped DEGs were enriched in regulation of the epithelial cell apoptotic process (GO:1904035) and muscle tissue development (GO:0060537). Component analysis detected that they were mainly located at the extracellular matrix (GO:0031012) and cell leading edge (GO:0031252). Molecular function analysis demonstrated that they participated in growth factor binding (GO:0019838) and toxic substance binding (GO:0015643) (Figure 4A and Table 1). Pathway analyses identified GEM-resistance was associated with the peroxisome proliferator-activated receptor (PPAR) signaling pathway (hsa03320), extracellular matrix (ECM)-receptor interaction (hsa04512), and focal adhesion (hsa04510) (Figures 4B,C and Table 1).

FIGURE 4

FIGURE 4. Enrichment analyses through 82 overlapped differentially expressed genes. (A) Gene ontology enrichment analysis. Biological process (BP) indicated they were enriched in regulation of epithelial cell apoptotic process (GO:1904035) and muscle tissue development (GO:0060537); cell component (CC) indicated they were enriched in extracellular matrix (GO:0031012) and cell leading edge (GO:0031252); molecular function (MF) indicated they were enriched in growth factor binding (GO:0019838) and toxic substance binding (GO:0015643); (B,C) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway and gene-concept network analysis. They were enriched in the peroxisome proliferator-activated receptor (PPAR) signaling pathway (hsa03320), extracellular matrix (ECM)-receptor interaction (hsa04512), and focal adhesion (hsa04510).

TABLE 1

TABLE 1. Functional enrichment analysis results

Metascape identified the interactions of the main 19 clustered enrichment terms (Supplementary Figure S3 and Supplementary Table S2). Table 2 shows that negative regulation of cell proliferation (GO:0008285) and its relevant enrichment terms were significantly associated with immune cells including endothelial cells (GO:2000351), T cells (GO:0050870), and granulocytes (GO:0071621), which confirmed the close relationship between GEM-resistance and immune system.

TABLE 2

TABLE 2. Immune-related enrichment terms associated with immune cells proliferation

PPI network and selection of hub genes

Finally, PPI network enrolled 37 nodes and 49 edges (Figure 5A). Among the interactions,CAV1, COL6A2, FABP4, FBLN1, PCOLCE, CSPG4, CCL2, THBS1, CFH, and COL6A1 with the highest degree scores were considered as the top 10 genes. Figure 5B showed the key module constituted by the 10 genes.

FIGURE 5

FIGURE 5. Protein–protein interaction (PPI) network and selection of hub genes. (A) PPI network of DEGs; (B) Hub module; (C) Six prognostic genes (CAV1, COL6A2, FABP4, FBLN1, PCOLCE, and CSPG4) through univariate Cox regression. Bold genes meant prognosis-related genes.

Ultimately, univariate Cox regression analysis detected that six genes (CAV1, COL6A2, FABP4, FBLN1, PCOLCE, and CSPG4) were identified to be prognosis-related (p < 0.05) (Table 3) (Figure 5C).

TABLE 3

TABLE 3. The six hub genes with high degree scores

Immunohistochemistry and validation of hub genes by TCGA and GTEx datasets

We collected and analyzed immunohistochemistry of BCa tissues and normal bladder tissues from THPA. Expression levels of the six hub genes (CAV1, COL6A2, FABP4, FBLN1, PCOLCE, and CSPG4) were evaluated at protein level. Figure 6 reveals the six hub genes were down-regulated in BCa tissues.

FIGURE 6

FIGURE 6. Immunohistochemistry and validation of six hub genes by the TCGA BLCA and GTEx datasets. Immunohistochemistry from THPA (left part in each subfigure) and bladder cancer (BCa) tissues (n = 404) and normal tissues (n = 28) from the TCGA BLCA and GTEx datasets (right part in each subfigure) indicated that the six selected genes (CAV1, COL6A2, FABP4, FBLN1, PCOLCE, and CSPG4) were down-regulated in BCa. (A) CAV1; (B) CSPG4; (C) FBLN1; (D) COL6A2; (E) FABP4; (F) PCOLCE.

In addition, we also validated the six hub genes between 404 BCa tissues and 28 normal bladder tissues from the TCGA BLCA and GTEx datasets. Box plots indicated that they were all down-regulated in BCa samples (p < 0.05), which was in accord with the results of immunohistochemistry (Figure 6).

Survival analysis and predictive value of the six hub genes in immunotherapy

KM curves suggested that the six hub genes could influence the OS time (Supplementary Figure S4). Lower expression levels of FABP4, FBLN1, and COL6A2 indicated worse OS time (p < 0.05). However, higher expression levels of CSPG4, CAV1, and PCOLCE might be indicators for worse OS time (p < 0.05) for BCa patients. In addition, we identified that FABP4, CSPG4, COL6A2, and PCOLCE were related to DFS time of BCa (p < 0.05) (Supplementary Figure S4).

IMvigor210 cohort indicated that COL6A2 (HR > 1), FABP4 (HR > 1), and FBLN1 (HR < 1) could predict the OS after immunotherapy with atezolizumab (p < 0.05) (Figure 7).

FIGURE 7

FIGURE 7. IMvigor210 cohort indicated that COL6A2, FABP4 and FBLN1 could predict the OS after immunotherapy with atezolizumab. (A) FBLN1; (B) FABP4; (C) COL6A2; (D) CAV1; (E) PCOLCE; (F) CSPG4.

Tumor-infiltrating immune cells with hub genes in TCGA BLCA dataset

We used TIMER to clarify the effects of immune cells on OS of BCa patients in the TCGA BLCA dataset. We found that five immune cells (CD4⁺ T cells, CD8⁺ T cells, neutrophils, endothelial cells, and cancer associated fibroblast cells) were associated with OS of BCa. Patients with high abundance of CD4⁺ T cells and CD8⁺ T cells had longer OS time than those with low infiltration levels. However, patients with low infiltration levels of neutrophils, endothelial cells, and cancer-associated fibroblast cells had longer OS time compared with those with high infiltration levels of these immune cells (Supplementary Figure S5).

Pearson correlation analysis evaluated the correlations between the six hub genes (CAV1, COL6A2, FABP4, FBLN1, PCOLCE, and CSPG4) and abundance of main immune cells (Table 4 and Figure 8A) (Supplementary Figure S6). When we restricted the robust R-value to more than 0.400, we found CAV1 had strong correlation with dendritic cells. COL6A2 and PCOLCE had strong correlations with macrophages (Figure 8B).

TABLE 4

TABLE 4. Pearson correlation analysis indicated the six hub genes (CAV1, COL6A2, FABP4, FBLN1, PCOLCE and CSPG4) were associated with immune cells infiltration

FIGURE 8

FIGURE 8. Tumor-infiltrating immune cell analysis with the six hub genes (CAV1, COL6A2, FABP4, FBLN1, PCOLCE, and CSPG4) through Pearson correlation method. (A) Heat map showed correlations between immune cells and six hub genes; (B) Robust correlations (R-value more than 0.400) was identified between CAV1, COL6A2, and PCOLCE and immune cells.

In addition, GSVA and Pearson correlation analysis identified PCOLCE, CSPG4, COL6A2, and CAV1 were positively correlated with the score of GEM-resistance (R-value >0, p < 0.05) (Figures 9A–D). However, FBLN1 and FABP4 were negatively correlated with GEM-resistance (R-value <0, p < 0.05) (Figures 9E,F).

FIGURE 9

FIGURE 9. Genetic mutations analysis of the six hub genes (CAV1, COL6A2, FABP4, FBLN1, PCOLCE, and CSPG4) in the TCGA BLCA dataset. (A–F) Correlations between GEM-resistance score and the six hub genes; (G) Mutation frequencies of CAV1, COL6A2, FABP4, FBLN1, PCOLCE, and CSPG4 in the TCGA BLCA dataset; (H) Kaplan–Meier survival curves showed that genetic mutations of six selected genes (CAV1, COL6A2, FABP4, FBLN1, PCOLCE, and CSPG4) were not associated with overall survival (OS) based on the TCGA BLCA dataset; (I) CSPG4 mutation was associated with infiltration levels of CD4⁺ T cells and natural killer (NK) cells.

Genetic mutations of hub genes in TCGA BLCA dataset

Figure 9G illustrated the mutation frequencies of the six hub genes (CAV1, COL6A2, FABP4, FBLN1, PCOLCE, and CSPG4) in 404 BCa patients from the TCGA BLCA dataset. Among the six hub genes, the top three most frequently mutated genes were FBLN1 (5.0%), FABP4 (4.0%), and CSPG4 (4.0%). FBLN1 mutations included fusion mutation, amplification mutation, and deletion mutation. Most of FABP4 mutations were amplification mutations and most of CSPG4 mutations were missense mutations.

Furthermore, mutations of all the six hub genes didn’t influence the OS time in comparison with the wild-type group (p > 0.05) (Figure 9H). TIMER identified that CSPG4 mutations could elevate the abundance of CD4⁺ T cells (p = 0.009) and NK cells (p = 0.021) (Figure 9I), while the genetic mutations of the other five hub genes were not associated with abundance of immune cells.

Construction of miRNA-gene regulatory network

A total of 72 miRNAs might regulate the expression levels of six hub genes through the miRTarBase and Targetscan databases (Table 5). The miRNA-gene regulatory network is displayed in Figure 10.

TABLE 5

TABLE 5. Upstream microRNAs (miRNAs) of the six hub genes (CAV1, COL6A2, FABP4, FBLN1, PCOLCE and CSPG4)

FIGURE 10

FIGURE 10. microRNA-gene regulatory network.

As we could see in Figure 11A, hsa-miR-124-3p was the overlapped upstream miRNA of CAV1, COL6A2, and CSPG4; hsa-miR-26b-5p and hsa-miR-192-5p were overlapped upstream miRNAs of CAV1 and PCOLCE, which indicated the three miRNAs might be key miRNAs in regulating the hub genes. Survival analysis demonstrated that higher expression levels of hsa-miR-124-3p and hsa-miR-192-5p were significantly related to better prognosis in BCa patients from the TCGA BLCA dataset.

FIGURE 11

FIGURE 11. Identification of hub microRNAs (miRNAs) and the key miRNA-gene regulatory network. (A) hsa-miR-124-3p, hsa-miR-26b-5p, and hsa-miR-192-5p were overlapped upstream miRNAs of four hub genes (CAV1, COL6A2, PCOLCE, and CSPG4) and their survival analyses; (B) Pairwise correlation analyses of four hub genes (CAV1, COL6A2, PCOLCE, and CSPG4); (C) Key miRNA-gene regulatory network formed by the three overlapped miRNAs and the four hub genes.

Furthermore, a pairwise correlation analysis of CAV1, COL6A2, CSPG4, and PCOLCE was carried out. Expression of the four genes demonstrated a significant positive correlation (p < 0.05). The increase of one hub gene was strongly correlated with the increase of another one (Figure 11B). Hence, based on the overlapped upstream miRNAs and correlation analysis, we confirmed that the three key miRNAs (hsa-miR-124-3p, hsa-miR-26b-5p, and hsa-miR-192-5p) could regulate the four hub genes (CAV1, COL6A2, CSPG4, and PCOLCE) in GEM-resistance and tumor immune microenvironment. We also plotted the key regulatory network in Figure 11C. In addition, we detected DEGs between high score of GEM-resistance and low score of GEM-resistance. We found PCOLCE, CSPG4, COL6A2, and CAV1 were up-regulated in patients with high GEM-resistance scores, which further confirmed their key roles in GEM-resistance (Supplementary Figure S7).

Gene set enrichment analysis

GSEA was conducted to investigate the possible role of the six hub genes (CAV1, COL6A2, FABP4, FBLN1, PCOLCE, and CSPG4) involved in GEM-resistance. We identified CSPG4 (Figure 12) was obviously enriched in cancer-related pathways and functions including the bladder cancer pathway (hsa05219) and TGF-β signaling pathway (hsa04350). In addition, CSPG4 exerted a vital role in chemotherapy-related functions including drug metabolism of cytochrome P450 (hsa00982), drug metabolism of other enzymes (hsa00983), cellular response to drug (GO:0035690), and response to drug (GO:0042493). Further, CSPG4 was also enriched in pathways of immune response and immune cells including the B cell receptor signaling pathway (hsa04662), NK cell-mediated cytotoxicity pathway (hsa04650), and T cell receptor signaling pathway (hsa04660) (Supplementary Tables S3, S4). In addition, the other five genes were also enriched in cancer-related (hsa05200 and hsa05219), immune-related (hsa04660, hsa4650, and hsa04662), and chemotherapy-related (hsa00982 and hsa00983) pathways (Supplementary Table S5).

FIGURE 12

FIGURE 12. Gene set enrichment analysis of CSPG4. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis and Gene Ontology (GO) functional analysis identified that CSPG4 was enriched in cancer-related, chemotherapy-related, and immune-related functions. (A) Cancer-related KEGG pathway analysis; (B) Chemotherapy-related KEGG pathway analysis; (C) Immune-related KEGG pathway analysis; (D) Cancer-related GO enrichment analysis; (E) Chemotherapy-related GO enrichment analysis; (F) Immune-related GO enrichment analysis.

Selection of targeted drugs for GEM-resistant BCa

CMAP analysis indicated 75 drugs might have antagonistic or synergistic effects on GEM-resistance. According to the enrichment score, Table 6 displayed the top 10 drugs with antagonism and the top 10 drugs with synergism, respectively. The top 10 antagonistic drugs for GEM-resistant BCa were lisinopril, rifabutin, clonidine, prasterone, vorinostat, prednisone, nifenazone, alvespimycin, trichostatin A, and tanespimycin.

TABLE 6

TABLE 6. Connectivity map (CMAP) database analysis

Discussion

Accumulating evidence indicates that aberrantly expressed genes are significantly associated with GEM-resistance in BCa. It is reported that CSNK1D played a key role in the metabolism of GEM, and inhibition of CSNK1D could sensitize BCa cells to GEM treatment, which might be utilized as a therapeutic target for metastatic BCa (Udhaya Kumar et al., 2020; Vena et al., 2020). Xie et al. identified that circular RNA (circRNA) circHIPK3 was an independent prognostic predictor, and up-regulation of circHIPK3 promoted GEM-sensitivity in BCa (Xie et al., 2020). Another study based on BCa cells indicated that inhibition of GP130 could enhance the sensitivity to GEM and reduce viability and migration of tumor cells through regulating the PI3K/AKT/mTOR signaling pathways (Li X. et al., 2019). This study analyzed the GSE77883 dataset from the GEO and TCGA BLCA datasets to identify promising biomarkers for GEM-resistant BCa through high-throughput sequencing data and bioinformatics analyses.

We assessed immune cells and identified that both BCa development and GEM-resistance were immune-related. We found that 82 key DEGs were significantly related to both BCa development and GEM-resistance. Functional enrichment analyses found these key DEGs were enriched in immune-related items, especially in the regulation of immune cell proliferation. After construction of the PPI network and Cox regression analysis, we selected six hub genes (CAV1, COL6A2, FABP4, FBLN1, PCOLCE, and CSPG4) with the highest connectivity degrees and prognostic values for further analyses. We used immunohistochemistry from THPA and expression profiles from larger samples to confirm the down-regulation of six hub genes in BCa at protein and mRNA level, respectively. Survival analyses demonstrated that they were related to OS time. Down-regulation of CSPG4, CAV1, and PCOLCE might be related to elevated chemotherapy sensitivity and thus lower expression levels of them were associated with better OS. IMvigor210 cohort validated that COL6A2, FABP4, and FBLN1 could predict the OS after immunotherapy with atezolizumab. Pearson correlation analysis revealed CAV1, COL6A2, and PCOLCE had strong correlations with immune cells, such as dendritic cells and macrophages. Next, we constructed the key miRNA-gene regulatory network based on four key genes (CAV1, COL6A2, PCOLCE, and CSPG4) and three key miRNAs (hsa-miR-124-3p, hsa-miR-26b-5p, and hsa-miR-192-5p).

CAV1 is the chief component of the caveolae plasma membranes in most human cells and participates in immune response and cancer progression (Shi et al., 2020). CAV1 in prostate cancer could induce epithelial–mesenchymal transition through activating cancer immune evasion, and CAV1 in cancer-derived exosomes was able to induce chemoresistance in recipient cells (Lin et al., 2019), which conformed to our results revealed in GEM-resistant BCa. In addition, CAV1 might also promote systemic lupus erythematosus through regulating pathways of T cell costimulation, lymphocyte costimulation, and B cell receptor signaling (Udhaya Kumar et al., 2020). CAV1 was pivotal in acute immune-mediated hepatic damage through driving RNS-mediated ferroptosis (Deng et al., 2020). The presently found CAV1 was the key gene regulated by three upstream miRNAs in the miRNA-gene regulatory network. Zhou et al. found that hsa-miR-124-3p and hsa-miR-192-5p suppressed the proliferation and invasion of tumor cells by targeting CAV1 (Zhou et al., 2018; Chen et al., 2019). Therefore, we hypothesized that CAV1 could facilitate tumor development and GEM-resistance via immune escape mechanism.

COL6A2 is one of the collagen family members and encodes one alpha chain of type VI collagen identified in most connective tissues (Hou et al., 2016). COL6A2 was reported to be up-regulated and to gather in the ECM-receptor interaction signaling pathway, which promoted the BCa progression (Zhu et al., 2019). Down-regulation of COL6A2 induced by decreased IDO1 could suppress host anti-tumor immune response through inhibiting immune-related pathways (Xiang et al., 2019). We found COL6A2 was positively correlated with most immune cells in a tumor-immune microenvironment, which could support the highly immunogenic nature of COL6A2 in BCa.

CSPG4 is a kind of transmembrane proteoglycan, considered as a promising tumor-associated antigen (Cavallo et al., 2007). Previous investigations have identified CSPG4 as a key gene in soft-tissue sarcoma, melanoma, and glioblastoma (Benassi et al., 2009; Wang et al., 2011). We found CSPG4 exerted a vital role in BCa prognosis, and both the expression and mutation of CSPG4 might influence the immune cells in BCa. In addition, Rolih et al. systematically summarized the evidence of CSPG4 in tumor biology and suggested that CSPG4 and anti-CSPG4 vaccination strategy had the potential to be an attractive target for anti-tumor immunotherapy (Rolih et al., 2017). GSEA detected that CSPG4 contributed to cancer-related pathways, immune system process, and drug metabolism, which further confirmed its value in drug-resistance and immunotherapy of BCa (Song et al., 2020). Nevertheless, further investigations are demanded to verify its mechanisms in BCa.

PCOLCE is a glycoprotein that elevates the activity of procollagen C-proteinase (Pulido et al., 2018). PCOLCE was up-regulated in osteosarcoma and promotes the distant metastasis (Wang et al., 2019). We found PCOLCE was down-regulated and was also associated with prognosis in BCa. The difference of PCOLCE expression between the two tumors could be explained by the miRNA-gene regulatory network; hsa-miR-26b-5p and hsa-miR-192-5p regulated the co-expression of PCOLCE and CAV1 in the network of BCa. Besides, PCOLCE was revealed to be involved in platelet and endothelial function and immune activation in human immunodeficiency virus (HIV) patients after pitavastatin treatment, which bound PCOLCE to immune related functions (DeFilippi et al., 2020).

The identified key miRNA-gene regulatory network indicated that hsa-miR-124-3p, hsa-miR-192-5p, and hsa-miR-26b-5p were key miRNAs in regulating the above genes. It is reported that higher expression level of hsa-miR-124-3p suppressed tumor proliferation and indicated better BCa prognosis in BCa through targeting downstream genes (Wang et al., 2015; Zhou et al., 2018; Zo and Long, 2019). Besides, hsa-miR-124-3p was frequently and tumor-specifically methylated in primary BCa, indicating that epigenetic silencing of hsa-miR-124-3p may also participate in BCa development (Shimizu et al., 2013). In vitro studies demonstrated that hsa-miR-192-5p was a suppressor for BCa cells by cell cycle regulation and clinical studies identified hsa-miR-192-5p as an independent prognostic marker based on multivariate COX regression (Jin et al., 2015; Hu et al., 2017). Bioinformatics analysis combined with in vitro experiments demonstrated that hsa-miR-26b-5p was a critical regulator in BCa progression by targeting the proliferation-related gene PDCD10 (Wu et al., 2018). With regard to the relationship between the three miRNAs and immune system, there is evidence that they played roles in modulating immune escape and immune cells (Li K. et al., 2019; Ji et al., 2019; Amoruso et al., 2020).

The CMAP database analysis identified the potential drugs with synergism or antagonism to GEM-resistance. We identified that two histone deacetylase inhibitors (HDACIs), trichostatin A and vorinostat, had antagonistic effects on GEM-resistance, indicating that the two drugs might circumvent the GEM-resistance and enhance the sensitivity to GEM. In vivo studies based on BCa cells revealed that trichostatin A may synergistically enhance GEM-mediated cell cycle arrest and apoptosis through inhibiting the Raf/MEK/ERK pathway (Jeon et al., 2011; Lin et al., 2018), which provided HDACIs as promising treatment methods to improve GEM-resistant BCa patients in future clinical practice. It is also worth noticing that oxybuprocaine and benzocaine had the synergistic effects to GEM-resistance and the two anesthetics might aggravate GEM-resistance in BCa. Furthermore, ascorbic acid, also called vitamin C, was identified to be associated with aggravate GEM-resistance. However, a phase I clinical trial based on pancreatic cancer patients indicated that ascorbic acid combined concurrently with GEM was well tolerated and could reduce adverse events, which is inconsistent with our results (Welsh et al., 2013).

In order to reduce the potential bias caused by one single method and to verify the relationship between key biomarkers in GEM-resistant BCa and tumor immune microenvironment, the present study used two methods including MCPcounter (Becht et al., 2016) and TIMER (Li et al., 2016; Li et al., 2017) for the calculating of the immune cells. By using these two methods separately, the results indicated that GEM-resistance and the key genes were closely related to immune cells, which further verified that the results were stable and convincing. However, several limitations existed. Even if immunohistochemistry was used to confirm the down-regulation of hub genes, and cell lines instead of in vivo models were used to perform CMAP analysis, we didn’t verify their actual molecular mechanisms. Therefore, the above hypotheses should be verified through experimental methods in future studies. Since most BCa tissues from the TCGA BLCA dataset were MIBC tissues, our findings were more applicable to GEM-resistant MIBC patients.

Conclusion

We identified both BCa development and GEM-resistance were immune-related. CAV1, COL6A2, FABP4, FBLN1, PCOLCE, and CSPG4 are hub genes in GEM-resistant MIBC. They could serve as potential prognostic predictors and immunotherapy targets for MIBC. In addition, the key miRNA-gene regulatory network suggested three key miRNAs (hsa-miR-124-3p, hsa-miR-26b-5p, and hsa-miR-192-5p) might also be implicated in GEM-resistance. Ultimately, CMAP analysis identified HDACIs (trichostatin A and vorinostat) might circumvent the GEM-resistance and enhance the sensitivity to GEM.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author contributions

YS, TX, YD, CQ, and HL: conception and design. YS, WY, and YD: collection and analysis of data. YS, TX, JL, MD, and JH: interpretation and manuscript writing.

Funding

This work was supported by National Key Research and Development Program of China (2018YFA0902802), Natural Science Foundation of Beijing, China (7202219), National Natural Science Foundation of China (81802533), and Beijing Municipal Science and Technology Commission (Z191100006619010).

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/fcell.2021.809620/full#supplementary-material

References

Agarwal, V., Bell, G. W., Nam, J.-W., and Bartel, D. P. (2015). Predicting Effective microRNA Target Sites in Mammalian mRNAs. Elife 4, e05005. doi:10.7554/eLife.05005

PubMed Abstract | CrossRef Full Text | Google Scholar

Amoruso, A., Blonda, M., Gironi, M., Grasso, R., Di Francescantonio, V., Scaroni, F., et al. (2020). Immune and central Nervous System-Related miRNAs Expression Profiling in Monocytes of Multiple Sclerosis Patients. Sci. Rep. 10 (1), 6125. doi:10.1038/s41598-020-63282-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Antoni, S., Ferlay, J., Soerjomataram, I., Znaor, A., Jemal, A., and Bray, F. (2017). Bladder Cancer Incidence and Mortality: A Global Overview and Recent Trends. Eur. Urol. 71 (1), 96–108. doi:10.1016/j.eururo.2016.06.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Becht, E., Giraldo, N. A., Lacroix, L., Buttard, B., Elarouci, N., Petitprez, F., et al. (2016). Estimating the Population Abundance of Tissue-Infiltrating Immune and Stromal Cell Populations Using Gene Expression. Genome Biol. 17 (1), 218. doi:10.1186/s13059-016-1070-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Benassi, M. S., Pazzaglia, L., Chiechi, A., Alberghini, M., Conti, A., Cattaruzza, S., et al. (2009). NG2 Expression Predicts the Metastasis Formation in Soft-Tissue Sarcoma Patients. J. Orthop. Res. 27 (1), 135–140. doi:10.1002/jor.20694

CrossRef Full Text | Google Scholar

Benjamini, Y., and Hochberg, Y. (1995). Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. J. R. Stat. Soc. Ser. B (Methodological) 57 (1), 289–300. doi:10.1111/j.2517-6161.1995.tb02031.x

CrossRef Full Text | Google Scholar

Benjamini, Y., and Hochberg, Y. (2000). On the Adaptive Control of the False Discovery Rate in Multiple Testing with Independent Statistics. J. Educ. Behav. Stat. 25 (1), 60–83. doi:10.3102/10769986025001060

CrossRef Full Text | Google Scholar

Bergman, A. M., Pinedo, H. M., and Peters, G. J. (2002). Determinants of Resistance to 2′,2′-difluorodeoxycytidine (Gemcitabine). Drug Resist. Updates 5 (1), 19–33. doi:10.1016/s1368-7646(02)00002-x

CrossRef Full Text | Google Scholar

Bray, F., Ferlay, J., Soerjomataram, I., Siegel, R. L., Torre, L. A., and Jemal, A. (2018). Global Cancer Statistics 2018: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA: A Cancer J. Clinicians 68 (6), 394–424. doi:10.3322/caac.21492

CrossRef Full Text | Google Scholar

Cao, J., Wang, Q., Wu, G., Li, S., and Wang, Q. (2018). miR-129-5p Inhibits Gemcitabine Resistance and Promotes Cell Apoptosis of Bladder Cancer Cells by Targeting Wnt5a. Int. Urol. Nephrol. 50 (10), 1811–1819. doi:10.1007/s11255-018-1959-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Cavallo, F., Calogero, R. A., and Forni, G. (2007). Are Oncoantigens Suitable Targets for Anti-tumour Therapy? Nat. Rev. Cancer 7 (9), 707–713. doi:10.1038/nrc2208

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, C., Li, A., Sun, P., Xu, J., Du, W., Zhang, J., et al. (2020). Efficiently Restoring the Tumoricidal Immunity against Resistant Malignancies via an Immune Nanomodulator. J. Controlled Release 324, 574–585. doi:10.1016/j.jconrel.2020.05.039

CrossRef Full Text | Google Scholar

Chen, P., Feng, Y., Zhang, H., Shi, X., Li, B., Ju, W., et al. (2019). MicroRNA-192 I-nhibits C-ell P-roliferation and I-nduces A-poptosis in H-uman B-reast C-ancer by T-argeting C-aveolin 1. Oncol. Rep. 42 (5), 1667–1676. doi:10.3892/or.2019.7298

PubMed Abstract | CrossRef Full Text | Google Scholar

Chou, C.-H., Shrestha, S., Yang, C.-D., Chang, N.-W., Lin, Y.-L., Liao, K.-W., et al. (2018). miRTarBase Update 2018: a Resource for Experimentally Validated microRNA-Target Interactions. Nucleic Acids Res. 46 (D1), D296–D302. doi:10.1093/nar/gkx1067

PubMed Abstract | CrossRef Full Text | Google Scholar

Choueiri, T. K., and Raghavan, D. (2008). Chemotherapy for Muscle-Invasive Bladder Cancer Treated with Definitive Radiotherapy: Persisting Uncertainties. Nat. Rev. Clin. Oncol. 5 (8), 444–454. doi:10.1038/ncponc1159

CrossRef Full Text | Google Scholar

Coen, J. J., Zhang, P., Saylor, P. J., Lee, C. T., Wu, C.-L., Parker, W., et al. (2019). Bladder Preservation with Twice-A-Day Radiation Plus Fluorouracil/Cisplatin or once Daily Radiation Plus Gemcitabine for Muscle-Invasive Bladder Cancer: NRG/RTOG 0712-A Randomized Phase II Trial. Jco 37 (1), 44–51. doi:10.1200/JCO.18.00537

PubMed Abstract | CrossRef Full Text | Google Scholar

DeFilippi, C., Toribio, M., Wong, L. P., Sadreyev, R., Grundberg, I., Fitch, K. V., et al. (2020). Differential Plasma Protein Regulation and Statin Effects in Human Immunodeficiency Virus (HIV)-Infected and Non-HIV-infected Patients Utilizing a Proteomics Approach. J. Infect. Dis. 222 (6), 929–939. doi:10.1093/infdis/jiaa196

PubMed Abstract | CrossRef Full Text | Google Scholar

Deng, G., Li, Y., Ma, S., Gao, Z., Zeng, T., Chen, L., et al. (2020). Caveolin-1 Dictates Ferroptosis in the Execution of Acute Immune-Mediated Hepatic Damage by Attenuating Nitrogen Stress. Free Radic. Biol. Med. 148, 151–161. doi:10.1016/j.freeradbiomed.2019.12.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Ebrahimi, H., Amini, E., Pishgar, F., Moghaddam, S. S., Nabavizadeh, B., Rostamabadi, Y., et al. (2019). Global, Regional and National Burden of Bladder Cancer, 1990 to 2016: Results from the GBD Study 2016. J. Urol. 201 (5), 893–901. doi:10.1097/JU.0000000000000025

PubMed Abstract | CrossRef Full Text | Google Scholar

Giannopoulou, A., Velentzas, A., Konstantakou, E., Avgeris, M., Katarachia, S., Papandreou, N., et al. (2019). Revisiting Histone Deacetylases in Human Tumorigenesis: The Paradigm of Urothelial Bladder Cancer. Ijms 20 (6), 1291. doi:10.3390/ijms20061291

PubMed Abstract | CrossRef Full Text | Google Scholar

Goebell, P. J., Legal, W., Weiss, C., Fietkau, R., Wullich, B., and Krause, S. (2008). Multimodale Therapien zum Blasenerhalt bei High-grade-Blasentumoren. Urologe 47 (7), 838–845. doi:10.1007/s00120-008-1715-4

CrossRef Full Text | Google Scholar

Gu, J., Zhang, J., Huang, W., Tao, T., Huang, Y., Yang, L., et al. (2020). Activating miRNA-mRNA Network in Gemcitabine-Resistant Pancreatic Cancer Cell Associates with Alteration of Memory CD4+ T Cells. Ann. Transl Med. 8 (6), 279. doi:10.21037/atm.2020.03.53

PubMed Abstract | CrossRef Full Text | Google Scholar

Hänzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: Gene Set Variation Analysis for Microarray and RNA-Seq Data. BMC bioinformatics 14, 7. doi:10.1186/1471-2105-14-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Hou, T., Tong, C., Kazobinka, G., Zhang, W., Huang, X., Huang, Y., et al. (2016). Expression of COL6A1 Predicts Prognosis in Cervical Cancer Patients. Am. J. Transl Res. 8 (6), 2838–2844.

Google Scholar

Hu, Y., Cheng, C., Hong, Z., and Shi, Z. (2017). Independent Prognostic miRNAs for Bladder Urothelial Carcinoma. Oncol. Lett. 14 (3), 3001–3005. doi:10.3892/ol.2017.6471

PubMed Abstract | CrossRef Full Text | Google Scholar

Humphrey, P. A., Moch, H., Cubilla, A. L., Ulbright, T. M., and Reuter, V. E. (2016). The 2016 WHO Classification of Tumours of the Urinary System and Male Genital Organs-Part B: Prostate and Bladder Tumours. Eur. Urol. 70 (1), 106–119. doi:10.1016/j.eururo.2016.02.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Jeon, H. G., Yoon, C. Y., Yu, J. H., Park, M. J., Lee, J. E., Jeong, S. J., et al. (2011). Induction of Caspase Mediated Apoptosis and Down-Regulation of Nuclear Factor-Κb and Akt Signaling Are Involved in the Synergistic Antitumor Effect of Gemcitabine and the Histone Deacetylase Inhibitor Trichostatin A in Human Bladder Cancer Cells. J. Urol. 186 (5), 2084–2093. doi:10.1016/j.juro.2011.06.053

CrossRef Full Text | Google Scholar

Ji, C., Guo, X., Ren, J., Zu, Y., Li, W., and Zhang, Q. (2019). Transcriptomic Analysis of microRNAs-mRNAs Regulating Innate Immune Response of Zebrafish Larvae against Vibrio Parahaemolyticus Infection. Fish Shellfish Immunol. 91, 333–342. doi:10.1016/j.fsi.2019.05.050

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, Y., Lu, J., Wen, J., Shen, Y., and Wen, X. (2015). Regulation of Growth of Human Bladder Cancer by miR-192. Tumor Biol. 36 (5), 3791–3797. doi:10.1007/s13277-014-3020-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Kaufman, D. S., Shipley, W. U., and Feldman, A. S. (2009). Bladder Cancer. The Lancet 374 (9685), 239–249. doi:10.1016/S0140-6736(09)60491-8

CrossRef Full Text | Google Scholar

Kirkali, Z., Chan, T., Manoharan, M., Algaba, F., Busch, C., Cheng, L., et al. (2005). Bladder Cancer: Epidemiology, Staging and Grading, and Diagnosis. Urology 66 (6), 4–34. doi:10.1016/j.urology.2005.07.062

PubMed Abstract | CrossRef Full Text | Google Scholar

Lamb, J., Crawford, E. D., Peck, D., Modell, J. W., Blat, I. C., Wrobel, M. J., et al. (2006). The Connectivity Map: Using Gene-Expression Signatures to Connect Small Molecules, Genes, and Disease. Science 313 (5795), 1929–1935. doi:10.1126/science.1132939

PubMed Abstract | CrossRef Full Text | Google Scholar

Lamb, J. (2007). The Connectivity Map: a New Tool for Biomedical Research. Nat. Rev. Cancer 7 (1), 54–60. doi:10.1038/nrc2044

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, B., Severson, E., Pignon, J.-C., Zhao, H., Li, T., Novak, J., et al. (2016). Comprehensive Analyses of Tumor Immunity: Implications for Cancer Immunotherapy. Genome Biol. 17 (1), 174. doi:10.1186/s13059-016-1028-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, K., Chen, Y., Li, A., Tan, C., and Liu, X. (2019a). Exosomes Play Roles in Sequential Processes of Tumor Metastasis. Int. J. Cancer 144 (7), 1486–1495. doi:10.1002/ijc.31774

CrossRef Full Text | Google Scholar

Li, T., Fan, J., Wang, B., Traugh, N., Chen, Q., Liu, J. S., et al. (2017). TIMER: A Web Server for Comprehensive Analysis of Tumor-Infiltrating Immune Cells. Cancer Res. 77 (21), e108–e110. doi:10.1158/0008-5472.CAN-17-0307

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X., He, S., Tian, Y., Weiss, R. M., and Martin, D. T. (2019b). Synergistic Inhibition of GP130 and ERK Signaling Blocks Chemoresistant Bladder Cancer Cell Growth. Cell Signal. 63, 109381. doi:10.1016/j.cellsig.2019.109381

PubMed Abstract | CrossRef Full Text | Google Scholar

Liberzon, A., Birger, C., Thorvaldsdóttir, H., Ghandi, M., Mesirov, J. P., and Tamayo, P. (2015). The Molecular Signatures Database Hallmark Gene Set Collection. Cell Syst. 1 (6), 417–425. doi:10.1016/j.cels.2015.12.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Lightfoot, A. J., Breyer, B. N., Rosevear, H. M., Erickson, B. A., Konety, B. R., and O'Donnell, M. A. (2014). Multi-institutional Analysis of Sequential Intravesical Gemcitabine and Mitomycin C Chemotherapy for Non-muscle Invasive Bladder Cancer. Urol. Oncol. Semin. Original Invest. 32 (1), e15–35. doi:10.1016/j.urolonc.2013.01.009

CrossRef Full Text | Google Scholar

Lin, A., Wei, T., Liang, J., Qi, C., Li, M., Luo, P., et al. (2021). CAMOIP: A Web Server for Comprehensive Analysis on Multi-Omics of Immunotherapy in Pan-Cancer. bioRxiv. doi:10.1101/2021.09.10.459722

CrossRef Full Text | Google Scholar

Lin, C.-J., Yun, E.-J., Lo, U.-G., Tai, Y.-L., Deng, S., Hernandez, E., et al. (2019). The Paracrine Induction of Prostate Cancer Progression by Caveolin-1. Cell Death Dis 10 (11), 834. doi:10.1038/s41419-019-2066-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, W.-C., Hsu, F.-S., Kuo, K.-L., Liu, S.-H., Shun, C.-T., Shi, C.-S., et al. (2018). Trichostatin A, a Histone Deacetylase Inhibitor, Induces Synergistic Cytotoxicity with Chemotherapy via Suppression of Raf/MEK/ERK Pathway in Urothelial Carcinoma. J. Mol. Med. 96 (12), 1307–1318. doi:10.1007/s00109-018-1697-7

CrossRef Full Text | Google Scholar

Mariathasan, S., Turley, S. J., Nickles, D., Castiglioni, A., Yuen, K., Wang, Y., et al. (2018). TGFβ Attenuates Tumour Response to PD-L1 Blockade by Contributing to Exclusion of T Cells. Nature 554 (7693), 544–548. doi:10.1038/nature25501

PubMed Abstract | CrossRef Full Text | Google Scholar

Nagy, Á., Lánczky, A., Menyhárt, O., and Győrffy, B. (2018). Validation of miRNA Prognostic Power in Hepatocellular Carcinoma Using Expression Data of Independent Datasets. Sci. Rep. 8 (1), 9227. doi:10.1038/s41598-018-27521-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Oh, K. S., Soto, D. E., Smith, D. C., Montie, J. E., Lee, C. T., and Sandler, H. M. (2009). Combined-Modality Therapy with Gemcitabine and Radiation Therapy as a Bladder Preservation Strategy: Long-Term Results of a Phase I Trial. Int. J. Radiat. Oncology*Biology*Physics 74 (2), 511–517. doi:10.1016/j.ijrobp.2008.08.021

CrossRef Full Text | Google Scholar

Pulido, D., Sharma, U., Vadon-Le Goff, S., Hussain, S.-A., Cordes, S., Mariano, N., et al. (2018). Structural Basis for the Acceleration of Procollagen Processing by Procollagen C-Proteinase Enhancer-1. Structure 26 (10), 1384–1392. doi:10.1016/j.str.2018.06.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies. Nucleic Acids Res. 43 (7), e47. doi:10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

Rolih, V., Barutello, G., Iussich, S., De Maria, R., Quaglino, E., Buracco, P., et al. (2017). CSPG4: a Prototype Oncoantigen for Translational Immunotherapy Studies. J. Transl Med. 15 (1), 151. doi:10.1186/s12967-017-1250-4

CrossRef Full Text | Google Scholar

Shi, Y.-B., Li, J., Lai, X.-N., Jiang, R., Zhao, R.-C., and Xiong, L.-X. (2020). Multifaceted Roles of Caveolin-1 in Lung Cancer: A New Investigation Focused on Tumor Occurrence, Development and Therapy. Cancers 12 (2), 291. doi:10.3390/cancers12020291

PubMed Abstract | CrossRef Full Text | Google Scholar

Shimizu, T., Suzuki, H., Nojima, M., Kitamura, H., Yamamoto, E., Maruyama, R., et al. (2013). Methylation of a Panel of MicroRNA Genes Is a Novel Biomarker for Detection of Bladder Cancer. Eur. Urol. 63 (6), 1091–1100. doi:10.1016/j.eururo.2012.11.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, Y., Jin, D., Chen, J., Liang, W., and Liu, X. (2020). Effects of Arsenic (+3 Oxidation State) Methyltransferase Gene Polymorphisms and Expression on Bladder Cancer: Evidence from a Systematic Review, Meta-Analysis and TCGA Dataset. Toxicol. Sci. 177 (1), 27–40. doi:10.1093/toxsci/kfaa087

PubMed Abstract | CrossRef Full Text | Google Scholar

Sternberg, C. N., Bellmunt, J., Sonpavde, G., Siefker-Radtke, A. O., Stadler, W. M., Bajorin, D. F., et al. (2013). ICUD-EAU International Consultation on Bladder Cancer 2012: Chemotherapy for Urothelial Carcinoma-Neoadjuvant and Adjuvant Settings. Eur. Urol. 63 (1), 58–66. doi:10.1016/j.eururo.2012.08.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Subramanian, A., Kuehn, H., Gould, J., Tamayo, P., and Mesirov, J. P. (2007). GSEA-P: a Desktop Application for Gene Set Enrichment Analysis. Bioinformatics 23 (23), 3251–3253. doi:10.1093/bioinformatics/btm369

PubMed Abstract | CrossRef Full Text | Google Scholar

Szklarczyk, D., Morris, J. H., Cook, H., Kuhn, M., Wyder, S., Simonovic, M., et al. (2017). The STRING Database in 2017: Quality-Controlled Protein-Protein Association Networks, Made Broadly Accessible. Nucleic Acids Res. 45 (D1), D362–D368. doi:10.1093/nar/gkw937

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, Z., Li, C., Zhang, K., Yang, M., and Hu, X. (2017). GE-mini: a mobile APP for Large-Scale Gene Expression Visualization. Bioinformatics w775, btw775. doi:10.1093/bioinformatics/btw775

PubMed Abstract | CrossRef Full Text | Google Scholar

Tooker, P., Yen, W.-C., Ng, S.-C., Negro-Vilar, A., and Hermann, T. W. (2007). Bexarotene (LGD1069, Targretin), a Selective Retinoid X Receptor Agonist, Prevents and Reverses Gemcitabine Resistance in NSCLC Cells by Modulating Gene Amplification. Cancer Res. 67 (9), 4425–4433. doi:10.1158/0008-5472.can-06-4495

PubMed Abstract | CrossRef Full Text | Google Scholar

Udhaya Kumar, S., Thirumal Kumar, D., Siva, R., George Priya Doss, C., Younes, S., Younes, N., et al. (2020). Dysregulation of Signaling Pathways Due to Differentially Expressed Genes from the B-Cell Transcriptomes of Systemic Lupus Erythematosus Patients - A Bioinformatics Approach. Front. Bioeng. Biotechnol. 8, 276. doi:10.3389/fbioe.2020.00276

PubMed Abstract | CrossRef Full Text | Google Scholar

Uhlen, M., Oksvold, P., Fagerberg, L., Lundberg, E., Jonasson, K., Forsberg, M., et al. (2010). Towards a Knowledge-Based Human Protein Atlas. Nat. Biotechnol. 28 (12), 1248–1250. doi:10.1038/nbt1210-1248

PubMed Abstract | CrossRef Full Text | Google Scholar

Vena, F., Bayle, S., Nieto, A., Quereda, V., Aceti, M., Frydman, S. M., et al. (2020). Targeting Casein Kinase 1 Delta Sensitizes Pancreatic and Bladder Cancer Cells to Gemcitabine Treatment by Upregulating Deoxycytidine Kinase. Mol. Cancer Ther. 19 (8), 1623–1635. doi:10.1158/1535-7163.MCT-19-0997

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, D., and Lippard, S. J. (2005). Cellular Processing of Platinum Anticancer Drugs. Nat. Rev. Drug Discov. 4 (4), 307–320. doi:10.1038/nrd1691

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Svendsen, A., Kmiecik, J., Immervoll, H., Skaftnesmo, K. O., Planagumà, J., et al. (2011). Targeting the NG2/CSPG4 Proteoglycan Retards Tumour Growth and Angiogenesis in Preclinical Models of GBM and Melanoma. PLoS ONE 6 (7), e23062. doi:10.1371/journal.pone.0023062

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, S., Zhong, L., Li, Y., Xiao, D., Zhang, R., Liao, D., et al. (2019). Up-regulation of PCOLCE by TWIST1 Promotes Metastasis in Osteosarcoma. Theranostics 9 (15), 4342–4353. doi:10.7150/thno.34090

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X., Wu, Q., Xu, B., Wang, P., Fan, W., Cai, Y., et al. (2015). miR-124 Exerts Tumor Suppressive Functions on the Cell Proliferation, Motility and Angiogenesis of Bladder Cancer by fine-tuning UHRF1. Febs J. 282 (22), 4376–4388. doi:10.1111/febs.13502

PubMed Abstract | CrossRef Full Text | Google Scholar

Welsh, J. L., Wagner, B. A., van’t Erve, T. J., Zehr, P. S., Berg, D. J., Halfdanarson, T. R., et al. (2013). Pharmacological Ascorbate with Gemcitabine for the Control of Metastatic and Node-Positive Pancreatic Cancer (PACMAN): Results from a Phase I Clinical Trial. Cancer Chemother. Pharmacol. 71 (3), 765–775. doi:10.1007/s00280-013-2070-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, K., Mu, X. Y., Jiang, J. T., Tan, M. Y., Wang, R. J., Zhou, W. J., et al. (2018). miRNA-26a-5p and miR-26b-5p I-nhibit the P-roliferation of B-ladder C-ancer C-ells by R-egulating PDCD10. Oncol. Rep. 40 (6), 3523–3532. doi:10.3892/or.2018.6734

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiang, Z., Li, J., Song, S., Wang, J., Cai, W., Hu, W., et al. (2019). A Positive Feedback between Ido1 Metabolite and COL12A1 via MAPK Pathway to Promote Gastric Cancer Metastasis. J. Exp. Clin. Cancer Res. 38 (1), 314. doi:10.1186/s13046-019-1318-5

CrossRef Full Text | Google Scholar

Xie, F., Zhao, N., Zhang, H., and Xie, D. (2020). Circular RNA CircHIPK3 Promotes Gemcitabine Sensitivity in Bladder Cancer. J. Cancer 11 (7), 1907–1912. doi:10.7150/jca.39722

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Y., Cheng, Z., Jia, X., Shi, N., Xia, Z., Zhang, W., et al. (2019). Mortality Trends of Bladder Cancer in China from 1991 to 2015: an Age-Period-Cohort Analysis. Cmar 11, 3043–3051. doi:10.2147/CMAR.S189220

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, G., Wang, L.-G., Han, Y., and He, Q.-Y. (2012). clusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters. OMICS: A J. Integr. Biol. 16 (5), 284–287. doi:10.1089/omi.2011.0118

CrossRef Full Text | Google Scholar

Yu, W. D., Wang, H., He, Q. F., Xu, Y., and Wang, X. C. (2018). Long Noncoding RNAs in Cancer‐immunity Cycle. J. Cell Physiol 233 (9), 6518–6523. doi:10.1002/jcp.26568

CrossRef Full Text | Google Scholar

Zhou, W., He, L., Dai, Y., Zhang, Y., Wang, J., and Liu, B. (2018). MicroRNA-124 I-nhibits C-ell P-roliferation, I-nvasion and M-igration by T-argeting CAV1 in B-ladder C-ancer. Exp. Ther. Med. 16 (4), 2811–2820. doi:10.3892/etm.2018.6537

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, Y., Zhou, B., Pache, L., Chang, M., Khodabakhshi, A. H., Tanaseichuk, O., et al. (2019). Metascape Provides a Biologist-Oriented Resource for the Analysis of Systems-Level Datasets. Nat. Commun. 10 (1), 1523. doi:10.1038/s41467-019-09234-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, H., Chen, H., Wang, J., Zhou, L., and Liu, S. (2019). Collagen Stiffness Promoted Non-muscle-invasive Bladder Cancer Progression to Muscle-Invasive Bladder Cancer. Ott 12, 3441–3457. doi:10.2147/OTT.S194568

CrossRef Full Text | Google Scholar

Zo, R. B., and Long, Z. (2019). MiR‐124‐3p Suppresses Bladder Cancer by Targeting DNA Methyltransferase 3B. J. Cell Physiol 234 (1), 464–474. doi:10.1002/jcp.26591

PubMed Abstract | CrossRef Full Text | Google Scholar