3.1 Expression change of CRGs and consensus clustering analysis of ovarian cancer patients

In this work, 10 classical CRGs were confirmed in the literature. Based on the TCGA and Genotype-Tissue Expression (GTEx) datasets, the heatmap revealed the differential mRNA expression of 10 cuproptosis-related genes (CRGs) in OC and normal samples. The results indicated that 6 genes (LIAS, LIPT1, DLD, PDHA1, PDHB, and CDKN2K) had significant differential expression (Fig. 1A). Our further summary analysis of the incidence of somatic mutations in these 10 CRGs revealed that 12 patients (2.75%) in the OC cohort had CRG mutations (Fig. 1B). Among them, the mutation frequency of DLAT was the highest, followed by MTF1, LIAS, LIPT1, CDKN2A, and GLS, while the other 4 CRGs had no mutation. Subsequently, we explored somatic copy number changes in these CRGs and detected common copy number changes in all 10 CRGs. Among them, DLD, MTF1, GLS, and LIPT1 had a high degree of copy number gain, while LIAS and PDHA1 had copy number loss (Fig. 1C). The locations of the copy number variations (CNVs) in CRGs on their respective chromosomes are shown in Fig S1A.

Fig. 1

Expression and genetic alteration of CRGs in OC: A the expression of 10 CRGs in OC and normal ovary tissues; B, C the CNV and mutation frequencies and classifications of 10 CRGs in OC patients from TCGA cohort. D Unsupervised clustering of 10 cuproptosis-related genes in TCGA OC cohort. Cuproptosis cluster, live status, age, recurrence, stage, grade, and BRCA mutation were utilized as patient annotations. E PCA analysis of two cuproptosis subtypes in TCGA OC cohort. *p < 0.01, ***p < 0.001

To further illuminate the clinical features or biological behaviors of these CRGs, we performed an unsupervised clustering algorithm and divided the TCGA OC cohort into 2 groups (k = 2) based on the expression profiles of 10 CRGs; there were 268 patients in cuproptosis-related Cluster 1 and 95 patients in cuproptosis-related Cluster 2 (Fig S1A, B, C). Cluster 1 had a survival advantage over Cluster 2 according to Kaplan–Meier survival analysis (log-rank test, p = 0.048; Fig S1D). We also explored the clinicopathological characteristics of these patients. The heatmap showed no significant difference in BRCA mutation and recurrence between the 2 clusters (Fig. 1D). Principal component analysis (PCA) demonstrated distinct separation in the cuproptosis transcriptional profiles between the 2 clusters (Fig. 1E).

3.2 Identification immune infiltration characteristics and differentially expressed genes between different clusters

As shown in Fig. 2A-D, the ESTIMATE algorithm was used to estimate the ratio of the immune matrix component of the tumor microenvironment (TME) based on the TCGA expression profiles; TME features were reflected by the immune score, stromal score, and estimate score, which are positively correlated with the infiltration level of immune cells and the content of matrix components [31]. The results showed that the immune score and estimate score was significantly higher in Cluster 1 than in Cluster 2, while no significant differences were shown between the two clusters in terms of stromal score. We further explored the relationship between clusters and immune cell infiltration established in the TCGA cohort by using the ssGSEA algorithm (Fig. 2E). The great majority of immune cells exhibited significantly higher infiltration in Cluster 1 samples than in Cluster 2 samples; these cells included memory B cells, dendritic cells, macrophages, CD4 T cells, CD8 T cells, natural killer cells, follicular helper T cells, and Tregs.

Fig. 2

Identification immune infiltration characteristics and DEGs in the two clusters. A–D Stromal score, immune score, estimate score, and purity between the OC clusters. E The GSA scores for each TME infiltrating cell in the two OC clusters. The line in the box represents the median value. F, G The volcano plot and heat map show the unsupervised clustering of DEGs in the TCGA OC cohort. H GO and KEGG enrichment analysis of DEGs. GO, Gene Ontology. KEGG, Kyoto Encyclopedia of Genes and Genomes. DEGs differentially expressed genes. TME, tumor microenvironment; ssGSEA, single sample gene set enrichment analysis. *p < 0.01, ***p < 0.001

