Different cell types in TNBC were identified by single-cell sequencing

The overall flow chart of our study is illustrated in Fig. 1. In this study, a total of 11,811 cells were derived in the quality control procedures using Seurat, including 8566 cells from primary TNBC samples (TNBC1, TNBC2, TNBC3) and 3245 cells from brain metastasis TNBC sample (TNBC4). All cells in the scRNA-seq datasets were clustered into 24 different cell clusters through PCA, and the resolution was set to 0.8. The high-quality cells were visualized using graph-based dimensionality reduction. As shown in Fig. 2A, different cell clusters were significantly distinguished in space. Then, we annotated all cells and distinguished eight different cell types, including B/plasma cells, endothelial cells, fibroblasts cells, luminal cells, myeloid cells, myoepithelial cells, T/NK cells, and unknown (Fig. 2B). The cell markers utilized were sourced from the Cell Marker database (Fig. 2C). These analyses could sufficiently discern the different cell types, for example, KRT8, KRT18, and KRT19 were markers of luminal cells, and thirteen cell clusters (0, 2, 3, 4, 5, 6, 7, 9, 12, 15, 16, 21, and 22) were identified as luminal cells. As shown in Fig. 2D, we described the proportion of the aforementioned eight cell types in inclusive TNBC samples. In summary, we successfully identified eight principal cell types for further exploration.

Fig. 1

Schematic design of the study. GSE, Gene Expression Omnibus series; TNBC, triple-negative breast cancer; TCGA, The Cancer Genome Atlas; NENF, Neuron derived neurotrophic factor

Fig. 2

Different cell types in TNBC were identified by single-cell sequencing. A tSNE and UMAP scatter plots displayed 24 different cell clusters in TNBC samples. B tSNE and UMAP scatter plots displayed 8 different cell types in TNBC samples. C Dot plots showed expression levels of marker genes used to note 8 different cell types. D The proportion of 8 different cell types in TNBC samples

Malignant cells in TNBC were identified by inferCNV analysis

It is obviously difficult to distinguish between benign and malignant epithelial cells at the single cell level through inherent cell markers. Therefore, we used inferCNV to assist in the identification of malignant cells. We first used the CNVs of myeloid cells as a reference to infer the CNVs of luminal cells. As shown in Fig. 3A-C, there was significant amplification or deletion of CNVs in luminal cells in TNBC1, TNBC2 and TNBC3 samples. Then, we accurately isolated malignant luminal cells through CNV correlation and CNV score. As demonstrated in Fig. 3D, malignant cells, defined by inferCNV, exhibited high heterogeneity in gene expression, which was significantly distinguished from control myeloid cells and normal luminal cells. Luminal cells in TNBC4 derived from metastatic sample were considered malignant tumor cells. Furthermore, 53 normal luminal cells were removed and a total of 7629 confident malignant cells were identified. All malignant luminal cells were performed further dimensionality reduction and clustering with the resolution at 0.2, identifying a total of nine different clusters of malignant cells (Fig. 3E). Taken together, the malignant cells were successfully identified for further exploration.

Fig. 3

Malignant cells in TNBC were identified by inferCNV analysis. Chromosomal landscape of inferred CNVs among luminal cells in TNBC1 (A), TNBC2 (B), and TNBC3 (C). D Scatter plot showed the CNV correlation and CNV score of TNBC luminal cells in TNBC1, TNBC2, and TNBC3. E tSNE and UMAP scatter plots displayed 9 different clusters in TNBC luminal cells

Cell trajectory and characteristics of various clusters of TNBC tumor cells

The tendency for metastasis is a characteristic of TNBC, the gene expression patterns during tumor metastasis exhibit temporal heterogeneity. Therefore, thoroughly explore the evolutionary trajectory of TNBC from the primary cells to the metastatic cells will contribute to understand potential biological processes preferably. Trajectory analysis performed for malignant cells uncovered the three states and three branches in the cell trajectory (Fig. 4A-B). Malignant cells of state 3 were presented at the beginning of the trajectory. As shown in Fig. 4C, cell clusters including cluster 2, 3, and 5 evolved to cluster 0, 1, 4, 6, 7, and 8. To further investigate the heterogeneity of malignant cells, we scored each cluster by using Hallmark gene sets. Accompanied by evolutionary trajectories, the expression of MARCKSL1, RBP7, and STMN1 gradually increased, IFI6 and RPL11 displayed as continuous expression. Conversely, the expression of ERRFI1, the crucial negative regulator of EGFR, was significantly decreased during (Fig. 4D). As shown in Fig. 4E, the canonical tumor malignant phenotype related signaling pathways, such as EMT, E2F targets, NOTCH signaling, PI3K/AKT/mTOR signaling were significantly enriched in cluster 0, 4, and 7. The analysis results indicated that the three clusters mentioned above possibly exhibit stronger levels of malignancy. In addition, to identify clusters with poor prognosis for TNBC, we used CIBERSOFTX to infer the abundance of TNBC patients in various clusters in the METABRIC dataset, and combined the cumulative overall survival (OS) information to perform survival analysis on different clusters. Interestingly, the results showed that only the cluster 0, and 4 showed significant correlation with poor prognosis in TNBC patients (Fig. 4F-G), while cluster 7 presented a significant association with favorable prognosis (Supplementary Figure S1A). The prognostic values of the other clusters were illustrated in Supplementary Figure S1B-G. In brief, cluster 0, and 4 displayed distinct malignant features of tumors and were significantly associated with poor prognosis in TNBC patients. These findings revealed a high degree of heterogeneity among different clusters of TNBC and the dynamic evolution process from primary TNBC cells to metastatic TNBC cells.

