ORIGINAL RESEARCH article

Front. Immunol., 01 July 2025

Sec. Cancer Immunity and Immunotherapy

Volume 16 - 2025 | https://doi.org/10.3389/fimmu.2025.1618193

This article is part of the Research TopicUnraveling Breast Cancer Complexity: Insights from Single-Cell Sequencing and Spatial TranscriptomicsView all 8 articles

Identification and validation of genes related to stem cells and telomere maintenance mechanisms as biomarkers for breast cancer

Shuang Zhen&#x;Shuang ZhenLifeng Huang&#x;Lifeng HuangQiannan ZhuQiannan ZhuRui ChenRui ChenJue WangJue WangXiaoming Zha*&#x;Xiaoming Zha*‡
  • Department of Breast Surgery, Department of General Surgery, The First Affiliated Hospital with Nanjing Medical University, Nanjing, China

Background: Stem cell-related genes (SCRGs) and telomere maintenance mechanism-related genes (TMMRGs) are pivotal in breast cancer (BC) pathogenesis by facilitating tumor cell proliferation and self-renewal. This study employed integrated transcriptomic and single-cell RNA sequencing (scRNA-seq) analyses to investigate SCRGs and TMMRGs as potential biomarkers for BC and to elucidate their underlying cellular mechanisms.

Methods: Total RNA was extracted from eight BC tumor samples and eight matched adjacent non-tumorous tissues. Differential expression profiling, protein-protein interaction (PPI) network construction, and Molecular Complex Detection (MCODE) were conducted. Biomarker candidates were identified using the least absolute shrinkage and selection operator (LASSO) algorithm, followed by pathway enrichment and immunological analyses. Publicly available scRNA-seq datasets were utilized to delineate BC cell types, with emphasis on cellular subsets exhibiting differential biomarker expression. Heterogeneity, communication, and pseudo-temporal analyses of key cells were examined. Biomarker expression was further validated by reverse transcription-quantitative polymerase chain reaction (RT-qPCR).

Results: JUN, NFKB1, and SP1 were significantly downregulated in BC, potentially modulating disease progression through mechanisms involving extracellular matrix (ECM) remodeling, intracellular signaling, oxidative stress response, and translational regulation. Activated B cells and natural killer (NK) cells demonstrated elevated infiltration levels, accompanied by increased expression of immune checkpoint molecules CD200, CD274, TIGIT, TNFRSF25, and TNFSF15. Nine distinct cellular lineages were annotated, among which mesenchymal cells exhibited pronounced biomarker expression differences and enhanced differentiation potential, designating them as key cellular mediators. Interactions between mesenchymal subpopulations (MSC1, MSC2, MSC3) and other cell types were markedly reduced in BC, despite an overall expansion in mesenchymal cell numbers during disease progression. MSC1 emerged as the predominant subtype. RT-qPCR analyses corroborated the downregulation of JUN, NFKB1, and SP1 in BC tissues.

Conclusion: JUN, NFKB1, and SP1 were identified as potential biomarkers for BC. These findings highlight the critical role of mesenchymal cells in tumor biology and suggest potential therapeutic targets.

1 Introduction

Breast cancer (BC) remains one of the most prevalent malignancies among women globally, with incidence rates continuing to escalate and posing a significant public health burden (1). Established risk factors—including age, sex, ethnicity, family history, genetic mutations, timing of menarche and menopause, and hormone replacement therapy—contribute to the diagnostic and therapeutic complexity of BC, adversely impacting patient survival (14). Although therapeutic advances have been made, five-year survival outcomes are still largely contingent upon the timeliness and accuracy of disease detection (5). The multifaceted nature of BC and its unpredictable clinical course highlight the urgent need for novel biomarkers to enhance early diagnosis, clarify underlying mechanisms, and support individualized treatment approaches, ultimately improving therapeutic outcomes and patient well-being (4, 6, 7). Consequently, intensified efforts in biomarker research are essential not only for advancing the biological understanding of BC but also for informing preventive, diagnostic, and therapeutic innovations—crucial for optimizing global healthcare resource allocation.

Telomeres, composed of repetitive nucleotide sequences at chromosomal termini, preserve genomic stability by preventing degradation and end-to-end fusion during mitosis. Progressive telomere shortening with successive cell divisions ultimately triggers replicative senescence or apoptosis upon reaching a critical threshold (8). Studies have shown that the shortening of telomeres is closely related to the senescence of immune cells, which can lead to the exhaustion of T cell function and the decline of proliferation ability, thereby weakening the immune response of the body (911). Telomere maintenance mechanisms (TMM), predominantly regulated by telomerase (12), counteract this attrition, enabling sustained cellular proliferation (13).The maintenance of telomere length is crucial for the normal function of cells, and telomere dysfunction is associated with various diseases (14). In cancer research, the activation of TMM is strongly linked to the occurrence and development of tumors. For instance, TMM activation in malignant neuroblastoma is closely associated with high-risk diseases and poor prognosis (15). Studies on hepatocellular carcinoma (HCC) have also shown that telomere maintenance plays a key role in the occurrence mechanism of HCC and helps define the clinical characteristics of patients (16). Stem cells, defined by their capacity for self-renewal and multilineage differentiation, are pivotal in tissue homeostasis, regeneration, and development (17, 18). Due to the unique immune escape mechanism of stem cells, they are crucial in the occurrence and development of tumor and immune therapy resistance (1921). For example, tumor stem cells can down-regulate the expression of major histocompatibility complex class I (MHC-I) molecules, therefore reducing antigen presentation and making it difficult for cytotoxic T lymphocytes to recognize and kill tumor cells (22). The interplay between TMM, stem cells, and BC is particularly significant (23, 24). Dysregulated TMM activity, frequently observed in malignancies including BC, facilitates the bypass of senescence and confers unlimited replicative potential to tumor cells (25). Numerous studies have shown that telomerase reactivation is a common strategy employed by cancer cells to stabilize telomeres and circumvent programmed cell death. Additionally, cancer stem cells—a distinct subset with tumor-initiating and self-renewing properties—exhibit robust telomere maintenance activity and engage in oncogenic signaling cascades that drive tumor progression and metastasis (26, 27). TMM supports the proliferative and regenerative capacity of these cells, thereby contributing to oncogenesis and therapeutic resistance. Despite substantial evidence linking TMM and stem cells to BC pathogenesis, the precise molecular pathways mediating their effects remain incompletely elucidated (28).

Single-cell RNA sequencing (scRNA-seq) can reveal the dynamic changes, functional status and intercellular interactions of tumor-infiltrating immune cells at the single-cell level, providing key insights in clarifying the mechanism of anti-tumor immune responses (29, 30). At present, many studies have been conducted on tumors through scRNA-seq sharing. Liu et al. revealed the immunosuppressive effect of APOE+ macrophages in immune checkpoint inhibitor therapy through scRNA-seq analysis (31). Xie et al. elucidated the landscape of BC brain metastases through scRNA-seq analysis and identified ILF2 as a potential therapeutic target (32). Meanwhile, they also clarified the contribution of intercellular communication of the tumor microenvironment cells guided by histone chaperones to BC metastasis through scRNA-seq analysis (33). These studies indicate that single-cell sequencing technology can systematically analyze the heterogeneity of the tumor immune microenvironment, providing pivotal data for revealing treatment resistance mechanisms and discovering new therapy targets.

BC is orchestrated by intricate genetic and cellular networks, wherein stem cell-related genes (SCRGs) and TMM-related genes (TMMRGs) play pivotal roles in sustaining tumor cell renewal and unchecked growth. Elucidating the function of these genes at both the transcriptomic and single-cell resolution offers a promising strategy for uncovering actionable biomarkers and therapeutic targets. This study integrated transcriptomic profiling and scRNA-seq to characterize SCRGs and TMMRGs in BC. Tumor and matched adjacent normal tissue samples from patients with BC, alongside publicly available scRNA-seq datasets, were analyzed to assess gene expression patterns and functional relevance. By delineating key biomarkers and mapping their expression across diverse cellular subsets, this investigation aims to clarify their roles in BC progression. The resulting insights enhance the molecular understanding of BC and inform the development of precision therapies aimed at improving clinical outcomes.

2 Materials and methods

2.1 Data source

Tumor tissue samples were collected from eight patients with BC at The First Affiliated Hospital of Nanjing Medical University, with matched paracancerous tissues from eight healthy individuals serving as controls. All specimens underwent RNA-seq and stringent quality control (QC) to form the whole transcriptome RNA-seq dataset. The study was approved by the Ethics Committee of The First Affiliated Hospital of Nanjing Medical University, and informed consent was obtained from all participants.

To enhance analytical robustness, additional BC-related datasets were sourced from public repositories. Specifically, RNA-seq data comprising 1,082 BC tumor samples and 113 paracancerous samples were retrieved from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) for expression validation. Concurrently, the scRNA-seq dataset GSE245601, based on the GPL18573 platform, was downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/gds), including 10 BC tumor samples and 2 normal tissue samples for cellular-level expression analysis (34).

