Exploring key biomarkers associated with axillary lymph node metastasis in breast cancer using single-cell RNA sequencing and Mendelian randomization

Abstract

Background

Axillary lymph node metastasis (ALNM) serves as a critical prognostic determinant in breast cancer, yet the molecular drivers governing lymphatic dissemination remain poorly characterized. Integrating single-cell transcriptomic profiling with Mendelian-randomization (MR)-based genetic prioritization may help reveal cell type-specific mechanisms underlying metastatic progression.

Methods

We analyzed the GSE195861 single-cell RNA sequencing dataset encompassing six invasive ductal carcinoma (IDC) samples and paired ALNM specimens. t-distributed Stochastic Neighbor Embedding-based clustering and SingleR annotation delineated cellular heterogeneity, while differential expression analysis identified metastasis-associated genes in epithelial compartments. MR analysis employing five robust methods (inverse variance-weighted, weighted median, MR-Egger, simple/weighted mode) integrated genome-wide association study data (GCST90018799) to establish causal gene-breast cancer associations. CellChat reconstructed ligand-receptor networks across nine annotated cell types.

Results

Unsupervised clustering resolved 27 cell clusters into nine lineages, revealing ALNM-specific expansion of monocytes, pre-B cells, and CD34+ hematopoietic stem cells (HSCs). Epithelial cells exhibited 2421 differentially expressed genes (DEGs) between IDC and ALNM, including 12 genes whose genetically predicted expression showed significant associations with breast cancer risk in MR analysis (P < 0.05). CD53 (odds ratio (OR) = 1.110, 95% confidence interval (CI) = 1.019–1.209, P = 0.017) and TCDD-inducible poly-ADP-ribose polymerase (TIPARP) (OR = 1.153, 95% CI = 1.032–1.288, P = 0.012) were prioritized as candidate genes, as their genetically predicted expression was associated with increased breast cancer risk in weighted median MR. Cell–cell communication analysis implicated macrophage-derived midkine-nucleolin signaling and B-cell-orchestrated macrophage migration inhibitory factor-(CD74 + CXCR4) axis in metastatic crosstalk. Functional enrichment linked DEGs to extracellular matrix remodeling and MAPK/PI3K-Akt activation.

Conclusion

This multi-omics integration prioritizes CD53 and TIPARP as ALNM-associated candidate genes with genetically supported associations with breast cancer risk, with macrophage-epithelial and B-cell-HSC interactions serving as potential therapeutic targets. Our findings provide a roadmap for developing metastasis-interceptive strategies through precision targeting of the ALNM-associated tumor microenvironment.

Keywords

breast cancer axillary lymph node metastasis single-cell RNA sequencing Mendelian randomization biomarker

Introduction

Breast cancer (BC) remains the most prevalent malignancy in women globally, with axillary lymph node metastasis (ALNM) serving as a pivotal determinant of disease staging, therapeutic decision-making, and patient survival.^1–3 Despite advancements in multimodal therapies, approximately 20–30% of patients with early-stage invasive ductal carcinoma (IDC) develop lymphatic metastases, as evidenced by large-scale clinical cohorts, underscoring persistent gaps in understanding the molecular drivers of metastatic dissemination.^4,5 The metastatic cascade—from local invasion to lymphovascular infiltration and colonization of distant niches—involves intricate crosstalk between neoplastic epithelial cells and stromal/immune components within the tumor microenvironment (TME).⁶ However, traditional bulk transcriptomic approaches obscure cellular heterogeneity, limiting insights into metastatic cell states and intercellular signaling networks.

The emergence of single-cell RNA sequencing (scRNA-seq) has significantly advanced our ability to dissect the TME at single-cell resolution, uncovering previously unrecognized epithelial cell subpopulations with metastatic potential and complex immune–stromal interactions.^7–9 In BC, recent scRNA-seq studies have identified epithelial cell subsets exhibiting hybrid epithelial–mesenchymal transition (EMT) features, as well as immunosuppressive myeloid cells that contribute to the formation of metastatic niches.^10–12 However, a key knowledge gap remains: the mechanisms by which spatially and temporally regulated cell–cell communication within primary tumors promotes the selection and expansion of tumor cell populations with a predisposition for lymph node metastasis are still poorly understood.

