Gene expression, prognostic, clinical, and pathological data of gliomas were downloaded from The Cancer Genome Atlas (TCGA) (gene expression: 701 glioma and 5 adjacent samples; clinical: 1,104 samples; mutation: 984 samples) and Gene Expression Omnibus (GEO) (50 glioma samples) databases. The datasets were managed by “Combat” algorithm to eliminate batch effects and then combined.

A consensus clustering algorithm based on Partitioning Around Medoid (PAM) algorithm was be used to categorize patients into cuproptosis subtypes. An analysis of patients’ clinicopathological characteristics and prognosis was carried out in order to determine the clinical features of the two subtypes. Cohort information included age, sex, glioma grade, and cuproptosis-related gene expression. We plotted the Kaplan-Meier curve with “survival” and “survminer” packages. Using “pheatmap” package, we generated a heatmap which would intuitively show the patient characteristics. GSVA and ssGSEA analyses revealed the TME differences between the two subtypes. ESTIMATE algorithm was used to calculate the immune and stromal scores of the cohort patients.

We used “limma” package to analyze differentially expressed genes between the two cuproptosis subtypes, with a standard of logFCfiter = 0.585 and adjusted p-value fliter = 0.05. To determine the potential functions of these genes, we used “clusterProfiler” and “enrichplot” to perform functional annotation of GO and analysis of KEGG pathway enrichment.

“Caret” package was used to randomly assigned the cohort patients to a training and test group (ratio of 1:1). We then analyzed the differentially expressed genes using LASSO and multivariate Cox analyses to select prognostic genes related to the gene subtypes. “Glmnet” package was used to avoid over-fitting. C R G − s c o r e = ∑ i c o e f   i   ∗   g e n e   i   e x p r e s s i o n We calculated the CRG-score for each patient. Patients were divided into a low- and high-risk group based on the median value of CRG-score.

The predictive capacity of the CRG-score was verified with multiple R packages. We used “ggpubr” package to identify the risk score difference between cuproptosis subtypes and gene subtypes. Survival analysis was performed by “survival” and “survminer” packages. The receiver operating characteristic (ROC) analysis was based on “timeROC” package. We also used the “pheatmap” package to reveal CRG-score model genes expression in the low- and high-risk group. Sankey diagram was based on “ggplot2” and “ggalluvial” packages, which showed an entire network of the subtypes, risk of CRG-score, and survival status.

Using ESTIMATE algorithm, stromal and immune scores of the low- and high-risk groups were evaluated. Additionally, we used CIBERSORT algorithm to calculate human immune-cell subsets of each patient in the two risk groups. To analyze mutation and tumor mutation burden (TMB), we used “maftools” and “ggpubr” packages. Curve of cancer stem cell analysis was made by “ggpubr,” “ggplot2” and “ggExtra” packages.

A nomogram was established by the “rms,” “regplot” and “survival” packages. The patients’ age, sex, glioma grade, and risk cluster were given a specific score, and the total score of each patient was the sum of these various factor scores. Based on a patient’s total score, the nomogram predicted 1-, 3- and 5-year survival accurately, which has been evaluated by calibration plots.

Ninety five glioma tissue samples were donated by 83 patients who had undergone surgery at Xiangya Hospital, including 27 grade 2, 6 grade 3, 50 grade 4, and 12 adjacent tissues. This study was approved by the Ethics Committee of Xiangya Hospital of Central South University.

TRIzol (Invitrogen, 15596018) was used to extract RNA from glioma tissue after grinding in liquid nitrogen. Then we converted RNA to cDNA by the RevertAid RT Reverse Transcription Kit (Thermo Scientific, K1691). The RT-qPCR data were analysed by the 2^−ΔΔCT strategy and values were relative to the housekeeping gene Glyceraldehyde-3-Phosphate Dehydrogenase (GAPDH). A complete list of the primers used is provided in in Supplementary Table S1. The 95 glioma tissue samples were manufactured using a tissue microarray. The following antibodies were used for IHC: FDX1 (Proteintech, 12592-1-AP), PDHA1 (ABclonal, A1895), and DLST (ABclonal, A13297). The grading rules as follows: The intensity was scored as follows: 0, negative; 1, weak; 2, moderate; and 3, strong. The frequency of positive cells was defined as follows: 0, less than5%; 1, 5%–25%; 2, 26%–50%; 3, 51%–75%; and 4, greater than 75%. The final grade of staining was determined by multiplying the score for staining intensity with the score for the frequency of positive cells (values, 0–12). This grading rules can reduce the grading bias caused by different cell density.

