Reprogramming of metabolism has received strongly increasing attention within the last decade [37], and several studies have shown that altered mitochondrial metabolism is considered to be one of the major emerging mechanisms of therapeutic resistance [38–41]. We constructed a risk signature composed of 8 MRGs (MRPL13, BCL2A1, ME3, SLC1A1, TP53AIP1, MTHFD2, ALDH2, and ACSL1) and fitted a nomogram model based on this signature according to corresponding clinical features. A brief description of the related functions of the 8 MRGs is shown in Table 2. Functionally, all 8 MRGs have been reported to be associated with cancer, and most of them are related to immunity and cell death, which arouses interest for further investigation.
We also tested the accuracy of the risk signature by ROC and AUC, in which the total AUC was 0.647. Furthermore, we calculated the AUC of the breast cancer patients for 1, 3, 5, and 10 years, and its trend was gradually strengthened (from 0.58 to 0.76), which means that this model also has a certain accuracy in the long-term prognosis of BC (Fig. 3b-c). Combining the riskScore with clinical features for Cox regression, it was found that the HR of the riskScore reached 3.496 with statistical significance in univariate Cox regression, and multivariant Cox regression also proved that the riskScore was also an independent risk factor for prognosis (Fig. 3d-e).
DCA is a method for evaluating nomograms that can meet the practical needs of clinical decision-making [42], and a calibration curve is another method to evaluate the consistency between the predicted and true values. The above methods proved the effectiveness of the nomogram model based on the risk signature. In addition, the prognosis of high-risk patients was significantly worse, and a higher riskScore was positively correlated with more aggressive clinical features (including higher PT, PN, AJCC stage, and more aggressive subtype). The corresponding results were also verified in the METABRIC cohort (Supplementary Fig. 1–2).
To be specific, the prognosis with high ME3 and TP53AIP1 was significantly better, while the prognosis with high MRPL13 was worse, and combined with the gene expression, it suggested that TP53AIP1 and ME3 may act as tumor suppressors and MRPL13 as an oncogene in BC. Similar to what we inferred, one study showed that TP53AIP1 has been proven to be a new tumor suppressor gene for BC and may become an effective target gene for therapy [43]. Another study reported that high MRPL13 is a poor prognostic factor for BC, and it can be used as a molecular marker for prognosis judgment and as a potential therapeutic target [24, 44]. It is worth noting that MRGs are associated with the immune response and cell death [45, 46], which also proves the reliability and potential to be used as biomarkers of our risk signature.
To explore the underlying reasons, we analyzed DEGs between risk groups and performed functional enrichment analyses. The results showed that in GO, most BP terms were related to the immune response, and KEGG analysis showed that Th17-cell differentiation, the NF-kappa B signaling pathway and most terms were closely related to the immune response and oncogenesis. GSEA also demonstrated similar results; in the low-risk group, the top 10 terms were all related to the immune response. These findings indicate that a low riskScore may be accompanied by an abundant immune response.
The ESTIMATE algorithm can effectively evaluate the proportion of tumor, infiltrating immune, and stromal cells in the tumor microenvironment (TME). Similar to our findings, the low-risk patients exhibited higher scores of ESTIMATE, immune and stromal, while lower tumor purity, that is, presented a higher degree of immune infiltration and showed a “warmer” condition.
Furthermore, CIBERSORT could evaluate the infiltration degree of 22 types of immune cells [36]. Among them, the riskScore was significantly negatively correlated with CD8 T cells, which are critical for tumor destruction [47], and the correlation coefficient was also the highest (R = -0.26, P < 2.2e-10). Moreover, in M2 macrophages, which tend to promote angiogenesis and neovascularization [48], the riskScore showed a significant positive correlation (R = 0.19, P = 2e-10). Notably, M2 macrophages can cause stromal activation and remodeling [49], which are endowed with a repertoire of tumor-promoting capabilities involving immunosuppression and result in poor prognosis of BC patients [50]. Meanwhile, the riskScore showed a negative correlation with M1 macrophages, which could enhance antitumor inflammatory reactions and act as major players in proinflammatory responses [51].
Meanwhile, based on multiple control modalities of intracellular and extracellular networks, TCGA cohort was divided into 5 subtypes [22], in which C1 was higher in high-risk (38% vs. 26%), and C3 was higher in low-risk (16% vs 7%). Among these subtypes, C1 (Wound Healing immune) has high expression of angiogenic genes and Th2 cells, which have high tumor cell proliferation and high intratumoral heterogeneity with a less favorable outcome. C3 (inflammatory), defined by elevated Th17 and Th1 genes, is associated with low to moderate tumor cell proliferation and minimal intratumoral heterogeneity, which has the best prognosis compared to other subtypes. After ranking the riskScore from low to high, the immune infiltration of individual cases was shown by ssGSEA. Consistent with the above, a low riskScore was accompanied by increased immune cell infiltration, making the immune infiltration of low-risk patients look “warmer”. Therefore, it is reasonable to speculate that the risk signature can divide patients with different immune responses, and this heterogeneous status of tumor immunity may account for the difference in prognosis.
Moreover, we explored the mutational landscape between the risk groups. Among them, TP53 was the major mutation in high-risk patients (37%), and the incidence of PIK3CA mutation was dominant in low-risk patients (42%). TP53 mutation was associated with frequencies of mutations in TNBC and HER2-positive subtypes of 80% and 70%, respectively, and in luminal A and B types of 10% and 30%, respectively [7, 52–54]. PIK3CA mutation has been widely recognized as a genomic marker of BC [55], and PIK3CA mutation rates were lower in TNBC (16%) than in HR+/HER2-(42%) and HER2+ (31%) BC [56]. These findings partly explained the variation in riskScore among molecular subtypes. In addition, the same gene mutation status had a significantly worse prognosis in high-risk patients but not in low-risk patients. This phenomenon may be due to the mutually exclusive nature of PIK3CA mutations and TP53 mutations in the low-risk group, so in the low-risk group dominated by PIK3CA-mut, TP53 tends to be the wild type, which has a positive effect on the prognosis of BC [57, 58]. In contrast, in the high-risk group, TP53 and PIK3CA mutations were not mutually exclusive, and they could be co-occurring, reflecting the difference in prognosis. However, the higher TMB in high-risk patients may be due to the more TP53 mutations and the higher proportion of triple-negative breast cancer [59].
After WGCNA, it was found that the turquoise module gene trend was the most consistent with our speculation. GO analysis of the turquoise module genes also revealed that immune regulation may account for the differences. In terms of treatment options, the high-risk cohort had a worse prognosis in three treatments regardless of TCGA or METABRIC cohort, which also reflected the reliability of our risk signature. The results from drug sensitivity analysis showed that low-risk patients were more sensitive to antitumor drugs, such as docetaxel, epirubicin and paclitaxel, which are the regular chemotherapy regimens for BC.
In this work, we established an 8-MRG-based risk signature for BC prediction through joint analysis of gene expression datasets from TCGA and METABRIC, which also effectively stratified BC patients. The signature was constructed based on transcriptome data, and its feasible clinical translation ability was demonstrated by clinical validation methods. Moreover, this method might also be suitable to explore the prognostic effects of mitochondria-related genes in other malignant tumors. However, there are several limitations in our study. First, the risk signature was established based on prognosis, so its application is mainly reflected in the survival of BC patients and screening of high-risk populations but has little effect on early diagnosis and screening. Second, our analysis only used data from public databases, and more real-world data are needed to further confirm our findings. Third, although several basic experimental studies have validated the function of certain genes in our risk signature in breast cancer (such as ACSL1 [60, 61], MTFHD2 [62], MRPL13 [24, 44, 45], TP53AIP1 [43] and BCL2A1 [63]), the function of other genes and the underlying mechanism of these 8 MRGs need to be further investigated, and we will conduct corresponding experimental studies in the future.