Cuproptosis-related gene expression is associated with immune infiltration and CD47/CD24 expression in glioblastoma, and a risk score based on these genes can predict the survival and prognosis of patients

Introduction Glioblastoma (GBM) is the most invasive type of glioma, is insensitive to radiotherapy and chemotherapy, and has high proliferation and invasive ability, with a 5-year survival rate of <5%. Cuproptosis-related genes (CRGs) have been successfully used to predict the prognosis of many types of tumors. However, the relationship between cuproptosis and GBM remains unclear. Methods Here, we sought to identify CRGs in GBM and elucidate their role in the tumor immune microenvironment and prognosis. To that aim, changes in CRGs in The Cancer Genome Atlas (TCGA) transcriptional and Gene Expression Omnibus (GEO) datasets (GEO4290 and GEO15824) were characterized, and the expression patterns of these genes were analyzed. Results A risk score based on CRG expression characteristics could predict the survival and prognosis of patients with GBM and was significantly associated with immune infiltration levels and the expression of CD47 and CD24, which are immune checkpoints of the “don't eat me “signal. Furthermore, we found that the CDKN2A gene may predict GBM sensitivity and resistance to drugs. Discussion Our findings suggest that CRGs play a crucial role in GBM outcomes and provide new insights into CRG-related target drugs/molecules for cancer prevention and treatment.


Introduction
Glioma is the most common malignant primary brain tumor, with glioblastoma (GBM) being the most aggressive type, constituting >50% of all gliomas, and having an incidence rate of 3-5 per 100,000 population (1,2). GBM shows the worst prognosis, with a median age at diagnosis of 65 years (3). Those under 70 years of age who do not undergo treatment have a median survival of approximately 3-4.5 months (4). Precision medicine, which combines molecular biomarkers and targeted therapies, has become increasingly important in modern cancer treatment (5). Therefore, it is particularly important to screen for molecular markers and target genes of GBM.
Copper (Cu) is an important cofactor in all organisms, but if the concentration exceeds the threshold maintained by the evolutionarily conservative homeostasis mechanism, it can induce a form of cell death named cuproptosis (6). Cu ions in mitochondria directly bind to the fatty acylated components of the tricarboxylic acid (TCA) cycle, resulting in the abnormal aggregation of fatty acylated proteins and loss of iron thiocluster proteins, which causes protein toxic stress reactions and eventually cell death (7). Cuproptosis-related genes (CRGs) [false discovery rate (FDR) < 0.01] include FDX1, LIAS, LIPT1, DLD, DLAT, PDHA1, PDHB, MTF1, GLS, and CDKN2A. CRGs have aroused interest in studying the regulation of mitochondrial copper homeostasis during normal and physiological changes and as a target for cancer therapy. CRGs successfully predict the prognosis of soft tissue sarcoma (8), hepatocellular carcinoma (9, 10), melanoma (11), and renal clear cell carcinoma (12,13). However, the role of CRGs in GBM remains unclear. Therefore, the aim of our study was to identify CRGs in GBM and elucidate their role in the tumor immune microenvironment and prognosis. This will provide insights not only into the signaling pathways and molecular mechanisms of cuproptosis in GBM, but also into the effects of immunotherapy on patients with GBM.

