Cation homeostasis-related prognostic genes uncovered by transcriptomic analysis in breast cancer

Introduction

Breast cancer (BC) is a malignant tumor originating from the epithelial tissues of the mammary ducts or lobules (1). According to data from the World Health Organization (WHO), there were 2.3 million new cases of BC and 670,000 BC-related deaths worldwide in 2022 (2). BC originates from the malignant transformation of mammary gland or ductal epithelial cells, and its occurrence is the result of the combined effect of genetic and environmental factors, including female gender, aging, obesity, alcohol consumption and previous history of chest radiotherapy, etc. (3,4). A notable feature of BC is its high heterogeneity, that is, tumors in different patients show significant differences in molecular structure, behavior and treatment response (5). Although current treatment strategies for BC (such as surgery, chemotherapy, radiotherapy, endocrine therapy, and targeted therapy) have achieved significant progress, advanced or metastatic disease still faces challenges such as drug resistance, adverse effects, and high treatment costs, with limited improvement in prognosis (6-8). Therefore, identifying novel prognostic biomarkers and potential therapeutic targets is crucial for achieving precision diagnosis and treatment of BC.

Cation homeostasis refers to the dynamic equilibrium process through which the body maintains relatively constant concentrations of various cations (such as sodium, potassium, calcium, and magnesium) within and outside cells, regulated by the coordinated activity of ion channels, transporters, and pumps—including Na⁺/K⁺-ATPase, calcium pumps, and sodium-calcium exchangers (9). As the core mechanism for maintaining basic physiological processes such as membrane potential, signal transduction, proliferation and apoptosis, its homeostasis imbalance has been confirmed to be closely related to the pathogenesis of BC and various malignant phenotypes (10). In BC, calcium signaling disorder is particularly significant, often driven by specific molecular alterations, such as overexpression of TRPV6 calcium channels and Orai1 calcium release activation channels, which are respectively associated with enhanced proliferation and poor prognosis of tubular type and BC (11). Similarly, the TRPM7 channel that regulates the flow of Ca²⁺ and Mg²⁺ regulates the migration and invasion of metastatic BC cells through the mitogen-activated protein kinase (MAPK) pathway (12). Although a large number of studies have confirmed that the abnormal expression of specific ion channels such as TRPV6 and TRPM7 is closely related to the proliferation, metastasis and drug resistance of BC, there are significant gaps in current research: most of the focus is on the isolated function of a single ion channel, lacking a systematic analysis of the co-regulatory network of cation homeostasis associated genes (CHRGs) and its overall value in the prognosis assessment of BC.

Unlike currently used clinical prognostic markers for BC, which primarily focus on “consequential” features such as tumor proliferation index (e.g., MKI67), specific molecular subtypes, or immune infiltration characteristics, the regulation of cation homeostasis operates at a more upstream level of cell fate determination (13). It is not an isolated biological event but rather a process that maintains fundamental physiological functions—including membrane potential, second messenger signaling, and energy metabolism—through the coordinated regulation of ion channels, transporters, and pumps (14). In tumor cells, the disruption of cation homeostasis often precedes or accompanies the emergence of malignant phenotypes and serves to integrate diverse oncogenic factors such as genetic mutations, metabolic reprogramming, and microenvironmental stress (15,16). Therefore, systematically dissecting the co-regulatory network of cation homeostasis-related genes can not only reveal the core electrophysiological mechanisms underlying tumor biological behavior but may also uncover a novel prognostic evaluation dimension independent of traditional clinicopathological parameters and individual molecular markers. This offers a more integrated theoretical basis for achieving precise classification and intervention in BC.

Based on the above background, this study aims to systematically explore the value of cation homeostasis-related genes in BC prognosis assessment. Through bioinformatics analysis, we constructed and validated a prognostic signature composed of eight cation homeostasis-related genes. This model effectively distinguishes between high- and low-risk groups of BC patients and was further integrated with clinicopathological features to develop a nomogram with good predictive performance. Our work not only confirms the important role of cation homeostasis in BC prognosis stratification but also provides new directions for understanding its tumor biology mechanisms and developing potential targeted intervention strategies. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1-2814/rc).

Methods

Data acquisition and preprocessing