Mendelian randomization (MR), a genetic instrumental variable approach, offers a powerful framework to infer causal relationships between molecular traits and disease outcomes while minimizing confounding biases.^13,14 When synergized with scRNA-seq-derived cell type-specific expression profiles, MR can help prioritize candidate genes identified from single-cell analyses by providing orthogonal genetic support for their relevance to BC biology.¹⁵

Here, we present a multi-omics investigation that integrates scRNA-seq of paired IDC and ALNM samples with MR analysis of genome-wide association study (GWAS) data to dissect the cellular and genetic landscape of BC lymphatic metastasis. Focusing on epithelial cells—the central drivers of metastatic spread—we explore how the TME is remodeled during lymph node colonization, unravel ligand–receptor signaling networks that mediate metastatic communication, and identify causal gene-disease associations through MR analysis. By bridging single-cell biology with genetic causal inference, our study offers novel insights and potential therapeutic targets for intercepting lymphatic metastasis.

Methods

Data source

The single-cell dataset GSE195861 was downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo). This dataset includes samples from seven cases of ductal carcinoma in situ (DCIS), six cases of IDC, six cases of ALNM, and one normal breast tissue sample. According to the GEO record, each of the six IDC cases was accompanied by a corresponding lymph node metastasis specimen collected from the same patient (Supplementary Table 1). Therefore, the comparative analyses in this study were based on six matched IDC–ALNM pairs. Gene-specific eQTL exposure datasets were obtained from the IEU OpenGWAS eqtl-a series, corresponding to the eQTLGen resource. According to the source database, these datasets represent blood-derived/whole-blood gene expression data from individuals of predominantly European ancestry and are aligned to HG19/GRCh37. The eQTLGen resource includes up to 31,684 samples from 37 cohorts, although the effective sample size may vary across specific gene-level exposure datasets available in OpenGWAS. BC GWAS summary statistics (ebi-a-GCST90018799) were used as the outcome dataset.¹⁶

scRNA-seq data processing and normalization

scRNA-seq is a transformative technology that provides insights into intercellular heterogeneity, high-resolution disease mechanisms, and potential clinical applications.¹⁷ In GSE195861, the Seurat package was used to create a Seurat object for scRNA-seq data.¹⁸ The PercentageFeatureSet function was applied to calculate the percentage of mitochondrial genes. Cells with fewer than 200 detected genes and genes expressed in fewer than three cells were filtered out. The data were then normalized, and genes with the highest cell-to-cell variation coefficients were selected for downstream analysis.

Dimensionality reduction, unsupervised clustering, and visualization

PC analysis was employed for dimensionality reduction of the integrated dataset.¹⁹ The variance explained by each principal component (PC) was ranked, and the top 20 PCs (dims = 20) were selected for downstream analysis. Graph-based clustering was performed using modularity optimization, and t-distributed Stochastic Neighbor Embedding (t-SNE) was used for visualization.²⁰ The SingleR package was utilized for cell type annotation, and the Monocle package was used for cell trajectory analysis.

Cell–cell communication analysis

Cell–cell communication was evaluated by quantifying ligand-receptor interactions between different cell types. The CellChat R package (www.cellchat.org) was used to infer intercellular communication, leveraging the CellChatDB database.²¹ This approach integrates expression data with known ligand-receptor interactions, including cofactors, to model cell–cell communication networks.

Identification of differentially expressed genes

Differential expression analysis between IDC and axillary lymph node metastatic epithelial cells was conducted using the FindMarkers function in the Seurat package. The filtering criteria for differentially expressed genes (DEGs) were adjusted P-value (P.adj) < 0.05 and |log2 fold change (log2FC)| > 1.

Functional enrichment analysis

To explore the functional roles of DEGs, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed.^22,23 These analyses identified common biological functions and pathways associated with the DEGs.

MR analysis

