Tumor senescence leads to poor survival and therapeutic resistance in human breast cancer

Abstract

Background:

Breast cancer (BRCA) is the most common malignant tumor that seriously threatens the health of women worldwide. Senescence has been suggested as a pivotal player in the onset and progression of tumors as well as the process of treatment resistance. However, the role of senescence in BRCA remains unelucidated.

Methods:

The clinical and transcriptomic data of 2994 patients with BRCA were obtained from The Cancer Genome Atlas and the METABRIC databases. Consensus clustering revealed senescence-associated subtypes of BRCA patients. Functional enrichment analysis explored biological effect of senescence. We then applied weighted gene co-expression network analysis (WGCNA) and LASSO regression to construct a senescence scoring model, Sindex. Survival analysis validated the effectiveness of Sindex to predict the overall survival (OS) of patients with BRCA. A nomogram was constructed by multivariate Cox regression. We used Oncopredict algorithm and real-world data from clinical trials to explore the value of Sindex in predicting response to cancer therapy.

Results:

We identified two distinct senescence-associated subtypes, noted low senescence CC1 and high senescence CC2. Survival analysis revealed worse OS associated with high senescence, which was also validated with patient samples from the National Cancer Center in China. Further analysis revealed extensively cell division and suppression of extracellular matrix process, along with lower stromal and immune scores in the high senescence CC2. We then constructed a 37 signature gene scoring model, Sindex, with robust predictive capability in patients with BRCA, especially for long time OS beyond 10 years. We demonstrated that the Sene-high subtype was resistant to CDK inhibitors but sensitive to proteosome inhibitors, and there was no significant difference in paclitaxel chemotherapy and immunotherapy between patients with different senescence statuses.

Conclusions:

We reported senescence as a previously uncharacterized hallmark of BRCA that impacts patient outcomes and therapeutic response. Our analysis demonstrated that the Sindex can be used to identify not only patients at different risk levels for the OS but also patients who would benefit from some cancer therapeutic drugs.

Introduction

Breast cancer (BRCA) is the most common malignant tumor and the leading cause of cancer-related death in women, making it a serious threat to human health (1). Due to the multifaceted nature of this complex disease, patients with BRCA are often classified into luminal A, luminal B, human epidermal growth factor receptor 2 (HER2)-enriched, and basal-like subtypes based on molecular expression by pathology (2, 3). Moreover, the Tumor, Nodes, Metastases (TNM) staging system, based on anatomical abnormalities including tumor size, lymph node involvement, and distant metastatic status, is the most frequently used tool for outcome prediction (4). However, there is still considerable heterogeneity in the treatment responses and the clinical outcomes among patients with similar clinical and pathologic backgrounds (5), highlighting the need for novel predictive markers other than clinical stages and pathohistological classifications.

As next-generation sequencing technologies can provide more information at a now acceptable cost, classification of patients with BRCA based on transcriptomic profiles has come to use, including the Prediction Analysis of Microarray 50 (PAM50), the 21-gene assay (OncotypeDX) and the 70-gene signature (MammaPrint) (6–8). Beside these classification systems based on the whole expression profile, there are also several signatures based on key biological features of the tumor, such as autophagy-related signatures (9, 10), N6-methyladenosine (11, 12), immune cell infiltration (13, 14), etc. These feature-specific signatures bring forth the hope for precision medicine, but many important biological processes remain uninvestigated.

Cellular senescence is the growth arrest of cells that have been intrinsically and/or extrinsically damaged by factors including oncogenic activation, mitochondrial dysfunction, radiation damage, oxidative and genotoxic stress, and chemotherapeutic agent-induced damage (15). Accumulating evidence suggests that cellular senescence is a double-edged sword in human cancer. On one hand, cellular senescence is considered tumor-suppressive by inhibition of cell division and tumor expansion (16–18). On the other hand, senescence-associated cellular plasticity and stemness reprogramming may be critical for treatment-resistance in many cancer types, including BRCA (19–21). Moreover, secretion of senescent cells (inflammatory cytokines, chemokines, growth factors), known as the senescence-associated secretory phenotype (SASP), influences cancerous, stromal, and immune cells in the tumor microenvironment (TME) (22–24). Despite numerous preclinical studies in cellular and animal models, our understanding on the variable effects of senescence in different cancers is still rudimentary, and the characteristics of senescence in patients with BRCA are particularly complex (25–27). Therefore, it is essential to elucidate the value of tumor senescence as a biomarker to guide clinical prognosis and treatment for better cancer prevention and therapy.

In this study, we comprehensively analyzed the senescent features in multiple BRCA cohorts and classified BRCA patients into two subtypes with distinct senescent status, microenvironment composition, mutation frequency, and signaling pathway activation. We also constructed the Senescence Index (Sindex) model, a tumor senescence scoring system that not only predicted the survival outcome of patients with BRCA, but also suggested potentially effective treatment strategies for the high-risk, high-senescence group of patients (Figure 1).

Figure 1

Materials and methods

Data collection and processing

