TIMER2 database (http://timer.cistrome.org/) was used to show the expression difference of NR1H3 between tumor and adjacent normal tissues in different cancer types of the TCGA project (Li et al., 2020). The “Expression analysis-Box Plots” module of the GEPIA2 web server (http://gepia2.cancer-pku.cn/#analysis) was used to obtain box plot of the expression difference between the breast tumor tissues and the corresponding normal tissues of the GTEx database (Tang et al., 2019). Additionally, the NR1H3 expressions in different pathologicals and clinical stages were obtained using the UALCAN database (http://ualcan.path.uab.edu/analysis-prot.html) (Chandrashekar et al., 2017). The Oncomine database (http://www.oncomine.org) was used to validate the expression of NR1H3 in breast cancer (Rhodes et al., 2007).

Paired human breast cancer and adjacent non-tumor paraffin tissue microarrays were purchased from Shanghai Zuocheng Biotech (Shanghai, China). The sections were subjected to antigen retrieval and incubated with primary antibodies against NR1H3 (ab41902, abcam) and CD68 (ab955, abcam) at 4°C in a humid chamber overnight. The next day, the sections were incubated with biotinylated secondary antibody for 60 min. Protein levels of NR1H3 and CD68 were evaluated as follows: the slides were appraised for the intensity of the staining (0–3) and the percentage of positively stained cells (0–4). Index of protein levels was calculated as the intensity of the staining × the percentage of positively stained cells. Therefore, slices were divided into 4 groups: negative (score 0), low expression (score 1–4), medium expression (score 5–8) and high expression group (score 9–12).

The subtypes of breast cancer for sub-group analysis are divided based on the 2013 StGallen criteria using the expression of HER2, ESR1 and MKI67, including basal (ESR1-/HER2-), luminal A (ESR1+/HER2-/MKI67 low), luminal B (ESR1+/HER2+/MKI67 high) and HER2(HER2+/ESR1-).

GEPIA2 and GSCA (http://bioinfo.life.hust.edu.cn/GSCA/#/) databases were used to reveal the correlation between NR1H3 expression and overall survival (OS), disease-free survival (DFS) or progress free survival (PFS) of breast cancer patients (Tang et al., 2019). Kaplan-Meier Plotter (https://kmplot.com/analysis/) was used to assess the effect of NR1H3 on OS, relapse-free survival (RFS), distant metastasis-free survival (DMFS) and post-progression survival (PPS) in breast cancer (Lanczky et al., 2016). Hazard ratios (HRs) with 95% confidence intervals (CI) and log-rank p-values were calculated. Additionally, we constructed univariate and multivariate Cox proportional hazard models. Multivariate analysis comprised seven variables, including the expression of NR1H3 gene, macrophage level, age, tumor stage, gender, race and tumor purity. The survival curves, featuring patterns of NR1H3 gene expression and macrophage level were shown on the diagram. The association between each immune cell type and OS was displayed under the low or high expression of NR1H3.

The correlation between NR1H3 expression and immune infiltration was determined using the TISIDB, TIMER and TIMER2. TISIDB (http://cis.hku.hk/TISIDB/index.php) was used to show the relations between NR1H3 expression and abundance of 28 tumor-infiltrating lymphocytes (TILs) types, immunoinhibitors, immunostimulators, MHC molecules, chemokines and chemokine receptors (Li T. et al., 2017; Ru et al., 2019). The TIMER2 online tool (http://timer.cistrome.org/) was used to analyze the correlation of NR1H3 with the infiltration level and prognostic value of immune cells, including macrophages, CD4⁺ T Cells, CD8⁺ T Cells, monocytes, B Cells, dendritic cells (DCs), neutrophils and natural killer (NK) cells (Li et al., 2020). We also used the TIMER2 to explore the immune infiltration distribution between different somatic copy number changes of NR1H3, and analyze the correlation between the expression of NR1H3 with monocyte markers (CD86, CD115/CSF1R, CD14), macrophage markers (CCL2, CD68, IL10, CD80), M1 macrophage markers (IRF5, INOS/NOS2, COX2/PTGS2), M2 macrophage markers (CD163, VSIG4, MS4A4A) and immune checkpoint molecules (PD-1/CD274, PD-L1/PDCD1, PD-L2/PDCD1LG2 and CTLA-4).

The scRNA-seq database TISCH (http://tisch.comp-genomics.org) was used to show the detailed cell-type annotation at the single-cell level in breast cancer (Sun et al., 2021). Sub-expression analysis of GEPIA 2021 (http://gepia2021.cancer-pku.cn/) visualized the NR1H3 expression in each immune cell type (B Cells, CD4⁺ T Cells, CD8⁺ T Cells, NK cells and macrophages) available in TCGA/GTEx sub-datasets.

The muTarget database (http://www.mutarget.com) is a cancer biomarker/target discovery tool that can identify mutations resulting in expression change. We used the database to predict the mutant genes that affect the expression of NR1H3 gene.

The gene ontology (GO) term enrichment analysis was performed by the LinkedOmics database (http://www.linkedomics.org/) pathway analysis. Gene Set Enrichment Analysis (GSEA) was used to search for Kyoto Encyclopedia of Genes and Genomes (KEGG) and Reactome pathways enrichment analysis. The network neighborhoods of NR1H3 were visualized by ConsensusPathDB-human (http://consensuspathdb.org). These data originate from currently 32 public resources for interactions (Kamburov et al., 2009). The GeneMANIA (http://genemania.org/), an online tool for investigation into associated or similar genes for target genes, was used to validate the gene interaction network results and conduct functional enrichment analysis (Franz et al., 2018).

The ROC Plotter platform (http://www.rocplot.org/) was used to identify NR1H3 whether predicts benefit from endocrine therapy and chemotherapy (Fekete and Győrffy, 2019). The platform integrates multiple gene expression datasets at transcriptome level and contains 3,104 breast cancer patients with treatment and response data. The ROC Plotter is a validation tool for predictive biomarkers.

The Kaplan‐Meier plotter, GSCA and GEPIA2 databases were used for generating survival plots, with data including either HR and p‐values or p‐values derived from a log‐rank test. The Cox proportional hazards regression model was used for univariate and multivariate analyses to evaluate the independence of NR1H3 in predicting prognosis. The correlation of gene expression was assessed by Spearman’s correlation analysis. p-values <0.05 were considered as statistically significant.

To determine the expression pattern of NR1H3 in breast cancer, we analyzed the NR1H3 expression profile based on multiple public databases. As shown in Figures 1A,B, expression of NR1H3 was significantly lower in breast invasive carcinoma (BRCA) compared with adjacent normal tissues. Three breast datasets in the Oncomine were adopted for the validation of lower NR1H3 expression in breast cancer (Supplementary Figure S1A–C). Immunohistochemistry analysis using the tissue microarray (including 83 paired breast cancer and adjacent normal breast tissues) showed that the protein level of NR1H3 was significantly downregulated in breast cancer compared to adjacent normal tissues (Figures 1C,E, p < 0.0001). In addition, according to the clinical data of these cancer cases, we found that the lower expression of NR1H3, the more lymph node metastasis (Figure 1D, p > 0.05). Then we analyzed the expression of NR1H3 in BRCA based on tumor subclasses using the UALCAN database. Luminal (p < 1.0e-12), HER2Pos (p = 4.1e-10), TNBC Basal-like 1 (TNBC-BL1) (p = 5.2e-03), TNBC Basal-like 2 (TNBC-BL2) (p = 1.54e-5), TNBC mesenchymal stem-like (TNBC-MSL) (p = 6.26e-5), TNBC Mesenchymal (TNBC-M) (p = 4.4e-16), TNBC unspecified (TNBC-UNS) (p = 1.0e-4) showed lower NR1H3 expression compared with normal tissues (Supplementary Figure S1D). No statistical difference was found in NR1H3 mRNA expression among different tumor stages (p > 0.05) (Supplementary Figure S1E). Moreover, the expression of NR1H3 was also significantly lower in other types of cancers compared with the corresponding normal tissues, such as colon adenocarcinoma (COAD), kidney chromophobe (KICH), lung adenocarcinoma (LUAD) and thyroid carcinoma (THCA) tissues (Figure 1A).

To evaluate the value of NR1H3 in predicting the prognosis of breast cancer patients, the association between NR1H3 expression and clinical prognosis of OS, DFS and RFS was analyzed in TCGA cohort. Kaplan-Meier survival curves showed that OS was shorter in BRCA patients with low NR1H3 expression in the GEPIA2 and GSCA databases (Supplementary Figure 2A,B). Then we used the Kaplan-Meier plotter approach to conduct a group of survival analyses using gene probe 203920_at. Similarly, the poor prognosis in breast cancer (OS p = 0.007; RFS p = 1.9e-8; DMFS p = 0.004; PPS p = 0.011) was shown to correlate with lower NR1H3 expression (Figure 2A; Supplementary Figure S2C). It is well known that breast cancer is a heterogeneous tumor and is divided into different subtypes based on ER/PR and HER-2 expression (Colombo et al., 2011). As shown in Figure 2; Supplementary Figure 2, OS (p = 0.015), DMFS (p = 0.004) and PPS (p = 0.023) were shorter in basal breast cancer patients with low NR1H3 expression, but not in luminal A, luminal B and HER2+ breast cancer patients (p > 0.05). Moreover, low NR1H3 expression was correlated with poorer prognosis of RFS in basal (p = 5.3e-6), luminal A (p = 3.0e-4), luminal B (p = 0.002) and HER2+ breast cancer patients (p = 0.0003) (Figure 2).

We further explored the prognostic characteristics of NR1H3 under different ER, PR, and HER-2 status. ER-positive subtype had shorter RFS (p = 0.023) in breast cancer with low NR1H3 expression (Figure 2E). Low NR1H3 expression was only correlated with worse PPS in PR-positive subtype (p = 0.03) (Supplementary Figure S2D). ER-negative (RFS p = 8e-04; DMFS p = 0.006) and PR-negative (RFS p = 6.3e-04; DMFS p = 0.045) subtypes were also statistically associated with clinical prognosis of RFS and DMFS, but only a trend towards poor survival without statistical significance of OS and PPS in NR1H3-low breast cancer (Figures 2F,G, Supplementary Figure S2F,G). Compared with high NR1H3 expression, low expression of NR1H3 indicated a worse survival prognosis of OS (p = 0.047), RFS (p = 3e-04) and PPS (p = 0.033) in HER2-positive breast cancer (Figure 2H; Supplementary Figure S2H). Among HER2-negative, only RFS (p = 5.8e-06) and DMFS (p = 0.023) showed statistical survival differences (Figure 2I and Supplementary Figure S2I). In addition, the correlation of NR1H3 expression with clinical and pathological features from Kaplan-Meier Plotter was integrated in Table 1. For instance, RFS was shorter in lymph node positive breast cancer patients with low NR1H3 expression (p = 0.015), but not in lymph node negative breast cancer patients (p = 0.058). For grade 3 breast cancer patients, low expression of NR1H3 indicated a worse survival prognosis of OS (p = 0.01), RFS (p = 0.025) and DMFS (p = 0.035). These results suggest that low NR1H3 expression may be a risk factor for a poor prognosis in breast cancer patients.

GO enrichment analysis revealed adaptive immune response and immune cells activation process were correlated with the expression of NR1H3 in breast cancer (Figures 3A,B). Additionally, the signaling pathways were significantly enriched of NR1H3 by KEGG analysis and Reactome analysis were presented in Figures 3C,D. Tumor-infiltrating immune cells, as prominent components of the TME, are closely linked to the initiation, progression or metastasis of cancer (Fridman et al., 2011; Steven and Seliger, 2018). Here, we found the relations between abundance of 28 TIL types and expression of NR1H3 were strongly correlated across different human cancer types (Supplementary Figure S3). Specifically, NR1H3 expression was closely related with infiltrating levels of TIL in BRCA. Next, we analyzed the correlation between NR1H3 expression and 6 types of infiltrating immune cells (B Cells, CD4⁺ T Cells, CD8⁺ T Cells, neutrophils, macrophages and DCs) in BRCA using TIMER database. Consistently, Figure 4 showed that NR1H3 expression level had significantly positive correlations with infiltrating levels of B Cells (r = 0.178, p = 2.18e-8), CD8⁺ T Cells (r = 0.108, p = 7.62e-4), CD4⁺ T Cells (r = 0.36, p = 7.41e-31), macrophages (r = 0.076, p = 1.68e-2), neutrophils (r = 0.302, p = 1.79e-21), and DCs (r = 0.324, p = 1.22e-24) in BRCA and with negative correlation with tumor purity (r = -0.332, p = 5.56e-27). Moreover, the same trend results were found in each subtype (Figure 4). These findings strongly suggest that NR1H3 is correlated with immune cells infiltration in breast cancer.

Next, TIMER, CIBERSORT, CIBERSORT-ABS, QUANTISEQ, XCELL, MCPCOUNTER and EPIC algorithms were further used to validate the potential relation between the expression of NR1H3 and the infiltration level of 8 types of immune cells (B Cells, CD4⁺ T Cells, CD8⁺ T Cells, monocytes, macrophages, DCs, neutrophils and NK cells) in diverse cancer types of TCGA. As shown in Supplementary Figure S4, CD8⁺ T Cells and macrophages were two immune cell types most strongly correlated with NR1H3 expression in BRCA.

The association between NR1H3 copy number variations (including deep deletion, arm-level deletion, diploid/normal, arm-level gain and amplification) and immune infiltrates in BRCA was investigated using different algorithms of TIMER2. The immune infiltration distribution by the somatic copy number alterations (sCNA) status of NR1H3 across TCGA cancer types was demonstrated in Figure 5A. Then, six of significant relationships between the changes in NR1H3 copy number variations and immune infiltrates in BRCA using TIMER2 were presented (Figure 5B). In particular, arm-level deletion (p = 3e-07), arm-level gain (p = 1.4e-09) and high amplification (p = 0.021) of NR1H3 had significant correlation with CD4⁺ Th2 cell infiltration level using XCELL algorithms. High amplification of NR1H3 was associated with low M2 macrophage infiltration level, compare with the “diploid/normal” status (p = 0.039). By CIBERSORT and CIBERSORT-ABS algorithms, high amplification of NR1H3 had high B Cell and low DC cell infiltration (Figure 5B). However, no statistical difference was found in CD8⁺ T Cell, neutrophil, and NK cell infiltration from TIMER2 (data not shown). These findings indicate the potential mechanism by which NR1H3 alterations affect immune infiltration distribution.

Considering the role of NR1H3 in immune infiltrates in BRCA and its prognostic impact, we used six BRCA data sets (BRCA_SRP114962, BRCA_GSE143423, BRCA_GSE138536, BRCA_GSE136206_mouse_aPD1aCTLA4, GSE114727-inDrop, BRCA_GSE114727_10X and BRCA_GSE110686) in the TISCH platform to analyze the expression of NR1H3 at the single-cell level. The results showed higher NR1H3 expressions in immune cells, mainly in monocyte/macrophage, compared with malignant cells (Supplementary Figure S5). Then we analyzed the GSE114727-inDrop dataset, which is divided into 12 types of cells. Figures 5A–D showed the number of cells in each cell type, with the distribution and number of various TME-related cells presented. In this data set, CD4⁺ T Cells were the most abundant immune cells (n = 5,413), whereas NR1H3 was highly expressed in monocyte/macrophage (Figure 6D). GEPIA2021 platform also revealed the consistent results that NR1H3 is highly expressed in macrophages in BRCA tumor/BRCA normal from TCGA and breast tissue from GTEx using EPIC algorithm (Figure 6E).

We further analyzed the correlation of NR1H3 expression and monocyte/macrophage markers in tumor tissues using TIMER2. We adjusted these results based on tumor purity, and revealed significant correlations between NR1H3 expression and monocyte markers (CD86, CD115/CSF1R, CD14), macrophage markers (CCL2, CD68, IL10, CD80), M2 macrophage markers (CD163, VSIG4, MS4A4A) and M1 macrophage markers (IRF5), whereas the M1 macrophage markers INOS/NOS2 and COX2/PTGS2 showed no correlation with NR1H3 expression (Supplementary Figure S6A).

In order to verify the findings from the database, we detected the protein level of NR1H3 and macrophage marker CD68 in paraffin tissue microarrays from breast cancer patients by IHC. The slides showed that NR1H3 and CD68 protein were expressed in interstitial cells of breast tumor tissues. A typical staining pattern is shown in Figure 7A. However, we did not find a linear relationship between NR1H3 and CD68 expression levels.

As mentioned above, we observed a statistical positive correlation between the immune infiltration of macrophages and NR1H3 expression in BRCA. Then we evaluated the prognostic efficiency of the combination of infiltrated macrophages and NR1H3 expression patterns for breast cancer (Supplementary Figure S6B). The low expression of NR1H3 accompanied by a high level of infiltrated macrophages was associated with poor prognosis in BRCA. However, there was no significant relations between the B Cells/CD4+ T Cells/CD8+ T Cells/neutrophils/NK cells/DCs and prognosis under the low expression level of NR1H3 based on most algorithms (Supplementary Figure S7). Specifically, under low NR1H3 expression, higher macrophage infiltration level had a worse outcome in BRCA using the TIMER algorithm (HR = 1.72, p = 0.0311), compared with lower macrophage infiltration level. On the contrary, the low M2 macrophage infiltration level predicted favorable prognosis under the high expression of NR1H3 using the CIBERSORT algorithm in BRCA (HR = 2.31, p = 0.0135), compared with the high M2 macrophage infiltration level. In BRCA-LumA, under low NR1H3 expression, higher monocyte level had a worse outcome (HR = 3.12, p = 0.0279). The statistically different scatterplot data of the above tumors produced using different algorithms was presented in Figures 7B–E; Supplementary Figure S8.

Additionally, to evaluate whether NR1H3 expression level and macrophage infiltration are independent risk factors for prognosis of BRCA, we conducted the univariate and multivariate analysis included seven variables: macrophage infiltration level, age, stage, gender, race, tumor purity and expression of NR1H3 (Table 2). The results showed that macrophage infiltration (HR = 6.20, p = 0.002), stage 3 (HR = 3.11, p = 0), stage 4 (HR = 13.17, p = 0) and NR1H3 expression (HR = 0.75, p = 0.018) were prognostic variables for the prognosis of OS in BRCA patients. After adjustments of age, stage, gender, race, and tumor purity, the level of macrophage infiltration (HR = 8.44, p = 0.002) and the expression of NR1H3 (HR = 0.73, p = 0.044) were independent prognostic factors in BRCA. These results suggest that NR1H3 is an independent prognostic biomarker and combining its expression level with the macrophage would help to play a more effective role in the prognosis prediction of BRCA.

The TISIDB database was used to infer the correlations between expression of NR1H3 and immunomodulators/chemokines across human cancers. As shown in Supplementary Figure S9, the relations between immunoinhibitors, immunostimulators, MHC molecules, chemokines and chemokine receptors and expression of NR1H3 were strongly correlated. Furthermore, NR1H3 was also positively associated with immune checkpoint molecules (PD-1/CD274, PD-L1/PDCD1, PD-L2/PDCD1LG2, and CTLA-4) in TIMER2 database (Supplementary Figure S10). These results suggest that NR1H3 is closely related to the immune status of human cancers.

To understand the biological function of NR1H3, ConsensusPathDB was used to integrate interaction network of NR1H3 in Homo sapiens. The network defined the neighborhood-based entity set centered by NR1H3 and containing 19 interaction nodes and 22 physical entity nodes (Figure 8A; Supplementary Figure S11A).

A gene interaction network was constructed using the GeneMANIA. Twenty NR1H3-associated genes were observed in the interaction network, functions of which focused on macrophage derived foam cell differentiation, regulation of macrophage derived foam cell differentiation, foam cell differentiation, regulation of interferon-gamma-mediated signaling pathway and regulation of inflammatory response (Figure 8B).

To further investigate the molecular mechanism of the NR1H3 in tumorigenesis, we attempted to screen out the targeting NR1H3-binding proteins and the NR1H3 expression-correlated genes for a series of pathway enrichment analyses. Based on the STRING tool, we obtained a total of 152 NR1H3-binding proteins, which were supported by experimental evidence. We used the GEPIA2 tool to combine all tumor expression data of TCGA and obtained the top 300 genes that correlated with NR1H3 expression. An intersection analysis of the above two groups showed two common members, ITGB2 and ITGB7 (Supplementary Figure S11B–D).

We also analyzed a gene interaction network of NR1H3 and common member ITGB2 using the GeneMANIA. The functions of phagocytosis, toll-like receptor signaling pathway, regulation of innate immune response, macrophage activation, regulation of toll-like receptor 4 signaling pathway, positive regulation of immune effector process and positive regulation of innate immune response were significantly related (Supplementary Figure S12).

In order to reveal the relationship between gene mutation status and NR1H3 expression in breast cancer, we screened mutations resulting in NR1H3 expression change using muTarget tool. The results showed that mutations of TMPRSS15, TBC1D4, ERCC5, ANKRD30A, SPINK5, TNXB, PHF8, FBXW7, ZEB, and RGS22 would lead to the alteration of NR1H3 expression (Figure 9A). Then we verified the above genes in the TIMER2 database and found that the FBXW7 mutation was significantly associated with high NR1H3 expression and high macrophage infiltration (Figure 9B, Figure 9C, Supplementary Figure S13A,B). Higher infiltration of M1 macrophage and lower infiltration of M2 macrophage were shown in the mutant group, compared with the wild type group (Figure 9C; Supplementary Figure S13C).

We used ROC Plotter to identify whether NR1H3 predicted benefit from endocrine therapy and chemotherapy. ROC Plotter showed that NR1H3 was upregulated in responders of luminal A (AUC = 0.564, p = 2.9e-02) and grade 1 subtype breast cancer patients with Tamoxifen treatment (AUC = 0.822, p = 1.1e-05) based on relapse-free survival (RFS) at 5 years (Supplementary Figure S14A). For pathological response, high NR1H3 expression predicted benefit from Anthracycline treatment in TNBC, luminal A, HER2 negative, ER negative, grade 1 and nodal positive subtype patients (Supplementary Figures S14B,C). NR1H3 was upregulated in responders of HER2 negative, ER negative, grade 1, nodal positive subtype patients with Taxane treatment (Supplementary Figure S15A). For patients treated with FAC (Fluorouracil, Adriamycin, Cytoxan), high NR1H3 expression predicted pathological response in luminal A and HER2 negative breast cancer (Supplementary Figure S15B). In addition, NR1H3 was highly expressed in pathological responders receiving FEC (Fluorouracil, Epirubicin, Cyclophosphamide) treatment in grade 1 breast cancer patients and Ixabepilone treatment in HER2 negative patients (Supplementary Figures S15C,D). These data indicated that NR1H3 may function as a predictor of chemoresponsiveness in breast cancer.