The Cancer Genome Atlas (TCGA) website (https://portal.gdc.cancer.gov/) provided the clinical phenotypes and TCGA-BRCA transcriptome data, which included 113 normal samples and 1118 tumor samples. For validation, the GSE199633 dataset, containing microarray data from 637 primary BRAC samples, was obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). This validation dataset was annotated using the GPL15048 platform.

To identify genes with altered expression related to sodium overload, we searched the GeneCards platform (https://www.genecards.org/) (16) for “sodium overload”. This query retrieved a comprehensive set of 2052 genes, including not only those directly involved in sodium transport (e.g., ion channels, exchangers, and pumps) but also a wide array of downstream effectors and genes implicated in the physiological consequences of altered sodium homeostasis. After processing the gene expression data, we performed differential gene expression analysis. Genes exhibiting an absolute log2 fold change (Log₂FC) > 0.585 and a p < 0.05 were defined as differentially expressed. Heatmaps and volcano plots were generated using the ‘pheatmap’ and ‘ggplot2’ R packages, respectively.

To elucidate the biological activities of all identified DESORGs, “clusterProfiler” R package was performed for KEGG and GO enrichment analyses. Additionally, GSEA was conducted on GO and KEGG gene sets to reveal overall functional enrichment patterns across different experimental groups (17). Hallmark pathway scores were computed per sample using the ‘GSVA’ R package (18).

Initial screening via univariate Cox (uni_Cox) regression identified prognostic DESORGs. Subsequently, 101 different modeling approaches were explored by combining ten machine learning algorithms with 10-fold cross-validation using the TCGA-BRCA for our training set. The selected model’s predictive accuracy was independently validated using GSE199633. The model exhibiting the greatest average C-index was selected as an ideal one for further investigation.

SHAP assigns an important value to each gene for every prediction by calculating its average contribution across all possible combinations of genes in the model (19). SHAP values were calculated using the Kernelshap R package to quantify each gene’s impact on the model’s output. The resulting SHAP values were illustrated utilizing shapviz R program to enhance comprehension of each gene’s influence on the predicted outcome.

Time-dependent ROC analysis assessed the DESORG-based risk model’s prognostic performance, with AUC values quantifying predictive accuracy. Kaplan–Meier survival curves were generated, and differences between groups were evaluated using log-rank tests. Multivariate Cox (mul_Cox) regression analysis was performed to determine independent predictors of prognosis. A nomogram was established to display how the risk score and clinical variables jointly predict survival outcomes. Finally, calibration curves were plotted.

MR analysis was performed using five different MR methods with “TwoSampleMR” R package. Single nucleotide polymorphisms (SNPs) used in this analysis were stringently selected based on the following criteria: strong association with the exposure factor (p < 5e-08), absence of linkage disequilibrium (r² threshold below 0.001 across a 10,000 kb genomic region.), and an F-statistic greater than 10. The IVW method served as the primary approach for inferring causality.

Human BC cell lines MCF7 and MDA-MB-231 were maintained at 37 °C in an atmosphere containing 5% CO₂. Cells were cultured in DMEM (Gibco, cat. #11965084) enriched with 10% FBS (Gibco, cat. #10091155) and antibiotic solution (Gibco, cat. #15140163). Small interfering RNAs (siRNAs) targeting human NR1H3, along with a scrambled control siRNA (siRNA-NC), were manufactured by Shenggong Co., Ltd. (Shanghai, China). The siRNA sequences for anti-human NR1H3 were siNR1H3#1: 5’-GCAUCCAGAUAUCUACAAA-3’; siNR1H3#2: 5’-CCACUUCAUGCUGUUGGAA-3’; siNR1H3#3: 5’- GGAAUGCAGCUUCAAGAUG-3’. MCF7 cells were transiently transfected with siRNAs with Lipofectamine RNAiMAX (ThermoFisher Scientific, cat. # 13778030). MDA-MB-231 cells underwent transfection with pENTR221-NR1H3 plasmid (Addgene, cat. # 79514) using Lipofectamine 2000 (ThermoFisher Scientific, cat. # 11668027).

Total RNA was extracted from MCF7 and MDA-MB-231 cells using an RNA miniprep kit (Zymo Research, cat. # R1054) and reverse-transcribed into cDNA with PrimeScript 1st strand cDNA Synthesis Kit (Takara, cat. # 6110A). The reverse transcription reaction was carried out at 37 °C for 15 min, then 85 °C for 5 sec to inactivate the enzyme. Subsequently, qRT-PCR was performed with SYBR Green PCR master mix (ThermoFisher Scientific, cat. # A46012) and qPCR was performed with these primers: NR1H3: F 5’-AATGCTGGGGAACGAGC-3’, R 5’-CGGCATTTGCGAAGCCGAC-3’ and β-ACTIN (control): F 5’-ACCATTGGCAATGAGCGGT-3’, R 5’-GGTCTTTGCGGATGTCCAC-3’. Reactions were set up in triplicate for each biological sample. Amplification was carried out on a real-time PCR instrument under the following cycling conditions: Initial denaturation: 95 °C for 2 min; 40 cycles of: 95 °C for 15 sec, 60 °C for 30 sec (annealing/extension); Followed by a melting (dissociation) curve from 65 °C to 95 °C. Gene expression was quantified via 2^-ΔΔCT analysis.

Total proteins were extracted from cells using RIPA buffer, separated by SDS-PAGE, transferred to PVDF membranes, and blocked with 5% skimmed milk for 1 hour, then subsequently incubated at 4 °C in 5% mike containing primary antibodies NR1H3 (1:5000, proteintech, cat. # 14351-1-AP), and GAPDH (1:10000, proteintech, cat. # 60004-1-Ig) overnight. The membranes were then washed with PBS and incubated with the appropriate peroxidase-conjugated secondary antibodies (1:10000). The images were performed using Immobilon Classico Western HRP substrate (MilliporeSigma, cat. # WBLUC0500) and analyzed with ImageJ software and GraphPad Prism 6. Protein expression levels were standardized to GAPDH expression levels. All Western blot experiments were performed in three independent biological replicates.

Cells (5×10³/well) were plated in 96-well plates and incubated for one week. Cell proliferation was evaluated using CCK8 (MeilunBio, cat. # MA0218) at 37 °C for 2 hours to detect absorbance at 450 nm (OD450) with a microplate reader at different time points.

A total of 1,000 cells were seeded into each well of a 6-well plate and incubated under standard conditions until distinct colonies became visible. The colonies were washed with PBS, fixed in methanol, and stained with 0.1% crystal violet for 20 minutes. Rinse wells gently with water to remove excess stain and allow plates to air-dry. ImageJ software (v1.53) was employed for automated counting. To ensure accuracy, a size threshold of 50 μm in diameter was established to exclude cellular debris from analysis.

Cells (1×10⁵/well) were plated in 12-well plates and cultured for 24 h to form monolayers. A sterile 200 µL pipette tip was used to create uniform scratches, followed by PBS washing to remove debris. Fresh medium was then added, and wound closure was monitored by imaging at 0 h and 24 h post-scratching. The migration rate was quantified by measuring the remaining wound area at both time points.

For migration assays, use uncoated 8 μm pore Transwell inserts (Corning, cat. # 3422). For invasion assays, add 40–50 μL of diluted Matrigel to the upper chamber of each insert, and incubate at 37 °C for 60 minutes to solidify. For both migration and invasion assays, 600 μL of complete medium supplemented 10% FBS was added to the lower chamber, while 200 μL of a cell suspension containing 5 × 10⁴ cells were seeded into the upper chambers of each insert. After 24 hours, remove non-migrated cells from the upper side of the membrane using a cotton swab and fix migrated or invaded cells with 70% ethanol for 10 minutes. Cells were stained with 0.1% crystal violet and then visualized under a phase contrast microscope at 200× magnification in multiple fields to obtain an average.

Potential drugs targeting NR1H3 were identified using the DSigDB database (https://dsigdb.tanlab.org/DSigDBv1.0/). Drug molecule structures were acquired from PubChem (https://pubchem.ncbi.nlm.nih.gov/), while NR1H3’s 3D structure came from the PDB (https://www.rcsb.org/). Subsequently, protein-ligand blind docking was carried out via the CB-Dock2 platform (https://cadd.labshare.cn/cb-dock2/index.php). This process utilized the CurPocket algorithm, which detects surface curvature-based cavities to predict potential binding sites on the NR1H3 protein, followed by performing blind docking of the selected drug molecules to these identified regions.

Data were performed with GraphPad Prism 6 (GraphPad Software) and R 4.2.2 (R Foundation). Data are presented as mean ± SD. For in vitro experiments, two-group comparisons used unpaired t-tests; multi-group comparisons employed ANOVA. Three independent biological replicates were performed, and p < 0.05 was considered statistically significant.

2052 genes related to the sodium overload pathway were retrieved from the GeneCards database. Differential expression analysis revealed a total of 753 genes exhibited significantly differential expression in tumor tissues relative to norma samples. Among the differentially expressed SORGs, 370 genes showed down-regulated and 383 exhibited up-regulated in tumor samples. A heatmap of DESORGs was shown in Figure 1A . The volcano plot in Figure 1B highlights both top 10 of down-regulated and up-regulated DESORGs with the most significant false discovery rate (FDR) values. Among the downregulated genes, the most significantly changed were VEGFD, CAVIN2, SCARA5, CA4, CAV1, MME, DMD, ADRB2, SLC2A4, and NPR1. Conversely, the top upregulated genes included CDKN3, INHBA, AURKB, MMP13, EZH2, KIF23, LMNB1, NME1, CCL11, and GFUS. To elucidate the functional implications of all DESORGs, KEGG was conducted and revealed significant enrichment of DESORGs in pathways including PI3K-Akt signaling, calcium signaling, lipid and atherosclerosis, as well as MAPK signaling ( Figure 1C ). Further, GO analysis ( Figure 1D ) indicated that the DESORGs participate in biological processes such as the positive regulation of phosphorylation, response to steroid hormone, and response to oxidative stress. In terms of cellular components, DESORGs were involved in the sarcomere, myofibril, and contractile muscle fiber. Regarding molecular function, the DESORGs showed enrichment in cytokine activity, growth factor activity, and protease binding.

Univariate Cox regression analysis ( Figure 2A ) revealed 67 DESORGs significantly linked to breast cancer prognosis, including 36 potentially protective genes (hazard ratio [HR] < 1), while 31 genes were associated with increased risk (HR > 1). To develop a robust prognostic model, we evaluated 101 models generated from combinations of 10 different machine learning algorithms ( Figure 2B ). The Ridge regression model outperformed other models, achieving greatest average C-index of 0.692 across all evaluated models. Specifically, it yielded a C-index of 0.775 in the TCGA training cohort and 0.609 in the GEO validation cohort, indicating consistent prognostic ability across different datasets. The Ridge model stratified TCGA patients into high/low-risk groups, with poorer survival in high-risk cases (p < 0.001; Figure 2C ). This finding was consistently observed in the GEO validation cohort, where high-risk patients also had significantly worse OS (p = 0.001; Figure 2D ). Furthermore, high-risk patients in testing cohort also exhibited shorter progression-free survival (PFS) (p < 0.001; Figure 2E ). In conclusion, these results underscore the predictive capability of the Ridge regression model for both OS and PFS in breast cancer patients.

To illustrate the interpretability of the Ridge model, SHAP values were used to explain feature importance and model predictions. The bar chart in Figure 3A illustrates the top-ranked features based on their mean absolute SHAP values, reflecting the average impact of each feature on the model’s predictions. Among these, IFNG is shown as the most important feature with a mean |SHAP value| of 0.329, followed by TFF1 (0.271), TRPM2 (0.259), RPA3 (0.254), SRD5A2 (0.252), and others. SHAP summary plot shown in Figure 3B provided a more detailed view of feature effects. Each row corresponds to a feature, ordered by global importance. Features where high values are predominantly associated with positive SHAP values include IFNG, TFF1, TRPM2, RPA3, SRD5A2, PGK1, TAGLN2, ADAMTSL1, EGR2, RACGAP1, NRG1, ABCD2, NFE2. Conversely, features like SOCS3 and ALDH3A1 show an opposite trend: high values are predominantly associated with negative SHAP values, indicating a tendency to decrease the model’s predicted outcome. SHAP waterfall plot ( Figure 3C ) shows how the model arrived at a prediction f(x) = 1.9 starting from a baseline expected value E[f(x)] = 2.77 (the average prediction over the dataset). Features contributing positively to the model’s prediction are highlighted in yellow with their corresponding positive SHAP values (e.g., KRT14 = 11.5 contributes +0.411, RACGAP1 = 3.72 contributes +0.457, EGR2 = 6.43 contributes +0.481). Features contributing negatively are shown in purple/maroon with their negative SHAP values (e.g., “20 other features” collectively contribute -1.36, SOCS3 = 7.77 contributes -0.458, ALDH3A1 = 1.79 contributes -0.405).

We aimed to create a practical clinical tool for overall survival (OS) prediction, and integrated the risk score derived from the Ridge regression model with relevant clinical variables from the TCGA dataset, including age, pathological stage, tumor size (T classification), lymph node (N classification), and metastasis (M classification). Univariate Cox proportional hazards regression analysis ( Figure 4A ) revealed that all considered variables were associated with increased risk, with the “riskScore” exhibiting the strongest association (HR = 3.491, p < 0.001). In the multivariate Cox regression analysis ( Figure 4B ), the “riskScore” remained a significant independent predictor of OS (HR = 3.168, p < 0.001), even after adjusting for other clinical factors. Age (HR = 1.030, p < 0.001) and pathological stage (HR = 1.609, p = 0.040) also were independent predictors. ROC analysis compared the ‘Risk’ model’s predictive accuracy against individual clinical variables, with AUC values calculated for quantitative evaluation ( Figure 4C ). The “Risk” model achieved a superior AUC of 0.845 compared to Age (0.611), Stage (0.722), T stage (0.631), M stage (0.578), and N stage (0.650), indicating its enhanced ability to discriminate between patients with different survival outcomes. Time-dependent ROC analyses further demonstrated the model’s consistent predictive performance over time, with AUC values of 0.845, 0.807, and 0.779 at 1, 3, and 5 years, respectively ( Figure 4D ). Additionally, the risk score consistently exhibited a higher C-index compared to traditional clinical variables across different time points, further supporting its superior prognostic value ( Figure 4E ). Furthermore, a nomogram was developed that integrates clinical variables with the risk score to predict OS probabilities ( Figure 4F ). To assess the nomogram’s accuracy, we constructed a calibration curve, which demonstrated satisfactory performance. The nomogram also exhibited good discriminative ability, with a C-index of 0.815 (95% CI: 0.777-0.853; Figure 4G ).

To assess the robustness of our risk stratification model, its prognostic power within various clinically defined subgroups were further evaluated. Figures 5A-E demonstrates the risk score’s robust prognostic value, with low-risk patients showing superior survival versus high-risk counterparts across all clinical strata (age, pathological stage, TNM classification; all p<0.001). Similarly, the risk score demonstrated consistent prognostic ability within molecularly defined subgroups, including molecular subtypes (Luminal, HER2, and TNBC), ER status, PR status, and HER2 status (all p < 0.001) ( Figures 5F-I ). These findings underscore the robust and independent prognostic value of our risk score across a wide spectrum of established clinical and molecular prognostic factors, consistently identifying patient groups with differing survival probabilities regardless of these other variables.

MR analysis was performed to assess causal links between the 67 prognostic DESORGs from the Ridge regression model and breast cancer risk. Figures 6A-C showed the estimated effect of each individual genetic variant (SNP, listed on the y-axis) used as an instrumental variable on breast cancer risk. In Figures 6D-F , these scatter plots where each point represents an instrumental SNP. The x-axis shows the SNP’s effect on breast cancer (β_XY), and the y-axis shows the inverse of the standard error of this effect (1/SE(β_XY)), indicating precision. Leave-One-Out sensitivity analysis ( Figures 6G-I ) shown how the overall MR estimate for the effect of ADAM15, HLA-F and NR1H3 on breast cancer changes by sequentially excluding each single nucleotide polymorphism (SNP). Figures 6J-L displayed the relationship between each SNP’s effect on the exposure (e.g., “SNP effect on ADAM15”) and its corresponding effect on breast cancer risk. Various MR methods were employed ( Figure 7A ). Weighted median, IVW and Simple mode methods show ADAM15 has a statistically significant odds ratio. All MR methods indicate HLA-F has a strong and statistically significant causal risk effect, with ORs substantially greater than 1 (e.g., IVW OR = 3.115, 95% CI: 2.859-3.393; Weighted mode OR = 2.781, 95% CI: 1.366-5.659). Most MR methods suggest NR1H3 has a strong and statistically significant protective causal effect, with ORs substantially less than 1 (e.g., IVW OR = 0.132, 95% CI: 0.060-0.289; Weighted mode OR = 0.105, 95% CI: 0.080-0.136). A circos plot designed to visualize genomic information, including the chromosomal locations of genes ( Figure 7B ): ADAM15 is located on chromosome 1, HLA-F is on chromosome 6, and NR1H3 is on chromosome 11. Compared to normal tissues, NR1H3 expression is significantly downregulated in tumor samples, whereas HLA-F and ADAM15 expression are significantly higher in tumor tissues ( Figure 7C ). ROC curves were generated to evaluate the ability of each gene to discriminate tumor and normal tissues. NR1H3 has the highest AUC of 0.817, ADAM15 has an AUC of 0.768, HLA-F has an AUC of 0.598 ( Figure 7D ), suggesting NR1H3 expression has the best discriminatory power among the three for distinguishing tumor from normal tissue. Kaplan-Meier survival analyses showed that only NR1H3 is associated with OS (p = 0.02) ( Figure 7E ). Higher expression of NR1H3 predicted a significantly better OS. Although lower HLA-F expression showed a trend toward better OS, this did not reach statistical significance (p = 0.068) ( Figure 7F ). No significant association was observed between ADAM15 expression and OS (p = 0.092) ( Figure 7G ).

To further explore the function of NR1H3 expression in breast cancer cell lines, we overexpressed NR1H3 in MDA-MB-231 cells and knockdown NR1H3 in MCF7 cells. The mRNA and protein levels were verified by qRT-PCR and western blot, respectively ( Figures 8A, B ). siNR1H3#2 exhibiting the best knockdown efficiency was selected for further experiments. Knockdown of NR1H3 (siNR1H3#2) in MCF7 cells significantly increased cell proliferation over 7 days compared to the control, while overexpressed NR1H3 in MDA-MB-231 had opposite effect ( Figure 8C ). Similarly, reducing NR1H3 expression increased the number of colonies formed in MCF7 cells, overexpression of NR1H3 in MDA-MB-231 cells resulted in a marked reduction in colony formation ( Figure 8D ). The wound healing assay demonstrated that NR1H3 downregulation increases while NR1H3 upregulation decreases the rate of wound closure (cell migration) at 24 hours compared to the control group ( Figure 8E ). Consistently, NR1H3 siRNA#2 significantly increased both cell migration and invasion through transwell membrane. In contrast, NR1H3 overexpression led to a substantial decrease in these invasive behaviors ( Figure 8F ). Collectively, these findings suggest that NR1H3 functions as a tumor suppressor in breast cancer by inhibiting cell proliferation, colony formation, migration and invasion.

To explore the pathways that NR1H3 involved in breast cancer, GSEA analysis was conducted. Among pathways enriched in samples with high NR1H3 expression, the top five were primarily linked to immune responses: GOBP_ANTIGEN_PROCESSING_AND_PRESENTATION_OF _EXOGENOUS_ ANTIGEN, GOCC_ MHC_PROTEIN_COMPLEX, GOMF_ANTIGEN_BINDING and GOMF_PEPTIDE_ANTIGEN_BINDING ( Figure 9A ). The enriched pathways in low NR1H3 group were more diverse and include: GOBP_AEROBIC_RESPIRATION, GOCC_PRESYNAPTIC_ACTIVE_ZONE_CYTOPLASMIC_ COMPONENT, GOBP_SENSORY_PERCEPTION_OF_TASTE and GOMF_TASTE_RECEPTOR_ ACTIVITY, which indicates that low NR1H3 levels are associated with alterations in cellular respiration and some neuronal or sensory-related pathways ( Figure 9B ). GSVA for KEGG pathways showed their correlation with NR1H3 expression (indicated by t-values). Pathways positively correlated with NR1H3 are again heavily involved in immune processes, such as: KEGG_CYTOSOLIC_DNA_SENSING_PATHWAY, KEGG_SYSTEMIC_LUPUS_ERYTHEMATOSUS, and KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION, while pathways negatively correlated with NR1H3 include processes like: KEGG_PROTEIN_ EXPORT and KEGG_UBIQUITIN_ MEDIATED_PROTEOLYSIS ( Figure 9C ). Similarly, GSVA for GO pathways show that pathways positively correlated with NR1H3 further confirm the strong association with immune activation, while pathways negatively correlated with NR1H3 are related to chromosomal organization and morphogenesis ( Figure 9D ).

DSigDB_All_detailed data was downloaded from DSigDB database. 155 drugs were identified associated with NR1H3. Among them, Cephaeline and Emetine were identified as potential upregulators of NR1H3 via the Connectivity Map (CMAP) database. This database contains gene expression data from cell lines treated with various compounds. CMAP data suggested that both Cephaeline and Emetine treatment lead to an upregulation of NR1H3 expression in MCF7 and HL60 cells. Following this lead from the gene expression data, we then performed molecular docking to investigate a plausible mechanism. Molecular docking shows that there are five potential binding pockets (C1 through C5) for the drug Cephaeline on the NR1H3 protein. The predicted binding affinities (Vina scores) for Cephaeline range from -9.2 kcal/mol (strongest binding, pocket C1) to -7.1 kcal/mol (weakest among those listed, pocket C5) ( Table 1 ). The binding pocket C1 was shown in Figure 10A displays visualizations of protein-ligand interactions between NR1H3 and Cephaeline. The ligand is shown interacting with 39 specific amino acid residues of NR1H3 protein in chain A, chain B and chain C. There are also five potential binding pockets (C1 through C5) for the drug Emetine on the NR1H3 protein. The predicted binding affinities (Vina scores) for Emetine range from -9.1 kcal/mol (strongest, pocket C1) to -7.3 kcal/mol (pocket C3) ( Table 2 ). The ligand is shown interacting with 31 specific amino acid residues of NR1H3 protein in chain A, chain B and chain C ( Figure 10B ). The favorable binding affinities (Vina scores of -9.2 kcal/mol for Cephaeline and -9.1 kcal/mol for Emetine) suggest a strong and stable interaction is possible.