To further investigate the conceivable biological processes of the cuproptosis-related clusters in the OC cohort, we utilized the R package “DEseq2” to identify the differentially expressed genes (DEGs) among the 2 clusters; 783 genes were highly expressed in Cluster 1 and 595 genes were highly expressed in Cluster 2 (Fig. 2F, G). GO enrichment analysis of the DEGs demonstrated that these DEGs were mainly related to some immune and carcinogenesis biological processes, such as T-cell activation, negative regulation of immune system process, regulation of T cell, signaling receptor activator activity, DNA − binding transcription activator activity, and immune receptor activity. Moreover, in the KEGG pathway enrichment analysis, these DEGs were mainly related to cytokine − cytokine receptor interaction, natural killer cell-mediated cytotoxicity, Th1 and Th2 cell differentiation, and cell adhesion molecules (Fig. 2H).

3.3 The cuproptosis-related scoring model exhibits great prognostic prediction ability

Using univariate Cox regression analysis, 33 prognostic cuproptosis-related genes were identified among the 1378 DEGs (Fig S2A). The cuproptosis scoring model was constructed depending on these prognostically significant DEGs. Least absolute shrinkage and selection operator (LASSO) Cox regression was applied to screen out the 33 prognostic DEGs to determine the optimal penalty parameter λ, and 21 genes were finally selected. A cuproptosis-related prognostic model was then established on the basis of the expression profiles of the 21 prognostic genes, and this model was named the CuRScore model (Fig. 3A–C).

Fig. 3

Construction of the CuRS model (A, B) Cvfit and lambda curves illustrating the LASSO regression model produced using tenfold cross-validation. C Coefficients of the 21 prognostic molecules in the Cox regression model. D The heatmap representing differences in the clinical information and genes expression between the high and low GuRS clusters. E Kaplan–Meier plots of OS for the two GuAS clusters in the TCGA data. F 1-, 2-, and 3-year ROC of GuAS prognostic performance in the TCGA data. G Ranked dot and scatter plots for the relationship between GuAS and prognosis status in the TCGA data. H, I, J Kaplan–Meier plots, ROC curve, and scatter plots of GuAS prognostic performance in GSE32062 data. K GuRScore differences among different clinical features in the TCGA data. L An alluvial diagram of the distribution of cuproptosis cluster, gene cluster in two risk groups, as well as survival outcomes.*p < 0.01, ***p < 0.001

Based on the median CuRScore, the samples in the TCGA and GSE32062 OC cohorts were further stratified into high-risk or low-risk subgroups (Fig. 3D). Patients in the high-risk subgroup showed substantially shorter overall survival than those in the low-risk subgroup according to the Kaplan‒Meier method in both cohorts (Fig. 3E, H). The area under the curve (AUC) was generated to compare prognostic ability within different models, and analysis indicated that CuRS possessed relatively high AUC values for 1-, 2-, and 3-year overall survival (0.655, 0.719, and 0.706, respectively) (Fig. 3F, G). Additionally, the AUCs for risk variables were 0.679, 0.641, and 0.584 in the GSE32062 OC cohort, which demonstrated a relatively high excellent predictive ability for OC patients (Fig. 3I, J). Moreover, we found that older age (> 60), advanced stage (stage III), death, and BRCA mutation were closely correlated with a high CuRS, whereas grade and recurrence revealed no relationship with CuRS (Fig. 3K). Moreover, considering the clinicopathological characteristics, CuRS was confirmed as an independent predictive factor for OC in the TCGA and GSE32062 cohorts by univariate and multivariate Cox regression analyses (Fig S3 A–D). Figure 3L illustrates the distribution of patients in the two cuproptosis clusters and two CuRScore subgroups. The majority of patients in cuproptosis cluster 2 were also in the low-risk group, which had a good prognosis, similarly, cluster 1 with high-risk score patients had the worst survival.

