Single-cell and bulk transcriptome analyses unravel cellular heterogeneity in TNBC
This study reconstructed a single-cell landscape of TNBC based upon scRNA-seq data from four primary TNBC specimens. Firstly, single cells with empty droplets or low quality were of removal, with 2599 / 3267 cells in GSM5457199 sample, 3872 / 4161 cells in GSM5457205 sample, 3755 / 4064 cells in GSM5457208 sample, 6233 / 7521 cells in GSM5457211 sample being retained in our analysis (Supplementary Fig. 1A-L). Next, the retained scRNA-seq data were scaled based upon PCA, with PC = 9 (Supplementary Fig. 2A-D). By use of UMAP approach, the selected single cells were clustered into 14 clusters, with remarkable cellular heterogeneity between TNBC and non-TNBC (Supplementary Fig. 3A-C). Also, marker genes in each cell cluster were determined (Supplementary Fig. 3D, E). In combination with the known marker genes of cell types, nine cell populations were classified, composed of B cells (n = 698), dendritic cells (n = 387), endothelial cells (n = 495), epithelial cells (n = 9984), fibroblasts (n = 350), macrophages (n = 732), monocytes (n = 32), plasmablasts (n = 1412), and T cells (n = 2369) (Fig. 1A). The marker genes were specifically expressed in the corresponding cell populations: MS4A1 for B cells, CD1C, and FCER1A for dendritic cells, PECAM1, VWF, CDH5, SELE, and CD34 for endothelial cells, EPCAM, CDH1, and KRT18 for epithelial cells, COL1A1 and PDGFRB for fibroblasts, APOE, CD68, MRC1, MSR1, and CXCL2 for macrophages, FCN1, LILRA5, and S100A8 for monocytes, JCHAIN for plasmablasts, and CD3D, CD3E, CD3G, and CD2 for T cells (Fig. 1B). These cell populations were notably different between TNBC and non-TNBC, with the higher cell ratio of B cells, dendritic cells, fibroblasts, macrophages, plasmablasts, and T cells, and the lower cell ratio of endothelial cells and epithelial cells in TNBC versus non-TNBC (Fig. 1C). We also determined novel marker genes of each cell population (Fig. 1D). We also gathered bulk transcriptome profiling of TNBC specimens from the TCGA dataset. By executing CIBERSORTx, the reference matrix of cell markers was established based upon scRNA-seq results, and the relative cell ratios of cell populations were estimated in bulk tissues (Fig. 1E). Consistent with scRNA-seq findings, the higher cell ratio of macrophages was found in bulk TNBC than non-TNBC tissues (Fig. 1F, G). Thus, macrophages were active in the TNBC microenvironment.
Cell-cell interactions in the TNBC and non-TNBC microenvironment
Next, we evaluated cell-cell interactions based upon ligand-receptor pairs in non-TNBC and TNBC, respectively. In comparison to non-TNBC, more active cell-cell interactions were investigated in TNBC, especially macrophages with other cell populations (Fig. 2A, B).
Transcriptional activity of macrophages in TNBC
Totally, 104 transcription factors presented differential expression in TNBC macrophages in comparison to non-TNBC macrophages (Fig. 2C; Supplementary table 1), which might transcriptionally regulate macrophage activity in TNBC. In addition, 110 DEGs were determined between TNBC and non-TNBC macrophages (Fig. 2D, E), of which 63 were up-regulated (Supplementary table 2) and 47 were down-regulated (Supplementary table 3) in TNBC macrophages versus non-TNBC macrophages. The DEGs were notably in relation to signaling transduction, cell surface receptor signaling pathway, response to cytokine, cellular response to cytokine stimulus, etc. (Fig. 2F). In addition, immune-related pathways were notably enriched by the DEGs, especially antigen processing and presentation (Fig. 2G, H). Altogether, the functional enrichment analysis unraveled the significance of the DEGs in modulating transcriptional activity of macrophages in the microenvironmental niche.
Construction of a novel macrophage-relevant prognostic signature for TNBC
The differentially expressed transcription factors and DEGs related to TNBC macrophages were included for univariate-cox regression analysis. Consequently, eight genes exhibited significant associations with TNBC prognosis (p ≤ 0.05), composed of C12orf60, CTSD, CTSL, ELK4, FCGR2A, FOLR2, HSPA8, and XRCC4. The genes were utilized for construction of a LASSO model. We randomized TCGA TNBC samples into training or test cohort. In the training cohort, LASSO analysis was executed for selecting signature genes with regression coefficient ≠ 0. Under the lambda minimum = 0.0267. As a result, five signature genes were eventually determined, comprising CTSD, CTSL, ELK4, HSPA8, and XRCC4 (Fig. 3A, B). The macrophage-based prognostic signature was constructed based upon the formula: risk score = 0.859575100676907 * CTSD expression + 0.0210700891980921 * CTSL expression + (-0.644138418956012) * ELK4 expression + 0.307340530719732 * HSPA8 expression + 1.31660312733179 * XRCC4 expression (Fig. 3C, D). Patients with risk score > median risk score were defined as high risk, while those with risk score ≤ median risk score were defined as low risk (Fig. 3E). High-risk group owned more dead or recurred/progressed status versus low-risk group (Fig. 3F, G). Survival analysis demonstrated the significantly shorter OS time for high-risk patients in the training cohort (Fig. 3H). Such survival difference was proven in the test and entire cohorts (Fig. 3I, J). In addition, ROCs were plotted for appraising the predictive efficacy of the model. In the training (Fig. 3K), test (Fig. 3L), and entire (Fig. 3M) cohorts, the model was excellently predictive of patients’ one-year survival (AUC > 0.9).
Associations of the macrophage-based prognostic model with more advanced status
Further analysis was conducted for evaluation of the correlations between the macrophage-based prognostic model and TNBC clinicopathological traits. TCGA patients with more advanced T, N, M stages, and pathological stages presented the notably higher risk score (Fig. 4A-D). In addition, the risk score displayed a negative association with tumor purity (Fig. 4E). Both in the GSE96058 and GSE45255 cohorts, patients with more advanced histological grades exhibited the prominently higher risk score (Fig. 4F, G). Altogether, the macrophage-based prognostic model was in relation to more advanced status of TNBC patients.
External verification of the macrophage-based prognostic model
The GSE96058 and GSE45255 cohorts were adopted for independently proving the efficacy of the macrophage-relevant signature in patient survival. In the GSE96058 cohort, it was proven that high-risk cases owned worse OS outcomes versus low-risk cases (Fig. 4H). Also, the model enabled to accurately predict one-year survival (Fig. 4I). In the GSE45255 cohort, shorter DFS time was investigated in high-risk cases, with the reliable efficacy in DFS prediction (Fig. 4J, K). These findings proved the generalizability of the macrophage-based prognostic model.
Heterogeneous somatic mutations between low- and high-risk TNBC patients
The study quantified TMB score across TNBC samples, with the median TMB of 1.34/MB (Fig. 5A). Overall, low-risk samples occurred relatively higher TMB score versus high-risk samples (Fig. 5B). TP53 had the highest mutated frequency both in low- and high-risk samples (Fig. 5C). OBSCN, UTP20, and KMT2D had significantly higher mutated frequency in high- versus low-risk cases, while FLG occurred significantly lower mutated frequency in low- versus high-risk cases (Fig. 5D, E). In addition, more frequent co-occurred mutations were investigated in high- than low-risk samples (Fig. 5F, G).
Molecular mechanisms underlying the macrophage-based prognostic model
For biological processes, mitotic spindle assembly checkpoint, and negative regulation of mitotic metaphase anaphase transition and mitotic sister chromatid separation were notably enriched in high-risk samples, while positive regulation of antigen receptor-mediated signaling pathway, and antigen processing and presentation of endogenous antigen, and positive T cell selection were prominently enriched in low-risk samples (Fig. 6A). For cellular components, preribosome, large subunit precursor, and preribosome were notably enriched in high-risk group, with the significant enrichment of immunoglobulin complex and T cell receptor complex in low-risk group (Fig. 6B). For molecular functions, 3’-5’ DNA helicase activity, single-stranded DNA helicase activity, and WNT-activated receptor activity were enriched in high-risk cases, with the prominent enrichment of peptide antigen binding, C-C chemokine binding, and immunoglobulin receptor binding in low-risk cases (Fig. 6C). Moreover, for KEGG pathways, DNA replication, RNA polymerase, and mannose type O-glycan biosynthesis presented the notable enrichment in high-risk group, with the notable enrichment of type I diabetes mellitus, allograft rejection, and graft-versus-host disease in low-risk group (Fig. 6D).
High-risk patients with higher response to immunotherapy
High-risk group exhibited the higher T-cell inflamed score and lower TIDE score in comparison to low-risk group (Fig. 7A, B). In addition, the higher ratios of responders to immunotherapy were investigated in high- versus low-risk group (Fig. 7C). In Fig. 7D, CTA number was remarkably higher in low- than high-risk group. Most immune checkpoints comprising CD80, CD86, IDO1, LAG3, LAIR1, PDCD1, HAVCR2, and LGALS3 presented the notably higher expression in high- versus low-risk group (Fig. 7E). Above findings proved that high-risk patients might respond to immunotherapy.
Heterogeneous sensitivity to drugs between low- and high-risk TNBC patients
Low-risk samples presented the notably lower IC50 of BI-2536, and UMI-77, indicative of higher sensitivity to the drugs (Fig. 7F). Meanwhile, high-risk individuals presented stronger sensitivity to 5-Fluorouracil, AZD2014, GSK2606414, AZD1332, Temozolomide, Oxaliplatin, Epirubicin, Taselisib, Entospletinib, Camptothecin, Dabrafenib, BMS-536924, PF-4708671, Topotecan, AZD8055, Rapamycin, Mitoxantrone, Dactolisib, Vinorelbine, Pictilisib, Foretinib, Ribociclib, GNE-317, Dactinomycin, Buparlisib, Teniposide, Irinotecan, JAK1_8709, Alisertib, AZD5363, AZ960, Palbociclib, LGK974, Gemcitabine, VX-11e, Uprosertib, KU-55933, Trametinib, SCH772984, JAK_8517, and Gallibiscoquinazole.
Experimental verification of the expression of the signature genes
We validated the expression of the signature genes in mammary epithelial cells (MCF10A) and two breast cancer cells (MDA-MB-231, and MCF-7). RT-qPCR results showed that CTSL expression was significantly lower in MCF-7 not MDA-MB-231 cells versus MCF10A cells (Fig. 8A). No significant difference in CTSD expression was detected between MCF10A cells and MDA-MB-231/MCF-7 cells (Fig. 8B). In comparison to MCF10A cells, lower ELK4 expression was found in MDA-MB-231, and MCF-7 cells (Fig. 8C). XRCC4 expression was remarkably higher in MDA-MB-231 not MCF-7 cells than MCF10A cells (Fig. 8D). In addition, HSPA8 was found to be remarkably lower in MCF-7 not MDA-MB-231 cells versus MCF10A cells (Fig. 8E). Western blot results demonstrated that CTSL, and CTSD expression was notably down-regulated and XRCC4 expression was notably up-regulated both in MDA-MB-231, and MCF-7 cells versus MCF10A cells (Fig. 8F). HSPA8 expression was modestly down-regulated in MDA-MB-231/MCF-7 cells versus MCF10A cells. However, ELK4 expression was not detected in the three cell lines.
Knockdown of XRCC4 attenuates migration of MCF10A, MDA-MB-231, and MCF-7 cells
For observing the role of XRCC4 in tumor cell migration, MCF10A, MDA-MB-231, and MCF-7 cells were transfected with specific siRNAs of XRCC4. Consequently, XRCC4 expression was significantly lowered in MCF10A, MDA-MB-231, and MCF-7 cells, notably si-XRCC4#1 and si-XRCC4#2 (Fig. 9A-C). Wound scratch assay results showed that MCF10A, MDA-MB-231, and MCF-7 cells with si-XRCC4#1 and si-XRCC4#2 had the significantly lower migration rate in comparison to those with si-NC transfection (Fig. 9D-I). Therefore, silencing XRCC4 attenuated migrative capacities of MCF10A, MDA-MB-231, and MCF-7 cells.
Overexpressed CTSL, CTSD, and HSPA8 mitigate migration of MCF10A, MDA-MB-231, and MCF-7 cells
CTSL, CTSD, and HSPA8 were significantly overexpressed by the corresponding gene overexpression plasmids in MCF10A, MDA-MB-231, and MCF-7 cells (Fig. 10A-C). As illustrated in wound scratch assay results, overexpressed CTSL, CTSD, and HSPA8 notably attenuated the migrative rate of MCF10A, MDA-MB-231, and MCF-7 cells (Fig. 10D-I). Thus, overexpressed CTSL, CTSD, and HSPA8 mitigated migrative capacities of MCF10A, MDA-MB-231, and MCF-7 cells.