We used Student’s t-test, Log-rank test, Cox regression model, GSVA analysis, ssGSEA analysis, PCA analysis, Least absolute shrinkage and selection operator (LASSO), multivariate Cox analysis, Spearman test and Wilcoxon test. All the statistical analyses were based on R (version 4.1.0). And the cutoffs for the Kaplan-meiers in this study was found by using the R package “survminer”.

19 cuproptosis-related genes were identified based on previous research. We analyzed the transcriptome profiling data of TCGA-GBM and TCGA-LGG cohorts to identify differences in mRNA expression between glioma and adjacent tissues. The expression of GLS (p < 0.001), LIPT1, FDX1, SLC31A1 (p < 0.01), DLST, CDKN2A, PDHA1, ATP7A, ATP7B, and NFE2L2 (p < 0.05) was different between glioma and adjacent tissues (Figure 1A). A summary analysis of simple nucleotide variations in the GBM and LGG cohorts showed that of the 984 samples, 45 (4.57%) had mutations in the cuproptosis-related genes. Only NLPR3, ATP3A, CDKN2A, ATP3B, MTF1, and GLS had an approximately 1% mutation frequency (Figure 1B).

Next, we analyzed copy number variations (CNV) and discovered that all 19 cuproptosis-related genes had CNVs. Among these genes, CDKN2A, ATP7B, and DLST had widespread CNV loss, especially CDKN2A, which had a CNV loss frequency of >40% (Figure 1C). The comparison of 19 cuproptosis-related genes CNV with all CNV was showed in Supplementary Table S2. Figure 1D shows location of 19 cuproptosis-related genes on the chromosomes and their CNV alterations. However, the CNV changes do not adequately explain differences in mRNA expression between the adjacent and glioma tissues. This suggests that not only CNVs affecting the expression of cuproptosis-related genes, but other factors, such as DNA methylation and transcription, may affect the mRNA expression of cuproptosis-related genes. There is an association between CDKN2A mutations and the development and recurrence of gliomas (Varn et al., 2022). It could be an important indicator for molecular diagnosis, such as IDH mutation or MGMT methylation.

We collected 3 eligible glioma cohorts (TCGA-LGG, TCGA-GBM, and GSE43378) to further study the expression patterns of cuproptosis-related genes that partake in tumorigenesis. A total of 725 patients were integrated from the 3 cohorts for further analysis (31 samples were removed because of the lack of complete clinical or gene expression information). Detailed patient information is presented in Supplementary Table S3. “Survival” and “survminer” packages were used to process the clinical data. The results of Kaplan-Meier analysis revealed the prognosis of cuproptosis-related genes, with p < 0.05. 14 of the 19 genes were thought to be associated with the prognosis of gliomas. Higher expression of ATP7A, CDKN2A, DLAT, DLD, DLST, FDX1, GLS, LIPT1, NFE2L2, NLRP3, and SLC31A1 predicted poorer survival. ATP7B, LIAS, and PDHA1 showed the opposite trend (Figure 2A). The Cox regression models of 19 cuproptosis-related genes were constructed to improve the survival analyses (Supplementary Table S4). The cuproptosis-related gene interactions, prognosis, and regulator connections are comprehensively illustrated in the cuproptosis network (Figure 2B).

With the aim of examining the role of cuproptosis-related genes in glioma, a consensus clustering algorithm based on the PAM algorithm was used to categorize patients with glioma by analyzing their cuproptosis-related gene expression conditions. According to the algorithm results, it appears that when k = 2 the glioma cohort was optimally sorting into subtypes A and B (Figure 3A). The other consensus matrix, where k is valued from 3 to 9 were showed in Supplementary Figure S1. Subtype A included 333 patients, and subtype B included 392 patients. We used Kaplan-Meier analysis to compare the 2 subtypes, and the curves showed that compared to patients with subtype B, those with subtype A had a better survival probability (p < 0.001, Figure 3B). PCA analysis also revealed that there were differences between the 2 subtypes on cuproptosis-related gene expression profiles (Figure 3C). The heatmap showed that both cuproptosis-related gene expression and pathological characteristics are clearly different in the 2 subtypes (Figure 3D). More patients with grade 4 glioma were classified as subtype B, whereas more patients with grade 2 glioma were classified as subtype A. The heatmap also showed that most of cuproptosis-related genes were highly expressed in subtype B, while LIAS, ATP7B, PDHA1, and CDKN2A were expressed higher in subtype A. Also subtype B had more patients aged >65. These characteristics may explain the prognostic differences between the 2 cuproptosis subtypes.

GSVA analysis was used to explore the TME differences between the 2 cuproptosis subtypes (Figure 3E). Subtype A converged on lipid metabolism, particularly linoleic acid and retinol. Subtype B converged on pathways associated with the metabolism and damage repair of hereditary substances, including mismatch repair, homologous recombination, and nucleotide excision repair significantly. In addition, cell cycle and protein metabolism pathways are worthy of attention in subtype B. GSVA analysis revealed that the cell cycle and metabolism may cause different prognoses between the 2 cuproptosis subtypes.