To genetically prioritize epithelial DEGs identified from the IDC-versus-ALNM comparison, we performed two-sample MR using eQTLs as instrumental variables and BC risk (ebi-a-GCST90018799) as the outcome. Single nucleotide polymorphisms (SNPs) associated with gene expression at P < 5 × 10⁻⁸ were selected as candidate instruments. To ensure independence among instruments, linkage disequilibrium clumping was performed using a window of 10,000 kb and an r² threshold of 0.001. Instrument strength was assessed by calculating R² and F-statistics, and SNPs with F ≤ 10 were excluded as weak instruments. Exposure and outcome datasets were harmonized using the standard procedure implemented in the TwoSampleMR package. Palindromic SNPs were handled according to allele frequency information, and ambiguous variants were excluded when strand alignment could not be reliably determined. The inverse variance-weighted (IVW) method was used as the primary analysis, complemented by weighted median, MR-Egger, simple mode, and weighted mode methods. Cochran's Q test was used to assess heterogeneity, the MR-Egger intercept was used to assess directional pleiotropy, leave-one-out analysis was used to evaluate robustness, and MR-PRESSO was additionally performed when the number of retained instrumental variables was sufficient. A two-sided P value < 0.05 was considered statistically significant²⁴ (Supplementary Tables 2 and 3).

Gene–gene interaction network

