Identification of a copper metabolism‐related gene signature for predicting prognosis and immune response in glioma

Abstract Background Gliomas are highly refractory intracranial cancers characterized by genetic and transcriptional heterogeneity. However, therapeutic options are limited. In the last years, copper‐induced cell death is becoming a prospective treatment strategy for gliomas and other solid tumors, but copper metabolism‐related genes associated with cancer development remain unclear. Methods We first collected gene expression data from The Cancer Genome Atlas (TCGA) to identify significantly differentially expressed copper metabolism‐related genes in gliomas. Using these genes, we performed COX regression and Least Absolute Shrinkage and Selection Operator (LASSO) regression to construct the prognostic model. The prognostic value of the model was further validated by CGGA testing set. Subsequently, functional analyses were carried out, including gene set enrichment analysis (GSEA), immune infiltration analysis, and mutation analysis. Finally, the expression levels of these genes were verified by immunohistochemical analysis. Results The prognostic model consisted of 7 genes: CDK1, LOXL2, LOXL3, NFE2L2, SLC31A1, SUMF1 and FDX1. According to this prognosis model, glioma patients could be split into the high‐risk group or low‐risk group, and the low‐risk group showed significantly better prognostic survival (p < 0.001). Moreover, the high‐risk group had higher levels of immune cell infiltration, immune checkpoint genes expression, and higher tumor mutational burden (TMB), which indicates that they might benefit more from immunotherapy. Finally, we confirmed the expression level of FDX1, SUMF1, and SLC31A1 protein as significantly different in glioblastoma, lower‐grade glioma, and non‐tumor brain tissues by immunohistochemical analysis, and the high expression of FDX1 and SLC31A1 protein was related to poor survival in glioma patients. Conclusions Our study could contribute to the prognosis prediction and decision‐making in patients with gliomas.


| INTRODUCTION
Gliomas are highly refractory and the most prevalent primary intracranial cancers, comprising 80% of all malignancies in the brain. 1 Diffuse gliomas consist of lowergrade gliomas (LGGs), including WHO grade II and III tumors, and glioblastoma (GBM). 2 The median OS time of patients with diffuse gliomas is 78.1 months for WHO grade II, 37.6 months for WHO grade III, and 14.4 months for WHO grade IV gliomas. 3 The main treatment of glioma is the combination of surgery and radiotherapy or chemotherapy. However, most of the patients experience limited benefit from them. 4 Thus, diffuse glioma is still a heavy burden for human health, and new therapies are needed to improve the prognosis of patients.
An important reason hindering the development of novel treatment for gliomas is that gliomas are highly heterogeneous in all molecular subgroups. 5 Histology alone is insufficient to give the accurate diagnosis and prognostic stratification needed for targeted treatment. Combination of histological assessment with molecular markers facilitates more accurate patient diagnosis and risk stratification. For example, isocitrate dehydrogenase 1 (IDH1) mutation has been identified as a typical favorable prognostic marker. 6 However, more effective prognostic molecular markers are particularly needed for improving prognosis and long-term life quality of patients with diffuse gliomas.
In detecting molecular markers related to prognosis of cancers, recent studies have indicated the importance of copper metabolism. 7 Copper, as an essential cofactor for enzymes and a cellular signaling agent, mediates a range of biochemical processes, including mitochondrial respiration, antioxidation, and hormone and neuropeptide biogenesis. However, excessive copper ion induces cytotoxicity; therefore, intracellular copper is restricted to low concentrations by homeostatic mechanisms to prevent the accumulation of free copper and copper-induced cytotoxicity. [8][9][10] Unlike normal tissues, it was proved that tumor cells have a relatively higher concentration of copper, which may contribute to the pathogenesis of it. [11][12][13][14][15][16][17] Tsvetkov et al. found abnormal elevated copper ions may cause a previously unidentified way of cell death, which induced by aggregation of lipoylated TCA cycle proteins and termed this distinct copper-dependent cell death mechanism "cuproptosis." 7 Copper-induced cell death is becoming a prospective treatment strategy for glioma and other solid tumors. [18][19][20][21][22][23][24] Studies showed that copper oxide-related cytotoxicity can cause severe cell death in C6 glioma cells, which can be prevented by adding copper chelators. 25 Zheng et al. showed a potential combinatory therapeutic strategy of gliomas, which contains copper. 26 Moreover, Tsvetkov et al. identified several genes involved in cuproptosis, whose altered expressions are highly correlated across human tumors, making tumors more sensitive to cuproptosis. 7 Copper metabolism-related genes thus emerge as promising therapeutic targets and potential prognostic markers for gliomas. However, there are still insufficient studies in constructing prognostic models in diffuse glioma based on copper metabolism, which hindered the clinical application of the copper metabolism-related genes in diagnosis and therapy design in glioma patients.
In this study, first, we performed a bioinformatic analysis based on mRNA expression and survival data from TCGA, and constructed a prognostic model for glioma patients using copper metabolism-related genes. Then we validated the model in the CGGA dataset. 27 Finally, GSEA, immune cell infiltration analysis, and mutation landscape analysis were performed to investigate the underlying mechanisms of the signature in the immune microenvironment and its guiding value for immunotherapy of glioma. Our study could be helpful for the prognosis prediction and decision-making in patients with gliomas.

