CXCL12, A Potential Biomarker and A Vital Modulator of Tumor Immune Microenvironment (TIME) of Bladder Cancer: From A Comprehensive Analysis of TCGA Database

Background: Tumor immune microenvironment (TIME) played a signicant role in the initiation and progression of bladder cancer (BC). However, there are few researches regarding the association between immune-related genes (IRGs) and tumor-inltrating immune cells (TICs) in TME of BC. Methods: We calculated the proportion of immune/stromal component and TICs in TME of 414 BC samples and 19 normal samples downloaded from The Cancer Genome Atlas (TCGA) database with the help of ESTIMATE and CIBERSORT algorithms. Differentially expressed genes (DEGs) were obtained from the comparison between Stromal and Immune Score and further analyzed by GO and KEGG enrichment analysis, as well as protein–protein interaction (PPI) network and COX regression analysis. CXC chemokine ligand-12 (CXCL12) was overlapping from the above analysis. Afterwards, single gene analysis of CXCL12 was carried out through difference analysis, paired analysis and Gene Set Enrichment Analysis (GSEA). The association between the expression of CXCL12 and TICs was assessed by difference analysis and correlation analysis. Results: The results indicated that immune and stromal component in TME of BC were associated with patients’ clinic-pathological characteristics. We further conrmed that 284 DEGs were primarily enriched in immune-associated activities and CXCL12 was the most signicant gene, which shared the leading nodes in PPI network and closely related with BC patients’ survival. Single gene analysis revealed that CXCL12 was down-regulated in BC samples and signicantly related with the clinic-pathological characteristics of patients. Further analysis indicated that CXCL12 greatly participated in immune-associated activities through closely communicating with TICs in TIME of BC. Conclusions: CXCL12 might be a potential biomarker and a vital modulator of TIME through communicating with multiple TICs, which might provide an extra insight for the immunotherapy of BC.


Introduction
Over the last 30 years, few advances have been made in the treatment of bladder cancer (BC) since the application of cisplatin-based chemotherapy in the mid-1980s [1] . Fortunately, immunotherapy has emerged as a novel potential therapy recently. Several clinical trials, such as IMvigor 210 study [2] and CheckMate 275 study [3] , have revealed that BC patients greatly bene ted from the treatment of immunecheckpoint blockade (ICB). Although multiple biomarkers have been associated with the prediction of immunotherapy effect, such as the expression of PD-L1 and TMB, it is still less than satisfactory to select BC patients who are likely to bene t most from immunotherapy [4] .
In recent years, the tumor microenvironment (TME) has caught our attention as its role in modulating the initiation and progression of cancers, including BC [5][6][7][8][9][10] . TME is composed of nonmalignant cells, vessels, lymphoid organs or lymph nodes, nerves, intercellular components and metabolites [11] . In brief, stromal component and immune component constitute the TME [12] . Furthermore, accumulating researches indicated that tumor immune microenvironment (TIME) had great potential to in uence tumor initiation, predict immunotherapeutic responsiveness and reveal new therapeutic targets [13] . Up to date, numerous studies revealed that immune-related genes (IRGs), obtained from TME, functioned to predict the survival of cancer patients, including breast cancer [14] , endometrial cancer [15] , liver cancer [16] , gastric cancer [5] , bladder cancer [17] and so on. For example, Q. Ding, et al. [18] found that a nine-gene signature was closely related with immune in ltration in TME and the survival of patients with ovarian cancer. Besides, tumor-in ltrating immune cells (TICs) in TIME, such as tumor-in ltrating lymphocytes (TILs), tumor-associated macrophages (TAMs) and so on, have the potential to be biomarkers and predict the prognosis of multiple cancers [19][20][21] . For instance, elevated level of TILs were associated with better overall survival (OS), higher pathological complete response (pCR) rate and fewer recurrence, as well as greatly bene ted from trastuzumab treatment in breast cancer [22][23][24][25] . De cient CD4 + T cell help suppressed the response of cytotoxic T lymphocytes (CTLs), which were able to establish e cient and durable antitumor activity [26] . TAMs acted to inhibit T cell recruitment and function as well as modulate the immunity of various tumors, thus affecting the response of immunotherapy [27] . Galectin-9 + TAMs predicted prognosis and response to adjuvant chemotherapy in BC patients [28] . However, there are few researches regarding the association between IRGs and TICs in the TME of BC. The exploration of the relationship between IRGs and TICs in TIME might provide us a new sight into the progression and immunotherapy of BC.
Fortunately, with the rapid development of transcriptome pro ling based on functional genomics analysis, comprehensive analysis of IRGs and TICs in the TIME of BC has become possible. In our study, we applied ESTIMATE and CIBERSORT algorithms to calculate the ratio of immune and stromal component, as well as TICs proportion of BC samples from The Cancer Genome Atlas (TCGA) database.
Next, we started with differentially expressed genes (DEGs) acquired by the comparison of immune and stromal component in BC samples, and found out that CXC chemokine ligand-12 (CXCL12) acted to be a potential biomarker of BC and a vital regulator of immune-associated activities in TIME of BC via communicating with multiple TICs.

