Spatially and functionally distinct subclasses of breast cancer-associated fibroblasts revealed by single cell RNA sequencing

Abstract

Cancer-associated fibroblasts (CAFs) are a major constituent of the tumor microenvironment, although their origin and roles in shaping disease initiation, progression and treatment response remain unclear due to significant heterogeneity. Here, following a negative selection strategy combined with single-cell RNA sequencing of 768 transcriptomes of mesenchymal cells from a genetically engineered mouse model of breast cancer, we define three distinct subpopulations of CAFs. Validation at the transcriptional and protein level in several experimental models of cancer and human tumors reveal spatial separation of the CAF subclasses attributable to different origins, including the peri-vascular niche, the mammary fat pad and the transformed epithelium. Gene profiles for each CAF subtype correlate to distinctive functional programs and hold independent prognostic capability in clinical cohorts by association to metastatic disease. In conclusion, the improved resolution of the widely defined CAF population opens the possibility for biomarker-driven development of drugs for precision targeting of CAFs.

Fig. 1

Unbiased clustering of fibroblast single cell transcriptomic data reveals four populations. a Schematic representation of negative selection strategy removing CD31⁺, CD45⁺, NG2⁺, and EPCAM⁺ cells. b Gating strategy and quantification of flow cytometry for single cell sequencing. After gating out doublets and DAPI⁺ dead cells, EpCAM⁻CD31⁻CD45⁻NG2⁻ CAFs made up 2.5% of the cells. FSC forward scatter, SSC side scatter. For single marker staining see also Supplementary Figure 1. c Violin plot of detected genes in 784 sorted fibroblasts. d t-SNE layout of CAFs (n = 716) by RPKM-normalized transcriptomic data. Colors represent clusters assigned by DBSCAN. e Expression plots on t-SNE layout. log₂(RPKM + 1) levels of CAF marker genes in individual cells. f Cell size and granularity as determined by forward-scattered light (FSC) and side-scattered light (SSC) of different CAF populations

Fig. 2

Distinct gene profiles in CAF subpopulations account for functional differences. a Heatmap of top 18 differentially expressed genes in each subpopulation estimated by ROTS. b Enrichment of the 150 most SDE in gene ontology (GO) terms. Gene ratio is determined by the number of detected genes within a GO term compared to the total number of genes. c Expression plots based on t-SNE layout. log₂(RPKM + 1) levels of endothelial cell genes in individual cells

Fig. 3

vCAF and mCAF marker genes can be used to trace back subpopulations in tissue sections. a Violin plots of selected vCAF differentially expressed genes in log₂(RPKM + 1). Violin colors represent mean expression of each population. Genes were sorted based on gene ontology terms. b Immunohistochemistry (IHC) staining of desmin on MMTV-PyMT tumor sections (6 µm). Images were acquired from the leading edge and the tumor center. Yellow boxes (left) indicate 2× magnified area (right). c IF staining of Nidogen-2 (green) and CD31 (magenta) or PDGFRα (red) on MMTV-PyMT tumor sections (5 µm) from mice of age 8 weeks, 12 weeks, and 15 weeks (top to bottom). Nuclei were counterstained with DAPI. d Immunofluorescence (IF) staining of Nidogen-2 on human tumor tissue (5 µm). Nuclei were counterstained with DAPI. e Violin plots of selected mCAF differentially expressed genes in log₂(RPKM + 1). f IHC staining of fibulin-1 and PDGFRα in MMTV-PyMT tumor sections (6 µm). Images were acquired from the leading edge and the tumor center. IHC staining of fibulin-1 (g) and PDGFRα (h) in human tumor tissue sections (6 µm). FACS-sort of MMTV-PyMT tumor (i) and mammary gland (j) tissue. Gating on single, living CD45⁻CD31⁻NG2⁻EPCAM⁻ cells followed by gating on PDGFRα⁺ cells (blue box) or PDGFRα^- cells (red box). k Violin plots of cell cycle gene expression in log₂(RPKM + 1). l IF staining for Nidogen-2 (red) and Ki-67 (grey) on sections (5 µm) from PDGFRα-EGFP (green) reporter mice. Arrows indicate Nid2⁺Ki67⁺ cCAF. Scale bar: 50 µm

Fig. 4

dCAF express the oncogenic driver gene and share gene expression with the tumor epithelium. a Violin plots of selected dCAF differentially expressed genes in log₂(RPKM + 1). Violin colors represent mean expression of each population. Genes were sorted based on gene ontology terms. b IF staining of SCRG1 (green) and EPCAM (magenta) on MMTV-PyMT tumor sections (5 µm) from mice of age 8 weeks, 12 weeks, and 15 weeks (left to right). Nuclei were counter stained with DAPI (blue). Yellow boxes indicate the area magnified in the lower panel. c Expression plots based on t-SNE clustering. log₂(RPKM + 1) levels of the virus-derived oncogenic driver gene in individual cells. d IF staining of Nidogen-2 (cyan) and SCRG1 (red) on MMTV-PyMT tumor tissue derived from PFGFRα-EGFP (green) reporter mice. Nuclei are stained with DAPI (blue). White arrowheads and stars indicate mCAF and dCAF, respectively. e IF staining of Nidogen-2 (cyan) and SCRG1 (red) on MMTV-PyMT tumor tissue derived from PFGFRα-EGFP (green) reporter mice. The image was acquired by 2-photon microscopy. Collagen fibers were detected by second harmonic generation (SHG, magenta). Dotted lines separate malignant tissue (T) from stroma (S). Scale bar: 50 µm

Fig. 5

CAF gene profiles correlate in human bulk RNA-sequencing data. Pearson correlation of genes from vCAF and mCAF profiles in TCGA datasets of a breast cancer, b pancreatic ductal adenocarcinoma, c lung adenocarcinoma, and d renal clear cell carcinoma. Correlation of the vCAF profile to e an endothelial metagene and f a matrix metagene. Correlation of the mCAF profile to g an endothelial metagene and h a matrix metagene in TCGA breast cancer data. i Correlation of vCAF and mCAF profiles to functional metagenes in the nested case–control study dataset of breast cancer patients. j Quantification of transwell invasion assay. The average percentage of area covered by invaded cells on the bottom of the membrane was quantified from 4 representative images in n = 3 biological repeats; Data depicted as mean ± s.d. *P = 0.026, **P = 0.0045; two-sided, unpaired Students t-test