CAF-mediated proliferation and migration of PCa cells

Primary NFs and CAFs derived from 3 PCa patients were isolated, and their identity was verified by western blot using their specific biomarkers (α-SMA, PDGFRα/β and FAP) (Fig. 1A). The culture supernatant of CAFs and NFs were collected for cellular experiments. As shown in Fig. 1B and C, the CAF-derived supernatant significantly facilitated PC-3 and DU145 cell migration (Fig. 1B) and proliferation (Fig. 1C).

Fig. 1

CAF heterogeneity and function in prostate cancer (PCa). (A) Western blot analysis for FAP, PDGFRα/β, α-SMA and GAPDH in the NFs (Normal fibroblasts) and CAFs (cancer-associated fibroblasts). (B) Transwell assay of PC-3 and DU145 cell migration mediated by CAF culture supernatant. (C) CCK-8 assessment of PC-3 and DU145 cell proliferation mediated by CAF culture supernatant. (D) Heatmap of differentially expressed genes (DEGs) in paired NFs and CAFs. (E) T-SNE plot of single cells from CAFs in PCa tissues. (F) Heatmap shows the marker genes of s distributed in the 12 clusters.

CAF heterogeneity in PCa tumors

Given the findings above, we sought to probe potential mechanisms underlying CAF-driven PCa cell migration/proliferation. Towards this end, we first analyzed transcriptomic profiles of NFs and CAFs in GSE85606 and GSE68164 PCa cohorts to identify differentially expressed genes (DEGs) between them. A total of 43 and 63 DEGs (|LogFC| > 1 and P < 0.05) were found in GSE85606 and GSE68164 cohorts, respectively (Fig. 1D, and Tables S1 and S2). There were only 4 overlapping DEGs (KRT7, IGFBP2, CPXM2 and TINAGL1) in both cohorts, among which KRT7 expression showed opposite trends. We further analyzed scRNA-seq data of CAFs and observed that CAFs could be divided into 11 clusters (Fig. 1E). Each of these 11 clusters contained unique top DEGs (Fig. 1F and Table S3). The differences in the DEGs among each cluster demonstrate the heterogeneity and plasticity of CAFs (Figure S1). Likely, those identified DEGs render CAFs stimulatory effects on PCa cell prolifreation and migration. Alternatively, the DEGs mark oncogenic CAF subpopulations.

Establishment of the CAF subtype score to predict patient PFI

To gain insights into CAF-driven PCa aggressiveness in more depth and more broadly, we further focused on the molecularly classified CAF subtypes. Based on the molecular heterogeneity of CAF populations obtained from scRNA-seq in solid tumors, CAFs have recently been stratified into the following 6 categories [22]: CAFmyo, CAFinfla, CAFadi, CAFendMT, CAFpn, and CAFap (Fig. 2A). To determine the effect of each CAF subtype on PCa progression and outcomes, we applied this CAF classification system to the bulk RNA-seq profiled PCa tumors (TCGA cohort) by using the top 30 expressed genes in each CAF subtype (Table S4), and CAFs in these PCa tumors were successfully categorized into the identical 6 subtypes, too. Univariate Cox regression analysis was first used to examine the association between patient PFI and the expression levels of top 30 genes in each CAF subtype (Figure S2A), but not all those 30 genes could predict PFI, indicating the role of the subtype rather than gene expression as per. Nevertheless, the CAF scores of each category were constructed according to the expression level of those 30 top genes. Using the median score as the cutoff, the CAF score was associated with PFI in all subtypes (Fig. 2B and C). Among the 6 different CFA categories, the area under the ROC curve (AUC) of the scores in all 6 categories was the largest in 7 years (Figure S2B). We further identified the top 10 DEGs in each subtype (Fig. 2D), and subsequent GSEA analyses showed both different and overlapping pathway enrichments among 6 subtypes with high CAF scores (Fig. 2E). The enriched pathways in the high-score groups mainly include cell proliferation, ECM, EMT, angiogenesis, inflammation and immune responses, which are intimately associated with tumor progression. The co-expression network between CAF typing and CAF-related genes was analyzed by Sankey plot (Figure S3A). The regulatory network of CAF-related genes in different subtypes represents the study of expression correlation and tumor progression in PCa patients (Figure S3B).

