MRS2 and mitochondrial gene networks in endometrial cancer: mechanisms, biomarkers, and therapeutic implications

Lactic acid increases ROS production by activating MRS2

It is well-established that tumor cells frequently exhibit a preference for aerobic glycolysis due to the Warburg effect, leading to significant lactate accumulation within the tumor microenvironment. Furthermore, emerging research indicates that lactate enrichment stimulates tumor proliferation, invasion, and immune evasion (Zong et al., 2024). Consequently, monitoring lactate content serves as an effective indicator for assessing tumor proliferation levels.

To explore the potential effects of lactate levels on cancer cell metabolism in response to LPS stimulation and to validate the functional role of MRS2 in endometrial cancer, we performed a series of in vitro experiments, including lactate detection, ROS measurement, cell proliferation assays, and transwell migration assay. Our findings revealed that increased lactate levels in endometrial cancer cells (KLE) after LPS stimulation, along with upregulated MRS2 expression (Fig. 1A). Immunofluorescence staining showed that MRS2 (green fluorescence) exhibited a weak signal in KLE cells, which was significantly enhanced in the KLE + LPS group cells. The merged image, a superposition of MRS2 and nuclear staining, clearly demonstrated that LPS stimulation could promote the expression of MRS2 in KLE cells (Fig. 2A). Flow cytometry analysis revealed that LPS-stimulated KLE cells (10 µg/ml, 48 h) produced higher levels of ROS (Figs. 2B–2C). Moreover, CCK assays (Fig. 1B) and Transwell migration tests (Figs. 2D–2E) confirmed that these cells exhibited stronger proliferative and invasive tendencies (p < 0.05). Finally, transmission electron microscopy (TEM) analysis revealed that KLE cells treated with LPS for 24 h exhibited a significant increase in mitochondrial number compared to the control group (Fig. 1C). Collectively, our observations indicate a potential association between lactate enrichment and MRS2 upregulation in KLE cells. This MRS2 induction appears to impact mitochondrial function, correlating with increased ROS levels and enhanced tumor proliferative/invasive activity. These findings support a possible mechanistic link whereby lactate may promote tumor progression through MRS2-mediated mitochondrial dysfunction and subsequent ROS generation, ultimately influencing cancer cell proliferation and invasion.

MRS2 drives endometrial cancer malignancy

To further investigate the functional role of MRS2 in endometrial cancer progression, we generated an MRS2-knockdown KLE cell model using siRNA. Interestingly, MRS2 knockdown did not significantly alter lactate levels in the cell supernatant, suggesting that LPS—but not MRS2—is the primary inducer of lactate production in endometrial cancer cells, as supported by earlier findings in ‘Lactic acid increases ROS production by activating MRS2’ (Fig. 3A).

Subsequent flow cytometry analysis revealed a notable increase in ROS levels in MRS2-knockdown cells compared to controls (Fig. 3B). Furthermore, CCK-8 assays demonstrated that cell viability was significantly reduced in the siRNA-MRS2 group at both 24 and 48 h relative to the control group (p < 0.05) (Fig. 3C). These results collectively indicate that MRS2 expression is positively associated with proliferative capacity in endometrial cancer cells.

Elevated MRS2 expression in endometrial cancer correlates with poor prognosis