3.4 Analysis of CuRS model using single‑cell RNA sequencing

After quality control, we utilized signal cell RNA-seq data (GSE17368 and GSE15893) to acquire gene expression profiles of at least the 21 signature genes for 6426 cells from 7 OC samples. Based on these malignant epithelial cells, a total of 15 clusters were identified, and 8 subtypes were confirmed based on CNVs (Fig S4A). Moreover, according to the median CuRS value, the 15 clusters of single cells were classified into high and low groups (Fig S4B). The allocation of the CuRScore was assessed and mainly existed between − 0.1 and 1 (Fig. 4A). Sorting of these cells according to the gene expression of each cell enabled further inference of the cell development pathway by the Monocle 3 algorithm. Consequently, a lineage from low CuRS OC cells to high CuRS OC cells was identified (Fig. 4B), and this lineage was consistent with the evolutionary pathway of the CuRS cluster. Copy number variation is a gradual development process, and according to the chronological order, the tumor cells that are closest to the copy number profile of normal cells should appear at an earlier time. As shown in Fig. 4, the CuRScore is positively correlated with tumor cell progression. Moreover, cluster 4, as the subtype with the lowest CNV frequency, was more likely than other clusters to be normal cells at represent the beginning of the progression of carcinogenesis (Fig S5A).

Fig. 4

Single-cell clustering analysis and pseudotime analysis in OC cohort. A The frequency distribution of the CuRscoring in OC cells. B The pseudotime trajectory in OC cells

3.5 Transcription factor differentially activated in high and low CuRS groups

Furthermore, we investigated the activation of transcription factors in two CuRS groups of OC cells. Because transcription factors mutually cross-regulate the expression of some genes, we grouped them into different modules (M1, M2, M3, M4, M5) (Fig S6A). Each module represented a group of collaborative transcription factors. The regulon activity score (RAS) of the respective module was computed according to the Connection Specificity Index defined in a previous study [32]. It appeared that M3 had a higher RAS than other modules in both CuRS groups, while the RAS of M5 was relatively low. Moreover, the RAS of M2 was negatively correlated with CuRS. There was no available expression diversity in the RAS of M4 and M1 between the high and low CuRS groups (Fig. 5A). The scores of the activity of regulons of these modules were classified with the AUCell algorithm and then shown by t-SNE plots. The M2 module was more activated in the low CuRS group (Fig S6B). We also analyzed the top 10 transcription factors in the high and low CuRS groups using the regulon specificity score (RSS).

Fig. 5

Transcription factor activation difference in two GuAS cluster cells. A The correlation between the scoring system and those modules relies on RNA velocity. B Top 10 differential activated TFs in the two GuRS clusters respectively

Then, we selected the top 10 transcription factors from the high (ELK3, ETS1, ETV4, NR2C2, PML, GABPB1, POLR2A, ELK4, HCFC1, and ETS2) and low (EGR3, EGR1, ATF3, HOXD9, DDIT3, FOS, JUN, JUNB, CHD2, and ELF3) CuAS groups for further analysis (Fig. 5B). The figure S6C, D shows the top 5 transcription factors enrichment distribution map in the high and low CuRS groups. Meanwhile, we performed GO and KEGG pathway enrichment analyses on the corresponding target genes based on the top 5 differential transcription factors. The transcription factors of high CuRS were significantly enriched in immune-related pathways, while TFs of cells with low CuRS were significantly enriched in pathways such as apoptosis. Moreover, tumor transcriptional dysregulation pathways were significantly enriched in both high and low CuRS cells (Fig S7).

3.6 Immune related pathways selectively activated in high and low CuRS groups