| Collection of copper metabolism-related genes
Copper metabolism-related genes were downloaded from the Molecular Signatures Database (MSigDB) version 7.5.1. 28,29 We also added the previously reported regulatory gene of copper-induced cell death. 7 Thus, a set of 111 genes related to copper metabolism was obtained (Table S1).  31 The values of gene expression profiles in these datasets were converted by log 2 for subsequent analysis.

| Differential analysis
Differential analysis was performed between the TCGA glioma samples and normal brain samples by "DESeq2" package in RStudio (version 4.1.3). Differentially expressed genes (DEGs) were filtered by a false discovery rate (FDR) threshold <0.05 and absolute log2-foldchange>1. From the DEGs, we selected the copper metabolism-related genes for further analysis.

| Building and validating the prognostic model
LASSO regression analysis was applied to construct a prognostic model, which has reduced the overfitting by selecting high-dimensional prognostic genes using the "glmnet" package. The risk score was calculated using the following formula: where expr genei is the normalized gene expression value and coefficient genei is its coefficient.
Through this formula, each LGG/GBM patient can be calculated with a risk score, and patients were grouped into the high-risk or low-risk group according to the median risk score. (Detailed procedures were described in the Supplementary Methods.)

| The construction of a nomogram
To better predict the prognosis in both the training and testing cohorts, a nomogram was constructed using the "rms" package. Calibration analysis and ROC curve were used to assess the accuracy of the nomogram. (Detailed procedures were described in the Supplementary Methods.)

| GSEA
The "DESeq2" package was used to obtain the DEGs between the high-and low-risk groups in the TCGA set. Then, GSEA 32 was performed using the "clusterProfiler" package. 33

| Immune cell infiltration analysis
We downloaded the immune infiltration estimation results for GBM/LGG patients in the TCGA database from the Timer2.0 database. 34

| Mutation landscape
Mutational data of patients (n = 681) with LGG/GBM from TCGA were obtained from cBioPortal. 35 The mutational data was analyzed by "maftools" R package 36 (Supplementary Methods).

| Protein expression in HPA database
The Human Protein Atlas (HPA) database (www.prote inatl as.org) was used for comparison of protein expression of prognostic genes between LGG and GBM samples by immunohistochemistry image. 37

| Tissue specimens and patient information
From Department of Neurosurgery, Sun Yat-sen Memorial hospital, we collected 91 surgical specimens, which had been histopathologically diagnosed gliomas by Department of Pathology of our hospital according to WHO criteria between 2012 and 2022. Among these, 60 surgical specimens with follow-up survival data were made into tissue microarrays (TMAs), while the rest (n = 31) were tissue slides. These samples include 55 GBM specimens, 31 LGG specimens and 5 normal brain tissues (from surgical specimens of patients with cerebrovascular disease) (Supplementary Methods).

