Identification and validation of a novel estrogen-related model for breast cancer to predict the prognosis

Introduction

Breast cancer (BRCA), as one of the most common malignant tumors in women globally, has continuously posed a significant challenge to public health due to its high incidence and mortality rates (1). Despite the advancements in BRCA prevention, diagnosis and treatment in recent years, it remains a major threat to women’s health worldwide (2). Traditional treatment approaches for BRCA primarily include surgical resection, radiation therapy, chemotherapy, and endocrine therapy (3). These treatment methods have greatly improved patient survival rates and extended their lives, but the survival rate of BRCA patients remains dissatisfactory (4). With the advancement of molecular biology, molecular targeted therapy has emerged as a promising treatment strategy in cancers, owing to its specificity and minimal side effects (3).

Estrogen, as a crucial hormone in the female body, not only regulates the reproductive system and bone health but is also closely associated with the occurrence and development of BRCA (5,6). By binding to estrogen receptors (ERs), estrogen activates several downstream signaling pathways, including PI3K-Akt, MAPK, and TGF-β, thereby regulating cell proliferation, differentiation, and apoptosis (7,8). Aberrant activation of the estrogen signaling pathway can lead to uncontrolled cell proliferation, metastasis, and increased therapeutic resistance, particularly in ER-positive BRCA (9-11). At present, various prognostic models for BRCA have been established, most of which are based on clinical features or a limited number of molecular biomarkers. However, these models generally exhibit limited predictive accuracy and insufficient applicability across different populations, which restricts their clinical utility (12-14). Recent studies have shown that estrogen-related genes (ERGs) possess potential prognostic value in several cancers, including prostate cancer and corpus endometrial carcinoma (15,16).

Therefore, in-depth investigation of the role of ERGs in BRCA prognosis is of significant importance in guiding clinical treatment and improving patient outcomes. This article aims to explore the predictive value of ERGs in BRCA prognosis and how these genes influence treatment strategies and patient survival. By delving into the study of these genes, more personalized and precise treatment approaches can be developed, further enhancing the survival rates and quality of life for BRCA patients. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1409/rc).

Methods

Data acquisition and identifying estrogen-related differentially expressed genes (ERDEGs) in BRCA

