Construction of a prognostic model based on efferocytosis for predicting survival and immunotherapy efficacy in cervical cancer

Introduction

Cervical cancer (CC) represents a malignancy originating in the cervix and ranks as the second most prevalent gynecological malignancy, and is the foremost contributor to cancer-related mortality among women (1). Persistent infection with high-risk human papillomavirus (HPV) is a significant risk factor for CC. While the immune system typically clears most HPV infections within a few months, certain high-risk strains, such as HPV16 and HPV18, evade immune clearance, resulting in the continued expression of viral oncogenes E6 and E7. This persistent oncogene expression induces heightened genomic instability, the accumulation of somatic mutations, and the integration of HPV DNA into the host genome, all of which contribute to CC (2).

CC can be effectively prevented through HPV vaccination, which serves as primary prevention and is primarily targeted at prepubertal and adolescent girls. In addition, cervical screening provides secondary prevention in women. Current treatment options for CC include surgical resection, radiotherapy, chemotherapy, and immunotherapy (3). CC-related mortality rates in developed countries have declined owing to the continual advancements in medical technology and the widespread implementation of CC vaccination programs (4). Nevertheless, the incidence of CC continues to rise globally, particularly among younger patients (5). It has been reported that 5% to 26% of patients with early-stage disease experience relapse following standardized treatment, and the 5-year survival rate for those with recurrent disease varies greatly, ranging from 15.0% to 50.0% (6). The delayed onset of symptoms until the advanced stages of the disease poses significant challenges to treatment efficacy (7). Hence, identifying and investigating genes associated with CC prognosis are crucial for improving patient outcomes.

Conventional prognostic indicators for CC currently used include age, pathological staging, histological type, and lymph node metastasis. Despite their utility in assessing prognostic risk to a certain extent, their predictive accuracy remains limited. Therefore, there is a critical need to identify novel biomarkers that can complement and enhance the predictive power of traditional clinical indicators.

High-throughput omics technologies provide a more comprehensive and detailed understanding of tumor biology by analyzing a large number of genes and proteins, including potential markers and signaling pathways related to prognosis, thereby enabling a more holistic prognostic risk assessment. Recent advances in high-throughput data analysis methods have revealed a correlation between CC prognosis and the tumor immune microenvironment (TIME). The TIME represents a complex ecosystem where tumor cells interact with their surroundings, playing a crucial role in the onset and progression of tumors as well as affecting clinical outcomes in various immunotherapy approaches. Efferocytosis is one of the most frequent events within the tumor microenvironment (TME) but its role has not yet been fully explored in relation to CC prognosis.

It is widely recognized that tumor development results from mutations in tumor-related genes, epigenetic alterations, or the accumulation of both in normal cells, and CC is no exception. Efferocytosis refers to the process of efficiently clearing apoptotic cells (ACs) without triggering an inflammatory response, which is a fundamental physiological mechanism for maintaining tissue development and homeostasis (8). The specific mechanism involves the conversion of methionine derived from ACs into S-adenosylmethionine (SAM) during efferocytosis, which serves as a substrate for DNA methylation mediated by recombinant DNA methyltransferase 3A (DNMT3A). This process plays a crucial role in tissue homeostasis by affecting the methylation level of relevant genes (9).

In addition, cytokines released during efferocytosis can promote genomic mutations to a certain extent, thereby promoting tumor progression and affecting treatment outcomes. Studies have shown that cell debris generated from the apoptosis of tumor cells during efferocytosis can stimulate tumor-associated macrophages (TAMs) to release cytokines (e.g., MMP, EGF, TGF-β, and so on), which induce epithelial-mesenchymal transition (EMT) to promote tumor progression. Simultaneously, efferocytosis can promptly remove damage-associated molecular patterns (DAMPs) released by damaged cells, thus promoting the immune escape of tumor cells.

Immunosuppressive cytokines released during efferocytosis can further inhibit the release of pro-inflammatory cytokines (such as IL-2, INF-γ, and TNF, among others), thus enhancing immune escape of tumor cells by altering the TME. In addition, cytokines secreted during efferocytosis can also promote the invasion and metastasis of tumor cells to distant tissues. For example, M2 macrophages induced by efferocytosis can stimulate the high expression of tumor-invasive cytokines such as MFG-E8 and NOX2 and secrete chitinase 3-like 1 (CHI3L1) to promote metastasis in breast and gastric cancers (10).

In CC, primary fibroblasts are capable of transforming HPV-positive CC cells through efferocytosis, thereby enabling cellular immortality (11). Furthermore, efferocytosis-related molecules and pathways have been identified as potential targets for cancer therapy, including phosphatidylserine and calreticulin, Tyro3, Axl, and Mer tyrosine kinase (MerTK), receptors of tyrosine kinase, indoleamine-2,3-dioxygenase 1, annexin V, CD47, TGF-β, IL-10, and macrophage phenotype switching. Combination therapies targeting these molecules alongside conventional therapies have demonstrated enhanced efficacy in cancer treatment (12). However, the potential of these mechanisms in promoting the progression, invasion, and metastasis of CC remains to be fully understood.

In this study, the prognostic predictive ability of ERGs in patients with CC was systematically evaluated using regression analysis methods based on publicly available CC databases. A risk model incorporating ERGs was developed and validated by assessing the expression of key genes in clinical samples using reverse transcription-quantitative polymerase chain reaction (RT-qPCR). The primary objective of this study was to effectively differentiate prognostic outcomes among patients with CC and optimize the predictive performance of current prognostic indicators. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1146/rc).

Methods

Data acquisition

