An Aging-Related lncRNA Signature Establishing for Breast Cancer Prognosis and Immunotherapy Responsiveness Prediction

Abstract

Purpose

Emerging evidence demonstrates the vital role of aging and long non-coding RNAs (lncRNAs) in breast cancer (BC) progression. Our study intended to develop a prognostic risk model based on aging-related lncRNAs (AG-lncs) to foresee BC patients’ outcomes.

Patients and Methods

307 aging-related genes (AGs) were sequenced from the TCGA project. Then, 697 AG-lncs were identified by the co-expression analysis with AGs. Using multivariate and univariate Cox regression analysis, and LASSO, 6 AG-lncs, including al136531.1, mapt-as1, al451085.2, otud6b-as1, tnfrsf14-as1, and linc01871, were validated to compute the risk score and establish a risk signature. Expression levels of al136531.1, mapt-as1, al451085.2, tnfrsf14-as1, and linc01871 were higher in low-risk BC patients, whereas otud6b-as1 expression was higher in high-risk BC patients. In the training and testing set, high-risk patients performed shorter PFI, OS, and DFS than low-risk patients.

Results

Our risk signature had the highest concordance index among other established prognostic signatures and displayed ideal predictive ability for 1-, 3- and 5-year patient OS in the nomogram. Additionally, BC patients with different risk score levels showed different immune statuses and responses to immunotherapy via GSEA, ssGSEA, ESTIMATE algorithm, and TIDE algorithm analysis. Of note, the qRT-PCR analysis validated that these 6 AG-lncs expressed quite differentially in BC tissues at various clinical stages.

Conclusion

The risk signature of 6 AG-lncs might offer a novel prognostic biomarker and promisingly enhance BC immunotherapy’s effectiveness.

Keywords:

• risk model
• senescence
• immune
• cancer therapy

Introduction

Globally, the number of newly diagnosed cases (11.7%) of female breast cancer (BC) is about 2.3 million, surpassing lung cancer as the most prevalent malignancy.¹ BC possesses metastatic heterogeneity and plasticity, with varying colonization preferences in distinct organs, such as brain, lung, liver, and bone.² Although there are tremendous advances in therapeutic strategy, given the highly invasive characteristic, it continues to be confronted with undesirable prognoses during recent decades.³ Therefore, the efforts to screen prognostic indicators and build accurate models for anticipating the prognosis are still critical for BC control.

With the update and advancements in genome sequencing technology, various newly emerging molecules have been screened out as significant participants in BC that were not previously appreciated.⁴ Tumor evolution commonly results in abnormally regulated epigenetic events, such as aberrant non-coding RNA (ncRNA) expression levels and DNA methylation. These altered ncRNAs in BC cells can accurately or partially reflect the information originating from primary and metastatic sites.⁵ Thus, it poses the potential prognostic implications of ncRNAs in BC outcomes. For instance, circ_0069094 was upregulated in paclitaxel‐resistant BC tissues and cells, involving in tumor growth and invasion through competitively sponging miR-136-5p.⁶

Pathogenesis and progression of malignant breast tumors have been linked to long noncoding RNAs (lncRNAs). LncRNAs function as key molecular mechanisms to regulate genomic activity, acting as tumor suppressors or oncogenes.⁷ LncRNAs can not be translated into proteins, but directly mediate gene modification, transcription, and translation, performing as competitive endogenous RNAs.⁸ In addition, lncRNAs are involved in multiple signaling pathways, which modify tumor immunity, cellular metabolism, cellular proliferation, and angiogenesis. LncRNAs extensively participate in the oncogenesis and progression of BC. For instance, BC lung metastasis dormancy was induced by the modulation of NR2F1 and Np63 by lncRNA nr2f1-as1, underlining the importance of this regulator in BC stem-like cell plasticity and metastatic dormancy.⁹ In addition, lncRNA SNHG1 could facilitate the capabilities of proliferation, migration, invasion, and EMT in BC cells by activating the miR-641/RRS1 axis.¹⁰ In the xenograft mouse model, SNHG1 intervention suppressed the growth of BC tumors. Ge et al revealed that the expression of metabolism-related lncRNA C6orf99 was remarkably higher in BC cell lines, and the knockdown of C6orf99 inhibited cell proliferation, migration, and invasion in MCF-7, indicating that C6orf99 contributed to BC growth and metastasis.¹¹

With growing interest in the relationship between lncRNAs and BC progression, emerging literature suggests that lncRNAs could serve as potential biomarkers of BC prognosis. LncRNAs could modulate the malignant behavior of BC cells via diverse complicated mechanisms, including binding to proteins, regulating downstream genes, blocking molecules, and competitively binding miRNA with target genes.¹² Stratification of BC risk can be efficiently achieved by integrating polygenic risk scores into prediction models, which improve surveillance and treatment effectiveness.¹³ For example, the upregulated lncRNA afap1-as1 was a well-established biomarker in BC tissue, which was also related to the poor malignancy status of patients.¹⁴ Similarly, Zhao et al pointed out that the developed prognostic risk-related signature, comprising ac004585.1, linc01235, al031316.1, and acta2-as1, accurately distinguished the high- and low-risk patients, showing the capacity of lncRNAs as indicators in BC prognosis and overall survival (OS).¹⁵ HOTAIR is an oncogenic lncRNA that has been recognized as a pivotal modulator in BC development, metastasis, and drug resistance.¹⁶ Interestingly, the serum expression level of HOTAIR exhibited a significant elevation in PR+/ ER+ BC patients, and was notably higher in stage III than in stages I and II, which validated the potential of HOTAIR as a prognostic tool for BC patients.¹⁷