Expression profile and clinical data of BRCA were downloaded from The Cancer Genome Atlas (TCGA, https://xenabrowser.net/hub/). An external validation dataset GSE42568 was acquired from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/). Information on 299 ERGs was extracted from a previous literature (17). The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

The differential gene expression analysis between normal and BRCA samples in TCGA data was conducted using the ‘limma’ package in R. Raw count data were normalized with TMM and log2-transformed. PCA was used for quality control to exclude outlier samples. A linear model in ‘limma’ compared two groups, with empirical Bayes moderation enhancing statistical reliability. Significantly differentially expressed genes (DEGs) were identified with P<0.05 and |log fold change (FC)| ≥1. Then, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed on the DEGs using the R package clusterProfiler. The GO terms contain biological process (BP), cellular component (CC), and molecular function (MF). Finally, DEGs were intersected with ERGs to identify ERDEGs using a Venn diagram.

The construction and validation of prognostic model

Based on the ERDEGs, univariate Cox regression analysis was conducted to identify genes significantly associated with the prognosis of BRCA. Subsequently, the least absolute shrinkage and selection operator (LASSO) regression analysis was used to prevent over-fitting using the glmnet package. The candidate genes selected by LASSO were then subjected to multivariate Cox regression analysis to determine the hub genes and their corresponding coefficients. A prognostic risk model was constructed based on the expression and coefficient value of hub genes, and individual risk scores were calculated accordingly. Patients were divided into low- and high-risk groups based on the optimal cutoff value of risk score determined by the maxstat package. Finally, the Kaplan-Meier (K-M) and time-dependent receiver operating characteristic (ROC) curves at 1-, 3- and 5-year intervals were performed to evaluate the predictive performance and stability of the risk model using the survival and qROC packages.

Immune infiltration analysis

The stromal score, immune score, ESTIMATE score, and tumor purity were acquired using the ESTIMATE algorithm of the IBOR R package. The infiltration levels of immune cells in the high- and low-risk groups were detected using the CIBERSORT algorithm of the IBOR package. Additionally, the Tumor Immune Dysfunction and Exclusion (TIDE) scores were predicted on the TIDE website (tide.dfci.harvard.edu/) and compared between the high- and low-risk groups.

Mutation analysis

After the mutation data of BRCA samples were obtained from the TCGA dataset, we analyzed the mutation profile through the R package “maftools”.

Survival analysis of hub genes

The K-M curve with hazard ratio (HR) and 95% confidence interval (CI) was used to evaluate the prognosis value between high- and low-expression of hub genes. Firstly, we calculated the optimal truncation value of the hub genes using the maxstat of the R package, and the cases of BRCA were divided into two groups according to the optimal truncation value. Then we assessed the difference between the two groups in the prognosis using the survfit function of the R package.

Gene set enrichment analysis (GSEA)

GSEA was performed by utilizing the R package and clusterProfiler to show the differential pathways between the high- and low-expression groups of hub genes. Gene set permutations were performed 1,000 times for each analysis. The pathways significantly enriched were adopted as adjusted P value of <0.05, False discovery rate q-value of <0.25 and normalized enrichment score >1.0. All reference gene sets were obtained in the Molecular Signatures Database (MSigDB, https://www.gseamsigdb.org/gsea/downloads.jsp).

Drug sensitivity analysis

Gene set cancer analysis (GSCA) is an integrated platform for the analysis of gene expression, immune response, mutation and drug sensitivity. We assessed the relationship between gene expressions and drug sensitivity in GSCA according to the Cancer Therapeutics Response Portal (CTRP) database.

Cell culture

The breast epithelial cell lines MCF 10A and human BRCA cell line MDA-MB-231 used in this study were purchased from Wuhan Pricella Biotechnology Co., Ltd (Wuhan, China) and routinely maintained in the DMEM/high glucose medium (Gibco-BRL, Rockville, IN, USA) with 10% fetal bovine serum (FBS, Haoyang Biological Manufacture, Tianjin, China), and 1% penicillin-streptomycin at a 37 ℃ cell culture incubator with 5% CO₂.

Cell transfection and real-time polymerase chain reaction (RT-PCR)

Transfection experiments were performed using Lipofectamine 2000 (Invitrogen) according to the manufacturer’s instructions. Briefly, cells were seeded at a density of 5×10⁵ per well in 6-well plates one day before transfection. Once the cells reached approximately 80% confluence, the plasmid DNA and Lipofectamine 2000 were separately diluted in Opti-MEM I Reduced Serum Medium. The diluted DNA and transfection reagent were then combined and incubated at room temperature for 20 min to allow complex formation. The resulting mixture was gently added dropwise to the cells. After incubation for 24–48 h, the transfected cells were harvested and prepared for subsequent experiments. Total RNA was extracted using TRIzol Reagent (Invitrogen) as the manufacturer’s instructions. First-strand cDNA synthesis was performed using the Revert Aid First Strand cDNA Synthesis Kit (Thermo Scientific). RT-PCR was conducted using SYBR Green qPCR Master Mix (ACCURATE BIOTECHNOLOGY, Changsha, China) on a Roche 480II system as follows: initially, 95 ℃ for 2 min for denaturation, followed by 40 cycles at 95 ℃ for 15 s and 60 ℃ for 1 min, then 95 ℃ for melting, and finally 50 ℃ for 30 s for cooling. The primer sequences were as follows:

BATF: F: CCCAGAAGAGCCGACAGA; R: TGCTTGATCTCCTTGCGTAG;
CLDN7: F: GGGGAGACGACAAAGTGAA; R: TGCCCAGCCAATAAAGAT;
TFPI2: F: GCCAATGTGACTCGCTAT; R: ATCCTCCCTGCTAACAAA;
TH: F: CACAGTCTTTGCTCCCACA; R: CCGCTAACCAGGATTTCAT;
GAPDH: F: AGAAGGCTGGGGCTCATTTG; R: AGGGGCCATCCACAGTCTTC.

Colony-formation assay

CLDN7 overexpression cells were digested by trypsin and seeded in a six-well plate at a density of 1,000 cells/well. MDA-MB-231 cells were cultured for 15 days. Then, the clones were washed with PBS, fixed with methanol for 5 min, and stained with 0.1% crystal violet.

Transwell assay

Cell migration ability was evaluated by transwell assay using a Matrigel-precoated Transwell chamber with a pore size of 8.0 µm according to the instructions of the manufacturer. MDA-MB-231 cells (1×10⁵) were suspended in serum-free medium and plated in the upper wells. The medium containing 20% FBS was added to the lower chamber as a chemoattractant. MDA-MB-231 cells were cultured for 24 h. After being fixed with methanol for 5 min, the chambers were stained with 1% crystal violet solution for 5 min. Then, the cells in the lower chamber were observed under an inverted microscope.

Scratch assay

After being transfected with selected genes, MDA-MB-231 cells were seeded on a 24-well plate at a density of 1×10⁵/well. A straight-line cell-free “scratch’’ was created by pipette tips, and a horizontal line at the back of the plate was drawn as a reference point to guarantee the same area of image acquisition. After washing with PBS to remove the debris, the plate was incubated in 5% CO₂ at 37 ℃. The migration speed was measured by calculating the difference in the distances between the two edges of the scratch at 0 and 24 h.

Western blotting

The proteins of the distinctively treated cells were lysed in a lysis buffer in the presence of protease inhibitors. After centrifugation at 12,000 rpm for 20 min at 4 ℃, the supernatant was collected. Next, 30 µg of protein was separated by 10% SDS-PAGE and transferred (100 V, 2 h) onto polyvinylidene fluoride membranes (Millipore, Bedford, MA, USA). After blocking with 5% nonfat milk for 1 h, the membranes were incubated overnight at 4 ℃ with the primary antibodies. The following primary antibodies were used: GAPDH poloclonal antibody (10494-1-AP, Proteintech, Wuhan, China), CLDN7 mouse monoclonal antibody (M1412-2, Huabio, Hangzhou, China), TH (J-AB499173, SanjingBio, Hangzhou, China), TFPI2 (SJ-AB477748, SanjingBio, Hangzhou, China), BATF (13507-1-AP, SanjingBio, Hangzhou, China), Phospho-STAT3 (p-STAT3) recombinant rabbit monoclonal antibody (ET1603-40, Huabio, Hangzhou, China), Phospho-Smad2/3 (p-Smad2/3) recombinant rabbit monoclonal antibody(ET1609-41, Huabio, Hangzhou, China), E-Cadherin recombinant rabbit monoclonal antibody (ET1607-75, Huabio, Hangzhou, China) and Vimentin recombinant rabbit monoclonal antibody (ET1610-39, Huabio, Hangzhou, China). After washing with TBS-T, the membrane was labeled with the secondary antibody for 1 h. The following secondary antibodies were used: HRP Conjugated Goat anti-Rabbit IgG polyclonal Antibody (HA1001, Huabio, Hangzhou, China) and HRP Conjugated Goat anti-Mouse IgG polyclonal Antibody (HA1006, Huabio, Hangzhou, China). Protein spots were visualized by ECL. GADPH was used as the endogenous control.

Statistical analysis

In this study, all data were analyzed in SPSS 25, R Studio and GraphPad Prism software. The independent sample t-test was used to calculate the difference between the two groups. The survival difference was evaluated by the log-rank test, and the correlation analysis was identified by the Pearson test. For cell experiments, all procedures were conducted in triplicate. The data were shown as mean ± standard deviation, and a comparison between the two groups was performed using t-test. P<0.05 was considered a significant difference.

Results

Identification of ERDEGs in BRCA

After the mRNA expression profile of tumor and control samples were obtained from TCGA, the differentially expressed analysis was performed and 7,382 DEGs were identified (Figure 1A). Then, KEGG enrichment analysis exhibited that DEGs were mainly involved in the PI3K-Akt signaling pathway (hsa04151), MAPK signaling pathway (hsa04010), focal adhesion (hsa04510), cellular senescence (hsa04218) and autophagy (hsa04140) (Figure 1B). Besides, the GO analysis results illustrated that DEGs mainly existed in the cell-substrate junction (GO:0030055), focal adhesion (GO:0005925), mitochondrial matrix (GO:0005759), chromosomal region (GO:0098687) and spindle (GO:0005819) (Figure 1C). Furthermore, these DEGs were enriched in small GTPase mediated signal transduction (GO:0007264), proteasome-mediated ubiquitin-dependent protein catabolic process (GO:0043161), cellular component disassembly (GO:0022411), mitotic cell cycle phase transition (GO:0044772) and positive regulation of cell adhesion (GO:0045785) in terms of BP (Figure 1C). In the aspect of MF, DEGs were mainly involved in cadherin binding (GO:0045296), DNA-binding transcription factor binding (GO:0140297), actin binding (GO:0003779), GTPase binding (GO:0051020) and ubiquitin-like protein ligase binding (GO:0044389) (Figure 1C). After crossing 7,382 DEGs with 299 ERGs, 113 ERDEGs were identified (Figure 1D).

Figure 1 Identification of ERDEGs in BRCA. (A) The volcano plot exhibited the DEGs between BRCA and normal samples. (B) The KEGG enrichment analysis of DEGs. (C) The GO enrichment analysis of DEGs. (D) The Venn diagram showed the overlap between DEGs and ERGs to obtain ERDEGs. BP, biological process; BRCA, breast cancer; CC, cellular component; DEGs, differentially expressed genes; ERDEGs, estrogen-related differentially expressed genes; ERGs, estrogen-related genes; MF, molecular function; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Estrogen-related prognostic risk model for BRCA

Then, the Cox analysis was used to determine the correlation between ERDEGs and prognosis. The result of univariate Cox analysis showed that seven genes, namely TH, BATF, SIAH2, TFF1, CLDN7, TFPI2 and ASCL1, were significantly related to the prognosis of BRCA (Figure 2A). Then, the result of LASSO Cox analysis identified that there were no collinearity and confounding factors among these seven genes (Figure 2B,2C). Thus, all of these genes were input to the multivariate Cox analysis, and four hub genes (BATF, CLDN7, TFPI2 and TH) were identified to independently predict the prognosis of BRCA (Figure 2D). Then, a four-gene risk model based on these hub genes was constructed, and according to the coefficient and expression of hub genes, the risk score was calculated: risk score = −0.015*BATF + 0.011*CLDN7-0.017*TFPI2 + 0.206*TH. The K-M curve showed that the patients with high risk scores had a lower survival probability than those with low risk scores (P<0.001, Figure 2E). Moreover, the ROC curve showed that AUC values of the risk score for predicting 1-, 3- and 5-year survival probability were 0.74, 0.71 and 0.69, respectively (Figure 2F).

Figure 2 Construction and validation of the prognostic risk model. (A) Univariate Cox analysis was performed on 113 ERDEGs. (B,C) The LASSO Cox analysis for seven genes identified by univariate Cox analysis. (D) Four hub genes were identified from the above seven genes through multivariate Cox analysis. (E) The Kaplan-Meier curves of high-risk and low-risk groups. (F) The ROC curves revealed the prognostic value of the risk model at 1-, 3- and 5-year in the TCGA dataset. AUC, area under the curve; CI, confidence interval; ERDEGs, estrogen-related differentially expressed genes; LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic.

The immune microenvironment and immunotherapy between high- and low-risk groups

As shown in Figure 3A,3B, ESTIMATE results indicated that the stromal score, immune score and ESTIMATE score were all significantly decreased in high-risk group (P<0.001, Figure 3A), while tumor purity was significantly increased in high-risk group (P<0.001, Figure 3B). Besides, the CIBERSORT result indicated that the levels of CD8 T cells, Tregs, activated natural killer (NK) cells and M1 macrophages were significantly decreased, but M0 macrophage and M2 macrophage were significantly increased in high-risk group (P<0.05, Figure 3C). Moreover, the TIDE score was lower in high-risk group compared with low-risk group (Figure 3D), whereas microsatellite instability (MSI) was obviously higher in high-risk group (Figure 3E).

Figure 3 Immune infiltration and TIDE analysis in different risk groups. (A,B) The discrepancies of immune score, stromal score and estimate score between high-risk and low-risk groups. (C) The infiltration level of immune cells between high-risk and low-risk groups. (D) The TIDE score between high-risk and low-risk groups. (E) The MSI score between high-risk and low-risk groups. The blue indicated low-risk group, while red indicated high-risk group. MSI, microsatellite instability; NK, natural killer; TIDE, Tumor Immune Dysfunction and Exclusion.

Mutation landscape between high- and low-risk groups

According to the MSI score, the gene mutation may be important for the development of BRCA. Thus, the mutation analysis was performed. As shown in Figure 4A,4B, it was obvious that the mutation of TP53 was remarkably downregulated in the high-risk group.

Figure 4 Mutation landscape between high-risk and low-risk groups. (A) The mutation landscape in high-risk group. (B) The mutation landscape in low-risk group. TMB, tumor mutational burden.

Expression and prognostic value of hub genes

From Figure 5A, BATF, CLDN7, and TH expressions were upregulated in tumor samples compared with normal samples in the TCGA dataset, while TFPI2 was downregulated in tumor samples (P<0.001). Consistent with the bioinformatics results, the qRT-PCR results illustrated that BATF, CLDN7 and TH were highly expressed in BRCA cells MDA-MB-231 compared to normal cells MCF 10A, while TFIP2 was down-regulated in MDA-MB-231 cells (P<0.001, Figure 5B). The expression of hub genes was then validated in an external dataset GSE42568. Similarly, BATF and CLDN7 expressions were both upregulated in tumor samples, and TFPI2 was downregulated in tumor samples (P<0.001, Figure 5C). However, there was no significant difference in TH expression between normal and tumor samples (Figure 5C). Besides, the K-M curve showed that high expression of CLDN7 and TH was closely related to the poor prognosis of BRCA, but the low expression of BATF and TFPI2 was closely related to the poor prognosis of BRCA (Figure 5D-5G).

Figure 5 Expressions and survival curves of hub genes. (A) The expression of hub genes in the TCGA cohort. (B) The expression of hub genes in MCF 10A and MDA-MB-231 cell lines. (C) The expression of hub genes in the GEO cohort. (D) The survival curve of CLDN7 in the TCGA cohort. (E) The survival curve of BATF in the TCGA cohort. (F) The survival curve of TH in the TCGA cohort. (G) The survival curve of THPI2 in the TCGA cohort. GEO, Gene Expression Omnibus; OS, overall survival; TCGA, The Cancer Genome Atlas.

Single-gene GSEA functional enrichment analysis

To further explore the potential mechanism of hub genes, the single-gene GSEA was performed based on the hallmark dataset. The results showed that high expression of BATF was enriched in IL2-STAT5 and P53 signal pathways and apoptosis (Figure 6A). High-expressed TH was mainly involved in glycolysis, oxidative phosphorylation and PI3K-AKT-mTOR signal pathway (Figure 6B). Besides, down-regulated CLDN7 was closely related to IL2-STAT5, IL6-JAK-STAT3 and TGF-β signal pathways (Figure 6C). In addition, high expression of TFPI2 was enriched in TNFα signal pathway, but the low expression of TFPI2 was enriched in Notch signal pathway and epithelial mesenchymal transition (Figure 6D).

Figure 6 Single-gene GSEA results of hub genes. (A) The single-gene GSEA results of BATF. (B) The single-gene GSEA results of TH. (C) The single-gene GSEA results of CLDN7. (D) The single-gene GSEA results of TFPI2. GSEA, gene set enrichment analysis.

Drug sensitivity analysis

To assess the clinical value of hub genes in chemotherapy, the drug sensitivity analysis was performed. As shown in Figure 7, it was obvious that TFPI2 and CLDN7 were positively and closely related to most drugs in CTRP database, while TH and BATF were negatively and closely related to most drugs.

Figure 7 Drug sensitivity analysis of four hub genes based on CTRP database. CTRP, Cancer Therapeutics Response Portal; FDR, false discovery rate.

Functional validation of four hub genes in RBCA cells

In the above analyses, CLDN7 was upregulated in the BRCA samples and cells, and its upregulation indicated a poorer prognosis for BRCA patients. Therefore, cell experiments were conducted on CLDN7 to explore its functions in BRCA. The knockdown and overexpression efficiency of the CLDN7 was verified by RT-PCR and Western blotting. When MDA-MB-231 cells were transfected with si-CLDN7, the mRNA and protein expression levels of CLDN7 were significantly decreased (P<0.05, Figure 8A,8B). At the same time, when MDA-MB-231 cells were transfected with vector-CLDN7 plasmid, the mRNA and protein expression levels of CLDN7 were significantly increased (P<0.05, Figure 8A,8B). In MDA-MB-231 cells, the knockdown of CLDN7 significantly reduced cellular proliferation, migration, and invasion (Figure 8C-8E). Conversely, CLDN7 overexpression significantly promoted the migration and invasion of MDA-MB-231 cells (Figure 8D,8E), although its effects on cell proliferation showed no statistical significance (Figure 8C). Additionally, the effects of CLDN7 on epithelial-mesenchymal transition (EMT)-related proteins were then analyzed. As shown in Figure 8F, CLDN7 overexpression significantly reduced the levels of E-cadherin while elevating the expression of Vimentin (P<0.05). Additionally, single-gene GSEA showed that CLDN7 was closely related to IL2-STAT5, IL6-JAK-STAT3 and TGF-β signal pathways, which were validated in vitro. The results demonstrated that CLDN7 overexpression significantly reduced expression levels of these pathway-related proteins, including p-STAT5, p-STAT3 and p-smad2/3 (P<0.05, Figure 8G).

Figure 8 Functional validation of CLDN7 in MDA-MB-231 cells. The transfection efficiency of the CLDN7 plasmid was verified at both the RT-PCR (A) and western blotting (B) analyses. (C) Effect of CLDN7 on cell proliferation was tested by colony formation assay. (D) Effect of CLDN7 on cell migration was tested by scratch assay; scale bar =500 µm. (E) Effect of CLDN7 on cell migration was tested by Transwell assay in MDA-MB-231 cells; scale bar =200 µm. *, P<0.05 compared with the si-NC; ^#, P<0.05 compared with the vector. (F) Protein expression levels of EMT-related proteins E-cadherin and Vimentin. (G) Protein expression levels of p-STAT5, p-STAT3 and smad2/3. *, P<0.05. Each experiment was performed at least three times. EMT, epithelial-mesenchymal transition; RT-PCR, real-time polymerase chain reaction.

To further verify the robustness of the hub gene network, BATF, TFPI2, and TH were subjected to parallel functional assays. First, the transfection efficiencies of BATF, TFPI2, and TH were verified. As shown in Figure S1A-S1C, siRNA-mediated knockdown significantly reduced the mRNA and protein expression levels of each target gene compared with the corresponding si-NC groups, whereas plasmid-mediated overexpression markedly increased their expression levels relative to vector controls (all P<0.01). These results confirmed the effectiveness and reliability of gene silencing and overexpression strategies. Consistent with the effects observed for CLDN7, silencing of TH significantly suppressed clonogenic ability, migration, and invasion of MDA-MB-231 cells, while BATF or TFPI2 knockdown enhanced cell proliferation, migration and invasion; overexpression of these genes exerted opposite effects (P<0.01, Figure S2A-S2I). These findings indicate that multiple hub genes identified in this study converge on regulating proliferative and metastatic behaviors of BRCA cells.

Discussion

In this study, 7,382 DEGs between BRCA and normal samples were identified, and these genes were enriched in PI3K-Akt and MAPK signal pathways. Previous studies demonstrated that the co-activation of PI3K-Akt and MAPK signal pathways induced tumor malignant features through promoting tumor growth and metastasis (18-20). It was worth noting that the cytokines of these signal pathways also played essential roles in the immune microenvironment through interacting with inflammatory cytokines (21). It indicated that the immune response may be closely related to the development of BRCA. Besides, four hub genes, namely TFPI2, CLDN7, TH and BATF, were screened out through Cox analysis, and a risk model was constructed by these hub genes. The ROC results indicated this risk model exhibited moderate predictive performance, suggesting potential value for prognostic risk assessment.

An increasing number of studies have demonstrated that elevated estrogen levels play a crucial role in the development and progression of BRCA (22). Estrogen signaling has been shown to regulate the transcription of genes encoding non-coding RNAs and to participate in angiogenesis and EMT in BRCA cells (23). Moreover, targeting the estrogen regulatory system has been reported to enhance ferroptosis and endocrine sensitivity in BRCA (24). Emerging evidence also suggests that estrogen signaling can modulate the immune response and influence the efficacy of immunotherapy in BRCA (10). In this study, we constructed a prognostic risk model based on four hub genes and stratified TCGA-BRCA patients into high- and low-risk groups. The high-risk group exhibited significantly lower immune, stromal, and ESTIMATE scores, accompanied by higher tumor purity, suggesting reduced infiltration of immune and stromal cells and a predominance of tumor cells within the tumor microenvironment. Such an “immune-cold” microenvironment is often associated with immune evasion, increased invasiveness, and poor prognosis. Further CIBERSORT analysis showed that CD8⁺ T cells, Tregs, activated NK cells, and M1 macrophages were markedly reduced, whereas M0 and M2 macrophages were significantly elevated in the high-risk group. These findings indicate that the tumor microenvironment of high-risk patients is dominated by immunosuppressive and tumor-promoting cell populations. M2 macrophages can secrete factors such as IL-10 and TGF-β to promote immune tolerance and tumor angiogenesis, thereby impairing the efficacy of immunotherapy (25). Additionally, we observed a significantly lower TP53 mutation frequency in the high-risk group. As a classic tumor suppressor gene, TP53 mutations typically lead to uncontrolled cell proliferation and impaired DNA repair. However, some studies have reported that TP53-mutated BRCA subtypes exhibit higher immunogenicity and increased T-cell infiltration (26,27). Therefore, the reduced TP53 mutation frequency observed in the high-risk group may align with its “immune-cold” phenotype, further suggesting lower immune activity and stronger immune evasion capacity, which may underlie the poor prognosis of these patients. Collectively, these findings suggest that estrogen signaling may influence BRCA progression by orchestrating the immune microenvironment through the synergistic regulation of the four hub genes.

According to the previous studies, the four key genes used to construct the risk model may play vital roles in the development of cancers. TFPI2 is a Kunitz-type serine proteinase inhibitor and prevents the extracellular matrix of cancer cells from degradation and tumor invasion (28). Studies have shown that estrogen 17β-estradiol (E2) induces the expression of TFPI2 through ERα (29). In BRCA, downregulation of TFPI2 has been associated with poor patient prognosis (30). Our analysis further confirmed that TFPI2 is downregulated in BRCA samples and predicts unfavorable outcomes. In addition, TFPI2 can inhibit the EMT process by suppressing TWIST1 expression, thereby reducing BRCA cell migration and invasion (31). Moreover, it has been proven that the overexpressed TWSIT1 could activate the Notch signal pathway to induce the EMT progression (32,33). In our study, we found that downregulated TFPI2 was enriched in the Notch signal pathway and EMT, indicating that TFPI2 may suppress BRCA malignant progression by inhibiting Notch signaling and EMT.

BATF belongs to the activator protein 1 (AP-1) transcription factor family, and it could negatively regulate the AP-1/ATF axis, which controls the differentiation of immune cells (34,35). Besides, some studies identified that the high expression of BATF was involved in the development of tumors by inducing the proliferation and migration of tumor cells (36,37). However, Seo et al. found that the overexpression of BATF in CD8 T cells could induce the production of effector cytokines and improve the anti-tumor function of T CAR T cells through interacting with IRF4 (38). Similarly, our results showed that BATF was also upregulated in BRCA based on bioinformatics analysis and qRT-PCR, and the high expression of BATF was closely related to the good prognosis of BRCA. It indicated that BATF may inhibit the growth of BRCA by enhancing the anti-tumor function of immune cells. However, the infiltration of CD8 T cells was decreased in the high-risk group, suggesting that BATF may not improve the prognosis of BRCA through regulating the immune response. In addition, BATF can act as an inhibitor of AP-1 activity, and the AP-1 transcription factor family regulates cell proliferation, differentiation, and survival (37,39). Combined with the GSEA results and mutation landscape, we inferred that BATF may inhibit the development of BRCA through inducing apoptosis and activating p53 signal pathway.

TH is an enzyme in the family of non-heme iron- and tetrahydrobiopterin (BH4)-dependent aromatic amino acid hydroxylases, which promotes the conversion of L-tyrosine to L-3,4-dihydroxyphenylalanine. Dysregulation of TH is associated with the death of dopaminergic neurons (40). The expression of TH is regulated by estrogen (41). Studies have found that neural elements are present in 85% of BRCA and exhibit cancer-type-specific functions (42). Previous research has shown that activation of the sympathetic nervous system drives BRCA progression in preclinical BRCA models. TH+ sympathetic nerves were detected at the tumor periphery of spontaneous orthotopic BRCA in mice (43). In this study, TH was upregulated in BRCA samples and associated with poor prognosis. Further GSEA indicated that this may be related to the activation of oxidative phosphorylation and the PI3K-AKT-mTOR signaling pathway. Activation of oxidative phosphorylation and PI3K-AKT-mTOR drives BRCA metastasis (44,45). Therefore, we speculate that TH promotes BRCA progression by activating oxidative phosphorylation and the PI3K-AKT-mTOR signaling pathway.

CLDN7 is a member of the claudin family and is important for constituting cell–cell tight junctions (46). As a key protein in epithelial cells, the abnormal expression of CLDN7 usually plays an essential role in the malignance of tumors through mediating EMT progression (47-49). Nguyen et al. identified that the overexpression of CLDN7 could cause the esophageal epithelial barrier dysfunction through regulating the TGF-β signal pathway (50). It indicated the interaction between TGF-β and CLDN7. Moreover, the TGF-β could also participate in the development of tumors through mediating EMT and immune response (51,52). In our results, CLDN7 was remarkably upregulated in BRCA and closely related to the poor prognosis, which is consistent with the findings of Fan et al. (53). The GSEA results also found that CLDN7 was negatively related to the IL2-STAT5, IL6-JAK-STAT3 and TGF-β signal pathways. In vitro experiments validated that CLDN7 overexpression enhanced the aggressive behavior of BRCA cells by promoting invasion, migration and EMT, and it modulated the expression of proteins related to the STAT5, STAT3, and TGF-β pathways. Previous studies have shown that the crosstalk between STAT5 and estrogen signaling influences the proliferation and metastasis of BRCA cells (54). The IL6/STAT3 pathway can hijack a subset of ER enhancers to drive a distinct transcriptional program, thereby promoting BRCA metastasis (55). Moreover, TGF-β plays a crucial role in estrogen-regulated EMT processes (56-58). Therefore, the correlations of CLDN7 with the IL2-STAT5, IL6-JAK-STAT3, and TGF-β signaling pathways suggest that CLDN7 may indirectly participate in the regulation of estrogen signaling, thereby influencing the hormonal responsiveness and invasive potential of BRCA cells.

In addition, drug sensitivity analysis revealed that TFPI2 and CLDN7 were positively correlated with the majority of drugs, whereas TH and BATF were negatively correlated with drug sensitivity. As a tight junction protein, CLDN7 regulates intercellular adhesion and drug permeability, and its aberrant expression may alter drug uptake efficiency. BATF, which is involved in T-cell differentiation and immune regulation, may contribute to drug resistance by modulating the tumor immune microenvironment. Therefore, these genes not only serve as potential prognostic biomarkers but may also act as predictors of chemotherapy and immunotherapy responses. The risk model constructed based on these genes provides a promising tool for predicting patient prognosis and may help clinicians develop individualized therapeutic strategies. By comprehensively analyzing gene expression profiles from the TCGA database and constructing a robust risk model, this study offers important insights into the potential mechanisms underlying BRCA progression.

However, there are notable limitations in this study. The risk model was validated only using data from the TCGA database, and its generalizability in clinical settings remains to be tested.

Conclusions

In conclusion, this study identified four hub genes—TFPI2, CLDN7, TH and BATF—that are closely associated with the prognosis of BRCA. A risk model incorporating these genes demonstrated moderate predictive performance, with potential values for clinical decision-making and personalized medicine. The involvement of these genes in critical signaling pathways, including TGF-β, underscores their potential as therapeutic targets. While this model offers preliminary prognostic and therapeutic insights, further validation with diverse cohorts and functional studies is necessary to fully understand the molecular mechanisms underlying BRCA progression and to confirm the clinical applicability of the risk model.

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-1409/rc

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

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

Funding: None.

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1409/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. Xiong X, Zheng LW, Ding Y, et al. Breast cancer: pathogenesis and treatments. Signal Transduct Target Ther 2025;10:49. [Crossref] [PubMed]
2. Gheller A, Tuttle R, Oxenberg J. Breast Cancer Risk, Screening, and Risk Reduction in Young Females. Am Surg 2025;91:1178-87. [Crossref] [PubMed]
3. Bai J, Gao Y, Zhang G. The treatment of breast cancer in the era of precision medicine. Cancer Biol Med 2025;22:322-47. [Crossref] [PubMed]
4. Habibi S, Bahramian S, Zare Jalise S, et al. Novel strategies in breast cancer management: From treatment to long-term remission. Crit Rev Oncol Hematol 2025;211:104715. [Crossref] [PubMed]
5. Adachi M, Yanagizono K, Okano Y, et al. Estradiol protects hair cells from cisplatin-induced ototoxicity via Nrf2 activation. Redox Rep 2023;28:2161224. [Crossref] [PubMed]
6. Lei D, Hong T, Li L, et al. Isobaric tags for relative and absolute quantitation-based proteomics analysis of the effect of ginger oil on bisphenol A-induced breast cancer cell proliferation. Oncol Lett 2021;21:101. [Crossref] [PubMed]
7. Clusan L, Ferrière F, Flouriot G, et al. A Basic Review on Estrogen Receptor Signaling Pathways in Breast Cancer. Int J Mol Sci 2023;24:6834. [Crossref] [PubMed]
8. Kim J, Munster PN. Estrogens and breast cancer. Ann Oncol 2025;36:134-48. [Crossref] [PubMed]
9. Beste ME, Kaunitz AM, McKinney JA, et al. Vaginal estrogen use in breast cancer survivors: a systematic review and meta-analysis of recurrence and mortality risks. Am J Obstet Gynecol 2025;232:262-270.e1. [Crossref] [PubMed]
10. Corti C, Binboğa Kurt B, Koca B, et al. Estrogen Signaling in Early-Stage Breast Cancer: Impact on Neoadjuvant Chemotherapy and Immunotherapy. Cancer Treat Rev 2025;132:102852. [Crossref] [PubMed]
11. Sun X, Zang H, Li S, et al. EFHD1 promotes breast cancer progression through estrogen signaling. Exp Cell Res 2025;450:114620. [Crossref] [PubMed]
12. Xu T, Chu C, Xue S, et al. Identification and validation of a prognostic signature of drug resistance and mitochondrial energy metabolism-related differentially expressed genes for breast cancer. J Transl Med 2025;23:131. [Crossref] [PubMed]
13. Huang G, Yu Y, Su H, et al. Integrating RNA-seq and scRNA-seq to explore the prognostic features and immune landscape of exosome-related genes in breast cancer metastasis. Ann Med 2025;57:2447917. [Crossref] [PubMed]
14. Liu H, Bao H, Zhao J, et al. Establishment and verification of a prognostic immune cell signature-based model for breast cancer overall survival. Transl Cancer Res 2024;13:5600-15. [Crossref] [PubMed]
15. Li M, Wang S, Huang H, et al. Reliable estrogen-related prognostic signature for uterine corpus endometrial carcinoma. Comput Biol Chem 2024;113:108216. [Crossref] [PubMed]
16. Zhang H, Fan MD, Hu Y, et al. Prognostic Assessment and Analysis of Underlying Biological Mechanisms of Prostate Cancer Based on Estrogen-Related Genes. Horm Metab Res 2025;57:273-85. [Crossref] [PubMed]
17. Zeng Y, Ma W, Li L, et al. Identification and validation of eight estrogen-related genes for predicting prognosis of papillary thyroid cancer. Aging (Albany NY) 2023;15:1668-84. [Crossref] [PubMed]
18. Shorning BY, Dass MS, Smalley MJ, et al. The PI3K-AKT-mTOR Pathway and Prostate Cancer: At the Crossroads of AR, MAPK, and WNT Signaling. Int J Mol Sci 2020;21:4507. [Crossref] [PubMed]
19. Stefani C, Miricescu D, Stanescu-Spinu II, et al. Growth Factors, PI3K/AKT/mTOR and MAPK Signaling Pathways in Colorectal Cancer Pathogenesis: Where Are We Now? Int J Mol Sci 2021;22:10260. [Crossref] [PubMed]
20. Morgos DT, Stefani C, Miricescu D, et al. Targeting PI3K/AKT/mTOR and MAPK Signaling Pathways in Gastric Cancer. Int J Mol Sci 2024;25:1848. [Crossref] [PubMed]
21. 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]
22. Larnder AH, Manges AR, Murphy RA. The estrobolome: Estrogen-metabolizing pathways of the gut microbiome and their relation to breast cancer. Int J Cancer 2025;157:599-613. [Crossref] [PubMed]
23. Mangani S, Piperigkou Z, Koletsis NE, et al. Estrogen receptors and extracellular matrix: the critical interplay in cancer development and progression. FEBS J 2025;292:1558-72. [Crossref] [PubMed]
24. Cao J, Zhou T, Wu T, et al. Targeting estrogen-regulated system x(c)(-) promotes ferroptosis and endocrine sensitivity of ER+ breast cancer. Cell Death Dis 2025;16:30. [Crossref] [PubMed]
25. Xu J, Ding L, Mei J, et al. Dual roles and therapeutic targeting of tumor-associated macrophages in tumor microenvironments. Signal Transduct Target Ther 2025;10:268. [Crossref] [PubMed]
26. Liu Z, Jiang Z, Gao Y, et al. TP53 Mutations Promote Immunogenic Activity in Breast Cancer. J Oncol 2019;2019:5952836. [Crossref] [PubMed]
27. Kim JY, Jung J, Kim KM, et al. TP53 mutations predict poor response to immunotherapy in patients with metastatic solid tumors. Cancer Med 2023;12:12438-51. [Crossref] [PubMed]
28. Ota Y, Uomoto M, Koizume S, et al. Tissue factor pathway inhibitor-2 inhibits integrin β1 activation and focal adhesion formation and suppresses peritoneal ovarian cancer dissemination in mice. Biochem Biophys Res Commun 2024;736:150890. [Crossref] [PubMed]
29. Andresen MS, Ali HO, Myklebust CF, et al. Estrogen induced expression of tissue factor pathway inhibitor-2 in MCF7 cells involves lysine-specific demethylase 1. Mol Cell Endocrinol 2017;443:80-8. [Crossref] [PubMed]
30. Hou J, Liu W, Yan M, et al. Unveiling heterogeneity and prognostic markers in ductal breast cancer through single-cell RNA-seq. Cancer Cell Int 2024;24:266. [Crossref] [PubMed]
31. Zhao D, Qiao J, He H, et al. TFPI2 suppresses breast cancer progression through inhibiting TWIST-integrin α5 pathway. Mol Med 2020;26:27. [Crossref] [PubMed]
32. Fahim Y, Yousefi M, Izadpanah MH, et al. TWIST1 correlates with Notch signaling pathway to develop esophageal squamous cell carcinoma. Mol Cell Biochem 2020;474:181-8. [Crossref] [PubMed]
33. Benjakul N, Prakobphol N, Tangshewinsirikul C, et al. Notch signaling regulates vasculogenic mimicry and promotes cell morphogenesis and the epithelial-to-mesenchymal transition in pancreatic ductal adenocarcinoma. PLoS One 2022;17:e0279001. [Crossref] [PubMed]
34. Fu Y, Koh B, Kuwahara M, et al. BATF-Interacting Proteins Dictate Specificity in Th Subset Activity. J Immunol 2019;203:1989-98. [Crossref] [PubMed]
35. Zhang R, Miao J, Zhai M, et al. BATF promotes extramedullary infiltration through TGF-β1/Smad/MMPs axis in acute myeloid leukemia. Mol Carcinog 2024;63:1146-59. [Crossref] [PubMed]
36. Feng Y, Pan L, Zhang B, et al. BATF acts as an oncogene in non-small cell lung cancer. Oncol Lett 2020;19:205-10. [Crossref] [PubMed]
37. Zhang Z, Lin M, Wang J, et al. Calycosin inhibits breast cancer cell migration and invasion by suppressing EMT via BATF/TGF-β1. Aging (Albany NY) 2021;13:16009-23. [Crossref] [PubMed]
38. Seo H, González-Avalos E, Zhang W, et al. BATF and IRF4 cooperate to counter exhaustion in tumor-infiltrating CAR T cells. Nat Immunol 2021;22:983-95. [Crossref] [PubMed]
39. Murphy TL, Tussiwand R, Murphy KM. Specificity through cooperation: BATF-IRF interactions control immune-regulatory networks. Nat Rev Immunol 2013;13:499-509. [Crossref] [PubMed]
40. Pajarillo E, Rizor A, Son DS, et al. The transcription factor REST up-regulates tyrosine hydroxylase and antiapoptotic genes and protects dopaminergic neurons against manganese toxicity. J Biol Chem 2020;295:3040-54. [Crossref] [PubMed]
41. Zhou L, Poon CC, Wong KY, et al. Prenylflavonoid Icariin Induces Estrogen Response Element-Independent Estrogenic Responses in a Tissue-Selective Manner. J Endocr Soc 2020;4:bvz025. [Crossref] [PubMed]
42. Li D, Hu LN, Zheng SM, et al. High nerve density in breast cancer is associated with poor patient outcome. FASEB Bioadv 2022;4:391-401. [Crossref] [PubMed]
43. Szpunar MJ, Belcher EK, Dawes RP, et al. Sympathetic innervation, norepinephrine content, and norepinephrine turnover in orthotopic and spontaneous models of breast cancer. Brain Behav Immun 2016;53:223-33. [Crossref] [PubMed]
44. Cheng S, Wan X, Yang L, et al. RGCC-mediated PLK1 activity drives breast cancer lung metastasis by phosphorylating AMPKα2 to activate oxidative phosphorylation and fatty acid oxidation. J Exp Clin Cancer Res 2023;42:342. [Crossref] [PubMed]
45. Yu X, Zhang Y, Xiong S, et al. Omics analyses of a somatic Trp53(R245W/+) breast cancer model identify cooperating driver events activating PI3K/AKT/mTOR signaling. Proc Natl Acad Sci U S A 2022;119:e2210618119. [Crossref] [PubMed]
46. West JJ, Golloshi R, Cho CY, et al. Claudin 7 suppresses invasion and metastasis through repression of a smooth muscle actin program. J Cell Biol 2024;223:e202311002. [Crossref] [PubMed]
47. Li Y, Gong Y, Ning X, et al. Downregulation of CLDN7 due to promoter hypermethylation is associated with human clear cell renal cell carcinoma progression and poor prognosis. J Exp Clin Cancer Res 2018;37:276. [Crossref] [PubMed]
48. Ji H, Ding X, Zhang W, et al. Claudin-7 Inhibits Proliferation and Metastasis in Salivary Adenoid Cystic Carcinoma Through Wnt/β-Catenin Signaling. Cell Transplant 2020;29:963689720943583. [Crossref] [PubMed]
49. Hämäläinen L, Soini Y, Pasonen-Seppänen S, et al. Alterations in the expression of EMT-related proteins claudin-1, claudin-4 and claudin-7, E-cadherin, TWIST1 and ZEB1 in oral lichen planus. J Oral Pathol Med 2019;48:735-44. [Crossref] [PubMed]
50. Nguyen N, Fernando SD, Biette KA, et al. TGF-β1 alters esophageal epithelial barrier function by attenuation of claudin-7 in eosinophilic esophagitis. Mucosal Immunol 2018;11:415-26. [Crossref] [PubMed]
51. Su J, Morgani SM, David CJ, et al. TGF-β orchestrates fibrogenic and developmental EMTs via the RAS effector RREB1. Nature 2020;577:566-71. [Crossref] [PubMed]
52. Peng H, Shen J, Long X, et al. Local Release of TGF-β Inhibitor Modulates Tumor-Associated Neutrophils and Enhances Pancreatic Cancer Response to Combined Irreversible Electroporation and Immunotherapy. Adv Sci (Weinh) 2022;9:e2105240. [Crossref] [PubMed]
53. Fan X, Qi A, Zhang M, et al. Expression and clinical significance of CLDN7 and its immune-related cells in breast cancer. Diagn Pathol 2024;19:113. [Crossref] [PubMed]
54. Oueslati M, Bettaieb I, Ben Younes R, et al. STAT-5 and STAT-6 in Breast Cancer: Potential Crosstalk With Estrogen and Progesterone Receptors Can Affect Cell Proliferation and Metastasis. J Clin Med Res 2022;14:416-24. [Crossref] [PubMed]
55. Siersbæk R, Scabia V, Nagarajan S, et al. IL6/STAT3 Signaling Hijacks Estrogen Receptor α Enhancers to Drive Breast Cancer Metastasis. Cancer Cell 2020;38:412-23.e9. [Crossref] [PubMed]
56. Fan DM, Qi PW, Gao SG, Chen YW, Cheng XL. TGF-β1 mediates estrogen receptor-induced epithelial-to-mesenchymal transition in some tumor lines. Tumour Biol 2014;35:11277-82. [Crossref] [PubMed]
57. Yoriki K, Mori T, Kokabu T, et al. Estrogen-related receptor alpha induces epithelial-mesenchymal transition through cancer-stromal interactions in endometrial cancer. Sci Rep 2019;9:6697. [Crossref] [PubMed]
58. Liu BW, Yu ZH, Chen AX, et al. Estrogen receptor-α-miR-1271-SNAI2 feedback loop regulates transforming growth factor-β-induced breast cancer progression. J Exp Clin Cancer Res 2019;38:109. [Crossref] [PubMed]