The GSE9750 dataset, comprising 33 CC samples and 24 healthy controls (HC), was retrieved from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) to identify differentially expressed genes (DEGs). Additionally, data from The Cancer Genome Atlas (TCGA) database (https://tcga-data.nci.nih.gov/)—specifically the TCGA-CC cohort—was accessed to provide information on 293 tumor samples, including clinical data, gene expression profiles, and survival outcomes, which was crucial for the construction of the prognostic risk model. To further validate this model, the Cancer Genome Characterization Initiative (CGCI): HIV+ Tumor Molecular Characterization Project (HTMCP)-CC dataset, containing survival information and gene expression profiles of 117 CC cases, was also used. A total of 67 ERGs were acquired based on previous research (13,14). Additionally, genes related to invasion, EMT, and angiogenesis genes were downloaded from the dbEMT2 and cancer single-cell state atlas (CancerSEA) databases, as well as the Molecular Signatures Database (MSigDB), respectively.

Identification of differentially expressed ERGs (DE-ERGs) from the GSE9750 dataset

The “limma” package was used to identify DEGs between CC and HC samples in the GSE9750 dataset (15), with a threshold of P.adjusted <0.05. Subsequently, the DEGs were overlapped with previously identified ERGs to extract DE-ERGs. The g:Profiler tool (https://biit.cs.ut.ee/gprofiler/gost) was utilized for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of these DE-ERGs. The outcomes were visualized using bubble plots and Sankey diagrams generated using the R “ggplot2” package. Additionally, pathways associated with DE-ERGs enriched in the HALLMARK and KEGG datasets were further analyzed through gene set variation analysis (GSVA) (16).

Development of a prognostic risk model for CC based on characteristic genes

Prognosis-related genes were identified through univariate Cox regression analysis of DE-ERGs within the TCGA-CC dataset with a significance threshold of P<0.05. Subsequently, the least absolute shrinkage and selection operator (LASSO) and multivariate Cox regression analyses were performed to select the most predictive characteristic genes (17). Based on the median risk score, individuals in the TCGA-CC cohort were stratified into two risk groups. The differences in expression of characteristic genes and overall survival (OS) between the high- and low-risk groups were analyzed utilizing the Wilcoxon test and Kaplan-Meier survival curves. The receiver operating characteristic (ROC) curves, generated using the “survROC” package (18), were used to assess the predictive accuracy of the risk model. The robustness of the prognostic risk model was further validated using the external CGCI-HTMCP-CC dataset.

Correlation analysis of clinicopathological attributes and risk scores in the TCGA-CC cohort

The following clinicopathological characteristics of individuals in the TCGA-CC cohort were evaluated: age, TNM (Tumor, Node, Metastasis) stage, tumor stage, G stage, invasion markers, and radiation therapy. The association between risk scores and these clinicopathological variables was assessed using the Wilcoxon and Kruskal-Wallis tests. Additionally, survival differences among patients with varying clinical features were examined using Kaplan-Meier survival analysis.

Analysis of independent prognostic factors in the TCGA-CC Cohort

Univariate and multifactorial Cox regression analyses were conducted to determine whether clinicopathological features and risk scores served as independent prognostic factors for individuals with CC. The diagnostic efficacy of individual or combined independent prognostic factors was evaluated using ROC curves. A nomogram predicting survival probability based on independent prognostic factors was constructed utilizing the “rms” (v6.3-0) (https://CRAN.R-project.org/package=rms). Calibration curves and decision curve analysis (DCA) were employed to validate whether the nomogram could serve as an optimal model for clinical decision-making.

Relationship between immunotherapy response of patients with CC and TIME

To specifically validate the predictive value of our risk model for immunotherapy response, we utilized the TCGA-CC cohort as the immunotherapy cohort. The potential response of patients to immune checkpoint inhibitors was inferred using the Tumor Immune Dysfunction and Exclusion (TIDE) score (comprising TIDE, dysfunction, and exclusion scores) and the immunophenoscore (IPS). A lower TIDE score and a higher IPS indicate a better potential response to immunotherapy. The estimation algorithm was utilized for comparing the stromal, immune, and ESTIMATE scores, as well as tumor purity, between the two risk groups. Pearson analysis was employed to explore the relationship between the TME and the risk score. The “immunedeconv” package was used to assess immune cell infiltration across all TCGA-CC samples (19), followed by Pearson analysis to evaluate correlations between immune cell populations. Additionally, the Wilcoxon test was performed to assess variations in immune cell infiltration between the two risk groups, while Spearman correlation analysis was utilized to examine the relationship between characteristic genes and immune cells. Additionally, the TIDE score was used to predict the likelihood of response to immunotherapy in CC, and its association with the risk score was evaluated using Pearson correlation.

Correlation of risk scores with metastasis/invasion indicators

The association between risk scores and genes related to invasion, EMT, and angiogenesis was analyzed using Pearson correlation analysis. The mRNA stemness indexes (mRNAsi) for patients with CC were sourced from an earlier study (20). Additionally, the invasion, neovascularization, EMT, and mRNAsi characteristics for the TCGA-CC cohort were evaluated using single-sample gene set enrichment analysis (ssGSEA).

Investigation of the regulatory mechanisms of characteristic genes in CC

Methylation regulators associated with m⁶A, m¹A, m⁵C, and m⁷G modifications were sourced from prior studies (Table 1). The relationship between these methylation regulators and the characteristic genes was analyzed through Pearson correlation analysis. Predicted interactions between the characteristic genes and their corresponding microRNAs (miRNAs) and lncRNAs (long non-coding RNAs) were identified using mirtarbase and tarbase databases, with predictions conducted with the R package “multiMiR” (21). The lncRNA-miRNA-mRNA network was subsequently constructed and visualized using Cytoscape software.

Validation of characteristic gene expression

Ten cancerous tissue samples were collected from patients diagnosed with HPV16-positive cervical squamous cell carcinoma (stages IA–IIA) at the Obstetrics and Gynecology Clinic of the hospital between May 2022 and January 2023. This sample cohort consisted of 3 cases at stage IA2, 6 cases at stage IB, and 1 case at stage IIA1. As per the pathological classification, there were 6 cases of well-differentiated cancer, while 4 cases were moderately differentiated. All tissue specimens were subjected to pathological diagnosis. The ages of patients ranged from 40 to 60 years, with a mean age of 49.5±2.39 years. Post-surgery, the samples were stored in an −80 ℃ ultra-low temperature freezer. Exclusion criteria included pregnancy, concurrent malignancies, prior preoperative radiotherapy or chemotherapy, and a history of cervical or vaginal lesion treatments. Clinical data and research results were kept confidential, and their use was limited to this study. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by the ethics committee of the Second Hospital of Shanxi Medical University of Shanxi Province, China (No. 2022-YX-226) and informed consent was obtained from all individual participants.

RT-qPCR was used to validate the expression of characteristic genes in all specimens. Total RNA was isolated from 10 CC and 10 adjacent tissue samples using the TRIzol reagent, and reverse transcription was performed using the SureScript First-Strand cDNA Synthesis Kit (Servicebio, China). Subsequently, RT-qPCR was conducted using Universal Blue SYBR Green qPCR Master Mix under the following thermal cycling conditions: initial denaturation at 95 ℃ for 60 seconds, denaturation at 95 ℃ for 20 seconds, annealing at 55 ℃ for 20 seconds, followed by amplification for 40 cycles. Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) served as the housekeeping gene. Primer sequences are presented in Table 2. The 2^-ΔΔCT method was employed to calculate gene expression levels, with GAPDH serving as the reference gene for normalization.

Statistical analysis

Continuous variables were tested for normality using the Kolmogorov-Smirnov test and presented as mean ± standard deviation or median values. Between-group comparisons used independent t-tests or Wilcoxon rank-sum tests, while categorical data were analyzed with chi-square tests. Repeated-measures analysis of variance (ANOVA) assessed within-group changes over time. Survival differences were evaluated with Kaplan-Meier analysis, and prognostic significance was assessed using univariate and multivariate Cox regression. ROC curves and area under the curve (AUC) metrics quantified predictive accuracy. Statistical significance was defined as P<0.05.

Results

Identified DE-ERGs

A total of 4,578 DEGs were identified in the GSE9750 dataset, comprising 2,293 upregulated DEGs and 2,285 downregulated DEGs, when comparing CC samples to HC (Figure 1A,1B). From these, 19 DE-ERGs were identified using a Venn diagram analysis (Figure 1C). These DE-ERGs demonstrated enrichment in 281 GO entries and 4 KEGG pathways (Figure 1D).

Figure 1 Differential expression analysis in the GSE9750 dataset. (A) Volcano plot and (B) heatmap of 4,578 DEGs between CC sample and HC samples. (C) Venn diagrams for the 19 differentially expressed efferocytosis-related genes (DE-ERGs) in CC. (D) A bubble chart for visualizing the number of enrichment entries of the 19 DE-ERGs. Sankey diagrams for the enriched MF (E), BP (F), CC (G) entries, and KEGG pathways (H) associated with the 19 DE-ERGs. BP, biological processes; CC, cellular component; DE-ERGs, differentially expressed efferocytosis-related genes; DEGs, differentially expressed genes; FC, fold change; GO, Gene Ontology; HC, healthy control; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function.

In terms of molecular function (MF), the DE-ERGs were involved in protein binding (Figure 1E), while biological processes (BP) analysis revealed their involvement in response to stimuli and cell communication (Figure 1F). Cellular components (CC) analysis indicated that these genes were associated with cellular anatomical entities (Figure 1G). Furthermore, DE-ERGs exhibited enrichment in arginine and proline metabolic pathways (Figure 1H).

Further analysis showed that the 19 DE-ERGs were predominantly enriched in 28 HALLMARK entries (Figure 2A) and 88 KEGG pathways (Figure 2B). Within HALLMARK gene sets, fatty acid metabolism, Notch signaling, and angiogenesis emerged as the main enriched pathways for DE-ERGs. In KEGG, the natural killer cell-mediated pathway exhibited the strongest association with the DE-ERGs.

Figure 2 Gene set variation analysis of the 19 differentially expressed efferocytosis-related genes. Bar chart illustrating the pathways enriched with DE-ERGs in HALLMARK (A) and KEGG (B) databases. DE-ERGs, differentially expressed efferocytosis-related genes; FC, fold change; KEGG, Kyoto Encyclopedia of Genes and Genomes.

EPO and UCP2 identified as key prognostic risk factors

Among the 19 DE-ERGs initially examined in the TCGA-CC dataset, four prognostic-related DE-ERGs were identified through univariate Cox analysis (Figure 3A). Subsequent LASSO regression analysis was performed to eliminate false positive genes (Figure 3B,3C), followed by multifactorial Cox analysis, which narrowed the selection to two key prognostic ERGs: EPO and UCP2 (Figure 3D).

Figure 3 Screening of two characteristic genes used to develop a prognostic risk model in CC. (A) Univariate Cox analysis to screen four prognostic-related differentially expressed efferocytosis-related genes in TCGA-CC dataset. (B) Two characteristic genes were identified using the LASSO model. (C) Cross-validation for tuning parameter selection in the LASSO model. (D) Forest plot of multivariate Cox analysis. (E) Distribution of risk curves and survival states for CC in the TCGA dataset. (F) Kaplan-Meier curves of distinct risk groups. (G) ROC curves assessing the predictive model for survival prediction in individuals with CC at 1-, 3-, and 5-year. ***, P<0.001. AUC, area under the curve; CC, cervical cancer; LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.

In the GSE9750 dataset, no significant difference was found in the expression of UCP2 between the CC and HC samples. However, EPO expression was higher in healthy samples compared to tumor samples (Figure S1A). Validation using RT-qPCR analysis further confirmed that both EPO and UCP2 expression levels were significantly higher in healthy tissues than in tumor tissues (Figure S1B,S1C).

A prognostic risk score was calculated using the following formula: Risk score = −0.136868067 × EPO − 0.398364215 × UCP2.

Based on the median risk score of 1.042299, patients were categorized into high-risk and low-risk groups (Figure 3E). Patients in the low-risk group exhibited a significantly better OS compared to those in the high-risk group (Figure 3F). The ROC curve for OS showed strong predictive capability of the risk signature, with AUC values exceeding 0.60 at 1, 3, and 5 years, highlighting the model’s robust prognostic power (Figure 3G). Furthermore, the predictive power of the risk model was further validated in the external CGCI-HTMCP-CC dataset, confirming its efficacy across datasets (Figure S2A-S2D).

Significant differences in clinicopathological characteristics

Significant differences in risk scores were observed among patients based on T-stage, tumor stage, and radiotherapy status, with all showing statistical significance (P<0.05; Figure 4A-4C). Stratified survival analysis based on clinical data revealed that patient survival differed significantly in specific clinical subgroups, including T1/T2, N0 stage, M0 stage, age <60 years, and tumor stage I/II (Figure 4D), highlighting the prognostic relevance of these clinical factors in CC.

Figure 4 Correlation of risk score and clinicopathological attributes. (A-C) Boxplots for risk score in distinct clinicopathological subgroups of TCGA-CC dataset, including (A) T stage, (B) tumor stage, and (C) radiation therapy. (D) Stratified survival analysis of different clinicopathological characteristics, including age <60 years, tumor stage I/II, M0, N0, and T1/T2. TCGA-CC, The Cancer Genome Atlas-cervical cancer.

Prognostic efficacy based on functional assessment of the nomogram model

Lymphovascular infiltration and the risk score were identified as independent prognostic factors for CC through univariate and multifactorial Cox analyses (Figure 5A). A nomogram model incorporating these independent prognostic factors was developed (Figure 5B), demonstrating that an increase in the overall score corresponded with a decline in survival rate. The model’s accuracy was confirmed using a calibration curve with a C-index of 0.942, indicating high reliability in survival predictions (Figure 5C). DCA results further supported the efficacy of the nomogram as an optimal tool for clinical decision-making (Figure 5D). The combined diagnostic model, which integrated risk score and lymphovascular infiltration, exhibited higher AUC values compared to diagnostic models relying solely on individual risk factors. The nomogram exhibited enhanced diagnostic accuracy and effectiveness in predicting survival outcomes (Figure 5E).

Figure 5 Clinical utilization of the risk score-based nomogram in CC. (A) Univariate Cox analysis was utilized to identify independent prognostic factors in TCGA-CC dataset. (B) A nomogram was developed as per lymphovascular infiltration and risk score. (C) Calibration curve of nomogram. (D) The clinical impact curve compares nomogram and other prognostic factors. (E) Decision curve analysis of nomogram, lymphovascular infiltration, and risk score in CC. AUC, area under the curve; OS, overall survival; TCGA-CC, The Cancer Genome Atlas-cervical cancer.

Immune cell variations between the high- and low-risk groups

Significant differences in immune cell profiles were observed between the high- and low-risk groups. The high-risk group showed notable correlations with immune score, ESTIMATE score, stromal score, and tumor purity, while no such correlations were observed in the low-risk group (Figure S3A). Additionally, there were varying abundances of 22 immune cells between both groups (Figure 6A), with their correlation indicating an association among immune cells in CC (Figure S3B). Notably, nine immune cells exhibited considerable variations between the two groups (Figure 6B). The strongest positive correlation was noted between EPO and resting mast cells, while the strongest negative correlation existed between UCP2 and activated mast cells (Figure 6C,6D).

Figure 6 Immune correlation analyses of two characteristic genes in CC. (A) Histogram of immune cell infiltration in TCGA-CC dataset. (B) Boxplot for 22 immune cells proportions in high- and low-risk samples. (C) Lollipop plot illustrating the link between nine differential immune cells and EPO. (D) Lollipop plot illustrating the link between nine differential immune cells and UCP2. *, P<0.05; **, P<0.01; ****, P<0.0001. TCGA-CC, The Cancer Genome Atlas-cervical cancer.

Furthermore, with respect to immunotherapy suitability in risk groups, individuals in the low-risk group appeared more likely to be responsive to immunotherapy (Figure 7A), as evidenced by a strong negative correlation between TIDE scores and the high-risk group (Figure 7B).

Figure 7 Immunotherapy exploration in high- and low-risk cervical cancer groups. (A) Boxplots of TIDE Score and IPS in the high- and low-risk groups. (B) Scatter plots of the correlations of dysfunction, exclusion, and TIDE scores with risk scores. IPS, International Prognostic Score; TIDE, Tumor Immune Dysfunction and Exclusion.

Risk score correlation with invasion, EMT, and angiogenesis-related genes

Our analysis revealed that only angiogenesis and EMT scores differed between the two groups (Figure S4A-S4D). Furthermore, these scores showed a clear correlation with risk scores (Figure 8A-8C).

Figure 8 Correlation analysis of risk scores and metastasis/invasion indicators. The Lollipop plot depicts the correlation between risk scores and (A) angiogenesis-related genes, (B) epithelial-to-mesenchymal transition-related genes, and (C) invasion-related genes.

Characteristic gene regulation via methylation-related regulators and 167 non-coding RNAs

Most methylation-related regulators were correlated with the characteristic genes, with the exceptions of EIF3A, IGF2BP2, IGF2BP3, and NSUN2 (Figure 9A-9D). In both databases, a total of 37 miRNAs were predicted to target the two characteristic genes (EPO and UCP2) (Table S1). Subsequently, building on these 37 miRNAs, an additional 130 lncRNAs were identified as potential regulators (online supplementary table: https://cdn.amegroups.cn/static/public/tcr-2025-1146-1.xlsx). The constructed lncRNA-miRNA-mRNA network (Figure 9E) revealed that hsa-miR-125b-5p specifically regulated EPO, while the lncRNAs HOTTIP, AC022966.1, MIRLET7BHG, SNHG32, and AC092279.1 were identified as key regulators of hsa-miR-125b-5p.

Figure 9 Methylation correlation and regulatory network targeting two characteristic genes. Correlation heatmap illustrates the association between characteristic genes and methylation regulators of (A) m⁶A, (B) m¹A, (C) m⁵C, and (D) m⁷G. (E) LncRNA-miRNA-mRNA regulatory network diagram targeting two characteristic genes. Green represents lncRNA, orange represents miRNA, and red represents characteristic genes. *, P<0.05; **, P<0.01; ***, P<0.001.

Discussion

CC remains a prevalent malignancy affecting women worldwide, particularly in developing countries where efforts to reduce its incidence have been less effective. For patients diagnosed in advanced stages, prognosis is often poor, highlighting the need for early identification of prognostic factors for CC. There is increasing evidence now to suggest that efferocytosis could play a pivotal role in both cancer therapy and in predicting the prognosis of cancer patients. The link between efferocytosis and ACs in HPV-positive CC has been demonstrated in several studies. For instance, Gaiffe et al. found evidence of apoptotic human CC-derived cell engulfment by human primary fibroblasts (HPFs) (22). However, the precise relationship between efferocytosis and the prognosis of CC remains largely unexplored.

In this study, we systematically explored how CC and efferocytosis are connected and established a prognostic risk model based on two ERGs in CC patients, namely, EPO and UCP2. To understand the role of cell death in the progression of CC, the pathways and MFs associated with DEGs involved in cell death processes were analyzed in this study.

Our analysis revealed significant enrichment in key pathways related to fatty acid metabolism, Notch signaling, angiogenesis, and NK cell activity. Notably, elevated levels of circulating saturated fatty acids, such as palmitate, were identified as disruptors of normal cell death processes (23). Fatty acid metabolism was found to be intricately linked to the malignant progression of CC, emerging as a novel therapeutic target for the disease (24). Kim et al. identified that the Phosphatidylserine receptor brain-specific angiogenesis inhibitor 1 (BAI1)-Rac1 signaling pathway, triggered by Notch1 signaling, promoted efferocytosis by cancer-associated fibroblasts (25). Activation of Notch signaling has also been found to contribute to squamous differentiation in CC (26). Another crucial finding by Di Carlo et al. in cancer mouse models was that macrophage efferocytosis and polarization could drive pathological angiogenesis and immunosuppression by slow-cycling ADAM12 + PDGFRα + mesenchymal stromal cells, located at the tumor margins (27). Therefore, cell death in CC appears to be closely related to metabolic and immune pathways, both of which significantly influence the development and progression of the disease.

Erythropoietin (EPO) is a cytokine secreted by juxtaglomerular cells in the kidney. It regulates the survival, proliferation, and differentiation of erythroid progenitor cells by activating its transmembrane receptor, EPOR. It has been found to promote the M2 phenotype macrophages to ameliorate apoptosis and efferocytosis (28). Recent studies have revealed the oncogenic potential of EPO and EPOR, showing that they can promote the onset and progression of various tumors through anti-apoptotic and angiogenic effects, ultimately influencing tumor prognosis. For instance, Vukelic et al. observed a significant positive expression of EPO and EPOR in patients with laryngeal squamous cell carcinoma, with a notable correlation between EPO expression and patient survival outcomes (29). In the context of CC, research has indicated a significant elevation in serum EPO levels among patients with early-stage disease. This elevation has been linked to reduced sensitivity of CC cells to cisplatin and has been suggested as a marker for assessing disease progression (30).

UCP2 is a transporter in the mitochondrial inner membrane that reduces mitochondrial membrane potential, thereby affecting the phagocytic capacity during efferocytosis. Besides, UCP2 has been shown to promote the growth of malignant tumors by enhancing glycolysis and altering cellular signaling pathways (31), making it a new target for cancer treatment and a key predictor of prognosis. Research has highlighted that UCP2 expression may predict the efficacy of platinum-based chemotherapy in ovarian serous carcinoma (32). Results from a study in China have highlighted the elevated expression of UCP2 in ovarian cancer cells; specifically, patients with low UCP2 expression demonstrated improved OS outcomes compared to those with high UCP2 expression (33).

The expression of UCP2 is notably elevated in CC compared to high-grade cervical intraepithelial neoplasia. Functional studies have shown that UCP2 knockdown in SiHa and HeLa cells halts apoptosis during the G2 phase of the cell cycle. This knockdown also leads to the downregulation of protein levels of extracellular signal-regulated kinases (ERKs), as well as the mRNA levels of Ras, matrix metalloproteinase 2 (MMP-2), and matrix metalloproteinase 9 (MMP-9). This phenomenon results in the deceleration or attenuation of the proliferative, migratory, and invasive capabilities of SiHa and HeLa cells (34).

Moreover, UCP2 assumes a pivotal role in the onset and progression of CC and serves as a predictor of chemotherapy outcomes. Imai et al. explored the association between UCP2 expression and the efficacy of neoadjuvant chemotherapy (NAC) in individuals with locally advanced CC. Their findings highlighted that low UCP2 expression was associated with greater sensitivity to NAC. This observation was further verified through cell experiments, suggesting that UCP2 expression levels may be a predictive marker of the efficacy of NAC in patients with locally advanced CC, thereby helping to guide treatment decisions and improving patient prognosis (35).

Given the current emphasis on immunotherapy, targeting apoptosis within the TME offers a promising approach for cancer treatment. Specifically, strategies focusing on efferocytosis present a potential avenue for novel therapeutic interventions in addressing tumorigenesis and cancer management. Combining traditional therapies with efferocytosis-targeted approaches or other forms of immunotherapy holds the potential to enhance treatment efficacy and ultimately improve patient outcomes (36).

The influence of the TME on the prognosis of CC is substantial. In addition to tumor cells, the TME comprises tumor-associated fibroblasts, endothelial cells, and infiltrating immune cells, forming intricate regulatory networks that influence tumor behavior and patient outcomes (37-39). Our findings in this study revealed a negative correlation between the risk score of the high-risk group and the stromal score, as well as a positive correlation with tumor purity, indicating potentially lower stromal cell infiltration in the high-risk group. Furthermore, immune infiltration analysis revealed significant differences between the high-risk and low-risk groups in initial B cells, resting NK cells, and mast cells. These variations in immune cell infiltration are crucial, as they are critically involved in tumorigenesis and tumor progression.

B cells exhibit a dual role in tumor immunity, acting as both promoters of antitumor immune responses and regulators that help control excessive inflammation (40). In an HPV-mediated CC mouse model, B cells were found to promote tumor progression (41). NK cells play a core role in combating virus-infected cells and tumors. NK cells can recognize and destroy HPV-infected cells, triggering an inflammatory response. However, HPV possesses mechanisms to evade both innate and adaptive immune responses, reducing the efficacy of NK cells in eliminating infected cells. Notably, a significant focus of CC immunotherapy research involves harnessing the potential of NK cells to improve tumor clearance (42).

Mast cells, on the other hand, play a crucial role in allergic reactions, pathogen immune responses during infection, angiogenesis, and innate and adaptive immune regulation (43). Derived from hematopoietic stem cells in the bone marrow, mast cells are characterized by basophilic granules containing a variety of active mediators and are widely distributed in mammals. Upon activation, mast cells produce a multitude of mediators and cytokines, participating in various physiological and pathological processes such as reproduction, immunity, allergic reactions, and cancer. The accumulation of mast cells within tumors is often associated with poor prognosis in CC (44), which is closely related to their critical role in promoting tumor angiogenesis and progression. In CC, activated mast cells release potent pro-angiogenic factors, such as VEGF, bFGF, TGF-β, TNF-α, and IL-8, which regulate the growth and migration of vascular endothelial cells. Proteases and heparin stimulate stromal cells to secrete matrix metalloproteinases and regulate their expression and function, degrading the extracellular matrix to create space for new blood vessels. Additionally, histamine and lipid-derived mediators can modulate the function and structure of nascent vessels, increasing their permeability and thereby facilitating tumor infiltration and metastasis. Our research results also highlight the strong correlation between key CC-related genes and mast cell presence. Mast cells may contribute to the invasiveness and metastatic potential of CC cells through the release of histamine and cannabinoids, facilitating tumor spread (45). In addition, mast cells are also known to be actively involved in the regulation of CC angiogenesis.

Recent evidence suggests that UCP2 could regulate degranulation and inhibit activation of mast cells (46,47). Similarly, EPO has been shown to act as an anti-inflammatory signal in murine mast cells (48). In the current study, we assessed the predictive potential of EPO and UCP2 and established an effective risk model to predict the clinical outcomes of patients with CC. We observed significant differences in prognosis between high-risk and low-risk groups, specifically in terms of T1/T2, N0, M0, age <60 years, and tumor stage I/II, underscoring the model’s ability to distinguish between patients with differing outcomes. Importantly, the inclusion of risk score and lymphatic vessel invasion in the model strengthened its predictive accuracy. Lymphatic vessel invasion is a well-established prognostic marker in CC, frequently associated with poor prognosis (49,50). For instance, Xu et al. reported that large tumor volume, lymph node metastasis, distant metastasis, and parauterine infiltration are usually linked to an unfavorable prognosis in patients with CC (51).

The infiltration of vascular and lymphatic vessels, along with diverse cell types within the tumor stroma, in the solid component of tumors is associated with the clinicopathological characteristics of early-stage CC. Given these associations, we speculated that both model genes—EPO and UCP2—may be involved in the inflammatory response mediated by mast cells during the progression of CC. To further elucidate the mechanism of action of EPO and UCP2 in CC, combined with the complex regulatory function of miRNAs in macrophages, we analyzed the miRNA regulatory mechanisms governing these two genes in this study.

miRNAs are non-coding RNAs that regulate gene expression at the post-transcriptional level, playing a significant role in tumorigenesis, development, and progression through various epigenetic modifications, including DNA methylation and histone modification. In this study, we found that hsa-miR-125b-5p specifically regulates EPO. Previous studies, such as that by Notley et al., showed that miR-125b promotes the progression of CC by regulating HMGA1. However, the degree of DNA methylation in ACs can determine their plasticity in immune regulation (52). Notably, the demethylation of AC DNA negates their protective effects, while re-methylation reverses the effects of activation and restores their ability to inhibit inflammation. These findings suggest that DNA methylation in ACs acts as a molecular switch that determines their immune function.

Efferocytosis is a crucial process for clearing dead cells and preventing inflammation. In tumors, it regulates tumor progression and prognosis by influencing the function of immune cells within the TME. Previous research has indicated that the methylation level of ACs is a key regulator of their immune modulatory function. Based on this, we hypothesized that the ERGs EPO and UCP2 might be associated with methylation mechanisms. To test this hypothesis, we analyzed the association between EPO/UCP2 and RNA methylation regulators and found significant correlations with the majority of these factors. This suggests that in CC, EPO and UCP2 may regulate key processes such as immunity, inflammation, and apoptosis at a deeper level by influencing the methylation status of miRNAs, ultimately affecting patient outcomes.

However, there are also some limitations in this study. The primary limitation of this study lies in its small sample size, particularly in the PCR validation cohort (n=10), which limited the statistical power for robust prognostic analysis. Although survival analysis using the survival R package showed no significant prognostic differences for EPO and UCP2, this is likely attributable to the small cohort size and limited follow-up data. As such, the current validation primarily confirms the detectability and relative expression of the characteristic genes rather than their prognostic value. The clinical relevance and predictive robustness of the risk score model require further verification in a larger, multi-center, and prospectively enrolled cohort with long-term follow-up to ensure generalizability. Moreover, the retrospective design and reliance on public databases may introduce selection bias, warranting cautious interpretation. Future studies should also incorporate additional clinical variables to comprehensively assess the prognostic significance of efferocytosis-related genes. As a result, extensive prospective studies, as well as in vitro and in vivo experimental investigations, are necessary to validate our findings. Particularly, further population-based validation is crucial for gaining deeper insights into the relationship between efferocytosis-related genes and CC. Despite these limitations, this study retains a high degree of reliability. The exploration of risk models based on efferocytosis-related genes in CC will remain a primary focus of future research endeavors.

Conclusions

In conclusion, a robust prognostic risk model based on two characteristic genes associated with CC was successfully established in this research, and their molecular mechanisms involved in the progression of this disease were further explored. The model proved effective in distinguishing prognostic differences among patients with CC, optimizing the efficacy of current prognostic indicators, and allowing for effective risk stratification of patients. While our preliminary validation cohort included patients with HPV16-positive, early-stage (IA–IIA) disease, the model’s generalizability to this specific subgroup requires future prospective validation. If successfully validated, it holds promise for guiding post-operative risk stratification and personalizing adjuvant therapy decisions in this clinical context, meriting further investigation.

Acknowledgments

We thank the TCGA, GEO, MSigDB and other databases for providing invaluable data for statistical analyses.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1146/rc

Data Sharing Statement: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1146/dss

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1146/prf

Funding: This study was funded by the Natural Science Foundation of Shanxi Province (No. 202203021221273), Shanxi Province Science Foundation for Youths (No. 2019047), Shanxi Province Key National Science and Technology Cooperation Projects (No. 202304041101013), and Shanxi Province Youth Science Research Project (No. 202403021212222).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1146/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by the ethics committee of the Second Hospital of Shanxi Medical University of Shanxi Province, China (No. 2022-YX-226) and informed consent was obtained from all individual participants.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Small W Jr, Bacon MA, Bajaj A, et al. Cervical cancer: A global health crisis. Cancer 2017;123:2404-12. [Crossref] [PubMed]
2. Pal A, Kundu R. Human Papillomavirus E6 and E7: The Cervical Cancer Hallmarks and Targets for Therapy. Front Microbiol 2019;10:3116. [Crossref] [PubMed]
3. Gopu P, Antony F, Cyriac S, et al. Updates on systemic therapy for cervical cancer. Indian J Med Res 2021;154:293-302. [Crossref] [PubMed]
4. Kusakabe M, Taguchi A, Sone K, et al. Carcinogenesis and management of human papillomavirus-associated cervical cancer. Int J Clin Oncol 2023;28:965-74. [Crossref] [PubMed]
5. Vu M, Yu J, Awolude OA, et al. Cervical cancer worldwide. Curr Probl Cancer 2018;42:457-65. [Crossref] [PubMed]
6. Cibula D, Dostálek L, Jarkovsky J, et al. Post-recurrence survival in patients with cervical cancer. Gynecol Oncol 2022;164:362-9. [Crossref] [PubMed]
7. Baezconde-Garbanati L, Agurto I, Gravitt PE, et al. Barriers and innovative interventions for early detection of cervical cancer. Salud Publica Mex 2019;61:456-60. [Crossref] [PubMed]
8. Doran AC, Yurdagul A Jr, Tabas I. Efferocytosis in health and disease. Nat Rev Immunol 2020;20:254-67. [Crossref] [PubMed]
9. Ampomah PB, Cai B, Sukka SR, et al. Macrophages use apoptotic cell-derived methionine and DNMT3A during efferocytosis to promote tissue resolution. Nat Metab 2022;4:444-57. [Crossref] [PubMed]
10. Zhao T, Su Z, Li Y, et al. Chitinase-3 like-protein-1 function and its role in diseases. Signal Transduct Target Ther 2020;5:201. [Crossref] [PubMed]
11. Hermetet F, Jacquin E, Launay S, et al. Efferocytosis of apoptotic human papillomavirus-positive cervical cancer cells by human primary fibroblasts. Biol Cell 2016;108:189-204. [Crossref] [PubMed]
12. Elder AM, Stoller AR, Black SA, et al. Macphatics and PoEMs in Postpartum Mammary Development and Tumor Progression. J Mammary Gland Biol Neoplasia 2020;25:103-13. [Crossref] [PubMed]
13. Mehrotra P, Ravichandran KS. Drugging the efferocytosis process: concepts and opportunities. Nat Rev Drug Discov 2022;21:601-20. [Crossref] [PubMed]
14. Trzeciak A, Wang YT, Perry JSA. First we eat, then we do everything else: The dynamic metabolic regulation of efferocytosis. Cell Metab 2021;33:2126-41. [Crossref] [PubMed]
15. Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015;43:e47. [Crossref] [PubMed]
16. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 2013;14:7. [Crossref] [PubMed]
17. Friedman J, Hastie T, Tibshirani R. Regularization Paths for Generalized Linear Models via Coordinate Descent. J Stat Softw 2010;33:1-22.
18. Hebert PD, Cywinska A, Ball SL, et al. Biological identifications through DNA barcodes. Proc Biol Sci 2003;270:313-21. [Crossref] [PubMed]
19. Sturm G, Finotello F, Petitprez F, et al. Comprehensive evaluation of transcriptome-based cell-type quantification methods for immuno-oncology. Bioinformatics 2019;35:i436-45. [Crossref] [PubMed]
20. Malta TM, Sokolov A, Gentles AJ, et al. Machine Learning Identifies Stemness Features Associated with Oncogenic Dedifferentiation. Cell 2018;173:338-354.e15. [Crossref] [PubMed]
21. Doosti Z, Ebrahimi SO, Ghahfarokhi MS, et al. Synergistic effects of miR-143 with miR-99a inhibited cell proliferation and induced apoptosis in breast cancer. Biotechnol Appl Biochem 2024;71:993-1004. [Crossref] [PubMed]
22. Gaiffe E, Prétet JL, Launay S, et al. Apoptotic HPV positive cancer cells exhibit transforming properties. PLoS One 2012;7:e36766. [Crossref] [PubMed]
23. Li S, Sun Y, Liang CP, et al. Defective phagocytosis of apoptotic cells by macrophages in atherosclerotic lesions of ob/ob mice and reversal by a fish oil diet. Circ Res 2009;105:1072-82. [Crossref] [PubMed]
24. Ping P, Li J, Lei H, et al. Fatty acid metabolism: A new therapeutic target for cervical cancer. Front Oncol 2023;13:1111778. [Crossref] [PubMed]
25. Kim HJ, Yang K, Kim K, et al. Reprogramming of cancer-associated fibroblasts by apoptotic cancer cells inhibits lung metastasis via Notch1-WISP-1 signaling. Cell Mol Immunol 2022;19:1373-91. [Crossref] [PubMed]
26. Li W, Zhang X, Gao T, et al. Jagged1 contained in MSC-derived small extracellular vesicles promotes squamous differentiation of cervical cancer by activating NOTCH pathway. J Cancer Res Clin Oncol 2023;149:18093-102. [Crossref] [PubMed]
27. Di Carlo SE, Raffenne J, Varet H, et al. Depletion of slow-cycling PDGFRα(+)ADAM12(+) mesenchymal cells promotes antitumor immunity by restricting macrophage efferocytosis. Nat Immunol 2023;24:1867-78. [Crossref] [PubMed]
28. Govindappa PK, Elfar JC. Erythropoietin promotes M2 macrophage phagocytosis of Schwann cells in peripheral nerve injury. Cell Death Dis 2022;13:245. [Crossref] [PubMed]
29. Vukelic J, Dobrila-Dintinjana R, Marijic B, et al. New insights into erythropoietin and erythropoietin receptor in laryngeal cancer tissue. Medicine (Baltimore) 2021;100:e23943. [Crossref] [PubMed]
30. Vázquez-Mellado MJ, Cortés-Ballinas LG, Blanco-Flores IC, et al. Erythropoietin promotes expression of survivin via STAT3 activation and reduces sensitivity to cisplatin in cervical cancer cells. Oncol Rep 2019;41:1333-41. [Crossref] [PubMed]
31. Lee S, Moon H, Kim G, et al. Anion Transport or Nucleotide Binding by Ucp2 Is Indispensable for Ucp2-Mediated Efferocytosis. Mol Cells 2015;38:657-62. [Crossref] [PubMed]
32. Kawanishi M, Fukuda T, Shimomura M, et al. Expression of UCP2 is associated with sensitivity to platinum-based chemotherapy for ovarian serous carcinoma. Oncol Lett 2018;15:9923-8. [Crossref] [PubMed]
33. Chen Z, Ho M, Chau PH. Prevalence, Incidence, and Associated Factors of Possible Sarcopenia in Community-Dwelling Chinese Older Adults: A Population-Based Longitudinal Study. Front Med (Lausanne) 2021;8:769708. [Crossref] [PubMed]
34. Zhu S, Lin W, Chen S, et al. The correlation of muscle thickness and pennation angle assessed by ultrasound with sarcopenia in elderly Chinese community dwellers. Clin Interv Aging 2019;14:987-96. [Crossref] [PubMed]
35. Imai K, Fukuda T, Wada T, et al. UCP2 expression may represent a predictive marker of neoadjuvant chemotherapy effectiveness for locally advanced uterine cervical cancer. Oncol Lett 2017;14:951-7. [Crossref] [PubMed]
36. Gheibi Hayat SM, Bianconi V, Pirro M, et al. Efferocytosis: molecular mechanisms and pathophysiological perspectives. Immunol Cell Biol 2019;97:124-33. [Crossref] [PubMed]
37. Ding J, Wang C, Sun Y, et al. Identification of an Autophagy-Related Signature for Prognosis and Immunotherapy Response Prediction in Ovarian Cancer. Biomolecules 2023;13:339. [Crossref] [PubMed]
38. Chang J, Wu H, Wu J, et al. Constructing a novel mitochondrial-related gene signature for evaluating the tumor immune microenvironment and predicting survival in stomach adenocarcinoma. J Transl Med 2023;21:191. [Crossref] [PubMed]
39. Yang J, Wang YP, Yu WG, et al. Comparative evaluation of antitumor Chinese patent medicine injections based on evidence-based medicine and Markov model. World J Tradit Chin Med. 2024;10:220-9.
40. Downs-Canner SM, Meier J, Vincent BG, et al. B Cell Function in the Tumor Microenvironment. Annu Rev Immunol 2022;40:169-93. [Crossref] [PubMed]
41. Tang A, Dadaglio G, Oberkampf M, et al. B cells promote tumor progression in a mouse model of HPV-mediated cervical cancer. Int J Cancer 2016;139:1358-71. [Crossref] [PubMed]
42. Gutiérrez-Hoya A, Soto-Cruz I. NK Cell Regulation in Cervical Cancer and Strategies for Immunotherapy. Cells 2021;10:3104. [Crossref] [PubMed]
43. Jiang L, Hua Y, Shen Q, et al. Role of mast cells in gynecological neoplasms. Front Biosci (Landmark Ed) 2013;18:773-81. [Crossref] [PubMed]
44. Wang J, Li Z, Gao A, et al. The prognostic landscape of tumor-infiltrating immune cells in cervical cancer. Biomed Pharmacother 2019;120:109444. [Crossref] [PubMed]
45. Rudolph MI, Boza Y, Yefi R, et al. The influence of mast cell mediators on migration of SW756 cervical carcinoma cells. J Pharmacol Sci 2008;106:208-18. [Crossref] [PubMed]
46. Sismanopoulos N, Delivanis DA, Alysandratos KD, et al. Mast cells in allergic and inflammatory diseases. Curr Pharm Des 2012;18:2261-77. [Crossref] [PubMed]
47. Tagen M, Elorza A, Kempuraj D, et al. Mitochondrial uncoupling protein 2 inhibits mast cell activation and reduces histamine content. J Immunol 2009;183:6313-9. [Crossref] [PubMed]
48. Wiedenmann T, Ehrhardt S, Cerny D, et al. Erythropoietin acts as an anti-inflammatory signal on murine mast cells. Mol Immunol 2015;65:68-76. [Crossref] [PubMed]
49. Senchukova MA, Makarova EV, Shurygina EI, et al. Morphological Characteristics and Clinical Significance of Different Types of Tumor Vessels in Patients with Stages I-IIA of Squamous Cervical Cancer. J Oncol 2020;2020:3818051. [Crossref] [PubMed]
50. Yan W, Qiu S, Ding Y, et al. Prognostic value of lymphovascular space invasion in patients with early stage cervical cancer in Jilin, China: A retrospective study. Medicine (Baltimore) 2019;98:e17301. [Crossref] [PubMed]
51. Xu C, Ma T, Sun H, et al. Markers of Prognosis for Early Stage Cervical Cancer Patients (Stage IB1, IB2) Undergoing Surgical Treatment. Front Oncol 2021;11:659313. [Crossref] [PubMed]
52. Notley CA, Jordan CK, McGovern JL, et al. DNA methylation governs the dynamic regulation of inflammation by apoptotic cells during efferocytosis. Sci Rep 2017;7:42204. [Crossref] [PubMed]