1. Abundance of CAFs is a poor prognostic factor associated with the progression of BLCA.
We first 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 significant 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 significantly 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 significant impact on BLCA grade as shown in Figure 1D (p<0.001), stage as highlighted in Figure 1E (p<0.001), T classification 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 verified 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 findings from TCGA. Although CAFs had no significant 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. Identification 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 significant 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 significant association with a correlation level of 0.8, while the correlation between gene significance (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 identified 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. Identification of three key genes related to CAFs in BLCA.
Univariate cox regression was first 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 significantly 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 identified CALD1, COL18A1 and TNC as the three key genes related to CAFs and further influenced OS in BLCA (Figure 4F). Notably, all these genes were significantly 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 significantly associated with progression(18). CALD1 was also defined 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 confirmed 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 significant difference between the high and low CALD1 expression groups (p<0.001). Additionally, the predictive value of CALD1 in cancer progression was confirmed through the ROC curve with an AUC of 0.679. Moreover, the expression levels of CALD1 differed significantly between different stages, T and N classifications (Figure 5A). Furthermore, the study observed a trend of increasing CALD1 levels with cancer metastasis, although no statistical significance 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 significantly 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 findings also confirmed 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 infiltrated 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 confirmed 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 efficacy 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.
6. 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 significantly 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. Significant differences in the expression levels of CALD1 were also observed between stage as well as T and N classification (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, significantly higher CTLA4, LAG3, PDL1, PDL2 and TIM3 expression were observed in patients with high expression of CALD1 (Figure 7D, E), further confirming 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 classification, 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 confirmed the association of CALD1 with CAFs and macrophages, which may further lead to the progression of BLCA.