Single-cell landscape of intratumoral heterogeneity and tumor microenvironment remolding in pre-nodal metastases of breast cancer

Abstract

Background: The metastasis of cancer cells is influenced by both their intrinsic characteristics and the tumor microenvironment (TME). However, the molecular mechanisms underlying pre-nodal metastases of breast cancer remain unclear.

Methods: We integrated a total of 216,963 cells from 54 samples across 6 single-cell datasets to profile the cellular landscape differences between primary tumors and pre-nodal metastases.

Results: We revealed three distinct metastatic epithelial cell subtypes (Epi1, Epi2 and Epi3), which exhibited different metastatic mechanisms. Specifically, the marker gene KCNK15 of the Epi1 subtype exhibited increased gene expression along the cell differentiation trajectory and was specifically regulated by the transcription factor ASCL1. In the Epi3 subtype, we highlighted NR2F1 as a regulator targeting the marker gene MUCL1. Additionally, we found that the Epi2 and Epi3 subtypes shared some regulons, such as ZEB1 and NR2C1. Similarly, we identified specific subtypes of stromal and immune cells in the TME, and discovered that vascular cancer-associated fibroblasts might promote capillary formation through CXCL9⁺ macrophages in pre-nodal metastases. All three subtypes of metastatic epithelial cells were associated with poor prognosis.

Conclusions: In summary, this study dissects the intratumoral heterogeneity and remodeling of the TME in pre-nodal metastases of breast cancer, providing novel insights into the mechanisms underlying breast cancer metastasis.

Keywords: Breast cancer; Pre-nodal metastasis; Single-cell analysis; Tumor microenvironment remodeling.

Fig. 1

Single-cell transcriptomic analysis of primary tumor and pre-nodal metastasis in breast cancer. a Workflow depicting the data collection, processing, and analysis in this study. b t-SNE plot displaying the distribution of 216,963 cells, categorized by primary tumor and pre-nodal metastasis. c Bubble plot showing the average expression of canonical marker genes for seven cell types. The size of each bubble represents the fraction of cell types expressing the gene. The color intensity of bubbles represents the scaled expression level of genes. d Radar plot showing the relative cell proportions of seven cell types in each group. The numbers surrounding the panel indicate the maximum range of cell proportions, starting from zero

Fig. 2

Reclustering analysis of epithelial cells. a t-SNE plot displaying the distribution of 67,784 epithelial cells, colored by clusters. b t-SNE plots displaying the distribution of Scissor-selected epithelial cells, categorized by primary tumor and pre-nodal metastasis. The red and blue dots represent cells associated with the metastatic and non-metastatic phenotypes, respectively. c t-SNE plot displaying the distribution of 67,784 epithelial cells, colored by subtypes. d Heatmap showing the expression levels of representative DEGs across epithelial cell subtypes. e Barplots showing the GO enrichment of DEGs in Epi1-3 subtypes

Fig. 3

Single-cell copy number variation analysis and trajectory analysis of epithelial cell subtypes. a The hierarchical heatmap showing the CNV profiles of epithelial cell subtypes. b Violin plots showing the CNV scores of epithelial cell subtypes. c Clonality trees showing the developmental course of epithelial cells, categorized by primary tumor and pre-nodal metastasis. The length of each branch is scaled according to percentage of cells in the subclone. d–g Monocle plots showing the differentiation trajectories of epithelial cells, colored by states (d), epithelial cell subtypes (e), EMT scores (f), and CtyoTRACE scores (g). h Heatmap showing the expression patterns of 37 branch-dependent genes across 3 states, which were identified by branched expression analysis modeling (BEAM) using the top 5 DEGs of each epithelial cell subtype. i, j Two-dimensional plots showing the expression patterns of two representative genes in cells of cell fate 1 (blue) and cell fate 2 (pink), respectively, along the pseudotime

Fig. 4

Transcriptional regulatory analysis of metastatic epithelial cell subtypes. a Heatmap showing the top 5 enrichment of transcription factors (TFs) in each epithelial cell subtype. b Network showing the TFs and their predicted target genes. The red nodes represent the differentially expressed genes (DEGs) in Epi1-3. The colors of the regions correspond to the transcriptional regulatory network specific to the Epi1 (red), Epi2 (yellow), and Epi3 (pink) subtypes. c Heatmap showing TF modules identified by Pearson correlation based on TF activities in epithelial cell subtypes. d t-SNE plots displaying the distribution of epithelial cell subtypes and lymph node epithelial cells. LN, lymph node. e Heatmap showing the correlation of TFs between epithelial cell subtypes and lymph node epithelial cells. Pearson correlation analysis was applied to calculate the correlation coefficients. f, g Violin plots showing the Epi scores of Epi1-3 subtypes in two breast cancer bone metastasis samples (f) and one breast cancer brain metastasis sample (g)

Fig. 5

Remodeling of the pre-nodal metastasis microenvironment. a t-SNE plots displaying the distribution of T cell subtypes between primary tumor and pre-nodal metastasis. b Heatmap showing the average expression of canonical marker genes for T cell subtypes. Barplot above the heatmap indicates the relative proportions of T cell clusters in primary tumor and pre-nodal metastasis. c Scatterplots showing the relative proportions of four T cell subtypes between primary tumor and pre-nodal metastasis. The Wilcoxon rank-sum test was applied to determine statistical significance. Each point represents a sample. d Bubble plot showing the expression pattern of signature genes for four T cell subtypes between primary tumor and pre-nodal metastasis. The size of each bubble represents the fraction of T cell subtypes expressing the gene. The color intensity of bubbles represents the scaled expression level of genes. e t-SNE plots displaying the distribution of fibroblast subtypes between primary tumor and pre-nodal metastasis. f Heatmap showing the average expression of canonical marker genes for fibroblast subtypes. Barplot above the heatmap indicates the relative proportions of fibroblast clusters in primary tumor and pre-nodal metastasis. g Line plots showing the differences in cell proportions of specific cell subtypes between primary tumor and pre-nodal metastasis. Each point represents a specific cell subtypes. The top circular plots depict the relative cell proportions

Fig. 6

Cellular interactions in the breast cancer microenvironment. a, b Heatmaps showing the overall interaction strength between specific cell subtypes in primary tumor (a) and pre-nodal metastasis (b). c Bubble plots showing the differences of specific ligand-receptor interactions from fibroblast subtypes to CXCL9⁺ macrophages between primary tumor and pre-nodal metastasis. d Bubble plots showing the differences of specific ligand-receptor interactions from CXCL9⁺ macrophages and FOLR2⁺ macrophages to endothelial cell subtypes between primary tumor and pre-nodal metastasis

Fig. 7

The clinical implications of metastatic epithelial cell subtypes. a-d Kaplan–Meier curves showing the survival differences between the high-score and low-score groups classified by metastatic epithelial cell subtypes in TCGA (a), GSE20685 (b), GSE7390 (c), and GSE3143 (d). e, f Barplots showing the differences of gene mutations between the high-score and low-score groups classified by metastatic epithelial cell subtypes in TCGA (e) and METABRIC (f). g Barplots showing the differences of clinical characteristics between the high-score and low-score groups classified by metastatic epithelial cell subtypes. h, i Ternary plots showing the relative Epi scores of Epi1-3 subtypes before and after treatment for patients who only received anti-PD1 treatment (h) and patients who received neoadjuvant chemotherapy followed by anti-PD1 treatment (i)