The RNA-seq data were downloaded from UCSC XENA (https://xenabrowser.net/datapages/) and subjected to uniform analysis. The normality test was performed on the FPKM data of the training set, yielding a p-value of 0.06726. As this value was greater than 0.05, the data were considered to be normally distributed. Analysis of TCGA data revealed significant upregulation of MRS2 in 21 tumor types, including endometrial cancer. High MRS2 expression correlated with poor tumor differentiation and worse prognosis (Table 1). MRS2 was upregulated across tumor tissues (Fig. 4A); in endometrial cancer, there was a significant difference in its expression between cancerous tissues and normal tissues (Fig. 4B). K-M curve analysis showed that high MRS2 expression was associated with low survival (Fig. 4C). ROC curve analysis yielded an area under the curve (AUC) of 0.705 for MRS2 (Fig. 4D), demonstrating favorable diagnostic efficacy and suggesting its potential as a diagnostic biomarker for endometrial cancer (Fig. 4D).

Identification of MRS2—associated mitochondrial genes in endometrial cancer

To explore the functional relevance of MRS2 in endometrial cancer, we identified differentially expressed genes (DEGs) associated with mitochondrial function. Differential expression analysis using DESeq2 identified 6,880 DEGs in endometrial cancer (|logFC| ≥ 1, adj.p < 0.05), including 4,161 upregulated and 2,719 downregulated genes (Fig. 5A). The differences in gene expression are visualized in the heatmap shown in Fig. 5B.

By analyzing the association of the MRS2 gene with genes related to mitochondrial function, 148 related genes were finally identified. The correlation results (Fig. 5C) exclusively demonstrated positive correlations, with the highest correlation coefficient being 0.6564 (MRPL3) and the lowest being 0.4008 (AFG1L).

The intersection of 6,880 DEGs and 148 MRS2-related mitochondrial function-related genes (MRS2_MFRGs) was obtained, as illustrated in Fig. 5D, yielding genes, including DARS2, DNA2, FAM136A, MGME1, MRPL13, MRPL15, MTFR2, MTHFD1L, MTHFD2, PDK1, PDSS1, PGAM5, SLC25A13, SLC25A33, VDAC1, and YARS2.

To comprehend the relevant biological functions and pathways of the candidate genes, visualization was conducted using the ggplot2 package, as demonstrated in Figs. 5E and 5F. GO analysis revealed that a total of 45 terms were enriched, including 16 cellular components (CC), seven molecular functions (MF), and 22 biological processes (BP). KEGG pathway analysis revealed that a total of seven functional pathways were significantly enriched, with most of these pathways related to mitochondrial function.

Identification of MRS2—associated mitochondrial genes most closely linked to prognosis

Univariate Cox regression analysis identified seven genes (MRPL15, MTHFD2, MTFR2, MRPL13, PDSS1, DNA2, and DARS2) significantly associated with prognosis (p < 0.05) (Fig. 6A). These genes were further analyzed to identify key prognostic markers.

The glmnet package in R was employed to conduct Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis on the seven genes identified above. The analysis selected the model with the minimum cross-validation error (lambda.min = 0.016), identifying three candidate key genes (MRPL15, MTHFD2, and MTFR2) (Fig. 6B).

Using the survival package in R, based on the expression levels of candidate key genes, EC samples were stratified into high and low expression groups (distinguished by the median), and survival analysis was conducted. The results (Fig. 6C) revealed significant differences in survival rates over time between the high and low expression groups for all three genes (p < 0.05). Therefore, MRPL15, MTHFD2, and MTFR2 were identified as the key genes in this study.

GSEA of key genes

Gene Set Enrichment Analysis (GSEA) showed that the three genes were significantly enriched in cell cycle, DNA replication, spliceosome and vascular smooth muscle contraction pathways (Fig. 7), suggesting that they may be involved in related biological processes and play a key regulatory role in the development of endometrial cancer.

Association of MRS2—related genes with immune infiltration and scoring in tumors

Single-sample gene set enrichment analysis revealed significant differences in immune cell infiltration between endometrial cancer and normal tissues (Fig. 8A). Among 28 immune cell types, 20 showed significant differences in enrichment scores (p < 0.05) (Fig. 8B), highlighting the role of immune microenvironment alterations in endometrial cancer progression. Next, we investigated the correlation between the expression levels of three key genes and MRS2 and the infiltration levels of 28 types of immune cells and visualized the results (Fig. 8C). The immune cells showing the highest positive correlation with these genes were activated CD4 T cells, effector memory CD4 T cells, and type 2 T helper cells.

Next, we evaluated EC samples in the training set using immune scores and stromal scores calculated based on the ESTIMATE algorithm. A total of three scores were predicted, namely Stromal Score (stromal score), Immune Score (immune score), and ESTIMATE Score (comprehensive score, the sum of the former two). We observed varying degrees of negative correlation between the four key genes and these three scores (p < 0.05). Among them, MRPL15 exhibited the strongest negative correlation with all three scores, while MTFR2 showed the weakest negative correlation (Fig. 9).

Transcription factors of MRS2—associated mitochondrial genes

To further explore the regulatory pathways of MRS2-related genes in endometrial cancer, we employed the ChIP-X Enrichment Analysis 3 transcription factor enrichment analysis tool (ChEA3) to predict the transcription factors associated with the three key genes. The composition of the top 10 transcription factors (CEBPG, E2F8, E2F6, CENPA, FOXM1, PA2G4, ZNF581, PRMT3, MTERF3, GTF3A) based on mean rank is illustrated in Fig. 10A. The correlation between the top 10 transcription factors and the four key genes is depicted in Fig. 10B, revealing varying degrees of positive correlation between the 10 transcription factors and the key genes. Except for the non-significant correlation between ZNF581 and MTFR2 and MTHFD2, all other relationships between the transcription factors and key genes were statistically significant (p < 0.05). The strongest correlation was observed between CENPA and MTFR2 (0.8235).

Next, we analyzed the differential expression of the top 10 transcription factors between the Normal and EC groups using the Wilcoxon test. The results (Fig. 10C) revealed significant differences in the expression of all transcription factors between the groups, except for PRMT3, which did not exhibit significant differences in expression between the two groups (p < 0.05).

Tumor drug sensitivity

A model was constructed using data from the Genomics of Drug Sensitivity in Cancer (GDSC) database (version 2), which includes data on 805 cell lines, 17,419 genes, and 198 drugs as the training set. The correlation between each drug’s IC₅₀ values and key genes was calculated, and based on the criteria of —r— > 0.3 and p < 0.05, 13 drugs with higher correlations with key genes were selected. The results are presented in Fig. 11A.

Subsequently, we stratified EC samples into high and low expression groups based on the median expression levels of key genes. We observed differences in drug IC₅₀ values between the high and low expression groups and generated a scatter plot depicting the relationship between key genes and drug IC₅₀. The results (Fig. 11B) indicated that the drug BMS_754807_2171 exhibited the strongest positive correlation with MRS2, MTFR2, MTHFD2, and MRPL15.

Protein expression validation

Subsequently, we queried the expression profiles of these four pivotal genes at the protein level in both Normal and EC samples within the Human Protein Atlas (HPA) database, comparing their disparities. The results (Fig. 12) revealed that the protein expression levels were higher in EC samples compared with normal samples, consistent with the gene expression patterns elucidated in our previous analysis.

Cell culture

Human KLE endometrial carcinoma cells (KeyGEN, Jiangsu, China) were cultured in PRMI-1640 containing 1% penicillin–streptomycin (Gibco, Billings, MT, USA) and 10% fetal bovine serum (FBS) (Gibco, Billings, MT, USA) at 37 °C in a humidified atmosphere (5% CO₂, 95% air).

Quantification of lactate content

Lactate concentration was measured enzymatically using the VITROS^® 5600 (Johnson & Johnson, New Brunswick, NJ, USA). Blood samples were immediately placed on ice after collection and centrifuged at 4 °C within 15 min to separate plasma. The enzymatic, solid-phase colorimetric method employed the VITROS Chemistry Products LAC Slide. Briefly, a 10 µL sample aliquot was dispensed onto the slide, and the subsequent enzymatic reaction, involving lactate oxidase and peroxidase, generated a colorimetric change. The rate of this change, measured by reflectance spectrophotometry at 540 nm, was kinetically monitored and quantified against a stored calibration curve.

Cell counting kit-8

The cell suspension was prepared and inoculated into a 96-well plate with six repetitive wells. After adhering to the wall, LPS (catalog numbers: MB5198, Meilunbio, Dalian, China) was added to each well and cultured in the incubator at 37 °C for a corresponding time. A total of 10 ul Cell counting kit-8 (CCK8) enhanced solution (Meilunbio, Dalian, China) was added to each well subsequently. After further incubation for 0.5–4 h, the absorbance was measured at 450 nm.

Flow cytometry detection for ROS

Cells were seeded and cultured to a confluence of 50–70% at the time of detection. A 10 µM final concentration of DCFH-DA (Meilunbio, Dalian, China) was prepared by diluting it in serum-free culture medium at a ratio of 1:1,000. The cell culture medium was then removed and a suitable volume of the diluted DCFH-DA working solution was added. The cells were incubated in the dark at 37 °C for 20–30 min. Following incubation, the cells were washed three times with serum-free culture medium to remove any uninternalized DCFH-DA. After collection, the cells were analyzed using a BD Canto II flow cytometer (Becton, Dickinson and Company, America). The fluorescence intensity was measured in real-time or at specific time points using an excitation wavelength of 488 nm and an emission wavelength of 525 nm.

Transwell migration assay

Cells were resuspended in serum-free medium at a concentration of 5 ×10⁴ cells/mL. A total of 600 µL of medium containing 20% fetal bovine serum was added to the bottom wells of a 24-well plate. A Transwell insert with an 8-µm pore membrane was placed in each well, and 200 µL of the cell suspension was added to the top chamber. The device was incubated at 37 °C with 5% CO₂ for 12–24 h. Non-migrated cells were removed, and migrated cells were fixed with 4% paraformaldehyde for 10 min. Cells were then stained with 1% crystal violet for 25 min. After washing and drying, the migrated cells were counted under a microscope to assess migration. Stained cells were imaged at 200 × magnification using an inverted microscope (Olympus, Tokyo, Japan).

Cell immunofluorescence

A 12-well plate was prepared with 5µl of culture medium, and cell coverslips were placed for subsequent fixation. Cells in the logarithmic growth phase were harvested, digested with trypsin, and then resuspended in medium for counting. Following this, the cells were inoculated into a 12-well plate. The supernatant was discarded, and the cells were washed three times with PBS. Fixation was performed by adding 4% paraformaldehyde for 15 min at room temperature, followed by three washes with PBS, each lasting 5 min. The cells were then permeabilized with 0.5% Triton X-100 for 20 min at room temperature and washed again three times with PBS for 5 min each. To block nonspecific binding, an appropriate amount of 10% goat serum blocking solution (diluted with PBS) was added, and the cells were incubated for 2 h at room temperature. After the blocking solution was removed, the primary antibody, MRS2 (catalog numbers: DF14799, 1:200, Affinity Biosciences, Cincinnati, OH, USA), was added, and the cells were incubated overnight at 4 °C. Following three washes with PBS, the secondary antibody (catalog numbers: YM2006, AbFluor 488, 1:1,000, Immunoway, San Jose, CA, USA) was added. The cells were incubated in the dark for 60 min at room temperature before the secondary antibody was removed and the cells were washed three times with PBS for 5 min each. Incubation with DAPI was performed in the dark for 5 min, followed by the removal of the DAPI and three final washes with PBS, each lasting 5 min. The coverslip was then removed, blotted dry with filter paper, and inverted onto a slide containing anti-fade reagent. The coverslip was mounted using a sealing agent. Finally, images were captured using a fluorescence microscope (Nikon, Tokyo, Japan) at a wavelength of 555 nm.

Transmission electron microscope

For transmission electron microscope (TEM) analysis, cultured cells were first fixed with 2.5% glutaraldehyde at room temperature for 5 min. Following fixation, cells were gently scraped in a unidirectional manner using a cell scraper to prevent membrane damage. The resulting cell suspension was transferred to a centrifuge tube and pelleted by low-speed centrifugation (≤3,000 rpm) for 2 min to obtain a compact cell aggregate approximately 3–4 mm in diameter. After replacing the initial fixative with fresh electron microscopy-grade fixative, the cell pellet was carefully resuspended and maintained in solution. Ultrastructural examination was subsequently performed using a HITACHI HT7700 transmission electron microscope (HITACHI, Tokyo, Japan).

Data source

Training Set: RNA-seq data (count and FPKM values) and clinical information for EC and normal samples were downloaded from the TCGA-UCEC dataset using the TCGAbiolinks package in R. After removing duplicates, 574 samples (539 EC and 35 normal) were retained for analysis. In addition, the FPKM data of the training set were tested for normality, and the obtained p-values were used to determine whether the data conformed to a normal distribution. When the p-value was greater than 0.05, the null hypothesis was rejected, indicating that the data conformed to a normal distribution.

Validation set: The GSE63678 dataset was downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc= GSE63678), consisting of 35 samples, including 12 samples of EC and their controls (EC = 7, Normal = 5). This dataset was utilized to validate the expression levels of key genes.

