Breast cancer (BC) is the most common malignant cancer with the highest incidence and mortality rate (15% of all tumor deaths) among women [1, 2]. Meanwhile, BC is a highly heterogeneous cancer with five major intrinsic molecular subtypes, namely, Luminal (40%), Normal-like (2%–8%), Luminal B (20%), human epidermal growth factor receptor 2 (HER2)-enriched (10%–15%) and Triple Negative (15%–20%) [3]. Triple-negative breast cancer (TNBC), showing a lack of human epidermal growth factor receptor-2 (HER2), estrogen receptor (ER) and progesterone receptor (PR) immunohistochemistry expression, is defined as HER2-, ER- and PR- negative, respectively [3, 4]. Compared to non-TNBC subtypes, TNBC is characterized as the most aggressive BC with frequent p53 mutation, higher proliferative and distant metastases rates (10%–24% of all BC) and poor prognosis (only 60% of 5-year survival rate) [5, 6]. TNBC exhibits different clinical and biological characteristics from other BC subtypes and occurs in younger women at a higher frequency, commonly presenting as an invasive ductal carcinoma with advanced stage [7, 8]. Common therapeutic regimens for BC include chemotherapy, radiotherapy, surgery and endocrine therapy, such as endocrine tamoxifen and trastuzumab, which could be used for HER2-positive patients but not for TNBC patients due to a lack of the three receptors [9]. Currently, combined chemotherapy of taxanes and anthracycline is the most common choice for TNBC patients; however, anthracyclines have irreversible toxicity to the heart, which encourages researchers to further study the molecular mechanism of TNBC to develop safer and more effective chemotherapy regimens [10].

Intratumoral heterogeneity of TNBC is a driving force not only to pathogenesis, but also to metastasis, treatment resistance and poor clinical outcomes [11]. Yeo and Guan [12] found that a single TNBC tumor tissue might consist of multiple subtypes, the compositions of which will change under chemotherapy pressure [13]. Although most primary TNBC cases are responsive to pre-operative chemotherapy, viable tumor cells remaining in the breast (pathologic complete response) are still closely related to the poor clinical outcomes in TNBC but not to ER+ BC [14], suggesting the involvement of a minor TNBC subpopulation in metastatic dissemination [15]. The development of single-cell RNA sequencing (scRNA-seq) technology has improved the understanding of complex molecular pathways and provided novel insight into genetic heterogeneity at the single-cell level for cancer research [16]. For instance, Patel et al. [17] revealed the cell state diversity and functional features of glioblastoma. The malignant cells of BC originate from epithelial cells of the breast and contribute to the complex tumor microenvironment (TME) via interacting with immune cells and stroma cells [18]. Karaayvaz et al. [15] found that the TNBC is composed of different subpopulations of epithelial cells, in which one epithelial cell subpopulation with distinct signatures could be associated with the long-term prognosis of TNBC patients. Thus, exploring the subtypes of TNBC could help understand the progression of TNBC. Iacono et al. [19] performed a global regulatory program to establish the gene regulatory networks (GRNs) based on the scRNA-seq data, facilitating the identification of critical regulators in disease. In addition, epithelial–mesenchymal transition (EMT) promotes the metastasis and invasion of cancer cells, partly due to the fact that epithelial cells can acquire the properties of cancer stem cells [20]. Special signaling factors in the TME, such as interleukin-6 (IL-6), epidermal growth factor (EGF), transforming growth factor-$\beta$ (TGF-$\beta$) and hematopoietic growth factor (HGF), activate the intracellular signaling pathways of epithelial cells [21, 22]. Particularly, SNAIL, TWIST and Zinc-finger enhancer binding (ZEB), as the key transcription factors in EMT, have been extensively studied [23] and the TGF-$\beta$ signaling pathway plays a vital role in EMT induction [24]. These three distinct extracellular TGF-$\beta$ isoforms bind to their complex receptor TGF-$\beta$ receptor type 1 (TGF-$\beta$R1) and TGF-$\beta$R2 to further phosphorylate suppressor of Mother against Decapentaplegic (SMAD2) and SMAD3 and recruit SMAD4 [25]. Subsequently, the SMAD2–SMAD3–SMAD4 complex migrates to the nucleus where they act as transcription factors or helper proteins to regulate the expression of genes involved in cell growth, proliferation, invasion, motility, acquisition of mesenchymal state, angiogenesis and apoptosis [26]. TGF-$\beta$ pathway could activate EMT through several different mechanisms, and the SMAD complex further activates the mesenchymal genes (fibronectin and vimentin) and EMT-TFs SNAIL, TWIST, ZEB1 and SLUG to suppress E-cadherin. Moreover, the TGF-$\beta$ pathway also induces the expression of a group of miRNAs and lncRNAs to promote EMT in carcinogenesis and fibrosis [27, 28].