Fig. 2

The CAF subtype classification and scores for progression prediction in PCa. (A) The classification of CAF subtypes. (B) Univariate Cox regression analysis of CAF subtype scores and association with progression-free interval (PFI). (C) Kaplan-Meier analysis of the CAF subtype scores and association with PFI. (D) Heatmap of top 10 differentially expressed genes in CAF subtype score low and high groups. (E) GSEA results showing the activated signaling pathways in the CAF score high group. (F-G) The validation of the CAF subtype score model to predict PFI, as determined by Kaplan-Meier analysis in MSKCC and CPGEA PCa cohorts

We further analyzed PCa cohorts from the MSKCC, CPGEA and GEO datasets to validate the CAF score model as a prognostic factor observed in the TCGA PCa patients. For the MSKCC cohort, similar results were obtained. In the six subtypes of CAF score models, patients in the high score group had more rapid disease progression than those in the low score one (Fig. 2F). In the CPGEA cohort, patients in the high score group had shorter PFI, but statistical significances were reached only for CAFmyo, CAFendMT and CAFap subtypes, while at a board-line for CAFpn (Fig. 2G). The GSE70770 cohort analysis showed that the scores for CAFmyo, CAFendMT and CAFap subtypes were significantly associated with patient PFI (Figure S6A), as observed in the CPGEA cohort.

The total CAF score model as a predictor for PCa patient PFI and treatment response

To simplify the CAF subtype score system above for potential clinical application, we integrated six subtype scores and all CAF associated genes to construct a total CAF (tCAF) score model. Taking the median tCAF score as the cutoff, the analysis of TCGA, MSKCC, CPGEA and GSE70770 PCa cohorts showed that tCAF score had high accuracy in predicting PFI. The patient PFI in the tCAF high group was significantly worse than that in the low one (Figs. 3A and S6B). The ROC curves for each group, when the third, fifth, and seventh years were evaluated as the end time points, demonstrated the robust predictive power of the tCAF score model (Figs. 3B and S6C). GSVA pathway analysis unraveled that proliferation-related pathways were highly enriched in the tCAF score high group (Fig. 3C). Further assessments of the TCGA PCa cohort showed that there were significant differences between the tCAF score and age, Gleason score, T and N stages (Fig. 3D and E). The Gleason score, an important indicator in PCa, was significantly higher in the tCAF score high group. Next, we explored whether the tCAF score could be used in the selection of drug therapy (excluding ADT) in PCa patients. The IC50 value of each drug for each patient in the tCAF score low and high groups was calculated using the oncoppredict package. We computationally identified 14 drugs that were more effective in the tCAF score low group and 46 drugs that were more effective in the tCAF score high group (Fig. 3F and Figure S4). The 3D structural tomography of talazoparib, zoledronate, cediranib, gemcitabine, and savolitinib that could potentially be used to treat patients in the tCAF high group was searched in PubChem database (Fig. 3G).

Fig. 3

The tCAF score model and association with survival, clinical features and drug sensitivity in PCa. (A) Kaplan-Meier analysis of the tCAF score low and high groups in TCGA, MSKCC and CPGEA PCa cohorts. (B) ROC curves of tCAF scores at 3, 5, and 7 years in the TCGA cohorts. (C) GSVA enrichment analysis showing the activation states of biological pathways in the tCAF score low and high groups. (D-E) Differences in clinical characteristics between the tCAF score low and high groups in the TCGA PCa cohort. (F) The effect of tCAF scores on response to commonly used drugs. (G) The 3D structure tomographs of 5 candidate small-molecule drugs for tCAF score high groups in PCa.