We further performed an enrichment analysis of 29 immune-related pathways using the GSVA algorithm based on bulk RNA-seq data. Correlation analysis demonstrated that the CuRS with prognostic signatures was connected with most of the 29 immune-related pathways, and similar results were found in the scRNA-seq data (Fig. 6A). The KEGG enrichment analysis based on the GSVA analysis on OC patients from bulk RNA-seq data demonstrated that the high CuRS cluster was significantly related to the activation of immune-related pathways, including antigen processing and presentation, the intestinal immune network for IgA production, and Th17 cell differentiation. In addition, according to the scRNA-seq data, these results showed that the high CuRS cluster was significantly associated with the oxidative phosphorylation metabolic pathway, ribosome, and TNF signaling pathway (Fig. 6B). Furthermore, the proportion of infiltrating immune cells was analyzed using the CIBERSORT algorithm and was found to be significantly different between the high and low CuRS clusters for cells such as M0 and M1 macrophages, CD4 T cells, CD8 T cells, activated NK cells, and activated mast cells (Fig. 6C, D).

Fig. 6

Immune-related pathway analysis rely on bulk RNA-seq and scRNA-seq data in OC. A The heatmap represents differences in the clinical information and immune-related pathways in bulk RNA-seq analysis (left) and scRNA-seq analysis (right). B KEGG enrichment analysis based on the GSVA algorithm in bulk RNA-seq (left) and scRNA-seq analysis (right). C Stromal score, immune score, estimate score, and purity between the high and low GuAS clusters. D The TME infiltrating cell base on the CIBERSORT algorithm in the two GuAS clusters. The line in the box represents the median value. *p < 0.05, **p < 0.01, ***p < 0.001

3.7 Novel ligand–receptors pairs between immunocytes and the CuRS model

Then, we analyzed the intercommunication patterns between OC cells and immunocytes by using the CellChat algorithm based on single-cell RNA-seq analysis (Fig. 7A). As illustrated, there were specific communication patterns for both incoming and outgoing signaling in high and low CuRS cells, where the SEMA4 pathway had significant communication intensity differences between the high and low CuRS groups in terms of outgoing signaling patterns. Additionally, regarding incoming signaling patterns, the SELL pathway had a significant difference in communication intensity (Fig. 7B). The Notch pathway exhibited a significant difference in sender, receiver, and influencer patterns (Fig. 7C). Figure 7D illustrates the communication pattern distribution of the signaling pathways in OC cells and immunocyte cells. Furthermore, we utilized the CellTaker algorithm to identify ligand‒receptor pairs between immunocytes and CuRS cells. As shown in Fig. 7E, OC cells with high CuRS communicate with macrophages, T cells, CD34 + B cells, fibroblasts, and monocytes more actively than OC cells with low CuRS, which may explain the immune landscape difference.

Fig. 7

Novel ligand–receptor pairs difference between high and low CuRS clusters. A The cross-talk intercommunication among CuRS OC cells and different immunocytes. B The correlation between CuRS groups or immunocytes and different signaling pathways. C NOTCH signaling pathway network. D Incoming and outcoming communication patterns of target and secreting cells. E ligand–receptor pairs among different cells

3.8 Chemotherapy differential sensibility in high and low CuRS groups

