We collected the scRNA-seq data for 26 hormonal-sensitive prostate cancer (HSPC) samples and 8 CRPC samples from 5 public datasets to study CAFs heterogeneity in prostate cancer ( Table S1 ) (16–18, 23). First, we used Seurat V3 R package to import the raw counts data of each sample into R software (24). Next, we removed out cells from the raw counts data which met the following exclusion criteria: cells with less than 500 genes expressed, cells with more than 15% of genes derived from mitochondrial genes, or cells with more than 15% of genes derived from ribosomal genes. Since genes associated with mitochondria and ribosome can cause unexpected noise and such genes were not the focus of our study, we removed out the mitochondrial genes (genes beginning with MT-) and ribosomal genes (genes beginning with RPL or RPS) from the raw counts data. The DoubletFinder R package was utilized to remove the doublets from the raw counts data (25). In the remaining high-quality cells, Harmony R package was used to integrate multiple scRNA-seq datasets. After integration, the reduction was set as “harmony” for subsequent dimensionality reduction. We normalized and scaled the raw counts data using the NormalizeData and SclaeData function, respectively. To reduce dimensionality, we performed principal component analysis on the normalized data using the top 2000 highly variable features determined by the FindVariableFeatures function. The appropriate number of principal components was selected using the ElbowPlot function. The t-distributed stochastic neighbor embedding (t-SNE) and uniform manifold approximation and projection dimensionality reduction (UMAP) were then conducted by RunTSNE and RunUMAP function, respectively. The FindClusters function was implemented to identify potential cell clusters.

The scRNA-seq data for 7 bone metastatic prostate cancer samples and 8 normal bone samples was collected from GSE143791 (26). The processes of data quality control, multiple datasets integration, and standard data processing were consistent with the above.

We annotated the cell clusters from the prostate cancer samples by the average expression of the following well-recognized cell type markers: T cells (PTPRC, CD3D, CD3E, and CD3G), B cells (CD79A, IGKC, IGLC2, and IGLC3), macrophages (CD68, CD163, FCGR3A, LYZ, and CSF1R), mast cells (ENPP3, KIT, SLC18A2, and MS4A2), fibroblasts (αSMA, FN1, and FAP), endothelial cells (PECAM1, ENG, CDH5, and VWF), and epithelial cells (EPCAM, AR, KRT5, KRT14, KRT8, and KRT18). For the annotation of cell clusters from bone metastatic prostate cancer samples, the cell type marker genes were obtained from Youmna et al. (26).