In tumorigenesis and development of glioma, tumor immune microenvironment plays a critical part. Understanding how cuproptosis relates to the immune landscape in the tumor microenvironment may provide new avenues to treat cancer (Klemm et al., 2020). Therefore, we used ssGSEA to determine whether there was an immune infiltration difference between the cuproptosis subtypes (Figure 3F). A significant difference was found in the level of immune infiltration based on the analysis. Subtype B showed a higher infiltration level, especially in activated CD4⁺ and CD8⁺ T-cell, gamma delta T-cell, natural killer (NK) T-cell, and NK cells. We then used “estimate” package to compute the TME score of two subtypes (Figure 3G). In contrast to patients with Subtype A, those with Subtype B got higher stromal, immune, and ESTIMATE scores (p0.001). This suggests that the TME of patients with Subtype B has a higher proportion of stromal cells or immune cells.

The results showed an association between cuproptosis-related genes expression and glioma prognosis. We constructed 2 subtypes based on cuproptosis. The analysis revealed that the cell cycle, tumor metabolism, and immune infiltration may be related to cuproptosis in glioma.

We used the “limma” package to analyze differentially expressed genes in the 2 cuproptosis subtypes, with logFCfliter = 0.585 and adjusted p-value fliter = 0.05. We identified 2,486 differentially expressed genes, and we believe that these genes are related to potential biological differences between the cuproptosis subtypes. A volcano plot showed all differentially expressed genes and highlighted the cuproptosis-related genes (Supplementary Figure S2). Functional enrichment analysis was be conducted by a gene ontology database (Figure 4A). And those cuproptosis subtype-related genes were enriched for cytokine production, cell cycle, cell adhesion, and extracellular environment. An analysis of KEGG revealed the enrichment of proteoglycans in cancer, cell cycle, senescence, cell adhesion, and p53 signaling (Figure 4B). We could find that GO and KEGG are similar in the following ways: cell cycle, cell adhesion and extracellular environment (Figures 4A, B). We also believed that GO and KEGG confirmed GSVA results (Figure 3E), which shows different pathways between two cuproptosis subtypes such as KEGG_CELL_CYCLE, KEGG_PROTEIN_EXPORT and KEGG_MISMATCH_REPAIR. In addition, it was revealed that cuproptosis subtype-related genes may influence glioma cell characteristics and TME, and may contribute to glioma oncogenesis. To learn more about the significance of cuproptosis subtype-related genes in glioma prognosis, we used a consensus clustering algorithm based on the PAM algorithm. The 725 cohort patients were divided into 2 new subtypes: gene subtypes A and B with 459 and 266 patients, respectively (Figure 4C; Supplementary Figure S2). Kaplan-Meier analysis illustrated a significantly difference in prognosis between the 2 gene subtypes (p < 0.001, Figure 4D). The heatmap of gene subtypes showed that most elderly patients and patients with grade 4 glioma were classified as gene subtype B, explaining the prognostic differences in cuproptosis gene subtypes (Figure 4E). Significant differences were also found in cuproptosis-related gene expression between the 2 cuproptosis gene subtypes (Figure 4F).

The CRG-score was constructed using differentially expressed genes among the gene subtypes. Using "caret" package, we randomly distributed the cohort patients (4 patients were excluded due to incomplete survival and gene expression data) into training (n = 360) and test (n = 361) groups. Our next step was to select the appropriate prognostic genes related to the gene subtypes through LASSO analysis and multivariate Cox analysis. Following result of the former with the minimum partial likelihood deviance, we kept 34 genes (Supplementary Figures S3A, B). And result of the latter revealed 11 genes, including 2 low-risk genes (NOG and MKX) and 9 high-risk genes (NBPF8, TSKU, AURKB, SLC25A43, P2RY6, STEAP1, CDK4, RARRES1, and KCNN4). The CRG-score was established as follows using the multivariate Cox regression analysis (Supplementary Table S5). The coefficients were kept 4 decimal places. C R G − s c o r e = − 0.2885   *   N O G + − 0.2045   *   M K X + 0.2943   *   N B P F 8 + 0.3126 *   T S K U + 0.1503   *   A U R K B + 0.2898   *   S L C 25 A 43 + 0.2731   *   P 2 R Y 6 + 0.1496   *   S T E A P 1 + 0.1252   *   C D K 4 + 0.1361   *   R A R R E S 1 + 0.1671   *   K C N N 4

The CRG-score showed a significant difference between gene subtypes (Figure 5A). Patients with subtype A generally had a lower CRG-score, which means that these patients had a lower risk and better prognosis. Similar differences were found between the 2 cuproptosis subtypes (Figure 5B). This indicates that the CRG-score may have relevance to immune infiltration, metabolism, and cell cycle.

