The CAFs Related Gene CALD1 Is a Prognostic Biomarker and Correlated With Immune In ltration in Bladder Cancer

YiHeng Du Suzhou Kowloon Hospital Shanghai Jiao Tong University School of Medicine https://orcid.org/00000002-5092-6448 Xiang Jiang Suzhou Kowloon Hospital Shanghai Jiao Tong University School of Medicine Bo Wang Suzhou Kowloon Hospital Shanghai Jiao Tong University School of Medicine Jin Cao Suzhou Kowloon Hospital Shanghai Jiao Tong University School of Medicine Yi Wang Suzhou Kowloon Hospital Shanghai Jiao Tong University School of Medicine Jiang Yu Suzhou Kowloon Hospital Shanghai Jiao Tong University School of Medicine XiZhi Wang Suzhou Kowloon Hospital Shanghai Jiao Tong University School of Medicine HaiTao Liu (  haitao.liu@shgh.cn ) Shanghai First People's Hospital, Shanghai Jiaotong University School of Medicine, Shanghai 200080 https://orcid.org/0000-0002-3674-0081


Introduction
Bladder Cancer (BLCA) is one of the most common malignant tumors of the genitourinary system (1).
Based on the invasion depth, BLCA can be divided into Non-muscle Invasive (NMIBC) and Muscle Invasive Bladder Cancer (MIBC). Notably, NMIBC has a better prognosis and more treatment options compared to MIBC (2). Besides, intravesical instillation of BCG was rst reported to be effective in treating NMIBC by Morales et al. in 1976, with complete remission rates of up to 70% − 80% (3). This was the earliest application of immunotherapy for the management of BLCA and has been used up to date. Moreover, BLCA represents a growing number of cancers characterized by the in ltration of a signi cant number of immune cells in the TME (4,5), making it suitable for immunotherapy.
Tumorigenesis and tumor development are complex processes affected by many factors. In addition to genetic changes and epigenetic defects, the tumor microenvironment (TME), which consists of cellular and non-cellular components, also plays an essential regulatory role (6). Tumor stromal cells are the main cellular components of the TME. Furthermore, stromal cells in the TME, especially Cancer-Associated Fibroblasts (CAFs) (7), have multiple effects on cancer growth and maintenance.
CAFs interact with immune components by secreting various factors such as collagen, Matrix Metalloproteinases (MMPs), and chemokines (8). Additionally, Tumor-In ltrating Immune Cells (TIICs) in the TME were critically associated with cancer outcomes. Out of all the TIICs, macrophages are always among the most abundant in the TME, including in BLCA (9). Moreover, existing evidence indicates that CAFs interact with the M2 macrophages, referred to as Tumor-Associated Macrophages (TAMs), promoting immunosuppression and inducing the occurrence and progression of cancers (10).
Furthermore, bioinformatics methods are now widely used in cancer research. Through high-throughput sequencing, mechanisms underlying pathological processes, including cancers, can be revealed by comparing different genes' expression networks. Currently, the CIBERSORT (11) and MCP-COUNTER (12) algorithms for calculating the abundance of immune cells are widely used in studies of the TME. CAFs scores and relative levels of 22 TIICs can be computed using the MCP-COUNTER and CIBERSORT algorithms, which provide signi cant help in studying the relationship between CAFs and TIICs. Moreover, genes highly co-expressed in cancer can be identi ed through the Weighted Gene Co-expression Network Analysis (WGCNA) (13). Therefore, through these bioinformatics means, the present study uncovered and validated that Caldesmon 1 (CALD1), a key gene associated with CAFs, was crucial in regulating both the stromal and immune microenvironment of BLCA. Consequently, it may become a promising biomarker of BLCA progression.

Methods And Materials
The source of data

