Introduction

Breast cancer (BC) is the most commonly diagnosed female malignancy, with an increasing rate of about 0.5% annually1. As the second leading cause of death, many clinicopathological biomarkers have been utilized to predict, and therefore improve its survival2,3,4. However, BC is a highly heterogeneous cancer5, and conventional approaches such as TNM staging and molecular subtypes are insufficient. Fortunately, the recent development of molecular signatures extensively promoted the survival of BC patients and the development of targeted therapy6,7,8. It leads to the great importance of gene expression profiling in prognosis prediction and high-risk BC patient selection.

Genome instability fosters aberrant hallmark functions of cancer, one of which is mitochondrial energy metabolism and apoptosis9,10. Mitochondrial ribosomal proteins (MRPs), encoded by the nuclear MRP gene family11, are indispensable components in mitochondrial translation. In previous studies, more than 40 MRPs were reported to be overexpressed in BC, acting as promotors of BC cellular viability12,13. Recently, X Li et al. reported MRPL52 upregulation in BC and its oncogenic role in apoptotic resistance and metastatic promotion14. Some MRPs serve as tumor suppressors. For example, MRPL41 is downregulated in BC15. Its upregulation was associated with positive estrogen receptor status and could be induced by estradiol and trichostatin A16.

MRP family members are important regulators in BC development, but the prognostic values for BC patients have not been elucidated. In this study, we conducted a comprehensive bioinformatics analysis of the MRP family. The prognostic effect was investigated, and an MRP-based model was further established to predict the BC survival. We aimed to help clinicians to predict the risk of BC patients and provide researchers with new therapeutic candidates.

Material and methods

Differential expression analysis

Gene Expression Profiling Interactive Analysis (GEPIA2), a web server with RNA sequencing expression data17, was used to find out the differentially expressed genes (DEGs) in BC compared with normal tissue. Based on the analysis of variance (ANOVA), the Log2(Fold Change) cutoff was 1 and the P value cutoff was 0.001. 1085 BC samples were analyzed and matched with 291 samples from TCGA normal and GTEx data. The DEGs were screened for subsequent molecular analysis. Additionally, the differential expression of MRP genes in various primary tumors, as well as BC cell lines, were validated in the Cancer Cell Line Encyclopedia (CCLE) database18. Due to the high heterogeneity of BC, we further explored the DEGs in different BC subtypes relative to normal tissue using GEPIA2, with the same criteria and matched normal data mentioned above. Differential protein expression of BC was also explored based on the Clinical Proteomic Tumor Analysis Consortium (CPTAC) Confirmatory/Discovery dataset in UALCAN19,20.

Genomic alteration and methylation analysis

cBioPortal is a platform containing multidimensional genomics datasets21, and was utilized to visualize the genomic alteration of DEGs and to evaluate the mRNA expression correlation with DNA copy number variation. Copy number alteration was determined using GISTIC 2.0, in which a value of 2 is equivalent to amplification. Besides genetic mutation and copy-number alteration, DNA methylation also plays a vital role in cancer progression9. Therefore, the Spearman correlations between MRP methylation and RNA expression (RNA Seq V2 RSEM, log2(value + 1)) were conducted. Moreover, we explored the prognostic values of genetic alteration and methylation levels in BC by cBioportal and MethSurv, respectively. MethSurv is an online tool providing survival analysis based on DNA methylation data from the TCGA database22. We employed a “single CpG” analysis module, and all the available genomic regions and individual CpG sites were evaluated. BC patients were dichotomized by median methylation levels, and the log-likelihood ratio (LR) test P value < 0.05 was considered significantly different in prognosis.

Functional enrichment analysis and protein interaction visualization

Gene Ontology (GO) represents the biological function of candidate genes23, and the Database for Annotation, Visualization, and Integrated Discovery (DAVID) was applied to conduct GO enrichment analysis24,25. Since proteins form protein complexes and play biological roles synergistically, protein–protein interaction (PPI) network construction would be essential to understanding MRPs. String, a database integrating all known physical and functional protein interactions26, was employed to explore the MRP family members. Cytoscape was adopted to integrate and modify the network27. The above-mentioned Log2 (Fold Change) was used to represent the node size, indicating the significance of expression difference.

Data acquisition

