The matched primary and adjacent normal samples were collected from seven BC patients during surgery and immediately stored in the GEXSCOPE Tissue Preservation Solution at 2-8°C. Histological characterization was determined according to the criteria of the World Health Organization by pathologists from the hospital. Pathologic staging was performed according to the current International Union against Cancer tumor–lymph node metastasis classification. Experiments were reviewed and approved by the Institutional Review Board of The Second Hospital of Shenzhen, China, and were conducted in compliance with the Helsinki Declaration. Each patient provided written informed consent before sample collection.

Single-cell suspensions were at once processed for the scRNA-seq using the Chromium platform (10x Genomics). Single-cell capture, barcoding, and library preparation were performed by following the 10x Genomics Single Cell Chromium 3′ protocols, and the final libraries were sequenced on the Illumina HiSeq 2500 platform. Data was processed using the CASAVA pipeline (Illumina), and sequencing reads were converted to FASTQ files and UMI read counts using the CellRanger software (10x Genomics, v.2.1.1).

The single-cell sequencing reads from the 10x Genomics Chromium were demultiplexed and then aligned to the GRCh38.p12 human genome reference using the CellRanger pipeline (v.3.1.0, 10x Genomics) with the default parameters and UMI count matrices against the corresponding Ensembl gene annotations were generated. The reference genome and the gene annotation were downloaded from the UCSC genome browser (27).