The raw RNA sequencing (RNA-seq) data of The Cancer Genome Atlas-breast invasive carcinoma (TCGA-BRCA) were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) on June 20, 2025, and used as the training set. It included 1,081 BC tissue samples and 99 normal breast tissue samples. Clinical characteristics such as age, grade, stage, T stage, N stage, M stage, and overall survival (OS) were also retrieved. Among these, 1,000 tumor samples had complete survival information. The GSE20685 dataset was obtained from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/gds) as the validation set (Platform: GPL570), which consisted of 327 BC tissue samples, all with available survival data. Additionally, a total of 4,006 cation homeostasis-related genes (CHRGs) were acquired from the Molecular Signatures Database (MSigDB) database (https://www.gsea-msigdb.org/gsea/msigdb) (Table S1). The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Differential expression analysis

Differential expression analysis between BC samples and normal breast tissues in the TCGA-BRCA training set was performed using the “DESeq2” package (v1.42.1) (17). Differentially expressed genes (DEGs) were identified with the threshold |log₂fold change (FC)| >2 and adjusted P value (P adj) <0.05. The “ggplot2” (v3.5.1) (18) and “pheatmap” (v1.0.12) (19) packages were used to generate volcano plots and heatmaps of the top 10 significantly up- and down-regulated genes ranked by |log₂FC|.

Candidate gene screening

To identify cation homeostasis-related DEGs in BC, the “ggvenn” package (v0.1.10) (20) was used to find the intersection between DEGs and CHRGs. The overlapping genes were defined as candidate genes for subsequent analyses.

Functional enrichment and protein-protein interaction (PPI) network analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were conducted using the “clusterProfiler” package (v4.10.1) (21). GO covers three functional categories: molecular function (MF), cellular component (CC), and biological process (BP). Pathways with P<0.05 were considered significantly enriched. The top 10 GO terms and top 30 KEGG pathways were visualized based on P value.

A PPI network was constructed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (http://string-db.org) with an interaction score threshold ≥0.9. The “STRINGdb” package (v3.21) (22) was used for visualization.

Prognostic gene identification

Univariate Cox regression analysis was performed on candidate genes using the “survival” package (v3.8.3) (23) based on TCGA-BRCA samples with survival data [threshold: hazard ratio (HR) ≠1 and P<0.05]. The proportional hazards (PH) assumption was tested (P>0.05), and genes satisfying this condition were considered potential prognostic genes. To refine the model, least absolute shrinkage and selection operator (LASSO) regression was applied using the “glmnet” package (v4.1.8) (24). A five-fold cross-validation was performed to select optimal genes under minimum error.

Prognostic model construction and validation

Each BC patient with complete survival data was assigned a risk score calculated using the formula:

Riskscore=∑i=1ncoef (genei)×expr(genei)[1]

where “coeft” represents the weight of each gene in the model and “Expression” denotes the expression level of each gene. Based on the optimal cut-off value of the risk score determined from the training set, patients were stratified into high-risk and low-risk groups. The distribution of risk scores and their association with survival status were visualized using the “survival” R package (v3.8.3).

To evaluate the model’s performance, the Kaplan-Meier survival analysis with the log-rank test was first applied in the training set to assess the significance of OS differences between the two groups (P<0.05). Subsequently, time-dependent receiver operating characteristic (ROC) curves for 1-, 3-, and 5-year OS were generated using the “pROC” package (v1.18.5) (25). The predictive accuracy of the prognostic model was evaluated by the area under the curve (AUC), with AUC >0.6 considered acceptable.

In the GSE20685 validation dataset, the same risk score cut-off value derived from the training set was applied to stratify patients into high- and low-risk groups. Model reliability was validated from three perspectives: Kaplan-Meier curves were plotted to visualize survival differences between the groups; AUC values for the 1-, 3-, and 5-year ROC curves were calculated. The risk score was confirmed as a statistically significant independent prognostic factor (P<0.05).

Independent prognostic analysis and nomogram construction

Univariate and multivariate Cox regression analyses were performed incorporating the risk score and clinical features [e.g., gender, age, stage, tumor (T) stage, lymph node (N) stage]. Features with P<0.05 in both analyses were considered independent prognostic factors. A nomogram was constructed using the “rms” package (v6.8.1) (26) to predict 1-, 3-, and 5-year OS. Calibration curves and decision curve analysis (DCA) were generated using the “regplot” package (v1.1) (27). ROC analysis (AUC >0.7) was used to evaluate predictive accuracy. Subgroup analyses were conducted to validate the model’s stability. Within each subgroup, patients were stratified using the predefined risk cutoff. Kaplan-Meier curves and time-dependent ROC analyses for 1-, 3-, and 5-year OS were used to evaluate prognostic performance.

Expression validation, chromosomal localization, and correlation

Analysis of Prognostic Genes Wilcoxon test was used to compare expression levels of prognostic genes between BC and normal samples (P<0.05). The “OmicCircos” package (v1.30.0) (28) was employed to visualize the chromosomal locations of prognostic genes. Pearson correlation analysis among prognostic genes was performed using the “psych” package (v2.4.3) (29) (|correlation coefficient (cor)| >0.3, P<0.05).

Cell and culture

The cell lines used in this study included: MCF-10A, SK-BR-3, BT-474, MCF-7, MDA-MB-231, MDA-MB-468, and BT549, obtained from the Chinese Academy of Sciences’ Cell Bank in Shanghai, China. All cells were cultured at 37 ℃ in a humidified atmosphere containing 5% CO₂. Specifically, MCF-10A cells were cultured in DMEM/F12 medium supplemented with 5% horse serum, 20 ng/mL epidermal growth factor, 0.5 µg/mL hydrocortisone, 10 µg/mL human insulin, 1% non-essential amino acids, and 1% penicillin-streptomycin solution. MCF-7, MDA-MB-231, SK-BR-3, and BT-474 cells were cultured in RPMI-1640 medium. BT549 and MDA-MB-468 cells were cultured in DMEM medium. All media were supplemented with 10% fetal bovine serum (FBS) and 1% penicillin-streptomycin solution, with BT-474 medium additionally containing 10 µg/mL insulin. All cell lines used in the experiments were within 10 passages.

Quantitative real-time polymerase chain reaction

Total RNA was extracted from MCF-10A, SK-BR-3, BT-474, MCF-7, MDA-MB-231, MDA-MB-468, and BT549 cells using the RNA-Quick purification kit. Subsequently, 1 µg of total RNA was reverse transcribed into complementary DNA (cDNA) with the Fast all-in-one RT kit. The expression levels of CACNA1H, CEMIP, CLIC6, CXCL13, NTSR1, S100B, SAA1, and SLC1A1 genes were detected by RT-qPCR using 2× SYBR Green qPCR Master Mix and the Roche Light Cycler 480 QPCR instrument (Roche Diagnostics, Germany). Results were analyzed using the 2^–ΔΔCt method, and the primer sequences are listed in Table 1.

Gene set variation analysis (GSVA)

GSVA was conducted using the “GSVA” package (v1.50.5) (30) with gene sets “c2.cp.kegg.v7.5.1.symbols.gmt” and “h.all.v7.5.1.symbols.gmt” from MSigDB. Differences between high- and low-risk groups were analyzed using the “limma” package (v3.58.1) (31). Pathways with |t-value| >2 and P<0.05 were considered significantly enriched. The top 10 pathways were visualized.

Gene set enrichment analysis (GSEA)

Spearman correlation was computed between each prognostic gene and all other genes using the “psych” package. Genes were ranked by correlation coefficient, and GSEA was performed with “clusterProfiler” (v4.10.1) (32) using the c2.cp.kegg.v7.5.1.symbols.gmt gene set. Significantly enriched pathways were identified with |normalized enrichment score (NES)| >1, false discovery rate (FDR) <0.25, and P<0.05. The top five enriched pathways per gene were visualized using the “msigdbr” package (v7.5.1) (33).

Construction of molecular regulatory networks

Transcription factors (TFs) regulating prognostic genes were predicted using NetworkAnalyst (https://www.networkanalyst.ca/), and TF-mRNA networks were visualized with Cytoscape (v3.10.2) (34).

Upstream micro-RNAs (miRNAs) targeting prognostic genes were predicted using Diana-microT (http://diana.imis.athena-innovation.gr) and Miranda (http://www.microrna.org). Common miRNAs from both databases were used to construct miRNA-mRNA networks via Cytoscape.

Immune microenvironment analysis

The ESTIMATE algorithm was applied to calculate immune, stromal, and ESTIMATE scores. Wilcoxon test was used to compare these scores between high- and low-risk groups (P<0.05).

Immune infiltration analysis

To analyze the differences in immune cell infiltration between the high-risk and low-risk groups of triple-negative breast cancer (TNBC) patients in the TCGA-BRCA training set, the ssGSEA algorithm from the “GSVA” package (v1.50.5) was used to calculate the enrichment scores of 28 immune cell types (35) in each sample. The relative infiltration levels of these immune cells were then compared between the high-risk and low-risk groups using the Wilcoxon rank-sum test, with immune cells showing a significant difference (P<0.05) being defined as differentially infiltrated. To further elucidate the interrelationships among these differentially infiltrated immune cells and their associations with the eight prognostic genes, Spearman correlation analysis was performed using the “psych” package (v2.4.3), with correlations meeting the thresholds of |cor| >0.3 and P<0.05 considered statistically significant. Additionally, the risk score’s correlation with immune cell infiltration was further assessed. All resulting data were visualized using the “ggplot2” package (v3.5.1).

Drug sensitivity analysis

Drug sensitivity data were obtained from the Genomics of Drug Sensitivity in Cancer (GDSC, https://www.cancerrxgene.org/). The “oncoPredict” package (v1.2) was used to compute half-maximal inhibitory concentration (IC₅₀) values for common chemotherapeutic drugs. Wilcoxon test compared IC₅₀ values between risk groups (P<0.05) (36). The top 10 most significantly different drugs were visualized using “ggplot2”.

Statistical analysis

All bioinformatic analyses were performed using R (v4.4.3). A P value <0.05 was considered statistically significant.

Results

Identification of DEGs

Differential expression analysis was performed between BC samples and normal controls in the TCGA-BRCA training set, with screening criteria set at |log₂FC| >2 and adjusted P value (P adj) <0.05. A total of 1,848 DEGs were identified, including 1,121 up-regulated and 727 down-regulated genes in BC samples (Table S2). Compared with normal tissues, genes such as IBSP, CGA, CST4, MUC2, CST5, CSAG1, MAGEA12, CARTPT, MAGEA3, and MAGEA6 were significantly up-regulated in BC, while GLYAT, APOB, MYOC, GPD1, PLIN1, LEP, HEPACAM, CIDEC, AL845331.1, and CA4 were markedly down-regulated (Figure 1A,1B).

Figure 1 Identification of differentially expressed genes in breast cancer. Volcano plot (A) and heatmap (B) of DEGs. DEGs, differentially expressed genes.

Identification of candidate genes and functional enrichment and PPI analysis

To identify differential genes associated with CHRGs in BC, we intersected the 1,848 DEGs obtained above with the 4,006 CHRGs retrieved from MSigDB, yielding 477 overlapping genes (Figure 2A, Table S3), which were defined as candidate genes for subsequent analyses.

Figure 2 Screening and identification of DEGs and CHRGs. (A) Venn diagram analysis showing the overlap between DEGs and CHRGs. (B) Quantitative summary of the Venn diagram: 3,529 genes unique to DEGs (65.6% of total), 477 genes overlapping between DEGs and CHRGs (8.9%), and 1,371 genes unique to CHRGs (25.5%). (C) PPI network visualization of CHRGs. (D) Topological analysis of core gene modules within the PPI network. BP, biological process; CC, cellular component; CHRGs, cation-homeostasis-related genes; DEGs, differentially expressed genes; MF, molecular function; PPI, protein-protein interaction.

GO functional enrichment analysis was conducted on the above 477 candidate genes, and a total of 2,008 pathways were significantly enriched (P<0.05) (Table S4). Among them, in BP, they were mainly enriched in pathways such as regulation of monovalent ion transmembrane transport, calcium ion transport, and calcium ion transmembrane transport; in CC, they were mainly enriched in pathways such as monovalent ion channel complex, transmembrane transporter complex, and cation channel complex transmembrane transporter complex; and in MF, they were involved in pathways such as monovalent ion gated channel activity, gated channel activity, and divalent ion channel activity (Figure 2B). The results of GO enrichment analysis indicated that these candidate genes were widely involved in the BPs of ion channels.

KEGG analysis revealed 72 significantly enriched pathways (P<0.05) (Table S5), including calcium signaling pathway (hsa04020), neuroactive ligand-receptor interaction (hsa04080), and neuroactive ligand signaling (hsa04082) (Figure 2C). A PPI network was constructed with high-confidence interactions (confidence ≥0.9), resulting in 232 nodes and 387 edges. Key hub genes such as GNAI1 interacted with multiple proteins (e.g., FGF2, GRIA2, KCNJ10), highlighting functional connectivity among candidate genes (Figure 2D).

Construction and validation of a BC prognostic model

Univariate Cox regression and PH assumption testing identified 8 prognosis-related genes: NTSR1, CEMIP, CACNA1H, SAA1, CXCL13, CLIC6, SLC1A1, and S100B. Among these, NTSR1, CEMIP, and CACNA1H were risk factors (HR >1), while SAA1, CXCL13, CLIC6, SLC1A1, and S100B were protective (HR <1) (Figure 3A, Table 2). PH assumption testing confirmed the suitability of these genes for Cox modeling (Figure S1). LASSO regression further validated these 8 genes (λ_min =0.0025) as core prognostic markers (Figure 3B).

Figure 3 Construction and validation of a cation channel-based prognostic model. (A) Forest plot of the eight genes significantly associated with overall survival (HR with 95% CI). (B) LASSO coefficient profiles and cross-validation curve for the eight candidate genes. (C) Distribution of risk scores based on the eight-gene signature in the training cohort (n=1,000), with patients stratified into high- and low-risk groups using the optimal cutoff. (D) Kaplan-Meier survival analysis comparing the two risk groups. (E) ROC curves evaluating the model’s performance for predicting 1-, 3-, and 5-year overall survival in the training set. (F-H) Validation in the GSE20685 cohort (n=327): risk score distribution (F), Kaplan-Meier survival analysis (G), and ROC curve analysis (H) based on the optimal cutoff. AUC, area under the curve; CI, confidence interval; HR, hazard ratio; LASSO, least absolute shrinkage and selection operator; OS, overall survival; ROC, receiver operating characteristic; TCGA-BRCA, The Cancer Genome Atlas-breast invasive carcinoma.

Based on these genes, a prognostic model was built. Using an optimal risk-score cutoff (−0.57), 1,000 BC patients were stratified into high-risk (n=660) and low-risk (n=340) groups (Figure 3C). Kaplan-Meier analysis showed significantly worse survival in the high-risk group (P<0.0001) (Figure 3D). Time-dependent ROC analysis demonstrated AUCs >0.6 for 1-, 3-, and 5-year survival, with increasing predictive accuracy over time (Figure 3E). The model was validated in the GSE20685 dataset (cutoff =1.58), consistently separating patients into high-risk (n=222) and low-risk (n=105) groups with significant survival differences (P<0.0001) and AUC >0.6 (Figure 3F-3H). The expression patterns of the prognostic genes in the GSE20685 dataset were consistent with those in the TCGA-BRCA dataset. These results demonstrate that the constructed BC prognostic model possesses excellent predictive potential.

Independent prognostic analysis and clinical predictive value in BC

Correlation analysis between the risk score and clinical characteristics showed that the P values for the risk score (riskScore), gender, age, and N stage (Stage-N) were all greater than 0.05, indicating that these variables satisfied the PH assumption (Figure S2, Table 3). Univariate Cox regression analysis revealed that riskScore, age, Stage-N (I vs. 0), and Stage-N (II vs. 0) had P values less than 0.05, suggesting that these four clinical features possessed significant prognostic predictive value. Moreover, all four indicators had HRs greater than 1, identifying them as risk factors (Figure 4A). Subsequent PH assumption testing confirmed that age, Stage-N, and riskScore met the PH assumption (P>0.05) (Figure S3, Table 4). Multivariate Cox regression analysis further validated these three factors—age, Stage-N, and riskScore—as independent prognostic predictors (Figure 4B). A nomogram constructed based on these independent prognostic factors demonstrated good predictive performance (Figure 4C). DCA at 1, 3, and 5 years indicated that the model incorporating age, Stage-N, and riskScore outperformed other comparative curves, confirming the superior predictive ability of the nomogram (Figure 4D). Furthermore, the AUC values of the ROC curves all exceeded 0.7, further confirming the clinical practicability of the nomogram integrating the riskScore, age, and Stage-N indicators (Figure 4E).

Figure 4 Construction and validation of a prognostic model integrating risk score and clinical characteristics. (A) Univariate Cox regression analysis identified risk score, Age, and Stage-N as significant prognostic factors. (B) Multivariate Cox regression analysis confirmed Age, Stage-N, and riskScore as independent prognostic factors. (C) A clinical prognostic nomogram incorporating the three independent factors for predicting 1-, 3-, and 5-year overall survival probabilities. (D) DCA of the nomogram. (E) ROC curve analysis evaluating the nomogram’s predictive performance for 1-, 3-, and 5-year survival. ***, P<0.001. AUC, area under the curve; CI, confidence interval; DCA, decision curve analysis; OS, overall survival; ROC, receiver operating characteristic; Stage-N, lymph node stage.

Subgroup analysis validates the broad applicability of the prognostic signature

The prognostic value of established clinicopathological factors and the eight-gene signature was evaluated across different patient subgroups. Analysis of traditional factors confirmed that advanced nodal metastasis (N1–2) and older age (>60 years) were significantly associated with poorer OS (P<0.0001), with time-dependent ROC analysis showing AUCs of 0.6, 0.62, and 0.61 for N stage and 0.7, 0.58, and 0.59 for age at 1-, 3-, and 5-year OS, respectively (Figure 5A,5B). The risk model was then applied within these strata. In the N0 cohort, high-risk patients had significantly worse OS than low-risk patients (P<0.01), with AUCs of 0.62, 0.76, and 0.69 for 1-, 3-, and 5-year OS (Figure 5C). This was mirrored in the N1–2 cohort, where high-risk patients also exhibited inferior outcomes (P<0.0001) and AUCs of 0.72, 0.69, and 0.72 (Figure 5D). Similar results were observed across age groups. Among younger patients (≤60 years), the high-risk group demonstrated significantly poorer survival (P<0.0001), with AUCs of 0.77, 0.71, and 0.76 (Figure 5E). In the older cohort (> 60 years), high-risk patients also showed markedly worse OS (P<0.01), yielding AUCs of 0.64, 0.71, and 0.65 (Figure 5F). These findings confirm that the eight-gene signature provides robust and consistent risk stratification across all clinical subsets, underscoring its stability and broad applicability.

Figure 5 Subgroup analysis of the eight-gene prognostic signature in breast cancer. (A,B) Kaplan-Meier survival curves and time-dependent ROC analyses for (A) nodal status (N0 vs. N1–2) and (B) age (≤60 vs. >60 years). (C-F) Evaluation of the eight-gene signature within specific patient subgroups. Kaplan-Meier curves and corresponding time-dependent ROC analyses for 1-, 3-, and 5-year OS are shown for the (C) N0 cohort, (D) N1–2 cohort, (E) younger patient cohort (≤60 years), and (F) older patient cohort (>60 years). Patients were stratified into high- and low-risk groups based on the predefined risk score cutoff. P values were calculated by log-rank test. AUC, area under the curve; N, lymph node; OS, overall survival; ROC, receiver operating characteristic; TCGA-BRCA, The Cancer Genome Atlas-breast invasive carcinoma.

Expression validation, chromosomal localization, and correlation analysis of prognostic genes in BC

Analysis of the expression levels of the eight prognostic genes in the training set revealed significant differences (P<0.001) between the control and BC groups. Specifically, NTSR1, CEMIP, CACNA1H, CXCL13, and SLC1A1 were up-regulated in BC tissues, while CLIC6, SAA1, and S100B were down-regulated compared to the control group (Figure 6A). Subsequent analysis of chromosomal localization showed that these genes are distributed across multiple chromosomes: NTSR1 on chromosome 20, CEMIP on chromosome 15, CACNA1H on chromosome 16, SAA1 on chromosome 11, CXCL13 on chromosome 4, CLIC6 and S100B on chromosome 21, and SLC1A1 on chromosome 9 (Figure 6B), indicating their genetic independence. Correlation analysis among the eight genes demonstrated several significant relationships: a negative correlation between SAA1 and CEMIP (cor =−0.35, P<0.001); a positive correlation between CLIC6 and SLC1A1 (cor =0.33, P<0.001); negative correlations between S100B and both CACNA1H (cor =−0.32, P<0.001) and CEMIP (cor =−0.32, P<0.001); and a strong positive correlation between S100B and SAA1 (cor =0.76, P<0.001) (Figure 6C). These findings suggest potential functional interactions among these prognostic genes. We validated the mRNA expression levels of eight candidate prognostic genes (CACNA1H, CEMIP, CLIC6, CXCL, NTSR1, S100B, SAA1, and SLC1A1). The results demonstrated that compared with the normal mammary epithelial cell line MCF-10A, these genes collectively exhibited significant differential expression patterns across a panel of BC cell lines. Notably, several genes including CACNA1H, CEMIP, CXCL, NTSR1, and SAA1 showed significant upregulation, which collectively underscores the potential of this gene set as prognostic biomarkers for BC (Figure 6D).

Figure 6 Analysis of differential expression, chromosomal localization, and interaction network of prognosis-related genes. (A) Differential expression of eight prognostic genes between BC and control groups. (B) Chromosomal localization of the eight prognostic genes. (C) Correlation analysis of expression among the prognostic genes. (D) RT-qPCR validation of eight core gene expressions in BC and normal breast cells. *, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001 and ns indicates no statistical significance (P≥0.05). BC, breast cancer; RT-qPCR, reverse transcription quantitative polymerase chain reaction.

GSVA and GSEA analysis of prognostic genes

To explore biological function and pathway differences between high- and low-risk groups, GSVA was performed using KEGG gene sets. A total of 186 pathways were significantly enriched, primarily involving N-glycan biosynthesis, other glycan degradation, and O-glycan biosynthesis (Figure 7A, Table S6). GSVA analysis with hallmark gene sets identified 50 enriched pathways, mainly associated with TNF-α signaling via NF-κB, hypoxia, cholesterol homeostasis, and mitotic spindle (Figure 7B, Table S7).

Figure 7 Pathway enrichment and functional analysis of prognostic genes in high- and low-risk groups. (A) GSVA using KEGG gene sets. (B) GSVA using Hallmark gene sets. (C) GSEA showing enriched pathways associated with the eight prognostic genes. GSEA, gene set enrichment analysis; GSVA, gene set variation analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Further investigation of the biological functions and signaling pathways associated with prognostic genes in BC was conducted through GSEA. The results revealed: CACNA1H was enriched in 27 pathways, primarily cytokine-cytokine receptor interaction and ribosome pathways; CEMIP was enriched in 45 pathways, mainly ribosome and extracellular matrix (ECM)-receptor interaction pathways; CLIC6 was enriched in 32 pathways, predominantly Parkinson’s disease and oxidative phosphorylation pathways; CXCL13 was enriched in 45 pathways, chiefly cytokine-cytokine receptor interaction and chemokine signaling pathways; NTSR1 was enriched in 12 pathways, primarily cell cycle and DNA replication pathways; S100B was enriched in 63 pathways, mainly cytokine-cytokine receptor interaction and hematopoietic cell lineage pathways; SAA1 was enriched in 68 pathways, primarily cytokine-cytokine receptor interaction and hematopoietic cell lineage pathways; SLC1A1 was enriched in 17 pathways, predominantly Parkinson’s disease and oxidative phosphorylation pathways (Figure 7C, Table S8).

Regulatory network analysis of prognostic genes

To identify upstream regulatory factors of the prognostic genes, a TF-mRNA regulatory network was constructed. The results showed that: CEMIP was associated with 32 TFs (e.g., TCF21, AHR); NTSR1 with 37 TFs (e.g., TET1, ATF3); S100B with 13 TFs (e.g., FLI1, HOXC9); CLIC6 with 36 TFs (e.g., ERG, VDR); SLC1A1 with 29 TFs (e.g., AF4, SOX11); CACNA1H with 31 TFs (e.g., CTCF, TAF2); SAA1 with 24 TFs (e.g., CBP, SRY); CXCL13 with 11 TFs (e.g., OCT1, CHD1) (Figure 8A). Furthermore, a miRNA-mRNA regulatory network was established to explore upstream miRNAs of these prognostic genes in BC: NTSR1 was linked to 40 miRNAs (e.g., hsa-miR-92a-2-5p,hsa-miR-185-3p); S100B to 8 miRNAs (e.g., hsa-miR-1287-5p, hsa-miR-3133); CLIC6 to 22 miRNAs (e.g., hsa-miR-1253, hsa-miR-3155a); SLC1A1 to 52 miRNAs (e.g., hsa-miR-96-5p, hsa-miR-130a-5p); CACNA1H and SAA1 were each associated with one miRNA (hsa-miR-32-5p and hsa-miR-660-5p, respectively); CXCL13 was connected to 27 miRNAs (e.g., hsa-miR-1304-5p, hsa-miR-186-5p); no miRNAs were predicted for CEMIP (Figure 8B). These molecular regulatory networks help clarify potential upstream regulatory mechanisms of the prognostic genes in BC.

Figure 8 Analysis of upstream regulatory networks for the prognostic genes. (A) TF-mRNA regulatory network showing transcription factors potentially regulating the eight prognostic genes. (B) miRNA-mRNA regulatory network depicting potential miRNA interactions with the prognostic genes. miRNA, microRNA; TF, transcription factor.

Immune mechanism and drug sensitivity analysis

To investigate differences in immune molecules between the high- and low-risk groups, the tumor microenvironment (TME) score—including Immune Score, Stromal Score, and ESTIMATE Score—was calculated using the ESTIMATE algorithm. The results showed that both the Immune Score and ESTIMATE Score exhibited P values <0.05, indicating significant differences in the immune microenvironment between the two risk groups (Figure 9A).

Figure 9 Characterization of the tumor immune microenvironment and clinical therapeutic implications in high- and low-risk groups. (A) Tumor microenvironment scores. (B) Heatmap of immune cell infiltration levels. (C) Differential infiltration of immune cell types. (D) Correlation network of differentially infiltrated immune cells. (E) Correlation between prognostic genes and immune cells. (F) Correlation between risk score and immune cells. (G) Top candidate drugs with differential sensitivity. *, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001. Cor, correlation coefficient; ESTIMATE, estimation of stromal and immune cells in malignant tumor tissues using expression data; IC₅₀, half maximal inhibitory concentration; MDSC, myeloid-derived suppressor cell; ssGSEA, single-sample gene set enrichment analysis.

Further analysis of immune cell infiltration revealed differential levels of 28 immune cell types between the groups. Activated B cells, activated CD4⁺ T cells, activated CD8⁺ T cells, memory B cells, and immature B cells showed lower infiltration abundance in the high-risk group and relatively higher abundance in the low-risk group (Figure 9B). Significant differences were observed in 25 immune cell subsets between the groups, including activated B cells, macrophages, activated CD8⁺ T cells, activated dendritic cells, T follicular helper cells, type 1 T helper cells, and myeloid-derived suppressor cells (MDSCs) (P<0.05) (Figure 9C). Correlation analysis among differentially infiltrated immune cells showed that MDSCs were strongly positively correlated with effector memory CD8⁺ T cells (cor =0.88, P<0.001), activated dendritic cells (cor =0.83, P<0.001), and immature B cells (cor =0.80, P<0.001) (Figure 9D). Analysis of correlations between prognostic genes and differential immune cells revealed that: S100B was significantly positively correlated with activated B cells (cor =0.56), activated CD8⁺ T cells (cor =0.41), type 1 T helper cells (cor =0.52), T follicular helper cells (cor =0.43), and plasmacytoid dendritic cells (cor =0.51) (all P<0.05); CXCL13 showed strong positive correlations with activated B cells (cor =0.68), activated CD4⁺ T cells (cor =0.62), and activated CD8⁺ T cells (cor =0.66); CACNA1H was significantly negatively correlated with activated B cells (cor =−0.24), activated CD8⁺ T cells (cor =−0.18), type 1 T helper cells (cor =−0.21), T follicular helper cells (cor =−0.16), and plasmacytoid dendritic cells (cor =−0.08) (all P<0.05) (Figure 9E). In addition, correlation analysis between the risk score and immune cell abundance revealed that the risk score was significantly negatively correlated with plasmacytoid dendritic cells (cor =−0.3) and activated B cells (cor =−0.38) (P<0.05) (Figure 9F).

Drug sensitivity analysis identified 192 drugs with significantly different IC₅₀ values between the risk groups (P<0.05) (Table S9). The top 10 most significant agents included BDP.00009066, camptothecin, and CDK9 inhibitors (Figure 9G), providing a candidate basis for subsequent targeted intervention studies in BC.

Discussion

BC, the most common malignant tumor among women worldwide, is closely associated with metabolic reprogramming (37). Disruption of cation homeostasis may further promote tumor progression (38). Recent studies have revealed that metal ion homeostasis (such as copper, zinc, and calcium) and pH regulation collaboratively influence BC metastasis, colonization, microenvironment adaptation, and stemness by modulating the activity and expression of key proteins, providing a theoretical basis for combination therapies targeting ion metabolism (39-42). Based on transcriptomic data, this study identified eight cation homeostasis-related prognostic genes in BC: NTSR1, CEMIP, CACNA1H, SAA1, CXCL13, CLIC6, SLC1A1, and S100B. Further analyses were conducted to explore the biological pathways, immune microenvironment infiltration, and drug sensitivity associated with these genes. The findings shed light on the potential mechanisms of cation homeostasis in BC and offer new insights for its diagnosis and treatment.

This study found that neurotensin receptor 1 (NTSR1) is highly expressed in BC and is a risk factor, which is consistent with the conclusions of previous studies (43). Its activation not only promotes malignant phenotypes through classical signaling pathways such as PI3K/Akt (44), MAPK/ERK (45), and Src (46), but, more critically, functions as a Gq-protein-coupled receptor to trigger calcium release from the endoplasmic reticulum stores via the Gq-PLCβ-IP₃ pathway. This leads to a sharp increase in cytosolic Ca²⁺ concentration. The initial calcium signal can further activate non-selective cation channels such as TRPC, resulting in Na⁺ influx and membrane depolarization, establishing a positive feedback loop and indirectly modulating other ion channels, thereby collectively driving cell migration, invasion, and drug resistance (47,48). The significant upregulation of CACNA1H (encoding the T-type calcium channel Cav3.2) in BC tissues provides a key voltage-dependent calcium influx pathway for tumor cells (49). This study confirmed that the sustained calcium signaling brought about by its high expression can directly activate downstream pro-survival and proliferative signaling axes such as PI3K/Akt, thereby enhancing cell proliferation and epithelial-mesenchymal transition (EMT) capabilities (50,51). This mechanism directly explains the poor prognosis of patients with high expression of CACNA1H, and the anti-tumor effect demonstrated by its inhibitors in preclinical studies further confirms that it can serve as a potential intervention target for high-risk patient groups. S100B is lowly expressed in BC overall (48). The S100 calcium-binding protein beta subunit (S100B), as a protective factor, is generally lowly expressed in BC, and its functional loss is associated with poor prognosis of patients (52). In the nervous system, S100B is a key calcium signal regulator; In BC cells, the loss of its expression may lead to impaired calcium sensing ability. The failure of this “calcium sensing” function through the EF-hand domain may prevent cells from precisely regulating calcium-dependent downstream processes (such as inhibiting EMT) (53), thereby indirectly removing the constraints on the malignant evolution of tumors. This, together with the “calcium release” mediated by NTSR1 and the “calcium influx” mediated by CACNA1H, constitutes a complete and functionally complementary calcium signaling network, systematically driving tumor development.

The cell migration-inducing hyaluronan-binding protein (CEMIP), as a risk gene, is upregulated in BC. It not only reshapes the ECM, but recent evidence suggests that it can also interact with the endoplasmic reticulum chaperone protein BiP. It directly induces calcium leakage from the endoplasmic reticulum into the cytoplasm, activates calcium-sensitive signals such as PKCα, and thereby closely links the two pro-cancer processes of microenvironment remodeling and intracellular calcium homeostasis imbalance (54). The protective gene chloride intracellular channel 6 (CLIC6) in this model was significantly downregulated in BC. The chloride channel protein encoded by this gene is normally located in the intracellular intimal system. The loss of its expression may disrupt the chloride ion flow, thereby disrupting cell volume regulation and membrane potential homeostasis (55). A stable membrane potential is the basis for maintaining the normal gating of voltage-dependent calcium channels. The membrane potential disorder caused by the down-regulation of CLIC6 may indirectly interfere with the calcium-signal-dependent cancer-promoting program driven by risk genes such as CEMIP, thereby exerting its protective effect. The successful analysis of this network not only explains the structural basis of the model’s outstanding predictive ability at the system level, but also provides a novel theoretical basis for implementing multi-target ion intervention strategies for high-risk patients.

Our research further reveals that the influence of calcium signaling extends beyond the tumor cells themselves to the tumor immune microenvironment. For instance, the chemokine C-X-C motif ligand 13 (CXCL13), a key factor in B cell recruitment, acts as a ligand for CXCR5 and directly induces calcium mobilization in immune cells (56-58). Although CXCL13 is highly expressed in BC (59) both our study and recent clinical evidence indicate its association with improved patient prognosis (HR <1) (60). This suggests that CXCL13 may exert a complex immunomodulatory role within the TME through a dual mechanism: orchestrating immune cell chemotaxis and fine-tuning local calcium signaling. Although the direct effect of serum amyloid A (SAA1) on calcium ions remains unclear, its role as an acute-phase protein connecting systemic inflammation and cancer (61), along with its complex expression patterns and functional paradoxes, suggests it may be a potential regulatory node. Finally, carrier family 1 gene (SLC1A1) through its electrogenic transport mechanism of “symport of 3 Na⁺ and 1 H⁺ coupled with antiport of 1 K⁺” (62,63), directly maintains the Na⁺/K⁺ transmembrane gradient, which is fundamental to cellular electrophysiological balance, including calcium homeostasis.

In summary, the eight genes identified in this study may collectively participate in the regulation of cation homeostasis through direct or indirect mechanisms. Among them, NTSR1 and CACNA1H constitute key nodes of calcium signal input by mediating endoplasmic reticulum calcium release and voltage-dependent calcium influx, respectively. S100B, functioning as a calcium sensor protein, may have its downregulated expression attenuate the cell’s precise response to calcium signals. CEMIP links microenvironment remodeling with alterations in intracellular calcium homeostasis by inducing calcium leakage from the endoplasmic reticulum. CLIC6 and SLC1A1 provide the electrophysiological foundation for cation homeostasis by maintaining chloride ion homeostasis and the sodium/potassium ion gradient, respectively. CXCL13 connects cation signaling with the immune microenvironment by inducing calcium mobilization in immune cells. Meanwhile, SAA1 may act as a potential node linking inflammation and ion regulation, exerting an indirect effect. In conclusion, these eight genes form functionally complementary synergistic relationships across multiple dimensions—including cation membrane transport, intracellular release, signal sensing, electrophysiological maintenance, and microenvironmental interaction—collectively contributing to the malignant progression of BC.

The enrichment of the high-risk group in N-glycan and O-glycan pathways suggests that glycosylation modifications serve as a critical initiating step for malignant phenotypes. Aberrant glycosylation drives malignant progression through multiple avenues. It directly promotes tumor proliferation and metastasis by modifying key receptors like EpCAM, enhancing their stability and signaling efficiency (64). More notably, it can disrupt membrane potential homeostasis by altering the conformation and function of ion channel proteins. This aberrant glycosylation establishes the initial conditions for subsequent cation homeostasis imbalance, representing an upstream event in the pathogenic network (65). This creates the initial conditions for subsequent cation homeostasis imbalance and signaling pathway abnormalities, constituting an upstream event in the network. The significant enrichment of the cytokine-cytokine receptor interaction pathway places the key biomarkers we identified (e.g., CACNA1H, CXCL13, S100B, SAA1) at the core of the TME. This pathway acts as a powerful signaling hub: on one hand, it directly promotes proliferation, angiogenesis, and immune evasion by activating downstream pathways such as JAK/STAT and NF-κB through factors like TNF-α (66); on the other hand, it regulates ionic events including intracellular calcium signaling (67), serving as a bridge connecting immune responses and cellular electrophysiological activities. Particularly crucial is the activation of the TNF-α/NF-κB pathway within this network, which not only constitutes the core of inflammatory responses but also acts as a key driver of metastasis by inducing EMT (68). To support their rapid proliferation and invasion, high-risk tumor cells exhibit significant metabolic remodeling. Hypoxia signaling drives glycolytic and lipid metabolic reprogramming through HIF-1α, meeting the bioenergetic and biosynthetic demands of rapidly proliferating cells (69). Concurrently, activation of the mitotic spindle pathway directly reflects abnormally vigorous cell division activity, potentially associated with enhanced genomic instability (70), together constituting the terminal execution phase of the malignant phenotype.

During both the initiation and effector phases of adaptive immunity, calcium influx serves as a “common switch” for T and B cell activation. Particularly for CD8⁺ T cells, clonal proliferation, production of effector cytokines, and differentiation into cytotoxic T cells all depend on sustained calcium influx triggered by T cell receptor activation (71). Furthermore, calcium signaling plays a dominant role in functional polarization of macrophages. The balance between pro-inflammatory M1 (anti-tumor) and anti-inflammatory M2 (pro-tumor) tumor-associated macrophages is crucial in determining the nature of the immune microenvironment (72). Current research clearly demonstrates that the polarization process of macrophages is finely regulated by Ca²⁺ signaling (73). Therefore, disturbances in calcium homeostasis within the TME likely promote an immunosuppressive microenvironment by preferentially driving macrophage polarization toward the M2 phenotype, thereby facilitating tumor immune escape. This study further revealed that the risk score was significantly negatively correlated with the abundance of plasmacytoid dendritic cells and activated B cells. Combined with the observed reduction in CD8⁺ T cell infiltration and the expansion of MDSCs in the high-risk group, this suggests that cation homeostasis dysregulation may impair the functional maturation of antigen-presenting cells such as plasmacytoid dendritic cells and B cells (74,75). This impairment weakens their ability to effectively activate CD8⁺ T cells, thereby compromising adaptive anti-tumor immune responses and creating favorable conditions for MDSC-mediated innate immune suppression (76). This remodeling of the immune cell lineage driven by cation signaling disruption may represent a key mechanism underlying the poor prognosis observed in high-risk group patients. Cation homeostasis dysregulation driven by key genes (such as NTSR1 and CACNA1H) in BC cells not only directly promotes malignant behavior of tumor cells but also remotely reshapes the tumor immune microenvironment by “hijacking” the calcium ion-centered public signaling system.

It is particularly important to emphasize that cation homeostasis imbalance likely represents the common mechanism and “information highway” permeating all the aforementioned pathways. It functions both as a downstream effector of upstream glycosylation modifications and cytokine signaling, and as a core regulator that reciprocally controls tumor proliferation, cell death, invasion, and microenvironment homeostasis by influencing signal transduction, gene expression, and immune cell function.

Our drug sensitivity analysis revealed that high-risk patients, characterized by significant cation homeostasis imbalance and immunometabolic reprogramming in our model, demonstrated heightened sensitivity to multiple pharmacological agents (including BDP.00009066, camptothecin, and CDK9 inhibitors). This finding carries profound mechanistic implications: the sensitivity to topoisomerase I inhibitor camptothecin may be directly linked to the active spindle assembly pathway and enhanced genomic instability observed in the high-risk group, rendering their rapidly dividing cells more vulnerable to DNA damage (77). Furthermore, the efficacy of CDK9 inhibitors aligns with the aberrant cell cycle progression indicated by the model, providing a promising approach for targeting transcriptional dependency in these aggressive tumors (78). In summary, these computationally predicted drug sensitivities not only lay the foundation for subsequent targeted intervention studies in high-risk BC patients but also indicate the potential clinical translational value of the prognostic model.

Based on the immune exhaustion characteristics presented in the high-risk group, such as reduced CD8⁺ T cell infiltration and MDSC expansion, the theoretical basis for combination therapy can be further explored. Existing studies have shown that cation homeostasis imbalance, particularly calcium signaling disturbances driven by key genes such as NTSR1 and CACNA1H, may not only directly promote malignant tumor proliferation but also participate in the formation of an immunosuppressive microenvironment by interfering with the calcium-calcineurin-NFAT signaling pathway essential for T cell activation (16,79,80). Therefore, in theory, the combined application of ion channel inhibitors and immune checkpoint inhibitors may produce synergistic effects: the former directly inhibits tumor cell activity by regulating ion homeostasis and indirectly ameliorates the immunosuppressive state, while the latter aims to restore the effector function of exhausted T cells (81,82). This interventional approach, which combines the intrinsic driving mechanism of tumor cells with the extrinsic immune escape pathway, may provide a new research direction for improving the efficacy of treatment in high-risk BC patients. In conclusion, based on prognostic stratification, the model constructed in this study provides theoretical clues and candidate biomarkers for subsequent exploration of combination therapeutic regimens targeting cation homeostasis dysregulation and immune checkpoints.

It should be emphasized that although this study reveals the potential value of targeting ion channels in BC treatment, currently, ion channel inhibitors approved by the U.S. Food and Drug Administration (FDA) are mainly used in the treatment of cardiovascular and neurological diseases such as arrhythmias, hypertension, and epilepsy. Their application in oncology is mostly still in preclinical research or early clinical trial stages, and there is still a considerable gap before true clinical translation (83). Therefore, the clinical significance of the prognostic model based on cation homeostasis-related genes proposed in this study is currently mainly reflected in providing a new tool for the precise identification and stratification of high-risk patients, rather than directly suggesting the immediate use of existing ion channel inhibitors. This study provides a theoretical foundation and candidate biomarkers for the future development of novel cancer treatment strategies targeting specific ion channels, but its translation into clinical practice still requires substantial subsequent research for validation.

Conclusions

The BC prognostic model constructed in this study reveals a core biological paradigm: cation homeostasis imbalance driven by key genes (such as NTSR1 and CACNA1H) in tumor cells, particularly calcium signaling dysregulation, constitutes a common pathway connecting “cell-intrinsic malignancy” and “immune microenvironment remodeling”. On one hand, it directly promotes tumor proliferation, migration, and drug resistance by activating signaling axes such as PI3K/Akt; on the other hand, it systemically suppresses anti-tumor immunity by hijacking calcium signaling (e.g., the calcineurin-NFAT pathway) essential for immune cell activation (including T cells, B cells, and macrophages). Consequently, the high predictive power of this model stems from its simultaneous quantification of both the autonomous invasive potential of tumors and their ability to mediate immune escape. This provides a novel theoretical foundation for developing combination therapies that co-target ion channels and immune checkpoints.

The main limitation of this study lies in its heavy reliance on bioinformatic analyses without experimental validation. Although the systematic screening of cation homeostasis-related prognostic genes and construction of the predictive model were achieved through computational biology methods, the research remains at the data analysis level. The interpretation of the cation homeostasis regulatory network is largely based on database predictions, lacking functional experiments (such as in vitro cell assays or animal models) to confirm that these genes indeed regulate ion balance. Subsequent experimental studies are required to further validate the actual biological significance of these bioinformatic findings.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1-2814/rc

Data Sharing Statement: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1-2814/dss

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1-2814/prf

Funding: This study was supported by the National Natural Science Foundation of China (Grant No. 82404685), Zhejiang Science and Technology Department “vanguard” “leading goose” Research (grant No. 2023C03044), and the Suzhou Youth Science and Technology Program (KJXW2021067).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1-2814/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Smolarz B, Nowak AZ, Romanowicz H. Breast Cancer-Epidemiology, Classification, Pathogenesis and Treatment (Review of Literature). Cancers (Basel) 2022;14:2569. [Crossref] [PubMed]
2. Kim J, Harper A, McCormack V, et al. Global patterns and trends in breast cancer incidence and mortality across 185 countries. Nat Med 2025;31:1154-62. [Crossref] [PubMed]
3. Yan P, Jimenez ER, Li Z, et al. Midkine as a driver of age-related changes and increase in mammary tumorigenesis. Cancer Cell 2024;42:1936-1954.e9. [Crossref] [PubMed]
4. Zuckerman JE, Hsu YL, Grimes WJ. Expression of c-fos and c-myc in Saccharomyces cerevisiae. J Mol Biol 1991;219:253-61.
5. Mardamshina M, Karagach S, Mohan V, et al. Integrated spatial proteomic analysis of breast cancer heterogeneity unravels cancer cell phenotypic plasticity. Nat Commun 2025;16:10482. [Crossref] [PubMed]
6. Waks AG, Winer EP. Breast Cancer Treatment: A Review. JAMA 2019;321:288-300. [Crossref] [PubMed]
7. Li Y, Zhang H, Merkher Y, et al. Recent advances in therapeutic strategies for triple-negative breast cancer. J Hematol Oncol 2022;15:121. [Crossref] [PubMed]
8. Garcia-Martinez L, Zhang Y, Nakata Y, et al. Epigenetic mechanisms in breast cancer therapy and resistance. Nat Commun 2021;12:1786. [Crossref] [PubMed]
9. Iwadate Y, Golubeva YA, Slauch JM. Cation Homeostasis: Coordinate Regulation of Polyamine and Magnesium Levels in Salmonella. mBio 2023;14:e0269822. [Crossref] [PubMed]
10. Silver BB, Zhang SX, Rabie EM, et al. Substratum stiffness tunes membrane voltage in mammary epithelial cells. J Cell Sci 2021;134:jcs256313. [Crossref] [PubMed]
11. Xu X, Li N, Wang Y, et al. Calcium channel TRPV6 promotes breast cancer metastasis by NFATC2IP. Cancer Lett 2021;519:150-60. [Crossref] [PubMed]
12. Meng X, Cai C, Wu J, et al. TRPM7 mediates breast cancer cell migration and invasion through the MAPK pathway. Cancer Lett 2013;333:96-102. [Crossref] [PubMed]
13. Zeng J, Tian D, Zhang J, et al. GNAL-driven calcium signaling reshapes the spatiotemporal immune landscape in ER(+) breast cancer: causal insights and prognostic implications. Transl Cancer Res 2026;15:24. [Crossref] [PubMed]
14. Drew D, Boudker O. Ion and lipid orchestration of secondary active transport. Nature 2024;626:963-74. [Crossref] [PubMed]
15. Sun L, Zhang H, Gao P. Metabolic reprogramming and epigenetic modifications on the path to cancer. Protein Cell 2022;13:877-919. [Crossref] [PubMed]
16. Zheng S, Wang X, Zhao D, et al. Calcium homeostasis and cancer: insights from endoplasmic reticulum-centered organelle communications. Trends Cell Biol 2023;33:312-23. [Crossref] [PubMed]
17. Wang Y, Zhuang H, Jiang XH, et al. Unveiling the key genes, environmental toxins, and drug exposures in modulating the severity of ulcerative colitis: a comprehensive analysis. Front Immunol 2023;14:1162458. [Crossref] [PubMed]
18. Wang J, Wu N, Feng X, et al. PROS1 shapes the immune-suppressive tumor microenvironment and predicts poor prognosis in glioma. Front Immunol 2022;13:1052692. [Crossref] [PubMed]
19. Zhang X, Chao P, Zhang L, et al. Single-cell RNA and transcriptome sequencing profiles identify immune-associated key genes in the development of diabetic kidney disease. Front Immunol 2023;14:1030198. [Crossref] [PubMed]
20. Zhou W, Li H, Zhang J, et al. Identification and mechanism analysis of biomarkers related to butyrate metabolism in COVID-19 patients. Ann Med 2025;57:2477301. [Crossref] [PubMed]
21. Lo HJ, Tsai CH, Huang TW. Apoptosis-associated genetic mechanisms in the transition from rheumatoid arthritis to osteoporosis: A bioinformatics and functional analysis approach. APL Bioeng 2024;8:046107. [Crossref] [PubMed]
22. Song DM, Shen T, Feng K, et al. LIG1 is a novel marker for bladder cancer prognosis: evidence based on experimental studies, machine learning and single-cell sequencing. Front Immunol 2024;15:1419126. [Crossref] [PubMed]
23. Li CL, Yao ZY, Qu C, et al. Machine learning model reveals the risk, prognosis, and drug response of histamine-related signatures in pancreatic cancer. Discov Oncol 2025;16:155. [Crossref] [PubMed]
24. Shi H, Yuan X, Liu G, et al. Identifying and Validating GSTM5 as an Immunogenic Gene in Diabetic Foot Ulcer Using Bioinformatics and Machine Learning. J Inflamm Res 2023;16:6241-56. [Crossref] [PubMed]
25. Robin X, Turck N, Hainard A, et al. pROC: an open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinformatics 2011;12:77. [Crossref] [PubMed]
26. Li M, Wei X, Zhang SS, et al. Recognition of refractory Mycoplasma pneumoniae pneumonia among Myocoplasma pneumoniae pneumonia in hospitalized children: development and validation of a predictive nomogram model. BMC Pulm Med 2023;23:383. [Crossref] [PubMed]
27. Zhang Z, Cortese G, Combescure C, et al. Overview of model validation for survival regression model with competing risks using melanoma study data. Ann Transl Med 2018;6:325. [Crossref] [PubMed]
28. Hu Y, Yan C, Hsu CH, et al. OmicCircos: A Simple-to-Use R Package for the Circular Visualization of Multidimensional Omics Data. Cancer Inform 2014;13:13-20. [Crossref] [PubMed]
29. Goncharova IA, Nazarenko MS, Babushkina NP, et al. Genetic Predisposition to Early Myocardial Infarction. Mol Biol (Mosk) 2020;54:224-32. [Crossref] [PubMed]
30. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 2013;14:7. [Crossref] [PubMed]
31. Xu M, Zhou H, Hu P, et al. Identification and validation of immune and oxidative stress-related diagnostic markers for diabetic nephropathy by WGCNA and machine learning. Front Immunol 2023;14:1084531. [Crossref] [PubMed]
32. Wu T, Hu E, Xu S, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb) 2021;2:100141. [Crossref] [PubMed]
33. Gao H, Sun Z, Hu X, et al. Identification of glycolysis-related gene signatures for prognosis and therapeutic targeting in idiopathic pulmonary fibrosis. Front Pharmacol 2025;16:1486357. [Crossref] [PubMed]
34. Xiao C, Zhang X, Hou B, et al. CYFIP2: potential pancreatic cancer biomarker and immunotherapeutic target. Discov Oncol 2024;15:847. [Crossref] [PubMed]
35. Rattu P, Glencross F, Mader SL, et al. Corrigendum to “Atomistic level characterisation of ssDNA translocation through the E. coli proteins CsgG and CsgF for nanopore sequencing” Comput Struct Biotechnol J 2022;20:1027. [Comput Struct. Biotechnol. J. 19 (2021) 6417–6430].
36. Maeser D, Gruener RF, Huang RS. oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief Bioinform 2021;22:bbab260. [Crossref] [PubMed]
37. Lei P, Wang W, Sheldon M, et al. Role of Glucose Metabolic Reprogramming in Breast Cancer Progression and Drug Resistance. Cancers (Basel) 2023;15:3390. [Crossref] [PubMed]
38. Zamay TN, Zamay SS, Zamay GS, et al. Systemic Mechanisms of Ionic Regulation in Carcinogenesis. Cancers (Basel) 2025;17:286. [Crossref] [PubMed]
39. Zhang X, Jiang Q, Su Y, et al. AMPK phosphorylates and stabilises copper transporter 1 to synergise metformin and copper chelator for breast cancer therapy. Br J Cancer 2023;128:1452-65. [Crossref] [PubMed]
40. Vogel-González M, Musa-Afaneh D, Rivera Gil P, et al. Zinc Favors Triple-Negative Breast Cancer's Microenvironment Modulation and Cell Plasticity. Int J Mol Sci 2021;22:9188. [Crossref] [PubMed]
41. Amith SR, Fliegel L. Regulation of the Na+/H+ Exchanger (NHE1) in Breast Cancer Metastasis. Cancer Res 2013;73:1259-64. [Crossref] [PubMed]
42. Makena MR, Ko M, Mekile AX, et al. Secretory pathway Ca(2+)-ATPase SPCA2 regulates mitochondrial respiration and DNA damage response through store-independent calcium entry. Redox Biol 2022;50:102240. [Crossref] [PubMed]
43. Dupouy S, Viardot-Foucault V, Alifano M, et al. The neurotensin receptor-1 pathway contributes to human ductal breast cancer progression. PLoS One 2009;4:e4223. [Crossref] [PubMed]
44. Hung YH, Wang HC, Hsu SH, et al. Neuron-derived neurotensin promotes pancreatic cancer invasiveness and gemcitabine resistance via the NTSR1/Akt pathway. Am J Cancer Res 2024;14:448-66. [Crossref] [PubMed]
45. Qin R, Fan X, Huang Y, et al. G protein-coupled receptors: pivotal hubs in gastric cancer malignancy-from multidimensional crosstalk to precision therapeutics. J Transl Med 2025;23:879. [Crossref] [PubMed]
46. Moody TW, Ramos-Alvarez I, Jensen RT. Peptide G-Protein-Coupled Receptors and ErbB Receptor Tyrosine Kinases in Cancer. Biology (Basel) 2023;12:957. [Crossref] [PubMed]
47. Wang H, Cheng X, Tian J, et al. TRPC channels: Structure, function, regulation and recent advances in small molecular probes. Pharmacol Ther 2020;209:107497. [Crossref] [PubMed]
48. Ehrlich AT, Couvineau P, Schamiloglu S, et al. Visualization of real-time receptor endocytosis in dopamine neurons enabled by NTSR1-Venus knock-in mice. Front Cell Neurosci 2022;16:1076599. [Crossref] [PubMed]
49. Daniil G, Fernandes-Rosa FL, Chemin J, et al. CACNA1H Mutations Are Associated With Different Forms of Primary Aldosteronism. EBioMedicine 2016;13:225-36. [Crossref] [PubMed]
50. Kulkarni GC, Saha R, Peters CJ. Ion channel expression and function in glioblastoma multiforme (GBM): pathophysiological mechanisms and therapeutic potential. Biochim Biophys Acta Mol Cell Res 2025;1872:119982. [Crossref] [PubMed]
51. Fattahi S, Amjadi-Moheb F, Tabaripour R, et al. PI3K/AKT/mTOR signaling in gastric cancer: Epigenetics and beyond. Life Sci 2020;262:118513. [Crossref] [PubMed]
52. Yen MC, Huang YC, Kan JY, et al. S100B expression in breast cancer as a predictive marker for cancer metastasis. Int J Oncol 2018;52:433-40. [Crossref] [PubMed]
53. Uemura T, Green M, Warsh JJ. Chronic LiCl pretreatment suppresses thrombin-stimulated intracellular calcium mobilization through TRPC3 in astroglioma cells. Bipolar Disord 2016;18:549-62. [Crossref] [PubMed]
54. Domanegg K, Sleeman JP, Schmaus A. CEMIP, a Promising Biomarker That Promotes the Progression and Metastasis of Colorectal and Other Types of Cancer. Cancers (Basel) 2022;14:5093. [Crossref] [PubMed]
55. Loyo-Celis V, Patel D, Sanghvi S, et al. Biophysical characterization of chloride intracellular channel 6 (CLIC6). J Biol Chem 2023;299:105349. [Crossref] [PubMed]
56. Legler DF, Loetscher M, Roos RS, et al. B cell-attracting chemokine 1, a human CXC chemokine expressed in lymphoid tissues, selectively attracts B lymphocytes via BLR1/CXCR5. J Exp Med 1998;187:655-60. [Crossref] [PubMed]
57. Harrer C, Otto F, Pilz G, et al. The CXCL13/CXCR5-chemokine axis in neuroinflammation: evidence of CXCR5+CD4 T cell recruitment to CSF. Fluids Barriers CNS 2021;18:40. [Crossref] [PubMed]
58. Hussain M, Adah D, Tariq M, et al. CXCL13/CXCR5 signaling axis in cancer. Life Sci 2019;227:175-86. [Crossref] [PubMed]
59. Chen L, Huang Z, Yao G, et al. The expression of CXCL13 and its relation to unfavorable clinical characteristics in young breast cancer. J Transl Med 2015;13:168. [Crossref] [PubMed]
60. Ding Y, Zhou Q, Ding B, et al. Transcriptome analysis reveals the clinical significance of CXCL13 in Pan-Gyn tumors. J Cancer Res Clin Oncol 2024;150:116. [Crossref] [PubMed]
61. Gaiser AK, Bauer S, Ruez S, et al. Serum Amyloid A1 Induces Classically Activated Macrophages: A Role for Enhanced Fibril Formation. Front Immunol 2021;12:691155. [Crossref] [PubMed]
62. Kanai Y, Clémençon B, Simonin A, et al. The SLC1 high-affinity glutamate and neutral amino acid transporter family. Mol Aspects Med 2013;34:108-20. [Crossref] [PubMed]
63. Guo W, Li K, Sun B, et al. Dysregulated Glutamate Transporter SLC1A1 Propels Cystine Uptake via Xc(-) for Glutathione Synthesis in Lung Cancer. Cancer Res 2021;81:552-66. [Crossref] [PubMed]
64. Ščupáková K, Adelaja OT, Balluff B, et al. Clinical importance of high-mannose, fucosylated, and complex N-glycans in breast cancer metastasis. JCI Insight 2021;6:e146945. [Crossref] [PubMed]
65. Issa FA, Hall MK, Hatchett CJ, et al. Compromised N-Glycosylation Processing of Kv3.1b Correlates with Perturbed Motor Neuron Structure and Locomotor Activity. Biology (Basel) 2021;10:486. [Crossref] [PubMed]
66. Habanjar O, Bingula R, Decombat C, et al. Crosstalk of Inflammatory Cytokines within the Breast Tumor Microenvironment. Int J Mol Sci 2023;24:4002. [Crossref] [PubMed]
67. Hasiakos S, Gwack Y, Kang M, et al. Calcium Signaling in T Cells and Chronic Inflammatory Disorders of the Oral Cavity. J Dent Res 2021;100:693-9. [Crossref] [PubMed]
68. Pires BR, Mencalha AL, Ferreira GM, et al. NF-kappaB Is Involved in the Regulation of EMT Genes in Breast Cancer Cells. PLoS One 2017;12:e0169622. [Crossref] [PubMed]
69. Infantino V, Santarsiero A, Convertini P, et al. Cancer Cell Metabolism in Hypoxia: Role of HIF-1 as Key Regulator and Therapeutic Target. Int J Mol Sci 2021;22:5703. [Crossref] [PubMed]
70. Hosea R, Hillary S, Naqvi S, et al. The two sides of chromosomal instability: drivers and brakes in cancer. Signal Transduct Target Ther 2024;9:75. [Crossref] [PubMed]
71. Park YJ, Yoo SA, Kim M, et al. The Role of Calcium-Calcineurin-NFAT Signaling Pathway in Health and Autoimmune Diseases. Front Immunol 2020;11:195. [Crossref] [PubMed]
72. Qiu Y, Chen T, Hu R, et al. Next frontier in tumor immunotherapy: macrophage-mediated immune evasion. Biomark Res 2021;9:72. [Crossref] [PubMed]
73. Nascimento Da Conceicao V, Sun Y, Ramachandran K, et al. Resolving macrophage polarization through distinct Ca(2+) entry channel that maintains intracellular signaling and mitochondrial bioenergetics. iScience 2021;24:103339. [Crossref] [PubMed]
74. Chen S, Cui W, Chi Z, et al. Tumor-associated macrophages are shaped by intratumoral high potassium via Kir2.1. Cell Metab 2022;34:1843-1859.e11. [Crossref] [PubMed]
75. Jhunjhunwala S, Hammer C, Delamarre L. Antigen presentation in cancer: insights into tumour immunogenicity and immune evasion. Nat Rev Cancer 2021;21:298-312. [Crossref] [PubMed]
76. Lelis FJN, Jaufmann J, Singh A, et al. Myeloid-derived suppressor cells modulate B-cell responses. Immunol Lett 2017;188:108-15. [Crossref] [PubMed]
77. O'Connell BC, Adamson B, Lydeard JR, et al. A genome-wide camptothecin sensitivity screen identifies a mammalian MMS22L-NFKBIL2 complex required for genomic stability. Mol Cell 2010;40:645-57. [Crossref] [PubMed]
78. Zhang H, Luo X, Li S, et al. Targeting CDK9 With Harmine Induces Homologous Recombination Deficiency and Synergizes With PARP Inhibitors in Ovarian Cancer. Phytother Res 2025;39:3648-63. [Crossref] [PubMed]
79. Srikanth S, Woo JS, Sun Z, et al. Immunological Disorders: Regulation of Ca(2+) Signaling in T Lymphocytes. Adv Exp Med Biol 2017;993:397-424. [Crossref] [PubMed]
80. Chen X, Ma C, Li Y, et al. Trim21-mediated CCT2 ubiquitination suppresses malignant progression and promotes CD4(+)T cell activation in breast cancer. Cell Death Dis 2024;15:542. [Crossref] [PubMed]
81. Chen X, Li J, Zhang R, et al. Suppression of PD-L1 release from small extracellular vesicles promotes systemic anti-tumor immunity by targeting ORAI1 calcium channels. J Extracell Vesicles 2022;11:e12279. [Crossref] [PubMed]
82. Wu Y, Wu Y, Gao Z, et al. Revitalizing T cells: breakthroughs and challenges in overcoming T cell exhaustion. Signal Transduct Target Ther 2026;11:2. [Crossref] [PubMed]
83. Griffin M, Khan R, Basu S, et al. Ion Channels as Therapeutic Targets in High Grade Gliomas. Cancers (Basel) 2020;12:3068. [Crossref] [PubMed]