A prognostic gene model with psychological stress-related genes captures immune activation in breast cancer.

1
Introduction
Breast cancer is the most prevalent malignancy in women worldwide, with approximately 2.3 million new cases reported in 2022 [1]. Although the overall survival of breast cancer has significantly improved due to the advancement in early diagnosis and various therapeutic strategies, disease burden including recurrence and metastasis continues to pose substantial challenges [2].
Accumulating evidence highlights the critical role of psychosocial factors, especially chronic stress and depression, in modulating oncogenesis and disease progression [3]. Epidemiologic results indicate that depression is associated with an increased risk of breast cancer [4]. A large-scale cohort study also revealed that 42 % of breast cancer patients developed comorbid mental disorders within a 4-week period, which ranks first among the investigated solid tumors and correlates with poor clinical outcomes [3,5]. These findings suggest a bidirectional relationship between depression and breast cancer.
Mechanistically, chronic stress and depression trigger activation of the hypothalamic-pituitary-adrenal axis and sympathetic nervous system, resulting in sustained release of glucocorticoids and catecholamines [6]. These molecules can influence tumor microenvironment (TME), reducing CD8+ T cell infiltration while increasing M2 macrophage polarization, and activate downstream oncogenic pathways, such as MAPK/ERK and STAT3, thereby promoting proliferation, angiogenesis, metastasis, and immune evasion [[7], [8], [9]]. Notably, genetic polymorphisms and epigenetic modifications of psychological stress-related genes, such as GPR54 and TSC22D3, can regulate neuroendocrine responses and influence tumor behavior. Expression levels of these genes may also act as prognostic factors for cancer patients [10,11]. Collectively, previous studies indicate that psychological stress can affect the expression of critical genes within tumor cells, and these psychological stress-related genes may have good prognostic predictive potential.
Using data from The Cancer Genome Atlas (TCGA) and the Molecular Taxonomy of Breast Cancer International Consortium (METABRIC), we constructed a 18-gene signature of literature reported psychological stress-related genes, and this signature showed good prognosis predicting efficacy. This study provides novel insights into the interplay between psychological stress related tumor genes and breast cancer progression, providing strategies to breast cancer risk stratification and diseases prognosis.

2
Methods
2.1
Patient data acquisition
TCGA gene expression data (raw counts and transcript per million) were downloaded from Gene Expression Omibus (GEO, accession number GSE62944) [12]. Survival and PAM50 intrinsic subtypes of TCGA were retrived from Liu et al. and our previous study, respectively [13,14]. Gene expression data from the METABRIC were downloaded from the Synapse software platform (syn1688369, Sage Bionetwors, Seattle, WA, USA) [15]. Gene expression data of 107 TNBC and 12 immunotherapy treated metastatic TNBC were downloaded from GEO with accession number GSE58812 and GSE241876 [16,17]. H&E based estimation of tumor infiltrating lymphocytes were retrieved from Saltz et al. [18].

2.2
PubMed search of psychological stress related genes
A systematic literature search (Supplementary Fig. 1, Supplementary Table 1) in PubMed was conducted with the following search keywords: “stress, psychological”, “psychological stress”, “psychosocial stress”, “breast neoplasms”, “mammary neoplasms, experimental”. Time filter of last 10 years (from 2015/01/01 to 2024/12/31) was used in our search. The following keywords were also used to identify psychological stress-induced tumor gene expression signatures in other cancer types: “stress, psychological”, “psychological stress”, “psychosocial stress”, “neoplasms”, “gene expression”. Article which met all the following criteria was include in the study: (1) studying psychological stress and cancer; (2) reporting mental psychological stress-related gene expression signature; (3) the gene signature was tumor gene expression signature. Any article which met one or more of the following criteria was excluded from the study: (1) review, protocol, case report, comment, editorial and retracted articles; (2) no gene involved in the study; (3) not studying psychological stress or tumor; (4) studying the psychological stress caused by BRCA mutation test; (5) not reporting tumor gene expression signature induced by psychological stress; (6) reporting gene expression signature in PBMC/blood/saliva but not in tumors.