1056 Primary BC samples with complete transcriptome data and clinical traits in the TCGA database were acquired via the UCSC Xena platform28, after excluding cases with incomplete overall survival (OS) information, or unknown age and TNM staging. The TNM stage was decided by the American Joint Committee on Cancer Tumor Stage Code. The status of hormone receptor (HR) and human epidermal growth factor receptor 2 (HER2) were assessed by immunohistochemistry (IHC) or fluorescence in situ hybridization (FISH). The mRNA profiles were measured using the Illumina HiSeq 2000 RNA Sequencing platform and shown as log2(value + 1) transformed RSEM normalized count. 77 MRP family members were included in the gene expression profiles, except MRPL57 which was not covered in the TCGA dataset. Based on the median level of each MRP, the samples were classified into high- and low- expression groups.

Survival analysis and prognostic model establishment

We performed univariate Cox proportional hazards regression to determine the prognostic factors for OS among all 77 MRP family members. The significant prognostic effects of MRPs were plotted using Kaplan–Meier curves. The least absolute shrinkage and selection operator (LASSO) regression was performed to exclude the collinearity between variables and improve the accuracy of the model. Stepwise multivariate Cox regressions were carried out to select the optimal combination of covariates. Based on the combined results, we developed an MRP-based prognostic signature, evaluated by calibration curves and the area under the curve (AUC) of the time-dependent receiver operating characteristic (ROC) curves.

Statistical analysis

Statistical analysis was performed using R 4.0.2 software. All statistical tests were two-sided. A P value less than 0.05 was used as a cutoff for statistical significance.

Ethical approval

There was no need for ethical approval as all data analyzed were obtained from public databases.

Results

Identification of differentially expressed MRP family members

In the MRP family, 12 members were specifically upregulated in BC (P < 0.001) and thus identified as differentially expressed genes (DEGs), including MRPL3, MRPL13, MRPL14, MRPL17, MRPL24, MRPL42, MRPL47, MRPS23, DAP3, MRPS30, MRPS34, MRPS35 (Fig. 1). The most overexpressed gene was MRPL14, with a Log2 (Fold Change) of 1.336, followed by MRPS34, with a Log2 (Fold Change) of 1.333. Among various types of cancer, the DEGs expressions of BC were comparable to other primary tumors (Fig. S1). But the median expressions of MRPL3, MRPL13, MRPL14, MRPL24, MRPL47, DAP3, and MRPS34 were higher than the average. The mRNA levels of DEGs across different BC cell lines were then illustrated in Fig. S2. It was noteworthy that all DEGs’ expressions were the lowest on BC cell line SUM102PT. Low expressions of MRPL3, MRPL42, MRPL47, and low expressions of MRPL13, MRPL14, MRPS23 were detected in HDQP1 and SUM185PE cells, respectively. Moreover, we discovered that MRP genes were upregulated in a subtype-specific manner (Table 1). Half of the 12 DEGs (MRPL13, MRPL14, MRPL17, MRPL42, DAP3, and MRPS35) were highly upregulated in all the subtypes. Apart from the above-mentioned MRP genes, high expression of MRPS28 was specifically observed in the Luminal B subtype, MRPL27 in the HER2-enriched subtype, and MRPL9, MRPL37, MRPL41, MRPL48, MRPS12 in the basal-like subtype. These results suggested the oncogenic effects of MRP family genes on BC progression.

Figure 1
figure 1

Differential expression of MRP genes between breast cancer and normal tissues (GEPIA2). Red and grey box plots shown the expression of breast tumor and normal tissues respectively. The asterisk (*) indicated significant differences (P < 0.001) between the two groups.

Full size image
Table 1 Differential MRP gene expression in different breast cancer subtypes.
Full size table

Genetic and epigenetic features, and the correlation with MRP expression