GeneMANIA (https://genemania.org/) was used to construct a gene–gene interaction network to identify genes functionally related to the identified biomarkers.²⁵ This tool integrates publicly available biological datasets to infer gene relationships.

Statistical analysis

MR results are reported as odds ratios (ORs) with 95% confidence intervals (CIs). After harmonization of exposure and outcome datasets, the IVW method was used as the primary MR analysis, complemented by weighted median, MR-Egger, simple mode, and weighted mode approaches. Cochran's Q test was used to assess heterogeneity among instrumental variables. Directional pleiotropy was evaluated using the MR-Egger intercept, with funnel plots used as a visual aid to inspect potential asymmetry. Leave-one-out analysis was performed to determine whether the observed associations were driven by any single SNP. A two-sided P value < 0.05 was considered statistically significant. Group comparisons were conducted using the Mann–Whitney U test.

Results

Identification of cell clusters and composition in IDC and metastatic axillary lymph nodes

In this study, we analyzed six IDC samples and their paired six axillary lymph node metastases from the GSE195861 single-cell dataset. t-SNE clustering identified 27 distinct cell clusters, which were further classified into nine major cell types using SingleR annotation: monocytes, pre-B cells (CD34−), epithelial cells, granulocyte colony-stimulating factor-mobilized hematopoietic stem cells (G-CSF-HSC), macrophages, bone marrow (BM) cells, CD34+ hematopoietic stem cells (HSC CD34+), neutrophils, and stromal cells (Figure 1(a) and (b)).

Figure 1.

Cell clustering and trajectory analysis of single-cell RNA-Seq data. (a) t-SNE plot showing cell clustering of all samples in the single-cell RNA-Seq dataset. Each color represents a distinct cell cluster. (b) Annotation of cell types using the SingleR package, identifying nine major cell types. (c) Comparison of cell type distributions between IDC and matched metastatic axillary lymph nodes. (d) and (e) Cell trajectory analysis using Monocle. Black curves indicate inferred differentiation trajectories, and numbered nodes (1, 2) represent potential critical transition states. Pseudotime denotes the inferred temporal progression of cellular states.

Compared to the IDC group, the metastatic lymph node group exhibited a significantly higher proportion of monocytes, pre-B cells, epithelial cells, G-CSF-HSC, macrophages, bone marrow cells, and CD34+ hematopoietic stem cells (Figure 1(c)). Cell trajectory analysis revealed dynamic changes in the tumor microenvironment, with black curves representing differentiation trajectories and numbered nodes (1, 2) indicating potential critical transition states (Figure 1(d) and (e)). The overall trend suggested that epithelial cells might undergo specific differentiation or phenotypic changes during ALNM, which could be closely related to tumor progression and metastatic potential.

Cell–cell communication in IDC and metastatic axillary lymph nodes

To investigate intercellular communication, we performed cell–cell interaction analysis on the nine identified cell types. The results showed that macrophages and B cells had the highest number of ligand-receptor interactions, while monocytes and B cells exhibited the strongest communication intensity (Supplementary Figure 1).

At the signaling pathway level, we constructed an interaction network to depict cell–cell communications. Bubble plots indicated that the macrophage migration inhibitory factor (MIF)-(CD74 + CXCR4) interaction was enriched between B cells and G-CSF-HSC, whereas the MIF-(CD74 + CD44) interaction was predominant between B cells and monocytes. Additionally, the MDK–NCL interaction between macrophages and epithelial cells was particularly strong, suggesting that these interactions may play a crucial role in ALNM in BC (Figure 2).

Figure 2.

Bubble plot of cell–cell communication at the signaling pathway level. A bubble plot was constructed to visualize the intercellular communication network at the signaling pathway level. Each bubble represents a ligand–receptor interaction between specific cell types, with bubble size indicating interaction strength and color denoting pathway enrichment.

Identification of key differentially expressed genes and MR analysis

Given that epithelial cells are the primary source of tumor cells and are likely to be closely associated with ALNM, we conducted a differential gene expression analysis between epithelial cells from IDC and their paired metastatic lymph nodes. A total of 2421 DEGs were identified.

To further explore the role of these DEGs, we obtained gene-specific eQTL instruments for these DEGs from the IEU OpenGWAS eqtl-a resource and performed MR analysis using BC (ebi-a-GCST90018799) as the outcome. Twelve DEGs were significantly associated with BC (p < 0.05), including ACAP1, APOL6, ASRGL1, CD53, FSCN1, ISG15, PI3, PLEKHO1, RNF144B, SPART, TCDD-inducible poly-ADP-ribose polymerase (TIPARP), and TTC23 (Supplementary Tables 2 and 3, Figure 3).

Figure 3.

Forest plot of expression quantitative trait loci (eQTL)-based MR results for ALNM-associated epithelial DEGs using breast cancer (ebi-a-GCST90018799) as the outcome. eQTLs of candidate central genes were extracted from GWAS summary statistics and analyzed using Mendelian randomization (MR), with breast cancer (GCST90018799) as the outcome. Sensitivity analysis was performed by comparing causal effect estimates across multiple MR methods, including MR-Egger, penalized weighted median, simple mode, inverse variance weighted (IVW), and weighted mode. A P-value < 0.05 was considered statistically significant.

Among them, ACAP1, CD53, PLEKHO1, and TIPARP showed positive associations between genetically predicted expression and BC risk (OR > 1), while the remaining genes appeared to have protective effects. The best causal estimates showed that CD53 (OR (95% CI) = 1.110 (1.019–1.209), P = 0.017, weighted median) and TIPARP (OR (95% CI) = 1.153 (1.032–1.288), P = 0.012, weighted median) were significantly associated with an increased risk of BC. Leave-one-out sensitivity analysis demonstrated the robustness of these findings, supporting the potential relevance of these genes in BC biology and their possible involvement in ALNM-associated molecular alterations (Figure 4).

Figure 4.

Mendelian randomization analyses of CD53 and TIPARP, two candidate genes prioritized from ALNM-associated epithelial DEGs, using breast cancer risk as the outcome. (a) and (e) Scatter plots showing the associations between genetically predicted expression of CD53 and TIPARP and breast cancer risk. (b) and (f) Forest plots showing effect estimates from multiple Mendelian randomization (MR) methods for CD53 and TIPARP. (c) and (g) Funnel plots assessing potential asymmetry suggestive of directional pleiotropy in the MR analyses of CD53 and TIPARP. (d) and (h) Leave-one-out sensitivity analyses evaluating the robustness of the MR results for CD53 and TIPARP.

Functional enrichment and pathway analysis of key genes

Using GeneMANIA, we predicted 32 functionally related genes associated with the 12 identified DEGs, with RNF144B showing the highest correlation. The protein–protein interaction (PPI) network revealed various types of associations, including co-expression, co-localization, and genetic interactions.

GO enrichment analysis indicated that these genes were involved in gland development, response to lipopolysaccharide, and collagen-containing extracellular matrix organization. KEGG pathway analysis showed significant enrichment in MAPK signaling and PI3K-Akt signaling pathways, suggesting their potential involvement in the molecular mechanisms underlying BC ALNM (Figure 5).

Figure 5.

Functional enrichment and pathway analysis of key genes. (a) Gene–gene interaction network of identified biomarker genes constructed using the GeneMANIA database. (b) and (c) Gene Ontology (GO) enrichment analysis of key genes, including biological processes, cellular components, and molecular functions. (d) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of key genes.

Discussion

In this study, we employed a comprehensive approach combining scRNA-seq, functional enrichment analysis, and MR to investigate the molecular mechanisms driving ALNM in BC. The integration of scRNA-seq data with clinical and genetic information provided deep insights into the cellular dynamics and gene interactions involved in BC metastasis, highlighting several key findings that could shape future diagnostic and therapeutic strategies.

Our scRNA-seq atlas of IDC and matched ALNM samples reveals a profound reorganization of the TME during lymphatic dissemination. The enrichment of monocytes, G-CSF-HSC, and CD34+ hematopoietic progenitors in ALNM aligns with emerging evidence that metastatic niches recruit BM-derived cells to support colonization.^26–31 Notably, the expansion of pre-B cells in ALNM echoes recent reports of B cell-mediated immune tolerance in metastatic sites, where regulatory B cells (Bregs) may suppress anti-tumor T-cell responses through interleukin-10 and transforming growth factor-beta (TGF-β) secretion.^32,33 The observed epithelial cell trajectory bifurcation (nodes 1→2) suggests a metastable transitional state—possibly reflecting EMT-like plasticity—that warrants further investigation using pseudotime-velocity analyses.

The dominance of macrophage-epithelial interactions via MDK–NCL signaling provides mechanistic insights into niche preparation. MDK, a heparin-binding growth factor overexpressed in BC, promotes metastasis through PI3K/Akt activation and vascular leakiness, while its receptor nucleolin facilitates extracellular MDK internalization to drive pro-survival signaling.^34,35 Similarly, the MIF-(CD74 + CXCR4) axis between B cells and HSCs may orchestrate stemness programs in metastatic epithelial cells, as MIF is known to enhance cancer stem cell properties via CXCR4-mediated Wnt/β-catenin activation.^36,37 These findings position macrophage- and B cell-derived ligands as microenvironmental “fertilizers” for metastatic outgrowth.

Our MR analysis provided genetic support for prioritizing CD53 and TIPARP as candidate genes associated with BC risk, and these genes may also be relevant to ALNM-associated molecular alterations. CD53, a tetraspanin regulating immune synapse formation in lymphocytes, paradoxically emerges as a metastasis promoter—a duality potentially explained by its role in stabilizing integrin clusters that enhance epithelial cell adhesion to lymphatic endothelia.^38–40 TIPARP (TCDD-inducible poly-ADP-ribose polymerase), traditionally studied in xenobiotic metabolism, may foster metastasis through poly-ADP-ribosylation-dependent modulation of hypoxia-inducible factor 1-alpha (HIF-1α) stability under hypoxic stress.^41,42 The protective association of APOL6 and ISG15 aligns with their roles in ER stress-induced apoptosis and interferon-mediated immune surveillance, respectively, suggesting their loss in metastasis enables immune evasion.^43–46

The convergence of DEGs on MAPK and PI3K-Akt pathways underscores their dual roles in BC progression: while MAPK signaling (e.g. via ERK1/2 phosphorylation) drives EMT and matrix remodeling, PI3K-Akt activation promotes survival in detached circulating tumor cells.^47–49 Pharmacologically targeting these pathways has shown limited success due to compensatory feedback loops; however, our cell-type-specific expression data suggest that combinatorial inhibition of macrophage-derived MDK and epithelial Akt may disrupt metastatic resilience. Additionally, the association of RNF144B—an E3 ubiquitin ligase regulating epidermal growth factor receptor degradation—with extracellular matrix organization implies a novel link between growth factor receptor turnover and metastatic niche formation.^50–52

Despite these advances, our study has limitations. The modest cohort size (n = 6 pairs) necessitates validation in larger scRNA-seq datasets. Spatial transcriptomics could further resolve ligand-receptor interactions within microanatomical niches, while CRISPR-based perturbation of CD53/TIPARP in organoid models would clarify their functional roles. Future MR analyses should incorporate cis-eQTLs with pleiotropy-robust methods to refine causal estimates.

In conclusion, this study provides novel insights into the cellular and molecular mechanisms of ALNM in BC. The identification of key genes, cellular interactions, and signaling pathways opens new avenues for the development of targeted therapies aimed at preventing or treating lymph node metastasis. Future studies should focus on validating these findings and exploring the potential therapeutic implications of targeting these molecular drivers in clinical practice.

Supplemental Material

sj-docx-1-jbm-10.1177_03936155261443650 - Supplemental material for Exploring key biomarkers associated with axillary lymph node metastasis in breast cancer using single-cell RNA sequencing and Mendelian randomization

Supplemental material, sj-docx-1-jbm-10.1177_03936155261443650 for Exploring key biomarkers associated with axillary lymph node metastasis in breast cancer using single-cell RNA sequencing and Mendelian randomization by Limeng Qu, Jinyang Li, Shirong Ding, Qian Long and Wenjun Yi in The International Journal of Biological Markers

Footnotes

Abbreviations

Acknowledgments

Thanks to the anonymous peer-reviewers for their insightful suggestions and careful reading of the manuscript.

ORCID iD

Qian Long

Ethics approval

Not applicable.

Consent to participate

Not applicable.

Author contributions

Q-LM and L-JY: Writing—original draft. D-SR: Writing—original draft. L-Q: Writing—review & editing, Supervision. Y-WJ: Writing—review & editing, Conceptualization.

Funding

This study was supported by the National Natural Science Foundation of China [grant number: 82303526 and 82403073],the Natural Science Foundation of the Hunan Province of China [grant number: 2023JJ40842 and 2024JJ6590], the Innovation Platform and Talent Plan of Hunan Province [grant number: 2023SK4019], the Natural Science Foundation of Changsha City [grant number: kq2208309], the China Postdoctoral Science Foundation [grant number: 2023M733955 and 2023M743946], Open Funds of State Key Laboratory of Oncology in South China [grant number: HN2024-04, NH2024-07], and the Scientific Research Launch Project for new employees of the Second Xiangya Hospital of Central South University [grant number: QH20230256 and QH20230268].

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Supplemental material

Supplemental material for this article is available online.

References

Gong

Wang

, et al. Global trends and forecasts of breast cancer incidence and deaths. Sci Data 2023; 10: 334.

Giaquinto

Sung

Newman

, et al. Breast cancer statistics 2024. CA Cancer J Clin 2024; 74: 477–495.

Bartels

SAL

Donker

Poncet

, et al. Radiotherapy or surgery of the axilla after a positive sentinel node in breast cancer: 10-year results of the randomized controlled EORTC 10981-22023 AMAROS trial. J Clin Oncol 2023; 41: 2159–2165.

Park

Lee

Kong

S-Y

, et al. Use of adjuvant chemotherapy in hormone receptor-positive breast cancer patients with or without the 21-gene expression assay. Breast Cancer Res Treat 2018; 170: 69–76.

Donker

van Tienhoven

Straver

, et al. Radiotherapy or surgery of the axilla after a positive sentinel node in breast cancer (EORTC 10981-22023 AMAROS): a randomised, multicentre, open-label, phase 3 non-inferiority trial. Lancet Oncol 2014; 15: 1303–1310.

Harris

Savas

Virassamy

, et al. Towards targeting the breast cancer immune microenvironment. Nat Rev Cancer 2024; 24: 554–577.

Azizi

Carr

Plitas

, et al. Single-cell map of diverse immune phenotypes in the breast tumor microenvironment. Cell 2018; 174: 1293–1308.e36.

Al-Eryani

Roden

, et al. A single-cell and spatially resolved atlas of human breast cancers. Nat Genet 2021; 53: 1334–1347.

Qian

Olbrecht

Boeckx

, et al. A pan-cancer blueprint of the heterogeneous tumor microenvironment revealed by single-cell profiling. Cell Res 2020; 30: 745–762.

10.

Yang

Peng

, et al. Pan-cancer spatially resolved single-cell analysis reveals the crosstalk between cancer-associated fibroblasts and tumor microenvironment. Mol Cancer 2023; 22: 170.

11.

Luo

Yang

Mastoraki

, et al. Single-cell RNA sequencing identifies molecular biomarkers predicting late progression to CDK4/6 inhibition in patients with HR+/HER2- metastatic breast cancer. Mol Cancer 2025; 24: 48.

12.

Bao

Shi

Zhao

, et al. Integrated analysis of single-cell RNA-Seq and bulk RNA-Seq unravels tumour heterogeneity plus M2-like tumour-associated macrophage infiltration and aggressiveness in TNBC. Cancer Immunol Immunother 2021; 70: 189–202.

13.

Sundquist

Zhang

, et al. Mitochondrial related genome-wide Mendelian randomization identifies putatively causal genes for multiple cancer types. EBioMedicine 2023; 88: 104432.

14.

Zhu

Zhang

, et al. C-reactive protein and cancer risk: a pan-cancer study of prospective cohort and Mendelian randomization analysis. BMC Med 2022; 20: 301.

15.

Zhang

Zhu

, et al. The association between immune cells and breast cancer: insights from Mendelian randomization and meta-analysis. Int J Surg 2025; 111: 230–241.

16.

Manolio

Collins

Cox

, et al. Finding the missing heritability of complex diseases. Nature 2009; 461: 747–753.

17.

Butler

Hoffman

Smibert

, et al. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol 2018; 36: 411–420.

18.

Satija

Farrell

Gennert

, et al. Spatial reconstruction of single-cell gene expression data. Nat Biotechnol 2015; 33: 495–502.

19.

Ringnér

. What is principal component analysis? Nat Biotechnol 2008; 26: 303–304.

20.

Kleino

Frolovaitė

Suomi

, et al. Computational solutions for spatial transcriptomics. Comput Struct Biotechnol J 2022; 20: 4870–4884.

21.

Jin

Plikus

Nie

. Cellchat for systematic analysis of cell-cell communication from single-cell transcriptomics. Nat Protoc 2025; 20: 180–219.

22.

The Gene Ontology Consortium. Gene ontology consortium: going forward. Nucleic Acids Res 2015; 43: D1049–D1056.

23.

Kanehisa

Goto

. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res 2000; 28: 27–30.

24.

Carter

Sanderson

Hammerton

, et al. Mendelian randomisation for mediation analysis: current methods and challenges for implementation. Eur J Epidemiol 2021; 36: 465–478.

25.

Franz

Rodriguez

Lopes

, et al. GeneMANIA update 2018. Nucleic Acids Res 2018; 46: W60–W64.

26.

Qian

B-Z

Zhang

, et al. CCL2 Recruits inflammatory monocytes to facilitate breast-tumour metastasis. Nature 2011; 475: 222–225.

27.

Kim

Lee

I-H

Cho

, et al. Spatiotemporal evolution of the primary glioblastoma genome. Cancer Cell 2015; 28: 318–328.

28.

Reynaud

Pietras

Barry-Holson

, et al. IL-6 controls leukemic multipotent progenitor cell fate and contributes to chronic myelogenous leukemia development. Cancer Cell 2011; 20: 661–673.

29.

Pradeep

Huang

Mora

, et al. Erythropoietin stimulates tumor growth via EphB4. Cancer Cell 2015; 28: 610–622.

30.

Alderton

. Tumour metabolism: feeding the TCA cycle in vivo. Nat Rev Cancer 2016; 16: 198.

31.

Dong

, et al. R-2HG exhibits anti-tumor activity by targeting FTO/m6A/MYC/CEBPA signaling. Cell 2018; 172: 90–105.e23.

32.

Cook

Metheny-Barlow

Tallant

, et al. Angiotensin-(1-7) reduces fibrosis in orthotopic breast tumors. Cancer Res 2010; 70: 8319–8328.

33.

Satoguina

Adjobimey

Arndts

, et al. Tr1 and naturally occurring regulatory T cells induce IgG4 in B cells through GITR/GITR-L interaction, IL-10 and TGF-beta. Eur J Immunol 2008; 38: 3101–3113.

34.

Visvader

. Midkine links aging with breast cancer-A new predictor of cancer risk. Cancer Cell 2024; 42: 1815–1817.

35.

Yan

Jimenez

, et al. Midkine as a driver of age-related changes and increase in mammary tumorigenesis. Cancer Cell 2024; 42: 1936–1954.e9.

36.

Yan

Wang

, et al. Breast cancer stem cells secrete MIF to mediate tumor metabolic reprogramming that drives immune evasion. Cancer Res 2024; 84: 1270–1285.

37.

Zhang

Chen

Wang

, et al. Macrophage migration inhibitory factor promotes proliferation and neuronal differentiation of neural stem/precursor cells through Wnt/β-catenin signal pathway. Int J Biol Sci 2013; 9: 1108–1120.

38.

Dunlock

V-ME

Arp

Singh

, et al. Tetraspanin CD53 controls T cell immunity through regulation of CD45RO stability, mobility, and function. Cell Rep 2022; 39: 111006.

39.

Okochi

Kato

Nashiro

, et al. Expression of tetra-spans transmembrane family (CD9, CD37, CD53, CD63, CD81 and CD82) in normal and neoplastic human keratinocytes: an association of CD9 with alpha 3 beta 1 integrin. Br J Dermatol 1997; 137: 856–863.

40.

Yeung

Anderson

JML

Wee

, et al. Leukocyte tetraspanin CD53 restrains α3 integrin mobilization and facilitates cytoskeletal remodeling and transmigration in mice. J Immunol 2020; 205: 521–532.

41.

Zhang

Cao

Dong

, et al. TiPARP forms nuclear condensates to degrade HIF-1α and suppress tumorigenesis. Proc Natl Acad Sci U S A 2020; 117: 13447–13456.

42.

NarenYa, Feng

, et al. A preliminary study of effect of melatonin on inflammation and hypoxia-related factors in a mouse model of elastase-induced intracranial aneurysm. Brain Behav. 2025; 15: e70371.

43.

Pays

. The Janus-faced functions of apolipoproteins L in membrane dynamics. Cell Mol Life Sci 2024; 81: 134.

44.

Liu

Chen

, et al. Upregulation of apolipoprotein L6 improves tumor immunotherapy by inducing immunogenic cell death. Biomolecules 2023; 13: 415.

45.

Wardlaw

Petrini

JHJ

. ISG15: A link between innate immune signaling, DNA replication, and genome stability. Bioessays 2023; 45: e2300042.

46.

Perng

Y-C

Lenschow

. ISG15 In antiviral immunity and beyond. Nat Rev Microbiol 2018; 16: 423–439.

47.

Wang

, et al. Role of p38γ MAPK in regulation of EMT and cancer stem cells. Biochim Biophys Acta Mol Basis Dis 2018; 1864: 3605–3617.

48.

Wen

Hou

, et al. Cancer-associated fibroblast (CAF)-derived IL32 promotes breast cancer cell invasion and metastasis via integrin β3-p38 MAPK signalling. Cancer Lett 2019; 442: 320–332.

49.

Kallergi

Agelaki

Kalykaki

, et al. Phosphorylated EGFR and PI3K/Akt signaling kinases are expressed in circulating tumor cells of breast cancer patients. Breast Cancer Res 2008; 10: R80.

50.

Abad

Sandoz

Romero

, et al. The TP53-activated E3 ligase RNF144B is a tumour suppressor that prevents genomic instability. J Exp Clin Cancer Res 2024; 43: 127.

51.

Conforti

Yang

Piro

, et al. PIR2/Rnf144B Regulates epithelial homeostasis by mediating degradation of p21WAF1 and p63. Oncogene 2013; 32: 4758–4765.

52.

Wang

Dai

Jiang

, et al. RBR E3 ubiquitin ligases in tumorigenesis. Semin Cancer Biol 2020; 67: 131–144.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.45 MB