A Cellular Senescence-Related Signature Predicts Cervical Cancer Patient Outcome and Immunotherapy Sensitivity

Cervical cancer (CC) is one of the most prevalent gynecological malignancies. The rate of mortality and morbidity among patients with CC is high. Cellular senescence is involved in tumorigenesis as well as in the cancer progression. However, the involvement of cellular senescence in CC development is still unclear and requires further investigation. In this study, we retrieved data on cellular senescence-related genes (CSRGs) from the "CellAge" Database. We used the TCGA-CESC and CGCI-HTMCP-CC datasets as the training and validation sets, respectively. Finally, a signature was constructed using "univariate" and "Least Absolute Shrinkage and Selection Operator" (LASSO) Cox regression analysis, which contains eight CSRGs. Using this signature, we calculated the risk scores of all patients in the training and validation cohorts and categorized them into the low-risk group (LR-G) and the high-risk group (HR-G). Results showed that, compared to patients in the HR-G, those in the LR-G demonstrated a more positive clinical prognosis, more abundant immune cell infiltrations, and a more active immune response. The signature could also modulate the expression of SASP factors. In vitro studies showed an increased expression of SERPINE1 and IL-1α genes included in the signature in CC cells and tissues. Our findings help to deepen our insights into the etiology of CC, which could be beneficial for prognostic prediction and immunotherapy in clinical practice. Supplementary Information The online version contains supplementary material available at 10.1007/s43032-023-01305-w.


Introduction
Globally, cervical cancer (CC) is the second most common gynecological malignancy.In 2020, nearly 604,000 new incidences and 342,000 mortalities due to CC were reported [1].According to previous reports, the prognosis of CC patients primarily depended on the stage at the time of diagnosis.The survival rate of the early-stage patients is greater than 90%; however, the survival rate of patients with advanced and metastatic CC is less than 20% [2].Additionally, the therapeutic strategy also depends on the stage of the disease.For the early-stage or localized CC, the current treatment options include surgical resection, concurrent chemotherapy, and radiotherapy.For patients with metastatic CC, systemic therapies are widely used [3][4][5]; however, their therapeutic efficacy is moderate.Immunotherapy is a rapidly evolving field, that involves modifying and recruiting the immune system to target cancer cells more efficiently and accurately.Currently, immunotherapy has been used in the treatment of various tumors.Common immunotherapy methods mainly include cytokine treatment, immune cell Huijing Shao and Xia Li contributed equally to this work.therapy, tumor vaccines, and immune checkpoint inhibitors [6][7][8].Research has shown that patients with CC could benefit from immunotherapy.Pembrolizumab (an anti-PD-1) has been approved as an effective immunotherapy drug for treating patients with CC [9].However, the efficacy of immunotherapy differs from patient to patient and is largely dependent on their immune status.Given the serious threat of cervical cancer to women's health, attention should be paid to identify effective biomarkers that could facilitate early diagnosis, predict prognosis, and improve immunotherapeutic efficacy.
Cellular senescence is characterized by a stable and terminal cell cycle arrest, accompanied by complex changes at the molecular level.Research has demonstrated that cellular senescence has both a positive and a negative impact on the development and advancement of cancer.Senescent cells undergo persistent cell cycle arrest, which helps to maintain tissue homeostasis, prevent tumorigenesis by arresting cell proliferation and enhance immune surveillance of cancer cells [10][11][12].However, the failure to eliminate senescent cells by the immune system and their accumulation in the tumor microenvironment (TME) could lead to adverse effects.Senescence-associated secretory phenotype (SASP) is a characteristic phenotype of cellular senescence that includes multiple components, such as chemokines, cytokines, growth factors, interleukins, and proteases [13].These SASP components could remodel the tumor microenvironment (TME) and promote tumorigenesis.Furthermore, the advent of high-throughput technology and several publicly available databases, such as "SeneQuest" [14], "Cel-lAge" [15] and "Human Cellular Senescence Gene Database (HCSGD)" [16] could aid in identifying cellular senescenceassociated gene signatures in cancer.Using senescenceassociated genes, Luo et al. established a new model for predicting the survival outcomes of breast cancer patients [17].However, the contribution of cellular senescence to the development of CC is still unknown and no risk models based on cellular senescence have been established for CC.
Therefore, in this study, we screened for cellular senescence-related biomarkers.In addition, we developed a risk model to predict the prognosis, the tumour immune microenvironment (TIME), and the response of patients to immunotherapy in CC.This risk model could aid clinicians in designing more effective therapeutic strategies.