Compared with the median risk score, patients with a lower CRG-score were divided as the low-risk group (n = 356), and a higher score were divided as the high-risk group (n = 365) (Figure 5C). According to the distribution plot of the CRG-score risk, death rate from glioma rises while survival times decrease with mounting risk score (Figure 5D). The heatmap of the CRG-score-related gene expression in the 2 risk groups was also consistent with the formula of the CRG-score (Figure 5E). PCA analysis showed the dimensions of the low-risk and high-risk group (Figure 5F). Kaplan-Meier analysis indicated that patients in the high-risk group had a significantly poor prognosis (Figure 5G). The expression of cuproptosis-related genes varied widely, which confirmed the prospective connection between the CRG-score and cuproptosis (Figure 5H). Sankey diagram illustrates entire network of the subtypes, risk of CRG-score, and survival status (Figure 5I). The ROC curve represented the survival rate of the CRG-score with area under the curve (AUC) values of 0.874, 0.924, and 0.875, which contain respectively the 1-, 3-, and 5- year (Figure 5J).

To validate the prognostic performance of the CRG-score, we computed the CRG-score of an external validation group (GSE83300) and the test group which was randomly chosen from the cohort by “caret” packages. The distribution plot, heatmap, Kaplan-Meier analysis, and ROC curve analysis are shown in Supplementary Figures S5, S6. Consistent with the previous results, the low-risk group had better prognosis. The AUC values at the 1-, 3-, and 5-year were relatively high, which indicates that the CRG-score had a satisfactory capability to predict the survival of patients with glioma patients.

We used RT-qPCR to analyze the expression of 11 CRG-score-related genes in gliomas and adjacent tissues (n = 5). The results are shown in Supplementary Figure S7.

With the CIBERSORT algorithm, we further explored the potential role of immune infiltration by evaluating the association between the CRG-score and the abundance of immune cells. The CRG-score was positively correlated with Tregs, CD8⁺ T-cell, follicular helper T-cell, resting NK cells, neutrophils, M0, M1 and M2 macrophages. It presented negative correlation with resting memory CD4⁺ T-cell, plasma cells, activated NK cells, monocytes, activated mast cells, and eosinophils (Figure 6A). And a vioplot used Wilcoxon test was more significantly showed the relationship between immune infiltration and CRG-score (Figure 6B). Further, the CRG-score model genes were also correlated to the abundance of most immune cells significantly (Figure 6C). Using “estimate” package, we evaluated the TME score of low- and high-risk group (Figure 6D). A high CRG-score was associated with a high stromal and immune score, which indicated a close relationship between the CRG-score and TME.

The mutation data of the TCGA-LGG and GBM cohorts revealed a higher TMB index in the high-risk group (Figure 7A). Spearman correlation analysis indicated that the CRG-score was positively correlated with TMB (Figure 7B). We analyzed the somatic mutations in low- and high-risk group (Figures 7C, D). Differences in classical genes related to oncogenesis and development in gliomas, such as IDH1, TP53, ATRX, PTEN, EGFR, CIC, and PIK3CA, were observed. The low-risk group had higher mutation frequency of IDH1, TP53, ATRX, and CIC, in contrast to the high-risk group that had higher mutation frequency of PTEN, EGFR, MUC16, and PIK3CA. The somatic mutation results were consistent with the existing researches about prognosis of the gliomas (Rasheed et al., 1997; Cancer Genome Atlas Research et al., 2015; Chen et al., 2019; Bjorkblom et al., 2022; Ferrer, 2023). We also found that there was negative correlation between the CRG-score and the cancer stem cell index, indicating that gliomas with higher CRG-score had less stem cell characteristics and a higher degree of cell differentiation (Figure 7E).

Based on the CRG-score and clinicopathological characteristics, 1-, 3-, and 5-year survival could be predicted by a nomogram, which was constructed by “rms” package (Figure 8A). The calibration curve suggested that the nomogram was competent enough (Figure 8B).

83 glioma tissue samples, including 27 grade 2 samples, 6 grade 3 samples, 50 grade 4 samples, and 12 adjacent tissues, were made into tissue microarrays. FDX1, DLST and PDHA1 were chose for immunohistochemistry. FDX1 and lipoylation act as key regulators in cuproptosis. DLST is one of the enzymes where lipoylation occur. PDHA1 was one of the protein targets of lipoylation. So, we believed that these three indicators are representative in cuproptosis. The expression of FDX1, DLST, and PDHA1 was detected and graded (Supplementary Figure S8). The results showed a general difference between LGG, GBM, and adjacent tissues, which is consistent with the database analysis (Figure 9).