Using cBioPortal, genomic alteration analyses of DEGs in the MRP family were conducted. Overall, aberrant alterations occurred in 42% of the samples. As shown in Fig. 2A, MRPL13 (19%), MRPL24 (11%), and DAP3 (11%) were altered most frequently. Genetic amplification was the most common variation. A strong correlation was revealed between gene amplification and higher levels of the corresponding mRNA for most of the DEGs, such as MRPL13, MRPL14, and MRPS35 (Fig. 2B). As for epigenetic modification, all the DEGs’ expressions were negatively correlated with the methylation level of the MRP gene family in BC (Fig. 2C), except MRPL13 and MRPL47, whose methylation data was not included in cBioPortal. In addition, the genetic alterations of DEGs were not significantly associated with prognosis (Fig. S3), whereas the methylation levels of DEGs were remarkably correlated with OS (Fig. 2D). Hypermethylation of MRPL24 promoter, MRPL42 promotor, and CpG cg16002248 annotated to the gene DAP3 resulted in better OS in BC patients. By contrast, patients with high MRPS23 and MRPS35 methylation had worse OS.

Figure 2
figure 2

Gene alteration and methylation of MRP genes in breast cancer. (A) Summary of alteration patterns. The proportion of samples with genetic alteration was expressed as a percentage next to the gene, and genetic alteration types were indicated in bars with different colors. (B) The impact of gene amplification on the transcriptional level. Grey and red circles represented the mRNA expression in samples with diploid genotype and gene amplification. (C) Correlation between DNA methylation and mRNA expression in breast cancer. Each sample was indicated in a blue dot, and the fitted line was in red. The Spearman correlation coefficient and P value examined the association. (D) Kaplan–Meier plots showing overall survival in BC patients with hypermethylation (red curves) and hypomethylation (blue curves) for cg22493673-MRPL24, cg15127806-MRPL42, cg18503387-MRPS23, cg16002248-DAP3, and cg08925658-MRPS35. Log-likelihood ratio (LR) test P values were used to determine the prognostic significance and the hazard ratio (HR) was a relative prognostic measure of BC patients with hypermethylation compared with those with hypomethylation.

Full size image

MRP protein expression and functional enrichment analysis

At translational levels, among the 12 DEGs, MRPL3, MRPL13, MRPL17, MRPL47, MRPS23, MRPS30, DAP3, and MRPS35 exhibit significant differences in proteomic expression based on the CPTAC dataset. As demonstrated in Fig. 3A, the levels of these 8 proteins in BC were prominently higher than in normal tissues. To probe the biological function of MRPs, GO enrichment analysis of DEGs was performed. It was recognized that they were involved in mitochondrial translational termination, elongation, and translation, and were mainly located on the mitochondrial inner membrane, ribosome, and mitochondrial large ribosomal subunit (Fig. 3B). As for molecular function, they primarily played roles in the structural constituents of ribosomes and poly(A) RNA binding. The interplay among the DEGs was revealed using the PPI network (Fig. 3C). Remarkably, all of them interacted with each other, implying the MRP family functioned in its entirety.

Figure 3
figure 3

Proteomic and functional analysis of MRPs in breast cancer. (A) Differential proteomic levels of MRP genes between BC and normal tissues (CPTAC). (B) Gene Ontology enrichment analysis of MRP genes. (C) Protein–protein interaction network of MRP genes. The node size was proportional to Log2 (Fold Change).

Full size image

The prognostic value of MRP expression in breast cancer

Regarding BC prognosis, all 77 MRP genes were incorporated in univariate Cox regression (Table S1). As summarized in Table 2, the mRNA expression of 7 MRPs were revealed to be associated with OS in BC. Age (HR, 1.034; 95% CI, 1.021–1.047; P < 0.001), stage (III vs I: HR, 3.026; 95% CI, 1.706–5.368; P < 0.001; IV vs I: HR, 13.042; 95% CI, 6.430–26.453; P < 0.001) and IHC-based subtypes (HR-/HER2- vs HR + /HER2-,1.811; 95% CI, 1.125–2.914; P = 0.015) were also prognostic factors for BC. In Fig. 4, Kaplan–Meier survival curves were plotted to assess the effects of differentially expressed MRPs on OS. It was manifested that higher transcription levels of MRPL1 (HR, 1.518; 95% CI, 1.082–2.128; P = 0.016), MRPL13 (HR, 1.514; 95% CI, 1.079–2.123; P = 0.016), MRPS6 (HR, 1.440; 95% CI, 1.030–2.014; P = 0.033), MRPS18C (HR, 1.405; 95% CI, 1.005–1.963; P = 0.047), and MRPS35 (HR, 1.459; 95% CI, 1.044–2.039; P = 0.027) were correlated with deterioration of prognosis. Conversely, patients with downregulation of MRPL16 (HR, 0.630; 95% CI, 0.449–0.884; P = 0.007), and MRPL40 (HR, 0.643; 95% CI, 0.458–0.903; P = 0.011) had significantly shorter survival.

