Introduction

Gliomas are the most common and highly aggressive primary cancers of the central nervous system (CNS) cancer.¹ They primarily consist of adult astrocytomas, oligodendrogliomas, and glioblastomas (GBM). Astrocytomas and oligodendrogliomas are characterized by high heterogeneity, relatively slow growth, a propensity for recurrence, and the potential to progress into GBM.^2,3 Different types of gliomas, such as astrocytomas, originate from astrocytes, making them fundamental and unique in tumor formation.⁴ GBM, the most malignant form of astrocytoma, represents the most aggressive primary glioma of astrocytic origin, with a dismal median survival of less than 16 months.⁵ Despite the current standard treatment, which includes surgery, chemotherapy, and radiation, tumor recurrence is ultimately inevitable.³ The glioma microenvironment is highly heterogeneous, influencing tumorigenesis and structural evolution, and exhibits both inter-tumoral and intra-tumoral heterogeneity due to distinct phenotypic variations, tumor cell types, and cellular states.^6,7 Therefore, a comprehensive exploration of the heterogeneity of cellular subtypes within the glioma microenvironment is essential for understanding the pathogenesis of gliomas and for the development of new therapeutic strategies.

The tumor microenvironment (TME) of gliomas is composed of heterogeneous subsets of tumor cells, astrocytes, tumor-associated macrophages, and various immune cells that play distinct roles in tumorigenesis, progression, and response to clinical therapies.^8,9 Among these, astrocytes, which are the most abundant glial cells, interact with T cells, monocytes, and microglia, and are involved in the regulation of the TME, thereby influencing therapies that modulate T cell responses in gliomas.^10–12 These findings suggest that the study of the TME in gliomas should be conducted at the single-cell level. Single-cell RNA sequencing (scRNA-seq) offers significant advantages in elucidating the heterogeneity within gliomas by analyzing the expression data of individual cells.¹³ Furthermore, while astrocytes have been shown to either inhibit or promote the growth of brain tumors in mice, little is known about the potential roles of astrocyte subpopulations within the TME of gliomas.^14,15

Mitochondria are crucial for energy generation and metabolism, and the survival of brain cells, particularly astrocytes. Astrocytes play a key role in maintaining brain redox balance, bioenergetics, and neurotransmission.^16,17 Oxidative stress-induced lipid peroxidation, resulting from impaired mitochondrial metabolism, can lead to ferroptosis in astrocytes.¹⁸ Cuproptosis, a newly identified form of copper-induced cell death, facilitates the aggregation of mitochondrial proteins and disrupts the function of mitochondrial metabolic enzymes. This process is characterized by copper-dependent intracellular reactive oxygen species (ROS) production and excessive accumulation of lipid peroxides.¹⁹ Copper is an essential micronutrient that maintains normal cellular function and contributes to oncogenesis due to its critical biological roles in cell signaling, oxygen transport, and enzyme activity, particularly in promoting ROS production.^20,21 Copper chelators can induce immune checkpoint inhibition and reduce the expression of PD-L1,²² thereby enhancing anti-tumor immune responses and improving treatment efficacy.^23,24 Copper toxicity may be mediated by the synthesis of lipoylated enzyme genes, including FDX1, CDKN2A, and MPC1, which, in conjunction with lipoylated proteins, play vital roles in cancer progression.^19,25 Therefore, the cuproptosis of astrocytes, driven by copper toxicity and mitochondrial metabolic damage, may represent a novel therapeutic approach for gliomas. Additionally, the composition and cellular interactions within the TME are closely related to the prognosis of gliomas.²⁶ However, the potential role of cuproptosis-related genes (CRGs) in the TME and their implications for prognosis prediction in gliomas remain unclear.

In this study, we aimed to comprehensively identify the tumor microenvironment in gliomas and characterize the CRGs within astrocyte subpopulations. We conducted an analysis of single-cell RNA sequencing using samples from three in-house glioma patients, along with data from 14 previously published patients from the Chinese Glioma Genome Atlas (CGGA). Additionally, we validated the characteristics of CRGs using bulk RNA sequencing data from both the CGGA and The Cancer Genome Atlas (TCGA) datasets. Furthermore, we established and validated a CRG-score system in different independent cohorts to quantify the cuproptosis characteristics of individual tumors and predict the prognosis of gliomas.