This study performed a comprehensive analysis of scRNA-seq data from five primary BC samples to identify TNBC subtypes and determine critical genes for each molecular subtype. Firstly, the cell clusters were classified according to the annotated marker genes. Secondly, the pathological source of each cell cluster was classified, and we found that epithelial cells had a larger proportion in the TNBC pathotype. Pathway enrichment analysis revealed that the TGF-$\beta$ pathway played an essential role in TNBC progression, and we further identified TNBC-specific epithelial cell subsets and their gene markers. Ligand–receptor interaction analysis revealed the potential molecular mechanism of TNBC-specific epithelial cells in promoting cancer development. Finally, we developed the gene regulatory networks (GRNs), based on which c-fos gene (FOS) and X-box binding protein 1 (XBP1) were screened as two critical regulons in TNBC cell growth and proliferation.

TNBC is an aggressive subtype of BC with high intratumoral heterogeneity. In this study, we performed a comprehensive analysis of five primary BC samples using scRNA-seq data and identified nine cell subsets based on gene-expression profiles and marker genes, with the epithelial cells accounting for the largest proportion in nine subsets and TNBC pathotypes. Previous studies showed that malignant BC cells derive from epithelial cells [18, 37, 38]. We investigated the difference in epithelial cells in the TNBC and non-TNBC pathotypes, and classified four TNBC-specific epithelial cell subsets closely related to cell proliferation. Interaction analysis showed that these four subsets suppressed the immune response through different mechanisms to interact with T cells, playing vital roles in promoting tumor progression. In addition, the ANXA3 epithelial cell subset had the largest number of T cell interacted-receptor pairs and two key regulons involved in proliferation. Thus, we speculated that the ANXA3 epithelial cells might be a leading cause of TNBC progression, and the critical regulons genes could serve as potential targets for chemotherapy.

The epithelial cells in the TNBC pathotype exhibited higher TGF-$\beta$, G2M checkpoint, mTORC1, ROS, p53 and PI3K-mTOR-AKT pathway activities. Normally, the TGF-$\beta$ signaling is complex during cancer progression and it can have both tumor-promoting and tumor-suppressive functions [39, 40]. In organs such as the breast and skin, TGF-$\beta$ signaling actively regulates homeostatic growth, inhibits cell proliferation and transformation, and promotes apoptosis at the early stages of tumorigenesis [41]. However, the tumor cells could not only selectively bypass or adapt to the inhibitory functions of TGF-$\beta$ [42], but also make use of TGF-$\beta$ to obtain growth advantage and undergo EMT to allow epithelial cells to convert to migratory and invasive cells [43]. TNBC epithelial cells had higher TGF-$\beta$ pathway activity, suggesting that these epithelial cells might experience EMT to suppress anti-tumor immune reaction and promote cancer progression and metastasis. G2M checkpoint activation associated with DNA damage can further activate the p53 pathway to increase cell apoptosis [44]. Oncogenic mutations stimulate tumorigenesis by activating diverse growth signaling pathways [45]. Activation of the PI3K-mTOR-AKT pathway in the majority of tumors promotes cell growth and proliferation and changes the TME [46]. Continuous overactivation of mTORC1 signals could increase the levels of cell metabolism to promote cancer development [47]. Thus, the cell growth and proliferation of TNBC-epithelial cells were enhanced via multiple signaling pathways, although there existed an antagonism between the p53 pathway and the growth signaling pathway.