Table 2 Univariate Cox proportional hazards regression analysis in the TCGA patients, and multivariate analysis of the MRP-based nomogram.
Full size table
Figure 4
figure 4

Kaplan–Meier survival curves based on differential MRP gene expression levels in BC patients. Low and high expression groups were represented in green and red respectively. Black dashed lines shown the median survival time.

Full size image

Development of an MRP-based prognostic model

The prognostic effects of 7 MRPs were displayed in a forest plot (Fig. 5A) using HR (Hazard Ratio) and 95% confidence intervals (CI). Given the striking effect of MRPs on the BC survival outcome, an MRP-based prognostic model was developed to predict the OS. Firstly, the 7 MRPs were subjected to LASSO regression. The model fit the best when the penalty coefficient was 7 (Fig. 5B, C). Therefore, all MRPs were selected for stepwise multivariate Cox regression analyses. After excluding the variates without significant statistical difference and comparing the concordance index, 4 MRPs (MRPL16, MRPL40, MRPS18C, and MRPS35), age, and stage were integrated into the nomogram. Table 2 displayed the multivariate Cox regression analysis of the markers included in the nomogram. Patients with higher transcription levels of MRPS18C (HR, 1.412; 95% CI, 0.991–2.012; P = 0.056) and MRPS35 (HR, 1.48; 95% CI, 1.047–2.092; P = 0.027) indicated significantly shorter OS, while upregulated MRPL16 (HR, 0.663; 95% CI, 0.456–0.963; P = 0.031) and MRPL40 (HR, 0.567; 95% CI, 0.391–0.822; P = 0.003) were identified as independent protective factors. The nomogram was demonstrated in Fig. 6A. According to the nomogram scoring system, each variable was given a nomogram score on the point scale. They were summed up and a total nomogram score was obtained for each patient, acquiring a predicted probability of OS. Based on the median value of total scores, TCGA patients were stratified into high-risk and low-risk groups.

Figure 5
figure 5

Construction of an MRP-based nomogram. (A) Forest plot of univariate Cox regression analysis of MRP genes. (B) Least absolute shrinkage and selection operator (LASSO) coefficient profiles of the 7 MRPs determined by the optimal lambda. (C) Selection result of the optimal lambda in the model. The partial likelihood deviance (binomial deviance) curve was drawn versus log(λ), and vertical dotted lines were made at the optimal values.

Full size image
Figure 6
figure 6

MRP-based prognostic model and its evaluation. (A) A nomogram to predict the overall survival in BC patients. (B) Calibration curves of the nomogram to predict the probability of 3-year and 5-year overall survival. (C) Time-dependent receiver operating characteristic curves for nomogram-based 3-year and 5-year overall survival prediction. (D) Patient distribution based on different risk scores. (E) Kaplan–Meier survival curves of different risk groups in BC patients.

Full size image

MRP-based prognostic nomogram evaluation

Calibration plots (Fig. 6B) confirmed good consistency between the prediction of the nomogram and actual observations. The 3-year and 5-year AUC values of the nomogram for predicting OS were 0.799 (95% CI, 0.737–0.860) and 0.769 (95% CI, 0.705–0.83.4) respectively, compared with 0.708 (95% CI, 0.643–0.773, P = 0.001), and 0.650 (95% CI, 0.585–0.715, P < 0.001) for cancer stage (Fig. 6C). Hence, the MRP-based nomogram had superior predictive accuracy in BC survival outcomes. The distribution diagram of risk scores and survival status (Fig. 6D), and Kaplan–Meier curves of different risk groups (Fig. 6E) suggested the survival differences between the high- and low-risk groups.

Discussion

