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

To identify key biomarkers in gemcitabine (GEM)-resistant bladder cancer (BCa) and investigate their associations with tumor-infiltrating immune cells in a tumor immune microenvironment, we performed the present study on the basis of large-scale sequencing data. Expression profiles from the Gene Expression Omnibus GSE77883 dataset and The Cancer Genome Atlas BLCA dataset were analyzed. Both BCa development and GEM-resistance were identified to be immune-related through evaluating tumor-infiltrating immune cells. Eighty-two DEGs were obtained to be related to GEM-resistance. Functional enrichment analysis demonstrated they were related to regulation of immune cells proliferation. Protein–protein interaction network selected six key genes (CAV1, COL6A2, FABP4, FBLN1, PCOLCE, and CSPG4). Immunohistochemistry confirmed the down-regulation of the six key genes in BCa. Survival analyses revealed the six key genes were significantly associated with BCa overall survival. Correlation analyses revealed the six key genes had high infiltration of most immune cells. Gene set enrichment analysis further detected the key genes might regulate GEM-resistance through immune response and drug metabolism of cytochrome P450. Next, microRNA-gene regulatory network identified three key microRNAs (hsa-miR-124-3p, hsa-miR-26b-5p, and hsa-miR-192-5p) involved in GEM-resistant BCa. Connectivity Map analysis identified histone deacetylase inhibitors might circumvent GEM-resistance. In conclusion, CAV1, COL6A2, FABP4, FBLN1, PCOLCE, and CSPG4 were identified to be critical biomarkers through regulating the immune cell infiltration in an immune microenvironment of GEM-resistance and could act as promising treatment targets for GEM-resistant muscle-invasive BCa.


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 GEMresistance. 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 . In addition, bioinformatics analyses constructs a microRNA (miRNA)-gene regulatory network associated with alteration of memory CD4 + T cells in GEMresistant pancreatic cancer cells, which suggests the immune system is implicated in the microenvironment of GEMresistance (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.

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 Populationscounter) 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 GEMresistant 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 .

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.
To confirm the differential expression of the six prognosisrelated 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).
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 immunotherapytreated 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.
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. Figure 1 showed the workflow.

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 GEMresistance 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 GEMresistance 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).

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 downregulated 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 upregulated and 1,289 down-regulated genes were obtained in GEM-resistant BCa cells ( Figure 3A and Figure 3C).
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). Frontiers in Cell and Developmental Biology | www.frontiersin.org January 2022 | Volume 9 | Article 809620 6 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).
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.

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.

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.
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 downregulated 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   Figure S4). IMvigor210 cohort indicated that COL6A2 (HR > 1), FABP4 (HR > 1), and FBLN1 (HR < 1) could predict the  Frontiers in Cell and Developmental Biology | www.frontiersin.org January 2022 | Volume 9 | Article 809620 8 OS after immunotherapy with atezolizumab (p < 0.05) (Figure 7).

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 cancerassociated 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).

Genetic mutations of hub genes in TCGA BLCA dataset
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.
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.
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 GEMresistance (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 chemotherapyrelated 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).

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. Frontiers in Cell and Developmental Biology | www.frontiersin.org January 2022 | Volume 9 | Article 809620 The top 10 antagonistic drugs for GEM-resistant BCa were lisinopril, rifabutin, clonidine, prasterone, vorinostat, prednisone, nifenazone, alvespimycin, trichostatin A, and tanespimycin.

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 . 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 immunerelated 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. Downregulation 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 GEMresistant 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.  suppressed the proliferation and invasion of tumor cells by targeting CAV1 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 The bold values indicated statistical significance.
Frontiers in Cell and Developmental Biology | www.frontiersin.org January 2022 | Volume 9 | Article 809620 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 . 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  -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).
Frontiers in Cell and Developmental Biology | www.frontiersin.org January 2022 | Volume 9 | Article 809620 15 PDCD10 . 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 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.
Frontiers in Cell and Developmental Biology | www.frontiersin.org January 2022 | Volume 9 | Article 809620 16 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. Frontiers in Cell and Developmental Biology | www.frontiersin.org January 2022 | Volume 9 | Article 809620

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.