Bioinformatics and experimental approach reveal potential prognostic and immunological roles of key mitochondrial metabolism-related genes in cervical cancer

Background: Metabolic remodeling is the hallmark of cancer. In recent years, mitochondrial metabolism (MM) has been considered essential in tumorigenesis and cancer progression. Understanding the role of MM in cervical cancer (CC) can provide insights into disease progression and potential therapeutic targets.

Methods: Clinical data of CC patients was downloaded from the UCSC Xena dataset, and differentially expressed genes (DEGs) were identified between tumor and normal samples. MM-related genes (MMRGs) were screened from the MSigDB database. DEGs and MMRGs were then intersected to identify differentially expressed MMRGs. A prognostic risk model was constructed based on these intersecting genes through Cox regression analysis, and its association with the tumor microenvironment and immune checkpoint-related genes was evaluated. Hub genes’ expression was evaluated in cells through qRT-PCR. Additionally, drug sensitivity analysis was conducted to explore potential therapeutic drugs.

Results: We identified 259 overlapping genes between DEGs and MMRGs, with 55 being prognosis-related. Two molecular clusters were revealed, with C1 exhibiting poorer prognosis. A prognostic risk model comprising five genes (BDH1, MIR210, MSMO1, POLA1, and STARD3NL) was established, showing significant associations with survival outcomes of CC patients. Functional enrichment analysis revealed that DEGs between high- and low-risk groups were tightly associated with the immune system. Analysis of the immune microenvironment showed significant differences between different risk groups, with higher immune and ESTIMATE scores observed in the low-risk group. Additionally, expression levels of immune checkpoint-related genes were significantly correlated with the risk score. Drug sensitivity analysis identified potential therapeutic agents correlated with the expression of the five prognostic genes.

Conclusion: Our findings underscore the importance of MM in CC progression and provide potential therapeutic targets for CC.

1 Introduction

Cervical cancer (CC) is one of the most prevalent gynecologic cancers, with an incidence of approximately 6.5% of all female cancer cases worldwide (1, 2). It is primarily caused by persistent infection with high-risk human papillomavirus (HPV) types (3). Despite advancements in screening programs and HPV vaccination efforts, CC remains a leading cause of cancer-related deaths among women globally (4). For patients with early or locally advanced CC, the 5-year survival rate exceeds 50% following surgical or chemoradiotherapy interventions (5, 6). However, the metastasis or recurrence significantly reduces survival rates, with a 5-year survival rate of only 10% for patients under such circumstances (4, 5). Therefore, it is urgent to delve deeper into the biological mechanisms of CC progression and identify novel prognostic biomarkers to enhance therapeutic strategies and patient outcomes.

Mitochondria, responsible for cellular energy generation, provide energy through the tricarboxylic acid (TCA) cycle and oxidative phosphorylation (OXPHOS) (7). Cellular energy imbalance is a recognized hallmark of cancer (8). In the 1920s, Otto Warburg postulated that cancer cells preferentially utilize glycolysis over OXPHOS for ATP production (9). For a long time, the major metabolic feature of tumor cells was considered to be the Warburg effect (10). However, in recent years, increasing evidence suggests that mitochondrial metabolism (MM) and function play a crucial role in tumorigenesis and cancer progression. Dysregulated mitochondrial function in cancer cells leads to alterations in energy production, metabolism, and redox balance, facilitating their proliferation, survival, and metastasis (11, 12). Moreover, mitochondria actively regulate apoptosis, allowing cancer cells to evade cell death signals and promote resistance to therapy (13–15). Therapeutic methods that target diverse pathways within MM, such as inhibiting cellular constituents engaged in mitochondrial synthesis, decreasing metabolite accumulation, or preventing energy production within mitochondria, have demonstrated therapeutic efficacy in various cancers (16–18). Several regents have been reported to alleviate CC by modulating MM. For example, Ginsenoside Rh2 induces CC cell apoptosis by suppressing mitochondrial electron transfer chain complex (19). Butyrate inhibited mitochondria-dependent apoptosis in CC cells (20). Therefore, MM could be a new therapeutic strategy for CC. Understanding the intricate interplay between MM and cancer biology holds promise for identifying novel therapeutic targets and developing personalized treatment strategies to combat CC.

This study aimed to elucidate the role of MM in CC progression and prognosis. Differentially expressed mitochondrial metabolism-related genes (MMRGs) in CC patients were identified and a prognostic risk model based on these genes was constructed. Through comprehensive bioinformatics analyses, we explored the clinical significance of MMRGs, as well as their association with tumor immune microenvironment and drug sensitivity, providing insights into potential therapeutic strategies for CC.

2 Methods

2.1 Acquirement of differentially expressed MMRGs