Correlation of tCAF scores with genomic alterations in PCa tumors

We calculated the tumor mutation burden (TMB) for each patient in the TCGA PCa cohort. The TMB in the tCAF high score group was significantly higher (p < 0.001) (Fig. 4A). There was a highly positive correlation between TMB and tCAF score (Fig. 4B). The overall TMB was 53.19% and 69.42% in the tCAF score low and high groups, respectively (Fig. 4D). The prognosis of patients with both tCAF score and TMB low was much better (Fig. 4C). Figure 4E showed the mutual exclusivity and co-occurrence of mutations in tCAF score groups. We further examined the mutation frequencies of nine major oncogenic pathways in the tCAF score low and high groups. Nine major oncogenic pathways were detected in the tCAF score low group, while 10 major oncogenic pathways were detected in the tCAF score high group, mainly including RTK-RAS, WNT, NOTCH and Hippo pathway (Fig. 4F). Cancers differ from each other in their mutational patterns. We examined these differentially mutated genes by comparing the two cohorts of CAF score. The results showed that besides FLG2, TP53, NALCN, SACS, PTEN, OBSCN, RYR1 and FOXA1 were highly mutated in the CAF score high group (Fig. 4G). CNV alterations (mainly copy number deletions) occurred more frequently in all the CAF score subtype high groups (Figure S5).

Fig. 4

The mutational profiles in tCAF score low and high PCa tumors. (A-B) The relationship between TMB and tCAF scores. (C) Kaplan-Meier analysis of the tCAF score and TMB. (D) The frequencies of mutated genes among the tCAF score low and high groups. (E) Mutual exclusivity and co-occurrence of mutations in tCAF score low and high groups. (F) The mutation frequencies of common oncogenic pathways in two tCAF score subtypes. (G) Differentially mutated genes in two tCAF score subtypes

Effect of the tCAF score on immune status in PCa tumors

The PCa tumors in the TCGA cohort were scored using ssGSEA to quantify the activity, enrichment level and function of immune cells in each sample, and then grouped according to their tCAF scores (Fig. 5A). The immune status was more active in the tCAF high score group. Based on the expression profile, the ESTIMATE algorithm was used to calculate the stromal, immune and ESTIMATE scores of PCa. The results showed that the ESTIMATE, immune and stromal scores in the tCAF score high group were all higher than those in the low one (Fig. 5B). The cytolytic activity (CYT) score, an immunotherapy biomarker characterizing the antitumor immune activity of CD8+ cytotoxic T cells and macrophages, was significantly higher in the tCFA score high group (Fig. 5D). These results indicate the active immune status in the high tCAF score tumors. IPS is a quantitative index to evaluate the cancer-immunity cycle (CIC) efficacy. In the case of CTLA4 expression, the tCAF score low group had a better response to immunotherapy, while there was no difference between the two groups in the case of PD-L1 expression (Fig. 5C). Consistently, CTLA4 was significantly expressed in the tCAF score high group, while NECTIN2 was significantly expressed in the low group. The expression of PD-L1, PD-L2 and CCA did not differ significantly between the two groups (Fig. 5E). We further examined the relationship between tCAF scores and major histocompatibility complex (MHC). Except for TNFRSF14 and CD28, the expression level of MHC gene sets tended to be higher in the tCAF high score group (Fig. 5F).

Fig. 5

The tCAF score prediction of response to immunotherapy in PCa. (A) Twenty-nine immune-related gene sets were enriched in TCGA PCa cohort. (B) The stromal, immune and ESTIMATE scores between two tCAF score subtypes. (C-D) The IPS and CYT scores between two tCAF score subtypes. (E-F) The expression of immune checkpoint related genes and MHC gene set in two tCAF score groups

Association between the tCAF score and immune cell infiltration