Stem cell-associated gene data were curated by referencing established literature (35), and 26 stem cell gene sets were extracted from the StemChecker database (http://stemchecker.sysbiolab.eu/), yielding 4,419 SCRGs. Additionally, 218 TMMRGs were incorporated based on published literature (13). The overall analytical workflow of this study was illustrated in Figure 1.

Figure 1
Flowchart depicting analysis of tumor and control tissue samples. It starts with RNA-seq, dividing into differentially expressed genes, telomere-related genes, and stem cell-related genes. These lead to GO and KEGG analysis, protein interactions, and scRNA-seq analysis. MCODE plugin and machine learning identify biomarkers. Steps include cell clustering, key cell identification, and pseudo-time analyses, ending with immune microenvironment analysis, regulatory network, GSEA, and RT-qPCR.

Figure 1. The analysis flowchart of this study.

2.2 RNA-seq and data preprocessing

Total RNA from the 8 tumor and 8 paracancerous tissue samples was isolated using TRIzol reagent (Invitrogen, CA, USA). RNA integrity and concentration were evaluated using a NanoDrop ND-1000 spectrophotometer (Wilmington, DE, USA) and an Agilent Bioanalyzer 2100 system (Agilent, CA, USA). Sequencing libraries were prepared with the Hieff NGS Ultifaillumina Dual-mode mRNA Library Prep Kit to generate inserts of 300 ± 50 bp, followed by high-throughput paired-end sequencing (PE150) on the Illumina NovaSeq 6000 platform.

Raw sequencing data in FASTQ format underwent QC using fastp software (v0.23.4) (36) with default parameters, including adapter trimming, deduplication, and removal of low-quality reads. Cleaned reads were aligned to the human reference genome (GRCh38) using HISAT2 (v2.2.1) (37), and alignment outputs were stored in BAM format. Transcript assembly and quantification were performed using StringTie (v2.2.0) (38), and expression levels were normalized to FPKM, calculated as:

FPKM=total exon fragments/(mapped reads in millions×exon length in kilobases).

2.3 Construction of differential expression analysis and enrichment analysis

Differential expression analysis between tumor and normal tissues was conducted on the transcriptome-wide RNA-seq data using the DESeq2 package (v1.38.0) (39). Differentially expressed genes (DEGs) were defined by |log2 fold change (FC)| > 0.5 and P < 0.05. Visualization of DEGs was performed using a volcano plot and heatmap, generated with ggplot2 (v3.3.6) (40) and ComplexHeatmap (v2.14.0) (41), respectively. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using the clusterProfiler package (v4.7.1.3) (42), with significance set at P < 0.05. GO terms were classified into three categories: biological processes (BPs), molecular functions (MFs), and cellular components (CCs).

2.4 Identifying hub genes through protein-protein interaction network visualization

The intersection of DEGs, SCRGs, and TMMRGs was identified using the VennDiagram package (v1.7.1) (43) to isolate candidate genes implicated in both stemness and telomere regulation within the context of BC. These candidate genes were subsequently analyzed using the Search Tool for the Retrieval of Interacting Genes (STRING) database (http://string-db.org) with a confidence threshold set at > 0.7, and a protein–protein interaction (PPI) network was generated to elucidate potential molecular interactions at the protein level. The Molecular Complex Detection (MCODE) plugin was used to select the highest-scoring cluster and identify hub genes. A PPI network was subsequently constructed for these hub genes. Both networks were generated using Cytoscape software (v 3.9.1) (44).

2.5 Combining machine learning algorithms and expression profiling to screen biomarkers

To assess the diagnostic potential of the hub genes, Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis was performed using the glmnet package (v4.1.4) (45), configured with the ‘binomial’ family parameter. The optimal model was established at the point of minimal model error and the lowest lambda value, at which the genes as candidate biomarkers were selected for expression analysis. These genes were then subjected to expression analysis in both the in-house whole transcriptome RNA-seq dataset and the TCGA-BC dataset, comparing tumor and normal tissues. Expression differences were visualized via box plots generated with the ggplot2 package, and genes exhibiting consistent and statistically significant differential expression across both datasets (P < 0.05) were selected as biomarkers.

2.6 Gene set enrichment analysis

To investigate the functional relevance of these biomarkers, enrichment analyses were conducted on the whole transcriptome dataset. Using the Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb), the curated gene set collection c2.cp.v2023.2.Hs.symbols.gmt was employed. For each biomarker, Spearman correlation coefficients were computed with all other genes, and genes were ranked in descending order of correlation. GSEA was performed using the clusterProfiler package, with statistical thresholds set at P < 0.05 and False Discovery Rate (FDR) < 0.05. The top five significantly enriched pathways per biomarker were visualized using the enrichplot package (v1.18.4) (46), ranked by enrichment significance.

2.7 Analyzing immune microenvironment differences and immune checkpoint expression profiling

To further explore the differences in the immune microenvironment between patients with BC and normal individuals, a single-sample GSEA (ssGSEA) approach from the GSVA package (v1.42.0) (47) was applied to estimate the infiltration levels of 28 predefined immune cell types (48) within each sample of the whole transcriptome dataset. Comparative analysis revealed statistically significant differences in immune cell abundance between tumor and normal samples (P < 0.05). Additionally, a set of 40 immune checkpoints was compiled from previously published literature (49), and differences in the expression of these checkpoints between tumor and normal samples were analyzed within the same dataset (P < 0.05).

2.8 Drug prediction analysis and construction of regulatory network

Drug–biomarker interactions were identified using the Drug Gene Interaction Database (DGIDB, https://dgidb.org), with the top five candidates ranked by interaction score. The resulting drug–biomarker network was visualized using Cytoscape.

To investigate the underlying molecular mechanisms of the biomarkers, the DESeq2 package was employed to analyze differentially expressed miRNAs and lncRNAs (DE-miRNAs and DE-lncRNAs) from a whole transcriptome RNA-seq dataset, comparing tumor and normal tissues under the criteria of P < 0.05 and |log2FC| > 0.5. Volcano plots were generated to depict the distribution of DE-miRNAs and DE-lncRNAs.

Candidate miRNAs targeting the biomarkers were predicted via miRDB (https://mirdb.org/) and TargetScan (https://www.targetscan.org/) using a prediction score threshold of 500,000. Correspondingly, target lncRNAs for these miRNAs were identified through starBase (http://starbase.sysu.edu.cn/index.php), filtered by a clipExpNum > 4.

Subsequent analysis involved intersecting DE-miRNAs with the predicted miRNA targets and DE-lncRNAs with the predicted lncRNA targets to isolate key regulatory miRNAs and lncRNAs. A lncRNA–miRNA–mRNA regulatory network was then constructed in Cytoscape to elucidate the multifaceted regulatory landscape influencing biomarker expression.

2.9 ScRNA-seq analysis

To investigate the cellular mechanisms underlying BC and to characterize biomarker heterogeneity at the single-cell level, comprehensive analyses were conducted using the GSE245601 dataset. QC, clustering, and annotation were performed with the Seurat package (v5.0.1) (50). Key quality metrics—including gene count, cell count, and mitochondrial gene percentage—were calculated. Cells with fewer than 200 genes and genes covered by fewer than three cells were filtered out, as were cells expressing either less than 200 or more than 5000 genes. Genes expressed in fewer than three cells or with total counts below 500 or above 20,000 were filtered out. Mitochondrial content was capped at 10%. Distributions of nFeature_RNA, nCount_RNA, and percent.mt before and after filtering were visualized, and only cells meeting all criteria were retained for downstream analysis. Subsequently, data normalization was applied, and the top 2,000 highly variable genes were selected using the vst method. Principal component analysis (PCA) was conducted to evaluate variance across components, with results visualized via a PCA elbow plot. Statistically significant principal components, identified using the JackStrawPlot function based on gene-level P-values, were selected for further analysis. Unsupervised clustering was performed using the FindNeighbors and FindClusters functions, and clusters were visualized using Uniform Manifold Approximation and Projection (UMAP) at a resolution of 0.4. Cell type annotation was based on canonical marker genes reported in prior research (34), and a bubble plot was generated to illustrate the specificity of these markers across identified cell populations. In addition, the NF-κB activation, stemness-related signatures, and ERK-JUN signaling pathway scores of the annotated cell types were evaluated using the AddModuleScore method.

2.10 Identification and heterogeneity analysis of key cells

In the scRNA-seq analysis, the AddModuleScore function from the Seurat package was utilized to compute enrichment scores of specific gene sets for each sample in the GSE245601 dataset, thereby assessing metabolic pathway activity across cell populations. A curated background gene set, encompassing candidate genes associated with both stem cells and TMM in BC, served as the basis for this scoring. Scores were calculated across annotated cell types in tumor and normal tissues, and cell populations demonstrating significant differences (P < 0.05) were designated as candidate key cells. Subsequent single-cell-level analysis examined the expression profiles of biomarkers within these candidate key cells across tumor and normal tissues. Cells exhibiting consistent and significant differential expression across all biomarkers (P < 0.05) were defined as key cells relevant to both stemness and telomere maintenance in BC. Additionally, heterogeneity within these key cells was examined. UMAP-based reclustering was performed using FindNeighbors and FindClusters to delineate cellular subtypes, followed by assessment of biomarker expression across the resulting subclusters (P < 0.05).

2.11 Cell communication and pseudo-time analyses of key cells

To infer intercellular interactions, cell–cell communication analysis was conducted by evaluating receptor–ligand expression and pairing patterns among cell types. The distribution and proportions of key cell subtypes relative to other populations in tumor and normal tissues were quantified. Communication networks were constructed using the CellChat package (v1.6.1) (51), generating aggregated intercellular signaling maps and visualizing the contribution of each cell cluster to the overall communication landscape. Developmental trajectories and lineage dynamics of key cells and their subtypes were further explored through pseudo-time analysis using the Monocle package (v2.26.0) (52).

2.12 Expression validation of biomarkers

Tumor tissue samples were obtained from five patients diagnosed with BC, alongside five adjacent normal tissue samples from healthy individuals, all collected at The First Affiliated Hospital of Nanjing Medical University. RT-qPCR was performed for experimental validation of biomarker expression. Total RNA was extracted from the ten samples using TRIzol reagent (Ambion, Austin, USA), and cDNA synthesis was carried out using the SureScript First-Strand cDNA Synthesis Kit (Servicebio, Wuhan, China), following manufacturer protocols. qPCR amplification was performed with the 2 × Universal Blue SYBR Green qPCR Master Mix (Servicebio, Wuhan, China). GAPDH served as the internal control, and gene expression was quantified using the 2−ΔΔCt method (53). Primer sequences are provided in Supplementary Table 1.

2.13 Statistical analysis

All statistical analyses were conducted in R software (v4.2.2). Group comparisons were performed using the Wilcoxon test, with P < 0.05 considered statistically significant.

3 Results

3.1 Identification of candidate genes and associated pathways

Differential expression analysis of the whole transcriptome RNA-seq dataset identified 8,193 DEGs between tumor and normal samples, including 7,014 upregulated and 1,179 downregulated genes in BC (Figures 2A, B). These DEGs underwent enrichment analysis, yielding 649 significantly associated GO terms: 404 BPs such as “cell cycle,” “DNA damage response,” and “signal transduction” (Figure 2C); 112 CCs, including “plasma membrane,” “extracellular exosome,” and “cell junction” (Figure 2D); and 133 MFs, including “protein binding,” “protein kinase activity,” and “DNA-binding transcription factor activity” (Figure 2E). Additionally, 68 KEGG pathways enriched in DEGs were identified, including “pathways in cancer,” “PI3K-Akt signaling pathway,” “proteoglycans in cancer,” “MAPK signaling pathway,” and “focal adhesion” (Figure 2F). These pathways reflect essential mechanisms contributing to BC pathogenesis, reinforcing the reliability of differential expression results and providing functional insight into molecular alterations associated with BC. By intersecting the 8,193 DEGs with 4,419 SCRGs and 218 TMMRGs, 33 candidate genes related to both stemness and TMM in BC were identified (Figure 2G).

Figure 2
The image contains multiple panels with graphs and charts related to gene expression analysis. Panel (A) is a volcano plot showing the relationship between fold change and p-value, highlighting upregulated and downregulated genes. Panel (B) is a heatmap illustrating gene expression levels across various samples. Panels (C) to (F) are circular charts depicting the distribution of biological processes, cellular components, molecular functions, and pathways, each labeled with specific categories and percentages. Panel (G) is a Venn diagram showing the overlap between differentially expressed genes (DEGs), specific gene clusters (SCRGs), and thermogenesis-related genes (TMMRGs).

Figure 2. Identification of candidate genes and associated pathways. (A) Volcano plots of DEGs between BC samples and normal samples. In the figure, each dot represents a gene. Pink indicates up-regulation of gene expression, while blue indicates down-regulation of gene expression. (B) The expression heat map of DEGs. Red represents high expression and blue represents low expression. The darker the color, the higher/lower the expression. (C-E) GO analysis with up-regulated (the outer circle) and down-regulated (the inner circle) DEGs, including biological process(C), cellular component (D), and molecular function (E). (F) KEGG analysis with up-regulated (the inner circle) and down-regulated (the inner circle) DEGs. (G) 33 overlapping genes were identified as both stem cell and telomere maintenance mechanism related DEGs in BC.

3.2 Selection of JUN, NFKB1, and SP1 as biomarkers for BC

A PPI network encompassing 26 nodes and 77 edges was constructed based on these 33 candidate genes using a publicly available database (Figure 3A). Further analysis with the MCODE plugin identified the most prominent functional module, comprising 9 hub genes: JUN, HDAC1, HSP90AB1, SP1, ESR1, NFKB1, MYC, STAT3, and TP53. An additional PPI network was generated to visualize the interactions among these hub genes (Figure 3B).

Figure 3
(A) Network diagram showing interactions among various genes, with labels such as TERF1, PML, and EPAS1. (B) A simplified network highlighting key genes like HSP90AB1 and STAT3. (C) Line graph of coefficients versus log lambda values, displaying multiple lines. (D) Plot of binomial deviance over log(Lambda) with error bars. (E) Box plots comparing gene expression levels of JUN, NFKB1, and others between normal and tumor groups. (F) Additional box plots for JUN, NFKB1, and others showing expression in the same groups.

Figure 3. Screening of JUN, NFKB1, and SP1 as biomarkers for BC. (A) Protein–protein interaction network of 33 DEGs. (B) A significant module containing 9 hub genes was selected using MCODE. (C, D) The gene coefficient plot of LASSO analysis and the error plot of cross-validation. The horizontal coordinates were all log(Lambda), and the vertical coordinates were respectively the coefficient of the gene and the error of cross-validation. LASSO algorithm indicates the optimal model contained 5 genes. (E, F) The expression levels of these genes in tumor and normal samples were assessed in both whole transcriptome RNA-seq dataset (E) and TCGA-BC dataset (F) by Wilcoxon rank sum test. The tumor samples were colored red, and the normal samples were colored blue. ns represents P > 0.05, ** represents P < 0.01, *** represents P < 0.001, **** represents P < 0.0001.

Subsequently, the nine hub genes were subjected to LASSO regression analysis. At log(lambda.min) = −7.49352 and log(lambda.1se) = −5.725879, the optimal model retained five genes: JUN, SP1, NFKB1, STAT3, and TP53 (Figures 3C, D). Expression profiling across both the whole transcriptome RNA-seq dataset and the TCGA-BC cohort identified JUN, NFKB1, and SP1 as exhibiting statistically significant and consistent expression patterns (P < 0.01) (Figures 3E, F). These three genes were subsequently designated as biomarkers associated with both stem cell function and TMM in BC. Notably, all three demonstrated markedly reduced expression in tumor tissues compared to normal counterparts (P < 0.01).

3.3 Exploring the pathways by which biomarkers influence BC development

GSEA was performed to elucidate biological pathways associated with the identified biomarkers. Due to the lack of significant pathway enrichment for NFKB1, downstream analysis focused on pathways enriched by JUN and SP1. Enrichment analysis revealed that JUN was predominantly associated with pathways such as “ECM organization,” “RhoA regulation pathway,” “Nrf2 pathway,” “S1P–S1P3 pathway,” and “ECM proteoglycans” (Figure 4A). In contrast, SP1 enriched pathways included “selenoamino acid metabolism,” “SRP-dependent cotranslational protein targeting to membrane,” “translation initiation via medicus reference,” “ribosome,” and “eukaryotic translation elongation” (Figure 4B). These pathways are involved in key processes including ECM structure and tissue reconstruction, cell signaling regulation, oxidative stress responses, and protein synthesis and transport, all of which may contribute to BC pathogenesis through their interaction with these biomarkers.

Figure 4
Panel (A) shows a line graph of running enrichment scores for specific terms across an ordered dataset. Panel (B) displays a similar graph with different terms. Panel (C) is a stacked bar chart comparing immune cell types in normal versus tumor samples. Panel (D) presents a box plot of ssGSEA scores for various pathways in both groups. Panel (E) features a box plot comparing gene expression levels between normal and tumor groups across several genes.

Figure 4. Reveal of signaling pathways associated with biomarkers and potential interaction of immune system in BC patients. (A, B) GSEA enrichment results of JUN (A) and SP1 (B). Each polyline represents a pathway, and the peak of each polyline is the enrichment fraction of that pathway. (C) The relative abundance of 28 immune cells. (D) The differences in ssGSEA scores of 28 immune cell types between tumor samples and normal samples were compared by Wilcoxon rank-sum test. (E) The expression differences of immune checkpoints between tumor samples and normal samples were compared by Wilcoxon rank-sum test. ns represents P > 0.05, * represents P < 0.05, ** represents P < 0.01.

3.4 Revealing altered immune microenvironments and activated immune checkpoints in patients with BC

Comprehensive analysis of the immune microenvironment in BC revealed distinct immune infiltration patterns, as visualized in a stacked bar chart illustrating the relative abundance of 28 immune cell types across individual samples, computed via the ssGSEA algorithm (Figure 4C). Notably, activated B cells and NK cells demonstrated statistically significant differences between tumor and normal tissues (P < 0.05), with both cell types showing elevated infiltration in tumor samples (Figure 4D). This heightened immune presence likely reflects immune activation within the tumor microenvironment. Specifically, activated B cells may contribute to antibody-mediated responses, while NK cells are known for their direct cytotoxic activity against malignant cells, indicating that differential immune cell infiltration may be associated with immune escape, tumor growth, and the broader pathophysiology of BC.

Expression profiles were available for 34 out of the 40 immune checkpoints identified in the literature. An analysis of these immune checkpoints in tumor and normal samples revealed significant differences in the expression of CD200, CD274, TIGIT, TNFRSF25, and TNFSF15 between the two groups (P < 0.05) (Figure 4E). These five immune checkpoints also exhibited elevated expression levels in tumor samples (P < 0.05). This suggests that these checkpoints may be activated within the tumor microenvironment, potentially aiding tumor cells in evading immune surveillance and clearance, thereby positioning them as promising therapeutic targets. Targeting these immune checkpoints could inhibit tumor growth and metastasis.

3.5 Investigating targeted drugs and potential molecular mechanisms for biomarkers

Utilizing DGIdb, drug–biomarker interactions were predicted and ranked by interaction scores (Supplementary Table 2). Based on the selected candidates, a drug–biomarker interaction network was constructed, comprising 14 nodes and 11 edges (Figure 5A).

Figure 5
(A) A network diagram showing interactions of genes JUN, NFKB1, and SP1 with various chemicals. (B) A volcano plot illustrating differential expression of miRNAs, highlighting significant changes. (C) Another volcano plot labeled AVSP, depicting data points divided into upregulated, non-differential, and downregulated categories. (D) A network diagram displaying connections between various genes and miRNAs, including SP1 and JUN.

Figure 5. Investigating targeted drugs and potential molecular mechanisms for biomarkers. (A) The drug-biomarker network selected by DGIdb. Biomarkers were represented in red, and drugs were represented in orange. (B, C) The volcano plot of differentially expressed miRNAs (B) and lncRNAs (C) between tumor and normal samples. The color pink represented upregulation, while the color blue indicated downregulation. The horizontal axis corresponded to log2FC, denoting the logarithmic value of the gene expression change ratio, and the vertical axis represented the -log10 of the P-value. (D) Biomarker-miRNA-lncRNA regulatory network. Red represented biomarkers, cyan indicated miRNAs, and orange denoted lncRNAs.

Parallel differential expression analysis identified 272 DE-miRNAs (248 upregulated, 24 downregulated) and 12,287 DE-lncRNAs (11,899 upregulated, 388 downregulated) between tumor and normal samples (Figures 5B, C). Predictions of miRNAs targeting the biomarkers, along with their corresponding target lncRNAs, were obtained from public databases. Intersections between DE-miRNAs and predicted target miRNAs, as well as between DE-lncRNAs and predicted target lncRNAs, enabled the identification of key miRNAs and lncRNAs. These regulatory elements were integrated into a lncRNA–miRNA–mRNA network consisting of 18 nodes and 31 edges (Figure 5D), including 7 miRNAs, 9 lncRNAs, and 2 biomarkers (SP1 and JUN). Within this network, complex regulatory relationships were delineated; for example, XIST may regulate JUN via hsa-miR-139-5p, and potentially modulate SP1 through hsa-miR-199a-3p, hsa-miR-7-5p, hsa-miR-324-5p, and hsa-miR-299-5p.

3.6 All 9 annotated cell types exhibited higher counts in BC

Following transcriptomic analysis, single-cell studies provided insights into the unique expression patterns of biomarkers across various cell types and investigated the underlying cellular mechanisms of BC. Data QC was initially performed on the GSE245601 dataset, with metrics for nFeature_RNA, nCount_RNA, and percent.mt visualized before and after QC (Supplementary Figures 1A, B). Initially, the dataset contained 69,275 cells and 25,558 genes, which were reduced to 55,204 cells post-filtering, retaining the same number of genes. The top 2,000 highly variable genes were selected for further analysis (Supplementary Figure 1C). PCA demonstrated that cells from both tumor and normal samples were concentrated, with minimal variation beyond the 30th principal component, leading to the selection of the first 30 components (dims = 30) and a clustering resolution of 0.4 (Supplementary Figures 1D, E). UMAP clustering revealed 21 distinct clusters (Figures 6A, B). Marker gene analysis identified nine cell types: plasmablasts, mast cells, B cells, myeloid cells, basal cells, endothelial cells, mesenchymal cells, T cells, and luminal cells (Figure 6C). The specificity of these marker genes was validated using a bubble plot, enabling precise annotation of cell types (Figure 6D). The distribution of these nine cell types in tumor and normal samples was also visualized, revealing similar clustering patterns between groups but with higher cell counts in tumors across all types compared to normal samples (Supplementary Figure 2). These differences are likely attributed to the tumor microenvironment and the activation of immune responses in patients with BC. Additionally, basal cells had higher scores in the ERK-JUN signaling and stemness-related signatures pathways, while myeloid cells had the highest score in the NF-κB activation pathway (Supplementary Figure 3). This suggested that the proliferation, differentiation, and stemness regulation of basal cells, as well as the inflammatory response of myeloid cells, might play important roles in the biological processes of BC.

Figure 6
Four-panel image depicting UMAP plots and a dot plot analysis. Panel A and B show UMAP plots with clustering of cells labeled by sample identifiers. Panel C displays a UMAP plot with clusters labeled by cell type, including B cells, T cells, and others. Panel D is a dot plot representing gene expression across different cell identities, with varying dot sizes and colors indicating percent expression and average expression levels, respectively.

Figure 6. Cell clustering analysis and annotation. (A, B) UMAP clustering analysis identified 21 clusters. (C, D) The UMAP (C) and the bubble plot (D) of 9 distinct cell types annotated by marker genes.

3.7 Defining mesenchymal cells as a key cell population

Subsequent analyses employed the AddModuleScore function, using the previously identified 33 candidate genes associated with stem cells and TMM in BC as the background gene set. The scores for the nine annotated cell types in tumor and normal samples were calculated and compared. Noteworthy differences in AddModuleScore were observed among endothelial cells, luminal cells, mesenchymal cells, and myeloid cells (P < 0.01), indicating alterations in cell function and state within the tumor microenvironment. These cells were thus classified as candidate key cells (Figure 7A). Specifically, endothelial, luminal, and mesenchymal cells exhibited significantly lower scores in tumor samples (P < 0.0001), suggesting suppression of their stem cell properties and TMM within the tumor microenvironment. This suppression may be linked to tumor aggressiveness, metastatic potential, or treatment resistance. In contrast, myeloid cells displayed significantly higher scores in tumor samples (P < 0.01), potentially reflecting their active role in supporting tumor growth and maintaining the tumor microenvironment, such as by promoting inflammation, enhancing tumor cell survival, or suppressing immune responses.

Figure 7
Grouped data visualization comprising four panels: (A) Box plots comparing AddModuleScore for various cell types between normal and tumor groups, with statistical significance indicated. (B) Violin plots showing JUN expression across endothelial, luminal, mesenchymal, and myeloid cells in normal versus tumor groups, with significance levels noted. (C) Violin plots for NFKB1 expression in the same cell types and conditions. (D) Violin plots illustrating SP1 expression, again comparing normal and tumor groups across the four cell types with statistical annotations.

Figure 7. Analysis of key cell biomarkers. (A) AddModuleScore score of different cell types in tumor and normal samples. (B-D) The expression patterns of JUN (B), NFKB1 (C), and SP1 (D) in these candidate key cells in tumor and normal samples. ns represents P > 0.05, * represents P < 0.05, ** represents P < 0.01, *** represents P < 0.001, **** represents P < 0.0001.

Further investigation of the expression patterns of JUN, NFKB1, and SP1 in these key cell types revealed significant differences in both mesenchymal and luminal cells (P < 0.01) (Figures 7B–D).

3.8 Recognition of mesenchymal cell subtypes: MSC1, MSC2, and MSC3

Subsequent investigations examined the heterogeneity of mesenchymal cells, with UMAP clustering categorizing these cells into three subtypes: MSC1, MSC2, and MSC3 (Figure 8a). The distribution of these subtypes in tumor and normal tissues revealed that, compared to normal samples, MSC1 and MSC2 were more abundant in BC, while MSC3 was less prevalent (Figure 8b). A bubble plot further emphasized the specificity of marker genes across the three subtypes (Figure 8c). Biomarker expression was then evaluated across these subtypes. In tumor samples, JUN exhibited significantly lower expression in MSC1 and MSC2 (P < 0.05), while NFKB1 was significantly reduced in MSC2 and MSC3 (P < 0.01) (Figures 8d, e). SP1, however, showed no significant expression differences across the three subtypes (P > 0.05) (Figure 8f).

Figure 8
Panel a shows a UMAP plot depicting three MSC subtypes: MSC1 in blue, MSC2 in yellow, MSC3 in red. Panel b presents UMAP plots comparing normal and tumor samples with the same subtype color coding. Panel c displays a dot plot of gene expression levels across MSC subtypes with blue indicating expression intensity and dot size reflecting the percentage expressed. Panel d-f feature violin plots comparing JUN, NFKB1, and SP1 gene expressions in normal versus tumor samples for each MSC subtype, highlighting variations in distribution and expression levels.

Figure 8. Heterogeneity of mesenchymal cells. (a) Mesenchymal cells can be further divided into three cell subtypes. (b) The expression of these subtypes in tumor and normal tissues. (c) The specificity of marker genes across these 3 subtypes. (d-f) The expression levels of JUN (d), NFKB1 (e) and SP1 (f) in these subtypes. *:p<0.05; **:p<0.01; ***:p<0.001; ns:p>0.05

3.9 Analyzing communication relationships between mesenchymal cell subtypes and other cell types

The study further assessed the distribution and abundance of three mesenchymal cell subtypes, alongside eight other cell types, in both tumor and normal tissues, revealing significant variations in cell type proportions between these environments (Figure 9A). These differences suggested that BC may influence the interaction dynamics among various cell populations. Subsequently, communication analysis was conducted to evaluate interactions between MSC1, MSC2, and MSC3, as well as the other eight cell types, under both normal and tumor conditions, focusing on both the frequency and strength of these interactions. Under BC conditions, the number of interactions between the mesenchymal subtypes and other cell types was generally reduced compared to normal tissues (Figure 9B), implying that the tumor microenvironment may alter cell functionality and behavior by diminishing communication between mesenchymal cells and other cell types.

Figure 9
Two panels illustrate cell group data. Panel A shows bar charts with cell types on the y-axis and fractions and numbers on the x-axis, comparing normal and tumor groups. Panel B displays a circular network diagram of interactions among cell types, with lines indicating interaction numbers. Panel C features a similar network diagram, focusing on interaction weight/strength.

Figure 9. Communication relationships between mesenchymal cell subtypes and other cell types. (A) The distribution and quantity of different cell types in tumor and normal tissues. (B) The interaction numbers between these cell types in tumor and normal tissues. (C) The weight of interactions between these cell types in tumor and normal tissues.

In terms of interaction strength, compared to normal conditions, communication between MSC1 and endothelial cells was weakened in BC, while interactions with B cells and myeloid cells were enhanced. Similarly, MSC2 showed reduced interaction with endothelial cells but increased communication with basal and myeloid cells. MSC3 exhibited similar patterns, with decreased communication with endothelial cells and increased interaction with myeloid cells (Figure 9C). These findings may reflect alterations in immune regulation within the tumor microenvironment, potentially contributing to tumor invasiveness and environmental changes. The results suggest that BC could reorganize functional and interaction patterns among cells, particularly between mesenchymal and immune cells, which may be critical for understanding tumor progression and developing effective therapeutic strategies. Additionally, ligand-receptor interactions among the mesenchymal cell subtypes and other cell types under both normal and BC conditions were also examined (Figure 10), highlighting the importance of these interactions in regulating cellular behavior, maintaining tissue function, and responding to environmental changes.

Figure 10
Scatter plot displaying gene expression levels for various genes and conditions. Each subplot represents a dataset, with color indicating expression value intensity from blue (low) to red (high). Rows are labeled with gene names, and columns with experimental conditions. Legends are provided for color mapping and expression value categories.

Figure 10. Ligand-receptor interactions among these cell types under control and BC conditions.

3.10 Pseudo-time trajectory inference of mesenchymal cells

In the final stages of the study, pseudo-time analysis was performed on mesenchymal cells to trace their developmental trajectory. The analysis identified a clear starting point, from which cells underwent maturation as they progressed along the trajectory. Four differentiation nodes were detected, signaling points where individual cell clusters began to diverge into distinct cellular states or fates, underscoring the heterogeneity within mesenchymal cells that guides their progression through various developmental paths (Figure 11A). The mesenchymal cell trajectory was divided into nine stages, with cells evenly distributed across early, middle, and late developmental stages, highlighting the continuous and complex nature of mesenchymal cell development (Figure 11B).

Figure 11
Graphical representation of cell data shown in five panels (A-E).   A: Scatter plot with color gradient indicating pseudotime phases from zero to twenty on axes labeled Component 1 and Component 2.   B: Scatter plot with cells in nine different states, color-coded, on the same axes.   C: Dual scatter plots comparing normal cells and tumor cells, both color-coded by state, on similar axes.  D: Scatter plot distinguishing normal and tumor groups, color-coded.   E: Scatter plot categorizing cell types MSC1, MSC2, and MSC3, color-coded.   Indicated pathways demonstrate cell progression.

Figure 11. Pseudo-time trajectory inference of mesenchymal cells. (A) The trajectory plot illustrated the differentiation of mesenchymal cells over time. During the differentiation process, mesenchymal cells were divided into four distinct stages. Darker colors represented earlier time points in the differentiation timeline. (B) The trajectory was categorized into 9 stages. (C, D) The trajectory of mesenchymal cells in normal and tumor tissues respectively (C) and synthetically (D). (E) The distribution of 3 mesenchymal cell subtypes along the differentiation trajectory.

Notably, under BC conditions, the number of mesenchymal cells across all developmental stages significantly increased, corroborating previous findings (Figures 11C, D). This suggests that the BC microenvironment promotes an expansion of mesenchymal cells, enhancing their differentiation potential and possibly contributing to tumor progression and the remodeling of the tumor microenvironment. The study also depicted the distribution of the three mesenchymal subtypes—MSC1, MSC2, and MSC3—along the mesenchymal cell differentiation trajectory (Figure 11E). These subtypes were present throughout all stages, indicating the intrinsic heterogeneity and complex interactions among mesenchymal cells within the tissue. Interestingly, MSC1 was predominantly more abundant at every stage, possibly reflecting its central role in mesenchymal cell development and function. This suggests that MSC1 may possess greater biological activity or functionality, influencing both neighboring cells and the broader cellular community.

3.11 Verification of biomarkers expression

Previous analyses revealed significantly lower expression levels of JUN, NFKB1, and SP1 in BC samples from both the whole transcriptome RNA-seq dataset and the TCGA-BC dataset (P < 0.01) (Figures 3E, F). To further validate these findings at the clinical level, RT-qPCR was employed to assess the expression of these biomarkers in patients with BC (Figures 12A–C). RT-qPCR results confirmed that both NFKB1 and SP1 exhibited significantly reduced expression in BC samples (P < 0.05), aligning with initial observations. Although JUN demonstrated a trend toward decreased expression in BC, the difference was not statistically significant (P = 0.0139).

Figure 12
Bar graphs show the relative expression levels of JUN, NFKB1, and SP1 normalized to GAPDH in normal and tumor samples. Graph A (JUN) shows no significant difference (p=0.0813). Graph B (NFKB1) shows significant downregulation in tumors (p=0.0139). Graph C (SP1) shows significant downregulation in tumors (p=0.0058). Error bars indicate variability.

Figure 12. Expression level of biomarkers’ mRNA. (A-C) Expression levels of JUN (A), NFKB1 (B), and SP1 (C) in patients with BC by RT-qPCR. ns represents P > 0.05, * represents P < 0.05, ** represents P < 0.01. ns means 'not significant', and it represents p>0.05.

4 Discussion

BC is a multifactorial disease, influenced by various biological factors, with SCRGs and TMMRGs playing critical roles in tumor progression through the enhancement of tumor cell proliferation and self-renewal capabilities (23). Understanding the interactions between these mechanisms is crucial for identifying potential therapeutic targets and improving patient outcomes. This study highlights JUN, NFKB1, and SP1 as biomarkers with significantly reduced expression levels in BC tissues. These biomarkers are implicated in the progression of BC through various biological processes, such as ECM remodeling, cell signaling regulation, oxidative stress response, and mechanisms governing protein synthesis and transport. Notably, scRNA-seq analysis revealed high infiltration levels of activated B cells and NK cells in BC, alongside elevated expression of immune checkpoints, including CD200, CD274, TIGIT, TNFRSF25, and TNFSF15.

JUN promotes the self-renewal and stemness maintenance of cancer stem cells by activating stem cell-related pathways such as β-catenin and Notch and enhances tumor initiation ability and chemotherapy resistance (54, 55). Abnormal activation of telomerase reverse transcriptase (TERT) in cancer cells enables them to escape cellular senescence and thereby acquire unlimited proliferation ability. JUN may indirectly promote the expression of TERT, activate telomerase activity, compensate for telomere shortening during cell division, maintain chromosomal stability, and assist in the unlimited proliferation of cancer cells (56, 57). NFKB1 encoded a key component of the NF-κB transcription factor complex, which was critical in the occurrence and development of various diseases including cancer (58). Moreover, NFKB1 is involved in regulating the self-renewal, differentiation and treatment resistance of cancer stem cells (59). In addition, studies have shown that NFKB1 has a subunit specific role in the DNA damage response (60). This suggests that NFKB1 may indirectly affect the stability of telomeres by influencing the DNA damage repair pathway. SP1 has also been found to be related to the characteristics of cancer stem cells. For example, in triple-negative breast cancer (TNBC), SGCE promotes BC stemness by facilitating the transcription of FGF-BP1 by SP1 (61). Furthermore, Liu et al. found that SP1 mediated the inhibition of TERT after the overexpression of LKB1, thereby affecting the progression of lung adenocarcinoma. That is to say, SP1 might also affect TERT through a similar mechanism in BC and thereby influence the maintenance mechanism of telomeres. In conclusion, JUN, NFKB1 and SP1 may jointly maintain the characteristics of cancer stem cells by regulating the stemness maintenance and telomerase activity of cancer stem cells, thus influencing the progression, treatment resistance and recurrence of BC.

Moreover, JUN, SP1 and NFKB1 play multiple roles in the regulation of the immune microenvironment of BC. JUN, as a kernel factor of the AP-1 family, promotes the invasion of tumor cells by activating Matrix Metalloproteinases and at the same time down-regulates the expression of MHC-I class molecules to evade the recognition of cytotoxic T cells; the secretion of IL-6/IL-8 induced by it can also recruit immunosuppressive myeloid cells and inhibit T cell function (6264). SP1 regulates the GC-rich region of the PD-L1 promoter, enhances its transcriptional level, and mediates tumor immune escape (65, 66). Meanwhile, it can promote the expression of chemokine CCL2, recruit immunosuppressive macrophages, and affect the immune microenvironment (67). NFKB1 can continuously activate the NF-κB pathway. The CXCL12/CXCR4 axis NFBK1 regulates can also promote the infiltration of regulatory T cells (Tregs) and inhibit anti-tumor immunity (6870). They three jointly reshape the immune microenvironment through the transcriptional regulatory network and weaken the body’s anti-tumor immune response, which may become a potential mechanism for drug resistance to immunotherapy in breast cancer.

Immune infiltration analysis revealed a significant increase in the abundance of activated B cells and NK cells in BC tumor samples compared to normal tissue (P < 0.05). This elevated infiltration suggests that the immune system is activated within the BC microenvironment, potentially reflecting the body’s attempt to combat tumor progression. B cells are important components of adaptive immunity and can mediate the tumor cells killing by generating antibodies (71). B cells can also act as antigen-presenting cells (APC), presenting tumor antigens to T cells and activating the anti-tumor immune response of T cells (72). Multiple studies have shown that B-cell infiltration is associated with a favorable prognosis in cancer patients (73, 74). A higher level of activated B cell infiltration may indicate a stronger anti-tumor immune response and better clinical outcomes. NK cells are important components of the innate immune system, capable of recognizing and killing tumor cells without prior sensitization (7576). Moreover, studies have shown that activated NK cells may play a significant role in early tumor control (77), thereby improving the survival rate of patients. In conclusion, B cells and NK cells exert immune effects respectively by mediating antibody killing, antigen presentation and direct killing of tumor cells. Their activation degree and the infiltration level are significantly associated with the prognosis of cancer patients.

Furthermore, significant differences in the expression of immune checkpoints—specifically CD200, CD274 (PD-L1), TIGIT, TNFRSF25, and TNFSF15—were observed between tumor and normal samples. All five immune checkpoints showed elevated expression in tumor samples (P < 0.05), suggesting their activation within the tumor microenvironment, which likely aids tumor cells in evading immune surveillance. CD200 functions as an immunosuppressive molecule through its receptor CD200R, and its upregulation in BC may facilitate immune escape and contribute to tumor progression and metastasis. PD-L1, a key immune checkpoint, interacts with PD-1 on T cells when expressed on tumor cells, leading to T cell inhibition. High PD-L1 levels in BC are associated with poor prognosis, and blocking PD-L1/PD-1 interactions has emerged as an effective treatment for various BC subtypes. The elevated expression of these immune checkpoints highlights their potential as therapeutic targets in BC. Inhibiting these checkpoints could restore immune surveillance and enhance the efficacy of immune-based therapies (78). This immune analysis highlights the complex interplay between immune cells and the tumor microenvironment in BC. The increased infiltration of activated B cells and NK cells suggests an active immune response, while the upregulation of specific immune checkpoints points to mechanisms of immune evasion. Understanding these dynamics is essential for developing novel immunotherapeutic strategies aimed at boosting anti-tumor immunity while counteracting immune suppression. Future research should focus on elucidating the roles of these immune components and their interactions in BC, as well as exploring combination therapies that could synergistically improve treatment outcomes.

This study explored the lncRNA-miRNA-mRNA regulatory network in BC, revealing interactions between non-coding RNAs and key biomarkers (79, 80). For instance, XIST was predicted to regulate JUN through hsa-miR-139-5p, and may also influence SP1 via several miRNAs, including hsa-miR-199a-3p, hsa-miR-7-5p, hsa-miR-324-5p, and hsa-miR-299-5p. These interactions suggest complex regulatory mechanisms in which non-coding RNAs play pivotal roles in modulating the expression of these biomarkers. Given that this network was constructed based on differential expression analysis and database predictions, the reliability of these regulatory relationships is notably high. Previous research has demonstrated that the loss of OPA interacting protein 5 can inhibit BC proliferation via the hsa-miR-139-5p/NOTCH1 axis (81), and hsa-miR-7-5p has also been implicated in BC progression (82). Our findings further support the potential of targeting these non-coding RNA interactions for therapeutic intervention in BC. Specifically, the results suggest that XIST and the identified miRNAs (e.g., hsa-miR-139-5p and hsa-miR-7-5p) could jointly regulate SP1 and JUN, thereby influencing key cellular processes such as proliferation, apoptosis, and metastasis, which contribute to BC development. This highlights the potential co-regulatory roles of non-coding RNAs within the tumor microenvironment and their importance in BC progression.

Additionally, the study identified several drug candidates targeting these biomarkers, with Bardoxolone (83), Andrographolide (84), and Terameprocol (85, 86) emerging as potentially effective agents in BC treatment. These drugs target critical pathways involved in oxidative stress, apoptosis, and transcriptional regulation, positioning them as promising candidates for further investigation. However, additional in vivo and in vitro studies, along with large-scale clinical trials, are necessary to validate their efficacy and safety in BC.

This study further identified three mesenchymal cell subtypes—MSC1, MSC2, and MSC3—each displaying distinct expression patterns of BC biomarkers. JUN was significantly downregulated in MSC1 and MSC2, while NFKB1 was downregulated in MSC2 and MSC3. However, SP1 did not show significant expression differences across these subtypes. The downregulation of JUN and NFKB1 in mesenchymal subtypes may impair their ability to maintain tumor-suppressive functions, highlighting their potential role in driving BC progression. In addition, studies have shown that the MSC1 subtype may have pro-inflammatory properties and play a specific role in immune regulation (87). For example, MSCS activated by Toll-like receptor 4 ligand lipopolysaccharide may present the MSC1 phenotype and express pro-inflammatory cytokines such as IL-6 and IL-8 (87, 88). The MSC2 subtype may have immunosuppressive and tissue repair functions (87, 89). It is mainly produced in an anti-inflammatory environment and is characterized by high expression of immunosuppressive molecules such as IDO, PGE2 and TGF-β, which can inhibit the activity of immune cells and alleviate the inflammatory response (87, 89, 90). However, there are relatively few studies on the MSC3 subtype, and its specific functions remain to be further clarified. The heterogeneity of these mesenchymal subtypes was further explored through communication and pseudo-time analysis. Notably, MSC1 and MSC2 were more abundant in tumor tissues, suggesting that these subtypes may play an active role in tumor progression. The reduced interactions between these mesenchymal subtypes, along with enhanced communication with immune cells such as B cells and myeloid cells under BC conditions, suggest that the tumor microenvironment alters mesenchymal cell functionality. These changes may contribute to immune suppression, tumor cell survival, and sustained tumor growth (91, 92).

Pseudo-time analysis also revealed that mesenchymal cells undergo continuous differentiation across nine distinct stages, emphasizing their developmental plasticity. This trajectory analysis indicated that BC significantly increases the number of mesenchymal cells at all stages of development, further underscoring their critical role in tumor development. MSC1, in particular, was the most prominent subtype across all stages, potentially due to its central role in mesenchymal cell function and its influence on surrounding cell types within the tumor microenvironment (92, 93). Overall, mesenchymal cells appear to contribute to several aspects of BC pathogenesis, including the establishment of the tumor microenvironment, invasion, metastasis, and resistance to treatment. Their ability to interact with immune cells and modulate the ECM suggests that targeting these cells could disrupt cancer-supportive processes. Additionally, targeting epithelial-mesenchymal transition (EMT) and mesenchymal-epithelial transition (MET) pathways in mesenchymal cells may offer new therapeutic strategies to inhibit tumor metastasis and progression in BC.

In summary, this study identified JUN, NFKB1, and SP1 as important biomarkers associated with stem cells and telomere maintenance mechanisms in BC, highlighting their critical roles in tumor progression, alterations in the immune microenvironment, and potential therapeutic targets. Additionally, the exploration of mesenchymal cell heterogeneity and their dynamic communication with other cell types provides valuable insights into the complexity of the tumor microenvironment. However, the study has certain limitations, particularly the small sample size used for validation, which may restrict the generalizability of the results. To address this, we will actively expand the sample sources, collaborate with multiple centers, and incorporate more clinical samples for analysis and validation in the future. Meanwhile, future research should focus on integrating bioinformatic predictions with experimental validation to establish direct evidence for the involvement of JUN, SP1, and NFKB1 in telomere maintenance and stemness regulation.

5 Conclusions

In conclusion, this study provides valuable insights into the roles of mesenchymal and luminal cells in BC, while also underscoring the need for further research to validate and expand these findings. Addressing the limitations of the current study will be critical for enhancing our understanding of the underlying mechanisms and improving treatment outcomes for patients with BC.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics statement

The studies involving humans were approved by Ethics Committee of The First Affiliated Hospital with Nanjing Medical University. The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study.

Author contributions

SZ: Conceptualization, Data curation, Validation, Visualization, Writing – original draft, Writing – review & editing. LH: Data curation, Validation, Visualization, Writing – review & editing. QZ: Validation, Writing – review & editing. RC: Visualization, Writing – review & editing. JW: Conceptualization, Supervision, Writing – review & editing. XZ: Conceptualization, Project administration, Supervision, Writing – review & editing.

Funding

The studies involving humans were approved by Ethics Committee of The First Affiliated Hospital with Nanjing Medical University, and the approval number and date of approval are as follows: 2023-SR-760 and 2023-11-22. The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study.

Acknowledgments

We would like to express our sincere gratitude to all individuals and organizations who supported and assisted us throughout this research. Special thanks to the following authors: Shuang Zhen, Lifeng Huang, Qiannan Zhu, Rui Chen, Jue Wang, Xiaoming Zha. In conclusion, we extend our thanks to everyone who has supported and assisted us along the way. Without your support, this research would not have been possible.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2025.1618193/full#supplementary-material

Glossary

SCRGs: Stem Cell-Related Genes

TMMRGs: Telomere Maintenance Mechanism-Related Genes

BC: Breast Cancer

RNA-seq: RNA Sequencing

scRNA-seq: Single-Cell RNA Sequencing

DEGs: Differentially Expressed Genes

PPI: Protein-Protein Interaction

MCODE: Molecular Complex Detection

LASSO: Least Absolute Shrinkage and Selection Operator

GO: Gene Ontology

KEGG: Kyoto Encyclopedia of Genes and Genomes

GSEA: Gene Set Enrichment Analysis

ssGSEA: Single-Sample Gene Set Enrichment Analysis

TCGA: The Cancer Genome Atlas

GEO: Gene Expression Omnibus

RT-qPCR: Reverse Transcription Quantitative Polymerase Chain Reaction

ECM: Extracellular Matrix

NK cells: Natural Killer cells

MSC: Mesenchymal Stem Cells

EMT: Epithelial-Mesenchymal Transition

MET: Mesenchymal-Epithelial Transition

PCA: Principal Component Analysis

UMAP: Uniform Manifold Approximation and Projection

FPKM: Fragments Per Kilobase of transcript per Million mapped reads

BP: Biological Process

CC: Cellular Component

MF: Molecular Function

lncRNA: Long Non-Coding RNA

miRNA: MicroRNA

mRNA: Messenger RNA

DGIdb: Drug Gene Interaction Database

TLS: Tertiary Lymphoid Structures

IFN: Interferon

TNF: Tumor Necrosis Factor

PD-L1: Programmed Death-Ligand 1

PD-1: Programmed Cell Death Protein 1

TIGIT: T cell Immunoreceptor with Ig and ITIM domains

TNFRSF25: Tumor Necrosis Factor Receptor Superfamily Member 25

TNFSF15: Tumor Necrosis Factor Superfamily Member 15

CD200: Cluster of Differentiation 200

CD274: Cluster of Differentiation 274

Nrf2: Nuclear factor erythroid 2–related factor 2

RhoA: Ras homolog family member A

S1P: Sphingosine-1-phosphate

SRP: Signal Recognition Particle

NF-κB: Nuclear Factor kappa-light-chain-enhancer of activated B cells

STAT3: Signal Transducer and Activator of Transcription 3

TP53: Tumor Protein p53

TMM: Telomere maintenance mechanisms

TERT: Telomerase reverse transcriptase

TNBC: Triple-negative breast cancer

MDSCs: Myeloid-derived suppressor cells

HCC: Hepatocellular carcinoma

MHC-1: Major Histocompatibility Complex Class I

References

1. Xie J, Zou Y, Gao T, Xie L, Tan D, and Xie X. Therapeutic landscape of human epidermal growth factor receptor 2-positive breast cancer. Cancer Control. (2022) 29:10732748221099230. doi: 10.1177/10732748221099230

PubMed Abstract | Crossref Full Text | Google Scholar

2. Hosio M, Urpilainen E, Marttila M, Hautakoski A, Arffman M, Sund R, et al. Association of antidiabetic medication and statins with breast cancer incidence in women with type 2 diabetes. Breast Cancer Res Treat. (2019) 175:741–8. doi: 10.1007/s10549-019-05185-0

PubMed Abstract | Crossref Full Text | Google Scholar

3. Mokhtari E, Jamshidi S, Daftari G, Farhadnejad H, Teymoori F, Momeni SA, et al. The relationship between the insulinemic potential of diet and lifestyle and risk of breast cancer: a case-control study among Iranian adult women. Arch Public Health. (2023) 81:4. doi: 10.1186/s13690-022-01016-9

PubMed Abstract | Crossref Full Text | Google Scholar

4. Shao X, Xie N, Chen Z, Wang X, Cao W, Zheng Y, et al. Inetetamab for injection in combination with vinorelbine weekly or every three weeks in HER2-positive metastatic breast cancer: A multicenter, randomized, phase II clinical trial. J Transl Int Med. (2024) 12:466–77. doi: 10.1515/jtim-2024-0022

PubMed Abstract | Crossref Full Text | Google Scholar

5. Katsura C, Ogunmwonyi I, Kankam HK, and Saha S. Breast cancer: presentation, investigation and management. Br J Hosp Med (London England: 2005). (2022) 83:1–7. doi: 10.12968/hmed.2021.0459

PubMed Abstract | Crossref Full Text | Google Scholar

6. El Masri J and Phadke S. Breast cancer epidemiology and contemporary breast cancer care: A review of the literature and clinical applications. Clin Obstet Gynecol. (2022) 65:461–81. doi: 10.1097/GRF.0000000000000721

PubMed Abstract | Crossref Full Text | Google Scholar

7. Mokhtari-Hessari P and Montazeri A. Health-related quality of life in breast cancer patients: review of reviews from 2008 to 2018. Health Qual Life Outcomes. (2020) 18:338. doi: 10.1186/s12955-020-01591-x

PubMed Abstract | Crossref Full Text | Google Scholar

8. Robinson NJ, Morrison-Smith CD, Gooding AJ, Schiemann BJ, Jackson MW, Taylor DJ, et al. SLX4IP and telomere dynamics dictate breast cancer metastasis and therapeutic responsiveness. Life Sci Alliance. (2020) 3. doi: 10.26508/lsa.201900427

PubMed Abstract | Crossref Full Text | Google Scholar

9. Jose SS, Bendickova K, Kepak T, Krenova Z, and Fric J. Chronic inflammation in immune aging: role of pattern recognition receptor crosstalk with the telomere complex? Front Immunol. (2017) 8:1078. doi: 10.3389/fimmu.2017.01078

PubMed Abstract | Crossref Full Text | Google Scholar

10. Ji Y, Dang X, Nguyen LNT, Nguyen LN, Zhao J, Cao D, et al. Topological DNA damage, telomere attrition and T cell senescence during chronic viral infections. Immun Ageing. (2019) 16:12. doi: 10.1186/s12979-019-0153-z

PubMed Abstract | Crossref Full Text | Google Scholar

11. Cervantes GIV, Katz J, Lopez-Rosas A, Han Y, Lee-Chang C, Miska J, et al. IMMU-06. THE ROLE OF BRAIN TUMOR MICROENVIRONMENT ON T CELL TELOMERE SHORTENING. Neuro-Oncology. (2023) 25:v142–v. doi: 10.1093/neuonc/noad179.0538

Crossref Full Text | Google Scholar

12. Yuan X, Yuan H, Zhang N, Liu T, and Xu D. Thyroid carcinoma-featured telomerase activation and telomere maintenance: Biology and translational/clinical significance. Clin Trans Med. (2022) 12:e1111. doi: 10.1002/ctm2.1111

PubMed Abstract | Crossref Full Text | Google Scholar

13. Zhang Y, Luo S, Jia Y, and Zhang X. Telomere maintenance mechanism dysregulation serves as an early predictor of adjuvant therapy response and a potential therapeutic target in human cancers. Int J Cancer. (2022) 151:313–27. doi: 10.1002/ijc.34007

PubMed Abstract | Crossref Full Text | Google Scholar

14. Lv J, Zhao X, Zhao L, Gong C, Zheng W, Guo L, et al. The notable role of telomere length maintenance in complex diseases. Biomedicines. (2024) 12:2611. doi: 10.3390/biomedicines12112611

PubMed Abstract | Crossref Full Text | Google Scholar

15. Werr L, Rosswog C, Bartenhagen C, George SL, and Fischer M. Telomere maintenance mechanisms in neuroblastoma: New insights and translational implications. EJC Paediatr Oncol. (2024) 3. doi: 10.1016/j.ejcped.2024.100156

Crossref Full Text | Google Scholar

16. Xu Z, Wang H, Tian H, Wang W, and Hua D. Identification of telomere maintenance-driven molecular subtypes in hepatocellular carcinoma: implications for prognosis and targeted therapy via KPNA2. Discov Oncol. (2025) 16:516. doi: 10.1007/s12672-025-02311-x

PubMed Abstract | Crossref Full Text | Google Scholar

17. Butti R, Gunasekaran VP, Kumar TVS, Banerjee P, and Kundu GC. Breast cancer stem cells: Biology and therapeutic implications. Int J Biochem Cell Biol. (2019) 107:38–52. doi: 10.1016/j.biocel.2018.12.001

PubMed Abstract | Crossref Full Text | Google Scholar

18. Yang F, Xu J, Tang L, and Guan X. Breast cancer stem cell: the roles and therapeutic implications. Cell Mol Life Sci: CMLS. (2017) 74:951–66. doi: 10.1007/s00018-016-2334-7

PubMed Abstract | Crossref Full Text | Google Scholar

19. Tsuchiya H and Shiota G. Immune evasion by cancer stem cells. Regener Ther. (2021) 17:20–33. doi: 10.1016/j.reth.2021.02.006

PubMed Abstract | Crossref Full Text | Google Scholar

20. An J, Hu X, and Liu F. Current understanding of cancer stem cells: Immune evasion and targeted immunotherapy in gastrointestinal Malignancies. Front Oncol. (2023) 13 - 2023. doi: 10.3389/fonc.2023.1114621

PubMed Abstract | Crossref Full Text | Google Scholar

21. Chakraborty S, Mukherjee S, Basak U, Pati S, Dutta A, Dutta S, et al. Immune evasion by cancer stem cells ensures tumor initiation and failure of immunotherapy. Explor Immunol. (2023) 3:384–405. doi: 10.37349/ei

Crossref Full Text | Google Scholar

22. Agudo J and Miao Y. Stemness in solid Malignancies: coping with immune attack. Nat Rev Cancer. (2025) 25:27–40. doi: 10.1038/s41568-024-00760-0

PubMed Abstract | Crossref Full Text | Google Scholar

23. Robinson NJ, Taylor DJ, and Schiemann WP. Stem cells, immortality, and the evolution of metastatic properties in breast cancer: telomere maintenance mechanisms and metastatic evolution. J Cancer Metastasis Treat. (2019) 5:39. doi: 10.20517/2394-4722.2019.15

PubMed Abstract | Crossref Full Text | Google Scholar

24. Quaglino E, Conti L, and Cavallo F. Breast cancer stem cell antigens as targets for immunotherapy. Semin Immunol. (2020) 47:101386. doi: 10.1016/j.smim.2020.101386

PubMed Abstract | Crossref Full Text | Google Scholar

25. Zhang L, Chen W, Liu S, and Chen C. Targeting breast cancer stem cells. Int J Biol Sci. (2023) 19:552–70. doi: 10.7150/ijbs.76187

PubMed Abstract | Crossref Full Text | Google Scholar

26. Bai X, Ni J, Beretov J, Graham P, and Li Y. Cancer stem cell in breast cancer therapeutic resistance. Cancer Treat Rev. (2018) 69:152–63. doi: 10.1016/j.ctrv.2018.07.004

PubMed Abstract | Crossref Full Text | Google Scholar

27. Taurin S and Alkhalifa H. Breast cancers, mammary stem cells, and cancer stem cells, characteristics, and hypotheses. Neoplasia (New York NY). (2020) 22:663–78. doi: 10.1016/j.neo.2020.09.009

PubMed Abstract | Crossref Full Text | Google Scholar

28. Uprety B and Abrahamse H. Targeting breast cancer and their stem cell population through AMPK activation: novel insights. Cells. (2022) 11:576. doi: 10.3390/cells11030576

PubMed Abstract | Crossref Full Text | Google Scholar

29. Ren X, Zhang L, Zhang Y, Li Z, Siemers N, and Zhang Z. Insights gained from single-cell analysis of immune cells in the tumor microenvironment. Annu Rev Immunol. (2021) 39:583–609. doi: 10.1146/annurev-immunol-110519-071134

PubMed Abstract | Crossref Full Text | Google Scholar

30. Zhang Y, Wang D, Peng M, Tang L, Ouyang J, Xiong F, et al. Single-cell RNA sequencing in cancer research. J Exp Clin Cancer Res. (2021) 40:81. doi: 10.1186/s13046-021-01874-1

PubMed Abstract | Crossref Full Text | Google Scholar

31. Liu C, Xie J, Lin B, Tian W, Wu Y, Xin S, et al. Pan-cancer single-cell and spatial-resolved profiling reveals the immunosuppressive role of APOE+ Macrophages in immune checkpoint inhibitor therapy. Adv Sci (Weinh). (2024) 11:e2401061. doi: 10.1002/advs.202401061

PubMed Abstract | Crossref Full Text | Google Scholar

32. Xie J, Yang A, Liu Q, Deng X, Lv G, Ou X, et al. Single-cell RNA sequencing elucidated the landscape of breast cancer brain metastases and identified ILF2 as a potential therapeutic target. Cell Prolif. (2024) 57:e13697. doi: 10.1111/cpr.13697

PubMed Abstract | Crossref Full Text | Google Scholar

33. Xie J, Deng W, Deng X, Liang JY, Tang Y, Huang J, et al. Single-cell histone chaperones patterns guide intercellular communication of tumor microenvironment that contribute to breast cancer metastases. Cancer Cell Int. (2023) 23:311. doi: 10.1186/s12935-023-03166-4

PubMed Abstract | Crossref Full Text | Google Scholar

34. Kim H, Whitman AA, Wisniewska K, Kakati RT, Garcia-Recio S, Calhoun BC, et al. Tamoxifen response at single-cell resolution in estrogen receptor-positive primary human breast tumors. Clin Cancer Res. (2023) 29:4894–907. doi: 10.1158/1078-0432.CCR-23-1248

PubMed Abstract | Crossref Full Text | Google Scholar

35. Zheng H, Liu H, Li H, Dou W, Wang J, Zhang J, et al. Characterization of stem cell landscape and identification of stemness-relevant prognostic gene signature to aid immunotherapy in colorectal cancer. Stem Cell Res Ther. (2022) 13:244. doi: 10.1186/s13287-022-02913-0

PubMed Abstract | Crossref Full Text | Google Scholar

36. Chen S. Ultrafast one-pass FASTQ data preprocessing, quality control, and deduplication using fastp. Imeta. (2023) 2:e107. doi: 10.1002/imt2.107

PubMed Abstract | Crossref Full Text | Google Scholar

37. Kim D, Paggi JM, Park C, Bennett C, and Salzberg SL. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat Biotechnol. (2019) 37:907–15. doi: 10.1038/s41587-019-0201-4

PubMed Abstract | Crossref Full Text | Google Scholar

38. Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, and Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. (2015) 33:290–5. doi: 10.1038/nbt.3122

PubMed Abstract | Crossref Full Text | Google Scholar

39. Love MI, Huber W, and Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. (2014) 15:550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | Crossref Full Text | Google Scholar

40. Gustavsson EK, Zhang D, Reynolds RH, Garcia-Ruiz S, and Ryten M. ggtranscript: an R package for the visualization and interpretation of transcript isoforms using ggplot2. Bioinformatics. (2022) 38:3844–6. doi: 10.1093/bioinformatics/btac409

PubMed Abstract | Crossref Full Text | Google Scholar

41. Gu Z and Hübschmann D. Make interactive complex heatmaps in R. Bioinformatics. (2022) 38:1460–2. doi: 10.1093/bioinformatics/btab806

PubMed Abstract | Crossref Full Text | Google Scholar

42. Yu G, Wang LG, Han Y, and He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics. (2012) 16:284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | Crossref Full Text | Google Scholar

43. Chen H and Boutros PC. VennDiagram: a package for the generation of highly-customizable Venn and Euler diagrams in R. BMC Bioinf. (2011) 12:35. doi: 10.1186/1471-2105-12-35

PubMed Abstract | Crossref Full Text | Google Scholar

44. Liu P, Xu H, Shi Y, Deng L, and Chen X. Potential molecular mechanisms of plantain in the treatment of gout and hyperuricemia based on network pharmacology. Evid Based Complement Alternat Med. (2020) 2020:3023127. doi: 10.1155/2020/3023127

PubMed Abstract | Crossref Full Text | Google Scholar

45. Li Y, Lu F, and Yin Y. Applying logistic LASSO regression for the diagnosis of atypical Crohn’s disease. Sci Rep. (2022) 12:11340. doi: 10.1038/s41598-022-15609-5

PubMed Abstract | Crossref Full Text | Google Scholar

46. Wang L, Wang D, Yang L, Zeng X, Zhang Q, Liu G, et al. Cuproptosis related genes associated with Jab1 shapes tumor microenvironment and pharmacological profile in nasopharyngeal carcinoma. Front Immunol. (2022) 13:989286. doi: 10.3389/fimmu.2022.989286

PubMed Abstract | Crossref Full Text | Google Scholar

47. Hänzelmann S, Castelo R, and Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinf. (2013) 14:7. doi: 10.1186/1471-2105-14-7

PubMed Abstract | Crossref Full Text | Google Scholar

48. Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, et al. Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep. (2017) 18:248–62. doi: 10.1016/j.celrep.2016.12.019

PubMed Abstract | Crossref Full Text | Google Scholar

49. Pei S, Zhang P, Yang L, Kang Y, Chen H, Zhao S, et al. Exploring the role of sphingolipid-related genes in clinical outcomes of breast cancer. Front Immunol. (2023) 14:1116839. doi: 10.3389/fimmu.2023.1116839

PubMed Abstract | Crossref Full Text | Google Scholar

50. Satija R, Farrell JA, Gennert D, Schier AF, and Regev A. Spatial reconstruction of single-cell gene expression data. Nat Biotechnol. (2015) 33:495–502. doi: 10.1038/nbt.3192

PubMed Abstract | Crossref Full Text | Google Scholar

51. Jin S, Guerrero-Juarez CF, Zhang L, Chang I, Ramos R, Kuan CH, et al. Inference and analysis of cell-cell communication using CellChat. Nat Commun. (2021) 12:1088. doi: 10.1038/s41467-021-21246-9

PubMed Abstract | Crossref Full Text | Google Scholar

52. Du J, Yuan X, Deng H, Huang R, Liu B, Xiong T, et al. Single-cell and spatial heterogeneity landscapes of mature epicardial cells. J Pharm Anal. (2023) 13:894–907. doi: 10.1016/j.jpha.2023.07.011

PubMed Abstract | Crossref Full Text | Google Scholar

53. Livak KJ and Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods. (2001) 25:402–8. doi: 10.1006/meth.2001.1262

PubMed Abstract | Crossref Full Text | Google Scholar

54. Xie X, Kaoud TS, Edupuganti R, Zhang T, Kogawa T, Zhao Y, et al. c-Jun N-terminal kinase promotes stem cell phenotype in triple-negative breast cancer through upregulation of Notch1 via activation of c-Jun. Oncogene. (2017) 36:2599–608. doi: 10.1038/onc.2016.417

PubMed Abstract | Crossref Full Text | Google Scholar

55. Zou Y, Lin X, Bu J, Lin Z, Chen Y, Qiu Y, et al. Timeless-Stimulated miR-5188-FOXO1/β-Catenin-c-Jun Feedback Loop Promotes Stemness via Ubiquitination of β-Catenin in Breast Cancer. Mol Ther. (2020) 28:313–27. doi: 10.1016/j.ymthe.2019.08.015

PubMed Abstract | Crossref Full Text | Google Scholar

56. Liu M, Zhang Y, Jian Y, Gu L, Zhang D, Zhou H, et al. The regulations of telomerase reverse transcriptase (TERT) in cancer. Cell Death Dis. (2024) 15:90. doi: 10.1038/s41419-024-06454-7

PubMed Abstract | Crossref Full Text | Google Scholar

57. Yang R, Han Y, Guan X, Hong Y, Meng J, Ding S, et al. Regulation and clinical potential of telomerase reverse transcriptase (TERT/hTERT) in breast cancer. Cell Commun Signal. (2023) 21:218. doi: 10.1186/s12964-023-01244-8

PubMed Abstract | Crossref Full Text | Google Scholar

58. Cartwright T, Perkins ND, and C LW. NFKB1: a suppressor of inflammation, ageing and cancer. FEBS J. (2016) 283:1812–22. doi: 10.1111/febs.13627

PubMed Abstract | Crossref Full Text | Google Scholar

59. Kaltschmidt C, Banz-Jansen C, Benhidjeb T, Beshay M, Förster C, Greiner J, et al. A role for NF-κB in organ specific cancer and cancer stem cells. Cancers (Basel). (2019) 11:655. doi: 10.3390/cancers11050655

PubMed Abstract | Crossref Full Text | Google Scholar

60. Voce DJ, Schmitt AM, Uppal A, McNerney ME, Bernal GM, Cahill KE, et al. Nfkb1 is a haploinsufficient DNA damage-specific tumor suppressor. Oncogene. (2015) 34:2807–13. doi: 10.1038/onc.2014.211

PubMed Abstract | Crossref Full Text | Google Scholar

61. Qiu T, Hou L, Zhao L, Wang X, Zhou Z, Yang C, et al. SGCE promotes breast cancer stemness by promoting the transcription of FGF-BP1 by Sp1. J Biol Chem. (2023) 299:105351. doi: 10.1016/j.jbc.2023.105351

PubMed Abstract | Crossref Full Text | Google Scholar

62. Li Y, Samuvel DJ, Sundararaj KP, Lopes-Virella MF, and Huang Y. IL-6 and high glucose synergistically upregulate MMP-1 expression by U937 mononuclear phagocytes via ERK1/2 and JNK pathways and c-Jun. J Cell Biochem. (2010) 110:248–59. doi: 10.1002/jcb.22532

PubMed Abstract | Crossref Full Text | Google Scholar

63. Dong G, Chen P, Xu Y, Liu T, and Yin R. Cancer-associated fibroblasts: Key criminals of tumor pre-metastatic niche. Cancer Lett. (2023) 566:216234. doi: 10.1016/j.canlet.2023.216234

PubMed Abstract | Crossref Full Text | Google Scholar

64. Lu J, Luo Y, Rao D, Wang T, Lei Z, Chen X, et al. Myeloid-derived suppressor cells in cancer: therapeutic targets to overcome tumor immune evasion. Exp Hematol Oncol. (2024) 13:39. doi: 10.1186/s40164-024-00505-7

PubMed Abstract | Crossref Full Text | Google Scholar

65. Tao LH, Zhou XR, Li FC, Chen Q, Meng FY, Mao Y, et al. A polymorphism in the promoter region of PD-L1 serves as a binding-site for SP1 and is associated with PD-L1 overexpression and increased occurrence of gastric cancer. Cancer Immunol Immunother. (2017) 66:309–18. doi: 10.1007/s00262-016-1936-0

PubMed Abstract | Crossref Full Text | Google Scholar

66. Yang K, Li J, Zhu J, Chen Y, He Y, Wang J, et al. HOOK3 suppresses proliferation and metastasis in gastric cancer via the SP1/VEGFA axis. Cell Death Discov. (2024) 10:33. doi: 10.1038/s41420-024-01808-8

PubMed Abstract | Crossref Full Text | Google Scholar

67. Speranza FJ, Mahankali M, and Gomez-Cambronero J. Macrophage migration arrest due to a winning balance of Rac2/Sp1 repression over β-catenin-induced PLD expression. J Leukoc Biol. (2013) 94:953–62. doi: 10.1189/jlb.0313174

PubMed Abstract | Crossref Full Text | Google Scholar

68. Mortezaee K. CXCL12/CXCR4 axis in the microenvironment of solid tumors: A critical mediator of metastasis. Life Sci. (2020) 249:117534. doi: 10.1016/j.lfs.2020.117534

PubMed Abstract | Crossref Full Text | Google Scholar

69. Bocchi M, de Sousa Pereira N, de Oliveira KB, and Amarante MK. Involvement of CXCL12/CXCR4 axis in colorectal cancer: a mini-review. Mol Biol Rep. (2023) 50:6233–9. doi: 10.1007/s11033-023-08479-1

PubMed Abstract | Crossref Full Text | Google Scholar

70. Zhi Y, Lu H, Duan Y, Sun W, Guan G, Dong Q, et al. Involvement of the nuclear factor-κB signaling pathway in the regulation of CXC chemokine receptor-4 expression in neuroblastoma cells induced by tumor necrosis factor-α. Int J Mol Med. (2015) 35:349–57. doi: 10.3892/ijmm.2014.2032

PubMed Abstract | Crossref Full Text | Google Scholar

71. Laumont CM, Banville AC, Gilardi M, Hollern DP, and Nelson BH. Tumour-infiltrating B cells: immunological mechanisms, clinical impact and therapeutic opportunities. Nat Rev Cancer. (2022) 22:414–30. doi: 10.1038/s41568-022-00466-1

PubMed Abstract | Crossref Full Text | Google Scholar

72. Rastogi I, Jeon D, Moseman JE, Muralidhar A, Potluri HK, and McNeel DG. Role of B cells as antigen presenting cells. Front Immunol. (2022) 13:954936. doi: 10.3389/fimmu.2022.954936

PubMed Abstract | Crossref Full Text | Google Scholar

73. Otterlei Fjørtoft M, Huse K, and Rye IH. The tumor immune microenvironment in breast cancer progression. Acta Oncol. (2024) 63:359–67. doi: 10.2340/1651-226X.2024.33008

PubMed Abstract | Crossref Full Text | Google Scholar

74. Schaafsma E, Jiang C, and Cheng C. B cell infiltration is highly associated with prognosis and an immune-infiltrated tumor microenvironment in neuroblastoma. J Cancer Metastasis Treat. (2021) 7:55–66. doi: 10.20517/2394-4722.2021.72

PubMed Abstract | Crossref Full Text | Google Scholar

75. Roberti MP, Mordoh J, and Levy EM. Biological role of NK cells and immunotherapeutic approaches in breast cancer. Front Immunol. (2012) 3:375. doi: 10.3389/fimmu.2012.00375

PubMed Abstract | Crossref Full Text | Google Scholar

76. Portale F and Di Mitri D. NK cells in cancer: mechanisms of dysfunction and therapeutic potential. Int J Mol Sci. (2023) 24:9521. doi: 10.3390/ijms24119521

PubMed Abstract | Crossref Full Text | Google Scholar

77. Buque A, Bloy N, Petroni G, Kroemer G, and Galluzzi L. NK cells beat T cells at early breast cancer control. Oncoimmunology. (2020) 9:1806010. doi: 10.1080/2162402X.2020.1806010

PubMed Abstract | Crossref Full Text | Google Scholar

78. Hollern DP, Xu N, Thennavan A, Glodowski C, Garcia-Recio S, Mott KR, et al. B cells and T follicular helper cells mediate response to checkpoint inhibitors in high mutation burden mouse models of breast cancer. Cell. (2019) 179:1191–206.e21. doi: 10.1016/j.cell.2019.10.028

PubMed Abstract | Crossref Full Text | Google Scholar

79. Kong X, Duan Y, Sang Y, Li Y, Zhang H, Liang Y, et al. LncRNA-CDC6 promotes breast cancer progression and function as ceRNA to target CDC6 by sponging microRNA-215. J Cell Physiol. (2019) 234:9105–17. doi: 10.1002/jcp.27587

PubMed Abstract | Crossref Full Text | Google Scholar

80. Yang S, Wang X, Zhou X, Hou L, Wu J, Zhang W, et al. ncRNA-mediated ceRNA regulatory network: Transcriptomic insights into breast cancer progression and treatment strategies. Biomed Pharmacother. (2023) 162:114698. doi: 10.1016/j.biopha.2023.114698

PubMed Abstract | Crossref Full Text | Google Scholar

81. Li HC, Chen YF, Feng W, Cai H, Mei Y, Jiang YM, et al. Loss of the Opa interacting protein 5 inhibits breast cancer proliferation through miR-139-5p/NOTCH1 pathway. Gene. (2017) 603:1–8. doi: 10.1016/j.gene.2016.11.046

PubMed Abstract | Crossref Full Text | Google Scholar

82. Abd EAS, Ghanem HM, Swellam M, and Taha AM. Involvement of FAM170B-AS1, hsa-miR-1202, and hsa-miR-146a-5p in breast cancer. Cancer biomark. (2024) 39:313–33. doi: 10.3233/CBM-230396

PubMed Abstract | Crossref Full Text | Google Scholar

83. Warady BA, Pergola PE, Agarwal R, Andreoli S, Appel GB, Bangalore S, et al. Effects of bardoxolone methyl in alport syndrome. Clin J Am Soc Nephrol: CJASN. (2022) 17:1763–74. doi: 10.2215/CJN.02400222

PubMed Abstract | Crossref Full Text | Google Scholar

84. Wang XR, Jiang ZB, Xu C, Meng WY, Liu P, Zhang YZ, et al. Andrographolide suppresses non-small-cell lung cancer progression through induction of autophagy and antitumor immune response. Pharmacol Res. (2022) 179:106198. doi: 10.1016/j.phrs.2022.106198

PubMed Abstract | Crossref Full Text | Google Scholar

85. Smolewski P. Terameprocol, a novel site-specific transcription inhibitor with anticancer activity. IDrugs. (2008) 11:204–14. doi: 10.1002/hup.898

PubMed Abstract | Crossref Full Text | Google Scholar

86. Sun Y, Giacalone NJ, and Lu B. Terameprocol (tetra-O-methyl nordihydroguaiaretic acid), an inhibitor of Sp1-mediated survivin transcription, induces radiosensitization in non-small cell lung carcinoma. J Thorac Oncol. (2011) 6:8–14. doi: 10.1097/JTO.0b013e3181fa646a

PubMed Abstract | Crossref Full Text | Google Scholar

87. Ruhl T, Corsten C, Beier JP, and Kim BS. The immunosuppressive effect of the endocannabinoid system on the inflammatory phenotypes of macrophages and mesenchymal stromal cells: a comparative study. Pharmacol Rep. (2021) 73:143–53. doi: 10.1007/s43440-020-00166-3

PubMed Abstract | Crossref Full Text | Google Scholar

88. Shoshina OO, Kozhin PM, Shadrin VS, Romashin DD, Rusanov AL, and Luzgina NG. Phenotypic features of mesenchymal stem cell subpopulations obtained under the influence of various toll-like receptors ligands. Bull Exp Biol Med. (2021) 170:555–9. doi: 10.1007/s10517-021-05105-7

PubMed Abstract | Crossref Full Text | Google Scholar

89. Waterman RS, Henkle SL, and Betancourt AM. Mesenchymal stem cell 1 (MSC1)-based therapy attenuates tumor growth whereas MSC2-treatment promotes tumor growth and metastasis. PloS One. (2012) 7:e45590. doi: 10.1371/journal.pone.0045590

PubMed Abstract | Crossref Full Text | Google Scholar

90. Pleyer L, Valent P, and Greil R. Mesenchymal stem and progenitor cells in normal and dysplastic hematopoiesis-masters of survival and clonality? Int J Mol Sci. (2016) 17:1009. doi: 10.3390/ijms17071009

PubMed Abstract | Crossref Full Text | Google Scholar

91. Ruan H, Wang Z, Sun Z, Wei J, Zhang L, Ju H, et al. Single-cell RNA sequencing reveals the characteristics of cerebrospinal fluid tumour environment in breast cancer and lung cancer leptomeningeal metastases. Clin Trans Med. (2022) 12:e885. doi: 10.1002/ctm2.885

PubMed Abstract | Crossref Full Text | Google Scholar

92. Yu M, Bardia A, Wittner BS, Stott SL, Smas ME, Ting DT, et al. Circulating breast tumor cells exhibit dynamic changes in epithelial and mesenchymal composition. Sci (New York NY). (2013) 339:580–4. doi: 10.1126/science.1228522

PubMed Abstract | Crossref Full Text | Google Scholar

93. He P, Ma Y, Wu Y, Zhou Q, and Du H. Exploring PANoptosis in breast cancer based on scRNA-seq and bulk-seq. Front Endocrinol. (2023) 14:1164930. doi: 10.3389/fendo.2023.1164930

PubMed Abstract | Crossref Full Text | Google Scholar

Keywords: breast cancer, telomere maintenance mechanism, stem cell, whole transcriptome, single-cell RNA sequencing, biomarkers

Citation: Zhen S, Huang L, Zhu Q, Chen R, Wang J and Zha X (2025) Identification and validation of genes related to stem cells and telomere maintenance mechanisms as biomarkers for breast cancer. Front. Immunol. 16:1618193. doi: 10.3389/fimmu.2025.1618193

Received: 25 April 2025; Accepted: 17 June 2025;
Published: 01 July 2025.

Edited by:

April Kartikasari, RMIT University, Australia

Reviewed by:

Jindong Xie, Sun Yat-sen University Cancer Center (SYSUCC), China
Binbin Guo, Soochow University, China

Copyright © 2025 Zhen, Huang, Zhu, Chen, Wang and Zha. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Xiaoming Zha, bmp6aGF4bUBxcS5jb20=

These authors have contributed equally to this work and share first authorship

ORCID: Xiaoming Zha, orcid.org/0009-0007-0910-2539

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.