| Immunohistochemical (IHC) staining
Immunohistochemical analysis was performed using tissue slides or TMAs of human surgical specimens (Supplementary Methods). Figure S1 shows the flowchart of this study. The DEG analysis was performed between TCGA glioma samples

| Identification of prognostic copper metabolism-related DEGs in the TCGA-LGG/GBM cohort
expr genei × coefficient genei and normal brain samples, among them, 39 were identified as copper metabolism-related DEGs. Then, univariate Cox regression analysis of OS was performed on these DEGs, indicating 34 of them were correlated with OS. From these genes, 18 genes were excluded for further analysis because 8 of them was upregulated in tumor samples with their expression levels positively correlated with survival rate, and 10 genes was downregulated in tumor samples with their expression levels negatively correlated with survival rate. Of the remaining 16 copper metabolism-related DEGs, 12 were upregulated, and 4 were downregulated in glioma compared to normal tissue expression ( Figure 1A). Functional enrichment analyses indicated these genes were significantly (Pcorrection <0.01) enrichment in cellular response to copper ion and copper ion binding ( Figure 1B).

| Construction and validation of a copper metabolism-related prognostic model
In order to further refine selected genes mentioned above and construct a prognostic model, we performed LASSO regression analysis of the 15 genes mentioned above (AQP1 was not found in CGGA database, so it was excluded). The results showed that the model had the best fitting effect when the number of genes was 7 and λ was 0.0889 ( Figures 1C,D). The results of univariate Cox regression analysis and Lasso regression coefficient of these seven genes are summarized in Figure 1E, which were used for constructing the prognostic model.
According to the expression levels and regression coefficients of selected genes, the risk scores were calculated for each cases in both sets. Figure 2A shows the differences between the high-risk and low-risk groups in the living state, survival times and the distribution of the prognostic genes expression in TCGA sets, which indicated that patient with a higher risk score may be correlated with a poorer outcome in the TCGA training cohort. All these seven genes displayed higher expression in the high-risk group than in the low-risk group. The similar gene expression profiles were observed in the testing set ( Figure 2D).
To further evaluate the predictive performance of the gene signature in predicting the prognosis of glioma patients, we performed ROC curve analysis. The results showed that the area under curve (AUC) values of 1-year reached 0.858, 3-year reached 0.894, and 5-year reached 0.823 ( Figure 2B), and the OS time of the high-risk group was significantly shorter than that of the low-risk group in the TCGA set. (p < 0.001, Figure 2C). In the CGGA dataset, we found that the AUC of 1, 3 and 5 years were 0.713, 0.788 and 0.808, respectively ( Figure 2E). Similarly, in the CGGA dataset, relative to patients in the low-risk group, patients in the high-risk group have a significantly worse prognosis (p < 0.001, Figure 2F). These results showed that the gene signature has good performance in predicting the prognosis of patients in both groups.

| Clinical stratified analysis
To investigate the correlation between risk score and other clinical characteristics, stratified analysis of risk score was performed according to different clinical characteristics. As shown in Figures S2, in the TCGA set, patients with higher risk score tend to be older than 40 years, to have higher WHO grade, unmethylated MGMT promoter, glioblastoma, non-codeletion of 1p/19q and wild type IDH1 (all p < 0.001). In the CGGA set, stratified analyses of risk score for WHO grade, histological type, 1p/19q codeletion, and IDH1 mutation status were consistent with those in TCGA database (all p < 0.001, Figures S2). No significant difference was found in risk score distribution between different age levels or MGMT methylation status (Figures S2). No correlation was observed between risk score and gender in both sets (data not shown). Thus, these findings indicated that risk scores constructed by the seven-gene signature are able to reflect the clinical status of patients in both TCGA and CGGA sets.