Aging is a complicated biological and pathological phenomenon involving elaborated alterations at the epigenetic and genetic levels.¹⁸ Numerous methods have been proposed to speed up the natural aging process, including loss of telomeres, oxidative stress, DNA damage, and epigenetic alteration.¹⁹ Besides, aging is related to the decrease in immune system competence to combat diseases and prevent tumor formation.²⁰ Various evidence confirms transcriptomic tissue-specific peculiarity in linking human aging and cancer.²¹ Moreover, the potential biological processes underlying cancers are remarkably similar to aging, including exacerbated inflammation, abraded telomeres, and unstable genomes.²² Aging cancer cells can repress tumor immunity, trigger angiogenesis, facilitate proliferative signaling, and promote invasiveness.²³ Aging-related genes (AGs) or ncRNAs are proposed to have underlying predictive capabilities in cancer patients’ prognoses. For example, Xu et al utilized a 6 AGs-based signature (plcg2, mxd1, epor, ywhaz, apoc3, and h2afx) to develop a risk signature in predicting lung adenocarcinoma prognosis, which had a more precise prediction ability than simple clinicopathologic nomograms.²⁴ Also, 7 AGs, including mapt, cdk1, rae1, pola1, eef1e1, socs2, and hdac1, were successfully used to compute a risk score in hepatocellular carcinoma, and patients with high risk significantly possessed a worse prognosis, a lower tumor differentiation but a higher clinical stage.²⁵

Nevertheless, there is little research focusing on the integration of aging-related lncRNAs (AG-lncs) into predictive risk signatures for BC prognosis. Our study meticulously established a co-expression network of aging-associated lncRNAs and mRNAs, leveraging the extensive expression data from 1090 BC samples and 112 normal samples from the Cancer Genome Atlas (TCGA) project. Through univariate and multivariate Cox regression analysis and LASSO regression analysis, we have constructed a risk signature composed of 6 AG-lncs. Moreover, the prognostic value of the risk signature was validated through analytical techniques such as gene set enrichment analysis (GSEA) and single-sample gene set enrichment analysis (ssGSEA) to elucidate the differential immune states between high- and low-risk groups, providing insights into their potential responses to immunotherapy and chemotherapy. Then, the qRT-PCR analysis validated that these 6 AG-lncs expressed quite differentially in BC tissues at various clinical stages. Interestingly, the risk score represented an independent prognostic indicator and could stratify BC patients with diverse clinical outcomes into two risk levels. Notably, high- and low-risk patients exhibited distinct immune statuses and responsiveness to immunotherapy. The whole flowchart of this study is summarized in Figure 1. Our risk signature based on AG-lncs is expected to propose a novel and effective predictor for BC prognosis.

Figure 1 The flowchart of this study.

Materials and Methods

Data Collection and Processing