Fig. 4

Cell trajectory and characteristics of various clusters of TNBC tumor cells. The trajectory of primary TNBC cells evolved into metastasis TNBC cells were revealed by monocle analysis, visualized by pseudo-time (A), distribution of three cell states (B), and distribution of 9 clusters (C). D 6 most relevant genes were identified in the evolutionary process. E The hallmark pathway enrichment score of different tumor subpopulation cells were illustrated by heat map. The OS of TNBC patients with different abundant levels of C0 (F) and C4 (G) cluster, depicted by KM curves

Identification of gene co-expression modules among TNBC cells

As cluster 0 and 4 were closely correlated with poor prognosis in TNBC, it was necessary to explore the co expressed gene networks that exert important roles in these two subgroups. The scale-free network of cluster 0, and 4 were constructed for the best connectivity with soft threshold set at 14 (Fig. 5A-B). Finally, seven modules were identified with representative top 10 genes (Fig. 5C&E). As shown in Fig. 5D, module 6 and 7 were observed a certain extent of correlation, while other modules constituted another relevant group. On the other hand, the enrichment score of modules 6 was more concentrated in cluster 0 and 4 than that of module 7, and was nearly not enriched in other clusters. (Figures 3E and 5E). Consistently, harmonization module characteristic genes (hME) of module 6 was significantly increased in cluster 0 and 4 (Fig. 5F). In a word, these findings implied that genes of module 6 were co-expressed network of cluster 0 and 4, which could promote TNBC metastasis.

Fig. 5

Identification of gene co-expression modules among TNBC cells. A Weighed gene co-expression network analysis was constructed among malignant cells. B The hdWGCNA dendrogram of 7 modules. C The first 10 eigengenes of each module, ranked by eigengene-based connectivity. D The correlation among 7 modules illustrated by heat map. E UMAP scatter plots displayed the expression of module 1–7 among all malignant cells. F The module 6 score in 9 different clusters

NENF was identified metastasis related gene and upregulated in TNBC

To explore the essential genes driving tumor metastasis, we conducted differential analysis of expression profiles between metastatic and primary tumor cells, and identified 44 genes upregulated in metastatic tumor cells (Supplementary File Table S3). Then, we took the intersection of this gene set and module 6 gene set (Supplementary File Table S4). Totally, 16 genes were identified as candidate genes, including APOE, ATP6VOE2, C1QBP, CITED4, EFEMP1, FABP7, IGFBP2, KRT10, MIA, NENF, PCSK1N, SERP1, SMS, SPP1, UQCRFS1 and ZG16B (Fig. 6A). To investigate transcriptional levels of 16 candidate genes, we detected the copy-number variation of these candidate genes in the TCGA-TNBC dataset. In TCGA-TNBC specimens, NENF was the most significantly amplified candidate gene (88/114, 77.2%) (Fig. 6B). Consistent with Fig. 2E, NENF was mainly expressed in clusters 0, 4, and 7 (Fig. 6C). And NENF was highly expressed in metastatic TNBC at both single-cell and tissue level (supplementary Figure S2A-B). As predicted by METABRIC dataset, higher NENF expression level was associated with poorer prognosis of TNBC, both OS (Fig. 6D) and recurrence-free survival (RFS) (Fig. 6E). Moreover, NENF had observed higher expression in TNBC tissues and associated with advanced stages (Fig. 6F-G). RT-qPCR was used to detect NENF mRNA levels in 30 TNBC tissues and paired adjacent normal tissues, the analysis results show that the expression of NENF was remarkably increased in TNBC specimens (Fig. 6H). In addition, 20 TNBC tissues and paired adjacent normal tissues were collected to examine NENF expression by IHC staining. Comparing to the normal tissues, the expression of NENF in TNBC tissues was significantly up-regulated (Fig. 6I). Furthermore, the standard stain score of IHC is established as Fig. 6J. Taken together, these findings showed that high expression of NENF was associated with TNBC metastasis and poor patient prognosis, and it was upregulated in tumor tissues.