MRP gene family is implicated in biogenesis, metabolism, and apoptosis in BC cells. Defects in mitochondrial function underlying cancer growth may arise from abnormal gene profiling9,29. Simultaneously, they can modify gene expression and genomic stability30. To understand the aberrant MRP expression and alteration in BC, our research provided thorough bioinformatics analysis of the MRP gene family. According to the results, 12 MRPs (MRPL3, MRPL13, MRPL14, MRPL17, MRPL24, MRPL42, MRPL47, MRPS23, DAP3, MRPS30, MRPS34, MRPS35) were associated with the occurrence of BC. Consistent with previous research, higher expressions of MRPL13 and MRPS23 were revealed in BC tissues and associated with poor prognosis31,32,33,34. Interestingly, downregulation of MRPL13 restrained the proliferation and migration of BC35, whereas it enhanced the invasiveness of hepatoma via the reactive oxygen species (ROS)-Claudin-1 pathway36. We also discovered that MRP expression varied with different BC cell lines and subtypes. The lowest expressions of 12 MRP genes were observed on BC cell line SUM102PT. It was acquired from a minimally invasive human BC, characterized by overexpression of EGFR without gene amplification37. Basal-like BC, the most aggressive subtype, had the most up-regulated MRP genes. Upregulation of mRNA expression level was correlated with relevant gene amplification, the most frequent genetic alteration in MRPs, as well as DNA hypomethylation. L Liu et al. found that methylation of MRPS23 led to MRPS23 degradation and promote BC metastasis38. Our survival analysis further validated that hypermethylation of MRPS23 could result in poor survival outcomes.

In our study, the biological function of mitochondrial translational termination, elongation, and translation, as well as structural constituents of ribosome and poly(A) RNA binding were enriched in the MRP family network. But MRPs work further than that. It has been reported that MRPs constitute oxidative phosphorylation (OXPHOS) enzyme complexes, leading to reactive oxygen species (ROS) overproduction39. They can be stimulated by the mTOR signaling pathway, as well as inactivation of the TCA cycle enzyme and AMP-activated kinase (AMPK) activities40. The upregulation of MRPs in malignancy could probably be explained by the “Warburg effect”, a phenomenon that cancer cells adapt to hypoxia and display high rates of glycolysis41. Through enhanced glycolytic ATP supply, PTEN inactivation, and activated protein kinase B activity42, MRP can contribute to cancerous metabolic alterations. On the contrary, some MRPs can retard cell cycle progression and induce apoptosis by regulation of p21WAF1/CIP1, p27Kip1, and p5316,43,44. The impact and interplay of our MRP signatures deserve further investigation.

In terms of prognostic values, MRPs have been identified as effective predictive markers for metastasis and recurrence45. However, most of the previous studies were characterized by a low number of MRPs and BC samples, or by insufficient follow-up data. For the first time, our study covered nearly all MRP family members in Cox proportional hazards regression. 7 MRP genes were identified as potential biomarkers for the prognostic evaluation of BC. Upregulated MRPL1, MRPL13, MRPS6, MRPS18C, and MRPS35 brought unfavorable survival outcomes for patients with BC. On the contrary, low levels of MRPL16 and MRPL40 implied a worse prognosis. Later, the most optimal combination was selected to construct an MRP-based prognostic model for BC survival prediction. In the nomogram, low MRPL16 and MRPL40 expression, as well as high MRPS18C and MRPS35 expression were acquired with more risk scores. MRPL40 has been proven to be elevated in BC12, but their prognostic effects remain unknown. These 4 MRP genes have seldom been studied and could serve as potential therapeutic candidates for treating BC in the future.

Our study comprehensively analyzed the expression, genomic alteration, and prognostic values of the MRP family in BC. MRPL16, MRPL40, MRPS18C, and MRPS35 were integrated into a novel prognostic nomogram to accurately predict the OS in BC patients. However, there are still some limitations. Although the MRP-based nomogram indicated good accuracy and conformity, external validation was required to confirm our results. Additionally, the association between the MRP gene family and drug therapeutic effect has not been studied. Further experiments are necessary to investigate the molecular mechanism of MRPs, especially MRPL16, MRPL40, MRPS18C, and MRPS35.

Conclusion

MRP genes were primarily upregulated in BC, and MRP amplification was the most frequent alteration. They were associated with protein interactions in mitochondrial translational termination, elongation, and translation. MRPL16 and MRPL40 were positively associated with survival outcomes, while MRPS18C and MRPS35 were inversely correlated with OS. Besides thorough expression and prognosis analysis, an MRP-based nomogram was constructed. It yielded good discrimination and calibration, and therefore could accurately predict the survival outcomes for BC patients.