Count matrices were further analyzed using Scanpy (https://github.com/scverse/scanpy) (version 1.9.3) (28). Of the fourteen samples, sample 5 and sample 9 were removed because limited cells were detected and sample 4 was removed because the expressed genes were abnormally lower than other samples. The count matrices of the remaining eleven samples were merged. Cells with fewer than 200 genes expressed were filtered and genes detected in less than 3 cells were filtered. In addition, cells with higher than 20% mitochondrial or 50% ribosomal RNAs were also filtered as they represented low-quality cells. Further, for potential doublets, we did not include cells expressing more than 5000 genes and used scDblFinder (29) to further mark the remaining cells. We also calculated the cell cycle scores for each cell using sc.tl.score_genes_cell_cycle.

Then the matrix was normalized based on total UMI counts per cell (1e4) and then log-transformed. To reduce potential batch effects, a graph-based data integration algorithm batch-balanced KNN (bbknn) (30) was applied to the matrix. Ridge regression (31) was also performed on samples and Leiden clusters. The integrated data was then used for downstream analysis. Highly variable genes were identified using sc.pp.highly_variable_genes from scanpy with parameters min_mean=0.0125, max_mean=3, min_disp=0.5. Principal components were calculated using sc.pp.pca and sc.pl.pca_variance_ratio was used to determine the number of principal components used for clustering. The neighborhood graph was calculated using sc.pp.neighbors and then clustered using sc.tl.leiden with parameter resolution=0.4. Then the clusters were visualized using Uniform Manifold Approximation and Projection (UMAP) calculated by sc.tl.umap.

To annotate the major cell types, sc.tl.rank_genes_groups was applied to identify differentially expressed genes in each cluster using the Wilcoxon rank-sum test. The top significant differentially expressed genes were then manually reviewed. In addition, we checked each cluster using the known canonical markers such as EPCAM, KRT19, KRT14, ERBB2, ESR1 for epithelial cells, PECAM1, VWF for endothelial cells, DCN, COL1A1, COL1A2, COL3A1, CFD, and PRGFRB for CAFs, ACTA2, TAGLN, MCAM for perivascular cells, CD79A, CD79B for B cells, LYZ, IL1B, MSR1 for macrophages, JCHAIN, MZB1 for plasma cells, and CD3G, CD3D, IL7R, NKG7, GNLY, CD8A for T cells.

To investigate the heterogeneity of CAFs, cells from the subtype were extracted, and similar procedures were performed. In the clustering, we selected a high-resolution parameter (resolution_value=0.8) to obtain the fine-tuned subpopulations from this cell type. Differentially expressed genes were calculated for each CAF subpopulation against the remaining CAF subpopulations. The Wilcoxon rank-sum test was used for the statistical significance test as previously.

To infer the potential cell-cell communication network between clusters, CellChat (32) (version 1.6.1) was used to quantitatively measure networks following the tutorial. A CellChat object was created from the normalized expression matrix with the createCellChat function. After preprocessing data with identifyOverExpressedGenes, identifyOverExpressedInteractions, computeCommunProb, filterCommunication, and computeCommunProbPathway were used to calculate potential ligands–receptor interactions. Finally, the aggregated cell–cell communication network was calculated using the “aggregateNet” function.

The STRING database (https://string-db.org/) (33) was used to identify SFRP4-related protein-protein interactions. Gene set enrichment analyses (GSEA) were performed on the significantly differently expressed genes using GSEApy (https://github.com/zqfang/GSEApy) (34) with gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) and MSigDB Hallmark signatures gene sets. Three GO terms biological processes (BPs), cellular components (CCs), and molecular functions (MFs) were included. We used a false discovery rate (FDR) of 0.25 as the significance criterion.

The copy number variation of single cells was inferred using InferCNV (https://github.com/broadinstitute/infercnv). A set of reference ‘normal’ cells were randomly selected from immune cells.

To infer the pseudotime of CAF subclusters, diffusion maps were calculated using sc.tl.diffmap. A putative initial cell was selected after studying the individual diffusion components and identifying the most extreme diffusion component in one dimension. The pseudotime was then calculated using sc.tl.dpt.

Bulk sequencing data was deconvoluted using Scaden (https://github.com/KevinMenden/scaden) (35). Scaden is a deep-learning based algorithm for cell type deconvolution of bulk RNA-seq samples. scRNA-seq data and bulk sequencing data were prepared as its instructions, and data simulations, data preprocessing, training, and prediction were run sequentially.

To visualize the spatial distribution of SFRP4+ CAF, we analyzed BC spatial transcriptomics data (V1_Breast_Cancer_Block_A_Section_1) downloaded from 10x Genomics. Cells with less than 5,000 counts or more than 35,000 counts were filtered. Cells with mitochondrial percentage >20 were also filtered. Genes detected in less than 10 cells were filtered. Then the counts were normalized, and log transformed. Highly variable genes were selected using sc.pp.highly_variable_genes with parameters flavor=“seurat”, n_top_genes=2000. Then sc.pp.neighbors, sc.tl.umap, and sc.tl.leiden were run to cluster the cells. sc.pl.spatial was used to visualize the expressions.

Two datasets [GSE20194 (36), GSE20271 (37)] with neoadjuvant chemotherapy (NAC) response information were downloaded from the NCBI GEO database. GSE20194 had 278 samples and GSE20271 had 178 samples. The median expression of the target gene was used to split samples into high and low-expression groups in each cohort. The average gene expression in samples according to NAC response, pathological complete response (pCR) and residual disease (RD), was also compared.

Python package lifelines (https://lifelines.readthedocs.io/en/latest/, version 0.26.4) (21) was used to perform Kaplan–Meier curve analysis, log-rank test in METABRIC dataset, and p < 0.05 was considered as significant. To validate the result, we further check the prognosis value of genes using an online database Kaplan-Meier Plotter (38).

Public single-cell gene expression datasets GSE164898 and GSE113197 were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). Bulk sequencing data METABRIC was downloaded from cBioPortal. TCGA-BRCA including normal and tumor samples were acquired from the TCGA database (http://portal.gdc.cancer.gov/). The expression of BC cell lines was downloaded from the CCLE database.

To characterize cell expression heterogeneity in BC, we performed droplet-based single-cell transcriptome profiling (scRNA-seq, 10× Genomics Chromium system) on fresh tumors and adjacent non-tumor tissues from seven breast cancer (BC) patients after surgery ( Figure 1A ). The clinical and pathological information of the patients was summarized in Table 1 . An average of 5,752 cells per sample was sequenced with 221,292 average reads per cell. After quality control (Methods), eleven of fourteen samples and 71,773 cells were retained, of which 60.4% were from tumor samples.

The count matrix was normalized, the dimensionality reduction was performed by the principal component analysis, and the graph-based Leiden clustering was applied to classify cells into twelve transcriptionally distinct clusters and the clusters were visualized using UMAP ( Figure 1B ). We annotated the clusters into eight major cell types based on the expression of canonical cell type marker genes and the top differentially expressed genes of each cluster ( Figure 1C ). The sample level distribution of cells was shown in Figure 1D .

We observed heterogeneity in the cell type compositions across samples, even from the same individual ( Figure 1E ). The cell type markers included EPCAM, KRT19, KRT14, ERBB2, ESR1 for epithelial cells ( Figure 1F ), PECAM1, VWF for endothelial cells ( Figure 1G ), DCN, COL1A1, COL1A2, COL3A1, CFD, and PRGFRA for cancer-associated fibroblasts (CAFs) ( Figure 1H ), ACTA2, TAGLN, MCAM for perivascular cells ( Figure 1I ), CD79A, CD79B for B cells ( Figure 1J ), LYZ, IL1B, MSR1 for macrophages ( Figure 1K ), JCHAIN, MZB1 for plasma cells ( Figure 1L ), and CD3G, CD3D, IL7R, NKG7, GNLY, CD8A for T cells ( Figure 1M ). A matrix plot of expressions of these markers across the cell types is shown in Figure 1N . In addition to the CAF markers mentioned above, GSN, LUM, APOD, and FBLN1 were also highly expressed.

To investigate the genetic alteration in CAF, we analyzed the copy number variation (CNV) using InferCNV. No significant CNV was detected in CAFs ( Supplementary Figure 1 ), which was consistent with the previous report (39). Furthermore, we analyzed the cell-cell communications between major cell types using CellChat. We found complex interactions between CAF and other cell types ( Supplementary Figure 2 ). Three outgoing patterns were identified in the communications of secreting cells, and CAF was dominant in pattern 1, which included interactions of COLLAGEN, CXCL, LAMININ, et al. ( Figure 1O ). This is consistent with the functional of fibroblasts in cell matrix.

To explore the cellular heterogeneity within CAFs, we normalized, reduced the dimensionality, and clustered 11,831 CAF cells using a finely tuned pipeline, described in Methods. Eight CAF subtypes have been identified ( Figure 2A ). The expressions of the top five marker genes for each subtype are shown in Figure 2B . CAF subtypes 1, 3, 5, and 6 clustered together, while the remaining subclusters 0, 2, 4, and 7 formed another group. The percentages of each CAF subcluster across samples were shown in Figure 2C . We named each CAF subtype according to their most distinctive marker gene: BTG1+ CAF, corresponding to subtype 0; OGN+ CAF, corresponding to subtype 1; CFD+ CAF, corresponding to subtype 2; C1R+ CAF, corresponding to subtype 3; IGFBP7+ CAF, corresponding to subtype 4; MFAP5+ CAF, corresponding to subtype 5; SFRP4+ CAF, corresponding to subtype 6; and PTMA+ CAF, corresponding to subtype 7. The expression of the top 2 marker genes for each subtype were shown in Figure 2D . We found that SFRP4+ CAF was the most separated from the other subtypes in the UMAP plot ( Figure 2A ). SFRP4 was specifically expressed in this subtype, while CLU was also expressed in other cell types ( Figure 2E ).

We performed differential expression gene analyses between CAFs and other cell types to investigate the potential function of CAFs. The differentially expressed genes involved extracellular matrix organization, collagen fibril organization, transforming growth factor beta binding, and ECM-receptor interaction, according to GO and KEGG enrichment ( Figure 2F and Supplementary Figure 3 ). GSEA enrichment on MSigDB Hallmark signatures showed that they associated with estrogen response late, myogenesis, epithelial-mesenchymal transition, and cholesterol homeostasis pathways ( Figure 2G ). These findings suggest that CAFs may play a key role in modulating the tumor microenvironment and promoting BC progression and metastasis.

We then analyzed the cell-cell communication patterns between CAF subtypes and other cell types. We found that CAF subtypes were enriched in pattern 1 and pattern 4 while T cells enriched in pattern 6, B cell/macrophage enriched in pattern 2, endothelial cells in pattern 3 and epithelial cells in pattern 5 ( Figure 2H ). Pattern 1 involved signaling pathways of COLLAGEN, CD99, MK, CXCL et al. while pattern 4 involved CLDN, CD34, ncWNT, et al. These patterns imply that CAFs may influence the immune response and the epithelial-mesenchymal transition of tumor cells. BTG1+ CAF, OGN+ CAF, CFD+ CAF, IGFBP7+ CAF were in pattern 1, and C1R+ CAF, MFAP5+ CAF were in pattern 4 ( Supplementary Figure 4 ). Interestingly, SFRP4+ CAF had a weak signal in four patterns, suggesting that it may have a different communication mode from the other subtypes. A diffusion pseudotime analysis suggested that SFRP4+ CAF was in a distinct path from the other subtypes ( Figure 2I ). The distinct cell-cell communication patterns and different pseudo development trajectory indicate that SFRP4+ CAF may have a unique biological function and role in breast cancer.

We compared the percentages of each CAF subtype between the normal and tumor samples using scRNA-seq data. We found that CAF subtypes 0 (BTG1+ CAF), 2 (CFD+ CAF), and 4 (IGFBP7+ CAF) were more abundant in the tumor than the normal, while CAF subcluster 1 (OGN+ CAF), 3 (C1R+ CAF), 5 (MFAP5+ CAF), and 6 (SFRP4+ CAF) were more abundant in the normal than the tumor ( Figure 3A ). Most of the marker genes of SFRP4+ CAF (subcluster 6) were highly expressed in the normal samples ( Figure 3B ). Of the top 30 SFRP4+ CAF marker genes, eleven genes (11/30 = 36.7%) were significantly upregulated in the normal samples ( Supplementary Figure 5 , Student’s T test, p<0.001). The expression of SFRP4 in normal samples and tumor samples were shown in Figure 3C .

We then confirmed the expression specificity of SFRP4 in SFRP4+ CAF by reanalyzing several publicly deposited datasets. GSE75688 (40), a BC single-cell dataset generated by full-length single-cell RNA sequencing, also demonstrated that SFRP4 was only expressed in some CAF cells. GSE113197 (41), another BC scRNA-seq dataset generated from cells first sorted using epithelial cell surface markers, confirmed that SFRP4 was not expressed in epithelial cells. GSE164898 (42), another scRNA-seq dataset in the breast, confirmed the existence of SFRP4+ CAF. The eight cell types were clustered and annotated using the same pipeline for the in-house dataset ( Figure 3D left). Marker DCN was highly expressed in CAF cluster ( Figure 3D middle) and SFRP4 was expressed in some CAFs ( Figure 3D right). Further CAFs were fine clustered into five subpopulations and genes SFRP4, CLU, OGN were highly expressed in CAF subpopulation 3 ( Figure 3E ). The marker gene expressions of CAF subpopulations were shown in Figure 3F . We further confirm the existence of SFRP4+ CAF using a BC scRNA-seq data GSE225600. SFRP4 and CLU were highly expressed in the CAF subcluster 5 ( Figure 3G ). Then we accessed the expression distribution of SFRP4 in BC tissue using a spatial transcriptome dataset from Illumina. We found that the distribution of SFRP4 was aligned with DCN, a CAF marker gene, and separated from EPCAM, an epithelial cell marker gene ( Figure 3H ; Supplementary Figure 6 ). In addition, we compared the expression of SFRP4+ CAF marker genes in normal and tumor samples using data from the TCGA-BRCA cohort and found that most of the marker genes of SFRP4+ CAF were significantly down-regulated in the tumor samples ( Supplementary Figure 7 ), which was consistent with our scRNA-seq results.

We investigated the possible function of SFRP4+ CAF by constructing a protein interaction network from STRING-db. We found that SFRP4 could directly interact with several WNT components ( Figure 3I ). GO enrichment of the top marker genes of SFRP4+ CAF showed that they were involved in collagen-containing extracellular matrix and focal adhesion in cellular component items, and regulation of cell migration and negative regulation of cell population proliferation in biological process items ( Figure 3I ). GSEA enrichment of genes found the enrichment of the WNT signaling pathway and frizzled binding. The MSigDB hall markers also found enrichment of the Wnt-beta catenin signaling pathway ( Figure 3I ). In the bulk sequencing data, we stratified the tumor samples into two groups according to the median expression of the marker genes (SFRP4, CLU, OGN). We found that samples with low expression of these genes have significantly higher cell cycle scores, DNA damage response scores, and cell proliferation scores ( Figure 3J ). These findings suggest that SFRP4+ CAF may play a key role in modulating the tumor microenvironment and inhibiting tumor progression.

We explored the interaction of SFRP4+ CAF with WNT signaling by checking the expression level of the WNT genes in our scRNA-seq data and CCLE database. We found that WNT5A was highly expressed in our scRNA-seq dataset using all epithelial cells ( Supplementary Figure 8 ) and BC cell lines ( Supplementary Figure 9 ). Therefore, we hypothesized that SFRP4+ CAF may inhibit BC progression through WNT5A in the WNT pathway.

We hypothesized that SFRP4+ CAFs secrete SFRP4 and inhibit the WNT signal pathway in BC cells. SFRP4 is a secreted frizzled-related protein that interferes with WNT-Fz interactions and thus inhibits the WNT pathway (43, 44). Our analysis suggested that SFRP4+ CAFs may contribute to cell proliferation and migration. To assess our hypothesis, we silenced SFRP4 in cancer-associated fibroblasts (CAFs) using three siRNAs and examined the effects on BC cells. We confirmed the reduced expression of SFRP4 in CAFs ( Figure 4A ). ELISA experiment demonstrated that the level of conditioned media SFRP4 protein in cell-culture dish decreased when siSFRP4 silencing ( Figure 4B ). We then treated SK-BR-3 cells with the conditioned media of CAFs-siSFRP4 and observed an up-regulation of CTNNB1 and GSK3B mRNA, which are involved in WNT signaling ( Figure 4C ). This was consistent with a weak negative correlation between SFRP4 and GSK3B mRNA levels using the GEPIA database ( Figure 4D ). The protein from SK-BR-3 and H1937, treated with conditioned media of CAFs-siSFRP4, were used for detecting the activation of WNT signaling. In addition, SFRP4 silencing in CAFs induced EMT phenotypic transition, as evidenced by increased vimentin and decreased E-cadherin expression. However, we did not observe any changes in WNT5a. The level of GSK3β phosphorylation and β-catenin accumulation were significantly up-regulated in SK-BR-3 and H1937 BC cells after treatment with CAFs-siSFRP4 conditioned media. Moreover, we detected an activation of AKT and ERK1/2 pathways ( Figures 4E, F ). The SK-BR-3 and HCC1937 cells had a higher cell proliferation rate when co-culture with siSFRP4-CAFs compared with the control group ( Supplementary Figure 10 ). These results indicate that SFRP4 silencing in CAFs activates WNT signaling pathway. Furthermore, the healing and migration assay demonstrated that the capability of invasiveness and migration of BC cells were significantly increased when treatment with the conditioned media of CAFs-siSFRP4 ( Figures 4G–J ). It is known that WNT-β-catenin activation enhanced the transcriptional function of TCF, so we also measured the mRNA expression of several downstream targets of TCF, such as MMP-7, CD44, MYC, COX2, FN, and SLUG, and found that they were generally increased after SFRP4 silencing in CAFs ( Figure 4K ). In conclusion, our study showed that SFRP4+ CAFs inhibit the WNT-β-catenin pathway in BC cells by secreting SFRP4. Silencing SFRP4 in CAFs induced EMT and increased the metastatic potential of BC cells.

Since SFRP4+ CAF can modulate the WNT pathway, it may influence the prognosis of BC. We show the associations by performing Kaplan–Meier curve survival analysis and log-rank test in the METABRIC dataset. From the top fifteen marker genes, nine of them - SFRP4, OGN, ADIRF, MGP, COL14A1, CST3, PLAC9, SERPINF1 and ABI3BP - consistently predict better overall survival (OS) ( Figure 5A ) and relapse-free survival (RFS) ( Figure 5B ) while higher expressed. We then validate the prognosis of these genes in an online database Kaplan-Meier Plotter and find that except for ADIRF did not exist, and PLAC9 is not significant in the Kaplan-Meier Plotter, others are consistently significant in the Kaplan-Meier Plotter ( Supplementary Figure 11 ). Further, we deconvolute the METABRIC dataset using Scaden. We show that SFRP4+ CAF composition can also predict the OS ( Figure 5C ) and RFS ( Figure 5D ), the higher the percentage of SFRP4+ CAF, the better the OS and RFS.

The modified tumor microenvironment may also contribute to drug therapy. Therefore, we further check the gene expression of SFRP4 and the BC NAC outcomes. To our surprise, the pCR group has significantly lower SFRP4 expression than the RD group while the low SFRP4 group has a significantly higher pCR rate than the high SFRP4 group in TNBC subtypes and in all BC samples ( Figure 5E ). We further validate the result using an independent dataset ( Figure 5F ). This result suggests that SFRP+ CAF may have an impact on the NAC outcomes.