Fig. 6

NENF was identified metastasis related gene and upregulated in TNBC. A The intersection of genes in module 6 and upregulated genes in metastatic TNBC cells showed by Venn diagram. B The copy number changes of candidate genes in TCGA-TNBC were displayed by heat map. C tSNE scatter plots displayed the expression of NENF in TNBC malignant cells. The survival analysis of TNBC patients with high‑ or low‑NENF expression levels, including OS (D) and RFS (E), depicted by KM curves. F The NENF expression level in normal tissues and TNBC tissues. G The NENF expression level among TNBC stages. H The mRNA expression level of NENF in 30 paired TNBC tissues detected by RT-qPCR. I The protein expression level of NENF in 20 paired TNBC were detected by IHC and the scores of IHC staining. J The standard score of IHC staining. **p < 0.01, ***p < 0.001

NENF was required for cell invasion and migration through regulating EMT in TNBC

To investigate the role of NENF in TNBC progression, RT-qPCR and western blot were conducted to analyze the mRNA and protein levels of NENF in breast cancer cells line (MCF-7, T47D, SUM159, CAL-51 and MDA-MB-231). The results showed a significant discrepancy in the expression of NENF between hormone receptor positive (HR+) and triple-negative breast cancer cell lines. Obviously, the expression of NENF was significantly higher in triple-negative breast cancer subtype (SUM159, CAL-51 and MDA-MB-231), compared with HR + breast cancer subtypes (MCF-7 and T47D) (Fig. 7A-B). Then, we examine the effect of NENF deficiency on the breast cancer cell invasion and migration by using three specific siRNAs targeting NENF, and found two of siRNAs could efficiently reduce the mRNA and protein expression of NENF in CAL-51 and MDA-MB-231 cells (Fig. 7C-D). By Matrigel coated transwell assay, we found that depletion of NENF resulted in a substantial decrease in the rate of cell invasion (Fig. 7E). Meanwhile, depletion of NENF could reduce cell migration in the Matrigel non-coated transwell (Fig. 7F) and cell wound healing assays (Fig. 7G). Epithelial to mesenchymal transition (EMT) process is crucial for tumor metastasis, which serves as a driving factor for tumor cell invasion and migration [22, 23]. To investigated the effect of NENF expression in breast cancer EMT, western blot was used to detect the expression of epithelial and mesenchymal markers. Compared to control cells, the expression of E-cadherin (epithelial marker) was dramatically elevated, while decreasing the N-cadherin and Vimentin expression (mesenchymal marker) in NENF-depleted CAL-51 and MDA-MB-231 cells, demonstrating that NENF was positively correlated with EMT (Fig. 7H). In brief, these results indicated that depletion of NENF reduced cell invasion and migration through regulating EMT in TNBC.

Fig. 7

NENF was required for cell invasion and migration through regulating EMT in TNBC. The mRNA and protein expression levels of NENF in breast cancer cell lines detected by RT‑qPCR (A) and western blot (B). RT‑qPCR (C) and western blot (D) analysis of NENF expression levels in MDA-MB-231-siNENF cells and CAL-51-siNENF cells compared with siControl cells, respectively. E Cell invasion in cells as in D were detected by Matrigel coated transwell analysis, respectively. F Cell migration in cells as in D were detected by Matrigel non-coated transwell (F) and wound healing (G) analysis, respectively. H The protein expression levels of EMT‑related markers in cells as in D were detected by western blot, respectively. ***p < 0.001

Individualized therapy for TNBC patients based on NENF

Presently, the treatment of TNBC is still limited. Therefore, it is urgent to develop novel and effective therapeutic targets for TNBC. In this study, we normalized the gene matrix of METABRIC dataset for predicting the immune response to immune checkpoints on TIDE database (http://tide.dfci.harvard.edu/). As shown in Fig. 8A-B, low TIDE scores indicate strong response to immune checkpoint (ICI) therapy, and patients with NENF high expression exhibited weak response to ICI therapy. These results indicated that NENF probably be a promising predictor for ICI treatment. In addition, we sifted potential therapeutic drugs targeting high expression of NENF by online database (https://clue.io/, Supplementary File Table S5). Figure 8C revealed 46 molecular pathways targeted by 42 compounds in high NENF group. According to the most important mechanism of action for the high NENF group, including protein synthesis inhibitor, HDAC inhibitor, and ATPase inhibitor. These inhibitors had been conformed to effective for cancers [24,25,26]. In conclusion, our findings provide novel strategy for TNBC immunotherapy and individual treatment.

Fig. 8

Individualized therapy for TNBC patients based on NENF A The correlation between TIDE score and immune response status. B The correlation between different NENF group and immune response status. C The potential drugs identified for high-NENF patients