In the four TNBC-specific epithelial cell subsets, their marker genes were associated with the apoptotic, proliferation and cell growth processes. Lipids are involved in the proliferation, metastasis and therapeutic resistance in many tumors [48], and hypoxia can induce the formation of reactive oxygen species (ROS), leading to lipid peroxidation [49]. Stearoyl-CoA desaturase (SCD), a rate-limiting enzyme for unsaturated fatty acid synthesis, is highly expressed in lung cancer [50]. In tumors, SCD can reactivate fatty acid biosynthesis and improve the overall saturation level of membrane lipids, helping the tumor avoid damage from ROS [51]. MKI67 expression is indicative of cell proliferation status, and its overexpression is a prognostic indicator for cancer [52]. Annexin A3 (ANXA3) as a Ca²⁺-dependent membrane binding protein is involved in channel regulation, signal transduction, cell growth, proliferation, migration, apoptosis and anti-inflammation in tumor progression [53]. Research reported that silencing ANXA3 inhibits the proliferation, invasion and colony-forming abilities of HCC-1954 and MDA-MB 231 BC cells [54, 55]. AQP5 is an aquaporin that controls the water flow across the membrane and is also responsible for regulating osmotic pressure based on the intracellular glycerol concentration [56]. The regulation of water movement and osmotic pressure further influence cellular functions, such as proliferation, migration and adhesion [57]. AQP5 is involved in breast carcinogenesis [58]. AQP5 silencing or induction of hyperosmotic stress downregulates the expression of AQP5 and negatively affects tumor cell proliferation and migration [59], which might be caused by AQP5-mediated oxidative stress imbalance [60]. Another study reported that overexpression of AQP5 in the TNBC samples with high-expressed MKI67 is associated with a poor outcome, suggesting that the co-expression of the two genes could serve as an independent prognostic indicator for TNBC patients [61]. Thus, these four TNBC-specific epithelial cell subsets all contributed to TNBC progression but through different mechanisms, specifically, SCD+ and AQP5+ epithelial cell subsets were responsible for maintaining the oxidative stress homeostasis, and MKI67+ and ANXA3+ epithelial cell subsets were responsible for the cell cycle, proliferation and invasion regulation.

The ANXA3+ epithelial cell subsets interacted with T cells through different means, with the two regulons (FOS and XBP1) playing prominent roles in cell proliferation and angiogenesis. FOS has been regarded as an oncogene, and high levels of FOS in the nucleus can bind to other transcription factors (TFs) for downstream transcriptional activation [62]. A similar study showed that FOS is a representative apoptotic driver, as higher expressions of both FOS and apoptotic effector genes in FOS-targeted cells are correlated with better survival [63]. In addition, Nanamiya et al. [64] found that FOS is partially involved in regulating the expression of PVRL4 in BC cells, which in turn regulates the response of cancer cells to the immune system. XBP1 is a unique basic-region leucine zipper (bZIP) protein (TF) involved in EMT activation, cell growth and proliferation in a variety of cancers including in BC [65]. A previous study discovered that downregulating XBP1 expression through RNA interference can restore E-cadherin levels and promote the formation of cell–cell junctions in BC cells, which results in the inhibition of cell invasion and tumor growth [66]. These results further confirmed the reliability of the key regulatory genes we screened for TNBC, with FOS and XBP1 playing important roles in TNBC progression.

This study developed the single-cell profile of BC and identified four TNBC-specific epithelial cell subsets and comprehensively analyzed their roles in TNBC progression. In addition, the GRNs of the two regulons involved in cell growth and proliferation were constructed; however, their actual function and molecular mechanisms remained to be further elucidated in vivo and in vitro.

We performed a comprehensive analysis of the scRNA-seq data of five primary BC samples, and found a significant contribution of four TNBC-specific epithelial cell subsets to the TNBC metastasis and malignant progression. Two critical regulons in the ANXA3+ epithelial cell subsets exhibited complex regulation relationships to promote tumor development, which could help explore the underlying molecular mechanism of intratumoral heterogeneity and may be used as potential targets for TNBC.

TNBC, Triple-negative breast cancer; GEO, Gene Expression Omnibus; DAVID, Database for Annotation, Visualization and Integrated Discovery; SCENIC, Single-cell regulatory network inference and clustering; scRNA-seq, single cell RNA sequencing; GRNs, Gene regulatory networks; BC, Breast cancer; HER2, human epidermal growth factor receptor-2; ER, estrogen receptor; PR, progesterone receptor; EMT, epithelial-mesenchymal transition; IL-6, interleukin-6; EGF, epidermal growth factor; TGF-$\beta$, transforming growth factor-$\beta$; HGF, hematopoietic growth factor; PCA, Principal Component Analysis; UMAP, Uniform Manifold Approximation and Projection; MSigDB, Molecular Signatures Database.

The datasets generated and/or analyzed during the current study are available in the [GSE180286] repository, [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE180286]. The raw data used and/or analyzed during the current study are available from the corresponding author on reasonable request.

All authors contributed to this present work: HZ & YNS designed the study, XND acquired the data. HZ improved the figure quality. HZ & YNS drafted the manuscript, HZ & XND revised the manuscript. All authors contributed to editorial changes in the manuscript. All authors read and approved the final manuscript. All authors have participated sufficiently in the work and agreed to be accountable for all aspects of the work.

Not applicable.

Not applicable.

This work was supported by the Yancheng Engineering Center (YC2022808).

The authors declare no conflict of interest.