Cuproptosis characteristic gene consistent clustering to identify sample subgroups
The RNA-seq data of 169 TCGA-GBM samples were obtained, and the tumor samples were clustered into two groups based on cuproptosis genes (FDX1, LIAS, LIPT1, DLD, DLAT, PDHA1, PDHB, MTF1, GLS, and CDKN2A) through consistent cluster analysis (Fig. 1A, B, C). Significant differences were not observed in the survival between the two subgroups of samples, suggesting that these 10 genes alone may not be able to characterize the effect of cuproptosis mechanism on patient survival benefit (Fig. 1D). We also observed the expression patterns of these 10 characteristic genes in two subgroups of samples (Fig. 1E). However, significant differences were observed in the landscape of mutation, immune checkpoint expression level, and cancer hallmarks between the two subgroups. First, the SNV mutation frequencies of TP53 and other genes showed significant differences between the two types of samples (Fig. 2A). The map of CNA showed that both types of patients had significant amplification on chromosome 7 (Fig. 2B), and significant difference was not observed in total frequency of CNA (Fig. 2C). In addition, significant differences were observed in intratumoral heterogeneity between the two subgroups (Fig. 2D). Moreover, we observed a significantly different expression level of immune checkpoint genes PD-1, IDO1, and LAG3 in the two subgroups (Fig. 2E). Significant differences were observed in the cancer hallmarks of fatty acid metabolism, KRAS, P53, NOTCH, and PI3K/AKT/MTOR signaling pathway between the two subgroups (Fig. 2F).
Construction of cuproptosis activation scoring model based on differentially expressed genes between the sample subgroups
First, the DEGs between the two groups of samples were identified based on DESeq2 (Fig. 3A). The functions of these DEGs were enriched in the cell cycle related processes, protein kinase activity, P53 signaling pathway, and TGF-βsignaling pathway (Fig. 3B, C). In the TCGA-GBM sample set, we identified 14 candidate prognostic marker genes that were significantly associated with the OS of patients based on univariate Cox regression (Fig. 3E) and further filtered the redundant factors using LASSO Cox to obtain 11 prognostic marker genes (Fig. 3D). Based on the PCA of these 11 genes, their contribution to principal components 1 and 2 were used as coefficients (Fig. 3F) to construct CuAS.
The prognostic efficacy of CuAS
In the training set, we scored the samples based on these 11 genes and the results of the PCA (Fig. 4A). We found that patients with higher CuAS scores had significantly worse OS (Log-rank P < 0.0001, Fig. 4B) and 6 month AUC (95%CI) = 0.625 (0.608–0.643), 1 year AUC (95%CI) = 0.69 (0.646–0.732), 2 year AUC (95%CI) = 0.797 (0.735–0.852), 3 year AUC (95%CI) = 0.825 (0.765–0.885) (Fig. 4E). In addition, two validation sets from the CGGA database showed that CuAS had stable prognostic efficacy (Log-rank P = 0.0075, Fig. 4C, Log-rank P = 0.0043, Fig. 4D). Univariate Cox regression analysis was performed on CuAS and multiple clinical features to evaluate the independence of the prognostic efficacy of CuAS. The results showed that CuAS was significantly associated with patient prognosis (HR (95% CI) = 7.51 (3.75, 15.05)) (Fig. 4F). Consistent results were also observed in the independent validation sets (HR (95% CI) = 1.74 (1.133, 2.593)) (Fig. 4G). Furthermore, we constructed a multivariate Cox regression model for CuAS and multiple clinical features and found that CuAS could still serve as an independent prognostic factor (HR (95% CI) = 7.35 (3.23, 16.7)) (Fig. 4H), HR (95% CI) = 1.90 (1.12, 3.2)) (Fig. 4F)).
Epiregulin (EREG) was an oncogenic gene that can influence immunity and cuproptosis
The EREG mRNA expression levels were high in tissues and multiple glioma cell lines (Fig. 5 A, B). We used Western blot (Fig. 5C) and IHC staining (Fig. 5E) to detect the EREG protein expression levels in tumors and normal tissues, and we found that the protein expression levels of EREG in tumors were higher than that in normal tissues. Additionally, we found that the protein expression levels of EREG in glioma cell lines were higher than that in normal astrocytic cell lines (Fig. 5G). Subsequently, we constructed knockdown stable cell lines of EREG and verified the knockdown effects on mRNA (Fig. 5F) and protein levels (Fig. 5D). Functional experiments showed that EREG knockdown (KD) can significantly inhibit the proliferation detected by Edu exepriments (Fig. 7E), invasion (Fig. 6D), migration (Fig. 6E), and colony forming ability (Fig. 6C) of tumor cells. Additionally, flow cytometry cell cycle assays suggested that EREG KD significantly inhibited cell cycle progression from the G0/G1 phase to the S phase (Fig. 7F). To explore the relationship between EREG and immune infiltration, we detected the expression level of PDL1 in EREG-KD group and found that PDL1 also decreased (Fig. 6A). To explore the relationship between EREG and cuproptosis, we performed the different treatment gradients of Cu-Elesclomol(ES) (1:1) on U251 cell lines and found that cell viability decreased with increasing time, and the effect of ES-Cu required a specific concentration range (5–50 nM) (Fig. 7A, B). Subsequently, we treated tumor cells with the same concentration (30 nM) with ES-Cu, and observed cell viability at 0, 12, 24, 36, 48, 60, 72, 84, and 96 h (Fig. 7D). It was found that the proliferation rate of the treated cells decreased significantly compared with the cells that were not treated with ES-Cu. Then, the same treatment was performed on the shEREG and shNC groups and found that the proliferation rate of the two treated groups decreased significantly when compared with the shEREG and shNC groups that were not treated with ES-Cu; however, the reduction rate of the proliferation of the shEREG group was higher than that of the shEREG group with ES-Cu treatment, indicating that EREG can influence cell proliferation by affecting the process of cuproptosis (Fig. 7C, D). Therefore, we detected the protein expression level of FDX1 in the shEREG and shNC groups. The results showed that FDX1, the core regulatory protein in cuproptosis, was down-regulated in the shEREG group (Fig. 6B). Based on the above results, we believe that EREG is an oncogenic gene that can affect immunity by influencing the expression level of PDL1 and is closely related to the process of cuproptosis.
Single cell transcriptome analysis of CuAS patterns
Based on the downloaded single cell data (GSE173278, 29339 cells,10X Genomics platform), R Package Seurat was used to process the data. The expression profile was transformed by Log10, and 2,000 highly mutated genes were identified based on the VST method. Subsequently, principal component analysis and dimensionality reduction visualization were performed using UMAP. As a strong batch effect was observed, Harmony was used for batch correction (see Additional file 2: Fig. S1). Follow-up analysis was conducted based on the corrected data. The default parameters were used for clustering, and the meaning reference of analysis results based on known markers (SingleR was BP and HPCA) for cell type annotation (Fig. 8A). The cells in the GBM samples were divided into 7 categories, three malignant cell (OLIG1 + malignant, n = 11637; VEGFA + malignant, n = 6446; CENPF + malignant, n = 5363), microglia (n = 3219), fibroblasts (n = 1020), endothelial cells (n = 919), and oligodendrocytes (n = 735) (Fig. 8C). CuAS was calculated by using the previous model coefficients, but several cuproptosis characteristic genes were not detected in single cell data, and the expression level of many characteristic genes was undetectable (see Fig. 8E and Additional file 2: Fig. S2). A small number of cells (approximately 2800 cells) with a high CuAS score accounted for less than 10% of the whole cells and were distributed in multiple cell subpopulations. Most of the subpopulations contained less than 5% of cells with high CuAS, and 1132 cells with high CuAS were present in VEGFA + malignant cells (hypergeometric test, p value < 0.05). Trajectory inference of tumor cells are depicted (see S Additional file 2: Fig. S3). Moreover, the ancestor clone was determined based on CNV in combination with the idea of clone evolution, so as to determine the evolutionary relationship between cells more accurately. Based on CNV, we found that the OLIG1 + malignant cell may be the ancestor clone (Fig. 8D). Furthermore, in order to better explain the functional role of CuAS at the single cell level, the functions of VEGFA + malignant cells were observed and functional enrichment analysis was conducted based on specific up-regulated genes obtained by differential expression analysis. Pathways are mainly enriched in pathways related to hypoxia and oxidative stress (Fig. 8B).
Differential activation of transcription factors between high and low CuAS
Annotated files of human transcription factors were obtained from the RcisTarget database and the list of human transcription factors were downloaded. Transcription factor regulatory network was constructed using pySCENIC. Subsequently, the AUCell algorithm was used to calculate the activity of each transcription factor, and according to the CSI between the different transcription factors, four regulatory modules were identified (Fig. 9A). Module score was performed for each cell sample. We explored the association between cell type and module score, which revealed that the score of Module1 was significantly higher in VEGFA + malignant cells, while the score of Module2 was significantly higher in the endothelial cell subset. The score of Module3 was significantly higher in the microglia cell subset, while the score of Module4 was higher in the CENPF + malignant cells and partial OLIG1 + malignant cells, which reflected the differences of TF activated by different malignant cell subsets (Fig. 9B, C). Similarly, the results of the hTFtarget database showed that transcription factors such as FOSL2, JUND, NFIC and PBX3 were highly active in the VEGFA + malignant subgroup (see Additional file 2: Fig. S4). In particular, we previously observed that cells with high CuAS scores were concentrated in VEGFA + malignant subsets, showing the potential association between CuAS scores and Module1.
Correlation between CuAS and immune microenvironment
First, GO/KEGG enrichment analysis based on GSVA algorithm was performed on TCGA-GBM tumor samples, and immune-related pathways were differently enriched between the high and low CuAS groups, including T cells, NK Cell, B cell signal, chemokine signal, cytokine interaction, and other pathways (Fig. 10A). GO enrichment also showed that the activation differentiation and proliferation of T cells, NK cell proliferation, cytotoxic reaction, and other characteristics were highly enriched when CuAS scores were high (Fig. 10B). In addition, we identified differentially enriched signatures between the high and low CuAS groups based on GSEA enrichment analysis and also captured immune reaction processes such as leukocyte adhesion migration and T cell activation (Fig. 10C), indicating that the mechanism of cuproptosis is closely related to immune reaction process. Further, we calculated the immune and stromal components by using ESTIMATE, and it was observed that CuAS was significantly positively correlated with the stromal, immune, and ESTIMATE scores, while it was significantly negatively correlated with tumor purity (Fig. 10D–G), which suggested the association between cuproptosis and immunity. By using GSVA to calculate immune cell infiltration, we found that cuproptosis was significantly correlated with various types of immune cell infiltration, including activated DC and NK cells. (Fig. 12H). In addition, CIBERSORT and xCell methods were used to calculate various immune cell infiltrates, which revealed similar results (see Additional file 2: Fig. S5).
Specific cell communication was different between high and low CuAS groups
As the high CuAS cells were significantly enriched in VEGFA + malignant cells, we mainly analyzed the difference in communication and function between the VEGFA + malignant and other cells. Extensive cell communication was observed in each cell subpopulation (Fig. 11A). Furthermore, by distinguishing between the incoming and outgoing signals, we found that fibroblasts are the dominant signaler of outgoing signaling, and VEGFA + malignant cells are the signal receivers (Fig. 11B). Furthermore, we identified two patterns of cell subpopulations in outgoing signaling, in which VEGFA + malignant cells belonged to Pattern 2 and corresponding pathways included VEGF, FGF, CDAM, CD22, ADGRE5, and other malignant progression related pathways (Fig. 11C). Meanwhile, we analyzed the dominant signaling pathway of each cell and found that VEGFA + malignant cells are not only involved in the signaling pathway of VEGF, but also in the CD99 signaling pathway, which was not proposed in the non-negative Matrix Factorization (NMF) analysis (Fig. 11D). The CD99 signaling pathway plays an important role in tumor progression and transendothelial migration of cancer cells. VEGF and CD99 signaling pathways were further analyzed, and it was found that VEGFA + malignant cells are the dominant signalers of VEGF signals, and the cell subsets that were affected are mainly the endothelial cells and fibroblasts, both of which are important components of angiogenesis (Fig. 11 E, G). In addition, VEGFA + malignant cells are the dominant signaler, receiver, and influencer of CD99 signaling pathways, indicating that CD99 signaling pathways can occur as feedback loops (Fig. 11 F, H). Meanwhile, endothelial and fibroblast cells are also affected by CD99 signaling pathways, suggesting that VEGEA + malignant cells can influence transvascular endothelial migration (Fig. 12). The immunofluorescence detection of tissue samples with high and low CuAS showed that VEGFA and CD99 were also highly expressed in tissues with high CuAS. The results were opposite in tissues with low CuAS (Fig. 13E, F, G), which provided a new idea for the intervention of cuproptosis-related tumor cells.
CuAS is associated with prognosis of immunotherapy
Based on all the immunotherapy data searched, we observed the ability of the CuAS score in predicting the prognosis and efficacy in the immunotherapy cohort. Phs001493 (Renal cell carcinoma, Anti-PD1 therapy) and PRJEB23709_ipiPD1 (Melanoma,anti-CTLA4 & AMP; Anti-pd1 dual antibody therapy) were significantly associated with worse prognosis (see Additional file 2: Fig. S6B and D). For a patient’s Progression Free Survival (PFS), we found that NCT02684006 (kidney cancer, anti-PDL1 treatment) was significantly associated with worse prognosis (see Additional file 2: Fig. S6E). Therefore, high CuAS patients may benefit from immunotherapy.
Potential targeted drugs for high CuAS glioblastoma cells
The expression data of cell lines were extracted from three databases: GDSC, CCLE, and CTRP. A lower AUC value represents a higher sensitivity to drugs. Using the AUC data provided by these databases, multiple drugs with a negative correlation between the AUC and signature were found in GDSC, such as methotrexate, BMS -708163, YM201636, FR -180204, and NU − 7441 (Fig. 12A). A variety of drugs with positive correlations were also found, such as cyclopamine (Fig. 12B). These drugs showed significant differences in the AUC between the groups of high and low signature (Fig. 12 C, D). No drugs with an IC50 significantly correlated with signature were found in CCLE, while drugs with a significant AUC (KU -0063794, cytochalasin B, GDC − 0941, cabozantinib, CI − 976, SJ -172550, SGX − 523, BRD − K71935468, temozolomide, AT7867, BRD-K66532283, palmostatin B, GDC-0879, ETP-46464, and NVP-BEZ235) negatively correlated with signature were found in CTRP (Fig. 12E, F). These drugs were found in the AUC values of the high and low signature groups were significantly different in CTRP (Fig. 12 J, H). Therefore, high CuAS samples are likely to be sensitive to these compounds, and these compounds may be novel treatment options for GBM.
Experimental validation of model genes
The genes expression levels in the model were detected by qRT-PCR, and the results showed that they were highly expressed in 20 pairs of tumor and normal tissues (UNCX, SLC6A3, AGAP2-AS1, LINC00968, PTX3 and SBSPON), while ITPRID1, DCST2, ETV3L, and ENSG00000261327 were down-regulated (Fig. 13A). According to the corresponding PCR results, we divided the tissue samples into high and low CuAS groups. IHC staining was performed on UNCX, SLC6A3, and PTX3, and it was found that the protein expression of the high CuAS group was higher than that of the low CuAS group (Fig. 13B, C, D).