Antitumor immunity in tumor tissue can be interpreted as seven sequential processes, including release of cancer antigens (step 1), cancer antigen presentation (step 2), priming and activation (step 3), Tracking of immune cells to tumors (step 4), infiltration of immune cell into tumors (step 5), recognition of cancer cells by T cells (step 6), and killing of cancer cells (step 7). Although only step1 and step5 showed active status in the tCAF score high group, step2, step3, step4, step6 and step7 showed similar active status in both groups (Fig. 6A and B). Further analysis of infiltrated immune cells in tumor tissues showed that CD4 memory, CD8 effector, CD8 naive, B cells, NK cells and DC cells were more abundant in the tCAF score high group than the tCAF score low group (Fig. 6C and D). However, the degree of infiltration of Th cells, CD8 memory, Monocytes CD16, and pDC cells was reversed (Fig. 6D). Based on the immunological classification of solid tumors by Thorsson et al. [24], we further examined distributions of immune subtypes in two tCAF score groups. In the tCAF score low group, C1 (wound healing), C2 (IFN-gamma dominant), C3 (inflammatory), and C4 (lymphocyte depleted) accounted for 5%, 2%, 83%, and 10%, respectively. However, in the tCAF score high group, C1, C2, C3 and C4 accounted for 13%, 7%, 69% and 11%, respectively (Fig. 6E). The proportion of C1 and C2 was significantly higher, while the proportion of C3 was lower in the tCAF score high group.

Fig. 6

Immune cell infiltration in CAF score low and high PCa tumors. (A) Tumor immune cycle analysis. (B) Heatmap of the seven-step cancer immunity cycle and tCAF score. (C-D) Proportional plot of immune cell infiltration in two tCAF score groups. (E) The immunological classification of solid tumors with tCAF score

Shorter telomeres in the tCAF score high PCa tumors and association with unfavorable PFI

Telomere shortening or dysfunction occurs with aging, which drives inflammation [26]. We thus sought to determine whether increased immune activity observed in the tCAF score high tumors above was associated with altered telomere length. To this end, the TCGA PCa cohort was analyzed [25]. As shown in Fig. 7A and G, PCa tumors with high tCAF score had dramatically shorter telomeres compared to low score tumors. Moreover, the worst PFI was observed in patients with shortest telomere-bearing tumors (Fig. 7H).

Fig. 7

The relationship between telomere length and tCAF scores in PCa tumors. (A-G) The telomeres length in different CAF score subtypes. (H) Kaplan-Meier analysis of the telomere length in the TCGA PCa cohort

The tCAF score system as a prognostic factor in pan-cancer

Because Luo et al. showed similarity in CAF heterogeneity and transcriptomic profile across cancer types, we determined whether our tCAF score system could predict PFI and immune status in other solid tumors by analyzing the TCGA pan-cancer. We established separate tCAF scoring models for all solid tumors to improve the accuracy of CAF scoring model adaptation. Individual tCAF scores were established by performing univariate Cox regression analyses in each solid tumor (Fig. 8A and B). tCAF scores predicted PFI in pan-cancer (p < 0.05) (Fig. 8B). MK plots further showed significantly shorter PFI in those tumors with high tCAF score (Fig. 8D). TMB score and CYT score were analyzed simutaneously. According to the tCAF scores, TMB showed significant differences in BRCA, CESC, HNSC, KIRC, KIRP, LGG, LIHC, LUAD, PAAD, STAD and THCA (Fig. 8C). CYT score was more active in BLCA, COAD, GBM, KIRP, LGG, LIHC, LUSC, OV and STAD with tCAF score high tumors, while in BRCA, HNSC, PAAD, SKCM and UCEC with tCAF score low tumors (Fig. 8E).

Fig. 8

The tCAF score model in pan-cancer. (A) The CAF score in solid tumor. (B) Univariate Cox regression analysis of CAF score in pan-cancer. (C) The TMB of CAF score low and high groups in pan-cancer. (D) Kaplan-Meier analysis of the CAF score low and high groups in pan-cancer. (E) The CYT score of CAF score low and high groups in pan-cancer