Materials and Methods

Study Design and Data Acquisition

Published single-cell RNA sequencing unique molecular identifier (UMI) counts from 14 glioma patients, encompassing 6,148 cells, were obtained from the CGGA dataset.²⁷ The scRNA-seq was conducted on an Illumina 4000 platform. Bulk expression data and clinical information were sourced from 656 glioma patients in the CGGA dataset (http://www.cgga.org.cn/) and 628 gliomas patients from TCGA-GBM and TCGA-LGG (low-grade glioma) (https://portal.gdc.cancer.gov/repository). We applied batch correction to the mRNA sequencing data using the ComBat algorithm for both the CGGA and TCGA datasets. The clinical features for the CGGA and TCGA datasets included age, gender, tumor grade, 1p19q codeletion, IDH mutation, MGMT methylation status, follow-up time, and survival status.

Sample Collection and Preparation

Three samples (two GBM and one LGG with tumor-only samples) were collected for single-cell RNA sequencing, along with ten samples containing both tumor and adjacent normal tissues for bulk RNA sequencing. These samples were obtained from glioma patients who had undergone surgical resection at the First Affiliated Hospital of Harbin Medical University. Pathological diagnoses were conducted independently by at least two qualified pathologists. Fresh tissues were immediately dissected into fractions and enzymatically digested to obtain single cells. The study received ethical approval from the Institutional Ethics Committee of Harbin Medical University and was conducted in accordance with the Declaration of Helsinki. Written informed consent was obtained from each patient.

The glioma tissue was dissociated using a tumor dissociation kit (Miltenyi Biotech, Germany) following the manufacturer’s instructions. The dissociated tissue was then stored in Wash Cell MACS Tissue Storage Solution supplemented with RPMI 1640 and 10% fetal bovine serum (FBS).

Subsequently, the samples were sieved through a 40 µm cell strainer and centrifuged at 300 x g for 7 min. The supernatant was completely removed, and the pelleted cells were suspended in red blood cell lysis buffer (Solarbio) to lyse the red blood cells. The cells were resuspended in RPMI1640 + 10% FBS and re-filtered through a 35-mm cell strainer. Finally, resuspend cells as required were used for scRNA-sequencing.

ScRNA-Sequencing Library Construction and Processing

Each cell suspension was generated to 3′ single-cell RNA sequencing libraries using the Chromium Next GEM Chip G Single Cell Kit, Single Cell 3′ Library, and Gel Bead Kit V3.1 (10x Genomics). The 10x Genomics Chromium platform was utilized to construct a 10x barcoded cDNA library. The 10x libraries were sequenced on an Illumina Novaseq 6000 sequencing platform (Illumina, San Diego, CA) and subsequently mapped to the GRCh38 human reference genome using the Cell Ranger toolkit (version 4.0.0). The gene-barcode matrix was then generated following the Seurat (version 4.1.3) pipeline in R software (version 4.1.1). Low-quality cells, defined as those with a detected gene count of fewer than 200 or more than 6,000, or those exhibiting a high mitochondrial transcript ratio (>20%), were excluded from the analysis. The total number of transcripts in each cell was normalized to 10,000 using log transformation.

Identification of Glioma Cells and Subtypes

After normalization and scaling, we utilized Seurat to eliminate batch effects among patients using the mutual nearest neighbor algorithm. Additionally, we regressed out mitochondrial contamination and cell cycle effects. Principal component analysis (PCA) was performed to identify the top 2,000 highly variable genes. A subset of significant principal components (PCs) was selected for cluster analysis using the PC ElbowPlot function in Seurat. Cell clustering, sub-clustering, and visualization via t-distributed stochastic neighbor embedding (tSNE) were conducted separately using the FindClusters and RunTSNE functions in the Seurat package, with appropriate resolutions. Cell identity annotations for each cluster were defined using the SingleR package and subsequently verified manually. The known marker genes for cell types included epithelial cells, stromal cells, myeloid cells, T cells, B cells, plasma cells, and undefined astrocytes.^28–30

Cell-Cell Communication Analysis

The key signals and latent intercellular communication networks from the annotation of cell identify on different cell clusters were identified by the R package “CellChat” (version 1.0.0), which contained ligand-receptor interaction databases using a non-negative matrix factorization (NMF) method that has been successfully applied in pattern recognition.³¹ The number of ligand-receptor and strength of cell-cell communication were estimated from the target cell cluster to other cell clusters using the function of aggregateNet in CellChat. Major signaling pathways for target cell groups and intercellular communication patterns were depicted. Additionally, we identified both outgoing and incoming signals of different cell subtypes across all NMF clusters. By utilizing the computeNetSimilarity function, we classified signaling groups based on the similarity of their functions and structures.

Developmental Pseudo-Time Analysis

Trajectory analysis of astrocyte types in glioma was conducted using the Monocle R package to investigate potential lineage differentiation among astrocyte subpopulations.³² First, significant marker genes in each cluster were identified using the “differentialGeneTest” function within the Monocle package. The cell differentiation trajectory was then analyzed with two-dimensional reduction parameters via the DDRTree method, and astrocytes were ordered using the “orderCells” function. Subsequently, pseudo-time-dependent genes were identified through the “differentialGeneTest” function. Finally, differentially expressed genes between branches were visualized to identify co-expressed gene modules using the BEAM function.

Construction and Validation of the Cuproptosis-Related Genes Model

According to a false discovery rate (FDR) threshold of <0.05, we identified 27 cuproptosis-related genes from the study by Peter Tsvetkov et al.¹⁹ The intersection of cuproptosis-related genes and pseudotime-dependent genes was determined using the Venn function in the R package. To assess the overall prognostic impact of these CRGs in glioma patients, we constructed prognostic risk models using the training dataset from the CGGA glioma expression data (N=656). In the CGGA training dataset, the model for cuproptosis-related genes was developed based on the mRNA expression levels of the intersecting genes and their corresponding coefficients, as determined by multivariate Cox regression analysis.

where Expri expressed these cuproptosis-related intersection genes, the multivariate model calculated the Coefi coefficient, and the CRG-score represents a multimode weighted sum of the prognostic value for each sample. Using the median as the cut-off value, CGGA glioma patients were categorized into low- and high-risk groups. The Kaplan-Meier survival curves were used to assess the prognostic performance of the model. Additionally, TCGA glioma dataset was utilized to validate the performance of the risk-score model as an independent testing set (N=628). The relative abundance of 22 immune-infiltrating cells between the low- and high-risk groups was quantified using the CIBERSORT algorithm. The Human Protein Atlas (https://www.proteinatlas.org/) database was used to verify the protein expression levels of the cuproptosis-related genes in glioma and normal tissues.

Targeted mRNA Sequencing and Data Processing

Total RNA from glioma tissues for mRNA sequencing was extracted using TRIzol (Invitrogen) reagent, followed by the removal of genomic DNA contamination using DNase. The mRNA was isolated for RNA sequencing library preparation for Illumina (HaploX Genomics Center, Jiangxi, China). RNA sequencing was performed on the Illumina NovaSeq 6000 instrument using the paired-end 2×150 sequencing mode. Raw reads were assessed for quality, and FastQC was utilized to ensure high read quality. The clean reads were aligned to the human reference genome (hg19) using HISAT2 software, and the read counts for each gene were calculated with featureCounts (v1.6.3).

Statistical Analysis

All statistical analyses of scRNA-seq data were conducted using R software (version 4.1.1). A P-value of less than 0.05 was considered statistically significant. The analyses utilized several packages, including Seurat, CellChat, SingleR, and Monocle. Differences between the two groups were assessed using the t-test, while correlations were determined using either Pearson’s or Spearman correlation tests.

Results

Landscape of the Glioma Ecosystem Revealed by scRNA-Seq

The study design used in this study is depicted in Figure S1. To systematically interrogate the intra-tumoral heterogeneity of the glioma ecosystem, we generated scRNA-seq profiles from three glioma patients. A total of 32,982 cells were acquired from these samples, and 32,293 cells were retained after quality control following standard procedures. Firstly, we used the harmony function to normalize the glioma samples and then reduced the dimensionality of the scRNA-seq data by PCA. Then, we observed 24 different cell clusters by tSNE approaches (Figure S2 and Table S1). According to the annotation of the uniquely established expressed marker gene and checked manually, we classified 24 cell clusters into five cell types, including those from myeloid cells, stromal cells, T cells, NK cells, and astrocytes (Figure S3A and Table S2). To further validate the heterogeneity of our scRNA-seq of the glioma ecosystem, we similarly analyzed the scRNA-seq of 14 glioma patients in the CGGA dataset based on tSNE, and found 17 clusters of cells annotated as myeloid cells, T cells, B cells, and astrocytes based on preferentially expressed marker genes (Figure S3B and Figure S4A and B). Astrocytes are considered as the most abundant glioma cell type in the central nervous system. In our in-house and CGGA datasets, we further divided into five subclusters and six subclusters of astrocyte, respectively (Figure S5). Finally, we identified nine similar cell types using tSNE in our in-house and CGGA samples (Figure 1A and B), including myeloid cells, T cells, and astrocytes.

Figure 1 Analysis of cuproptosis-related genes across various cell types in glioma. (A) tSNE visualization of our glioma single-cell transcriptional profiles. (B) tSNE visualization of CGGA glioma single-cell transcriptional profiles. Each dot represents an individual cell, with different colors indicating distinct cell types. (C) Whole-genome CRISPR-Cas9 positive selection screen utilizing two copper ionophores (Cu-DDC and elesclomol-copper) in cells. Overlapping hits with an FDR value of less than 0.05 were analyzed. (D) Bubble plot illustrating cuproptosis-related genes for each cell type in glioma. Abbreviations: tSNE, t-Distributed Stochastic Neighbor Embedding; FDR, false discovery rate.

Cuproptosis played an important role in tumor growth, angiogenesis, and immune response in gliomas.²⁰ A total of 27 cuproptosis genes were used in this study (Figure 1C). To further elucidate the role of these cuproptosis genes in various cell types within the tumor microenvironment, we calculated the proportion of different cell types for cuproptosis-related genes in each cell subset. The expression levels of cuproptosis genes on different cell types were depicted in Figure 1D, including CDKN2A, MPC1, RPL3, etc. Additionally, we evaluated the expression levels of cuproptosis-related genes in astrocyte subpopulations and found that seven genes of CDKN2A, REXO2, MPC1, RPL3, CAPRIN1, SOX2, and AHR have higher expression levels in the tSNE and violin plots (Figure S6 and S7).

Intercellular Communication Patterns Predict Potential Signaling Pathways

To further elucidate the plasticity of various cell subtypes within the tumor microenvironment, we used CellChat to quantitatively analyze intercellular communication networks among ligands, receptors, and their connections derived from our in-house glioma scRNA-seq data. In Figure 2A, we observed the aggregated cellular communication network, revealing significant changes in interaction strength between myeloid cells and astrocyte subpopulations. We further identified potential signaling pathways within these cell types and categorized these pathways based on their functional and structural similarities (Figure 2B). Four groups of pathways were determined based on functional similarity grouping. Group 1 consists of EGF, FGF, and PTN pathways, which represent promiscuous signaling. Group 2 includes IGF, CD80, and OSM, primarily representing signaling from astrocyte subpopulations. Group 3 includes PDGF, APP, and JAM, and largely represents signaling between stromal cells and T cells. Group 4 is dominated by inflammatory pathways (eg, IL, IFN, TGFβ, CCL). Meanwhile, structural similarity grouping also identified four groups of signaling pathways. CellChat analysis can simply predict putative function by grouping well-known signaling pathways. To further identify the detailed communication patterns for the key individuals in different cell types, we analyzed communication patterns that coordinate cell groups with signaling pathways, which determined five outgoing signaling patterns and three incoming signaling patterns (Figure 2C). The outgoing communication pattern indicated that all outgoing signaling from astrocytes is characterized by patterns 1 and 5, which represent pathways involving TGFβ, CXCL, IL, IFN, and HGF. Conversely, the incoming signaling to astrocytes is predominantly represented by pattern 1, which includes multiple pathways such as TGFβ, IFN, and CXCL.

Figure 2 Interactions in Cell-Cell Communication in Glioma. (A) A landscape of cell-cell network communications. The thickness of the edges represents the number (Left) or the weight (Right) of interactions among cell types. (B) Categorization of signaling pathways based on cell functionality (Left) and structure (Right). Each dot represents a distinct signaling pathway. (C) Outgoing (Left) signaling patterns of target cells and incoming (Right) signaling patterns of secreting cells.

To further elucidate the complex interactions within intracellular networks among astrocyte subpopulations and other immune cells, CellChat identified 95 signaling pathways across nine cell types, including TGFβ, BMP, TNF, CXCL, CCL, and MIF pathways, as demonstrated by heatmap analysis (Figure S8A). Our findings indicated that the TGFβ signaling pathway exhibited a higher relative strength in astrocyte subpopulations. The TGFβ signaling network revealed that myeloid cells were the predominant sources of TGFβ ligands, which operated in an autocrine way among myeloid cells and in a paracrine way from immune cells to diverse astrocyte subpopulations, particularly in astrocyte cluster 1 (Figure S8B). In addition, we observed that myeloid cells, as a primary source of ligands, functioned in both autocrine and paracrine manners from immune cells to astrocyte clusters 3 and 4 in the CGGA dataset (Figure S8C). Within the TGFβ signaling pathways, we systematically explored the predicted signals from astrocytes to immune cells and identified five ligand-receptor pairs that were active between different astrocyte subpopulations and myeloid cells during the process of astrocyte specification, predominantly involving the TGFB3 ligand and its multimeric ACVR1B/TGFBR2 receptor (Figure S8D).

Pseudotime-Dependent Trajectory Analysis Revealed Alterations in Astrocyte Subpopulations and Cuproptosis-Related Genes in Glioma

To further analyze the developmental progress of all astrocyte subpopulations in glioma, we constructed potential developmental trajectories for five astrocyte clusters using the Monocle method (Figure 3A). We observed that three distinct branches could assign pseudotime-dependent progression states to the astrocytes (Figure 3B and C). Astrocyte cluster 1 and a portion of astrocyte cluster 4 were primarily located at the beginning of the developmental trajectory. In contrast, astrocyte cluster 2 and astrocyte cluster 3 were situated at opposite ends of the trajectory. Additionally, we identified a major pathway with distinct branches to validate the potential developmental trajectories of six astrocyte subpopulations in the CGGA dataset (Figure S9A). Astrocyte cluster 2 was predominantly positioned at the start of the timeline trajectory, whereas astrocyte cluster 3 and astrocyte cluster 5 were located at the end of the timeline trajectory (Figures S9B and C).

Figure 3 Pseudotime analysis of different astrocyte subpopulations. (A–C) Distribution of astrocyte cells along the pseudotime trajectory based on astrocyte cell subpopulations: (A) astrocyte cell subpopulations, (B) state, and (C) pseudotime. (D) Heatmap analysis of the expression of cuproptosis-related genes along the pseudotime, categorized by cell cluster. (E and F) Variability in the expression levels of seven genes along the pseudotime (E) and within each astrocyte cell subpopulation (F), based on cell states in our scRNA-seq data. (G and H) Variability in the expression levels of seven genes along the pseudotime (G) and within each astrocyte cell subpopulation (H), based on cell states in the CGGA scRNA-seq dataset. Abbreviation: CGGA, Chinese Glioma Genome Atlas.

In addition, we further analyzed the expression patterns of cuproptosis-related genes along the trajectory of glioma development. The pseudotime heatmap of CRGs derived five distinct clusters of pseudotime patterns (Figure 3D). The expression of cuproptosis genes gradually increased with pseudotime in astrocyte cluster 1, while it gradually decreased with pseudotime in astrocyte clusters 2 and 3. These results illustrate the pseudotime-dependent trajectory of astrocyte subpopulations and the time-dependent changes in cuproptosis-related genes in glioma.

To investigate the overlapping genes between cuproptosis and pseudotime-related genes, we identified seven genes common to both categories: AHR, CARPIN1, CDKN2A, MPC1, REXO2, RPL3, and SOX2 (Figure 3E, F, and Table S3). Furthermore, we validated similar findings in the CGGA single-cell RNA sequencing data, which included CDKN2A, MPC1, and SOX2 (Figure 3G, H, and Table S4). We also observed dynamic changes in the expression patterns of overlapping genes along the trajectory of glioma development and pseudotime in both our in-house and CGGA datasets.

Development and Validation of Cuproptosis-Related Prognosis Signature in Glioma

According to the results of cuproptosis and pseudotime-related genes in our in-house and CGGA datasets, we further constructed a cuproptosis-related signature using three prognosis-related cuproptosis genes. The risk score for CRGs is calculated as follows: CRGs risk score = (−0.09 × expression level of CDKN2A) + (−0.08 × expression level of MPC1) + (0.134 × expression level of SOX2). In the CGGA training dataset, 656 glioma patients with high-risk scores exhibited poorer overall survival (OS) compared to those with low-risk scores (Log rank test, P < 0.001) (Figure 4A, C, and 4D). By incorporating the risk score and other clinical characteristics into univariate and multivariate Cox regression analyses, we identified the risk score as an independent prognostic factor for glioma (Table S5). To verify the applicability and generalizability of the risk score, we further applied it to independent external cohorts from TCGA (n = 628) and obtained similar results (Log rank test, P = 0.037) (Figure 4B, Figure S10A and B). Additionally, we constructed a nomogram that included risk scores and key clinical features. The calibration analysis indicated that the prognostic nomogram is an effective model for predicting the 1-, 3-, and 5-year OS of the 656 glioma patients (Figure 4E and F).

Figure 4 Construction, validation, and evaluation of a risk score model based on cuproptosis-related genes. (A and B) Kaplan–Meier curve analysis for OS of glioma patients, comparing high-risk and low-risk groups in the training (A) and testing (B) sets. (C) Overview of survival in the high-risk and low-risk groups in the training set. (D) Distribution of risk scores among the high-risk and low-risk subgroups in the training set. (E) Nomogram for predicting the 1-, 3-, and 5-year OS of all glioma patients. (F) Calibration curves of the nomogram for predicting the 1-, 3-, and 5-year OS of all glioma patients. (G) Abundance of 22 infiltrating immune cell types in high-risk and low-risk groups. (H) Ranked dot and scatter plots illustrating the correlation between risk score and 22 infiltrating immune cell types in the training dataset. (I and J) Gene sets enriched in the high-risk (I) and low-risk (J) subgroups in the training dataset. *, P < 0.05; **, P < 0.01; ***, P < 0.001. Abbreviations: OS, overall survival.

According to the CIBERSORT algorithm, we further observed significant differences in immune infiltrating cells between high- and low-risk scores, particularly in B cells memory (P < 0.05), T cells follicular helper (P < 0.05), monocytes (P < 0.001), and macrophages M0 (P < 0.05) in the training dataset of 656 glioma patients (Figure 4G and Table S6). Furthermore, we found that B cells memory (P<0.05) and monocytes (P < 0.01) were correlated with risk scores (Figure 4H). Meanwhile, immune infiltrating monocytes that significantly correlated with risk scores were also identified in TCGA testing dataset (P < 0.001) (Figure S10C and S10D). To further elucidate the biological functions associated with the high- and low-risk groups, we conducted a Gene Set Enrichment Analysis (GSEA) and found that the gene sets of the high-risk group were enriched in tumor signaling pathways (Figure 4I), while the gene sets of the low-risk group were enriched in metabolic and immune networks (Figure 4J). Similar results also were validated in the TCGA validation dataset (Figure S10E-S10F). The detailed results of the GSEA are summarized in Table S7.

Experiment Validation for Expression Levels of the Three Cuproptosis-Related Genes

We further evaluated the differential expression levels of the CDKN2A, MPC1, and SOX2 genes through targeted mRNA sequencing in glioma and adjacent normal tissues. Our analysis revealed that CDKN2A (P < 0.05) and SOX2 (P < 0.05) were highly expressed, while MPC1 (P < 0.01) exhibited decreased expression in glioma tissue (Figure 5A). To further investigate the clinical and functional role of these genes, we examined the protein expression levels of tumor and normal tissues based on the Human Protein Atlas database. The expression levels of these three genes (CDKN2A, MPC1, and SOX2) were consistent with our in-house results (Figure 5B).

Figure 5 Experimental validation of three cuproptosis-related genes. (A) mRNA expression levels of the three genes in glioma tissues compared to adjacent normal tissues. (B) Immunohistochemical analysis of the three genes in normal and glioma tissues sourced from the Human Protein Atlas database. *, P < 0.05; **, P < 0.01.

Discussion

Glioma is the most common primary brain tumor and exhibits the highest mortality rate in grade IV glioma – glioblastoma.³ Over the past decades, the prognosis for glioma patients has remained unsatisfactory. Unlike other malignancies, the TME of glioma is highly immunosuppressive, influenced by the characteristics of glial, neuronal, and immune cell populations in the brain.^33,34 Therefore, a deeper understanding of the cellular components within the glioma TME is essential for identifying potential therapeutic targets for glioma patients. Cuproptosis is a recently identified copper induced cell death characterized by copper-dependent intracellular ROS production and excessive accumulation of lipid peroxides.¹⁹ Despite recent increase in the exploration of cuproptosis, a better understanding of cuproptosis in cellular heterogeneity within the TME is still needed. In this study, we comprehensively characterized the role of cuproptosis-related genes in the immune microenvironment of gliomas using single-cell transcriptome profiling and established a cuproptosis-related prognostic model for glioma.

Gliomas exhibit significant cellular and molecular heterogeneity within the TME.³⁵ Illustrating tumor heterogeneity is crucial for the diagnosis and treatment of gliomas.³⁶ Previously established prognostic models related to gliomas have typically relied on bulk RNA sequencing data,^37,38 which reflects the average gene expression levels across entire cell populations. This approach dilutes the contributions of rare cell populations to the overall gene expression patterns³⁹ and fails to precisely illustrate tumor heterogeneity and ignores tumor changes itself in the process of tumor development. Instead, scRNA-seq has proven to enhance our understanding of glioma heterogeneity and can more effectively depict alterations in cell expression patterns and cell-cell interactions.^40,41 Here, we for the first time developed a cuproptosis-related prognostic model by integrating scRNA-seq and bulk mRNA sequencing data to clearly understand the heterogeneity of gliomas, which is critical for uncovering therapeutic targets for glioma patients.

In our in-house and CGGA scRNA-seq datasets, we identified nine distinct cell types that included astrocyte subpopulations. Astrocytes are the most abundant glial cells in the CNS and secrete a variety of immunomodulatory molecules, such as cytokines, chemokines, and metabolites, which play significant roles in inflammation and immune regulation.^42,43 Heiland et al utilized single-cell sequencing data analysis to uncover heterogeneity in the functional responses of astrocyte subpopulations. Their findings suggest that astrocyte subpopulations that facilitate CNS repair and mitigate inflammation may also promote tumor growth and contribute to the generation of an immunosuppressive TME in GBM.¹² Recent studies have explored the roles of astrocyte subpopulations in CNS malignancies and their potential as therapeutic targets.^11,12,44 Similar to macrophages, these reports describe various subpopulations of astrocytes but do not delineate their specific functional roles. Consequently, we systematically identified cuproptosis-related genes in astrocyte-related subpopulations and analyzed their developmental trajectories and intercellular communications.

To further investigate the potential cell cross-talk that may contribute to the development and progression of gliomas, we conducted the analysis of intercellular communications using scRNA-seq. The complex correlations between astrocytes and immune cells shape the TME in gliomas. In the scRNA sequencing data analysis, astrocytes were the most abundant cell types in gliomas with higher relative strength in the TGFβ signaling pathway. The role of TGFβ in glioma biology is multi-dimensional. TGFβ can coordinate the interplay of glioma cell proliferation and migration, immunosuppression, and angiogenesis.⁴⁵ TGFβ plays a crucial role in glioma evasion of host immune responses and can exert immunosuppressive effects on almost all infiltrating immune cells.⁴⁶ Recent studies have shown that TGFβ is closely related to the oncogenic function of gliomas, TGFβ has high activity in GBM, which predicts poorer prognosis of patients.^47,48 Gong et al also revealed that Cu chelator obviously suppresses the activation of TGFβ,⁴⁹ indicating that cuproptosis may be correlated with the development of glioma by regulating the TGFβ signaling pathway.

Current evidence increasingly suggests that cuproptosis may serve as a novel therapeutic target for tumors. Recent studies have demonstrated that investigating cuproptosis in glioma could provide new insights into clinical treatment strategies.^50,51 Bao et al identified and developed a model based on cuproptosis-related differentially expressed genes in low-grade glioma,³⁸ which may represent a promising target for enhancing the efficacy of immunotherapy and cuproptosis treatment in gliomas. However, further research is necessary to determine whether cuproptosis-related genes along the trajectory of astrocyte progression are associated with the prognosis of gliomas. Therefore, we established a cuproptosis-related three-gene signature, comprising CDKN2A, SOX2, and MPC1, which demonstrates a better survival advantage in gliomas with low-risk scores compared to those with high-risk scores. Similar results were also observed in the TCGA testing set. Additionally, GSEA analysis revealed that high-risk gliomas were enriched in tumor signaling pathways, while low-risk gliomas were enriched in metabolic and immune networks. These findings indicate that cuproptosis is closely related to the progression of astrocytes in glioma and the efficacy of treatments. Our CRG-score signature may serve as a reliable indicator for accurately predicting the prognosis of glioma patients and could contribute to guiding personalized immunotherapy for gliomas.

Cyclin-dependent kinase inhibitor 2A (CDKN2A), recognized as a tumor suppressor, inhibits tumor cell proliferation and promotes apoptosis.⁵² Recent studies have demonstrated that CDKN2A is correlated with both ferroptosis⁵³ and cuproptosis,¹⁹ which are characterized by regulated cell death resulting from the deterioration of antioxidant function and the accumulation of lipid ROS. A recent study reported that decreased CDKN2A expression significantly induces tumor recurrence and progression following immune checkpoint blockade therapy in patients with non-small cell lung cancer patients.⁵⁴ Mitochondrial pyruvate carrier 1 (MPC1) is involved in the tumor metabolite and the process of biological synthesis.⁵⁵ Previous study revealed that GBM patients with decreased MPC1 copy numbers segments were associated with poor prognosis and temozolomide resistance.⁵⁶ In addition, Peter et al found that MPC1 gene is a critical mediator of elesclomol-copper ionophore-induced cell death,¹⁹ suggesting the potential role of MPC1 in the induction of cuproptosis in glioma. SOX2 is a known stem cell transcription factor and its upregulated expression is closely correlated with the development and progression of various cancers.⁵⁷ The expression of SOX2 was mediated by the TGFβ signaling pathway, thereby inducing tumor cell self-renewal to prevent the differentiation of glioma-initiating cells.⁵⁸ Moreover, a recent study indicates that SOX2 may be a potential biomarker to predict the sensitivity of cancer cells to ferroptosis, and its high expression can resistant to ferroptosis, which may largely attenuate the potential therapeutic efficacy by inducing ferroptosis.⁵⁹ Similar to ferroptosis, oxidation of SOX2 in cuproptosis could be also a potential therapeutic target for cancer treatment. Previous studies have demonstrated these three genes play an important role in cuproptosis, which provides theoretical support for the reliability of our prognostic model and clinical treatments.

To the best of our knowledge, the present study is the first to systematically integrate single-cell RNA sequencing and bulk RNA sequencing to investigate the prognostic value of cuproptosis-related genes in astrocyte subpopulations within gliomas. Additionally, we comprehensively identify cell subsets and interactions among cells in the glioma TME in our in-house scRNA-seq data and the CGGA scRNA-seq dataset. Based on the similarity of pseudo-time-dependent trajectories of astrocyte subsets in gliomas and the time-dependent changes of cuproptosis-related genes, we developed and validated a novel cuproptosis-related prognostic risk model, which demonstrates significant predictive value. However, we acknowledge some potential limitations. Firstly, our prognostic model was constructed based on the CGGA dataset, which may introduce selection bias; nonetheless, we validated this model using TCGA dataset, confirming its reliability to a certain extent. Furthermore, our targeted mRNA sequencing also corroborates the involvement of these three genes. Secondly, our study is primarily descriptive, and experimental validation in vitro and in vivo of these genes is necessary to enhance our understanding of cuproptosis in gliomas in the future. Nevertheless, our study contributes to the understanding of cellular heterogeneity, intercellular communication, astrocyte subpopulations, and cuproptosis in gliomas.

Conclusion

In summary, we investigated the tumor microenvironment in gliomas and the characteristics of cuproptosis in astrocyte subsets using scRNA-seq data from three in-house glioma patients and 14 glioma patients from the CGGA dataset. We developed and validated a novel cuproptosis-related prognostic signature consisting of three genes (CDKN2A, SOX2, and MPC1) in a total of 1284 glioma patients from the CGGA and TCGA datasets. These results highlight that cuproptosis-related prognostic signature in unique astrocyte subpopulations for glioma patients, which may provide valuable insights for clinical decision-making and guide the development of appropriate therapeutic strategies.