Data collection
The sequencing data and clinical characteristics of GBM tissues were obtained from UCSC XENA (https://xenabrowser.net/ datapages/, accessed on 7 July 2022) (14) through the Toil process, which uniformly processes RNAseq data in transcripts per million reads (TPM) format of The Cancer Genome Atlas (TCGA) and the Genotype-Tissue Expression (GTEx) database.
The corresponding normal tissue data were also extracted. A total of 1,846 samples were obtained, including datasets of 689 tumor samples, five paracancerous datasets, and 1,152 normal tissue samples. The GSE4290 (15) and GSE15824 (16) datasets were derived from the Gene Expression Omnibus (GEO) (https:// www.ncbi.nlm.nih.gov/geo/, accessed on 12 July 2022). The ggplot2 R package was used for visualization.

Sample collection
The human glioblastoma cell lines A172 were purchased from the Cell Bank of the Servicebio (Wuhan, China). A172 were cultured in RPMI-1640 medium (Procell, China), with 10% fetal bovine serum (Gibco, USA). This cell were cultured at 37°C with 5% CO 2 . Three paired GBM and adjacent non-tumor tissues were collected in Xijing Hospital; written informed consent was obtained from the patients. All these GBM patients did not receive chemotherapy or radiotherapy prior to surgery. All tissue samples were snap-frozen and stored in liquid nitrogen (−80°C) until RNA extraction. The study was approved by Xijing Hospital of Fourth Military Medical University (KY20193098).

Copper assay
Tissue (0.1 g) was removed and irrigated in cold saline, and wiped dry of surface water; 0.86% saline was added according to the ratio of weight (g):volume (ml) = 1:9. It was then crushed for 1 min with an ultrasonic grinder (Ningbo Xinyi YJ92-IIN) and subjected to lowtemperature and low-speed centrifugation at approximately 2,500 rpm for 15 min; the supernatant was taken for analysis. Copper content was measured using a Copper Assay Kit (Nanjing Jiancheng Bioengineering Institute), according to the manufacturer's instructions. The absorbance was determined at 600 nm wavelength using a microplate reader (TECAN Infinite M200 Pro).

Plasmids, shRNAs, and transfection
Short hairpin RNA (shRNA) oligonucleotide sequences were designed and synthesized by Shanghai Genechem Company (Shanghai, China). The shRNA sequences for knockdown of CD24 and CD47 is shown in Table 1, and A172 was transfected with shRNA using Lipofectamine 3000 (Invitrogen, USA) according to the manufacturer's instructions. Cells were then used for assays 48 h post-transfection.

Quantitative PCR
Total RNA was extracted using the GeneJET RNA Purification Kit (Thermo Scientific, USA), and cDNA was acquired by reverse transcription using the PrimeScript ™ RT reagent kit (DIYI, China). Quantitative real-time PCR was performed with a TB Green Fast qPCR Mix (DIYI, China). The 2 −DDCt method was used to quantify the expression of each gene normalized to that of actin. Detailed information on the primer sequences for each gene is shown in Table 2.

Consensus clustering analysis of CRGs
Ten CRGs were genotyped in TCGA and GTEx data. The ConsensusClusterPlus R packet (v1.54.0) was used for consistency analysis; the maximum number of clusters was six and 80% of the total samples were drawn 100 times. A cluster heat map was created in R using the package ggplot2 (version 3.3.3). The statistical method used was Spearman's rank correlation coefficient. The cBioPortal database (http://www.cbioportal.org/, accessed on 19 July 2022) was used to obtain the CRG mutation profile in GBM.

Gene network and enrichment analysis of CRGs
To analyze the potential interactions of these genes, we performed gene network analysis using STRING (17). Furthermore, we performed pathway enrichment analysis of the CRGs using the Database for Annotation, Visualization, and Integrated Discovery (18). Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) were used as references, and enrichment analysis was performed using the R package "clusterProfiler." We applied the Benjamini-Hochberg method for multiple correction and an FDR < 0.05 was considered significant.

Differential expression analysis and validation
To verify the differential expression level of CRGs in GBM and normal tissues, we collected the data of GSE4290 (GPL570) and GSE15824 (GPL570) and used box charts and the R packet "ggplot2" to compare the expression of CRGs in different datasets. After calculating the change in log2 multiple and 95% confidence intervals, we conducted a meta-analysis of the results of differential expression to improve the statistical ability of our research.

Analysis of correlation with immune infiltration
The Tumor Immune Estimation Resource (TIMER2.0; http:// timer.cistrome.org/) (19) was used to investigate the relationship between the expression of CRGs and the abundance of six immune cells (B cells, CD8 + T cells, CD4 + T cells, neutrophils, macrophages, and dendritic cells). We also examined two immune checkpoints involved in the "don't eat me" signaling pathway, including CD24 and CD47, because the expression levels of immune checkpointrelated genes are related to the therapeutic response to immune checkpoint inhibitors (20). Pearson's correlation analysis was used to test the relationship between GBM and CRGs.

Prognostic significance of CRGs in GBM
We used the pROC package (21) to transform the RNA-seq data in TPM format into log2 to compare the expression between samples. The area under the curve (AUC) values were determined by performing receiver operating characteristic (ROC) analysis, and the predictive capabilities of the signature were assessed. The abscissa is the false positive rate (FPR), and the vertical coordinate is the true positive rate (TPR). Survminer (version 0.4.9) and the survival package (version 3.2-10) were used to draw survival line charts. The Cox proportional hazard models for 1-, 5-, and 10-year biochemical recurrence (BCR)-free probability were obtained using the "rms" R package, and the variables included age, gender, race, and CRGs. Based on the data reported by Ceccarelli et al. (22), we analyzed the co-expression of chromosome arms 1p and 19q (1p/19q codeletion) and isocitrate dehydrogenase (IDH) status and their relation to CRGs in GBM.

Statistical analysis
R (version 3.6.3) and related packages were used for all statistical analyses. The Wilcoxon test was used to compare two independent non-parametric samples. Statistical significance in qPCR and copper assay experiments were calculated using Student's t-test. We used mean ± standard deviation to describe continuous variables with normal distribution, and correlated variables without normal distributions were examined via Spearman's correlation analysis. Statistical significance was set at p < 0.05.
To validate the associations between differential expression levels of CRGs and GBM, we collected two independent validation GEO datasets (GSE4290 and GSE15824) and performed a meta-analysis to evaluate the overall effect. We found differences in the expression of LIPT1 and CDKN2A in different datasets. Therefore, we adopted the random-effects model for the meta-analysis and found that there was a significant difference in the expression level of LIPT1 (95% CI: 0.34-1.66, p = 0.003) in GBM tissues. Although the expression of CDKN2A in tumor tissues of the GSE4290 dataset and TCGA-GTEX dataset was significantly lower than that of normal tissues, after adding the data of the GSE15824 dataset, meta-analysis showed that there was no significant difference in the expression of CDKN2A between tumor tissues and normal tissues (Supplementary Figure S2). In addition, we investigated the correlation between different gene expression levels and found a strong correlation ( Figure 1B). For example, we found that DLAT expression was highly positively correlated with that of PDHA1 (r = 0.755, p < 0.001) ( Figure 1C).
We used OncoPrint from a query for alterations in the CRGs in GBM; 15.1% of the CRGs in the samples per patient were altered, and CDKN2A had the highest level of alteration (13%) ( Figure 1D). FDX1 and DLAT showed gene amplification (red bars), GLS showed homozygous deletions (green bars), and CDKN2A showed nonsynonymous mutations (blue bars).

Protein-protein interaction analysis and functional enrichment of CRGs
After setting up medium confidence (0.400) and no more than 50 interactors (1st shell), protein-protein interaction network analysis was performed to explore the interactions of CRGs, showing that DLD, PDHB, DLAT, and PDHA1 were hub genes (Figure 2A).
To further clarify the biological function of CRGs, related approaches were analyzed using the GO and KEGG databases. The main biological processes of the three CRGs in the GO analysis  were iron-sulfur cluster binding; oxidoreductase activity, acting on the aldehyde or oxo group of donors; and oxidoreductase activity, acting on the aldehyde or oxo group of donors with NAD or NADP as an acceptor. The cellular components were dihydrolipoyl dehydrogenase complex, oxidoreductase complex, and mitochondrial matrix. Molecular functions included the TCA cycle, acetyl-CoA biosynthetic process from pyruvate, and acetyl-CoA biosynthetic process. In addition, the TCA cycle, pyruvate metabolism, and glycolysis/gluconeogenesis were closely related to the CRGs in the KEGG analysis ( Figure 2B).

Correlation between CRG expression and immune infiltration in GBM
Tumor-infiltrating immune cells (TIICs) are an indication of the host immune reaction to tumor antigens (25). TIIC, tumor, and stromal form an ecosystem in the tumor microenvironment and have shown potential prognostic value (26,27). We used the  Figure S3E). CD47 (28) and CD24 (29) are dominant innate immune checkpoints and part of the novel "don't eat me" signal that promotes cancer immune escape. Our results also showed that CD47 and CD24 were highly expressed in GBM tissues (Supplementary Figure S4A To validate the results obtained in the bioinformatic target prediction analysis, we used two independent shRNAs to silence   (Table 1) in 172 cell lines. Subsequently, we observed the gene expression levels of CRGs and found a positive correlation between CRGs and CD24 (p < 0.01) (Supplementary Figure S4B) and a negative correlation between CRGs and CD47 (p < 0.01) (Supplementary Figure S4C). As determined by the Kaplan-Meier curves, patients with high expression of CDKN2A (log-rank test, p < 0.001), PDHA1 (log-rank test, p = 0.001), or LIAS (log-rank test, p = 0.009) had a longer overall survival than patients with low expression, while patients with low expression of FDX1 (log-rank test, p < 0.001), DLD (logrank test, p < 0.001), DLAT (log-rank test, p < 0.001), or LIPT1 (logrank test, p < 0.001) had a longer overall survival than patients with high expression ( Figure 5B) (Supplementary Figure S5).

Nomogram development and validation for GBM
To facilitate the clinical application of the prediction model, we created a nomogram to predict the 1-, 5-, and 10-year survival  probabilities based on the patient characteristics and CRGs. The top row of the nomogram is marked with "points" to calculate the points associated with 13 variables. Subsequent lines ("age" to "LIPT1") are the risk factors (variables) used in the model, reflecting a relatively excellent predictive performance of the nomogram. The data indicate that the low CDKN2A expression and high FDX1 expression are proportional to the survival probability of GBM patients ( Figure 6A). Model calibration was assessed using calibration curves that measured the relationship between the results predicted by the model and those observed in the cohort. The model made predictions that were close to the actual results ( Figure 6B). Finally, we evaluated the differences in CRG expression between other pathological conditions of glioma and GBM. The results showed that the expression of FDX1, LIPT1, and DLAT was highest in GBM (Supplementary Figure S6). IDH status was the molecular marker in the 2016 WHO classification of GBM (30); here, we analyzed the IDH mutation status of CRGs in patients with GBM and found that the status of mutant FDX1 (p < 0.0001), LIAS (p < 0.0001), PDHB (p < 0.0001), DLD (p < 0.0001), PDHA1 (p < 0.0001), CDKN2A (p < 0.0001), and LIPT1 (p < 0.0001) was significantly different from that of the wildtype genes (Supplementary Table 2). The 1p/19q codeletion highly benefits diagnosis and prognosis in gliomas (31); our data showed that the FDX1 (p < 0.0001), LIAS (p < 0.0001), DLAT (p < 0.0001), DLD (p = 0.032), PDHA1 (p < 0.0001), MTF1 (p < 0.0001), and GLS (p < 0.0001) of 1p/19q non-codeletion (non-codel) GBM were significantly different from those of 1p/19q codeletion (codel) GMB (Supplementary Table 3).

A CRG-drug/molecule-pathway network reveals novel treatment strategies
To predict potential drugs for the treatment of GBM, we systematically evaluated the relationship between CRGs and drug/ molecule responses of related pathways and constructed a multigroup integrated interaction network. As shown in Figure 7, in the CRG-drug/molecule-pathway interaction network, the most important gene for interaction was CDKN2A, and the most important pathway was the drugs/molecules of pyruvate To predict the drug sensitivity of mutated or wild-type CDKN2A, we systematically evaluated the relationship between CDKN2A in GBM and drug sensitivity and resistance in the GDSC database. Using the IC 50 values, we calculated the correlation between the IC 50 and CDKN2A scores of the drugs/ molecules. Combining the correlation results and drug treatment information, we obtained five drug sensitivity candidate groups and one drug resistance candidate. GBM with CDKN2A mutations was significantly sensitive to BDOCA000347e (p = 0.00745), dihydrorotenone (p = 0.00277), remodelin (p = 0.0177), OF-1 (p = 0.00145), and TAF1_5496 (p = 0.0151), and was significantly resistant to VX-11e (p = 0.000307) (Supplementary Figure S7).

Discussion
GBM is the most common primary malignant tumor of the central nervous system (32). Owing to the high proliferative and invasive ability of GBM and its insensitivity to radiotherapy and chemotherapy, the prognosis of patients with GBM is poor (33), with 5-year survival rates of less than 5% (34-36). Understanding the molecular basis of CRGs is a critical step in the continued evolution of precision medicine and cancer therapy. In this study, we examined the expression characteristics of CRGs in GBM tissues and performed functional enrichment and linear regression analyses of immune checkpoints. To the best of our knowledge, this is the first study to construct a prognosis score and target drug prediction for CRGs and GBM to provide an accurate treatment plan for GBM. Heavy metal ions are essential micronutrients, but insufficient or excessive metal content causes cell death (37). For example, ironinduced death has been defined as an iron-dependent form of nonapoptotic cell death (38). Studies have shown that intracellular Cu can induce a new form of regulatory cell death that is different from oxidative stress, which is known as cuproptosis (6). Mitochondria regulate cell death induced by copper ionophores through the TCA cycle. Our results showed that the primary pathway for CRG enrichment is the TCA cycle ( Figure 2B). Growing cells under hypoxic conditions (1% O 2 ) attenuates copper ionophore-induced cell death, whereas forced stabilization of the HIF pathway with the HIF prolyl hydroxylase inhibitor FG-4592 under normoxic conditions (21% O 2 ) does not (6). Consistent with this finding, our research data showed that oxygen ranks first among the related drugs/molecules predicted by CRGs; therefore, we believe that appropriate hyperbaric oxygen therapy can increase cuproptosis in tumor tissue.
GBM immunotherapy has recently attracted considerable attention (39), and understanding the mechanisms of GBM immunosuppression will help us to develop immunotherapy strategies (40). Today, immunotherapy options for GBM remain limited (41), and there is an urgent need for new and effective targets. The present study analyzed immune cells and checkpoints in CRGs of GBM to establish a potential strategy for GBM immunotherapy. Our research data showed that a variety of CRG molecules are highly correlated with macrophage infiltration and significantly correlated with the expression of t CD47 and CD24, which are immune checkpoints of the "don't eat me" signal. Cancer cell death is also usually accompanied by downregulation of "don't eat me" signals (42).
In the 1990s, the first phagocytosis checkpoint axis, CD47signal regulatory protein alpha (SIRPa), was identified (43). CD47 is the ligand of SIRPa, and is involved in inflammatory response and is recognized as an immune checkpoint for tumor evasion (44), which operates as a "don't eat me" signal (45). CD24 also operates as a "don't eat me" signal that promotes cancer immune escape (29). Previous studies have shown that CD24 and CD47 expression are upregulated in GBM cells (46)(47)(48). In this study, GBM tissue had significantly increased expression of CD24 and CD47 (Supplementary Figure 4A). CD47 expression levels were positively correlated with those of FDX1, DLAT, DLD, GLS, LIPT1, PDHA1, PDHB, MTF1, and LIAS and negatively correlated with those of CDKN2A in GBM ( Figure 4A). CD24 expression levels were positively correlated with those of LIPT1, A B FIGURE 6 Nomogram development and validation. (A) Univariate Cox regression of CRGs, age, gender, and race to predict the 1-, 3-, and 5-year survival probabilities. (B) Calibration curve for the CRG nomogram model in GBM. The abscissa is the survival probability predicted by the model, the ordinate is the actual observed survival probability, and the gray diagonal is the ideal line.
CDKN2A, PDHA1, MTF1, and LIAS and negatively correlated with those of GLS in GBM ( Figure 4B). These findings provide new insights for developing immunotherapies against GBM.
Primary GBM comprises 90% of the cases and is IDH wild-type, while secondary GBM develops from lower-grade glioma and carries mutations in IDH (49). In this study, the IDH status of GBM with mutant types of FDX1, LIAS, PDHB, DLD, PDHA1, CDKN2A, and LIPT1 was significantly different from that of GBM with the wild-type genes (Supplementary Table 2). Furthermore, in the CRG-drug/molecule-pathway interaction network, the most important gene for interaction was CDKN2A (Figure 7). The GDSC database is a resource for biomarker discovery for the development of therapeutics for cancer cells and contains information from 138 anticancer drugs across 696 cell lines (23, 50). We observed that the CDKN2A gene predicted sensitivity and resistance to drugs in most GBM cell lines in the GDSC database. For example, GBM with mutated CDKN2A was significantly sensitive to BDOCA000347e, dihydrorotenone, and remodelin, but significantly resistant to VX-11e (Supplementary Figure 7). These findings provide a basis for the selection of therapeutic drugs.
Our study has several limitations. First, although immunoassay and clinical prognostic scores focusing on the expression of CRGs in GBM performed well, other significant genes with potential predictive values were not included in this study. Second, given that the prognostic signature and target drug prediction were built and validated by exploiting data from public databases, although we conducted qPCR experiments in tissue and A172 cell line, further experimental validation is required.

Conclusion
In summary, this study systematically analyzed the landscape of molecular alterations and interactive genes involved in cuproptosis in GBM. Our study suggests that these CRGs may play a crucial role in GBM outcome. The risk score based on CRG expression characteristics can predict the survival and prognosis of GBM patients, and is significantly associated with immune infiltration levels and the expression of CD47 and CD24, the immune checkpoints of the "don't eat me signal." Our results provide new insights into CRG-related target drugs/molecules for cancer prevention and treatment. Future biological studies are required to confirm our findings.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors. The CRG-drug/molecule-pathway network reveals novel treatment strategies. The network, based on Genecard, Drugbank, and GDSC databases, shows the interaction network of CRGs (pink circle), pathways (yellow hexagons), and drugs/molecules. The size of the shape increases with the increase in the degree of interaction, which reflects the relationship with the CRG-drug/molecule-pathway.

Ethics statement
The study was approved by Xijing Hospital of Fourth Military Medical University (KY20193098). The patients/participants provided their written informed consent to participate in this study.

Author contributions
EL, XZ, XM, and BL were responsible for the main experimental concept and design. XM collected samples and provided ethical proof. qPCR was performed by JS and LL. Analyses were performed by HQ, QM, YH, and SL, and the manuscript was written by EL. All of the authors approved the final version. All authors have read and agreed to the published version of the manuscript.

Funding
This work was financially supported by the National Natural Science Foundation of China (No. 82174166) and the Innovation Capability Support Program of Shaanxi (No. 2021TD-45).