| Construction and validation of the nomogram
To further investigate the independence of the risk score as a risk factor for OS time in diffuse glioma patients, univariate and multivariate Cox regression analysis were performed by considering potential risk factors as parameters. As shown in Table 1, in multivariate Cox regression, the risk group was identified as independent risk factor (HR = 2.257, 95% CI = 1.413-3.606, p < 0.001). Besides, the WHO grade, age, IDH status, 1p/19q codeletion status were also considered to be independent risk factors related to prognosis, and were used to construct the nomogram.
To better evaluate the risk of glioma patients, a nomogram was constructed by combining the prognostic model and other independent risk factors. The C-index of the nomogram model was 0.847 (95% CI = 0.837-0.858). As shown in Figure 3A, each item in the nomogram was assigned a point according to the actual condition of the case. The total points were calculated for each case and used to predict the survival rates. Calibration analysis and ROC curve analysis were performed to assess the prediction performance of the nomogram. In TCGA set, the AUC reached 0.877 at 1 year, 0.937 at 3 years, and 0.879 at 5 years OS rates of nomogram ( Figure 3B). The calibration curves showed good fitting between the observed value and optimized value of 1-, 3-, and 5-year OS rates ( Figure 3D-F). In addition, the accuracy of our nomogram was also validated in the CGGA set. The AUC reached 0.821 at 1 year, 0.868 at 3 years and 0.878 at 5 years OS rates of nomogram ( Figure 3C). Calibration curves also showed good fitting between the observed and optimized survival rate in CGGA set ( Figure 3G-I). These results suggest that the nomogram has good predictive ability in both sets and can be used as a reference for clinical decision-making.

| GSEA
According to the above analyses, patient prognosis differed significantly between two risk groups. To explore the reasons for this difference, we performed GSEA on genes expressed significantly differently between two risk groups in the TCGA set. GSEA results suggested that genes highly expressed in the high-risk group were significantly enriched in G2M checkpoint, PD-1 signaling, epithelial mesenchymal transition, complement, IL-10 signaling, PI3KAKT signaling pathways ( Figure 4A-F). As the PD-1 signaling pathway and IL-10 signaling pathway are key pathways for immunosuppression, 38,39 these results suggest that immunosuppression and innate immune dysregulation may play an important role in the progression of diffuse glioma, which further suggests that diffuse glioma is a kind of tumor that highly associated with immune dysfunction.

| Immune infiltration analysis and mutation landscape
Based on the results of GSEA, the immunosuppressive tumor microenvironment (TME) may play a role in the progression of diffuse glioma, so the comparison of immune infiltration levels between two groups in the TCGA dataset was performed. As shown in Figure 5A, in the high-risk group, there were more immune cells infiltration, including CD8+ T cells, regulatory T cells, neutrophils, macrophages, myeloid dendritic cells, and higher stroma in tumor tissue. Next, we explored the expression of immune checkpoint genes. As shown in Figure 5B, the expression levels of most immune checkpoint genes, such as CD274, PDCD1, LAG3, CTLA4, and HAVCR2, are higher in the high-risk groups (p < 0.0001). The high expression of immune checkpoint genes may result in suppressed immune cell function and blunted immune response. Therefore, although the high-risk group was related to a higher level of immune cell infiltration, the tumor cells remodel the TME by regulating the immune checkpoint genes expression of tumor cells and immune cells, thus inhibit antitumor immunity and lead to tumor immune evasion, which can be one of the explainations for the different prognosis between two risk groups.
To further explore the potential role of the prognostic model in predicting the response of immunotherapy, we analyzed the mutation landscape in the TCGA dataset. In the low-risk group, the mutation incidence was 94.96% and the most significantly mutated gene was IDH1 ( Figure 5C). While the mutation incidence in the high-risk group was 92.64%, and the most significantly mutated gene was TP53 ( Figure 5D). Figure S3 shows the overall mutated genes with significantly different mutation incidence between the two groups. Furthermore, as TMB was found to be a predictive biomarker for the clinical efficacy to immunotherapy, 40 we calculated the correlation between risk score and TMB for patients in the TCGA-LGG/GBM samples. The result showed that TMB was positively correlated with risk score (R = 0.51, p < 0.0001, Figure 5E). In addition, the TMB was significantly higher in the high-risk group than in the low-risk group (p < 0.001; Figure 5F). Therefore, these results indicated that patients in the high-risk group may benefit more from immunotherapy.