2.3
Differential expression and pathway analysis
Differential expression analysis was performed using R package DESeq2 (version 1.44.0) and limma (version 3.62.1) for TCGA and METABRIC data, respectively [19,20]. Gene set enrichment analysis (GSEA) was performed using R package ClusterProfiler (version 4.14.3) with 50 hallmark gene sets downloaded from the Molecular Signature Database (MsigDB, version 7.5.1) [21]. GSEA default settings were used except the following parameters: eps = 0, seed = 12345, pvalueCutoff = 1. Pathway analysis in the 37 survival related genes was performed using the enrichGO function in ClusterProfiler and the enrichment of Gene Ontology (GO) biological processes were analyzed. T cell inflamed gene expression signature score for each tumor was calculated with R package GSVA (version 2.0.5).

2.4
Immune cell infiltration analysis
The proportions of different immune cells in TCGA and METABRIC datasets were estimated using CIBERSORTx (https://cibersortx.stanford.edu) with default settings [22]. TCGA breast cancer cellular composition predicted by Kassandra cell deconvolution tool were retrieved from Zaitsev et al. [23].

2.5
Connectivity map (CMap) analysis
Genes with absolute fold change>1.5 and adjusted p-value (FDR) < 0.05 were defined as differentially expressed genes and used as input in the CMap query (https://clue.io) to identify drugs which induce the most similar (or dissimilar) gene expression pattern [24]. Touchstone database gene expression data (L1000, version 2020-12-17) were queried, and the results in TNBC cell line MDA-MB-231 were reported in this study.

2.6
Statistical analysis
All statistical analyses were performed using R software (version 4.4.1). The logistic least absolute shrinkage and selection operator (LASSO) regression was conducted with R package glmnet (version 4.1–8). Z-score normalized gene expression of the 37 overall survival-related genes and the overall survival in METBARIC basal-like breast cancer were used as input. A 10-fold cross validation was used in the LASSO regression with seed 20241027. “Lamda.1se” was used to build the simplest model with comparable accuracy with the best model. 18 genes were selected to construct the final Cox proportional-hazards regression model using R package survival (version 3.7-0). R packages timeROC (version 0.4) and rms (version 6.8-2) were used to plot the time-dependent receiver operating characteristic (ROC) curves and the calibration plot. R package ggscidca (version 0.2.3) was used to plot the decision curves. Log-rank test was used for two-group survival comparison in Kaplan-Meier plot. All p-values were two-sided and p-value<0.05 was considered as significant in this study.

3
Results
3.1
Identification of survival associated genes which were affected by stress and depression
We first retrieved 374 genes whose expression in tumor was influenced by psychological stress (Supplementary Table 2) from the literature review [[25], [26], [27], [28]]. We then used TCGA and METABRIC dataset to identify genes related to the overall survival in TNBC and/or METABRIC (Fig. 1). 15 and 23 genes (Fig. 2A) were significantly related to the overall survival in TCGA and METABRIC basal-like breast cancers, respectively. Pathway analysis (Fig. 2B–Supplementary Table 3) with the 37 survival related genes showed that the top 3 enriched pathways were related to hypoxia, which shares overlapping physiological responses and molecular pathway changes with psychological stress [29].

3.2
Establishment and evaluation of an 18-gene psychological stress-related signature to predict overall survival
With the 37 survival related genes, we next performed LASSO regression in METABRIC basal-like breast cancers (Fig. 3A), and 18 genes were selected to construct the final Cox proportional-hazards regression model to predict overall survival (Supplementary Table 4). We carried out ROC curve, calibration curve and decision curve analyses to evaluate the performance of the prognosis prediction model. Patients with high 18-gene signature score had significantly worse prognosis in METABRIC and external validation datasets TCGA and GSE58812 (Fig. 3B, E and 3H). The AUCs were 0.764, 0.646 and 0.659 for predicting 5-year overall survival in METABRIC, TCGA and GSE58812, and were slightly lower than the 10-year AUC of 0.772, 0.703 and 0.768 (Fig. 3C, F, 3I). The calibration plot showed good agreement between the prediction by the model and the actual observed 5-year overall survival in METABRIC, TCGA and GSE58812 datasets (Fig. 3D, Supplementary Fig. 2), with a Hosmer-Lemeshow goodness-of-fit test p-value = 1. The decision curve analysis suggested that the 18-gene prediction model was of value in clinical decision making (Fig. 3G, Supplementary Fig. 2).

3.3
Immune-related pathways were activated in tumors with low 18-gene signature score
To further understand the connection between low 18-gene signature score and good prognosis, we performed differential expression (DE) and pathway analysis comparing the breast cancers with low and high 18-gene signature scores in METABRIC and TCGA (Supplementary Tables 5–8). GSEA showed that “Hallmark allograft rejection”, “Hallmark interferon alpha response” and “Hallmark interferon gamma response” were activated in the low-signature score group (Fig. 4A and C). “Hallmark estrogen response late”, which was found to be related to lower CD8+ T cell infiltration and worse response to immunotherapy in our previous study [30], was downregulated in the low-signature score group (Fig. 4B and D). Contrary to our expectation, we compared the levels of total tumor infiltration lymphocytes between the two groups and didn't observe difference (Fig. 4E). We further compared the proportion of detailed immune cell types between the low- and high-signature score groups. We found the low-signature score group had significantly higher CIBERSORT estimated CD8+ T cells proportion in both METABRIC and TCGA datasets, whereas the M2 macrophages were significantly lower in the low-signature score group (Fig. 4F and G). Using Kassandra cell deconvolution tool predicted TCGA cellular composition, we also identified significantly higher proportion of CD8+ T cells and NK cells (Supplementary Fig. 3) in the low-signature score group. Besides, we found the low-signature score group also had significantly higher T cell inflamed gene expression score (Fig. 4H), suggesting breast tumors with low-signature score had immune inflamed tumor microenvironment [31]. Multivariate analysis with the T cell inflamed gene expression score and the 18-gene signature score showed the established 18-gene signature score was prognosis to overall survival in METABRIC and TCGA dataset independently of the T cell inflamed score (Fig. 4I).

3.4
Low 18-gene signature score was related to better response to immunotherapy
As the low 18-gene signature score was related to activated immune pathways and higher proportion of CD8+ T cells, we next examined whether tumors with low 18-gene signature score had better response to immunotherapy. Four TME subtypes were previously identified by Bagaev et al. [32], of which the immune-enriched/fibrotic (IE/F) and immune-enriched/non-fibrotic (IE) subtypes had better response to immunotherapy than the firbotic (F) and desert (D) subtypes. We found that TCGA basal-like tumors in the IE/F and IE subtypes had significantly lower 18-gene signature score than that in the F and D subtypes (Fig. 5A). Besides, metastatic TNBC with complete response and partial response to immunotherapy also had significantly lower 18-gene signature score (Fig. 5B). These data suggest low 18-gene signature score was related to better response to immunotherapy.
To identify potential medications that could turn the tumor gene expression in the high 18-gene signature score group towards the pattern in the low-signature score group, we queried CMap Touchstone datasets with the DE genes between the low- and high-signature score groups in METABRIC. 13 drugs which could induce significantly (FDR<0.05) similar or opposing (Fig. 5C–Supplementary Table 9) gene expression patterns as the input DE genes were identified. Epinephrine, a short-term and long term regulator of stress [33], was identified to induce opposing gene expression patterns (negative normalized connectivity score) towards the low 18-gene signature score tumor group. Bopindolol, which blocks beta adrenergic receptors and further inhibits the effects of epinephrine and norepinephrine, could induce similar gene expression patterns (positive normalized connectivity score) like in the low 18-gene signature group.

4
Discussion
Stress and depression are frequently exhibited in breast cancer patients, owing to disease burden like mastectomy, therapeutic adverse effects, and economic problems. Beyond traditional biological drivers, psychosocial disorders—notably chronic stress and depression—are increasingly recognized as pivotal regulators of breast cancer progression [3]. Depression and anxiety are highly prevalent among breast cancer patients, and are correlated with accelerated metastasis and reduced overall survival, while studies on the depression and anxiety induced gene expression changes in tumors were limited [[3], [4], [5]]. Our study identified 374 psychological stress influenced tumor genes from literature review. Using the data in TCGA and METABRIC, we developed and validated a prognostic signature composed of 18 psychological stress influenced tumor genes to stratify patients with different risks and predict prognosis, with high AUC values especially for long-term survival prediction. Notably, the 10-year AUCs were higher than the 5-year AUCs, indicating the genes included in this model may be related to long-term survival mechanisms like immune escape or metabolic adaption. We furthermore performed pathway analysis and identified the 18 gene signature could capture the activation of immune-related pathways.
Currently, limited studies have focused on the role of psychological stress-related genes in cancer prognosis and therapeutic efficacy. Several studies built the breast cancer prognosis or diagnosis model with psychological stress-related gene list, while the genes they used were either from the blood sample or simply reported as depression-related [34,35]. The expression of these genes may not be influenced by psychological stress in tumors. Specifically, Wang et al. established a gene signature with psychological stress related genes that included 10 genes to predict breast cancer overall survival and therapeutic sensitivity, including chemotherapy and endocrine therapy. Their risk model was trained with all breast cancer subtypes. In TNBC, their model had AUCs of 0.730 and 0.609 in internal and external validation groups at 5 years [35]. Compared to their 10-gene risk model, our 18-gene signature was constructed specifically for TNBC and performed better in predicting prognosis (Fig. 3C, F, 3I). Additionally, they found significant different abundances of immunosuppressive cells in different risk groups [35], which is also similar to our results.
In our model, several included genes like CMTM6 [36], NFATC2 [37], and PPP1R15A [38], were proved to be associated with immune microenvironment or immunotherapy response. For instance, overexpression of CMTM6 inhibits the interaction of Hsc70 and PD-L1 that promoting its degradation, resulting in blockade of PD-L1 lysosomal degradation and influence PD-L1-dependent immune evasion [36]. Besides, NFATC2 expression can repress Id3 transcription and modulate PD-1 transcription [37]. These findings suggest our 18-gene score could serve as a potential biomarker for immunotherapy candidacy [39]. Stress and negative emotions can remodel TME through regulating hypothalamic-pituitary-adrenal axis, adrenergic receptor signaling, and epigenetic modification, inducing an immunosuppressive TME [[6], [7], [8],10,40,41]. Our model was constructed with 18 tumor genes whose expression was influenced by psychological stress. Our data here found tumors with low-signature scores exhibited significant activation of immune-related pathways, including allograft rejection and interferon alpha/gamma response (Fig. 4A and C). Previous studies have uncovered that these pathways are hallmarks of cytotoxic immune activation, usually associated with enhanced antigen presentation and T-cell-mediated tumor elimination. For example, interferon gamma, mainly released by activated T cells, can inhibit tumor cell proliferation and trigger apoptosis [42]. Activation of interferon gamma signaling in tumor cells increases the sensitivity to T-cell-mediated killing and immunotherapy [43]. Conversely, the estrogen response late pathway, related to immunosuppressive microenvironment and poor immunotherapeutic response [44], was downregulated in the low-signature score group (Fig. 4B and D). Furthermore, cell-type deconvolution revealed significantly elevated CD8+ T cells and NK cells in the low-score group, while M2 macrophages were significantly lower (Fig. 4F and G). CD8+ T cells and NK cells are crucial anti-tumor cytotoxic cell populations. Elevation of these cell groups can activate the granzyme pathway, leading to extensive apoptosis and triggering the anti-tumor immune response [45]. In contrast, M2 macrophages are modulators of immune remodeling, angiogenesis, T cell suppression, and tumor promotion [46]. Reduction of M2 macrophages could reshape the immunosuppressive TME and enhance the antitumor immune response [47]. These findings showed the 18-gene model successfully captured the immune inflamed tumor microenvironment in breast cancer, and the different activation of immune-modulation pathways may contribute to the distinct clinical outcomes in the high- and low-score groups. Of note, low 18-gene signature was associated with better overall survival independently of the T cell inflamed gene expression score. Multiple non-immune pathways like “Hallmark Epithelial Mesenchymal Transition” and “Hallmark Angiogenesis”, which could also be regulated by chronic stress [[48], [49], [50]], were downregulated in the low 18-gene signature group. This suggested that other mechanisms may also contribute to the better survival in the low 18-gene score group besides the activation of anti-tumor immune response. Further investigations are needed to reveal the underlying mechanisms fully.
The 18-gene signature score model developed in our study holds dual clinical value for basal-like breast cancer. On the one hand, it provides precise risk stratification. High-risk patients may be inclined to benefit from intensified surveillance or novel therapeutic combinations, while low-risk patients characterized by an immune-active TME may be optimal candidates for immunotherapy. Stepped-care psychological interventions based on this signature model could potentially reduce disease burden and healthcare expenditures. On the other hand, identifying potential agents capable of shifting high-signature score tumors toward low-signature score might provide hope for combination therapy with immunotherapy. Some of the identified drugs in our CMap analysis have been reported to modulate immune characteristics and tumor progression. Epinephrine activates adrenergic receptors and the downstream signaling pathways that induce inflammatory injury, immunodysregulation, and cancer progression [40,41]. Previous studies have demonstrated that β-adrenergic receptor signaling inhibition could reverse the immunosuppressive status and improve the sensitivity to immunotherapy [51,52]. Parallelly, we identified beta-adrenergic receptor blockade Bopindolol as a potential drug to improve breast cancer survival outcomes and immunotherapy response. Additionally, p38 MAPK signaling modulates PD-L1 expression levels, tumor associated macrophages reprogramming, and tumor growth in various cancers [[53], [54], [55]]. Doramapimod inhibits the p38 MAPK pathway and may thereby enhance immunotherapy response. Besides, another identified drug Rigosertib, a PLK1/PI3K inhibitor, can downregulate PD-L1 expression and induce cancer cells autophagy through AMPK-ULK 1 activation in colorectal cancer [56], indicating the synergetic effect to immunotherapy. The safety and therapeutic effect of combining Rigosertib with Nivolumab is now evaluating in a phase I/IIa clinical trial (NCT04263090). These agents, with established safety profiles, could be tracked into clinical trials targeting high-signature score tumors, potentially reversing immunosuppressive characteristics and improve clinical outcomes.
This study has some limitations. The three datasets used in our study used different sequencing platforms and had imbalanced data, which may lead to lower prediction ability of the 18-gene signature. This predictive model needs further validation through prospective clinical studies and large-scale retrospective data. Besides, the selection of the 374 psychological stress-related genes was based on the report in previous studies. Only a few studies reported psychological stress influenced tumor genes and thus the selection of the 374 genes may introduce selection bias. Future sequencing study in breast cancer with psychometric data is warranted to obtain psychological stress-related tumor genes more comprehensively, and to validate the direct connection between stress and this 18-gene signature. The 18-gene signature was derived from TNBC cohorts within three datasets which exhibit ethnic homogeneity. The METABRIC and GSE58812 cohorts comprise almost exclusively European-ancestry individuals. While TCGA exhibits greater ethnic diversity with notable (33.4 % in basal-like breast cancer) black and African individuals, white women remain the predominant demographic. This limited ethnic representation indicates that the 18-gene signature requires validation in more diverse population in future study. Additionally, in vivo experiment and studies on the detailed mechanism are needed before the application of the identified compound in clinical use. Future prospective study should also be performed to confirm the synergetic medications that can reverse the high-risk status and sensitize breast cancer cells to immunotherapy.

5
Conclusion
Collectively, the 18-gene signature score model can stratify breast cancer patients into different risk groups and capture the immune activation in the tumor microenvironment. This scoring system may predict immunotherapy response and enhance prognostic precision in breast cancer therapy.

CRediT authorship contribution statement
Tian Du: Writing – review & editing, Writing – original draft, Methodology, Formal analysis, Data curation. Jiayi Tang: Writing – original draft, Formal analysis, Data curation. Hanqi Zhang: Writing – original draft, Formal analysis, Data curation. Qianwen Liu: Writing – review & editing, Methodology, Conceptualization. Yanan Kong: Writing – review & editing, Funding acquisition, Conceptualization.

AI assistance declaration
GPT-4 (OpenAI) and DeepSeek were used exclusively for grammar checking and enhancing readability of non-technical sections. No AI contributed to data synthesis or claim formulation. The authors are responsible for all academic content and clinical insights herein.

Funding
This research was funded by the 10.13039/501100001809National Natural Science Foundation of China (No. 82403103, Yanan Kong). Tian Du was funded by the Science and Technology Project in Guangzhou (Project number 2024A04J4844).

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.