The FindMarkers function embedded in Seurat V3 R package was used to identify the differentially expressed genes (DEGs) between different cell clusters or cell types. Genes that met the following criteria were defined as DEGs: genes expressed in at least 20% of cells, the adjusted P value ≤ 0.01, and the |Fold Change| ≥ 2. Function enrichment analysis of DEGs was performed using the Metascape online tool (http://metascape.org).

We used monocle2 R package to investigate the cell lineage trajectory of CAFs (27). The dispersion Table function was used to identify significant genes that meet the criteria of mean_expression≥0.5 and dispersion_empirical≥1*dispersion_fit. Cells were then ordered using orderCells function. We used DDRTree function for reducing dimension and plot_cell_trajectory function for visualization. After the cell lineage trajectory was constructed, we used the differentialGeneTest function to identify DEGs along the pseudotime.

We used gene set variation analysis function embedded in gsva R package to calculate the CAF-C0 signature score and CAF-C1 signature score (28). The detailed gene lists associated with the above scores were described in Table S2 .

We used pySCENIC to identify the potential gene regulatory networks in the scRNA-seq data (29). The regulons activities of transcriptional factors were calculated using the AUCell function embedded in pySCENIC. The analysis was performed with the default parameters. The input expression matrix was prefiltered using the following criteria: cells with more than 500 genes expressed and genes expressed in at least 5 cells.

To explore potential cell-cell interactions between cells, we used CellPhoneDB to identified significant receptor-ligand interactions (30). The analysis was performed with default parameters. Valid receptor-ligand interactions were defined as the expression of specific receptors in one cell type and the expression of the corresponding ligands in another cell type.

The RNA-seq data and clinical information for TCGA-PRAD, TCGA-BLCA, TCGA-KIRC, TCGA-STAD, TCGA-LGG, TCGA-MESO, and TCGA-SKCM dataset were obtained from UCSC Xena (http://xena.ucsc.edu) (31). The gene expression matrix and clinical information for GSE116918 dataset were collected from the GEO database (32). The RNA-seq data and clinical information for metastatic urothelial carcinoma patients treated with anti-PD-L1 immunotherapy were obtained from http://research-pub.gene.com/IMvigor210CoreBiologies (33). The RNA-seq data and clinical information for metastatic melanoma patients treated with anti-PD1 immunotherapy were obtained from GEO database GSE78220 (34). The gene expression matrix for GSE32269, GSE2443, GSE6811, GSE31410, GSE21887, and GSE61379 was obtained from GEO database (35–40).

Gene Set Enrichment Analysis (GSEA) was applied using the GSEA software (Version: 3.0; http://software.broadinstitute.org/gsea/index.jsp). Significance was set at NES>1.0, P value<0.05, and FDR<0.25.

We performed hierarchical clustering on the normalized expression data of CAFs-C0 and CAFs-C1 signature components. Euclidean distance was chosen as the distance metric and Ward as the clustering method.

We collected 56 treatment-naive HSPC samples and 15 post-endocrine therapy CRPC samples from the Renji Hospital. Formalin-fixed and paraffin-embedded tumor tissues were collected. For antigen retrieval, 3 μm paraffin-embedded sections were unmasked in 1 x Tris-EDTA buffer for 20 minutes at 95°C. For immunohistochemistry staining, slides were incubated with αSMA (1:5000; Abcam; ab124964), FN1 (1:2500; Abcam; ab268020), or FAP (1:200; Abcam; ab207178) overnight at 4°C and then incubated with HRP-linked secondary antibody (1:500; Cell Signaling Technology; 7074S) for 1 hour at room temperature. DAB was used for visualization. Immunofluorescence staining was performed in 3 μm paraffin-embedded sections using the Three-color Fluorescence kit (Recordbio Biological Technology, Shanghai, China) based on the tyramide signal amplification technology according to the manufacture’s instruction. The antibody catalog numbers and staining concentrations used in immunofluorescence staining were: αSMA (1:5000; Abcam; ab124964), CAV1 (1:250; Cell Signaling Technology; 3238S), FN1 (1:2500; Abcam; ab268020), FAP (1:200; Abcam; ab207178), and STAT2 (1:200; Abcam; ab32367).

All statistics and graphs were completed using the R software (v3.6.3) and GraphPad Prism software (v8.0). We used two-tailed t-test or Wilcoxon rank-sum test to calculated the P value. Kaplan–Meier survival analyses were completed using survival and survminer R packages. Data were represented as mean values ± SD. Significant P value was set at <0.05.

To identify CAFs subtypes in prostate cancer microenvironment, we collected scRNA-seq data for 26 HSPC samples and 8 CRPC from 5 public datasets ( Table S1 ). After data quality control, multiple datasets integration, and standard data processing, we discovered seven cell types in the prostate cancer microenvironment, including epithelial cells, T cells, B cells, macrophages, mast cells, fibroblasts, and endothelial cells ( Figures S1A–C ). Reclustering of 8,316 fibroblasts generated eight CAFs subtypes ( Figure 1A ). We next assessed the expression of well-established fibroblasts marker genes among the eight subtypes (C0-C7). Several genes (VIM and S100A4) were highly expressed in all eight CAFs subtypes, whereas other genes (FN1, FAP, αSMA, CAV1, PDGFRA, and PDGFRB) were nonuniformly expressed among the eight CAFs subtypes ( Figure 1B ). As shown in Figures 1C, D , these eight CAFs subtypes could be distinguished from each other by the expression of subtype-specific marker genes. Specially, CAFs-C0 highly expressed αSMA, CAV1, and microvasculature signature genes such as MYH1, MCAM, and RGS5, while CAFs-C1 highly expressed FN1, FAP, and collagen molecules such as COL1A1, COL3A1, and COL1A2. Pathway enrichment analysis on specific marker genes of each CAFs subtype revealed the functional heterogeneities among CAFs subtypes ( Figures 1E , S1D ). Pathways related to microvasculature development were enriched in CAFs-C0, while CAFs-C1 had high levels of extracellular matrix organization signatures. Like CAFs-C1, the path-ways enriched in CAFs-C2 were also related to extracellular matrix organization. CAFs-C3 had significant enrichment of inflammatory response regulation signatures. Since scRNA-seq analysis revealed that αSMA⁺ CAV1⁺ and FN1⁺ FAP⁺ CAFs were the two most abundant subpopulations of CAFs in prostate cancer microenvironment, we thus stained prostate cancer samples with these four fibroblasts markers to verified our findings from the bioinformatic analysis. The staining results confirmed the co-expression of αSMA and CAV1 and the co-expression of FN1 and FAP ( Figure 1F ), with minimal overlap between αSMA and FN1 expression ( Figure 1G ). Recent scRNA-seq analysis of different tumors have revealed various CAF subtypes, notably myofibroblastic CAFs (myCAFs) and inflammatory CAFs (iCAFs). myCAFs generally highly express collagen molecules and αSMA has been widely recognized as the marker for myCAFs (41–44). However, our analysis revealed that FN1⁺ FAP⁺ CAFs-C1, but not αSMA⁺ CAV1⁺ CAFs-C0, had high expression levels of collagen molecules and contributed to extracellular matrix remodeling. Of note, while some studies, including our analysis, have revealed that αSMA⁺ and FAP⁺ CAFs are distinct subpopulations of CAFs (7, 42), other studies have reported the co-expression of αSMA and FAP on CAFs, which suggests the considerable heterogeneity among tumors (43, 45, 46). iCAFs, characterized by cytokine secretion, have been reported to be involved in tumor progression and treatment resistance (8, 11), but our analysis showed that the two most abundant CAFs subtypes in prostate cancer microenvironment did not have the feature of cytokine secretion.

The subtype composition of CAFs changed with disease progression: CAFs-C0 as the most abundant subtype in HSPC but CAFs-C1 as the most abundant subtype in CRPC ( Figure 2A ). To further investigate the dynamic transition of compositional and transcriptional profiles in CAFs during disease progression, we performed trajectory analysis of CAFs using Monocle2. Trajectory analysis revealed that CAFs-C0 were at the early stage of the trajectory path, whereas CAFs-C1 were at the late stage of the trajectory path ( Figure 2B ). We then analyzed the dynamics of CAFs-C0 and CAFs-C1 signature scores along the pseudotime. We defined these signatures by obtaining the overlapping genes between upregulated DEGs of fibroblasts and upregulated DEGs of CAFs-C0 or CAFs-C1 ( Figures S2A, B ; Tables S2, S3 ), only retaining genes specific for both fibroblasts and CAFs subtypes. Consistently, CAFs-C0 signature scores were downregulated, while CAFs-C1 signature scores were upregulated along the pseudo-time ( Figure 2C ). We also observed significant downregulation of CAFs-C0 marker genes (αSMA and CAV1) but significant upregulation of CAFs-C1 marker genes (FN1 and FAP) during these transitions ( Figure 2D ). We next analyzed the transcriptional profile changes associated with CAFs state transitions and divided the transcriptional profile changes into 3 phases ( Figure 2E ). CAFs-C0 and CAFs-C3 were phase 1 cells, CAFs-C1 and CAFs-C4 were phase 3 cells, and other CAFs subtypes were phase 2 cells ( Figure 2E ). Phase 1 cells had high expression levels of αSMA, CAV1, PDGFRB, TLN1, and FLNA, matching the specific marker genes of CAFs-C0 and CAFs-C3 ( Figure 2E ; Table S2 ). Phase 3 cells highly expressed extracellular matrix markers such as DCN, FBLN1, FN1, COL1A1, and COL1A2, matching the specific marker genes of CAFs-C1 ( Figure 2E ; Table S2 ). Pathway analysis suggested that phase 1 cells were implicated in vasculature development, whereas phase 3 cells were involved in extra-cellular matrix remodeling. These results indicated that CAFs had different properties along the trajectory path, with a CAFs-C0 phenotype in the early stage but a CAFs-C1 phenotype in the late stage. Furthermore, GSEA analysis reveals that CAF-C0 signature was enriched in HSPC samples and CAF-C1 signature was enriched in CRPC samples ( Figure 2F ).

To investigate the regulatory programs which may lead to the substantial differences in transcriptional profile among CAFs subtypes, we investigated the transcriptional regulatory networks for CAFs. pySCENIC inferred the top five activating transcriptional factors for each CAFs subtype ( Figure 3A ). Among these transcriptional factors, both STAT2 and PRRX2 had high expression levels and transcription activities in CAFs-C1 ( Figures 3A, B , S3A ). However, only the downstream target genes of STAT2 contained the specific marker genes for CAFs-C1 ( Figures 3C , S3B ). Therefore, STAT2, but not PRRX2, may be involved in the regulation of CAFs-C1 differentiation. To exclude the influences of other STAT family members in CAFs-C1, we examined the expression levels and transcription activities of other STAT family members in CAFs-C1. The expression levels and transcription activities of other STAT family members were considerably lower than those of STAT2 in CAFs-C1 ( Figures S3C, D ). Subsequently, we evaluated the relationship between STAT2 expression and CAFs subtypes by immunofluorescence staining. The staining results suggested that CAFs-C1 exhibited high expression of STAT2, whereas CAFs-C0 lacked expression of STAT2 ( Figures 3D, E ). Additionally, by mining the prostate cancer dataset from TCGA, we observed a significant positive correlation between the expression levels of STAT2 and CAFs-C1 marker FN1, while no significant correlation was observed between the expression levels of STAT2 and CAFs-C0 marker αSMA ( Figure S4 ). These results suggested that STAT2 may participated in the differentiation of CAFs-C1.

Using CellPhoneDB, we investigated the cell–cell communication networks between CAFs and other cell types. CAFs showed frequent interactions with other cell types and endothelial cells were the predominant cell type interacting with CAFs ( Figure 3F ). We next focused on the cell–cell interactions between CAFs-C0, CAFs-C1, and other cell types. The most frequent cell-cell interactions between CAFs-C0, CAFs-C1, and other cell types were related to extracellular matrix organization ( Figure 3G ).

Bone metastases are the primary cause of mortality in prostate cancer patients (47–49). Osteoblasts are responsible for the formation of bone matrix by secreting and organizing collagens and other proteins. A large body of evidence implies that prostate cancer cells frequently interact with osteoblasts in the bone metastatic microenvironment, which are crucial for the metastatic colonization and proliferation of prostate cancer cells (50, 51). To investigate the osteoblasts landscape in the bone metastatic microenvironment, we collected scRNA-seq data for 7 bone metastatic prostate cancer samples. Thirteen major cell types were identified in the bone meta-static microenvironment ( Figures S5A–C ). Reclustering of osteoblasts produced two subtypes, termed osteoblasts-C0 and osteoblasts-C1 subtype ( Figure 4A ). Osteoblasts-C0 and osteoblasts-C1 highly expressed specific marker genes of CAFs-C0 and CAFs-C1, respectively ( Figure 4B ). CAFs-C0 signature and CAFs-C1 signature scores exhibited non-uniform distribution among osteoblasts ( Figure 4C ). We next sought to identify the DEGs between osteoblasts-C1 and osteoblasts-C0 ( Figure 4D ). The upregulated genes of osteoblasts-C0 contained lots of CAFs-C0 signature genes, while osteoblasts-C1 was characterized by increased expression of CAFs-C1 signature genes. Pathway analysis on the DEGs of osteoblasts-C0 and osteoblasts-C1 showed that osteoblasts-C0 were related to blood vessel development and osteoblasts-C1 were involved in extracellular matrix remodeling ( Figure 4E ). Additionally, the osteoblasts-C1 was the most abundant subtype of osteoblasts in the bone metastatic microenvironment ( Figure 4F ). We also analyzed the scRNA-seq data for 8 normal bone samples and found that the distribution of osteoblasts was absent in normal bone samples ( Figure S6 ). These results revealed that osteoblasts in the bone metastatic lesions and CAFs in the primary tumor sites demonstrated similar phenotypic and functional properties, which suggested that CAFs-C0 and CAFs-C1 signatures were the two most predominant phenotypic signatures of CAFs and osteoblasts in prostate cancer. We also sought to investigated the cell–cell interaction networks between osteoblasts and other cell types. In the bone metastatic microenvironment, endothelial cells were the predominant cell type interacting with osteoblasts ( Figure 4G ), consistent with our earlier observations in the primary tumor microenvironment.

To further explore the role of CAFs heterogeneity in CRPC development, we extended our analysis to bulk-sequencing profiles. GSEA analyses revealed significant enrichment of CAFs-C1 signature in CRPC samples compared with HSPC samples ( Figure 5A ). However, there was no difference in enrichment of CAFs-C0 signature between HSPC and CRPC samples ( Figure S7A ). We next examined the protein expression of αSMA and FAP in HSPC and CRPC samples. Immunohistochemistry staining suggested that CRPC samples exhibited decreased αSMA expression but increased FAP expression in comparison to HSPC samples ( Figure 5B ).

Previous studies indicate that different CAFs subtypes can be associated with different clinical prognosis, we therefore evaluated the prognostic significance of heterogeneous CAFs subtypes in prostate cancer. The expression levels of αSMA and CAV1 were not correlated with patient prognosis in TCGA-PRAD cohort ( Figure S8A ). High CAV1 expression was correlated with favorable prognosis in GSE116918 cohort ( Figure S8B ). In contrast, high expression of FN1 and FAP was associated with poor prognosis in both TCGA-PRAD and GSE116918 cohort ( Figures 6A, B ). Our data consistently showed that expression of αSMA had no association with patient prognosis while high FN1 expression was associated with poor prognosis in Renji cohort ( Figures 6C , S8C ). To better evaluated the prognostic significance of heterogeneous CAFs subtypes, we classified prostate cancer patients into different subtypes by hierarchical clustering of normalized expression profiles of CAFs-C0 and CAFs-C1 signature genes. Survival analyses suggested that patients with higher CAFs-C1 signature had worse prognosis, whereas CAFs-C0 signature had no association with patient prognosis in TCGA-PRAD and GSE116918 cohort ( Figures 6D , S8D ). Of note, pathway analysis on the DEGs of high CAFs-C1 patients from TCGA-PRAD cohort revealed significant enrichments of pathways related to extracellular matrix remodeling ( Figure 6E ).

Since we identified that CAFs-C1 signature enrichment was an adverse prognostic factor in prostate cancer, we sought to investigate the clinical significance of CAFs-C1 signature in distinct cancer types. As shown in Figure 7A , high CAFs-C1 signature was associated with poor overall survival in bladder cancer (TCGA-BLCA), kidney renal clear cell carcinoma (TCGA-KIRC), stomach adenocarcinoma (TCGA-STAD), brain lower grade glioma (TCGA-LGG), mesothelioma (TCGA-MESO), and skin cutaneous melanoma (TCGA-SKCM). Specific subpopulation of CAFs could suppress antitumor immunity by secreting inflammatory cytokines and recruiting immunosuppressive cells (14, 52, 53), we thus studied whether CAFs-C1 signature was involved in treatment resistance to immunotherapy. To this end, we acquired two datasets which contained the transcriptional profiles and treatment information of urothelial carcinoma patients receiving anti-PD-L1 immunotherapy and melanoma patients receiving anti-PD1 immunotherapy. GSEA analyses showed that non-responders had significant enrichment of CAFs-C1 signature in comparison to responders, indicating that CAFs-C1 signature was involved in treatment resistance to immune checkpoint inhibitors Figures 7B, C .