| Aberrant expression of the identified genes in glioma samples detected by IHC and survival analysis
Subsequently, we verified the protein expression level of these identified genes by immunohistochemical analysis. First, we searched these genes in HPA database. Among them, LOXL3 and FDX1 were not detected in LGG or GBM samples. NFE2L2 protein highly expressed in LGG while moderately expressed in GBM. CDK1 and LOXL2 had low expression in GBM while were not detected in LGG ( Figure S4). Since the HR values of FDX1, SUMF1, SLC31A1 were the highest in the seven-gene signature, we were going to validate their protein expression in human glioma specimens. The aberrant protein expressions of these three genes in glioma samples were confirmed by IHC ( Figure 6A,B). FDX1 shows higher protein expression in GBM and LGG relative to normal brain tissue (Both p < 0.001). Notably, the expression of FDX1 protein in GBM is significantly higher than that in LGG (p < 0.05). The expression of SLC31A1 protein in GBM was significantly higher than that in LGG (p < 0.001). The expression of SUMF1 protein in GBM and LGG was significantly higher than that in normal brain tissue (Both p < 0.001).
As for survival analysis, the high expression of FDX1 and SLC31A1 protein was associated with poor survival in glioma patients (p < 0.05 for FDX1 and p < 0.001 for SLC31A1, Figure 6C). However, there was no significant correlation between the protein expression of SUMF1 and OS time in patients with gliomas.