We acquired publicly available datasets and information of patients with BRCA from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) and the Molecular Taxonomy of Breast Cancer International Consortium (METABRIC) databases (cBioPortal, https://www.cbioportal.org/). For the TCGA dataset, raw count expression data was downloaded from the TCGA Data Portal with the TCGAbiolinks R package (28), and then the raw counts were transformed into normalized count with the VST function from the Deseq2 R package for subsequent analysis (29). Microarray gene-expression data was downloaded from the METABRIC database (https://www.cbioportal.org/study/summary?id=brca_metabric). The resources used in this study were summarized in Supplementary Table 4. All the R packages used in this study run in the R software (version 4.1.2, https://www.r-project.org). Data from the clinical trial I-SPY2 and IMvigor 210 was deposited under GSE196096 and GSE145281 (30, 31).

Patient population and clinical information

In this study, we retrospectively collected 2994 patients with BRCA for subsequent analysis: the METABRIC cohort and TCGA cohort each included 1904 and 1090 patients with detailed clinical information, survival follow-up, and complete genomic data, respectively (Supplementary Table 5). In addition, we collected pathological tissue sections and clinical data of patients with BRCA from the National Cancer Center of China. This study was approved by the ethics committee of the Chinese Academy of Medical Sciences and Peking Union Medical College (Ref. 19/327-2111).

Identification of differentially expressed genes and GSEA analysis

The DEGs between BRCA and normal breast tissues were identified with the “Deseq2” R package based on raw count data from the TCGA dataset (29). Genes with a false discovery rate < 0.05 and the |fold change| > 1 were defined as significant DEGs. Gene set enrichment analysis (GSEA) was applied to elucidate biological differences between BRCA and normal adjacent tissue (NAT) (32).

Consensus clustering for the senescence subtypes

Based on BRCA-senescence related genes, patients with BRCA were clustered with “ConsensusClustering Plus” R package (33). Single-sample GSEA (ssGSEA) was performed for senescence related Reactome gene sets using “GSVA” R package (34), and the results were visualized with the R package “ComplexHeatmap” (35). Principal component analysis (PCA) plot was visualized using “PlotPCA” function of the Deseq2 R package (29).

Characterize biological effects of senescence

We first used GSEA to compare difference in biological pathway activities between CC1 and CC2; then calculated the activity of hallmark pathways with ssGSEA, and investigated each pathway’s correlation with senescence. The tumor purity scores, immune infiltration levels, and stromal contents in different samples were evaluated via the “ESTIMATE” algorithm (36). To evaluate stromal and immune cell infiltration characteristics of BRCA, we used the “Xcell” R package (37) to quantitatively analyze the infiltration levels of different immune cells and fibroblasts. The ssGSEA algorithm was used to quantify the activity of immunogenic cell death, apoptosis, TGF-b, extracellular matrix (ECM) assembly, and angiogenesis, all gene sets retrieved from GO-BP. For WES data, the R package “Maftools” was utilized to calculate tumor mutation burden (TMB) and compare mutation differences between patients with high and low senescent statuses (38).

Multiplex immunofluorescence immunohistochemistry and imaging

A tissue microarray (TMA) spotted with tumor samples from 74 BRCA patients (HBreD090PG01) was purchased from Shanghai Outdo Biotech Co. Ltd. A total of 90 cores on the slide consisted of 74 cases of BRCA tissue and paired normal breast tissue. All tissues were collected in accordance with the ethical standards with the donor being informed completely and with their consent, from National Human Genetic Resources Sharing Service Platform: 2005DKA21300. For multiplex immunofluorescence staining, the manual multiplex immunofluorescence immunohistochemistry kit was used according to the manufacturer’s instructions, and the molecules panel, which consisted of five antibodies including anti-PANCK, anti-α-SMA, anti-CD8, anti-CD68, anti-P16, was conducted on the same slide.

WGCNA and Sindex construction and validation

All model construction work was first established in the METABRIC dataset, and then tested in TCGA as external validation. We first used the “WGCNA” package in R (with default parameters) to identify gene modules associated with high senescence, defining the module with the highest absolute module significance as the key module (red). To develop a senescence-based signature of individual tumor with better clinical utility and practicability, we first used the “Ezcox” R package to apply univariate proportional hazards (Cox) regression analysis to preliminarily screen overall survival (OS)-related genes in the module. Then, to remove the multicollinearity among these genes, we applied least absolute shrinkage and selection operator (LASSO) regression with the optimal penalty parameter and a minimum 10-fold cross-validation to identify the most valuable prognostic genes (39). The final 37 genes were used to establish a final model with Cox regression, and the linear predictor values were defined as the Sindex, based on which survival probabilities were calculated with the equation from Cox regression. Time-dependent receiver operating characteristic (ROC) curves were made using survivalROC R package, conducting calibration at different timepoints to determine the robustness of the model. Moreover, the Sindex of patients were used to divided patients into high and low subgroups according to an optimal cutoff value selected by the R package “maxstat” (40), and generated survival curves via the Kaplan-Meier (K-M) method and log-rank test with the R packages “survival” and “survminer” (41). Collectively, area under the curve (AUC) of the ROC curve and K-M survival analysis were considered when evaluating the predictive ability of the model.

For nomogram construction, we first used univariant Cox regression to identify survival-related clinical characteristics (including Sindex). The independent predictive characteristics were screened out with multivariable Cox regression model, and visualized with a forest plot drawn by the survminer R package. The stepwise multivariate Cox regression was used to construct the final nomogram, which included lymph nodes, tumor size, Nottingham prognostic index, inferred menopausal state, age and Sindex. The nomogram was further evaluated for predictive probability of 5-, 10-, and 15-year survival in patients with BRCA using the “survival” R package (41).

Exploration of the significance of senescore in response to cancer therapy

To evaluate the potential value of Sindex in therapeutic suggestion for patients with BRCA, we calculated the half-maximal inhibitory concentration (IC50) of common therapeutic drugs based on the Genomics of Drug Sensitivity in Cancer (GDSC) databases through the Oncopredict algorithm (42). The IC50 of these chemotherapeutic drugs in high and low Sindex subgroups were compared by Wilcoxon test with the results exhibited in box diagrams using the “ggpubr” R package (43). I-SPY2 is an ongoing multicenter, phase II neoadjuvant platform trial for patients with BRCA that used a variety of treatments and combinations (30). Based on the clinical response and transcriptomic data, we calculated Sindex for each patient, and compared the average Sindex between pathologic complete response (pCR) and non-pCR patients treated with each treatment schedule. IMvigor210 is a phase III clinical trial of atezolizumab in urothelial cancer patients. We calculated the Sindex for each patient and compared the average Sindex between responders (R) [complete response (CR), partial response (PR)] and non-responders (NR) [stable disease (SD), progressive disease (PD)] to elucidate the effect of senescence on immunotherapy response.

Cell lines and reagents

In this study, we used the breast cancer cell line MDA-MB-436 obtained from ATCC and pancreatic cancer cell line SW1990 with P53 knockout kindly provided by Prof Ziwen Liu. Both cell lines were grown in the growth medium recommended by ATCC, Leibovitz’s L-15 medium (Thermo Fisher, 11415064) with1% penicillin and streptomycin (Thermo Fisher, 15070063), and 10% and 20% fetal bovine serum (Thermo Fisher, 10099), respectively. Cell cultures are regularly tested for mycoplasma infection with MycoAlert mycoplasma detection kits (Lonza, LT07).

Cloning of TP53 vectors

Wildtype and mutant TP53 vectors were amplified from SW1990 (TP53WT) and PANC1 (TP53MUT, c.818G>A, p.R273H) cells, respectively, as previously described by Y Liu and W Bodmer (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1327731/). The TP53 vectors are then cloned into pInducer20 (Addgene) through the Gateway method: briefly, 50 ng of pCR8 TP53 and 150 ng of pInducer20 plasmids (Addgene) were mixed with LR Clonase II (Invitrogen, 11791) and incubated at 25 °C for 1 hour before reaction termination by Proteinase K. The reaction products were then used to transform DH5a bacteria and positive clones were selected on ampicillin agar plate.

P53 re-expression

Cell lines with stable overexpression of p53 WT or MUT were generated with lentivirus transduction. Lentiviruses were produced by transfecting Lenti-HEK-293 cells with 2 μg of p53 construct, 1.5 μg of viral protein R (VPR), and 0.5 μg of vesicular stomatitis Indiana virus G protein (VSVG), using Lipofectamine 3000 transfection reagent (Thermo Fisher, L3000015). Viral particles were collected from the medium supernatant after centrifuging at 4000 rpm for 30 min. 48 h after infection with 1ml of lentivirus and 10 μg of Polybrene, puromycin was added to select and generate cell lines with stable expression. Stable cells were tested for expression of P53 WT or MUT through sequencing. Transient re-expression was achieved simply through Lipofectamine 3000 aided plasmid transfection according to the manufacturer’s instruction.

Quantification of mRNA expression

Total RNA was obtained with RNeasy kit (Qiagene 74004) following manufacturer’s instruction. cDNA was synthesized with the M-MLV reverse transcriptase kit (Thermo Fisher, 28025013). RT-qPCR reaction was performed in a final volume of 20 μ containing 12 μl TaqMan^® Universal PCR Master Mix, 5 μl H2O, 1 μl of forward and reverse primers, and 1 μl cDNA (approximately 10ng/μl). The reaction was put in an ABI PRISM^® 96-Well Optical Reaction Plate under the standard thermal cycling conditions by ABI PRISM^® 7000 Sequence Detection System (TaqMan^®): initial 50°C for 2 min and 95°C for 10 min followed by 40 cycles at 95°C for 15 sec and 60°C for 1 min were used. All reactions were performed in three duplicates. The primers used for qPCR reactions are shown in Table 1.

mIHC staining

Microarray tissue samples collected from breast cancer patients were used as experimental samples and tonsil tissues were used as controls (both as FFPE samples). All the tissues were cut and made as section slides with 2-μm thicknesses. The slides were deparaffinized in xylene for 10 mins and repeat three times, and rehydrated in absolute ethyl alcohol for 5 mins and repeat twice, 95% ethyl alcohol for 5 mins, 75% ethyl alcohol for 2 mins, sequentially. Then the slides were washed with distilled water 3 times. A microwave-oven is used for heat-induced epitope retrieval; during epitope retrieval, the slides were immersed in boiling EDTA buffer (Alpha X Bio, Beijing, China) for 15mins. Antibody Diluent/Block (Alpha X Bio, Beijing, China) was used for blocking. The mIHC staining part was performed and analyzed according to a 6-plex-7-color panel, and specifications (with primary antibodies used) are as the following: CD8 (ZA0508, ZSGB-BIO, CHINA), CD68 (ZM0060, ZSGB-BIO,CHINA), PANCK (ZM0069, ZSGB-BIO,CHINA), P16 (ab51243, Abcam, Cambridge, UK), a-SMA (ab7817, Abcam, Cambridge, UK) and CD31 (ab76533, Abcam, Cambridge, UK). All the primary antibodies were incubated for 1 hr at 37°C. Then slides were incubated with Alpha X Ploymer HRP Ms+Rb (Alpha X Bio, Beijing, China) for 10 mins at 37°C. Alpha X 7-Color IHC Kit (Alpha X Bio, Beijing, China) was used for visualization. The correspondences between primary antibodies and fluorophores are listed in below:

AlphaTSA 480 (CD8), AlphaTSA 520 (CD68), AlphaTSA 570 (PANCK), AlphaTSA 620 (P16), AlphaTSA 690 (CD31) AlphaTSA 780 (a-SMA). After each cycle of staining, heat-induced epitope retrieval was performed to remove all the antibodies including primary antibodies and secondary antibodies. The slides were counter-stained with DAPI for 5 mins and enclosed in Antifade Mounting Medium (I0052; NobleRyder, Beijing, China). Axioscan7 (ZEISS, Germany) was used for imaging the visual capturing. Data analysis was performed with Halo (3.4, Indica Labs, United States).

Statistics

Statistical analysis was performed based on R software v4.1.2 (https://www.r-project.org/) and GraphPad Prism v9.3.0 (https://www.graphpad.com/). Categorical variables were analyzed using χ2 test or Fisher’s exact test. Continuous variables for paired samples were analyzed using Student’s t test. Multiple groups of continuous variables were analyzed using one-way ANOVA. Survival analysis was performed based on the univariate and multivariate Cox regression. Pearson coefficient of correlation was calculated to measure the correlation between two variables. Unless stated otherwise, two-tailed p < 0.05 was regarded as statistically significant.

Results

Senescence-based subtyping of BRCA

We first compared the transcriptomic data of BRCA and paired NAT from the TCGA dataset. GSEA showed extensive cellular senescence in BRCA but not NAT, in terms of biological pathways including cellular senescence (NES = 2.1709, P < 0.001), DNA damage/telomere stress induced senescence (NES = 2.42879, P < 0.001), formation of senescence associated heterochromatin foci (NES = 1.759044, P < 0.001), oxidative stress induced senescence (NES = 2.244097, P < 0.001), senescence associated secretory phenotype (NES = 2.274088, P < 0.001) and senescence TP53 targets DN (NES = 2.123696, P < 0.001) (Figure 2A).

Figure 2

To further investigate the role of senescence in breast cancer, we reviewed available literatures and databases for senescence-related genes (SRGs) published in the Kyoto Encyclopedia of Genes and Genomes database, the Gene Ontology database, the Reactome database, the Molecular Signatures Database and the Csgene database. Collectively, we found 28 gene sets and 957 SRGs (Supplementary Table 1). We then used Deseq2 to obtain 6236 DEGs between BRCA and NAT, among which 3729 were upregulated and 2505 were downregulated (Figure 2B). To identify SRGs with key biological roles in breast cancer, we overlapped the DEGs and the SRGs to get 187 upregulated and 96 downregulated genes, making up the list of 283 BRCA-SRGs (Figure 2C).

We then investigated the value of senescence in breast cancer subtyping. Consensus clustering of the TCGA dataset based on the expression of BRCA-SRGs classified patients into two distinct senescence-associated clusters, noted CC1 and CC2 respectively (Figure 2D). ssGSEA analysis of widely used Reactome senescence-related gene sets suggested that senescence is activated in CC2 but suppressed in CC1, except for SASP (Figure 2E). We further investigated the role of other patient parameters in the above subtyping scheme, and found no differential distribution of age, menopausal stage, tumor stage, and pathological subtype between the two subtypes. Since senescence is often considered an age-related process, we further confirmed that there was no correlation between senescence and patient age (Supplementary Figure 1A). However, the PAM50 subtypes were unevenly distributed, with more Luminal B, HER2-enriched, and triple-negative subtypes in CC2, and predominantly Luminal A in CC1 (Figure 2E, Supplementary Figure 1B). Similar results were obtained from the METABRIC dataset (Supplementary Figures 1C–F).

Survival analysis between CC1 and CC2 revealed worse OS associated with CC2 in both TCGA and METABRIC (Figures 2F, G). Moreover, PCA based on transcriptomic data revealed clear separation, suggesting distinct molecular characteristics and biological behaviors between the two subtypes (Figure 2H, Supplementary Figure 1G).

Biological effect of senescence

We next elucidated the impact of senescence on the behaviors of BRCA, through which high senescence mediated poor survival. We performed GSEA enrichment analysis between CC1 and CC2 in both METABRIC (Figure 3A) and TCGA (Supplementary Figure 2A), revealing extensively cell division and suppression of ECM process in CC2. As to the hallmark characteristics of tumor, correlation analysis between senescence and hallmark pathways in both METABRIC (Figure 3B) and TCGA (Figure 3C) revealed positive correlation with cell division and replication, in line with previous results. Surprisingly, senescence was negatively correlated with estrogen response and UV response DN.

Figure 3

We next explored how senescence influence the tumor microenvironment (TME), employing ESTIMATE analysis for stromal and immune score, Xcell analysis for fibroblasts and immune cells, and ssGSEA analysis for ECM pathway activation, TGF-b signaling activation, angiogenesis score, immunogenic cell death and apoptosis. For the stromal compartment, CC2 had significantly suppressed stromal development in both TCGA (Figure 3D) and METABRIC (Supplementary Figure 2B). For the immune compartment, CC2 had lower immune score and lower infiltration of CD8⁺T cells and macrophage (Figure 3E, Supplementary Figure 2C). CC1 and CC2 had comparable TMB (Supplementary Figures 2D, E). As different forms of cell death have profound effects on the surrounding microenvironment, we determined that CC2 underwent less immunogenic cell death and apoptosis, potentially causing the dissimilar TME (Figure 3E, Supplementary Figure 2C). Interestingly, further analysis of paired whole-exon sequencing (WES) from TCGA and METABRIC revealed high mutation frequency of TP53 in CC2 (48%) but not CC1 (11%) (Figure 3F, Supplementary Figure 2F). Previous studies reported that the activated p53 modulates senescence with a dual effect of promoting or, in some cases, inhibiting the senescence program. Re-expression of P53^WT in the TP53 truncated cell line MDA-MB-436 breast cancer cell line significantly inhibited senescence, while P53^R175H slightly increased senescence compared to P53^null (Figure 3G, Supplementary Figure 2G). The results have been verified by transient re-expression of P53WT versus P53^R273H in the P53KO SW1990 cell line derived from human pancreatic cancer (Supplementary Figure 2H).

Validation of senescence at the protein level

To further validate the relationship between senescence and TME at translational level, a sophisticated method of multiplex immunofluorescence histochemistry, which allows simultaneous detection of multiple target proteins, was employed to analyze the protein expression of the five genes, PANCK, α-SMA, CD8, CD68, P16, identified above on a BRCA tissue microarray. The five proteins were the marker of epithelial cell, fibroblast, CD8⁺T cell, macrophage and senescence, respectively. High expression of P16 corresponds to low expression of PANCK, α-SMA, CD8, and CD68, suggesting low stromal development and low infiltration of immunocyte (Figures 4A–C).

Figure 4

To better understand the status of senescence and prognosis for BRCA patients, we included 30 BRCA patients in the National Cancer Center in China who have pathological tissue. According to the expression level of P16 determined by immunofluorescence, we divided the patients into high and low senescence groups by midline expression. Survival analysis showed higher OS in low senescence group (log-rank test, P = 0.029; median OS 87.0 months versus 36.9 months, Figure 4D).

Construction of a senescence scoring model

On a population-based view, senescence markedly affected the prognosis of BRCA patients. Therefore, we constructed a subtyping method, Sindex that allowed characterization of senescence status and subsequent outcome prediction for each patient. We first applied weighted gene co-expression network analysis (WGCNA) to all genes in the METABRIC dataset. After adjustments of WGCNA parameters, the starting genes were divided into 17 modules by average linkage hierarchical clustering (Figures 5A, B). The magenta module (506 genes) exhibited the highest correlation with CC2 (Pearson’s correlation coefficient = 0.71, P < 0.001) (Figures 5B, C), collectively termed CC2 or senescence-high related genes (SHRGs). Next, the 506 SHRGs were primarily screened by univariate Cox regression, leaving 278 OS-correlated SHRGs for further analysis (Supplementary Table 5). Enrichment analysis of the SHRGs demonstrated significant correlation with cellular senescence and cell division (Figure 5D). Finally, the LASSO regression algorithm screened out 37 genes for Sindex calculation (Figures 5E, F, Supplementary Table 2).

Figure 5

Validation of Sindex to effectively predict patient outcomes based on senescence

Sindex was first internally validated with the METABRIC dataset, yielding an impressive predicting ability for OS (5-year AUC = 0.680, 10-year AUC = 0.675, 15-year AUC = 0.680, 20-year AUC = 0.726, 25-year AUC = 0.789) as indicated by the time-dependent ROC curve analysis (Figure 6A). External validation in the TCGA dataset further confirmed the robustness of Sindex in predicting OS for breast cancer patients (5-year AUC = 0.651, 10-year AUC = 0.632) (Figure 6D). The best range of OS prediction was 5-25 years, suggesting that Sindex was especially good for long-term prognosis, making it useful for BRCA, which is a chronic disease.

Figure 6

Looking at the Sindex of individual patients, an optimal cutoff value was selected with the maximally selected rank statistics from the ‘maxstat’ R package to define high senescence and low senescence subtypes, termed Sene-high and Sene-low respectively (Figures 6B, E). Kaplan-Meier survival analysis demonstrated clear separation of the two survival curves and shorter OS of the Sene-high subtype in METABRIC (log-rank test, P < 0.0001; median OS 19.6 months versus 56.2 months, Figure 6C) as well as TCGA (log-rank test, P < 0.001; median OS 19.6 months versus 56.2 months, Figure 6F).

Furthermore, we used univariant Cox regression to identify all OS-related characteristics with significant regression coefficients and p-values (Supplementary Table 3), and then used multivariate Cox regression to determine that Sindex is an independent and robust prognostic factor (HR, 2.097; 95% CI, 1.751–2.511; P < 0.001; Figure 7A). Based OS-related characteristics, we constructed a nomogram (Figure 7B) that could accurately predict the probability of 5-, 10-, and 15-year survival for BRCA patients (Figures 7C–E). The Sindex, lymph nodes positivity, tumor size, Nottingham prognostic index, inferred menopausal state, and age were incorporated into the nomogram as related predictors of breast cancer patient survival.

Figure 7

Senescence and drug resistance

Using the Oncopredict algorithm for drug response prediction, we demonstrated that the Sene-high(subtype were resistant to CDK inhibitors (CDK9i-AZ5576 and Dinaciclib) (Figures 8A, B) and Epirubicine (Figure 8C), but there was no significant difference in sensitivity to proteasome inhibitors (MG132 and Bortezomib) between Sene-low and Sene-high (Figures 8D, E). Further analysis revealed that Bortezomib had the lowest IC50 and the most negative correlation with Sindex (Figure 8F), suggesting that Bortezomib may be the most effective drug for the group of high-senescence and high-risk BRCA patients. Interestingly, the Sene-high subtype was resistant to some anti-microtubule chemotherapeutic drugs (Supplementary Figures 3A, B), but there was no significant difference in paclitaxel sensitivity between the two groups, suggesting more complex effect-response-resistance mechanisms of these drugs (Figure 8G).

Figure 8

We further explored the real-world effect of senescence on drug response in transcriptomic data from I-SPY2, a clinical research of neoadjuvant treatment for BRCA patients. There was no significant difference in Sindex bewtween pCR and non-pCR patients receiving paclitaxel, in line with previous result (Figures 8G, H). However, the significantly higher Sindex in non-pCR patients receiving veliparib and carboplatin suggested resistance to this treatment strategy conferred by senescence (Figure 8I). Therapeutic response to Neratinib (TKI), MK2206 (AKT inhibitor), Ganitumab (anti-IGF1R-mAb), and Trebananib (neutralizing peptide against Ang1/2) are not significantly affected by senescence, but non pCR patients to the Hsp90 inhibitor Ganetespib have significantly higher Sindex (Supplementary Figures 4A, B, D–F). As to patients receiving TDM1 and pertuzumab, pCR patients had significantly higher Sindex (Figure 8J), but patients with pCR and non-pCR to pertuzumab had comparable Sindex (Supplementary Figure 4C), suggesting that the high-senescence and high-risk group of patients may benefit from T-DM1 treatment. We noted that high-senescence patients may not benefit from Pembrolizumab immunotherapy (Figure 8K), consistent with previous results suggesting lower immune infiltration into senescent tumors. This was further confirmed by data from urothelial cancer patients treated with Atezolizumab in IMvigor 210 (Figure 8L).

Discussion

In 2020, the incidence of BRCA reached 11.7%, surpassing lung cancer for the first time to become the most diagnosed cancer in the world (1). Despite widely used classification system by molecular markers, the heterogeneous nature of BRCA poses difficulty onto prognosis and therapeutic decision (2, 44). Therefore, exploring new methods of BRCA subtyping has become the most urgent task in bringing precision medicine into reality (45). In this study, we identified senescence as a hallmark of BRCA that is associated with distinctive TME and genomic alterations. Based on the senescent characteristics of BRCA, we developed a subtyping scheme Sindex that provides an insightful perspective on the biological property of tumors from a new perspective.

We detected that high senescence is associated with significantly worse survival outcome in BRCA patients, suggesting that the activation of cellular senescence pathways may promote malignancy and therapeutic resistance. Although senescence is often considered a defense mechanism against cancer due to its cell cycle arresting properties (46), senescence is a potentially reversible process of epigenetic and transcriptional changes, and senescent dormancy may offer protection by immune-evasion, drug resistance, and resistance to apoptosis or other forms of cell death (47). Sun et al. reported that senescence is associated with good survival in gastric cancer (48), suggesting multimodal function and different mechanisms of senescence in different cancers.

The complexity of senescence deserves further investigation on the potential biological effect of senescence in BRCA. Different methods of data analysis in two datasets both reached the same conclusions: (1) senescence is associated with extensive cell replication and suppression of stroma and immune development; (2) high senescence suppresses stromal development and immune infiltration in the TME. SASP is a well-investigated property of senescence generally considered to be pro-tumorigenic, as senescent tumor cells have specific secretory phenotypes marked by increased secretion of pro-inflammatory (but not anti-tumor immunity), bioactive molecules including IL-6, IL-1a, chemokines (CXCL8), VEGF, and proteases (49). However, our result did not reveal a significant difference in SASP pathway activity between the Sene-high and the Sene-low groups, suggesting other mechanisms by which senescence remodels the TME. We further investigated other key characteristics crucial for anti-tumoral immunity, as cell death and antigen release are crucial for the initiation of anti-tumoral immune response, we indeed found out that senescence suppresses immunogenic cell death, resulting in lower antigen-presenting cell infiltration.

To increase the clinical application value and create better clinical practicability, we successfully constructed a novel senescence-related scoring tool (Sindex) to determine the prognostic risk of BRCA. The Sindex effectively stratified patients with BRCA into high- and low-risk groups. Survival analysis revealed that the Sene-high group had shorter OS than Sene-low group, and ROC curves exhibited a great predictive capacity of Sindex for the 5-, 10-, 15-, 20- and 25-year survival of BRCA. The model is validated in two individual datasets, METABRIC and TCGA.

Previous studies have demonstrated that cellular senescence could increase the drug resistance and side effects of the chemotherapy (50–52). Using the Oncopredict algorithm and I-SPY2 dataset, we found out that indeed, most drugs have a higher IC50 value in the Sene-high group, especially CDK inhibitors, epirubicin, and olaparib. Such resistance is explainable as those drugs share common mechanisms in inhibition of cell cycle progression through cell death induced damage. On the other hand, patients with high senescence may benefit from metabolic drugs such as proteosome inhibitors and MTORC inhibitors, consistent with the fact that senescent cells are metabolically active. Using the I-SPY2 and IMvigor 210 data sets, we noted that high-senescence patients may not benefit from immune checkpoint therapies, consistent with previous results showing lower immune infiltration in Sene-high subtype.

Our study has the limitation that the high-throughput sequencing data sets for initial analysis were relatively insufficient as it was simply obtained from public databases. The corresponding results and conclusions remain to be further investigated through more external congeneric research and should be validated via functional experiments in vivo and in vitro. Furthermore, several conclusions of this study require further research to confirm its reproducibility, improve the clinical application of senescence-related clusters, and elaborate on the role of Sindex in predicting the response to cancer therapy for BRCA.

Conclusions

In conclusion, our study lead to the identification of senescence subtypes that are associated with significantly different TME and survival outcomes. Our analysis demonstrated that the Sindex can be used to identify not only patients at different risk levels for the OS but also patients who would benefit from some cancer therapeutic drugs. Nevertheless, the validation of our findings in a wide spectrum of patient cohorts, and the findings that the senescence features reflect biological and clinical characteristics associated with sensitivity or resistance to the therapy, would pave a way for developing more rational therapy recommendations and promoting personalized cancer therapy.

References

• 1

  SungHFerlayJSiegelRLLaversanneMSoerjomataramIJemalAet al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA: Cancer J Clin (2021) 71:209–49. doi: 10.3322/caac.21660

  • CrossRef
  • Google Scholar

• 2

  PerouCMSørlieTEisenMBvan de RijnMJeffreySSReesCAet al. Molecular portraits of human breast tumours. Nature (2000) 406:747–52. doi: 10.1038/35021093

  • CrossRef
  • Google Scholar

• 3

  SørlieTPerouCMTibshiraniRAasTGeislerSJohnsenHet al. Gene expression patterns of breast carcinomas distinguish tumor subclasses with clinical implications. Proc Natl Acad Sci United States America (2001) 98:10869–74. doi: 10.1073/pnas.191367098

  • CrossRef
  • Google Scholar

• 4

  ParkYHLeeSJChoEYLa ChoiYLeeJENamSJet al. Clinical relevance of TNM staging system according to breast cancer subtypes. Ann Oncol Off J Eur Soc Med Oncol (2019) 30:2011. doi: 10.1093/annonc/mdz223

  • CrossRef
  • Google Scholar

• 5

  JiangYLiTLiangXHuYHuangLLiaoZet al. Association of adjuvant chemotherapy with survival in patients with stage II or III gastric cancer. JAMA Surg (2017) 152:e171087. doi: 10.1001/jamasurg.2017.1087

  • CrossRef
  • Google Scholar

• 6

  SorlieTTibshiraniRParkerJHastieTMarronJSNobelAet al. Repeated observation of breast tumor subtypes in independent gene expression data sets. Proc Natl Acad Sci United States America (2003) 100:8418–23. doi: 10.1073/pnas.0932692100

  • CrossRef
  • Google Scholar

• 7

  PetkovVIMillerDPHowladerNGlinerNHoweWSchusslerNet al. Breast-cancer-specific mortality in patients treated based on the 21-gene assay: a SEER population-based study. NPJ Breast Cancer (2016) 2:16017. doi: 10.1038/npjbcancer.2016.17

  • CrossRef
  • Google Scholar

• 8

  CardosoFvan't VeerLJBogaertsJSlaetsLVialeGDelalogeSet al. 70-gene signature as an aid to treatment decisions in early-stage breast cancer. New Engl J Med (2016) 375:717–29. doi: 10.1056/NEJMoa1602253

  • CrossRef
  • Google Scholar

• 9

  WuQLiQZhuWZhangXLiH. Identification of autophagy-related long non-coding RNA prognostic signature for breast cancer. J Cell Mol Med (2021) 25:4088–98. doi: 10.1111/jcmm.16378

  • CrossRef
  • Google Scholar

• 10

  LinQGLiuWMoYZHanJGuoZXZhengWet al. Development of prognostic index based on autophagy-related genes analysis in breast cancer. Aging (2020) 12:1366–76. doi: 10.18632/aging.102687

  • CrossRef
  • Google Scholar

• 11

  ZhongXLiJWuXWuXHuLDingBet al. Identification of N6-Methyladenosine-Related LncRNAs for predicting overall survival and clustering of a potentially novel molecular subtype of breast cancer. Front Oncol (2021) 11:742944. doi: 10.3389/fonc.2021.742944

  • CrossRef
  • Google Scholar

• 12

  ZhangXShenLCaiRYuXYangJWuXet al. Comprehensive analysis of the immune-oncology targets and immune infiltrates of n (6)-Methyladenosine-Related long noncoding RNA regulators in breast cancer. Front Cell Dev Biol (2021) 9:686675. doi: 10.3389/fcell.2021.686675

  • CrossRef
  • Google Scholar

• 13

  WangSXiongYZhangQSuDYuCCaoYet al. Clinical significance and immunogenomic landscape analyses of the immune cell signature based prognostic model for patients with breast cancer. Briefings Bioinf (2021) 22. doi: 10.1093/bib/bbaa311

  • CrossRef
  • Google Scholar

• 14

  DenkertCvon MinckwitzGDarb-EsfahaniSLedererBHeppnerBIWeberKEet al. Tumour-infiltrating lymphocytes and prognosis in different subtypes of breast cancer: a pooled analysis of 3771 patients treated with neoadjuvant therapy. Lancet Oncol (2018) 19:40–50. doi: 10.1016/s1470-2045(17)30904-x

  • CrossRef
  • Google Scholar

• 15

  SiebenCJSturmlechnerIvan de SluisBvan DeursenJM. Two-step senescence-focused cancer therapies. Trends Cell Biol (2018) 28:723–37. doi: 10.1016/j.tcb.2018.04.006

  • CrossRef
  • Google Scholar

• 16

  SunYCoppéJPLamEW. Cellular senescence: The sought or the unwanted? Trends Mol Med (2018) 24:871–85. doi: 10.1016/j.molmed.2018.08.002

  • CrossRef
  • Google Scholar

• 17

  Pérez-ManceraPAYoungARNaritaM. Inside and out: the activities of senescence in cancer. nature reviews. Cancer (2014) 14:547–58. doi: 10.1038/nrc3773

  • CrossRef
  • Google Scholar

• 18

  ArrealLPivaMFernándezSRevandkarASchaub-CleriguéAVillanuevaJet al. Targeting PML in triple negative breast cancer elicits growth suppression and senescence. Cell Death differentiation (2020) 27:1186–99. doi: 10.1038/s41418-019-0407-5

  • CrossRef
  • Google Scholar

• 19

  Triana-MartínezFLozaMIDomínguezE. Beyond tumor suppression: Senescence in cancer stemness and tumor dormancy. Cells (2020) 9:346. doi: 10.3390/cells9020346

  • CrossRef
  • Google Scholar

• 20

  SmigielJMTaylorSEBrysonBLTamagnoIPolakKJacksonMW. Cellular plasticity and metastasis in breast cancer: a pre- and post-malignant problem. J Cancer metastasis Treat (2019) 5:47. doi: 10.20517/2394-4722.2019.26

  • CrossRef
  • Google Scholar

• 21

  MongiardiMPPellegriniMPalliniRLeviAFalchettiML. Cancer response to therapy-induced senescence: A matter of dose and timing. Cancers (2021) 13:484. doi: 10.3390/cancers13030484

  • CrossRef
  • Google Scholar

• 22

  CoppéJPPatilCKRodierFSunYMuñozDPGoldsteinJet al. Senescence-associated secretory phenotypes reveal cell-nonautonomous functions of oncogenic RAS and the p53 tumor suppressor. PloS Biol (2008) 6:2853–68. doi: 10.1371/journal.pbio.0060301

  • CrossRef
  • Google Scholar

• 23

  JacksonJGPantVLiQChangLLQuintás-CardamaAGarzaDet al. p53-mediated senescence impairs the apoptotic response to chemotherapy and clinical outcome in breast cancer. Cancer Cell (2012) 21:793–806. doi: 10.1016/j.ccr.2012.04.027

  • CrossRef
  • Google Scholar

• 24

  XuCXieNSuYSunZLiangYZhangNet al. HnRNP F/H associate with hTERC and telomerase holoenzyme to modulate telomerase function and promote cell proliferation. Cell Death differentiation (2020) 27:1998–2013. doi: 10.1038/s41418-019-0483-6

  • CrossRef
  • Google Scholar

• 25

  HeSSharplessNE. Senescence in health and disease. Cell (2017) 169:1000–11. doi: 10.1016/j.cell.2017.05.015

  • CrossRef
  • Google Scholar

• 26

  HerranzNGilJ. Mechanisms and functions of cellular senescence. J Clin Invest (2018) 128:1238–46. doi: 10.1172/jci95148

  • CrossRef
  • Google Scholar

• 27

  ChildsBGGluscevicMBakerDJLabergeRMMarquessDDananbergJet al. Senescent cells: an emerging target for diseases of ageing. nature reviews. Drug Discovery (2017) 16:718–35. doi: 10.1038/nrd.2017.116

  • CrossRef
  • Google Scholar

• 28

  ColapricoASilvaTCOlsenCGarofanoLCavaCGaroliniDet al. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res (2016) 44:e71. doi: 10.1093/nar/gkv1507

  • CrossRef
  • Google Scholar

• 29

  LoveMIHuberWAndersS. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol (2014) 15:550. doi: 10.1186/s13059-014-0550-8

  • CrossRef
  • Google Scholar

• 30

  WolfDMYauCWulfkuhleJBrown-SwigartLGallagherRILeePREet al. Redefining breast cancer subtypes to guide treatment prioritization and maximize response: Predictive biomarkers across 10 cancer therapies. Cancer Cell (2022) 40:609–623.e606. doi: 10.1016/j.ccell.2022.05.005

  • CrossRef
  • Google Scholar

• 31

  BalarAVGalskyMDRosenbergJEPowlesTPetrylakDPBellmuntJet al. Atezolizumab as first-line treatment in cisplatin-ineligible patients with locally advanced and metastatic urothelial carcinoma: a single-arm, multicentre, phase 2 trial. Lancet (London England) (2017) 389:67–76. doi: 10.1016/s0140-6736(16)32455-2

  • CrossRef
  • Google Scholar

• 32

  SubramanianATamayoPMoothaVKMukherjeeSEbertBLGilletteMAet al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci (2005) 102:15545–50. doi: 10.1073/pnas.0506580102

  • CrossRef
  • Google Scholar

• 33

  WilkersonMDHayesDN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinf (Oxford England) (2010) 26:1572–3. doi: 10.1093/bioinformatics/btq170

  • CrossRef
  • Google Scholar

• 34

  HänzelmannSCasteloRGuinneyJ. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinf (2013) 14:7. doi: 10.1186/1471-2105-14-7

  • CrossRef
  • Google Scholar

• 35

  GuZHübschmannD. Make interactive complex heatmaps in r. Bioinf (Oxford England) (2021) 38:1460–2. doi: 10.1093/bioinformatics/btab806

  • CrossRef
  • Google Scholar

• 36

  YoshiharaKShahmoradgoliMMartínezEVegesnaRKimHTorres-GarciaWet al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun (2013) 4:2612. doi: 10.1038/ncomms3612

  • CrossRef
  • Google Scholar

• 37

  NewmanAMLiuCLGreenMRGentlesAJFengWXuYet al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods (2015) 12:453–7. doi: 10.1038/nmeth.3337

  • CrossRef
  • Google Scholar

• 38

  MayakondaALinDCAssenovYPlassCKoefflerHP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res (2018) 28:1747–56. doi: 10.1101/gr.239244.118

  • CrossRef
  • Google Scholar

• 39

  RanstamJCookJA. Statistical models: an overview. Br J Surg (2016) 103:1047. doi: 10.1002/bjs.10240

  • CrossRef
  • Google Scholar

• 40

  SullivanLMMassaroJMD'AgostinoRBSr. Presentation of multivariate data for clinical use: The framingham study risk score functions. Stat Med (2004) 23:1631–60. doi: 10.1002/sim.1742

  • CrossRef
  • Google Scholar

• 41

  RanstamJCookJA. Kaplan-Meier Curve. Br J Surg (2017) 104:442. doi: 10.1002/bjs.10238

  • CrossRef
  • Google Scholar

• 42

  YangWSoaresJGreningerPEdelmanEJLightfootHForbesSet al. Genomics of drug sensitivity in cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Res (2013) 41:D955–961. doi: 10.1093/nar/gks1111

  • CrossRef
  • Google Scholar

• 43

  WhiteheadMJMcCanneyGAWillisonHJBarnettSC. MyelinJ: an ImageJ macro for high throughput analysis of myelinating cultures. Bioinf (Oxford England) (2019) 35:4528–30. doi: 10.1093/bioinformatics/btz403

  • CrossRef
  • Google Scholar

• 44

  KravchenkoJAkushevichISeewaldtVLAbernethyAPLyerlyHK. Breast cancer as heterogeneous disease: contributing factors and carcinogenesis mechanisms. Breast Cancer Res Treat (2011) 128:483–93. doi: 10.1007/s10549-011-1347-z

  • CrossRef
  • Google Scholar

• 45

  CardosoFKyriakidesSOhnoSPenault-LlorcaFPoortmansPRubioITet al. Early breast cancer: ESMO clinical practice guidelines for diagnosis, treatment and follow-up†. Ann Oncol Off J Eur Soc Med Oncol (2019) 30:1194–220. doi: 10.1093/annonc/mdz173

  • CrossRef
  • Google Scholar

• 46

  ColladoMSerranoM. Senescence in tumours: evidence from mice and humans. nature reviews. Cancer (2010) 10:51–7. doi: 10.1038/nrc2772

  • CrossRef
  • Google Scholar

• 47

  WangLLankhorstLBernardsR. Exploiting senescence for the treatment of cancer. nature reviews. Cancer (2022) 22:340–55. doi: 10.1038/s41568-022-00450-9

  • CrossRef
  • Google Scholar

• 48

  ZhouLNiuZWangYZhengYZhuYWangCet al. Senescence as a dictator of patient outcomes and therapeutic efficacies in human gastric cancer. Cell Death Discovery (2022) 8:13. doi: 10.1038/s41420-021-00769-6

  • CrossRef
  • Google Scholar

• 49

  CoppéJPDesprezPYKrtolicaACampisiJ. The senescence-associated secretory phenotype: the dark side of tumor suppression. Annu Rev Pathol (2010) 5:99–118. doi: 10.1146/annurev-pathol-121808-102144

  • CrossRef
  • Google Scholar

• 50

  DemariaMO'LearyMNChangJShaoLLiuSAlimirahFet al. Cellular senescence promotes adverse effects of chemotherapy and cancer relapse. Cancer Discovery (2017) 7:165–76. doi: 10.1158/2159-8290.cd-16-0241

  • CrossRef
  • Google Scholar

• 51

  SalehTTyutyunyk-MasseyLGewirtzDA. Tumor cell escape from therapy-induced senescence as a model of disease recurrence after dormancy. Cancer Res (2019) 79:1044–6. doi: 10.1158/0008-5472.can-18-3437

  • CrossRef
  • Google Scholar

• 52

  MantovaniFCollavinLDel SalG. Mutant p53 as a guardian of the cancer cell. Cell Death differentiation (2019) 26:199–212. doi: 10.1038/s41418-018-0246-9

  • CrossRef
  • Google Scholar