Calculation of CAFs scores
The MCP-counter algorithm (12,14), provided by TIMER 2.0 (http://timer.cistrome.org/), was used to calculate the CAFs score of patients recruited from both the TCGA and GEO databases.

Estimation of the tumor microenvironment
The ESTIMATE algorithm in the estimate package of the R software (version 4.0.2) was used. Three scoring forms, namely immune score, stromal scare and ESTIMATE score, were positively correlated with the proportion of immune and stromal components as well as the sum of both. Therefore, it is indicated that the higher the score, the more signi cant the ratio of the corresponding parts in the TME.

Screening for differentially expressed genes ( DEGs )
To screen for DEGs in both the TCGA and GEO cohorts, we divided the patients into two groups (High and low CAFs) according to the medium level of the CAFs score. After that, analysis of differential gene expression between these two groups was performed using the Bioconductor limma package of the R software (version 4.0.2). Genes with p 0.01 and |logFC| 1 were considered to be signi cantly differentially expressed.

WCGNA for the determination of key genes
The Weighted Gene Correlation Network Analysis (WGCNA) R-package was used for co-expression network analysis. CAFs related genes were selected from the most signi cant modules related to high levels of CAFs. Key genes were obtained from the intersection of CAFs related genes and DEGs in both the TCGA and GEO cohorts.

GO and KEGG enrichment analysis
The "clusterpro ler" R package was used to perform GO and KEGG enrichment analyses for hub genes.
GO terms or KEGG pathways with corrected P-values 0.05 were considered to be signi cantly enriched in the hub genes.

Gene Set Enrichment Analysis (GSEA)
The Hallmark and C2 Kegg gene sets v7.2 were used for GESA, which was performed using the GSEA software (version 4.1.0) obtained from the Broad Institute. Gene sets with NOM p < 0.05 and False Discovery Rate (FDR) q < 0.05 were considered to be signi cant.

Survival Analysis
Survival analysis was conducted using the "survival" and "survminer" packages of the R software.
Additionally, the Kaplan-Meier method with the best cut-off value was used to draw the survival curves. P values from the log-rank test that was less than 0.05 were considered to be signi cant

TIICs Pro le
The R software's CIBERSORT algorithm was used to determine the pro le of TIICs (including 22 immune cells) in all tumor samples of the TCGA cohort. 235 patients with p < 0.05 were selected for subsequent analysis.

Immunohistochemistry analysis
Here, we used ACTA2, a marker of CAFs, and CD206, a marker of macrophage M2, to discuss the relation between CALD1, CAFs, and macrophage M2. The expression of CALD1 (1:250, Abcam, ab32330), ACTA2 (1:250, Abcam, ab7817), and CD206 (1:4000, Abcam,ab252921) in tumor tissues was detected using the BenchMark GX automatic immunohistochemical staining system (Roche, Switzerland). After depara nization, the tissue sections were incubated with primary antibody for 32 minutes. Biotinylated anti-IgG antibody and horseradish peroxidase were used to show positive expression areas. Hematoxylin was used for counterstaining and Bluing Reagent for post counterstaining.

Statistics analysis
Univariate Cox regression was performed using the "survival" package of the R-software. Genes with p values of less than 0.05 were considered to be related to survival. Moreover, survival-related hub genes were identi ed from the intersection of survival-related genes in the TCGA and GEO cohorts. Correlation analysis was conducted through the Pearson correlation. Gene expression was transformed by log2 (FPKM + 1). Differential analysis between the high and low CALD1 expression groups was then conducted through the Wilcoxon test. The Kruskal test estimated statistical signi cance for variables of more than two groups. P-value < 0.05 was considered signi cant.

Results
1. Abundance of CAFs is a poor prognostic factor associated with the progression of BLCA.
We rst calculated the abundance of CAFs in the TCGA cohort using the MCP-COUNTER algorithm. The results indicated that CAFs were more abundant than any other cell types in the tumor microenvironment ( Figure 1A) and had a signi cant correlation with the stromal score (R=0.73), immune score (R=0.37) as well as the ESTIMATE score (R=0.59) ( Figure 1B). Moreover, high levels of CAFs were signi cantly associated with low survival in BLCA patients, as shown in Figure 1C (p=0.003). Furthermore, we compared the effects of different cell types in the microenvironment on the clinical characteristics of BLCA. Notably, CAFs were found to have a signi cant impact on BLCA grade as shown in Figure 1D (p 0.001), stage as highlighted in Figure 1E (p 0.001), T classi cation as indicated by Figure 1F (p 0.001), and lymph node metastasis as shown in Figure 1G (p 0.001). These results, therefore, suggested that the abundance of CAFs supported the progression of BLCA. These results were further veri ed by analyzing the role of CAFs in the GEO BLCA cohort. Similar results were obtained since CAFs were highly abundant in the TME, showed a strong correlation with the ESTIMATE score (Figure 2A), and lowered the Overall Survival (OS) of BLCA patients ( Figure 2B). Moreover, CAFs were closely associated with the stage of BLCA in patients and lymph node metastasis ( Figure 2C), which strongly validated the ndings from TCGA. Although CAFs had no signi cant effect on distant metastasis in BLCA in both cohorts due to the limited number of M1 patients, the study still observed a trend in which CAFs promoted distant metastasis. Therefore, the above results suggested that the abundance of CAFs is a poor prognostic factor and enhances the progression of BLCA.
2. Identi cation of 74 hub genes related to CAFs as well as their enriched functions and pathways.
We further categorized patients into the high and low CAFs groups then screened for DEGs in the TCGA and GEO cohorts. A total of 555 and 187 genes were differentially expressed between the low and high CAFs groups in the TCGA and GEO cohorts. The heatmap shows the top 50 DEGs in Figures 3A and B. Moreover, WCGNA was applied to screen for modules that had the most signi cant association with levels of CAFs in both the TCGA ( Figure 3C) and GEO cohorts ( Figure 3D). The yellow module in the TCGA cohort showed the most signi cant association with a correlation level of 0.8, while the correlation between gene signi cance (GS) and module membership (MM) was 0.96. Similarly, the correlation between the blue module and the abundance of CAFs was shown to be 0.52, while GS and MM's correlation was 0.56.
Additionally, the intersection of DEGs and genes in the most related modules identi ed 74 hub genes ( Figure 4A). Go functional analysis and KEGG enrichment analysis indicated that these genes were crucial in functions related to remodeling of the extracellular matrix. Notably, the following GO terms were enriched; extra matrix organization, collagen-containing extracellular matrix and extracellular matrix structure constituent, et al. (Figure 4B). On the other hand, the following KEGG pathways were enriched; focal adhesion, ECM-receptor interaction, et al. (Figure 4C).
3. Identi cation of three key genes related to CAFs in BLCA.
Univariate cox regression was rst conducted on both the TCGA and GEO cohorts based on the expression of hub genes. The results showed that 10 and 22 genes, respectively, were signi cantly related to patients' survival with p values less than 0.05. The genes' p values and hazard ratios are shown in forest plots separately ( Figure 4D and 4E). The intersection of survival-related hub genes in TCGA and GEO identi ed CALD1, COL18A1 and TNC as the three key genes related to CAFs and further in uenced OS in BLCA ( Figure 4F). Notably, all these genes were signi cantly correlated with markers of CAFs, including; ACTA2 (α-SMA), MFAP5, MMP2, PDGFRB, VIM, FN1, FAP, FOXF1 and ZEB1 ( Figure 4G) (15).
Additionally, TNC was reported to be a biomarker of CAFs (16) and is a well known independent risk factor for BLCA (17). COL18A1 was previously reported to be involved in a 12-gene progression score signi cantly associated with progression (18). CALD1 was also de ned as a poor prognostic factor in BLCA (19). In the present study, we selected CALD1 for further analysis. GESA analysis through the hallmarks gene sets con rmed that CALD1 was positively involved in pathways related to epithelium to mesenchymal transition and hypoxia, which are crucial for inducing immunosuppression of the TME ( Figure 4H). Besides, GSEA of KEGG gene sets indicated that CALD1 was involved in multiple microenvironment remodeling pathways such as adhesion molecules cams, ECM receptor interaction and focal adhesion. It was also enriched in immune-related pathways, including the chemokine signaling pathway and cytokine-cytokine receptor interaction ( Figure 4I). 4. Correlation between CALD1, OS, and clinical characteristics in the TCGA BLCA cohort and its involvement in the modulation of the TME In the TCGA BLCA cohort, CALD1 was shown to markedly impact BLCA patients' OS since there was a signi cant difference between the high and low CALD1 expression groups (p 0.001). Additionally, the predictive value of CALD1 in cancer progression was con rmed through the ROC curve with an AUC of 0.679. Moreover, the expression levels of CALD1 differed signi cantly between different stages, T and N classi cations ( Figure 5A). Furthermore, the study observed a trend of increasing CALD1 levels with cancer metastasis, although no statistical signi cance was obtained. We further compared CAFs, macrophages and ESTIMATE scores between the high and low CALD1 expression groups. Results showed that the high CALD1 group had signi cantly higher CAFs, macrophages, stromal, immune, and ESTIMATE scores than the low CALD1 group (Figure 5B,5C). These results, therefore, indicated that CALD1 is a detrimental factor in the progression of BLCA. The ndings also con rmed that CALD1 was involved in modulating both stromal and immune microenvironment, which was possibly achieved through CAFs and macrophages. 5. Involvement of CALD1 in the regulation of TIICs and the immune checkpoint pathway.
The CIBERSORT algorithm was further used to validate the effect of CALD1 on TIICs in BLCA. The proportion of each TCGA BLCA patient's in ltrated immune subsets was analyzed using the CIBERSORT algorithm ( Figure 6A). Notably, correlation analysis showed that CALD1 was positively associated with macrophages (M0, M2) and negatively related to CD8+ T cells ( Figure 6C). A comparison of the TIICs levels between the high and low expression of CALD1 also con rmed an elevated level of macrophages (M0, M2) and decreased CD8 + T cells in the high CALD1 expression group ( Figure 6D). Consequently, the study further examined whether CALD1 was correlated with immune checkpoints such as PD-L1, which was also crucial in predicting immunotherapy e cacy in BLCA. Immune-checkpoint related genes, including CTLA-4, LGALS9 (GAL9), LAG-3, PDCD1 (PD-1), PDCD1LG2 (PD-L2), CD274 (PD-L1), TIGHT and HAVCR2 (TIM-3) were therefore selected for further analysis. Interestingly, almost all the genes (CTLA-4, LAG-3, PD1, PDL2, PDL1, TIGIT and TIM-3) were up-regulated in patients with high expression of CALD1 ( Figure 6E, F). These results, therefore, highlighted the role of CALD1 in regulating TIICs and immune checkpoint pathways.

Validation of the immune regulatory role of CALD1 in the GEO cohort
To validate the results from the TCGA cohort, we further analyzed the effect of CALD1 on BLCA patients in the GEO cohort. The results showed that high expression of CALD1 was signi cantly associated with a shorter OS. Moreover, the ROC curve revealed that CALD1 had an AUC of 0.730 in predicting localized BLCA progression to metastatic BLCA. Signi cant differences in the expression levels of CALD1 were also observed between stage as well as T and N classi cation ( Figure 7A). Moreover, up-regulated CAFs, macrophages and ESTIMATE scores were shown in the high CALD1 expression group compared to the other category ( Figure 7B, C). Furthermore, signi cantly higher CTLA4, LAG3, PDL1, PDL2 and TIM3 expression were observed in patients with high expression of CALD1 ( Figure 7D, E), further con rming the role of CALD1 in regulating immune checkpoints.
7. Expression of CALD1 in clinical specimens in the validation cohort. 40 BLCA patients with different grades, stage and TNM classi cation, were recruited to validate the above results. Expression levels of CALD1 were examined in pathological sections after clinical treatment with TURBT or radical cystectomy. The results revealed high expression levels of CALD1 in patients with a higher grade ( Figure 8A) and stage ( Figure 8B). Moreover, High co-expression was found between CALD1, ACTA2 and CD206 in the tumor stroma, especially in patients with advanced BLCA. These results con rmed the association of CALD1 with CAFs and macrophages, which may further lead to the progression of BLCA.

Discussion
BLCA is among the cancers characterized by in ltration of abundant immune cells in the TME, con rmed by BCG's treatment e ciency. Moreover, recent advances in immune checkpoint inhibitor therapy for BLCA further demonstrate that BLCA is profoundly regulated by tumor immunity (20).
Tumor cells and the microenvironment are a whole functional unit where the cells are regarded as seeds while the microenvironment in which they are located is considered the soil (21). Therefore, tumor cells and their microenvironment interact with each other and evolve together to promote tumors. Moreover, the TME is a complex integrated system divided into the immune microenvironment and non-immune microenvironment. Stromal components dominated the non-immune microenvironment (22), especially CAFs (23). There have been numerous studies on the crosstalk between stromal and immune components of the TME (24,25). Notably, macrophages are usually the most abundant TIICs in the tumor microenvironment, including in BLCA (26). Macrophages consist of two groups with different phenotypes, namely M1 and M2. Macrophages M2 are associated with immunosuppressive functions, angiogenesis, and the extracellular matrix's degradation, contributing to cancer migration and metastasis (27).
Numerous studies have demonstrated that the interaction between CAFs and macrophages can further promote the progression of cancer. For instance, Mazur et al. reported that CAFs could increase macrophages' adhesive ability and promote cancer invasion and metastasis (28). Additionally, Betul et al. showed that CAFs could recruit and differentiate monocytes into M2 macrophages and exert their immunosuppressive role through the PD-1 axis (29). The macrophages could also reversely modulate the status of CAFs. Zhang et al. demonstrated that macrophages could turn umbilical cord mesenchymal stem cells into CAFs, further promoting cancer progression through Epithelial-mesenchymal Transition (EMT) (30). The above evidence demonstrates that CAFs are involved in tumor immunomodulation and promoted tumor progression by interacting with macrophages. In the present study, CAFs were shown to be abundant in BLCA. Furthermore, our results also indicated that CAFs conferred signi cant adverse effects on the progression of BLCA and OS. We also con rmed a close relationship between CAFs and macrophages. Therefore, combined with existing evidence, the interaction between these two cell types can signi cantly promote the progression of BLCA. Targeting these two cell types may be a potential strategy for the treatment of BLCA.
On the other hand, we identi ed three key genes related to CAFs by WCGNA and take CALD1 for further study. Previous researches have demonstrated that CALD1 was involved in cell proliferation and migration through actin cytoskeleton recombination (19). It was also recognized as a tumor-speci c splicing variant in colon, bladder, and prostate tissue samples. The mis-splicing of CALD1 is an independent epigenetic event that is related to the destruction of the tight junctions between epithelial cells, hence altering the stiffness of the extracellular matrix and promoting cancer invasion and metastasis (31). Studies also con rmed CALD1 to be a risk factor for the progression of BLCA, but its role with regard to CAFs and immune regulation in BLCA is yet to be reported.
In this study, the vital role of CALD1 in the progression of BLCA was demonstrated in three independent cohorts. As a key gene associated with CAFs, CALD1 showed strong correlations with stromal and immune scores, suggesting its dual regulation of stromal and immune components. Also, CALD1 exhibited essential involvements in the processes modulating the TME, including EMT (32), hypoxia (33), extracellular matrix remodeling (34) and cytokine regulation (35), as revealed by GSEA. Moreover, CALD1 represented a positive correlation with M2 macrophages and a negative correlation with CD8+ T cells. Also, high correlations between the expression levels of CALD1 and multiple immune checkpoint genes were observed. These results highlighted the critical function of CALD1 in inducing immunosuppression in BLCA. The increased expression of CALD1 may be correlated with high expression of macrophage M2 and checkpoints and the depletion of T cell CD8. In the clinical validation cohort, we well con rmed the co-expression of CALD1, ACTA2 and CD206 through immunohistochemistry, which demonstrated the close correlation between CALD1, CAFs and macrophages. At the same time, in tumor specimens with different stages and depth of tumor invasion, we found signi cantly differentially expressed CALD1 level, which further indicated that CALD1 has the potential to be used as a marker of BLCA progression From our study, we con rmed that CALD1 is a risk factor in the progression of BLCA. It was also approved for the rst time that CALD1 regulated the tumor microenvironment associating with CAFs and macrophages. Also, our results clearly demonstrate the importance of bioinformatics analysis in cancer researches. Through the bioinformatics means such as WCGNA, CIBERSORT and MCP-counter, combined with clinical veri cation, we can get information that other means cannot achieve. It is believed that with the continuous progress of the algorithms, bioinformatics can provide more signi cant help for clinical diagnosis and treatment.
Despite the insightful ndings, limitations still exist in our study. First, relationships among CALD1, CAFs, macrophages and immunosuppression were only veri ed by correlation analysis. Further veri cation from in vitro and in vivo experiments are required for exploring the exact mechanisms. Second, the potential of CALD1 to become a speci c marker of CAFs in BLCA still need validation since we did not discuss its expression in normal broblasts. Last, although signi cant co-expression was found between CALD1, ACTA2 and CD206 in BLCA sections, a more extensive validation cohort is still necessary to avoid the selection bias.

Conclusion
In conclusion, this study con rmed the pro-tumor function of CAFs and identi ed CAFs-related genes in BLCA through WCGNA and screening for DEGs in both the TCGA and GEO cohorts. Moreover, CALD1 was recognized as one of the key genes related to CAFs and outcomes in BLCA. Further analyses showed that CALD1 played a vital role in regulating the TME of BLCA. Furthermore, CIBERSORT and correlation analysis con rmed that CALD1 was related to the in ltration level of multiple immune cells in the TME, especially macrophages M2 and CD8 T cells. It was also shown that high expression of CALD1 might lead to an increased level of immune checkpoint-related genes, including PDL1. Therefore, CALD1 may be associated with the immunosuppression status of TME in BLCA, which further leads to tumor progression. Further studies on CALD1 may provide insights into the immune network in BLCA and offer new targets for cancer treatment.

Declarations
Ethics approval and consent to participate 40 BLCA spacemen were collected from Shanghai First People's Hospital, with patients' informed consent and the approval of the Medical Ethics Committee of Shanghai First People's Hospital (approval number: 2020KY143 ). CALD1 regulated the TME and promoted the progression of BLCA A. Kaplan-Meier curves of OS between high and low CALD1 group using Log-Rank test in TCGA BLCA cohort. ROC curve demonstrated the accuracy of CALD1 in predicting cancer progression, with an AUC of 0.679. The expression of CALD1 was signi cantly elevated with cancer progression. Three different bar plots showed the expression level of CALD1 in different tumor stages, T and N classi cation, respectively B. CAFs (p<0.001), macrophages (p<0.001), stromal score (p<0.001), immune score (p<0.001) and ESTIMATE score (p<0.001) in the CALD1 high expression group were signi cantly higher than those in the CALD1 high expression group C. Close relationships of CALD1 with CAFs (R=0.88,p<0.001), Macrophages (R=0.58,p<0.001), stromal score (R=0.84,p<0.001), immune score (R=0.51,p<0.001) and ESTIMATE score (R=0.72,p<0.001) were observed in the TCGA cohort. Validation of the TME regulating, Oncogenic promoting and immune checkpoint associating role of CALD1 in the GEO cohort. A. Kaplan-Meier curve validated the difference of OS between high and low CALD1 expression BLCA Patients in the GEO cohort, with a P-value <0.001. The AUC value of the ROC curve for CALD1 prediction of tumor progression in the GEO cohort was 0.73. A signi cantly higher level of CALD1 was also detected in patients with higher grade, stage, T and N classi cation in the GEO cohort. B-C. CAFs (p<0.001), macrophage (p<0.001), stromal score (p<0.001), immune score (p<0.001) and ESTIMATE score (p<0.001) were signi cantly different between the high and low CALD1 expression