Data acquisition
The transcriptome pro ling and clinical data of 414 BC samples and 19 normal samples were downloaded from TCGA database (https://portal.gdc.cancer.gov/). Difference analysis between Scores and the clinic-pathological characteristics Difference analysis was conducted to learn the correlation between Immune/Stromal/ESTIMATE Score and clinic-pathological characteristics, such as age, gender, grade, stage and TNM classi cation. It was also carried out to nd the association between the expression of CXCL12 and clinic-pathological characteristics as well as the association between the expression of CXCL12 and TICs. Wilcoxon rank sum and Kruskal-Wallis rank sum test were applied for the comparison. P < 0.05 was considered to be statistically signi cant.
DEGs acquisition 414 BC samples were grouped into to two subgroups, including high Immune/Stromal Score group and low Immune/Stromal Score group based on the comparison with the median. Package limma in R was applied for the analysis. A fold change (FC, log2 (high-score /low-score) ) > 2 and false discovery rate (FDR) <0.05 were used to search the DEGs. Pheatmap package in R was used to plot the heatmaps of DEGs.

Intersection analysis
VennDiagram package in R was used to plot the venn diagram of DEGs.
Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis GO and KEGG enrichment analysis of the above DEGs were further carried out by using the ClusterPro ler, enrichplot, and ggplot2 packages in R. P< 0.05 was considered to be statistically signi cant.

Protein-Protein Interaction (PPI) Network and Gene Set Enrichment Analysis (GSEA)
The preliminary PPI network of 284 DEGs was acquired from the Search Tool for Retrieval of Interacting Genes/Proteins (STRING) database (version 11.0) and further reconstructed in Cytoscape (3.6.1). The con dence of interactive relationship of the nodes was >0.95. GSEA (4.1.0) based on the different gene sets, was applied to learn the speci c functional pro le of CXCL12. P< 0.05 was considered to be statistically signi cant.

COX Regression Analysis
Package survival in R was applied to conduct univariate COX regression.

Difference analysis and paired analysis of CXCL12
Package Beeswarm in R was used to assess the expression of CXCL12 in bladder tumor samples and normal samples. Packages Ggpubr and BiocManager were applied to learn the expression of CXCL12 in bladder tumor samples and the paired normal samples. Wilcoxon test was used for the comparison. P < 0.05 was considered to be statistically signi cant.

Survival analysis
Packages survival and survminer in R were used to carry out the survival analysis. Kaplan-Meier plot and log-rank tests were conducted to learn the associations between the expression of CXCL12 and the survival of BC patients. P < 0.05 was considered to be statistically signi cant.

TICs Pro le
TICs pro le in BC samples was evaluated by CIBERSORT algorithm.

Analysis process of the study
We rst downloaded the transcriptome pro ling and clinical data of 414 BC samples and 19 normal samples from TCGA database and further analyzed the ratio of immune and stromal component in TME of each tumor sample through ESTIMATE algorithm. Difference analysis was conducted to nd out the association between immune/stromal component in TME and the clinic-pathological characteristics of BC patients. A total of 284 DEGs were further acquired based on the immune and stromal component in TME of BC. GO and KEGG enrichment analysis of these 284 DEGs were performed to learn their biological functions. Finally, CXCL12 was found out through the intersection analysis of PPI network and univariate COX regression. Then, we focused on the expression pro le of CXCL12 in BC samples and normal samples, and the association between the expression of CXCL12 and survival and clinic-pathological characteristics of BC patients. GSEA of different gene collections was also carried out to learn the functional pro les of CXCL12. TICs pro le in TIME of BC samples was calculated by CIBERSORT algorithm. Difference analysis and correlation analysis were applied to nd out the correlation between the expression of CXCL12 and TICs.
2. Immune and stromal component in TME were associated with the clinic-pathological characteristics of BC To learn the association between the proportion of immune/stromal component in TME and the clinicpathological characteristics of BC patients, we rst set up the immune/stromal component evaluating system of TME in BC samples. Immune, Stromal and ESTIMATE Score represented the proportion of immune, stromal and the total component in TME of each tumor sample, respectively. The higher Score suggested the larger proportion of immune or stromal component in TME. The clinic-pathological characteristics of BC patients, including age, gender, grade, stage and TNM classi cation, were concluded in Supplement Table 1. Difference analysis revealed that Stromal Score was associated with patients' age and gender, especially the grade (Figure 2a-c, p=0.018; 0.032; <0.001). The higher Stromal Score indicated the higher grade of BC (Figure 2c, p<0.001). It was also found that Stromal Score was related with T, N classi cation and clinical stage of BC (Figure 2d-g). Regarding Immune Score, it was signi cantly up-regulated in female patients (Figure 2i, p=0.037) and the higher Immune Score also suggested the higher grade of BC (Figure 2j, p<0.001). However, neither TNM classi cation nor clinical stage of BC was associated with Immune Score (Figure 2k-n, p>0.05). As for ESTIMATE Score, the results similarly showed that it was up-regulated in female patients (Figure 2p, p=0.028) and indicated the higher grade of BC (Figure 2q, p<0.001). In addition, ESTIMATE Score was connected with T classi cation and clinical stage of BC (Figure 2r-u). From above, we might conclude that the stromal and immune component in TME were signi cantly associated with the grade of BC, and partly related with TNM classi cation and clinical stage. This provided us a new sight to explore the underlying mechanisms of the development and progression of BC.

DEGs were primarily enriched in immune-associated activities
In order to further understand the underlying mechanisms of BC TME, we conducted difference analysis to acquire the DEGs pro le in TME. BC samples were rst divided into two groups, including high Stromal Score group and low Stromal Score group (or high Immune Score group and low Immune Score group). Heat map displayed the DEGs pro le in TME of BC (Figure 3a-b). To be speci c, there were 537 DEGs obtained from the comparison between high Stromal Score and low Stromal Score, among which 496 DEGs were up-regulated and 41 DEGs were down-regulated (Figure 3c, Supplement Table 2). Besides, there were 531 DEGs obtained from the comparison between high Immune Score and low Immune Score, among which 466 DEGs were up-regulated and 65 DEGs were down-regulated (Figure 3c, Supplement Table 3). In conclusion, a total of 284 DEGs were commonly up-regulated or down-regulated in Stromal or Immune Score group (Figure 3c, Supplement Table 4). Next, we carried out GO and KEGG enrichment analysis to assess the biological functions of these 284 DEGs. Go enrichment analysis revealed that 284 DEGs were primarily enriched in immune-associated activities, such as leukocyte proliferation, migration and chemotaxis, lymphocyte and mononuclear cell proliferation and so on (Figure 3d-f). Similarly, KEGG enrichment analysis suggested that these DEGs were mainly enriched in immune-associated activities, including complement and coagulation cascades, cytokine-cytokine receptor interaction, B cell receptor signaling pathway, chemokine signaling pathway and so on. Therefore, we considered that DEGs acquired from Stromal and Immune Score group of BC were signi cantly associated with immuneassociated activities in TME.

PPI network and univariate COX regression analysis of 284 DEGs
To learn the detailed reciprocity among these 284 DEGs, we further ploted PPI network through STRING database and cytoscape software. The potential interactions among DEGs were displayed in Figure 4a.
The top ten DEGs, which shared the leading nodes, were ranked in Figure 4b, including ITGAM, CXCL12, CXCL13 and so on (Supplement Table 5). Univariate COX regression of DEGs was conducted at the same time to nd out the speci c DEGs, which were signi cantly associated with BC patients' survival ( Figure  4c). MMP9, COMP, F13A1 and CXCL12 came to our eyes (Figure 4c, p<0.05; Supplement Table 6). Intersection analysis was nally carried out to search the DEGs who shared the leading nodes in PPI network and were signi cantly related with BC patients' survival ( Figure 4d). Fortunately, CXCL12 emerged.

CXCL12 was down-regulated in BC tissues and associated with the clinic-pathological characteristics
The CXC chemokine CXCL12 greatly participated in multiple physiological and pathological processes through interacting with its receptors CXC chemokine receptor 4 (CXCR4) and atypical chemokine receptor 3 (ACKR3) [29] . In our study, CXCL12 was signi cantly down-regulated in BC tissues compared with normal tissues (Figure 5a, p<0.001). Besides, the paired analysis also indicated that CXCL12 was obviously down-regulated in BC tissues compared with the paired normal tissues (Figure 5b, p<0.001). However, survival analysis suggested that the expression of CXCL12 was not signi cantly associated with the survival of BC patients (Figure 5c, p=0.050). As for clinic-pathological characteristics of BC patients, the results showed that the expression of CXCL12 changed with ages and higher expression of CXCL12 indicated the higher grade of BC (Figure 5d, f; p<0.001, p<0.001). In addition, the expression of CXCL12 was partly correlated with T, N classi cation and clinical stages of BC (Figure 5h-k). Hence, we might conclude that CXCL12 had the potential to participate in the progression of BC.

CXCL12 greatly participated in immune-associated activities
On account of the potential role of CXCL12 in the progression of BC, we further explored the underlying mechanisms of CXCL12. BC cases were rst divided into two groups, including high CXCL12 expression group and low CXCL12 expression group. GSEA suggested that for C5 collection, the gene ontology sets, genes in high CXCL12 expression group were enriched in cytokine binding, leukocyte migration and negative regulation of interleukin 6 production (Figure 6a, p<0.05). For C2 collection, the KEGG gene sets database, genes in high CXCL12 expression group were enriched in allograft rejection, antigen processing and presentation, B cell receptor signaling pathway, complement and coagulation cascades, and so on (Figure 6b, p<0.05). For hallmark gene sets, genes in high CXCL12 expression group were primarily enriched in allograft rejection, complement, IL2 stats signaling, IL6 jak stat3 signaling and so on. Finally, for C7 collection de ned by MSigDB, the immunologic gene sets, genes were found to be enriched in multiple immune-associated activities, which were related with CD4 T cell, CD8 T cell, B cell and so on (Figure 6d, p<0.05). As a result, according to GSEA, CXCL12 was signi cantly related with immuneassociated activities and might play an important role in regulating TIME of BC.

CXCL12 communicated with TICs in TIME
Since GSEA suggested that CXCL12 was signi cantly correlated with immune-associated activities, we speculated that there might be underlying connections between CXCL12 and TIME of BC. In consequence, we calculated the proportion of TICs in BC samples through CIBERSORT algorithm. The speci c ratio of 22 TICs in tumor samples was shown in Figure 7a. And the association among these TICs was exhibited in Figure 7b. Difference analysis was carried out to learn the association between the expression of CXCL12 and speci c TICs. The results showed that 10 kinds of TICs, such as Macrophages M2, B cells naïve, T cells gamma delta, T cells CD4 naïve, T cells follicular helper and so on, were signi cantly associated with the expression of CXCL12 (Figure 7c, Supplement Table 7). The correlation analysis suggested that the expression of CXCL12 were positively or negatively related with 8 kinds of TICs, among which B cells naïve, macrophages M2 and mast cells resting were positively associated, and dendritic cells activated, dendritic cells resting, NK cells resting, T cells CD4 naïve and T cells follicular helper were negatively associated (Figure 7e-l, p<0.05; Supplement Table 7). Finally, intersection analysis between difference analysis and correlation analysis suggested that 6 kinds of TICs were signi cantly related with the expression of CXCL12 (Figure 7d). In conclusion, CXCL12 was obviously related with various TICs in TIME of BC and might greatly participated in regulating the TIME of BC via closely communicating with multiple TICs.

Discussion
In this study, we focused on the associations between IRGs and TICs in the TME of BC. Firstly, we downloaded the transcriptome pro ling and clinical data of 414 BC samples and 19 normal samples from TCGA database, and calculated the proportion of immune and stromal component in TME of each BC samples with the help of ESTIMATE algorithm. Secondly, difference analysis found out that Immune/Stromal Score was associated with the clinic-pathological characteristics of BC. Thirdly, DEGs were obtained through the comparison between high Immune Score and low Immune Score (or high Stromal Score and low Stromal Score), and GO and KEGG enrichment analysis suggested that these DEGs were mainly enriched in immune-related activities. Fourthly, PPI network and univariate COX regression analysis indicated that CXCL12 shared the leading nodes in PPI network and was signi cantly related with BC patients' survival. Fifthly, signal gene analysis was conducted and found out that CXCL12 was down-regulated in BC samples compared with normal sample and signi cantly related with the clinic-pathological characteristics of BC. GSEA revealed that CXCL12 was signi cantly associated with immune-associated activities and might play an important role in regulating TIME of BC. Finally, CIBERSORT algorithm was applied to calculate the proportion of TICs in BC samples, and difference analysis as well as correlation analysis suggested that CXCL12 was obviously connected with multiple TICs in TME of BC and might modulate the TIME of BC via closely communicating with TICs.
CXCL12 was traditionally classi ed as a homeostatic CXC chemokine and took a great part in modulating kinds of physiological and pathological processes via binding to its receptors CXCR4 and ACKR3 [29,30] . CXCL12/CXCR4 axis had been proved to be associated with the progression and therapy of cancers. For example, CXCL12/CXCR4 axis advanced pancreatic cancer development, invasion, and metastasis through complex crosstalk with other pathways, and was correlated with the poor prognosis of patients [31] . Treatment targeting the CXCL12/CXCR4 pathway increased the e cacy of radiotherapy of locally advanced cervical cancer [32] . The role of CXCL12 in tumor development mainly depended on the speci c microenvironment of tumors [33] . In addition, it was found that the CXCL12-CXCR4/CXCR7 axis had a great in uence in gastrointestinal malignancies through immune resistance [34] . Further researches suggested that the mechanisms of immunotherapy resistance might be associated with the CXCL12/CXCR4 axis [35] . Up to date, only few documents explored the association between CXCL12 and the tumorigenesis and progression of BC [36] . Researches indicated that CXCL12 was down-regulated in BC tissues compared with normal bladder mucosal tissues, and positively associated with the differentiation degree and invasive depth of BC tissues [37] . Additionally, CXCL12/CXCR4 advanced invasion of BC cells through activating Stat3 [38] .
Furthermore, increasing studies have revealed that tumor-in ltrating immune cells played an important role in the progression and treatment of BC. For example, intratumoral TIGIT + CD8 + T-cell abundance functioned as a potential prognostic factor for patients' survival and a predictive biomarker for adjuvant chemotherapeutic effect [39] . CD19 + tumor-in ltrating B-cells activated CD4 + tumor-in ltrating T-cells in the TMB of BC and acted as an independent prognostic factor for post-surgery survival and adjuvant chemotherapy bene ts of BC patients [40] . Besides, DC-SIGN + TAM in ltration was signi cantly associated with a tumor-promoting TIME and might functioned as a prognostic indicator and therapeutic target in the immunotherapy of BC [41] . However, there are no researches focusing on the associations between CXCL12 and TICs in TIME of BC.
In our study, we con rmed that CXCL12 was signi cantly down-regulated in BC tissues and associated with patients' clinic-pathological characteristics. In addition, functional analysis revealed that CXCL12 greatly participated in immune-associated activities and might regulate TIME of BC through communicating with multiple TICs, such as macrophages M2, B cells naïve, T cells follicular helper, mast cells resting, dendritic cells resting and T cells CD4 naïve. Further researches should be carried out to clarify the accuracy of the above combined analysis, and focused on the underlying mechanisms of the communication between CXCL12 and TICs in TIME of BC.

Declarations
Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.

Availability of data and material
The datasets generated and analyzed in this research are available in TCGA database (https://portal.gdc.cancer.gov)

Competing interests
The authors declare that they have no competing interests.