Data Retrieval and Processing
The transcriptomic data of patients with CC was obtained from the "TCGA-TARGET-GTEx" database using the UCSC Xena (http:// xena.ucsc.edu/) browser.This database was used for differential expression analysis.Further, the data on gene expression, as well as clinical and somatic mutations in patients with CC, were retrieved from the "TCGA-CESC" (https:// portal.gdc.cancer.gov/) dataset, which served as the training set.Additionally, the data on RNA sequencing and survival of patients with CC were retrieved from the "CGCI-HTMCP-CC" (https:// portal.gdc.cancer.gov/) dataset, which served as the validation set.Samples without vital clinicopathological information or survival data were excluded from the study.Finally, we included 286 patients from the "TCGA-CESC" dataset and 117 patients from the CGCI-HTMCP-CC dataset.We identified 279 genes from the "CellAge" (https:// genom ics.senes cence.info/ cells/) database as cellular senescence-related genes (CSRGs).Survival analysis of a single gene was performed using the GEPIA database (http:// gepia.cancer-pku.cn/ index.html).

Identification and Functional Enrichment of DEGs
Utilizing the "DESeq2" R package, we identified differentially expressed genes (DEGs) with a threshold criterion of "P < 0.05" and "|log2-fold change|≥ 1".Next, we employed the "ClusterProfiler" R package [18] for performing "Gene Ontology (GO)" and "Kyoto Encyclopedia of Genes and Genomes (KEGG)" pathway enrichment analyses to determine the pathways that are significantly enriched by the DEGs (an adjusted P < 0.05).Finally, we constructed a volcano plot in the "ggplot2" R package [19] for visualizing DEGs and the results of GO and KEGG analyses.

Establishing and Validating the Prognostic Signature associated with Cellular Senescence
We employed the "VennDiagram" (version 1.7.3)R package to construct a Venn diagram to identify the overlapping genes between the DEGs and CSRGs.First, a univariate Cox regression analysis was used to determine the correlation between the 23 overlapping genes and the survival outcomes of patients.Then, the "glmnet" R package [20] was employed to conduct the "Least Absolute Shrinkage and Selection Operator (LASSO)" Cox regression analysis on significant genes (P < 0.05).Subsequently, we used tenfold cross-validation in a LASSO Cox regression analysis for constructing an optimal risk model.Finally, we used the risk model to calculate the risk scores based on the following formula: The risk scores of all patients in the training set were calculated using this formula.Next, we used the median value as the parameter to categorize all patients into the high-risk group (HR-G) and the low-risk group (LR-G). (1) Principal Component Analysis (PCA) was performed using the "prcomp" R package for subsequent clustering.Additionally, we performed Kaplan-Meier (KM) analysis using the "survival" and "survminer" R packages [21] to determine differences in the survival outcomes of the patients in the two groups.Finally, the receiver operating characteristic curve (ROC) curves were generated using the "timeROC" R package to assess the predictive performance of the risk model.Next, we verified the predictive efficiency of the risk model in patients in the validation set using the aforementioned formula.

Clinical Significance of the Signature and a Novel Prognostic Nomogram
We employed the "Univariate and Multivariate Cox Regression Analyses" for determining the independent clinical prognostic significance of the risk model.Next, the "rms" (Version 6.3.0)R package was employed for constructing a predictive nomogram using the results of univariate and multivariate Cox analysis.First, a nomogram was constructed for predicting the 1-, 3-, and 5-year overall survival (OS) rates of patients with CC.Subsequently, the calibration curves were plotted and Decision Curve Analysis (DCA) was conducted to study the validity of the nomogram in the clinical setting.

Gene Set Enrichment Analysis
We identified DEGs in patients in the two risk groups using the "DESeq2" R package based on the following threshold criteria: "P < 0.05" and "|log2-fold change|≥ 1".This was used to elucidate the mechanisms underlying the pathogenesis of CC.Subsequently, "Gene Set Enrichment Analysis" (GSEA) was performed on the DEGs."P < 0.05", "false discovery rate q-value < 0.25", and "normalized enrichment score > 1.5" were considered to be significantly enriched genes.

Evaluation of Immune Cell Infiltration and Tumor-Associated Gene Set Scores
We used the three most commonly used immune analysis methods for assessing the infiltration of immune cells in the TME of CC patients; (1) To determine differences in the TME of patients in both groups, we assessed the levels of stromal and immune cell infiltration using the "Estimation of Stromal and Immune Cells in Malignant Tumor Tissues using Expression Data" (ESTIMATE) R package.(2) Next, we determined the proportion of 22 immune cell types using the "Cell-Type Identification by Estimating Relative Subsets of RNA Transcripts" (CIBERSORT) R package.(3) Finally, we used "Single-Sample Gene Set Enrichment Analysis" (ssGSEA) for evaluating the percentage of 29 immunocytes and determining the enrichment scores of angiogenesis, epithelial-mesenchymal transition (EMT), and hypoxia-related genes.Angiogenesis-, EMT-, and hypoxia-related gene sets were retrieved from the Molecular Signatures Database.Table S1 shows the relevant marker genes.

Prediction of Response to Immunotherapy
We utilized the Wilcoxon test to compare the immune checkpoint gene expression among patients in the two groups.
The "IMvigor210" (http:// resea rch-pub.gene.com/ IMvig or210 CoreB iolog ies/) dataset contains patients with metastatic urothelial cancer who were treated with an anti-PD-L1 agent.We used the "IMvigor210 Core Biologies" R package [22] on this data to determine the correlation between immunotherapeutic efficacy and the risk model.We used the median risk score to categorize these patients into two risk groups.Finally, the immunotherapeutic efficacy and clinical outcomes of patients in these two groups were compared.

Analysis of Mutations and Drug Sensitivity
The "maftools" R package was used to convert the original files from TCGA to the mutation annotation format, for comparing the mutational landscapes in patients in the two groups [23].Finally, to examine the clinical response of patients to chemotherapy, the "pRRophetic" R package was used for calculating the half-maximal inhibitory concentration (IC50) values of 138 commonly used chemotherapeutic agents [24].

Sample Collection
We collected ten pairs of matched CC and normal tissues from the Putuo Hospital Affiliated to Shanghai University of Traditional Chinese Medicine between May and December 2022 and stored them at -80 °C.Next, we performed a real-time polymerase chain reaction (qRT-PCR) on these samples.This study was approved by the Ethics Committee of Putuo Hospital Affiliated to Shanghai University of Traditional Chinese Medicine.Written informed consent was obtained from all participants.

Cell Culture
HcerEpic, CaSki, HeLa, and SiHa cells were purchased from the Shanghai Institute of Cell Biology (Shanghai, China).HcerEpic and CaSki cells were cultured in RPMI-1640 medium.HeLa and SiHa were cultured in DMEM (Gibco) and α-MEM (Gibco), respectively.All media were supplemented with 10% FBS (Gibco), penicillin, and streptomycin.The cells were incubated at 5% CO2 and 37 °C.

Western Blotting
RIPA buffer (Epizyme, Shanghai, China) was used for cell lysis to extract total protein, which was then separated using gel electrophoresis at 120 V and transferred to PVDF membranes.After being blocked with 5% milk, the PVDF membranes were incubated with primary antibodies for an overnight period at 4 °C, followed by secondary antibodies for an additional two hours at room temperature.A ECL system was used to visualize the protein bands (NCM Biotech, Suzhou, China).The protein bands were visualised using an ECL system (NCM Biotech, Suzhou, China).The primary antibodies used for Western blotting were an anti-SERPINE1 (66261-1-Ig, Proteintech, IL, USA) and an anti-IL-1α (ab300501, Abcam) antibody.

Statistical Analysis
We performed all statistical analyses with the help of the R (version 4.1.3)software and Perl.We compared the patients' OS rates in both groups using KM analysis.The infiltration of immunocytes and the expression of immune checkpoint genes in patients from the two groups were compared using the Wilcoxon test, with a significance level of p < 0.05.

Differentially Expressed CSRGs in Patients with CC
Figure 1 shows the study design.We identified a total of 1,691 DEGs in the tissues of 306 patients with CC and 13 adjacent normal samples from the TCGA-TARGET-GTEX dataset and visualized using a volcano plot.Of these DEGs, 695 were downregulated and 996 were upregulated genes (Fig. 2a).GO and KEGG analysis revealed that the DEGs were associated with immune and inflammatory responses (Fig. S1).To identify cellular senescence-related DEGs, we intersected DEGs with 279 CSRGs retrieved from the "Cel-lAge" database.Finally, we identified 23 overlapping genes for further investigation (Fig. 2b).

Establishment of a CSRGs Prognosis Signature for CC
Univariate Cox regression analysis was performed on 23 overlapping genes to identify genes related to prognosis in CC.We performed the LASSO Cox regression analysis on eight genes with a p-value of less than 0.05 (Fig. 2c and Table S1).We then identified eight significant prognosisrelated genes using LASSO Cox regression analysis with tenfold cross-validation based on the optimal lambda value (Fig. 2d and e).We constructed the risk model based on the expressions and coefficients of these eight genes.A correlation was observed among the expression levels of these eight prognostic genes.We can observe correlations among the eight prognostic genes, such as a positive correlation between SERPINE1 and IL-1α expression and a negative correlation between SERPINE1 and E2F1 expression (Fig. 2f).

Assessing the Risk Model in Patients from the Training and Validation Sets
The median risk score served as the basis for distinguishing the patients from the training set into the HR-G and the LR-G.PCA revealed that the risk scores could efficiently distinguish patients in the HR-G from those in the LR-G (Fig. 3a). Figure 3b shows a risk plot ranking and categorizing patients into the HR-G and the LR-G.Interestingly, the number of deceased patients was higher in the HR-G (Fig. 3c).Furthermore, the heatmap in Fig. 3d displays the expression of the eight genes that constitute the risk model.High levels of AGT, CAVIN1, GNG11, IL-1α, and SERPINE1 were observed in patients in the HR-G.On the other hand, low expression of SOX2, E2F1, and IFNG were observed in patients in the LR-G.Additionally, KM survival curves showed that the prognosis of patients in the HR-G was poor compared to the LR-G (Fig. 3e).The "time-dependent ROC" analysis (Fig. 3f) showed that the Area Under the Curve (AUC) values for predicting the 1-year OS rate was 0.800, the 3-year was 0.6999, and the 5-year OS rate was 0.671.The results demonstrated the capability of the risk model in predicting the prognosis of CC patients.To validate the reliability of the risk model in estimating patients' prognoses, a risk score was calculated for each individual from the CGCI-HTMCP-CC dataset, and the entire cohort was also divided into the HR-G and LR-G.PCA showed significant differences between the two groups (Fig. 3g), further indicating that the risk signature could satisfactorily distinguish the prognoses of CC patients.Additionally, KM survival curves revealed that the survival rate of patients in the HR-G was poor (Fig. 3h).
The AUCs for predicting the 1-and 2-year OS rates were 0.579 and 0.588, respectively (Fig. 3i).Together, these results suggest that the risk model could distinguish between the prognoses of patients.The risk scores of the deceased (Fig. S2a) and elderly patients (Fig. S2b) were higher.However, no difference in the risk score was observed in patients with different tumor grades (Fig. S2c) or stages (Fig. S2d).Furthermore, univariate COX regression analysis revealed that the risk score of patients and the tumor stage were both reliable prognostic factors, which were subsequently included in the multivariate COX regression analysis (Fig. 4a).Multivariate Cox regression analysis showed that the risk score (P < 0.001) and tumor stage (P = 0.006) served as independent prognostic factors (Fig. 4b).Finally, we established a predictive nomogram by integrating the risk score and tumour stages of patients to determine the 1-, 3-, and 5-year OS rates of patients with CC (Fig. 4c).The calibration curves demonstrated that the actual and the nomogram-predicted OS rates were highly similar (Fig. 4d).Furthermore, we calculated the AUCs of the risk model, nomogram, and tumor stage for determining the 1-, 3-, and 5-year OS rates using the "timedependent ROC" analysis.Figure 4e-g revealed that the AUCs of the tumor stage, risk model, and nomogram were significant, thus indicating adequate predictive accuracy for these three factors.In addition, DCA revealed that the risk model, nomogram, and tumour stage demonstrated good performance in predicting the clinical net benefit (Fig. 4h-j).

Biological Processes Associated with the Risk Model and Correlation between the SASP and the Risk Model
The results of the survival analysis of CC patients in our study confirmed the predictive efficiency of the risk model, prompting us to explore the underlying mechanisms.We identified 1,058 DEGs in patients in the two groups in the training set (with "P < 0.05" and "|log2-fold change|≥ 1")."GO and KEGG pathway enrichment" analysis revealed that multiple immune and inflammatory pathways were altered (Fig. S3).Next, we compared the expression of different SASP factors in patients in the two groups.

Relationship Between the Signature and the Immune Landscape
SASP is characterized by its ability to induce inflammation.Given the pro-inflammatory nature of the SASP, it is likely that senescent cells can attract immune cells.Figure 6a shows a significant enrichment of multiple pathways related to immune responses, such as natural killer cell-mediated cytotoxicity, B and T cell pathways, and the Toll-like receptor signaling pathway, in patients in the LR-G.These results suggest that the risk model could aid in predicting the immune status of patients with CC.Subsequently, we used algorithms such as the "ESTIMATE", "CIBERSORT", and "ssGSEA" to assess the levels of immune cell infiltration in the TME.The "ESTIMATE" algorithm was used to perform immune analysis in patients with CC.In the LR-G, CC patients had higher immune scores and lower stromal scores than those in the HR-G (Fig. 6b).Similarly, Pearson correlation analysis showed a positive correlation between the patients' stromal and risk scores (r = 0.131, P = 0.027; Fig. 6c).Further, a negative correlation between the immune and risk scores was observed (r = -0.379,P = 0.001; Fig. 6d).Additionally, ssGSEA showed significant differences between the two groups in terms of their immune profiles.A significant increase in B cells, neutrophils, NK cells, CD8 + T cells, and Th1 and Th2 cell infiltration was observed in patients in the LR-G (Fig. 6e).Additionally, multiple immune function signatures were significantly activated in patients in the LR-G (Fig. 6f).Immune analysis performed using the "CIBERSORT" algorithm revealed an increase in the abundance of activated NK cells, naïve B cells, neutrophils, T follicular helpers, M1 macrophages, and CD8 + T cells, as well as activated mast cells in patients in the LR-G.However, a significant increase in the abundance of neutrophils, resting memory CD4 + T cells, and M0 macrophages was observed in patients in the HR-G (Fig. 6g).Together, these results showed an increased antitumor immune activity in low-risk patients.And, the risk model was an effective tool in assessing the immune status of CC patients with different degrees of risk.

Correlation Between the Risk Model and Immune Checkpoints, as well as its Potential to Predict Immunotherapy Response
Studies have shown the significant involvement of immune checkpoint genes in regulating the infiltration of immune cells.
Our results demonstrated that a correlation exists between the risk model and the tumour immune microenvironment of patients with CC.Therefore, we determined the expression of PD-1, CTLA-4, and PD-L1 in patients in the two groups to understand the complex correlation among the risk scores, immunocytes infiltration, and the expression of immune checkpoint genes.Figure 7a shows a significant increase in PD-1, CTLA-4, and PD-L1 expression in patients in the LR-G The survival of patients expressing high PD-1 levels in the LR-G was significantly longer compared to those expressing high PD-1 levels in the HR-G (P = 0.006; Fig. 7b).Additionally, the OS of patients expressing low PD-1 levels in the LR-G was significantly longer (P = 0.038; Fig. 7b).The survival patterns of patients divided into groups based on their risk scores, PD-L1 and CTLA-4 expressions in the training set were similar, as demonstrated in Fig. 7c and d.Furthermore, we assessed the ability of the risk model to predict the efficacy of immunotherapy in the IMvigor210 cohort.The response (P < 0.001; Fig. 7e) and OS (P = 0.025; Fig. 7f) of patients with low-risk scores to anti-PD-L1 drugs were significantly better compared to those in the HR-G.Further, in the LR-G, the response of patients (28.1%) to anti-PD-L1 therapy was better compared to that of the HR-G (17.4%;Fig. 7g).Results suggest that those with low-risk CC may gain more benefit from anti-PD-L1 treatment.

Correlation Between the Risk Model, Tumor Mutational Burden (TMB), and Drug Sensitivity
TMB is a reliable factor for evaluating the immune response against tumors.The results showed that the TMB scores of patients in the LR-G were higher than those in the HR-G (Fig. 8a), indicating that low-risk patients could benefit more from immunotherapy.Next, to assess the impact of TMB on clinical outcomes, we compared the survival rates of patients in the low and high TMB groups and found no significant difference between them (Fig. 8b).Additionally, we examined the frequency of somatic mutations in both groups.Figure 8c shows the TOP20 driver genes with the most mutations and their somatic mutation profiles.In both groups, the most frequently mutated genes were MUC16, PIK3CA, TTN, and KMT2C; however, the frequency of mutation of these genes was higher in patients in the LR-G compared to the HR-G.Currently, resistance to chemotherapy is a serious problem in cancer therapy.Therefore, we examined the therapeutic efficacy of 138 chemotherapeutic drugs in patients from the TCGA dataset.The responses of patients in the LR-G to 38 drugs, including 5-fluorouracil, gemcitabine, and ruxolitinib, were positive (Fig. 8d-f).In the HR-G, the response of patients to dasatinib, thapsigargin, WH-4-023, midostaurin, and TGX22 were positive (Fig. 8g-k).

Expression of Signature Genes in CC Cells and Tissues
The risk model comprised eight CSRGs.The survival of patients with CC expressing high levels of IL-1α and SERPINE1 was poor (Fig. 9a and c).These results were verified using data from the GEPIA database (Fig. 9b  and d).In addition, an increase in SERPINE1 and IL-1α expression was observed in patients in the HR-G (Fig. 9e  and f).These results suggest a close correlation between SERPINE1 and IL-1α expression and the occurrence and the development of CC.Therefore, we determined the expression of SERPINE1 and IL-1α in ten pairs of normal and CC tissues.An increase in SERPINE1 and IL-1α expression was observed in the CC tissues compared to the normal tissues (Fig. 9g and h).Additionally, compared to normal cervical epithelial cells, an increase in SERPINE1 and IL-1α expression was observed in CC cells (Fig. 9i and j) at mRNA as well as protein levels (Fig. 9k).

Discussion
Cellular senescence is a stable cell cycle arrest that remodels the TME through various cancer hallmarks, such as tumor proliferation, migration, invasion, angiogenesis, EMT, and tumor immune response, thereby affecting cancer patients' prognosis [25,26].However, the correlation between cellular senescence and CC is still unclear.In this study, we sought to explore the potential role of cellular senescence in the pathogenesis of CC, using bioinformatics and in vitro experiments, to better predict the prognosis of CC and provide novel targets for its treatment.Firstly, we explored the expression of senescence-associated genes in normal and CC samples using databases.We found a total of 23 genes associated with senescence that were differentially expressed in CC samples.These 23 genes include AGT, SOX2, E2F1, GNG11, SERPINE1, IL-1α, IFNG, CAVIN1, CDKN2B, CXCL1, HEPACAM, HSPB2, ID1, IGFBP5, ING2, KL, MMP9, PMVK, SFN, SIX1, TNFSF15, VENTX, and ZFP36.This result preliminarily illustrates the important role of cellular senescence in the occurrence and progression of CC.
With the development of high-throughput sequencing technology as well as computer algorithms, constructing gene sets for disease prediction has provided strong support for prognostic prediction of tumors, which is beneficial for optimizing treatment decisions in clinical practice.SASP is an important feature of cell senescence that includes cytokines, chemokines, growth factors, and proteases.Different SASP molecules serve different functions in the TME.We therefore analyzed the expression of SASP molecules in the HR-G and LR-G, and we found that multiple SASP molecules, including IL-6, IL-8, IL-1β, and VEGF-A, were increased in the HR-G.This is consistent with earlier findings.Multiple studies have shown that IL-6 is important for the cervical carcinogenesis.Pan et al. conducted an immunohistochemistry analysis of IL-6 expression in CC tissues, discovering a significantly elevated expression in the tumor tissues [29].Additionally, research has shown that IL-6 is abundantly expressed in invasive CC and is implicated in the pathogenesis of HPV-related CC [30].All the above findings suggest that IL-6 is a detrimental factor for the development of CC.IL-8 is a pro-inflammatory factor that promotes tumor growth.Fujimoto et al. [31] found that CC patients with high levels of IL-8 had an extremely poor prognosis, whereas those with lower levels had a better 24-month survival rate, indicating that IL-8 is a prognostic indicator of CC.In 2017, Jia et al. [32] found that IL-8 was associated with the tumorigenesis of CC, and exogenous IL-8 promoted the carcinogenic potential of HeLa cells.VEGF-A is a key factor in blood vessel formation, and previous studies have shown that the expression amount of serum VEGF-A is upregulated in CC, and targeting VEGF-A is beneficial for the treatment of CC [33,34].It therefore has the potential to be an effective treatment modality for cervical cancer by modulating SASP molecules in the TME.
Another important finding in our study is the significant correlation between CSRGs and the composition of the tumor-infiltrating immune cells.Many studies have shown that cellular senescence is associated with the TIME [35].Senescent cells secrete numerous cytokines and chemokines to induce immune cells and promote the body's immune response.Multiple studies have confirmed that most cervical cancers are HPV positive, and the body shows an antiviral immune response after infection with the HPV virus [35].Therefore, TIME is important for the development of cervical cancer, which deserves to be fully studied.In this work, GSEA analysis revealed the up-regulation of several pathways associated with immune responses, such as the natural killer cell-mediated cytotoxicity, the B cells, T cells, and Toll-like receptor signaling pathways, in patients in the LR-G.Our correlation analysis revealed that the stromal score had a positive correlation with the risk score, whereas the immune score had a negative correlation with the risk score.Based on this, we analyzed the composition of immunocytes, and we found that most immune cells were up-regulated in the LR-G, including B cells, CD8 + T cells, NK cells, and neutrophils.These results illustrate a more active and complex immune response in patients in the LR-G, which also lays the foundation for immunotherapy in CC.
Immunotherapy is a rapidly developing therapeutic strategy that holds tremendous potential in clinical settings.Immunotherapy targets and eliminates tumor cells by activating the patient's immune system.Studies have shown the efficacy of immunotherapy in treating various solid cancers, such as lung, breast, and renal [36][37][38][39].Additionally, studies have demonstrated the benefits of immunosuppressive agents targeting PD-1 and CTLA-4 or its primary ligand PD-L1 in treating patients with advanced and metastatic CC [40,41].Cellular senescencerelated immune remodelling could influence the efficacy of immune checkpoint blockade.Our results showed an increase in PD-1, CTLA-4, and PD-L1 expression in patients in the LR-G, thus indicating a higher sensitivity of these patients to immune checkpoint blockade therapy.Additionally, in the IMvigor210 cohort, patients with low-risk scores were highly sensitive to PD-L1 inhibitors.Therefore, the risk model could aid in screening patients who could benefit from a combination therapy.
We established a risk model based on eight CSRGs.KM analysis showed an independent association between SERPINE1 and IL-1α expression and the prognosis of patients with CC.SERPINE1 negatively regulates the pericellular proteolytic pathway.It has been observed that high levels of SERPINE1 expression are associated with a worse outcome and shorter disease-free survival in cancers such as breast and gastric [42,43].Hazelbag et al. used a multivariate Cox regression analysis and identified SER-PINE1 as a strong, independent prognostic factor for CC.Additionally, the study has shown an association between SERPINE1, poor survival, and disease recurrence in a subgroup of patients with CC without lymph node metastases [44].Interestingly, our results showed high SER-PINE1 expression in patients in the HR-G.The clinical outcome of patients with CC who expressed high levels of SERPINE1 was poor.In addition, SERPINE1 expression was higher in the CC cells and tissues.IL-1α is a crucial cytokine involved in inflammatory processes and promotes the pathogenesis of cancer.However, IL-1α exerts pro-and anti-cancer effects; hence its involvement in cancer progression is still controversial.Liu et al. [45] showed that IL-1α promotes breast cancer progression by increasing the activation of the NF-kB and STAT3 signaling pathways.Nevertheless, Dagenais et al. [46] showed that IL-1α suppresses breast cancer by inhibiting cell proliferation through the IL-1α signaling pathway.However, no study has reported the expression and functions of IL-1α in CC.Our results showed an increase in IL-1α expression in patients in the HR-G.The prognosis of patients expressing high levels of IL-1α was poor.Additionally, an increase in IL-1α expression was observed in the CC cells and tissues, thus indicating that IL-1α could promote malignant transformation of cells, thereby detrimental to the prognosis of patients with CC.In vitro as well as in vivo should be conducted to determine the involvement of SERPINE1 and IL-1α in CC.
However, our study has several limitations.First, we used data extracted from publicly available databases.Hence, prospective studies involving human subjects are required to validate our results.Additionally, cell-based and animal experiments should be performed to enhance our understanding of the mechanisms of CSRGs in the progression of CC.

Conclusion
We constructed a risk model comprising eight CSRGs to predict the prognosis of patients with CC.In addition, we used the risk model to determine the clinical outcomes and immune cell infiltration profiles of patients with CC.Finally, our risk model may help to design accurate and personalized therapeutic strategies for patients with CC.

Fig. 1
Fig. 1 Study design and flow diagram

Fig. 2
Fig. 2 Identification of prognosis-related DEGs.(a) The volcano plot shows DEGs in tissues from CC and normal samples from the TCGA-TARGET-GTEX dataset; red points indicate upregulation, and blue points indicate downregulation.(b) Venn plot shows 23 overlapping genes between the DEGs and the CSRGs.(c) The forest plot shows

Fig. 3
Fig. 3 Establishment of a risk model based on CSRGs.(a) PCA plot of patients from the TCGA dataset in the HR-G and LR-G (b) Disribution of the risk scores of patients with CC from the TCGA dataset.(c) The scatter plot shows the survival time and outcomes of patients with CC from the TCGA dataset.(d) Heatmap of eight genes in the prognostic signature.(e) KM analysis in two groups.(f) Time-

Fig. 4
Fig. 4 Clinical significance of the risk signature.(a) The forest plot shows the results of a univariate Cox regression analysis.(b) The forest plot shows the multivariate Cox regression analysis results of the prognostic signature and clinical variables.(c) Nomogram for assessing 1-, 3-, and 5-year OS rates.(d) The calibration curve for a nomogram.(e-g) Timedependent ROC curves of the nomogram, tumor stage, and the risk model for predicting the 1-, 3-, and 5-year OS rates.(h-j) DCA of the clinical stage, the risk scores, and the nomogram for predicting the 1-, 3-and 5-year OS rates

Fig. 5
Fig. 5 Correlation between the risk model and the SASP, as well as the enrichment scores of specific gene sets.(a) The expression of SASP factors in patients in the HR-G and LR-Gs.(b) GSEA revealed activation of the EMT, angiogenesis, and hypoxia pathways in

Fig. 6
Fig. 6 Evaluation of immune cell infiltration.(a) Multiple significant immune-related pathways were identified using GSEA.(b) Differences in the immune and stromal scores of patients in the HR-G and LR-G.(c) Stromal scores of patients in the two groups.(d) Immune scores of patients in both groups.(e) Infiltration of 16 immune cell types was measured using the "ssGSEA" algorithm.(f) Enrichment

Fig. 7
Fig. 7 Efficiency of the risk signature in predicting response to immunotherapy.(a) Expression of immune checkpoint genes like PD-1, PD-L1, and CTLA-4 in patients in both groups.(b-d) KM survival curves for assessing the survival of patients from the TCGA dataset stratified based on the risk scores and the expression of immune checkpoint genes like PD-1, PD-L1, and CTLA-4.(e) The

Fig. 8
Fig. 8 Correlation between the risk score, TMB, and chemosensitivity.(a) Comparison of TMB scores between two risk groups.(b) KM analysis in the two risk groups.(c) The waterfall plot shows somatic

Fig. 9
Fig. 9 SERPINE1 and IL-1α expression in CC cells and tissues (a, c) KM survival curves to determine the survival of patients categorized based on SERPINE1 and IL-1α expression.(b, d) KM survival curves to determine the survival of patients from the GEPIA database categorized based on SERPINE1 and IL-1α expression.(e, f) Relative