From TCGA database (https://portal.gdc.cancer.gov/), we obtained RNA sequence transcriptome data and relevant clinical information of 304 patients with CC and 3 normal adjacent tissue samples. From the GEO database (https://www.ncbi.nlm.nih.gov/geo/, GSE44001), we downloaded RNA sequencing data and clinical information of 300 patients for external validation.

The list of TCA cycle-related genes was obtained through literature mining (Arnold et al., 2022). Their mRNA expression levels between CC and normal adjacent tissue samples were compared according to TCGA cohort. The limma software package was used to identify the differentially expressed TCA cycle-related genes with the significance threshold (p < 0.05 and |log₂FC|>1.5), which were presented as a heatmap. The “corrplot” package was used to reveal correlations between DEGs associated with the TCA cycle. An interaction network of proteins among TCA cycle-related DEGs was constructed using STRING and visualized using Cytoscape 3.8.0.

The “ConsensuClusterPlus” R package was used for the analysis the comprehensive expression of the 18 differentially expressed TCA cycle-related genes to identify distinct subgroups of 302 CC samples. It was repeated 1,000 times to ensure classification stability (parameters: clustering algorithm, k-means; distance, Euclidean). The optimal k value was determined based on cumulative distribution function and delta area values. Principal Component Analysis (PCA) were performed by “Rtsne” R package to reduce the dimension of the 18 DEGs. The Kaplan-Meier method and log-rank test were used to evaluate the overall survival (OS) rate of patients with different subtypes. Chi-square test was used to analyze the distribution of age, tumor grade, tumor stage and histological type among different clusters.

To sort out TCA cycle related genes with potential prognostic value (p < 0.05), univariate Cox analysis was performed for OS. Next, using a least absolute shrinkage and selection operator (LASSO) regression model, the optimal value was determined to build a prognostic gene signature. We used R’s glment package to perform Cox regression analysis and LASSO. On the basis of the following formula: Risk score = ∑Coefgene × Expgenes, risk scores for every single patient was calculated, where Coefgene represents the coefficient of each prognostic gene and Expgenes represents the expression level of each gene. According to the median risk scores, patients were divided into high-risk and low-risk group. In addiction, we plotted the receiver operating characteristic (ROC) curves and Kaplan–Meier plots. To perform the validation of the prognostic model, GEO dataset (GSE44001) was analyzed the prognostic value with similar methods.

Gene set enrichment analysis (GSEA) (https://www.broadlnstitute.org/gsea/) was used to identify differential expressions of genes (gene sets) that were functionally related and whose enrichment in CC patient subgroups was significant. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway databases were downloaded from the molecular signatures database (MSigDB) as the functional enrichment reference set (http://www.gsea-msigdb.org). Finally, the significantly enriched KEGG pathways are shown centrally. Gene set variation analysis (GSVA) was implemented in the “GSVA” R package to investigate potential molecular characteristics that differed between high- and low-risk groups. Access to the MSigDB, the Hallmark gene set “c2. cp.kegg.v2022.1. Hs.symbols.gmt” was gotten to applied in GSVA. According to a threshold of |log₂FC| > 0.1 and p < 0.05, DEGs between high- and low-risk groups were screened and undergone to gene ontology (GO) and KEGG analyses using the “clusterProfiler” R package. Furthermore, we made a comparation the expression levels of cuproptosis-related genes between the high-risk group and the low-risk group, and visualized them by box plots.

Spearman correlation analysis was performed using “cor.test” in R to analyze the correlation between immune checkpoint expression with p < 0.05 as the cutoff for significance. The pRRophetic R package was used to predict chemosensitivity based on data from the Genomics of Drug Sensitivity in Cancer pharmacogenomics database. The half maximal inhibitory concentration (IC50) of clinically commonly used chemotherapy drugs was calculated to evaluate the role of predictive signatures in predicting the treatment response of CC. We compared the IC50 values between the high- and low-risk groups via Wilcoxon signed-rank test.

Both cervical cancer and adjacent non-cancerous tissues used in this study were obtained from postoperative patients with cervical cancer from 2019 to 2022 in Department of gynecology, Guangdong Women and Children Medical Hospital. All samples were obtained through review by the ethics committee, and the informed consent of CC patient was acquired. We extracted RNA from specimens by utilizing the TRIzol reagent (Ambion, United States), followed by reverse transcription into cDNA utilizing the QuantiTect Reverse Transcription Kit (Promega, United States). Quantitative PCR (qPCR) is a technique for measuring the amount of DNA present in a sampl in real time. With the aid of SYBR-Green (Vazyme, China), real-time qPCR assays were carried out, and expression levels were standardized to beta-actin levels. The sequences of primers are listed in Table 1.

All statistical analyses were performed with the use of R software (Version 4.2.1). Wilcoxon test was used to analyze the difference in the expression of TCA-related genes between normal and tumor tissues. Cox regression model was used for univariate and multivariate survival analysis to screen independent prognostic signature. The OS of patients in the high and low risk groups was analyzed by the Kaplan-Meier method and log-rank test. The ROC curve was drawn and the area under the curve (AUC) was determined by using the “survivalROC” software package. GraphPad Prism 9 program was used to draw scatter plots, and paired t-test was used to detect the differences in the expression of prognostic related genes between cervical cancer tissues and adjacent tissues. A p-value of less than 0.05 was considered to be statistically significant (p < 0.001 = ***, p < 0.01 = **, and p < 0.05 = *).

A list of 117 TCA cycle-related genes was identified (Supplementary Table S1), based on published data, and their RNA expression levels compared in TCGA data from 304 CC and 3 normal adjacent tissue samples. There were 18 differentially expressed TCA cycle-related genes identified, with a threshold of p < 0.05 and |log₂FC > 1.5|, of which 17 (PKM, GPI, IDH1, SHMT1, MTHFD1, SHMT2, ENO1, IDH2, ALDOA, DHFR, ELOVL3, SCD, TYMS, HK2, ALDOB, and PKLR) were upregulated and only one gene, ACAT1, was downregulated in tumor tissues (Figure 1A). Correlations among the mRNA expression levels of TCA cycle-related DEGs were analyzed by Pearson correlation analysis (Figure 1B). The results showed that all the TCA cycle-related DEGs had a positive correction with each other. In particular, FASN was significantly correlated with SCD (r = 0.66, p < 0.001) and PKM was significantly correlated with ENO1 (r = 0.63, p < 0.001). Construction of a PPI network revealed that the top 5 genes including PKM, IDH1, ENO1, PKLR and IDH2 were selected, based on their values of closeness, to be the hub nodes in the PPI network (Figure 1C).

To explore the relationships between CC subtypes and expression of the 18 TCA cycle-Related DEGs, consensus clustering analysis was performed to classify tumors according to expression levels of TCA cycle -related DEGs. Clustering variable (k) values from 2 to 9 were applied; when k = 2, intragroup correlations were low. Hence, patients with CC could be divided into two different subtypes, including 215 cases in cluster 1 and 87 cases in cluster 2 (Figures 2A, B). PCA was conducted to verify the ability of the model to group patients in the entire set and observed that patients in different clusters were dispersed in two directions (Figure 2C). There was a significant difference in OS time between the two clusters (p = 0.041) (Figure 2D). Further, the associations between the clustering and clinicopathological parameters were examined. The significant difference was found between cluster 1 and cluster 2, for the survival state (p < 0.05) and pathological type (p < 0.01). In contrast, other parameters such as age, tumor grade and clinical stage were no significant different (Figure 2E).

The TCA cycle-related genes were all chosen for the univariate Cox regression analysis, and we found that 12 genes were significantly associated with OS in TCGA CC (Figure 3A). A LASSO regression analysis was applied to establish a prognostic gene signature using the 12 genes mentioned above. Following LASSO analysis to minimize overfitting (Figures 3B,C), seven genes involving ACSL1, ALDOA, FOXK2, GPI, MDH1B, MDH2 and MTHFD1 were identified (Figure 3D).

The risk score of seven genes was also calculated for further univariate and multivariate Cox regression analyses. The risk score formula to predict OS was developed as follows: risk score = (−0.080 × ACSL1) + (−0.036 × ALDOA)+ (0.014 × FOXK2) + (−0.068 × GPI) + (0.102 × MDH1B) + (−0.216 × MDH2) + (−0.089 × MTHFD1). It is well-known that survival times vary among patients with different pathological types of CC. Thus, prognosis analysis of the seven genes in different pathological types of CC including cervical squamous cell carcinoma and cervical adenocarcinoma. Using this signature, patients were further classified into equal high- and low-risk groups, based on the median risk value (Figures 4A, E). As illustrated in the scatter diagram in Figures 4B, F, individuals in the high-risk score group had worse outcomes than those in the low-risk group. In addition, a significant difference in OS time was detected between the two groups by Kaplan-Meier analysis (p < 0.01, p < 0.05) (Figures 4C, G). ROC curve analysis was conducted to evaluate the sensitivity and specificity of the prognostic model, resulting in AUC values of the models for predicting 1-, 3-, and 5-year OS of 0.613, 0.663, and 0.736 in cervical squamous cell carcinoma, while 1-, 3-, and 5-year OS of 0.699, 0.663, and 0.633 in cervical adenocarcinoma respectively (Figures 4D,H).

To test the robustness of the gene signature model built from the TCGA data, data from 300 patients with CC in the GEO cohort, the GSE44001 dataset were also divided into high- and low-risk groups using a similar formula to TCGA data (Figure 5A). According to the uniform formula, the survival analyses found that patients with higher risk scores had poorer OS (p = 0.001) (Figures 5B, C). In the GSE44001 dataset, the AUC was 0.705 at one year, 0.701 at three years and 0.68 at five years (Figure 5D).

To explore the potential biological processes in high- and low-risk groups, we performed a GSEA. The KEGG pathway analysis showed that phototransduction, RNA polymerase, and steroid biosynthesis were mainly enriched in the low-risk group (Figure 6A), while allograft rejection, glycosaminoglycan biosynthesis—keratan sulfate and other glycan degradation were principally enriched in the high-risk group (Figure 6B). To further identify the expression difference of these two groups, the GSVA enrichment analysis revealed that cancer pathways, including Wnt, Notch and mTOR signaling pathway were highly expressed in low-risk group, compared with high-risk group (Figure 6C). These results suggest that Metabolic reprogramming modulates tumor proliferation, apoptosis, and cell cycle via these pathways. The box plot of Cuproptosis-related gene analysis illustrated that in the low-risk group, CDKN2A DLAT DLD GLS LIAS MTF1, and PDHA1 were significantly downregulated in the high-risk group (Figure 6D). The finding echoed the definition of the novel cell deathmodality “Cuproptosis” which is featured by disturbing specific mitochondrial metabolic enzymes (Tsvetkov et al., 2022).

Correlation assessment of the association between risk score and immune checkpoint-related genes found that PD-L1, 4-1BBL, OX40L, GITR, B7.1, and B7.2 had a negative correction with the risk score (Figure 7A). In other words, CC patients with higher risk scores had lower expression levels of these immune checkpoints. Our data suggest that patients in low-risk group may be more sensitive to immunotherapy. In addition to immunotherapy, we also analyzed the association between the predictive signature and the efficacy of general chemotherapy for CC. The results found that the IC50 of sunitinib, paclitaxel, cisplatin, and docetaxel was lower in the low-risk group (Figure 7B), which is helpful for exploring individualized treatment schemes suitable for high- and low-risk group patients.

To determine whether the seven prognostic genes were differentially expressed in CC tissues, a total of 19 paired clinical CC tissues and adjacent normal tissues were analyzed each gene expression using qRT-PCR. The findings illustrated that the expression levels of ACSL1, ALDOA, FOXK2, MDH2, and MTHFD1 in cervical cancer specimens were differentially expressed in contrast with those in normal specimens, whereas there was no significant difference of the expression levels of GPI and MDH1B (Figure 8).