Mitochondrial Function-Related Genes: A total of 1,136 mitochondrial function-related genes (MFRGs) were obtained from the human MitoCarta3.0 database (https://www.broadinstitute.org). Besides, 1,967 mitochondrial function-related genes were acquired from the Integrated Mitochondrial Protein Index database (IMPI, https://www.mrc-mbu.cam.ac.uk/). After taking the intersection of the two sets, 1,070 mitochondrial function-related genes were identified.

Differential expression analysis

Differential gene expression analysis was performed using the DESeq2 package (1.36.0) in R (Love, Huber & Anders, 2014), with significant genes identified based on |logFC| ≥ 1 and adj.p < 0.05. Gene expression patterns were visualized using heatmaps generated with the pheatmap package.

Candidate gene identification

To explore the role of MRS2 and mitochondrial function-related genes in the pathogenesis of EC, using the psych package in R, the relationship between the MRS2 gene and mitochondrial function-related genes was analyzed based on Pearson correlation. Mitochondrial function-related genes associated with the MRS2 gene were selected using the criteria of |r| > 0.4 and p < 0.05.

Functional enrichment analysis

Functional enrichment analysis of disease-associated differentially expressed genes, including gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, was performed using the clusterProfiler package (4.4.4) in R (Wu et al., 2021). Significantly enriched functional entries and pathways were identified with a significance threshold of p < 0.05.

Survival analysis

To screen for key genes significantly associated with survival, univariate Cox regression analysis (survival package (3.4-0) (MGrambsch, 2000) in R, p < 0.05) identified prognosis-related genes. LASSO regression analyses were performed using the glmnet software package (4.1-4) (Friedman, Hastie & Tibshirani, 2010), and the parameter with the smallest error was determined by 10-fold cross-validation (lambda.min = 0.016), and the results were screened accordingly to further screen candidate genes. To evaluate the prognostic significance of the key genes, the survival (3.4-0) (MGrambsch, 2000) and survminer packages (0.5.0) in R were used to classify the candidate key genes into high and low expression groups according to the median expression value. Kaplan–Meier survival curves were plotted, and genes with a significant correlation with survival (p < 0.05) were identified as the final key genes.

Gene set enrichment analysis

To understand the biological functions and signaling pathways associated with key genes, samples were divided into high and low expression groups based on the median expression level of marker genes, followed by differential expression analysis. GSEA was then conducted using the C2: KEGG gene set from the MSigDB database, obtained via the msigdbr R package, as the background reference, and enrichment analysis was conducted on the ranked differentially expressed genes in the background gene set using the GSEA function in R (adj.p < 0.05).

Analysis of immune infiltration

Immune cell infiltration levels were analyzed using ssGSEA (GSVA package in R), with enrichment scores calculated for 28 immune cell types. Wilcoxon tests were used to compare infiltration levels between EC and normal samples (A p-value less than 0.05 indicated statistical significance). The ggdotchart function of the ggpubr package (0.4.0) was used to visualize the correlation between three key gene expression levels and 28 levels of immune cell infiltration in a lollipop plot.

ESTIMATE score

The estimate package in R was used to estimate the following scores in the EC samples of the training set: StromalScore, ImmuneScore, and ESTIMATEScore, with the ESTIMATEScore being the sum of StromalScore and ImmuneScore. A p-value < 0.05 indicated that the difference in the level of immune infiltration between the groups was statistically significant.

Transcription factor identification and drug sensitivity screening

Transcription factors associated with the key genes were predicted using the ChEA3 tool (https://maayanlab.cloud/chea3/), and the top 10 significantly enriched were screened according to a significance threshold of p < 0.05 to analyze the differences in distribution between groups and the correlation between transcription factors and expression levels of key genes. This study utilized the oncoPredict software package (0.2) (Maeser, Gruener & Huang, 2021) in R to predict the IC₅₀ values of multiple drugs in EC samples. The research employed version 2 data from the Cancer Drug Sensitivity Genomics Database (GDSC) as the training dataset, which contains 805 cell lines, 17,419 genes, and 198 drugs. A correlation model was established to analyze the relationship between drug HIC values and key genes. Drugs showing significant correlations with key genes (as indicated by correlation coefficients |r| > 0.3 and p < 0.05) were identified through computational analysis.

Protein expression validation

To validate our findings at the protein level, immunohistochemical staining data for the key genes in endometrial tissue were retrieved from the Human Protein Atlas (HPA) database (https://www.proteinatlas.org/). The protein expression patterns were compared with transcriptional data to evaluate their consistency.

Statistical analysis

Statistical analyses were performed using IBM SPSS Statistics 20 (IBM Corp., Armonk, NY, USA) GraphPad Prism 8, and R 3.5.1. Data are presented as mean ± SEM. Statistical comparisons between three or more groups were performed using one-way ANOVA for normally distributed data or the Kruskal–Wallis test for non-normally distributed data. Pearson correlation was used for correlation analysis, with p < 0.05 considered statistically significant.