To further explore the correlation between the CuRS subgroups and the effectiveness and sensitivity of chemotherapy, the OC chemotherapy cohort with considerable absolute clinical data was investigated based on the imvigor210 database. Patients in the high CuRS group had a substantially worse prognosis in the imvigor210 OC cohort (Fig. 8A). However, there were no significant differences in the chemotherapy drug response rate, CR/PR rate, or SD/PD rate between the high and low CuRS groups (Fig. 8B, C). These results suggested that the CuRS could be utilized to predict the efficacy of OC patient chemotherapy survival. Additionally, we analyzed the differences in chemotherapy and targeted therapy drug resistance between the high and low CuRS groups (Fig. 8D). Based on the CCLE database, the IC50 values of cisplatin, topotecan, cyclophosphamide, docetaxel, and 5-fluorouracil were higher in patients with high CuRS, and there was a significant positive correlation between IC50 and CuRS (Fig. 8E). Remarkably, similar results were obtained from the GDSC database analysis. In the GDSC database, we identified 172 kinds of chemotherapy drugs with a significant positive correlation between IC50 and CuRS subgroups. Among them, 13 kinds of chemotherapy drugs had correlation coefficients greater than 0.3, and no chemotherapy drug IC50 had a significant negative correlation with CuRS (Fig S8). We also focused on the relationship between CuRS subgroups and resistance to the PARP inhibitors olaparib and niraparib. In both the CCLE database and GDSC database, in patients with high CuRS, the IC50 values of olaparib and niraparib were significantly increased (Fig. 8F, G). In conclusion, these results suggest that drug IC50 has an excellent positive correlation with CuRS and CuRS was associated with chemotherapy and targeted therapy drug resistance.

Fig. 8

The correlation between CuRS and drugs sensitivity. A Kaplan–Meier curve of OS with high and low GuAS clusters for IMvigor210 cohort. B, C Boxplot and Bar’s graph shows the therapeutic response [complete response (CR)/partial response (PR) and stable disease (SD)/progressive disease (PD)] to immunotherapy in the two GuAS clusters in IMvigor 210 cohort. D The relationship of drug sensitivity and CuRS score. E Drugs characterized by IC50 positively correlated with CuRS score in the CCLE database. F, G PARP inhibitors sensitivity positively correlated with CuRS score in the CCLE and GDSC database. ***p < 0.001

3.9 Knockdown of KIF26B increased the sensitivity of OC cells to cuproptosis agonists

According to the above results, we further explored whether the expression of cuproptosis model related genes contributed to platinum resistance in OC. The online kmplot database (http://kmplot.com/analysis/) analysis results suggested that high expression of KIF26B, which among cuproptosis-related genes, predicted poorer prognosis in ovarian cancer patients treated with cisplatin (Fig. 9A). The IHC showed that the protein level of KIF26B in ovarian cancer was significantly higher than that in normal ovarian epithelial tissues. Meanwhile, the level of KIF26B protein in platinum-resistant OC samples was significantly higher than that in platinum-sensitive OC samples (Fig. 9B). To clarify the biological function of KIF26B, we constructed a cisplatin resistant ovarian cancer cell line A2780 (A2780/DDP). The expression levels of KIF26B in OC cell lines and the efficiency of knockdown KIF26B were confirmed by quantitative real-time PCR (Fig. 9C). The CCK-8 assay and clonogenic assay showed that knockdown of KIF26B expression significantly reduced ovarian cancer cell proliferation (Fig. 9D, E). In line with these findings, the apoptosis assay was utilized to validate the above results (Fig. 9F, G). Moreover, under the treatment of Elesclomol-CuCl2, the CCK-8 assay results revealed that the sensitivity of OC cells to cuproptosis agonist was increased after knockdown of KIF26B (Fig. 9H). These results thus suggested that the knockdown of KIF26B facilities cuproptosis in OC cells.

Fig. 9

KIF26B is related to copper ionophore-induced cell death. A Kaplan–Meier analysis of the OS and PFS of OC platinum chemotherapy patients with KIF26B high or low expression level. B Representative IHC images showing KIF26B expression in normal ovarian epithelial tissues, platinum-sensitive, and platinum-resistant OC specimens. C the expression level in IOSE80, A2780, A2780/DDP and after KIF26B knockdown. D CCK-8 assay analyzes cell viability after KIF26B knockdown treatment with cisplatin. E Colony formation and quantification analysis after KIF26B knockdown treatment with cisplatin. F, G Representative apoptosis assay and quantification analysis after KIF26B knockdown treatment with cisplatin. H CCK-8 assay analyzes cell viability after KIF26B knockdown treatment with elesclomol-CuCl2. *p < 0.05