| DISCUSSION
Gliomas are highly heterogeneous and refractory tumors with high recurrence rates and resistance to chemo-and radiotherapy. 41 Copper is an important nutrient involved in many biological functions; however, due to its redox properties, abnormal copper levels can be toxic to cells. The connection between copper and copper-dependent diseases, especially cancer, has been noted, since numerous observations have shown that tumors have a greater demand for copper compared to healthy tissue. Furthermore, elevated copper levels have been proven to be directly related to cancer progression. 42 Copper metabolism was reported to take part in glioma proliferation, angiogenesis, 43 and the TME. 26 Copper metabolismrelated genes are abnormally expressed in tumor tissues relative to normal tissues, 44 and thus emerge as potential prognostic markers for gliomas. However, there are still insufficient studies in exploring the relationship between copper metabolism and clinical characteristics among glioma patients.
In the current study, we conducted a bioinformatic analysis using gene expression and survival data of from TCGA, and screened for 16 DEGs with prognostic value. Then, seven genes were screened to construct a prognostic model by LASSO regression analysis. The prognostic model consisted of seven genes: CDK1, LOXL2, LOXL3, NFE2L2, SLC31A1, SUMF1, and FDX1. CDK1, a kinase related to the cell cycle, was found to be up-regulated by temozolomide (TMZ) in a NF-κB dependent way in GBM, 45 and can be a target for glioma treatment by affecting the cell cycle. 46 LOXL2 and LOXL3, members of the lysyl oxidase (LOX) family, LOXL2 was found to promote glioma progression, and improve the TMZ resistance of glioma cells. 47 A study found that downregulation of LOXL3 decreased cell proliferation and invasion, and increased cell adhesion in U87 cells. 48 NFE2L2 is a transcription factor which protects cells from the oxidative damage and is associated with resistance to cancer treatments, 49 which was indicated as a prognostic indicator in gliomas. 50 SLC31A1 encodes copper transporter 1 (CTR-1), which is a transmembrane pump responsible for copper absorption. It is highly expressed in different types of cancer characterized by elevated copper levels, and high CTR-1 expression enhances PD-L1 expression by promoting cellular copper influx, aggravating PD-L1-driven tumor immune escape. 51 SUMF1, has been reported whose increased expression could improve relapse-free survival of breast cancer patients. 52 FDX1, which encodes a reductase that reduce Cu 2+ to its more toxic form, Cu 1+ , is a direct target of elesclomol (a copperloaded ionophore). 53 Peter Tsvetkov found through knockout studies that FDX1 deletion confers resistance to copper-induced cell death. Since cells with high levels of FDX1 are sensitive to copper-induced cell death, the authors suggest that copper ionophore treatment should be targeted at tumors with this metabolic profile.
Based on GSEA, genes highly expressed in the highrisk group were enriched in pathways related to tumor progression and immune response. It is reasonable to assume that copper metabolism may be closely linked to tumor progression and TME in glioma. Kenichi reported that copper ionophore NSC319726 induced DNA damage, causing cell cycle arrest in glioblastoma cells. However, it showed no substantial toxicity to differentiated and ungrown neurons, and this selective lethality of copper ionophore is an ideal feature as an anti-GBM drug. 54 PD-1 is an immunosuppressive receptor that causes T-cell dysfunction and apoptosis by binding to its ligand PD-L1 55 and plays a role in braking inflammatory responses and conspiring immune evasion of glioblastoma cells. 56 Epithelial-mesenchymal transition (EMT) is the main cause of recurrence and poor prognosis of glioma. 57 IL-10 is an anti-inflammatory cytokine produced mainly by M2 macrophages. 58 Vidhya M found that IL-10 release drives T-cell exhaustion in glioblastoma, thereby contributing to an immunosuppressive TME. 59 The PI3K/Akt/mTOR signaling pathway is one of the key pathways in tumor progression, and its molecular alteration is a typical hallmark of glioma. 60,61 We further explored the heterogeneity of the tumor immune environment between two groups to look for potential causes of different outcomes. The risk score was positively correlated with the level of immune cell infiltration and stromal infiltration, which was consistent with previous study indicating that tumor-surrounding stromal cells, such as tumor-associated macrophages (TAMs), myeloid-derived suppressor cells (MDSCs), and regulatory T cells (Tregs), could produce chemokines to promote tumor growth, resulting in poor prognosis in patients. 62 In addition, the higher expression levels of immune checkpoint proteins in patients in the high-risk group, where the more immunosuppressive TME predominated, provides a plausible explanation for their poorer survival outcome.
TMB is an emerging potential biomarker for immunotherapy. 63 A prevailing view in the field of immunotherapy is that tumors with increased TMB present more neoantigens and therefore are correlated with better clinical response to immunotherapy. 64 One of the major obstacles to effective immunotherapy for GBM is the low TMB in most GBM tumors. 65 Our study showed a positive correlation between risk scores and TMB in gliomas, which may indicate that the patient with higher risk score respond better to immunotherapy.
As cuproptosis was identified to be a novel mechanism of cell death, the role of copper metabolism in tumor progression and its guiding value for treatment have attracted much attention. Some studies have proved that copper ionophores may contribute to the treatment of glioma, 54 but no study has yet explored the value of the copper metabolism-related genes in predicting the prognosis and guiding the treatment of glioma. In this study, copper metabolism-related gene prognostic model for glioma was developed, and the clinical features of glioma patients were further integrated to develop a nomogram to better predict the survival of glioma patients. Compared with other published models, our model has better prediction efficiency (Table S2). In addition, we explored the potential mechanisms of different outcomes of glioma patients from the perspective of tumor immune microenvironment, and predicted the possible immunotherapy response of glioma patients. We further confirmed the abnormal expression of FDX1, SUMF1, and SLC31A1 proteins in glioma samples by immunohistochemistry, and the correlation between high expression of FDX1 protein and poor prognosis in glioma patients. Considering that FDX1 is a key regulator mediating cuproptosis, our study also provided a new perspective and direction for the future research on the treatment of glioma using the cuproptosis mechanism.
However, this study also has some limitations. First, immune microenvironment analysis and mutation analysis were not validated in the testing set due to lack of data. Second, in this study, only partial genes were verified for protein expression and survival analysis in clinical samples, without further experiments in vitro or in vivo to explore their roles and mechanisms in glioma.

| CONCLUSION
In conclusion, we constructed a copper metabolismrelated gene prognostic model in glioma for the first time. Using this model, we can well predict the prognosis and the possible immunotherapy response of glioma patients. We further confirmed the abnormal expression of FDX1, SUMF1, and SLC31A1, and found the correlation between high expression of FDX1 and SLC31A1 protein and poor prognosis in glioma patients, which provided a new direction for the future research on the treatment of glioma using the cuproptosis mechanism.