307 human AGs were downloaded from the human aging genome resource (HAGR) (http://genomics.senescence.info/genes/). The normalized RNA sequencing data (FPKM format) of AGs and the associated clinical characteristics of 1202 samples, including 1090 BC samples and 112 normal samples were obtained from the TCGA data portal (https://portal.gdc.cancer.gov/). Besides, the DNA mutation data and copy number variation (CNV) on all BC patients were also downloaded from the TCGA project. The criteria for inclusion were followed: 1) patients diagnosed as BC; 2) cases with intact clinical information and expression profiles. The criteria for exclusion were followed: 1) patients without survival information or with fewer than 30 d of survival; 2) patients without information on the pathological grade or clinical stage. Table S1 displays the clinical characteristics of the recruited BC patients. All 1090 BC samples were assigned to the training group, whereas 544 BC patients were chosen randomly for the testing group.

Identification of DEAGs

DEAGs (p < 0.05, |log2FC| > 1) were identified through the limma package of R. The protein-protein interaction (PPI) network based on DEGs was established via the Search Tool for the Retrieval of Interacting Genes (STRING) (http://string-db.org/) (threshold combined score set as ≥ 0.4). Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) analyses were performed with the clusterProfiler package of R.

Co-Expression Analysis and Identification of AG-Lncs

The correlation between 307 AGs from HAGR and AG-lncs was identified by Pearson correlation analysis. Filtering of AG-lncs was performed based on the criteria of correlation coefficient > 0.3 and p < 0.001. Candidate AG-lncs that were selected according to the above criteria were further analyzed.

Construction of the Risk Score Based on AG-Lncs

The Lasso Cox regression analysis used the “glmnet” software package. The risk score was calculated according to the following formula:

The risk score = Expression_mRNA1 × Coefficient_mRNA1 + Expression_mRNA2 × Coefficient_mRNA2 +. + Expression_mRNAn × Coefficient_mRNAn

R language v4.0.2 was conducted to make the survival analysis. The receiver operating characteristic (ROC) curve was plotted by the “survivalROC” R package to assess the signature’s specificity and sensitivity for 1-, 2-, and 3-year OS in BC patients.

The Development of Nomogram and GSEA

In the TCGA datasets, the risk score and the other independent risk factors were recruited to develop a nomogram to predict BC patients’ 1-, 2-, and 3-year survival rates by the “rms” R package. The calibration curve was applied to evaluate the effect and accuracy of the nomogram. The GSEA method was applied to determine the latent discrepancies in the biological function and signaling pathways by the clusterProfiler package of R (p < 0.05).

The Immune- and Tumor- Microenvironment Analysis

The ssGSEA was conducted to analyze the infiltrating scores of 16 immune cells and the reactivity of 13 immune-related pathways by the “gsva” R package. Using the ESTIMATE score, ESTIMATE algorithm, stromal score, tumor purity, and the Immune score for each BC sample were identified to evaluate the abundance of immune and stromal cells.

The Sensitivity of the Immune Therapy Response Analysis

The tumor immune dysfunction and exclusion (TIDE) algorithm (http://tide.dfci.harvard.edu/) was performed to evaluate the immunotherapy response. The National Comprehensive Cancer Network (NCCN) guidelines recommend the use of anti-tumor drugs, such as doxorubicin, etoposide, and rapamycin in the treatment of BC. Drug sensitivity was assessed with IC50 values of commonly used chemotherapeutic drugs in BC patients in the TCGA project. The ggplot2 and pRRophetic were used to visualize the results.

Tissue Collection

To evaluate the prognostic potential of the risk score, 20 BC patients at different clinical stages were collected from Sept. 2020 to Oct. 2021 in plastic and cosmetic surgery at Tongji Hospital of Huazhong University of Science and Technology. The corresponding clinicopathologic data of those 20 patients is provided in Table S2. All the participants involved in the study have signed and approved informed consent.

Quantitative Reverse Transcription Polymerase Reaction (qRT-PCR) Analysis

Total RNA was isolated from BC tissues by the TRIzol reagent kit (Invitrogen). The RNA concentration was evaluated by a K5800 spectrophotometer (Kaiao, Beijing, China). Then, using the PrimeScript RT kit (Takara, Japan), the RNA was reversely transferred into complementary DNA (cDNA) after 40 cycles at 103°C for 5 s, 37°C for 10 min, and 4°C for 15 min. The cDNA samples were subjected to qPCR by SYBR GreenTM Master Mix (Yeasen, China). Table S3 provides the primer sequences for the 6 lncRNAs detection. The GAPDH level was set as an internal standard control. All the gene expression levels were gathered and quantified by the 2^−ΔΔCt method. The results were obtained from three independent experiments.

Statistics Analysis

All statistical analysis was conducted by GraphPad Prism (version 8.0) and R version 4.0.5. The continuous variables between the two groups were compared by the independent t-test, and the proportion differences were compared by the χ2-test. The prognostic value was assessed by multivariate and univariate Cox proportional hazards regression analyses. The progression-free interval (PFI), disease-specific survival (DSS), and OS between different groups were compared by Kaplan-Meier analysis through the Log rank test. The TIDE score was compared between groups by the Wilcoxon test. The correlation between the two variables was assessed by Pearson correlation analyses. All statistical tests were performed as two-sided tests, and the values of p < 0.05 were indicated to be significant.

Results

Identification of DEAGs and Construction of a Risk Score Based on AG-Lncs

307 human AGs in total were collected from the HAGR (Table S4). Then, The TCGA project discovered roughly 73 DEAGs, including 33 downregulated and 40 upregulated DEAGs, comparing BC and non-tumor samples (|logFC| > 1 and FDR < 0.05) (Figure 2A). The 73 DEAG interactions were visualized via the PPI (Figure 2B). Furthermore, these DEAGs were considerably abundant in the biological aging process in the GO analysis (Figure 2C). For example, gland development, response to oxidative stress, and response to peptide hormones, which all participated in cellular senescence and cancer progression, were intensely correlated with DEAGs. Likewise, in the KEGG analysis, certain DEAGs were linked to many malignancies, particularly BC, and also associated with PI3K-AKT and JAK-STAT signaling pathways, which were key regulatory pathways in cancer development (Figure 2D).

Figure 2 The expression, PPI, and function analysis of DEAGs and risk model construction. (A) A heatmap of the expression level of 73 DEAGs between the normal (N, blue) and BC tissues (T, red). (B) PPI network visualizing the interactions of DEAGs (interaction score = 0.5). (C) Barplot graph for KEGG analysis of DEAGs. The longer the bar, the more enriched the genes; The deeper the red, the more obvious the differences (q-value: the adjusted p-value). (D) Bubble graph for GO analysis of DEAGs. The bigger the bubble, the more enriched the genes; The deeper the red, the more obvious the differences (q-value: the adjusted p-value). (E and F) LASSO regression analysis of the 6 prognostic AG-lncs. (G) The forest plot of the HR for the correlation between 6 AG-lncs and OS of BC.

The co-expression analysis confirmed 697 AG-lncs. Subsequently, the univariate Cox regression analysis identified 25 prognostic AG-lncs in BC patients’ survival (Table S5). Afterward, we performed the stepwise multivariate Cox regression analysis and LASSO regression analysis to cut down gene counts, thus generating 6 AG-lncs (al136531.1, mapt-as1, al451085.2, otud6b-as1, tnfrsf14-as1, and linc01871) to develop the risk signature (Figure 2E–G). The risk score was calculated by the multivariate Cox regression analysis coefficient. The formula was as follows:

The risk score = al136531.1_expression × (−0.3430) + mapt-as1_expression × (0.3604) + al451085.2_expression × (0.2184) + otud6b-as1_expression × (0.0467) + tnfrsf14-as1_expression × (0.2501) + linc01871_expression × (0.2600). The Kaplan-Meier survival curves validated the prognostic value of these 6 AG-lncs in BC survival (p < 0.001) (Figure S1). Specifically, in accordance with our risk score formula, al136531.1, mapt-as1, al451085.2, tnfrsf14-as1, and linc01871 were higher expressed in low-risk patients, whereas otud6b-as1 was higher expressed in high-risk patients. In other words, the higher expression of otud6b-as1 predicted a lower survival rate of BC patients.

Validation of the Risk Score Based on AG-Lncs

Next, in the training cohort, 1090 BC patients were classified into high-risk (545 patients) and low-risk groups (545 patients) based on the median risk score (Figure 3A). Likewise, 544 BC patients were also classified into high-risk (267 patients) and low-risk groups (277 patients) in the testing cohort (Figure 3B). The high-risk score endowed BC patients with a likely higher death probability. Meanwhile, BC patients in the high-risk group displayed shorter OS (Figure 3C–F), PFI (Figure 3D–G), and DSS (Figure 3E–H) than patients in the low-risk group both in the training cohort and testing cohort. Importantly, both in the training cohort and testing cohort, the multivariate and univariate Cox regression analyses discovered that the risk score was relevant to the poor OS, and was an independent risk predictor in the survival of BC patients, compared with the other established predictors (clinical stage, TNM stage, and age) (Figure S2A–D). Furthermore, the nomogram enrolling all prognostic factors predicted the 1-, 3- and 5-year survival possibility in the training and testing cohort (Figure 4A and B). The calibration plot corresponded well to the practical observation, reflecting the optimal performance of the nomogram (Figure 4C and F). Besides, the AUC values of 1-, 2- and 3-year OS were 0.710, 0.716, and 0.724 in the training cohort (Figure 4D), and 0.671, 0.706, and 0.710 in the testing cohort (Figure 4G). More notably, the AUC value of the risk score at 3-year OS was the highest compared to other prognostic factors, indicating the risk score’s potential to predict the long-term survival rate in BC patients (Figure 4E and H).

Figure 3 Identification of prognostic aging-related lncRNA signature in training and testing cohort. The distribution and survival status of BC patients with different risk scores in training cohort (A) and testing cohort (B). The green and red dots represented survival and death, respectively. The heatmap exhibited the expression levels of the 6 AG-lncs in the high- and low-risk groups. The Kaplan-Meier survival analysis showed the OS, the PFI, and the DSS of BC patients between high- and low-risk groups in the training cohort (C–E) and testing cohort (F–H).

Figure 4 The development of nomogram and ROC analysis. Construction of a nomogram containing the risk score and other independent clinical factors to predict BC’s 1-, 3-, and 5-year OS of BC in training cohort (A) and testing cohorts (B). *p < 0.05, **p < 0.01; ***p < 0.001. The calibration plots of the nomograms in terms of the agreement between nomogram-predicted and observed 1-, 3- and 5-year and 5-year survival outcomes of BC in training cohort (C) and testing cohort (F). The ROC curves of the risk score at 1-, 2-, and 3-year in training cohort (D) and testing cohort (G). The AUC for the risk score and other clinical features (age, stage, TNM stage) based on the ROC curves in the training cohort (E) and testing cohort (H).

The Differential Function Analysis and PCA

To evaluate the variation in enriched pathways, GSEA was performed in gene expression in two risk subgroups. In the high-risk group, biological occurrences such as homophilic cell adhesion and microtubule cytoskeleton organization were overall enriched, while in the low-risk group, immune-related pathways and biological processes were commonly abundant, such as activation of antigen receptor-mediated signaling pathway, immune response, adaptive immune response, and B cell activation (Figure 5A). Of note, the high-risk group was also abundant in specific energy-related biosyntheses, such as unsaturated fatty acid biosynthesis, steroid hormone biosynthesis, and steroid biosynthesis, which might promote malignant tumor progression (Figure 5B). Meanwhile, in the high-risk group, the extracellular matrix (ECM) receptor interaction and cell cycle were also enriched. Additionally, the GSEA found that low-risk patients were inclined to suffer from autoimmune diseases, such as autoimmune thyroid disease, systemic lupus erythematosus, allograft rejection, primary immunodeficiency, and viral myocarditis (Figure 5B). These findings provided the possibility of further exploring the distinction of immunotherapeutic treatments between the two risk subgroups. Moreover, the PCA exhibited three distinctive distribution patterns in BC patients. Compared to all AGs (Figure 5C) and all AG-lncs (Figure 5D), BC patients displayed two distinct risk stratifications based on 6 AG-lncs (Figure 5E).

Figure 5 Functional enrichment analysis between the high- and low-risk groups by GSEA. (A) The GO analysis of the significantly enriched biological process between the high- and low-risk groups. (B) The KEGG analysis of the significantly enriched signaling pathways between the high- and low-risk groups. The 3D scatterplot of the PCA showed the distribution of BC patients based on all AGs (C), all AG-lncs (D), and the 6 AG-lncs (E).

Comparing the Risk Signature with Other Established Models

Moreover, the predictive value was assessed in the external validation cohort through ROC analysis to verify the reliability of the AG-lnc signature. Compared with the other prognostic signature in BC established by Zhao et al,¹⁵ Zhang et al,²⁶ Ping et al,²⁷ and Luo et al,²⁸ the AUC value at the 5-year OS of the risk score was the highest, illuminating the optimal ability of our risk signature in prediction of the longer survival outcomes in BC patients (Figure 6A). Besides, the concordance index of our risk model was also the highest, further supporting the good representation of our risk signature in BC prognosis (Figure 6B). The Kaplan-Meier analysis so showed the different OS of low- and high-risk BC patients based on our risk signature and other risk signatures (Figure 6C). The restricted mean survival curve showed all established signatures’ appropriate HR and p-values (Figure 6D).

Figure 6 Verification of the optimal aging-related lncRNA signature. (A) The ROC curves of the risk score and other 4 signatures at 1-, 3-, and 5 years. (B) The concordance index of the risk score and the other 4 signatures. (C) The Kaplan-Meier survival analysis for OS of BC patients between high- and low-risk groups based on the risk score and other 4 signatures. (D) The HR of the risk score and the other 4 signatures.

Analysis of Clinical Correlations and Genetic Alterations

The heatmap displayed the clinical characteristics of 6 AG-lncs (Figure 7A). Among those 6 AG-lncs, al136531.1, mapt-as1, al451085.2, tnfrsf14-as1, and linc01871 were higher expressed in low-risk patients, while otud6b-as1 was higher expressed in high-risk patients (Figure 7B). Meanwhile, the distinct distributions of the risk score among different clinical stages were analyzed via the Wilcoxon signed-rank test (Figure 7C–H). BC patients got higher risk scores at age > 65, death stage, late TNM stages, and more advanced clinical stage, suggesting an underlying prediction of the risk score in different clinical outcomes of BC patients. Interestingly, the frequency of CNVs and profiles of gene somatic mutations were considerably distinct between the two risk subgroups (Figure 7I and J). Patients in the low-risk group had the most genetic alterations of pik3ca (35%), while patients in the high-risk group possessed the most genetic alterations of tp53 (37%).

Figure 7 The correlation of the risk score with the clinical features of BC patients. (A) The heatmap visualizing the distribution of the 6 AG-lncs in two risk groups and different BC clinical outcomes (age, gender, TNM stage, and state). (B) The expression of 6 AG-lncs between high- and low-risk groups. ***p < 0.001. The correlation of the risk score with various clinical characteristics of BC, including age (C), state (alive or dead) (D), clinical stage (E), and TNM stage (F–H). The mutation frequency and classification of all genes in BC patients between low-risk group (I) and high-risk group (J).

Landscape of Tumor Immunity

The correlation between the immunity and risk score was explored by analyzing the various immune cell abundances among two risk groups in 7 datasets (Figure S3). With the risk score increasing, the heatmap visualized dynamic changes in immune cell infiltration and expression profiles. Additionally, the CIBERSORT algorithm analyzed 22 different immune cell types, uncovering that memory B cells, naive B cells, plasma cells, follicular helper T cells, CD4 memory-activated T cells, CD8+ T cells, regulatory T cells, activated NK cells, and resting dendritic cells were enriched in the low-risk group, while M0 and M2 macrophages were enriched in the high-risk group (Figure 8A). Moreover, the ssGSEA quantified the enrichment scores of 13 immune cell-related functions, discovering that in the low-risk group, CCR, inflammation-promoting, check-point, APC co-inhibition, and cytolytic activity were remarkably activated (Figure 8B). These results demonstrated a highly immune-stimulated environment among low-risk patients. Furthermore, the violin plots visualized the differences in the tumor purity, ESTIMATE score (immune score combined with stromal score), immune score, and stromal score among two risk subgroups (Figure 8C). The high-risk patients displayed higher tumor purity but lower ESTIMATE and immune scores compared with low-risk patients. The high-risk score generated a higher tumor purity, accompanied by altering immune stromal cells in the tumor microenvironment. Among those immune cells correlating with BC patients’ survival, follicular helper T cells, memory B cells, M0 macrophages, and M2 macrophages were likely the risk factors, while the plasma cells, naive B cells, and memory resting CD4+ T cells were likely the protective factors (Figure S4). In addition, the heatmap further exhibited the relationship between the interactive immune cells (Figure 8D). In Pearson’s correlation analysis, follicular helper T cells, naive B cells, and plasma cells showed a negative, while M0 and M2 macrophages were positively correlated with risk scores (p < 0.05, | r | > 1) (Figure 8E).

Figure 8 The differential immune state between the two risk subgroups. Comparison of the ssGSEA scores of 16 types of immune cells (A) and 13 immune-related pathways (B) between low-risk (blue box) and high-risk (red box) groups in the TCGA project. *p < 0.05, **p < 0.01; ***p < 0.001. (C) The tumor purity, ESTIMATE score, Stromal score, and Immune score between the low- and high-risk groups. (D) The correlation between immune cells. (E) The Pearson correlation analysis of the risk score with the immune cells.

Profiles of Distinct Immunotherapy Responsiveness Between Two Risk Groups

TIDE is commonly an effective tool for predicting immunotherapy efficacy.²⁹ A higher TIDE prediction score represents a higher immune escape possibility and a worse immune checkpoint inhibitor (ICI) treatment efficacy.³⁰ In the present study, the patients in the high-risk group got a lower TIDE score and a lower T-cell dysfunction score but a higher T-cell exclusion score, confirming that our high-risk patients might possess better efficacy of ICI therapy (Figure 9A). The occurrences of microsatellite instability showed no obvious distinction between the two risk groups. Additionally, two risk groups exhibited differential expression profiles of immune checkpoint-related molecules (Figure 9B). These results provided novel insights into the risk score’s potential application in assisting in predicting clinical immunotherapy responsiveness. Intriguingly, m6A-related gene expression showed different distributions between the two risk groups, proposing a profound exploration of a possible connection between epigenetic patterns and immune landscape or immunotherapy effects (Figure 9C). Besides, the results of the IC50 analysis of chemotherapeutic drugs exhibited that high-risk BC patients had a better response to doxorubicin, etoposide, and rapamycin (p < 0.01) (Figure 9D–F).

Figure 9 The differential response to the immune therapy between the high- and low-risk groups. (A) The TIDE, microsatellite instability, and T-cell exclusion and dysfunction score in two risk subgroups. The scores between the two risk subgroups were compared through the Wilcoxon test. ns, no significance; ***p < 0.001. (B) The expression levels of immune checkpoint molecules between the high- and low-risk groups in BC. *p < 0.05, **p < 0.01; ***p < 0.001. (C) The expression levels of m6A-related genes between the high- and low-risk groups in BC. **p < 0.01; ***p < 0.001. The IC50 of chemotherapeutic drugs such as doxorubicin (D), etoposide (E), and rapamycin (F) was higher in the high-risk group.

Validating the Expression of 6 AG-Lncs in BC Tissues

20 BC samples with different clinical stages were selected as a new clinical validation cohort. Based on the expression of AG-lncs from qRT-PCR analysis, the risk score was re-calculated as: the risk score = al136531.1_expression × (−0.3430) + mapt-as1_expression × (0.3604) + al451085.2_expression × (0.2184) + otud6b-as1_expression × (0.0467) + tnfrsf14-as1_expression × (0.2501) + linc01871_expression × (0.2600). The coefficient was derived from the multi-Cox regression in the training set of BC patients from the TCGA. Based on the median re-calculated risk score, 20 BC patients were classified into low- and high-risk subgroups. Among them, al136531.1, mapt-as1, al451085.2, tnfrsf14-as1, and linc01871 expressed higher in low-risk patients, while otud6b-as1 expressed higher in high-risk patients (Figure 10A–F). This validation was in accordance with the sequencing outcome from the TCGA (Figure 7B), proposing a further validation or exploration of the biological function of these 6 AG-lncs in BC progression or clinical prediction.

Figure 10 Validation of the predictive ability of the risk score in an external cohort. The expression of 6 AG-lncs in BC patients at different clinical stages: (A) al136531.1, (B) mapt-as1, (C) al451085.2, (D) otud6b-as1, (E) tnfrsf14-as1, and (F) linc01871. *p < 0.05, **p < 0.01; ** p < 0.001.

Discussion

BC encompasses several biological entities with disparate molecular subtypes and outcomes that require priorities for identifying the prognostic indications. Constructing predictive models, such as a stratified machine learning system integrated with multi-omics data, can confer predictive efficacy for the diagnosis and prognosis of breast cancer. For instance, Tabl et al constructed a hierarchical machine learning system on the basis of gene expression data, clinical information, and treatment type from 347 BC patients, which could predict the five-year survival rate of BC patients.³¹ The Nottingham Prognostic Index (NPI) is an indicator for prognosis after BC surgery, which is based on lymph node status, tumor size, and histological grading.³² Zhou et al built a residual neural network classification model with gene expression, CNV, and mRNA from 1885 female BC patients by high-dimensional embedding to extract prognostic biomarkers to predict NPI categories with higher accuracy.³³ Similarly, Yan et al developed and validated a diagnostic nomogram according to histopathological outcomes, neoadjuvant chemotherapy cycle count, and tumor size, identifying BC patients with a high probability of pathological complete response.³⁴ Moreover, systematic analysis combining multiple databases could contribute to the screening of prognostic biomarkers in BC.³⁵ Apart from traditional tissue biopsy, non-invasive liquid biopsy is also a promising tool for BC prognosis prediction.³⁶

The current study successfully established an effective signature of 6 AG-lncs in predicting BC prognosis. Firstly, after selecting the expression of DEAGs between BC and normal from the TCGA project, 6 AG-lncs were used to establish a predictive signature in BC prognosis. Secondly, through the multivariate and univariate Cox regression, ROC curve analysis, and survival analysis in the training cohort and testing cohort, the AG-lnc model was apparently correlated with BC survival and had an independent predictive value compared with other clinics’ pathological features. Importantly, the risk score could effectively predict the long-term survival rate in BC patients as the augmented AUC value with the survival time increasing. Thirdly, it further found that high-risk patients exhibited more induced immune infiltration, and might gain more benefits from ICI treatment. Besides, m6A-related genes also exhibited differential expression patterns between the two risk groups, revealing the possible regulatory mechanisms of epigenetics. Notably, the external clinical cohort validated the differential expression of all 6 AG-lncs in BC tissues with different stages. Our results demonstrated that the 6 AG-lnc signature might complement BC patients’ clinical and pathologic stratification, contributing to BC patients’ diagnosis and rational therapy.

The natural aging process is a gradual degeneration and loss of homeostasis resulting from transforming constitutive substances and tissue structures.³⁷ Aging plays a multifaceted role in tumor biology, mediating tumor cell senescence or promoting growth and metastasis.³⁸ For example, gliomas are heterogeneous brain malignancies with poor prognoses.³⁹ Gnanavel et al constructed a miRNA signature that was related to the aging-related senescence in glioblastoma, intensively suggesting that senescence-linked molecular mechanisms were involved.⁴⁰ Coincidentally, Luo et al also chose 4 AGs (sstr3, pon1, tert, and lep), but not miRNAs, from TCGA and Chinese Glioma Genome Atlas to establish a risk model.⁴¹ It demonstrates the value of this model in predicting clinical outcomes in glioma patients. In our study, we successfully established a risk signature of BC prediction based on 6 AG-lncs, including al136531.1, mapt-as1, al451085.2, otud6b-as1, tnfrsf14-as1, and linc01871. Under this model, the risk score successfully stratified BC patients into two different survival subgroups, where patients in the high-risk subgroup had a shorter survival time, PFI, and DSS than low-risk patients. In marked contrast to previous studies, here we used AG-lncs as combinatorial features instead of AGs or miRNAs, also achieving good predictive results. The potential explanation may lie in the high tumor specificity of lncRNA and aging-related lncRNA expression patterns.

Aging-related risk signature is a reliable predictor in evaluating the OS, the severity as well as the features of immune profiles in cancer.⁴² Given the complex and dynamic nature of immune cell distribution in the tumor microenvironment, the immune response to therapy is still unpredictable Increasing evidence agrees with the hypothesis that immune cell infiltration patterns are crucial to BC progression.⁴³ For example, in head and neck squamous cell carcinoma (HNSCC), Yang et al analyzed the potential link among AGs, prognosis, inflammation, and tumor immunity, by constructing a risk signatue of 7 AGs.⁴⁴ This robust model indicated an independent prognostic index of the risk score, displaying the immunosuppressive and chronic inflammatory status in HNSCC patients. It provided pivotal clues for individualized immunotherapy in cancer patients. Zhang et al designed and verified a reliable signature based on 24 AGs to forecast gastric cancer patients’ prognosis in an independent Gene Expression Omnibus dataset.⁴⁵ The immune infiltration analysis verified that the amount of follicular helper T cell cells was much higher in low-risk patients, while the amount of monocytes had a reverse trend.

Similarly, our results verified that according to the risk score, there were remarkably higher expression patterns in various immune function scores including CCR, inflammatory promotion, checkpoint, APC co-inhibition, and cytolytic activity in the low-risk subgroup compared with the high-risk subgroup. Furthermore, the tumor purity was enhanced with the risk score increasing, as well as resulting in stromal immune cell changes, showing a consistent phenomenon that AG-lncs were strongly associated with the immune infiltration in BC. Deciphering these AGs and AG-lncs will benefit to design the effective immunotherapy strategies.

In exploring our risk signature, low-risk patients related to a highly immune-stimulated environment, while high-risk patients displayed lower immune scores with more abundant M2 and M0 macrophages. Moreover, al136531.1, mapt-as1, al451085.2, tnfrsf14-as1, and linc01871 expressed higher in low-risk patients, while otud6b-as1 expressed higher in high-risk patients. LncRNAs play a potentially important role in modulating immune cell function, reshaping tumor immune environment, and influencing BC progression. Interestingly, many other studies drew similar conclusions as we did. In HER2+ BC, linc00624 inhibited the presentation of MHC-I antigen, the infiltration of CD8+T cells, and the activation of the IFN pathway, leading to HER2-targeted therapy resistance.⁴⁶ Liu et al found that the lncRNA irena could facilitate the macrophage conversion from tumor-suppressing phenotype to tumor-promoting phenotype by triggering NF-κB signaling, which enhanced inflammatory cytokine production.⁴⁷ Furthermore, lncRNA xist competed with miR-101 to suppress macrophage M2 polarization, inhibiting BC cell migration and proliferation.⁴⁸ Diverse lncRNAs acted as suppressors or promotors of tumor immune responses in BC, regulating immune cell functions. Therefore, dysregulated lncRNAs exert a critical function in BC immune microenvironment and are promising biomarkers to indicate immunotherapy prognosis. In vitro and in vivo models will be conducted for the following in-depth investigations, such as the BC three-dimensional spheroid model for tumorigenesis and drug resistance assessment.⁴⁹

Nevertheless, some issues remain to be resolved in this study. BC is a frequently diagnosed cancer type with high heterogeneity and plasticity. BC possesses complicated molecular typing and tumor subtypes, including HER+ BC and triple-negative breast cancer (TNBC), which are accompanied by great differences in prognosis. Therefore, whether the risk model composed of AG-lncs could predict prognosis in BC subtypes with different molecular-driven characteristics. Specifically, future studies should validate the applicability of our AG-lncs-based risk model in different subtypes of hormone receptor-positive, HER2-positive, and TNBC. Such targeted analyses may reveal the specific role of AG-lncs in different subtypes and improve risk stratification and prognostic accuracy. In addition, the establishment of prospective cohort studies with comprehensive molecular characterization could provide the data needed to assess the prognostic value of AG-lncs across subtypes and to explore potential mechanistic pathways of how AG-lncs influence tumor biology in different subtypes. Combining AG-lncs expression profiles with detailed clinical data, including treatment response and outcomes, will enhance our understanding of the use of these biomarkers in clinical decision-making and personalized medicine strategies. In addition, functional studies are needed to elucidate the biological roles of identified AG-lncs in breast cancer, which may involve in vitro and in vivo models to examine the effects of AG-lncs on cancer cell proliferation, migration, immune escape, and response to therapy. Given the increasing role of immunotherapy in the treatment of breast cancer, especially in subtypes such as TNBC, it will be valuable to study the association between AG-lncs and the immune microenvironment, which may lead to the identification of new markers predictive of immunotherapeutic response. By addressing these questions, it will contribute to determining the prognostic evaluation of AG-lncs-associated models in different subtypes of breast cancer. It should also be mentioned that the distinct sensitivity of immune checkpoint inhibitors in different risk groups of BC was of value to conduct sample verification and discussion in vivo and in vitro experiments. The greatest strength of this study is that it utilizes a gene expression model that provides some reference for the clinical treatment and prognosis of BC.

However, it is undeniable that this study cannot replace the traditional gold standard approach. Emerging technologies and advanced bioinformatic tools, including single-cell sequencing and multi-omics approaches, could be utilized to enhance the precision and scope of our risk model. These tools will enable us to better understand the heterogeneity of the 6 AG-lncs in different BC subtypes and to explore their interactions with the tumor immune microenvironment. Specifically, the mechanism of AG-lncs in regulating immune cell infiltration and immune escape mechanisms is still a mystery, which is critical for the development of novel immunotherapeutic strategies, including personalized cancer vaccines or immune checkpoint inhibitors. Furthermore, we intend to integrate our findings with other established biomarkers to develop a more holistic diagnostic and prognostic model for BC in multiple dimensions, such as ALDH1,⁵⁰ a marker of breast cancer stem cells. In addition, it is important to investigate AG-lncs as potential therapeutic carriers or targets. Using nanotechnology and targeted delivery systems, it is feasible to construct novel therapeutic carriers to deliver small-molecule drugs, siRNAs, or CRISPR gene editing tools directly into tumor cells or specific microenvironmental cells. Such strategies can improve the precision of treatment and minimize the toxic side effects on normal tissues, providing safer and more effective therapy options for BC patients. Specifically, small molecule drugs, known as biological missiles, are a targeted therapeutic form with high specificity and few adverse effects.⁵¹ By constructing the “drug-ingredient-target-disease” network, network pharmacological analysis could be a valid path to identify and filter the active ingredients against BC.⁵² In summary, we are dedicated to advancing our work in a direction that facilitates personalized medicine and improves BC patient prognosis.

Conclusion

Collectively, an effective predictive risk model was constructed in BC prognosis based on 6 AG-lncs, including al136531.1, mapt-as1, al451085.2, otud6b-as1, tnfrsf14-as1, and linc01871. More intriguingly, high-risk BC patients had presumptive worse clinical outcomes, whereas they might be in an immune-desert state and might respond more positively to immunotherapy. This well-established prognostic signature under AG-lnc patterns might be conducive to BC prognosis recognition and rational individual therapy.

Abbreviations

AG, aging-related gene; AG-lncs, aging-related lncRNAs; AUC, area under curve; BC, breast cancer; CBC, contralateral breast cancer; cDNA, complementary DNA; DEAGs, differentially expressed AGs; DSS, disease-specific survival; ECM, extracellular matrix; GO, Gene Ontology; GSEA, gene set enrichment analysis; HNSCC, head and neck squamous cell carcinoma; HAGR, human aging genome resource; ICI, immune checkpoint inhibitor; KEGG, Kyoto Encyclopedia of Genes and Genomes; lncRNAs, long non-coding RNAs; ncRNAs, non-coding RNAs; OS, overall survival; PCA, Principle Component Analysis; PFI, progression-free interval; PPI, protein-protein interaction; ROC, receiver operating characteristic; ssGSEA, single-sample gene set enrichment analysis; TCGA, The Cancer Genome Atlas; TNBC, triple-negative breast cancer; TIDE, tumor immune dysfunction and exclusion.

Ethics Approval and Informed Consent

The study was performed after agreement from the ethics committee of Tongji Hospital of Huazhong University of Science and Technology the patients’ informed consent. The clinical samples in the study comply with the Declaration of Helsinki.

Consent for Publication

All authors have provided their consent for publication.

Author Contributions

All authors made a significant contribution to the work reported, whether that is in the conception, study design, execution, acquisition of data, analysis, and interpretation, or in all these areas; have drafted or written, or substantially revised or critically reviewed the article; have agreed on the journal to which the article will be submitted; reviewed and agreed on all versions of the article before submission, during revision, the final version accepted for publication, and any significant changes introduced at the proofing stage; and agree to take responsibility and be accountable for the contents of the article. Zi-Hui Yang and Hong Zeng are co-correspondence authors for this study.

Disclosure

The authors declare that they have no conflicts of interest to report regarding the present study.

Data Sharing Statement

All the datasets displayed in this study can be obtained in the online database. Further questions can be directed to the corresponding author.

Additional information

Funding

The authors received no specific funding for this study.