CC cohort was downloaded from the TCGA-UCSC Xena (http://xena.ucsc.edu). DEGs were identified between tumor samples (n = 305) and normal samples (n = 4) using the “limma” package in R software (version 4.1.3), with the criterion of |log2 fold change (FC)| > 1 and adjust p-value < 0.05. These DEGs were visualized in a volcano plot using “ggplot2”. MMRGs were acquired from the MSigDB database (gsea-msigdb.org). Overlap genes in the DEGs and MMRGs were visualized through “upsetR” and “VennDiagram” packages.

2.2 Functional enrichment analysis

Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using the “clusterProfiler”, with a criterion of p < 0.05. The results of GO and KEGG were visualized through the “Goplots” package in R.

2.3 Consensus clustering analysis

The MMRGs underwent univariate Cox regression analysis using the SPSS (p < 0.05). Consensus clustering analysis was then performed using R package ConsensusClusterPlus.

2.4 Establishment and evaluation of prognostic risk model

To mitigate overfitting, LASSO regression analysis was conducted on genes selected through univariate Cox regression utilizing the R package “glmnet”, with the penalty function lambda (0.5765) employed via cross-validation to identify and eliminate overfitting genes. Finally, a multivariate Cox regression analysis was performed on the retained genes using SPSS to establish a prognostic risk model for MMRGs in CC. The risk score was calculated using the obtained regression coefficients from the multivariate Cox regression analysis, with a formula as follows:

$$riskscore={\displaystyle \sum _{i=1}^{n}ex{p}_{RNAi}}\ast Coe{f}_{RNAi}$$

Subsequently, the risk score for each sample in the TCGA dataset was calculated. Based on their RiskScore values, samples were categorized into high- and low-risk groups, with the median RiskScore serving as the threshold. Kaplan-Meier (K-M) curves and receiver operating characteristic (ROC) curves were then generated for both groups using the R packages “survival” and “survminer” to assess the predictive performance of the model. The area under the ROC curve (AUC) for 1-, 3-, and 5-years overall survival was calculated. The expression of prognostic biomarkers in each sample was visualized in a heatmap using the “ggplot2” package in R.

2.5 Construction of nomogram

Clinical characteristics encompass age, TNM stage, and pathology stage. A nomogram was established based on these clinical features and risk scores. The calibration curve depicted the 45-degree dashed lines representing the best predictions of the nomogram. The nomogram and calibration curve was constructed using “rms”, “regplot”, and “survival” packages.

2.6 Evaluation of tumor immune microenvironment landscape

The stromal score, immune score, ESTIMATE score, and tumor purity were acquired using the ESTIMATE algorithm. They were then compared between different risk groups using the Wilcoxon rank-sum test. Immune cell abundance in the high- and low-risk groups was detected using the MCPcounter algorithm.

2.7 Immune checkpoint and immunotherapy response analysis

Immune checkpoint-related genes in CC were searched from the Checkpoint Therapeutic Target Database (CKTTD) (21). Expression levels of these genes in different risk groups were explored, and the correlation of them with risk score was visualized using the R package “ggplot2”. Tumor immune dysfunction and exclusion (TIDE) score was determined using the TIDE algorithm. p < 0.05 indicates statistical significance.

2.8 Drug sensitivity analysis

CellMiner (https://discover.nci.nih.gov/cellminer) encompasses 60 distinct cell lines originating from 9 types of malignancies, serving as essential screens in the development of novel anti-tumor medications. It includes 262 drugs, either FDA-approved or undergoing clinical trials (22). NCI60 drug response data were acquired from the CellMiner tool. Drug sensitivity between different risk groups was evaluated, and the association between drug sensitivity and prognostic gene expression was assessed using Pearson’s test.

2.9 Cell culture

Human cervical epithelial cells (Cat NO.: CP-H059) and cervical cancer cells Hela S3 (Cat NO.: CL-0350) were purchased from Wuhan Pricella Biotechnology Co., Ltd. All cells were cultured in a specialized medium containing Ham’s F-12K supplemented with 10% FBS and 1% P/S in an incubator at 37°C with 5% CO₂. The medium was changed every 2–3 days.

2.10 Quantitative real-time (qRT)-PCR

Total RNA was extracted using TRIzol reagent (Invitrogen) according to the manufacturer’s protocol. The concentration and purity of the extracted RNA were assessed using a NanoDrop spectrophotometer (Thermo Fisher Scientific, CA, USA). The isolated RNA was then reverse-transcribed into cDNA with the PrimeScript RT reagent kit (Takara, Dalian, China). The reaction conditions were as follows: 42°C for 15 min and 95°C for 3 min. For quantification, qRT-PCR was performed using the Hieff UNICON Universal Blue qPCR SYBR Green Master Mix (Yeasen, Shanghai, China) on the ABI7900HT System. The thermocycling program was set as follows: initial denaturation at 95°C for 30 s, followed by 40 cycles of 95°C for 3 s and 60°C for 20 s. The relative mRNA expression levels were determined using the 2^-ΔΔCt method, with GAPDH serving as the internal control. The primer sequences were designed by Gene Creat Bioengineering CO. (Wuhan, China) and were shown in Supplementary Table S1.

2.11 Statistical analysis

Statistical analyses were conducted using R software version 4.2.0 and GraphPad Prism 8.0.2. Student t-test was used to compare the differences between the two groups, and the Pearson method was used for correlation analysis. All experiments were performed in triplicate, and data were analyzed using GraphPad Prism software. Statistical significance was determined using a Student’s t-test. Statistical significance was defined as p < 0.05.

3 Results

3.1 Identification of prognosis-related differentially expressed MMRGs in patients with CC

Based on the TCGA dataset, 4353 DEGs were identified between the tumor and normal samples. Volcano plot showed 1976 upregulated genes and 2377 downregulated genes in the tumor group (Figure 1A). Additionally, 1234 MMRGs were screened from the MsigDB database (Figure 1B). Venn diagram revealed 259 overlapping genes in the DEGs and MMRGs (Figure 1B). Then, the overlapping genes were seeded on the univariate Cox regression, identifying 55 prognosis-related MMRGs (Figure 1C). GO enrichment analysis showed that these genes were associated with biological process, such as regulation of small GTPase mediated signal transduction, small molecule metabolic process, cellular lipid metabolic process, and small molecule catabolic process. They were also related to cellular components, including AP-type membrane coat adaptor complex, cytosol, and membrane coat. As for molecular functions, metallopeptidase activity, metalloendopeptidase activity, and catalytic activity were enriched (Figure 1D). The enriched KEGG pathways were related to environmental information processing, metabolism, and organismal systems, such as ABC transporters, fatty acid degradation, and bile secretion (Figure 1E).

Figure 1

Figure 1. Identification of prognosis-related differentially expressed MMRGs in patients with CC. (A) Volcano plot of DEGs between tumor and normal samples in the TCGA database. The blue dots represent downregulated genes and the red dots represent upregulated genes. (B) A total of 4353 DEGs and 1234 MMRGs were screened from public databases, and 259 intersecting genes were identified between the DEGs and MMRGs using the Venn diagram. (C) Prognostic features were selected using the univariate Cox regression based on 259 differentially expressed MMRGs; p < 0.05. (D, E) GO and KEGG enrichment analysis was performed on the genes identified by univariate Cox regression analysis. CC, cervical cancer; DEG, differentially expressed gene; MMRG, mitochondrial metabolism-related gene; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

3.2 Identification of two CC molecular subtypes

Subsequently, consensus clustering analysis was performed on the 55 prognosis-related MMRGs. The CDF curve showed that k = 2 is the optimal number of clusters (Figures 2A, B). Figure 2C displayed the consensus values for different k values, with the highest consensus value at k = 2. Therefore, all samples were clustered into two subtypes. The heatmap showed that the samples were well separated, with a clear distinction between the two subtypes (Figure 2D). K-M curves revealed that patients in C1 had significantly poorer OS than those in C2 (Figure 2E). Compared to C2, patients in C1 had higher pathological stages (Table 1). These results suggest that patients in C1 may have a worse prognosis and a higher degree of malignancy.

Figure 2

Figure 2. Two molecular subtypes of CC were identified. (A) Consensus CDF curves for different values of k (k = 2 to 10). (B) Relative change in the area under the CDF curve for different k values, indicating the optimal number of clusters. (C) Bar plot showing consensus values for different k values, with different colors representing varying cluster numbers. (D) Consensus clustering heatmap showing the division of samples into two clusters (C1 and C2) based on consensus matrix values. (E) K-M curve showed the overall survival time of CC patients between C1 and C2. CDF, cumulative distribution function; CC, cervical cancer; K-M, Kaplan-Meier; HR, hazard ratio; CI, confidence interval; OS, overall survival.

Table 1

Table 1. Differences in stages and grades between two different subtypes.

We then compared the expression levels of 55 prognosis-related MMRGs between C1 and C2. As shown in Figure 3A, 17 MMRGs were significantly overexpressed and 25 MMRGs were downregulated in C1. KEGG pathways analysis revealed that both overexpressed genes and downregulated genes in C1 were associated with metabolic pathways (Figures 3B, C). These results indicated that metabolic pathways may be involved in the imbalance of these differentially expressed prognosis-related genes to regulate the CC tumor microenvironment and progression.

Figure 3

Figure 3. Expression and functions of prognosis-related genes in C1 and C2. (A) Expression levels of 55 prognosis-related genes in C1 and C2. (B) KEGG pathways related to overexpressed genes in C1. (C) KEGG pathways related to downregulated genes in C1. *p < 0.05, **p < 0.01, ***p < 0.001, “-”means no significant.

3.3 Establishment of the MM-related risk model for patients with CC

To further explore the role of 55 MMRGs in CC prognosis, lasso regression was performed on these genes, and 10 MMRGs were identified (Figure 4A). The 10-round cross-validation was used to determine the optimal values of the penalty parameter (Figure 4B). Finally, using multivariate Cox regression, five genes, including BDH1, MIR210, MSMO1, POLA1, and STARD3NL were identified (p < 0.05, Figure 4C). Results revealed that BDH1 (HR = 0.71, 95%CI: 0.530-0.951) and POLA1 (HR = 0.581, 95%CI: 0.357-0.946) were protective features, while MIR210 (HR = 1.702, 95%CI: 1.326-2.186), MSMO1 (HR = 1.592, 95%CI: 1.158-2.189), and STARD3NL (HR = 1.776, 95%CI: 1.199-2.631) were harmful features (Figure 4C).

Figure 4

Figure 4. Identification of prognostic genes in CC. (A, B) Lasso Cox regression identified 10 genes related to the prognosis of patients with CC, and 10-round cross-validation was performed to detect the optimal values of the penalty parameter. (C) Multivariate Cox regression identified 5 prognostic genes based on the above 10 genes; p < 0.05. CC, cervical cancer; MMRG, mitochondrial metabolism-related gene.

Furthermore, based on the expression levels of the above 5 prognostic genes, the risk score of each sample in the TCGA dataset was calculated as follows: risk score = -0.343*BDH1 + 0.532*MIR20 + 0.465*MSMO1 + -0.543*POLA1 + 0.575*STARD3NL. CC patients in the TCGA dataset were divided into high- and low-risk groups based on the median risk score. The risk score distribution of patients was visualized in Figure 5A. Our data indicated that patients with high risk scores had worse prognoses (Figure 5A). The expression of BDH1 and POLA1 was higher in the low-risk group, while the expression of MIR210, MSMO1, and STARD3NL was higher in the high-risk group (Figure 5A). The K-M survival curve demonstrated superior survival outcomes for the low-risk group compared to the high-risk group (Figure 5B). The AUCs of the enrolled patients at 1-, 3-, and 5-year were 0.78, 0.77, and 0.75, respectively (Figure 5C).

Figure 5

Figure 5. Assessment of MM-related risk model for CC. (A) Risk score distribution, OS time of each patient in the TCGA cohort, and heatmap of the five prognostic genes. (B) Kaplan-Meier curve revealed the survival probability of patients in different risk groups. (C) The ROC curve of the risk model at 1-, 3-, and 5-year. Red represents 1-year, green represents 3-year, and blue represents 5-year. (D) Nomogram for predicting 1-, 3-, and 5-year OS time for CC patients in the TCGA database. (E) Calibration curves of the nomogram observed 1-, 3-, and 5-year outcomes. MM, mitochondrial metabolism; CC, cervical cancer; OS, overall survival; ROC, receiver operating characteristic.

A nomogram was then constructed to show the performance of risk score and clinical features on CC prognosis. As shown in Figure 5D, the M stage, pathology stage, and risk score ranked in the top three in terms of contribution to predicting CC, followed by T stage, G stage, N stage, and age. The calibration curves of the nomogram for the probability at 3- and 5-year indicated a good clinical value (Figure 5E).

3.4 Functional enrichment analysis of DEGs between high-risk and low-risk groups

Subsequently, DEGs between high- and low-risk groups were identified, and the function of these genes was explored using GO and KEGG enrichment analyses. The top five enriched GO biological process terms were immune system process, response to stress, regulation of response to stimulus, response to chemical, and system development (Figure 6A). The top five GO cellular component terms included extracellular region, vesicle, endomembrane system, intrinsic component of membrane, and integral component of membrane (Figure 6B). The GO molecular function terms were mainly enriched in signaling receptor binding, protein-containing complex binding, molecular function regulator, anion binding, and small molecule binding (Figure 6C). The KEGG pathways were major enriched in cell adhesion molecules, Th17 cell differentiation, IL-17 signaling pathway, T cell receptor signaling pathway, and MAPK signaling pathway (Figure 6D). These results revealed that DEGs between high- and low-risk cohorts were tightly associated with tumor immune microenvironment.

Figure 6

Figure 6. Functions of DEGs between high and low risk groups. (A–C) GO enrichment analysis was performed on DEGs in high and low risk groups, including biological process (A), cellular component (B), and molecular function (C). (D) KEGG pathways related to DEGs in high and low risk groups. DEG, differentially expressed gene; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

3.5 Tumor immune microenvironment analysis

Given the above results that differentially expressed prognosis-related genes were related to metabolic pathways and that DEGs in high risk vs. low risk were related to the immune system process, we then explore the immune microenvironment between different risk groups. Our data showed that immune and ESTIMATE scores in the low-risk group were significantly higher than those in the high-risk group (p < 0.05, Figure 7A), and the high-risk group had a higher tumor purity (Figure 7B). Additionally, we observed five immune cells were significantly infiltrated in the low-risk group, including T cells, CD8 T cells, cytotoxic lymphocytes, B lineage, and myeloid dendritic cells (Figure 7C).

Figure 7

Figure 7. Different immune profiles between different risk groups in the TCGA dataset. (A) Stromal, immune, and ESTIMATE scores between the low- and high-risk groups. (B) Tumor purity between the low- and high-risk groups. (C) Cell fraction between the low- and high-risk groups. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001, ns, no significant.

The effector function of CD8 T cells is regulated by immune checkpoints. Due to the differences in the levels of CD8 T cells between high and low risk groups, we further analyzed the differences in immune checkpoint gene expression between high and low risk groups. Among 13 checkpoint-related genes, 7 genes were observed to be significantly correlated to the risk score, and all of them were downregulated in the high-risk group, including CDK1, EZH2, ICOS, IDO1, PLK1, TIGIT, and TLR8 (p < 0.05, Figure 8A). Further research showed that the high-risk group had a higher TIDE score and exclusion score than the low-risk group (Figure 8B). These results suggest that the high risk group is associated with an immunosuppressive environment.

Figure 8

Figure 8. Analysis of immune checkpoints between the low- and high-risk groups. (A) Expression levels of immune checkpoint-related genes between the low- and high-risk groups. (B) TIDE score between the low- and high-risk groups. *p < 0.05, **p < 0.01, ***p < 0.001. TIDE, tumor immune dysfunction and exclusion; MSI, microsatellite instability.

3.6 Expression of hub genes

To further explore whether the differences in the microenvironment between different groups are related to the hub genes, we examined the expression of these hub genes among different groups. Compared to the low-risk group, BDH1 and POLA1 were downregulated in the high-risk group, while MIR210, MSMO1, and STARD3NL were upregulated (Figure 9A). Similarly, in the C1 group, BDH1 and POLA1 were downregulated, while MIR210, MSMO1, and STARD3NL were upregulated (Figure 9B). These results are consistent with previous findings, indicating a higher mortality risk and therefore a poorer prognosis in the C1 group. Furthermore, we conducted cell experiments to validate the expression of these genes. Compared to human cervical epithelial cells, BDH1 and POLA1 were downregulated, while MIR210, MSMO1, and STARD3NL were upregulated in the CC cell line Hela S3 cells (Figure 9C). These findings suggest that the expression of these five hub genes may be closely related to tumor immune-metabolic regulation.

Figure 9

Figure 9. Expression levels of 5 MMRGs in different groups. (A) Expression levels of 5 prognostic genes in the low- and high-risk groups. (B) Expression levels of 5 prognostic genes in C1 and C2. (C) Expression levels of 5 prognostic genes in human cervical epithelial cells and CC cells. *p < 0.05, **p < 0.01, ***p < 0.001. MMRG, mitochondrial metabolism-related gene; CC, cervical cancer.

3.7 Drug sensitivity analysis

The correlation between the expression levels of 5 prognostic genes and drug sensitivity was further investigated utilizing the CellMiner database. The top three drugs significantly correlated with the expression of each gene are shown in Supplementary Figure S1. BDH1 is positively associated with ciclosporin, Raloxifene, and Tamoxifen (Supplementary Figure S1A). MIR210 was positively associated with Cediranib, ergenyl, and Motesanib (Supplementary Figure S1B). Similarly, MSMO1 displayed a positive correlation with Amiodarone, uridin, and Zoledronate (Supplementary Figure S1C), and POLA1 was positively associated with Methylprednisolone, PX-316, and ZM-336372 (Supplementary Figure S1D). STARD3NL showed a positive association with JNJ-38877605 and Lovastatin while displaying a negative association with Fluorouracil (Supplementary Figure S1E). Among these drugs, the high-risk group was significantly sensitized to certain drugs, including ciclosporin, Raloxifene, Tamoxifen, Zoledronate, and Lovastatin (Supplementary Figure S2).

4 Discussion

Although advances in the treatment of CC have been made, its poor prognosis still poses a significant threat to women’s health (4). Recent studies suggest that MM is essential for tumor growth, and some clinical trials have demonstrated the feasibility of modulating MM to treat cancer (23). The bioinformatic research on MM in CC is limited. Therefore, this study constructed a prognostic model based on five MMRGs for CC and explored their association with the immune microenvironment.

Metabolic remodeling is one of the hallmarks of cancer. Current evidence indicates that MM-related pathways are reprogrammed in cancer, playing crucial roles in bioenergetics, biosynthesis, and redox homeostasis (24). The regulation of redox balance in tumor cells is influenced by their significantly increased glucose uptake, which produces TCA cycle metabolites. These metabolites supply electrons to the mitochondrial electron transport chain (ETC) (25). Inhibiting ETC-related genes could heighten the vulnerability of cancer cells to glucose depletion, consequently impeding tumor progression (26). MMRGs have been considered prognostic markers for various cancers, including breast cancer (27), osteosarcoma (28), and ovarian cancer (29). Based on the MMRGs. the present study identified two molecular subtypes of CC, and C1 showed shorter OS time than C2. Additionally, 5 key prognostic MMRGs in CC were identified, including BDH1, MIR210, MSMO1, POLA1, and STARD3NL. BDH1 is a key catalytic enzyme in ketogenesis, catalyzing the reversible conversion of acetoacetate to beta-hydroxybutyrate (30). Within the mitochondria, ketone bodies undergo oxidation via the TCA cycle, leading to the generation of acetyl-CoA and NADH. Downregulation of BDH1 is a prognostic marker in hepatocellular carcinoma (31). POLA1 encodes DNA polymerase, which facilitates DNA replication and repair, ensuring the maintenance of mitochondrial genome integrity. POLA1 has antitumor activity in inhibiting cancer cell proliferation and inducing apoptosis (32). A previous study demonstrates that POLA overexpression is associated with the poor prognosis of CC patients (33). MIR210 originates from mitochondria, its expression correlates with hypoxia gene signatures, and it could reduce the activity of proteins controlling MM (34). Nakada et al. points out that MIR210 induces energy metabolism shift from OXPHOS to glycolysis via acting on the mitochondrial inner membrane (35). MSMO1 catalyzes the demethylation of C4-methyl sterol, a critical step in cholesterol biosynthesis within mitochondria. Abnormal expression of MSMO1 would lead to CC (36). STARD3NL is involved in MM by mediating the transfer of cholesterol between membranes, potentially contributing to lipid metabolism and homeostasis within mitochondria (37). The risk score constructed using these five MMRGs exhibits good prognostic function for patients with CC. Patients within the low-risk group have longer survival time.

Furthermore, KEGG analysis of DEGs in the normal and tumor groups was enriched for several metabolic pathways. Overexpressed genes and downregulated genes in C1 were also associated with metabolic pathways. Metabolism and immunity are usually inextricably linked. After further analyzing the function of DEGs in high- and low-risk populations, we found that these DEGs are primarily associated with the immune system, including the Th17 cell differential, IL-17 signaling pathway, and MAPK signaling pathway. Th17 cells are a subset of T-helper cells that produce IL-17, a pro-inflammatory cytokine (38). IL-17 acts on tumor cells and various components of the tumor microenvironment to promote tumor growth, angiogenesis, and metastasis (39). It can also induce the production of other pro-inflammatory cytokines and chemokines, creating a pro-tumor inflammatory environment (40). The MAPK signaling pathway is a crucial intracellular signaling cascade involved in cell proliferation, survival, and differentiation (41). Activation of MAPK signaling can occur downstream of IL-17 receptor engagement. In breast cancer, the Th17/IL-17/MAPK cascade signaling pathway plays a multifaceted role in cancer progression (42). Previous studies have indicated that MM abnormalities can influence tumor cell antigen presentation and processing, thereby aiding tumor cells in evading recognition and attack by the immune system (43). Activation of the IL-17/MAPK signaling pathway can lead to an increase in immunosuppressive cells such as regulatory T cells and myeloid-derived suppressor cells (MDSCs), thereby inhibiting the anti-tumor immune response. Additionally, activation of the IL-17/MAPK signaling pathway may result in changes in tumor cell surface antigens, making tumor cells less recognizable and susceptible to clearance by the immune system (44). Th17 and the MAPK signaling pathway have been demonstrated to play a pro-econogenic role in promoting CC (45, 46). Therefore, we hypothesize that MMRGs may enhance tumor immune evasion by modulating the IL-17/MAPK signaling pathway in cervical cancer.

Given the crucial role of the immune microenvironment in cancer progression, we further investigated differences in immune cell infiltration and immune scores among individuals at different risk groups. The results indicated that individuals at low risk had higher immune scores, and correspondingly, we observed more immune cell infiltration in the low-risk group, including T cells, CD8 T cells, cytotoxic lymphocytes, B lineage cells, and myeloid dendritic cells. These findings demonstrate the correlation between MM-related prognostic models and immune infiltration in CC. To regulate immune responses, PD-L1 was expressed in these immune cells to maintain immune homeostasis and protect the body from foreign pathogens (47). Furthermore, analysis of immune checkpoint-related genes revealed higher expression of immune checkpoints in the low-risk group, suggesting better efficacy of immune therapy in this group compared to the high-risk group. Among the identified immune checkpoint genes, TIGIT showed the most significant difference between the two groups. TIGIT is a crucial target in tumor immunotherapy, as it can prevent NK cells from releasing tumor antigens, impair dendritic cell-induced T cell priming, or inhibit CD8+ T cell-mediated killing of cancer cells (48). This is consistent with the findings of Han et al., suggesting that TIGIT may kill cancer cells in the low-risk group by reducing dendritic cell-triggered T-cell priming (49). However, the specific mechanism of action of TIGIT in cervical cancer remains to be further elucidated. Additionally, Our results revealed that BDH1 expression is positively related to ciclosporin, a typical immunosuppressive drug with an anti-tumor effect (50). Importantly, the drug sensitivity of ciclosporin is higher in the high-risk group, indicating that CC patients have better outcomes with this drug. Research suggests that alterations in MM within cancer cells may also modulate drug sensitivity and resistance through the MAPK signaling pathway (51). Whether MM-modulated MAPK signaling influences the drug sensitivity of ciclosporin needs further research.

5 Conclusion

In conclusion, this study elucidated the significance of MM in CC progression and its interplay with the immune microenvironment. By constructing a prognostic model based on five MMRGs (BDH1, MIR210, MSMO1, POLA1, and STARD3NL) and exploring their association with immune infiltration, significant insights were gained into CC treatment. Furthermore, the study highlighted the potential involvement of the Th17/IL-17/MAPK signaling pathway in mediating immune escape in CC, possibly influenced by MMRGs. Additionally, the analysis of immune checkpoint-related genes suggested a potential for improved efficacy of immune therapy in the low-risk group, with TIGIT emerging as a significant target. These findings underscore the intricate relationship between MM, immune regulation, and therapeutic outcomes in CC, providing valuable insights for the development of novel prognostic markers and therapeutic strategies.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Ethics statement

Ethical approval was not required for the study involving humans in accordance with the local legislation and institutional requirements. Written informed consent to participate in this study was not required from the participants or the participants’ legal guardians/next of kin in accordance with the national legislation and the institutional requirements.

Author contributions

QH: Conceptualization, Data curation, Formal analysis, Writing – original draft. Y-FX: Investigation, Methodology, Writing – original draft. H-PL: Data curation, Methodology, Writing – original draft. TZ: Formal analysis, Supervision, Writing – original draft, Writing – review & editing.

Funding

The author(s) declare that no financial support was received for the research and/or publication of this article.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2025.1522910/full#supplementary-material.

Supplementary Figure 1 | Correlation between predicted drugs and prognostic genes. (A) Correlation of BDH1 with ciclosporin, Raloxifene, and Tamoxifen. (B) Correlation of MIR210 with Cediranib, ergenyl, and Motesanib. (C) Correlation of MSMO1 with Amiodarone, uridin, and Zoledronate. (D) Correlation of POLA1 with Methylprednisolone, PX-316, and ZM-336372. (E) Correlation of STARD3NL with JNJ-38877605, Lovastatin, and Fluorouracil.

Supplementary Figure 2 | Drug sensitivity between different risk groups. *p < 0.05, **p < 0.01, ns, no significant.

References

1. Jha AK, Mithun S, Sherkhane UB, Jaiswar V, Osong B, Purandare N, et al. Systematic review and meta-analysis of prediction models used in cervical cancer. Artif Intell Med. (2023) 139:102549. doi: 10.1016/j.artmed.2023.102549

PubMed Abstract | Crossref Full Text | Google Scholar

2. Ming C, Bai X, Zhao L, Yu D, Wang X, Wu Y. RPL24 as a potential prognostic biomarker for cervical cancer treated by Cisplatin and concurrent chemoradiotherapy. Front Oncol. (2023) 13:1131803. doi: 10.3389/fonc.2023.1131803

PubMed Abstract | Crossref Full Text | Google Scholar

3. Schubert M, Bauerschlag DO, Muallem MZ, Maass N, Alkatout I. Challenges in the diagnosis and individualized treatment of cervical cancer. Medicina (Kaunas). (2023) 59:925. doi: 10.3390/medicina59050925

PubMed Abstract | Crossref Full Text | Google Scholar

4. Han S, Wang S, Lv X, Li D, Feng Y. Ferroptosis-related genes in cervical cancer as biomarkers for predicting the prognosis of gynecological tumors. Front Mol Biosci. (2023) 10:1188027. doi: 10.3389/fmolb.2023.1188027

PubMed Abstract | Crossref Full Text | Google Scholar

5. Li C, Cang W, Gu Y, Chen L, Xiang Y. The anti-PD-1 era of cervical cancer: achievement, opportunity, and challenge. Front Immunol. (2023) 14:1195476. doi: 10.3389/fimmu.2023.1195476

PubMed Abstract | Crossref Full Text | Google Scholar

6. Lv X, Jia Y, Li J, Deng S, Yuan E. The construction of a prognostic model of cervical cancer based on four immune-related LncRNAs and an exploration of the correlations between the model and oxidative stress. Front Pharmacol. (2023) 14:1234181. doi: 10.3389/fphar.2023.1234181

PubMed Abstract | Crossref Full Text | Google Scholar

7. Wang SF, Tseng LM, Lee HC. Role of mitochondrial alterations in human cancer progression and cancer immunity. J BioMed Sci. (2023) 30:61. doi: 10.1186/s12929-023-00956-w

PubMed Abstract | Crossref Full Text | Google Scholar

8. Huang L, Wei B, Zhao Y, Gong X, Chen L. DYNLT1 promotes mitochondrial metabolism to fuel breast cancer development by inhibiting ubiquitination degradation of VDAC1. Mol Med. (2023) 29:72. doi: 10.1186/s10020-023-00663-0

PubMed Abstract | Crossref Full Text | Google Scholar

9. Foo BJ, Eu JQ, Hirpara JL, Pervaiz S. Interplay between mitochondrial metabolism and cellular redox state dictates cancer cell survival. Oxid Med Cell Longev. (2021) 2021:1341604. doi: 10.1155/2021/1341604

PubMed Abstract | Crossref Full Text | Google Scholar

10. Kubicka A, Matczak K, Labieniec-Watala M. More than meets the eye regarding cancer metabolism. Int J Mol Sci. (2021) 22:9507. doi: 10.3390/ijms22179507

PubMed Abstract | Crossref Full Text | Google Scholar

11. Dong L, Gopalan V, Holland O, Neuzil J. Mitocans revisited: mitochondrial targeting as efficient anti-cancer therapy. Int J Mol Sci. (2020) 21:7941. doi: 10.3390/ijms21217941

PubMed Abstract | Crossref Full Text | Google Scholar

12. Passaniti A, Kim MS, Polster BM, Shapiro P. Targeting mitochondrial metabolism for metastatic cancer therapy. Mol Carcinog. (2022) 61:827–38. doi: 10.1002/mc.23436

PubMed Abstract | Crossref Full Text | Google Scholar

13. Dabravolski SA, Nikiforov NG, Zhuravlev AD, Orekhov NA, Mikhaleva LM, Orekhov AN. The role of altered mitochondrial metabolism in thyroid cancer development and mitochondria-targeted thyroid cancer treatment. Int J Mol Sci. (2021) 23:460. doi: 10.3390/ijms23010460

PubMed Abstract | Crossref Full Text | Google Scholar

14. Sainero-Alcolado L, Liano-Pons J, Ruiz-Perez MV, Arsenian-Henriksson M. Targeting mitochondrial metabolism for precision medicine in cancer. Cell Death Differ. (2022) 29:1304–17. doi: 10.1038/s41418-022-01022-y

PubMed Abstract | Crossref Full Text | Google Scholar

15. Li W, Zhang L. Rewiring mitochondrial metabolism for CD8(+) T cell memory formation and effective cancer immunotherapy. Front Immunol. (2020) 11:1834. doi: 10.3389/fimmu.2020.01834

PubMed Abstract | Crossref Full Text | Google Scholar

16. Shrestha R, Johnson E, Byrne FL. Exploring the therapeutic potential of mitochondrial uncouplers in cancer. Mol Metab. (2021) 51:101222. doi: 10.1016/j.molmet.2021.101222

PubMed Abstract | Crossref Full Text | Google Scholar

17. Mamouni K, Kallifatidis G, Lokeshwar BL. Targeting mitochondrial metabolism in prostate cancer with triterpenoids. Int J Mol Sci. (2021) 22:2466. doi: 10.3390/ijms22052466

PubMed Abstract | Crossref Full Text | Google Scholar

18. Paragliola RM, Torino F, Papi G, Locantore P, Pontecorvi A, Corsello SM. Role of mitotane in adrenocortical carcinoma - review and state of the art. Eur Endocrinol. (2018) 14:62–6. doi: 10.17925/EE.2018.14.2.62

PubMed Abstract | Crossref Full Text | Google Scholar

19. Liu Y, Yu S, Xing X, Qiao J, Yin Y, Wang J, et al. Ginsenoside Rh2 stimulates the production of mitochondrial reactive oxygen species and induces apoptosis of cervical cancer cells by inhibiting mitochondrial electron transfer chain complex. Mol Med Rep. (2021) 24:873. doi: 10.3892/mmr.2021.12513

PubMed Abstract | Crossref Full Text | Google Scholar

20. Zhang K, Ji X, Song Z, Song W, Huang Q, Yu T, et al. Butyrate inhibits the mitochondrial complex Iota to mediate mitochondria-dependent apoptosis of cervical cancer cells. BMC Complement Med Ther. (2023) 23:212. doi: 10.1186/s12906-023-04043-3

PubMed Abstract | Crossref Full Text | Google Scholar

21. Zhang Y, Yao Y, Chen P, Liu Y, Zhang H, Liu H, et al. Checkpoint therapeutic target database (CKTTD): the first comprehensive database for checkpoint targets and their modulators in cancer immunotherapy. J Immunother Cancer. (2020) 8:e001247. doi: 10.1136/jitc-2020-001247

PubMed Abstract | Crossref Full Text | Google Scholar

22. Shankavaram UT, Varma S, Kane D, Sunshine M, Chary KK, Reinhold WC, et al. CellMiner: a relational database and query tool for the NCI-60 cancer cell lines. BMC Genomics. (2009) 10:277. doi: 10.1186/1471-2164-10-277

PubMed Abstract | Crossref Full Text | Google Scholar

23. Vasan K, Werner M, Chandel NS. Mitochondrial metabolism as a target for cancer therapy. Cell Metab. (2020) 32:341–52. doi: 10.1016/j.cmet.2020.06.019

PubMed Abstract | Crossref Full Text | Google Scholar

24. Boese AC, Kang S. Mitochondrial metabolism-mediated redox regulation in cancer progression. Redox Biol. (2021) 42:101870. doi: 10.1016/j.redox.2021.101870

PubMed Abstract | Crossref Full Text | Google Scholar

25. Oliveira GL, Coelho AR, Marques R, Oliveira PJ. Cancer cell metabolism: Rewiring the mitochondrial hub. Biochim Biophys Acta Mol Basis Dis. (2021) 1867:166016. doi: 10.1016/j.bbadis.2020.166016

PubMed Abstract | Crossref Full Text | Google Scholar

26. Sun L, Wan A, Zhou Z, Chen D, Liang H, Liu C, et al. RNA-binding protein RALY reprogrammes mitochondrial metabolism via mediating miRNA processing in colorectal cancer. Gut. (2021) 70:1698–712. doi: 10.1136/gutjnl-2020-320652

PubMed Abstract | Crossref Full Text | Google Scholar

27. Lin Y, Huang Z, Zhang B, Yang H, Yang S. Construction and analysis of a mitochondrial metabolism-related prognostic model for breast cancer to evaluate survival and immunotherapy. J Membr Biol. (2024) 257:63–78. doi: 10.1007/s00232-024-00308-1

PubMed Abstract | Crossref Full Text | Google Scholar

28. Yang S, Liu L, Liu X, Li X, Zheng Y, Ren Z, et al. The mitochondrial energy metabolism pathway-related signature predicts prognosis and indicates immune microenvironment infiltration in osteosarcoma. Med (Baltimore). (2023) 102:e36046. doi: 10.1097/MD.0000000000036046

PubMed Abstract | Crossref Full Text | Google Scholar

29. Fu Y, Huang Z, Huang J, Xiong J, Liu H, Wan X. Metabolism-related gene vaccines and immune infiltration in ovarian cancer: A novel risk score model of machine learning. J Gene Med. (2024) 26:e3568. doi: 10.1002/jgm.3568

PubMed Abstract | Crossref Full Text | Google Scholar

30. Zhang Z, Bi X, Lian X, Niu Z. BDH1 promotes lung cancer cell proliferation and metastases by PARP1-mediated autophagy. J Cell Mol Med. (2023) 27:939–49. doi: 10.1111/jcmm.17700

PubMed Abstract | Crossref Full Text | Google Scholar

31. Luo W, Wu S, Zhang F, Chen X, Ma Y, Mo Y. Decreased expression of 3-hydroxybutyrate dehydrogenase 1 is a prognostic marker and promotes tumor progression in hepatocellular carcinoma. Pathol Res Pract. (2022) 238:154111. doi: 10.1016/j.prp.2022.154111

PubMed Abstract | Crossref Full Text | Google Scholar

32. Han T, Goralski M, Capota E, Padrick SB, Kim J, Xie Y, et al. The antitumor toxin CD437 is a direct inhibitor of DNA polymerase alpha. Nat Chem Biol. (2016) 12:511–5. doi: 10.1038/nchembio.2082

PubMed Abstract | Crossref Full Text | Google Scholar

33. Yu L, Wei M, Li F. Longitudinal analysis of gene expression changes during cervical carcinogenesis reveals potential therapeutic targets. Evol Bioinform Online. (2020) 16:1176934320920574. doi: 10.1177/1176934320920574

PubMed Abstract | Crossref Full Text | Google Scholar

34. Chan SY, Zhang YY, Hemann C, Mahoney CE, Zweier JL, Loscalzo J. MicroRNA-210 controls mitochondrial metabolism during hypoxia by repressing the iron-sulfur cluster assembly proteins ISCU1/2. Cell Metab. (2009) 10:273–84. doi: 10.1016/j.cmet.2009.08.015

PubMed Abstract | Crossref Full Text | Google Scholar

35. Nakada C, Hijiya N, Tsukamoto Y, Yano S, Kai T, Uchida T, et al. A transgenic mouse expressing miR-210 in proximal tubule cells shows mitochondrial alteration: possible association of miR-210 with a shift in energy metabolism. J Pathol. (2020) 251:12–25. doi: 10.1002/path.5394

PubMed Abstract | Crossref Full Text | Google Scholar

36. Zheng G, Wang Z, Fan Y, Wang T, Zhang L, Wang M, et al. The clinical significance and immunization of MSMO1 in cervical squamous cell carcinoma based on bioinformatics analysis. Front Genet. (2021) 12:705851. doi: 10.3389/fgene.2021.705851

PubMed Abstract | Crossref Full Text | Google Scholar

37. Xu Y, Bao X, Chen X, Wu P, Chen S, Zhang B, et al. STARD3NL inhibits the osteogenic differentiation by inactivating the Wnt/beta-catenin pathway via binding to Annexin A2 in osteoporosis. J Cell Mol Med. (2022) 26:1643–55. doi: 10.1111/jcmm.17205

PubMed Abstract | Crossref Full Text | Google Scholar

38. Yasuda K, Takeuchi Y, Hirota K. The pathogenicity of Th17 cells in autoimmune diseases. Semin Immunopathol. (2019) 41:283–97. doi: 10.1007/s00281-019-00733-8

PubMed Abstract | Crossref Full Text | Google Scholar

39. Nalbant A. IL-17, IL-21, and IL-22 cytokines of T helper 17 cells in cancer. J Interferon Cytokine Res. (2019) 39:56–60. doi: 10.1089/jir.2018.0057

PubMed Abstract | Crossref Full Text | Google Scholar

40. Salazar Y, Zheng X, Brunn D, Raifer H, Picard F, Zhang Y, et al. Microenvironmental Th9 and Th17 lymphocytes induce metastatic spreading in lung cancer. J Clin Invest. (2020) 130:3560–75. doi: 10.1172/JCI124037

PubMed Abstract | Crossref Full Text | Google Scholar

41. Wei J, Liu R, Hu X, Liang T, Zhou Z, Huang Z. MAPK signaling pathway-targeted marine compounds in cancer therapy. J Cancer Res Clin Oncol. (2021) 147:3–22. doi: 10.1007/s00432-020-03460-y

PubMed Abstract | Crossref Full Text | Google Scholar

42. Shibabaw T, Teferi B, Ayelign B. The role of Th-17 cells and IL-17 in the metastatic spread of breast cancer: As a means of prognosis and therapeutic target. Front Immunol. (2023) 14:1094823. doi: 10.3389/fimmu.2023.1094823

PubMed Abstract | Crossref Full Text | Google Scholar

43. Guo Y, Luo C, Sun Y, Guo W, Zhang R, Zhang X, et al. Inhibition of mitochondrial fusion via SIRT1/PDK2/PARL axis breaks mitochondrial metabolic plasticity and sensitizes cancer cells to glucose restriction therapy. BioMed Pharmacother. (2023) 166:115342. doi: 10.1016/j.biopha.2023.115342

PubMed Abstract | Crossref Full Text | Google Scholar

44. Lin Z, Huang Q, Liu J, Wang H, Zhang X, Zhu Z, et al. Interleukin-17D promotes lung cancer progression by inducing tumor-associated macrophage infiltration via the p38 MAPK signaling pathway. Aging (Albany NY). (2022) 14:6149–68. doi: 10.18632/aging.204208

PubMed Abstract | Crossref Full Text | Google Scholar

45. Liu Y, Guo QF, Chen JL, Li XR, Hou F, Liu XY, et al. Correlations between alterations of T-helper 17 cells and treatment efficacy after concurrent radiochemotherapy in locally advanced cervical cancer (stage IIB-IIIB): a 3-year prospective study. Chin Med J (Engl). (2021) 134:954–62. doi: 10.1097/CM9.0000000000001475

PubMed Abstract | Crossref Full Text | Google Scholar

46. Sun HN, Xie DP, Ren CX, Guo XY, Zhang HN, Xiao WQ, et al. Ethyl beta-carboline-3-carboxylate increases cervical cancer cell apoptosis through ROS-p38 MAPK signaling pathway. In Vivo. (2022) 36:1178–87. doi: 10.21873/invivo.12817

PubMed Abstract | Crossref Full Text | Google Scholar

47. Sun C, Mezzadra R, Schumacher TN. Regulation and function of the PD-L1 checkpoint. Immunity. (2018) 48:434–52. doi: 10.1016/j.immuni.2018.03.014

PubMed Abstract | Crossref Full Text | Google Scholar

48. Liang JY, Wang DS, Lin HC, Chen XX, Yang H, Zheng Y, et al. A novel ferroptosis-related gene signature for overall survival prediction in patients with hepatocellular carcinoma. Int J Biol Sci. (2020) 16:2430–41. doi: 10.7150/ijbs.45050

PubMed Abstract | Crossref Full Text | Google Scholar

49. Ge X, Su Z, Wang Y, Zhao X, Hou K, Zheng S, et al. Identifying the intervention mechanisms of polydatin in hyperuricemia model rats by using UHPLC-Q-Exactive Orbitrap mass spectroscopy metabonomic approach. Front Nutr. (2023) 10:1117460. doi: 10.3389/fnut.2023.1117460

PubMed Abstract | Crossref Full Text | Google Scholar

50. Wu PJ, Hsin IL, Hung WL, Lee MS, Wang PH, Ko JL. Combination treatment with cyclosporin A and arsenic trioxide induce synergistic cell death via non-apoptotic pathway in uterine cervical cancer cells. Chem Biol Interact. (2022) 368:110177. doi: 10.1016/j.cbi.2022.110177

PubMed Abstract | Crossref Full Text | Google Scholar

51. Vu MH, Iswanto ABB, Lee J, Kim JY. The role of plasmodesmata-associated receptor in plant development and environmental response. Plants (Basel). (2020) 9: 216. doi: 10.3390/plants9020216

PubMed Abstract | Crossref Full Text | Google Scholar