Integrative analysis of crotonylation-associated genes reveals prognostic and therapeutic targets in gliomas
Bowen Yin, Zhuoyang Fan, Panpan Yu, Jin Li, Yilin Wang, Minfeng Shu
TL;DR
This study explores how crotonylation, an epigenetic modification, affects glioma progression and identifies genes that could help predict outcomes and guide treatment.
Contribution
The study introduces a novel risk score based on crotonylation-associated genes for glioma prognosis and treatment stratification.
Findings
Crotonylation levels were negatively correlated with glioma grade and patient survival outcomes.
High-risk patients showed increased mutations and an immunosuppressive tumor environment.
Reduced crotonylation suppressed CXCL1 and promoted GZMB expression in glioma cells.
Abstract
Crotonylation, an emerging epigenetic modification, has been implicated in various biological processes, including tumor progression. However, its role in glioma remains poorly understood. This study aims to investigate the prognostic and therapeutic implications of crotonylation-associated genes in glioma. Crotonylation levels were assessed by IHC in glioma tissues of varying grades. Key crotonylation-associated genes were identified and analyzed across five glioma datasets. A prognostic risk score was developed using machine learning algorithms and validated in multiple cohorts. Genomic alterations, immune landscapes, and therapeutic responses were examined in relation to the risk score. Single-cell dataset GSE131928 was analyzed to explore the relationship between the risk score and immune cell infiltration. After crotonate treatment of T98G cells, ChIP-seq and qPCR were performed…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsGlioma Diagnosis and Treatment · Cancer, Hypoxia, and Metabolism · ATP Synthase and ATPases Research
Introduction
Gliomas, the most common and aggressive primary brain tumors, are characterized by remarkable heterogeneity in their genetic, epigenetic, and immune profiles, posing significant challenges for effective treatment and prognosis (1–3). Among gliomas, glioblastoma (GBM) is the most malignant subtype and remains refractory to current therapeutic modalities (4). The standard treatment, known as the Stupp protocol, combines maximal surgical resection with radiotherapy and concurrent temozolomide (TMZ) chemotherapy (5). Despite these interventions, the prognosis for GBM patients remains dismal, with a median survival of approximately 14–18 months and a high likelihood of recurrence (6). The advent of immunotherapy, which has revolutionized the treatment landscape for cancers such as melanoma and lung cancer (7, 8), has also sparked interest in developing immune-based strategies for glioma. However, challenges such as the immunosuppressive tumor microenvironment and low immunogenicity of gliomas necessitate further exploration of novel biomarkers and therapeutic approaches to improve patient outcomes (9, 10).
Recent studies have uncovered crotonylation, a novel post-translational modification, as a crucial regulatory mechanism in cancer biology (11, 12). This dynamic acylation modification occurs on lysine residues of histones and non-histone proteins, influencing chromatin structure and transcriptional activity (13). Unlike acetylation, crotonylation is associated with transcriptionally active chromatin and is tightly linked to cellular metabolism, particularly the production of Crotonyl-CoA, a key substrate for this modification (14). Crotonylation has been implicated in diverse cellular processes, including DNA repair, cell proliferation, and immune responses, suggesting its potential role in oncogenesis (15, 16). In glioma, while other epigenetic modifications like methylation and acetylation have been extensively studied, the specific contribution of crotonylation to tumor progression, immune regulation, and therapy resistance remains underexplored. Understanding the role of crotonylation-associated genes and their regulatory networks in glioma may uncover novel biomarkers and therapeutic targets, offering new avenues for treatment (17, 18).
We systematically evaluated crotonylation-associated genes in gliomas, focusing on their prognostic and therapeutic significance. By analyzing multi-omics datasets from large glioma cohorts, we constructed a prognostic risk score based on crotonylation-related genes and assessed its relationship with genomic alterations, the immune microenvironment, and therapeutic response. This work highlights the potential of crotonylation-associated genes as biomarkers for glioma stratification and treatment, providing a foundation for future epigenetic studies in gliomas.
Result
Identification and prognostic impact of crotonylation-associated gene clusters in glioma
We first examined the levels of crotonylation modification in human glioma samples of varying grades (D79717: grade 2; D88307: grade 3; E22087: grade 4). The results revealed a progressive decline in overall crotonylation levels with increasing glioma grade (Figure 1A). To identify genes associated with crotonylation, we focused on enzymes involved in the production of Crotonyl-CoA and those catalyzing histone and non-histone crotonylation modifications. Key genes identified “included” GCDH and ECHS1, both of which play crucial roles in lysine metabolism. Specifically, GCDH catalyzes the conversion of glutaryl-CoA to Crotonyl-CoA, while ECHS1 facilitates the hydration of Crotonyl-CoA to produce 3-hydroxybutyryl-CoA (19–21). Additionally, ACSS2, ACADS, ACOX1, and ACOX3 are essential in short-chain fatty acid metabolism, where ACSS2 catalyzes the formation of Crotonyl-CoA from crotonate, and ACADS, ACOX1, and ACOX3 catalyze the conversion of butyryl-CoA to Crotonyl-CoA (19, 22, 23). Reported crotonylation “writers” include EP300, CREBBP, KAT8, and KAT2B, while “erasers” include HDAC1, HDAC2, HDAC3, SIRT1, SIRT2, and SIRT3. Recognized “readers” of crotonylation modifications include DPF2, MLLT3, and YEATS2 ( 11, 12, 24–32) (Supplementary Table S1).
Consensus clustering of crotonylation profiles in glioma. (A) IHC staining of Kcr in glioma samples with different grades (G2, G3, G4): D79717, D88307, and E22087. (B) PCA of raw expression profiles for the combined dataset. (C) PCA after COMBAT adjustment for combined expression profiles. (D) Unsupervised clustering of 20 crotonylation-related genes reveals two distinct clusters with patient annotations including cohort source, survival status, gender, and age. (E) PCA indicating two distinct clusters within the combined cohorts. (F) Kaplan-Meier analysis estimating overall survival differences between Cluster-A and Cluster-B.
To assess the impact of these genes on glioma patient outcomes, we analyzed five glioma datasets: TCGA (GBM + LGG), CGGA1, CGGA2, GSE42669, and GSE7696. Batch effects were minimized across these datasets using the “ComBat” function in the “sva” R package, resulting in reduced inter-dataset variability, thus enabling cohesive analysis (Figures 1B, C). To further explore the regulatory role of crotonylation-related genes in relation to glioma patient survival, we applied the ConsensusClusterPlus package in R. Unsupervised clustering revealed several distinct clusters, all with comprehensive survival data. Importantly, two independent clusters, labeled A and B, were identified based on a consensus matrix with k = 2 (Supplementary Figure S1; Figure 1D). Principal component analysis (PCA) indicated a significant distinction between these two clusters (Figure 1E). Survival analysis demonstrated that patients in cluster B were associated with significantly improved prognosis across TCGA (log-rank test, P < 0.001), CGGA1 (log-rank test, P = 0.041), CGGA2 (log-rank test, P < 0.001), GSE7696 (log-rank test, P = 84), GSE42669 (log-rank test, P = 0.011) and the combined cohort (log-rank test, P < 0.001) (Figure 1E; Supplementary Figure S2).
Construction of the risk score in glioma
To quantify the relationship between crotonylation-related genes and glioma patient survival, we developed a risk score system tailored to glioma patients. crotonylation-associated genes were incorporated into this system using a variety of machine learning algorithm combinations. Following the approach previously reported (33), we integrated 10 machine learning algorithms—including CoxBoost, Stepwise Cox, Ridge, RSF, GBM, Survival-SVM, Lasso, Enet, plsRcox, and SuperPC—to achieve a risk score with high accuracy and stability across different cohorts. In the TCGA-glioma cohort, over 100 predictive models were constructed, and the average C-index of each model was calculated across four independent validation cohorts (Supplementary Figure S3). Among these models, the combination of RSF and StepCox algorithms demonstrated the highest average C-index in the validation cohorts. Genes such as SIRT, MLLT3, CREBBP, EP300, KAT2B, GCDH, SIRT2, YEATS2, DPF2, SIRT3, KAT8, ECHS1, and ACOX1 were negatively correlated with the risk score, while others exhibited positive correlations (Figure 2A).
Construction of the risk score in glioma. (A) Heatmap showcasing the expression levels of 20 crotonylation-related genes alongside their associated risk scores in combined glioma cohorts. A bar chart to the left displays the correlations between these genes and the risk score. (B) Concordance index evaluating the predictive accuracy of the risk score across various cohorts. (C) Univariate Cox regression analyses depicting the association between the risk score and survival across different cohorts. (D) Time-dependent Area Under the Curve for the risk score across different cohorts. (E) Kaplan-Meier survival analysis comparing overall survival rates between low- and high-risk groups in a combined-cohort.
The C-index of the risk score was also assessed across individual cohorts, where it showed strong predictive performance (Figure 2B). Univariate Cox regression analyses confirmed the risk score as an independent prognostic biomarker for overall survival in glioma patients (Figure 2C). Furthermore, the time-dependent area under the curve (AUC) values underscored the risk score’s utility as a prognostic indicator in the TCGA, CGGA, and GEO datasets (Figure 2D). Based on the median risk score across combined cohorts, glioma patients were categorized into high- and low-risk groups. Kaplan–Meier survival analysis revealed significantly poorer overall survival for patients in the high-risk group across all cohorts (Figure 2E; Supplementary Figure S4). Collectively, these findings indicate that the risk score is both stable and robust across diverse independent cohorts, making it a valuable tool for glioma prognosis.
Genomic status of different risk groups
To characterize the genomic landscape of different risk groups within the TCGA-glioma dataset, we first examined the frequency of somatic mutations. A positive correlation was identified between the risk score and somatic mutation count, indicating that the high-risk group exhibited a higher frequency of somatic mutations, encompassing both synonymous and non-synonymous types (Figures 3A–F). To further explore the relationship between genomic alterations and risk score, we analyzed copy number variations (CNVs) across different risk groups. In the high-risk group, genes on chromosomes 7 and 12 showed a tendency toward amplification, while genes on chromosomes 1, 8, and 13 were more prone to deletion (Figure 3G). Genomic instability in gliomas drives somatic mutations and chromosomal copy number variations (CNVs), both of which shape tumor evolution and therapeutic responses. In high-risk gliomas, recurrent CNVs—such as chromosome 7/12 amplifications (EGFR, CDK4/MDM2) and chromosome 1/8/13 deletions (CHD5, DLC1, RB1)—play pivotal roles in promoting proliferation, immune evasion, and DNA repair defects. For instance, EGFR amplifications activate PI3K/AKT signaling to fuel tumor growth (3), while CDK4 gains bypass cell cycle checkpoints, and RB1 losses impair homologous recombination repair (3, 34, 35). These CNVs synergize with somatic mutations to amplify neoantigen burden, paradoxically creating immunogenic tumors that remain therapy-resistant due to compensatory immunosuppressive mechanisms.
Genomic states of different risk groups in glioma. (A-C) Boxplots illustrating comparisons of total mutation counts (A), synonymous mutation counts (B), and non-synonymous mutation counts (C) between low- and high-risk groups within the TCGA-glioma cohort. (D-F) Analysis of the correlation between risk scores (log transformed) and mutation counts for total mutations (D), synonymous mutations (E), and non-synonymous mutations (F) in the TCGA-glioma cohort. (G) Overview of copy number variations, showing gains and losses in groups categorized by low and high risk. (H) Boxplots detailing the differences in neoantigen counts between low- and high-risk groups in the TCGA-glioma cohort.
Neoantigens are critical determinants of tumor immunogenicity. In this study, we utilized a neoantigen gene set, defined by mutations capable of generating neoepitopes, to investigate immune response characteristics in glioma. As hypothesized, the high-risk group displayed a significantly higher neoantigen count compared to the low-risk group (Figure 3H). These findings indicate that an elevated risk score is linked to an increased mutation burden in glioma tumor cells, potentially contributing to a greater neoantigen load.
Distinct tumor microenvironment characteristics in high- and low-risk glioma groups
To characterize the biological features of tumors under the risk classification system, Gene Set Enrichment Analysis was conducted to determine pathways that were associated with differentially expressed genes in high- and low-risk groups. Initially, differential gene expression analysis between the high- and low-risk groups was performed using three R packages: DESeq2, edgeR, and limma (Supplementary Figure S5). The results from each package were then analyzed separately to identify the top 10 enriched pathways for each group. In the low-risk group, pathways such as adrenergic signaling in cardiomyocytes, aldosterone synthesis and secretion, cardiac muscle contraction, cortisol synthesis and secretion, long-term depression, and nicotine addiction were predominantly enriched. In contrast, the high-risk group showed enrichment of pathways associated with IL-17 signaling, inflammatory bowel disease, intestinal immune network for IgA production, and viral protein interaction with cytokine and cytokine receptors (Figures 4A, B; Supplementary Figure S6).
Identification of two distinct patterns in high-risk and low-risk sample groups. (A) GSEA analysis showing pathway enrichment of differentially expressed genes analyzed by the limma package in the low-risk group. (B) GSEA analysis showing pathway enrichment of differentially expressed genes analyzed by the limma package in the high-risk group. (C) Analysis of hallmark gene sets in both low- and high-risk groups.
To further investigate the cancer-related pathways underlying these groups, a heatmap was generated using hallmark gene sets from the msigdbr R package, revealing significant expression differences between high- and low-risk groups (Figure 4C). In the high-risk group, pathways related to immune and inflammatory responses—such as interferon-gamma response, IL6-JAK-STAT3 signaling, IL2-STAT5 signaling, inflammatory response, and TNF-alpha signaling via NF-kB—were markedly upregulated. In contrast, genes associated with pancreatic beta-cell function, hedgehog signaling, and KRAS signaling were highly expressed in the low-risk group. The methylation status of the MGMT (O6-methylguanine-DNA methyltransferase) promoter is a crucial biomarker in gliomas, particularly in predicting the response to temozolomide (TMZ) chemotherapy. Methylation silences MGMT expression, thereby impairing its ability to repair DNA damage induced by TMZ, leading to better therapeutic outcomes (36, 37). IDH (isocitrate dehydrogenase) mutations, predominantly IDH1 and IDH2, are another critical molecular hallmark, commonly associated with lower-grade gliomas and better prognoses. IDH mutations result in the production of oncometabolites that alter cellular metabolism and epigenetic states, including histone and DNA modifications (38, 39). We observed that patients in the low-risk group were more likely to have a methylated MGMT promoter, while patients in the high-risk group were more frequently characterized by IDH mutations. Moreover, glioblastoma patients, representing high-grade gliomas, were almost exclusively classified into the high-risk group. These findings suggest that the risk scoring system effectively differentiates distinct tumor microenvironments in glioma, highlighting the unique biological characteristics associated with each risk group.
Immune landscape of different risk groups in glioma
From the preceding analysis, we observed that the high-risk group exhibited a greater number of neoantigens and significant enrichment in immune response-related pathways. To further characterize the immune status within different risk groups, we employed four distinct methods—ESTIMATE, CIBERSORT, xCELL, and MCPcounter—to evaluate immune cell infiltration. Using the ESTIMATE algorithm, we found that the high-risk group had significantly higher stromal, immune, and ESTIMATE scores compared to the low-risk group (Figures 5A–C). Overall, gliomas exhibited limited infiltration of lymphocytes such as CD8+ T cells, CD4+ T cells, and B cells, as well as antigen-presenting cells like dendritic cells, consistent with the generally immunosuppressive nature of gliomas. Comparing the high- and low-risk groups, the high-risk group demonstrated increased infiltration of relatively immunosuppressive cell types, including macrophages (M0, M1, and M2 types), Th2 cells, fibroblasts, monocytes, and neutrophils. Conversely, the low-risk group showed a higher infiltration of NK cells (Figure 5D), reflecting distinct immune microenvironments between the two groups.
Immune landscape of different risk groups in glioma. (A-C) The correlation between the risk score and stromal score (A), immune score (B), and ESTIMATE score (C). (D) Heatmap showing the infiltration of immune cells between low- and high-risk groups in the TCGA-glioma cohort. (E) Regulation of Immunomodulators, illustrating from left to right: mRNA expression, expression versus methylation, amplification frequency, and deletion frequency (represented as amplifications) for 75 immunomodulatory genes by immune subtype. The colored blocks above, ranging from red, yellow, green, cyan, to blue from left to right, represent the division of the risk score into five equally spaced groups, ordered from smallest to largest. (F) Assessment of the risk score correlation with pathways related to anti-tumor immune response steps in TCGA glioma. (G) Correlations between the expression of risk score and immune checkpoint markers at mRNA levels.
Given that the risk score was constructed based on crotonylation-associated genes, we further explored the connection between crotonylation and immune-related regulation. Crotonylation, like DNA methylation, acts as an epigenetic regulatory mechanism, with both often working in concert to influence gene expression. Although DNA methylation and histone crotonylation operate on different molecular substrates (DNA versus histones), they are closely interlinked in gene regulation. DNA methylation typically recruits methyl-binding proteins (such as MeCP2) that subsequently recruit HDACs, leading to the removal of histone crotonylation. Since histone crotonylation is associated with an open chromatin structure, its removal results in tighter chromatin conformation, further repressing gene expression (40, 41). To investigate the role of immune-related molecules, we analyzed co-inhibitors, co-stimulators, ligands, receptors, cell adhesion molecules, and antigen presentation factors in glioma. We divided the TCGA-glioma samples into five groups based on their risk scores. Most of these molecules, except for HLA, exhibited relatively low expression levels in glioma. To account for this low expression potentially obscuring their intrinsic trends, we also examined their methylation levels and mutation status, finding that alterations in methylation and mutations were amplified to varying degrees (Figure 5E). In particular, IL-13, C10orf54, ENTPD1, TNFSF4, and HLA-DR showed significant correlations with the risk score. Further correlation analysis of expression levels with the risk score revealed a positive association for TNFSF4, HLA-DRA, and ENTPD1(Supplementary Figure S7). Additionally, we observed that the risk score negatively correlated with key processes such as T-cell priming and activation, recognition of cancer cells by T cells, and killing of cancer cells (Figure 5F). However, it positively correlated with immune checkpoint molecules (Figure 5G), suggesting a distinct immune profile linked to the risk score in glioma.
Predictive value of the risk score in immunotherapy and chemotherapy
To predict immune checkpoint blockade (ICB) response, the TIDE score, a widely recognized predictor, is commonly used to evaluate immune response. The risk score was positively correlated with TIDE and exclusion, while negatively correlated with dysfunction and microsatellite instability (MSI Expr sig) (Supplementary Figure S8A). To understand the impact of the risk score on the clinical efficacy of glioma treatments, we analyzed correlations between the risk score and the IC50 of candidate drugs in the Genomics of Drug Sensitivity in Cancer (GDSC) database (Supplementary Figure S8B). Over 50 candidate drugs (|Rs| > 0.5) were identified, most of which exhibited a positive correlation between IC50 and the risk score. These drugs targeted IGF-1R inhibitors, CDK4 and CDK6 inhibitors, EGFR-TK inhibitors, and other AMPK inhibitors. Additionally, six drugs exhibited a negative correlation between IC50 and the risk score, including mTOR inhibitors, Lck inhibitors, and Src/Abl inhibitors (Supplementary Figures S8C, D).
In summary, these findings suggest that the risk score plays a crucial role in mediating immune response and is associated with drug sensitivity. Thus, the risk score could serve as a potential biomarker for establishing appropriate treatment strategies.
Exploring the role of CXCL1 and tumor-associated macrophages
To achieve a higher resolution in detecting the relationship between risk scores and immune cell infiltration, we conducted a series of analyses on the single-cell dataset GSE131928, which includes data from nine glioma patients. First, we clustered the cells based on malignancy, categorizing them into malignant cells, immune cells, and other cells (Figure 6A). Then, using SingleR, we further subdivided the non-malignant cell populations and compared the immune cell subgroup distributions across the nine patients, with cluster identities validated by canonical marker genes (Supplementary Figure S9). We found that Tumor-Associated Macrophages (TAMs) constituted a relatively large proportion of the cell populations (Figure 6B). After calculating the risk scores for each patient, we discovered a positive correlation between the risk score and the infiltration of TAMs (Figure 6C), thereby validating our previous results at the single-cell level.
Exploring the role of CXCL1 and tumor-associated macrophages. (A) UMAP plot illustrating cell clustering in the single-cell dataset GSE131928. (B) UMAP plot displaying immune cell subgroup distributions across different patients. (C) Boxplot showing the proportion of different immune cell subgroups ordered by risk score from low to high for each patient in the GSE131928 dataset. (D) Schematic representation of the ChIP-seq procedure after crotonate treatment of T98G cells. (E) KEGG enrichment analysis highlighting peak differences between crotonate-treated T98G and control T98G cells. (F) MA plot depicting the distribution of peaks in control versus crotonate-treated groups. (G) UMAP plot demonstrating the relative expression of CXCL1 across different subgroups, with malignant cells circled by dashed lines. (H) qPCR bar graph comparing CXCL1 mRNA levels between control and crotonate-treated T98G cells (ND, Not Detected). (I) Immunofluorescence staining of PD-1 and GZMB in tissue samples E22087 and D79717. PD-1 is shown in yellow, GZMB in pink, and nuclei are stained with DAPI in blue.
According to reports, crotonylation modification can occur on histone lysine residues, thereby regulating the transcription levels of DNA in those regions (13). To further determine the relationship between TAMs and risk scores, we conducted a series of experiments using the human glioma cell line T98G. Crotonate, a crotonyl-CoA precursor, has been reported to increase intracellular crotonylation modification levels. After treating the cells with crotonate for 48 hours, we performed ChIP-seq analysis using a pan-kcr antibody (Figure 6D). KEGG pathway analysis showed that, compared to the crotonate-treated T98G cells, the peaks of chemokine signaling pathway was significantly upregulated in control group (Figure 6E), with CXCL1 exhibiting the greatest difference (Figure 6F). Umap plots indicated that CXCL1 is relatively high expressed in the malignant cell populations within gliomas (Figure 6G). Subsequently, we detected the mRNA levels of CXCL1 using qPCR. The results showed a significant reduction in CXCL1 expression in crotonate-treated T98G cells compared to the control group (Figure 6H). CXCL1 is a chemokine that mediates the migration of macrophages to tumor sites through binding with its receptor CXCR2 (42). Tumor cells and surrounding stromal cells secrete CXCL1, creating a chemotactic gradient that attracts monocytes from the circulation into the tumor tissue, where they further differentiate into TAMs. Within the TME, CXCL1 not only participates in the recruitment of macrophages but also influences their polarization (43). TAMs typically exhibit an M2-like phenotype, possessing pro-tumor, pro-angiogenic, and immunosuppressive functions. CXCL1 promotes the polarization of macrophages towards the M2 phenotype by regulating signaling pathways such as PI3K/Akt and STAT3 (44), thereby supporting tumor growth and metastasis. The interaction between CXCL1 and TAMs also involves immunosuppressive mechanisms. M2-TAMs secrete inhibitory cytokines such as IL-10 and express immune checkpoint molecules like PD-L1, which suppress the activity of effector T cells, helping tumors evade immune system surveillance and attack (45). Therefore, we examined the expression levels of PD-1 and GZMB in glioma samples with different levels of crotonylation modification. Consistent with our previous results, PD-1 expression was negatively correlated with crotonylation modification levels, while GZMB expression was positively correlated with crotonylation modification levels (Figure 6I).
Discussion
This study systematically investigates the prognostic and therapeutic significance of crotonylation-associated genes in glioma, revealing their critical role in tumor progression, immune regulation, and patient outcomes. By focusing on the expression and functional impact of crotonylation-related genes across large glioma cohorts, we constructed a robust risk score that effectively stratifies glioma patients into distinct prognostic groups. The results highlight crotonylation as a pivotal epigenetic modification influencing glioma biology, providing a foundation for potential therapeutic interventions.
Crotonylation is a post-translational modification first identified in histones, characterized by the addition of a crotonyl group to lysine residues (13). In the context of cancer, crotonylation has emerged as a significant epigenetic marker influencing tumor biology. Studies have shown that crotonylation plays a role in cancer cell proliferation, invasion, and immune evasion. For instance, crotonylation modifications have been implicated in colorectal cancer and liver cancer, where dysregulated metabolic pathways alter the availability of crotonyl-CoA, affecting chromatin dynamics and gene expression (41, 46, 47). Despite these advancements, research on crotonylation in gliomas remains sparse. Gliomas, particularly glioblastoma multiforme, are known for their heterogeneity and poor prognosis, yet the role of crotonylation in regulating their complex tumor microenvironment and progression is largely unexplored. This paucity of data underscores the need for studies like ours to investigate the potential regulatory mechanisms and therapeutic implications of crotonylation in glioma biology.
The crotonylation-associated risk score in this study was developed using multiple machine learning approaches to identify the optimal combination of prognostic genes (33). Methods such as random forest analysis, and stepwise Cox regression were employed to select genes with the most significant contribution to patient outcomes. Each method offered unique advantages: random forest analysis identified feature importance through ensemble learning, and stepwise Cox regression refined the model to ensure that only the most independent predictors were included. The use of combined machine-learning techniques allowed for a more nuanced and data-driven selection process, capturing complex relationships between gene expression and survival outcomes. Additionally, this combination of methods ensured that the risk score achieved a balance between simplicity and predictive accuracy.
We revealed a significant relationship between the crotonylation-associated risk score and neoantigen burden in glioma patients, suggesting that crotonylation-related genes may influence tumor immunogenicity. High-risk patients exhibited a higher neoantigen load, which is typically associated with increased tumor immunogenicity. However, despite this elevated neoantigen burden, the immune response in these tumors appeared to be suppressed, likely due to the immunosuppressive tumor microenvironment. Gene set enrichment analysis (GSEA) further supported these findings, as immune-related pathways, including those involved in T cell activation, antigen presentation, and cytokine signaling, were significantly enriched in the high-risk group.
Furthermore, our study revealed the predictive value of the risk score in immunotherapy and chemotherapy. High-risk patients exhibited higher TIDE scores, indicating resistance to immune checkpoint blockade (ICB) therapy. This resistance was supported by the positive correlation between the risk score and immune exclusion signatures. Further drug prediction analysis identified potential therapeutic agents with distinct efficacy profiles for high- and low-risk gliomas. These findings emphasize the utility of the risk score in guiding treatment decisions, enabling personalized therapeutic strategies for glioma patients.
Our integrated single-cell and bulk sequencing analyses revealed a significant correlation between tumor-associated macrophages (TAMs) and the crotonylation-associated risk score. TAMs, predominantly polarized into the M2 phenotype, drive tumor progression through multifaceted mechanisms, including immune suppression, hypoxia induction, and metastasis promotion (48). ChIP-seq analysis demonstrated that crotonate treatment globally elevates crotonylation levels but paradoxically reduces peaks at specific loci, such as the chemokine CXCL1. We hypothesize: first, although crotonate treatment globally upregulated crotonylation levels, different genes may respond differently to crotonylation modification. Some genes might, in a high crotonylation environment, regulate their crotonylation status through feedback mechanisms or the action of specific de-crotonylases (such as HDACs), leading to a decrease in their crotonylation peaks. Second, crotonylation modifications can affect chromatin openness and its interactions with various epigenetic regulatory proteins. Crotonate treatment may alter the chromatin structure of specific gene regions, reducing their sensitivity to crotonylation modification or changing the binding affinity of transcription factors associated with these modifications, thereby causing a decrease in the crotonylation peaks of these genes. CXCL1, suppressed under high crotonylation conditions, is critical for recruiting monocytes via CXCR2 binding, a process further amplified in hypoxic niches by damage-associated molecular patterns like HMGB1 (49–51). Malignant cells secrete M2-polarizing cytokines (e.g., IL-10, CCL5, CXCL12) to drive their differentiation into immunosuppressive M2 TAMs (52), which in turn secrete IL-10, TGF-β, and HLA-G while recruiting myeloid-derived suppressor cells (MDSCs). MDSCs suppress T-cell function through arginine depletion via arginase-1/iNOS and amplify immunosuppression by producing IL-10, which skews macrophages toward the M2 phenotype (53). M2 TAMs further enhance immune evasion by upregulating PD-L1 and B7-H4, immune checkpoints induced through direct tumor cell contact and correlated with poor prognosis in glioblastoma and hepatocellular carcinoma (54, 55). These checkpoints inhibit cytotoxic T-cell activity, as evidenced by restored CD4+/CD8+ T-cell infiltration following TAM depletion in preclinical models. Importantly, crotonylation is metabolically regulated by lysine-derived crotonyl-CoA, positioning lysine metabolism as a therapeutic lever to modulate epigenetic-immune crosstalk. Dietary or pharmacological interventions targeting lysine intake or HDAC activity could recalibrate crotonylation dynamics, potentially disrupting TAM recruitment and polarization to enhance immunotherapy efficacy in glioma. In parallel, our ongoing development of crotonylation-enhancing therapeutic peptides represents a complementary strategy to directly boost intracellular crotonyl-CoA production through targeted activation of metabolic enzymes. These two approaches, though mechanistically distinct, converge on the common goal of elevating global crotonylation levels – a biochemical modification that appears to exert tumor-suppressive effects through epigenetic reprogramming of both malignant cells and immune components within the glioma ecosystem.
Additionally, our results identified specific immune-related molecules—TNFSF4, HLA-DRA, and ENTPD1—that were significantly correlated with the crotonylation-associated risk score, further linking crotonylation pathways to the regulation of the immune microenvironment in gliomas. Mechanistically, crotonylation-mediated chromatin remodeling alters transcriptional programs, ultimately reshaping the tumor immune landscape and modulating immune exhaustion-related pathways. Among these molecules, ENTPD1 (CD39) stands out due to its well-documented role in tumor immunosuppression. ENTPD1 is an ectonucleotidase that hydrolyzes extracellular ATP into AMP, a precursor of immunosuppressive adenosine. This enzymatic activity contributes to the establishment of an immunosuppressive microenvironment by reducing pro-inflammatory ATP signaling and increasing adenosine-mediated T cell suppression (56–58). Elevated expression of ENTPD1 has been associated with poor prognosis in multiple cancers, including gliomas, where it promotes T cell exhaustion and supports the expansion of regulatory T cells (Tregs). Subsequently, we examined the correlation between Kcr levels and glioma grade in human glioma samples and observed a correlation between Kcr and the expression levels of PD-1 and GZMB in the glioma microenvironment.
Despite these significant findings, our study has some limitations. Gliomas are characterized by profound intra- and intertumoral heterogeneity, which manifests at genetic, epigenetic, and cellular levels. While our risk score effectively stratifies patients into prognostic groups using bulk transcriptomic data, this approach inherently averages subclonal variations. For instance, spatially distinct tumor regions may harbor subclones with divergent crotonylation profiles, potentially leading to underestimation of high-risk subpopulations in bulk analyses. Such heterogeneity could influence the reliability of the risk score in predicting outcomes for patients with mixed molecular subclones. Furthermore, gliomas with heterogeneous IDH mutation status or MGMT promoter methylation—critical molecular subtypes—may exhibit varying crotonylation dynamics. Finer stratification across subtypes could refine the score’s precision in future studies. First, functional studies are needed to elucidate the mechanistic roles of crotonylation-related genes in glioma. Second, the connection between crotonylation and specific immune processes warrants further exploration to fully understand its impact on glioma immunotherapy. Lastly, while the predictive performance of the risk score was robust, additional prospective studies are necessary to confirm its clinical utility. To overcome these limitations, we propose integrating single-cell RNA sequencing and spatial transcriptomics in subsequent research. scRNA-seq would enable resolution of crotonylation-associated gene expression at the single-cell level, identifying rare subclones or immune-stromal interactions masked in bulk data. Spatial techniques could map crotonylation patterns across tumor regions, clarifying their relationship to microenvironmental niches. Additionally, multi-region sampling of tumors in validation cohorts would help quantify regional variability in crotonylation and its prognostic relevance. Finally, validating the risk score within homogeneous molecular subgroups could enhance its clinical applicability.
In conclusion, this study underscores the critical role of crotonylation-associated genes in glioma progression, prognosis, and therapeutic response. The risk score based on these genes provides a valuable tool for stratifying glioma patients and tailoring treatment strategies. Further research into the mechanisms underlying crotonylation and its interplay with the TME may pave the way for novel therapeutic approaches in glioma management.
Method
Datasets
RNA sequencing datasets encompassing glioma patients were sourced from multiple databases. For TCGA cohorts, which included Lower Grade Glioma (LGG; WHO grade II–III) and Glioblastoma Multiforme (GBM; WHO grade IV), data were retrieved from the UCSC Xena platform (https://xenabrowser.net/datapage/). Accompanying clinicopathological details were acquired from the cBioPortal (http://www.cbioportal.org/). Additionally, RNA sequencing and clinical data for patients from the Chinese Glioma Genome Atlas (CGGA) were obtained from the CGGA’s official website (https://www.cgga.org.cn). Further, two RNA sequencing datasets (GSE42669, GSE7696 and GSE131928) along with their clinical information were downloaded from the Gene Expression Omnibus (GEO) database. Any patient records lacking survival data were systematically excluded from subsequent survival analyses. Pan-Kcr ChIP-seq data are accessible at GEO under accession codes GSM8532192, GSM8532193, and GSE297117.
Data processing and mutation analysis
Batch effects across different cohorts were corrected using the “ComBat” algorithm within the “SVA” package. Somatic mutation and copy number alteration (CNA) data for glioma were obtained from TCGA. Additionally, copy number variations associated with glioma were assessed using GISTIC 2.0 through the GenePattern platform.
Construction of the risk score
To construct the risk score, we implemented a robust methodological framework, evaluating 117 algorithmic combinations to optimize the predictive utility of crotonylation-related genes. This approach incorporated ten distinct algorithms: Lasso, Ridge, Stepwise Cox, CoxBoost, Random Survival Forest (RSF), Elastic Net (Enet), Partial Least Squares Regression for Cox (plsRcox), Supervised Principal Components (SuperPC), Generalized Boosted Regression Modeling (GBM), and Survival Support Vector Machine (Survival-SVM). Each predictive model was subjected to rigorous validation across multiple cohorts, with the Concordance Index (C-index) calculated to quantify predictive accuracy (Supplementary Table S2). The final crotonylation risk signature was derived using the RSF and Stepwise Cox algorithms, which yielded the highest average C-index, affirming their superior predictive reliability across the validation cohorts.
Prediction of response to immunotherapy and chemotherapy
To explore various tumor immune evasion mechanisms, we employed the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm. Chemotherapy response predictions were generated for each sample using the “pRRophetic” package, leveraging a ridge regression model based on cancer drug sensitivity genomics. Additionally, Spearman correlation and differential analyses were performed across risk groups to identify promising therapeutic agents.
Functional and pathway enrichment analysis
Gene Ontology and Kyoto Encyclopedia of Genes and Genomes analyses were conducted to explore the functions of differentially expressed genes, with significance criteria set at a corrected p-value < 0.05 and an absolute log fold change > 1. Enrichment analysis was performed using the “clusterProfiler” R package (v4.8.3), and results were visualized with the “ggplot2” R package (v3.4.4).
Estimation of TME cell infiltration
Tumor microenvironment characteristics, including immune, stromal, and ESTIMATE scores for tumor samples, were assessed using the “ESTIMATE” R package. Immune cell infiltration levels within the TME of glioma were quantified using the CIBERSORT, xCell, and MCPcounter algorithms.
Single-cell data analysis
The single-cell dataset GSE131928 was downloaded from the GEO database. Subsequent quality control, dimensionality reduction, and normalization were performed following the steps outlined on the Seurat website. Cell clustering was conducted using the SingleR package.
Quantitative real-time PCR
Total RNA of T98G cells was extracted with Trizol reagent (Invi- trogen, USA) following the manufacturer’s instructions.Reverse tran- scription was performed using the 1st strand cDNA Synthesis kit (Yeason, China).quantitative real-time PCR were performed using the Hieff^®^ qPCR SYBR^®^ Green Master Mix (Yeason, China). The relative expression of RNA was calculated according to the 2− ΔΔCT method. The qPCR primers are: H− CXCL1− F: ACAACAATTACGCGCTGCGT; H− CXCL1− R: GTTTCTTAACTATGGGGGATGC; H− ACTIN− F: CACTCTTCCAGCCTTCCTTC; H− ACTIN− R: GTACAGGTCTTTGCGGATGT.
Immunohistochemistry/immunofluorescence staining
Glioma tissues and paraffin sections were provided by the Department of Neurosurgery, Shanghai Changhai Hospital, Naval Medical University. Tumor tissues were sent to Servicebio for staining analysis.
Statistical analysis
All statistical analyses were conducted using R software (v4.3.1). Associations between continuous variables were assessed via Spearman’s correlation test, while differences between two groups were evaluated with the Wilcoxon test. Survival differences were analyzed using Kaplan-Meier curves with the log-rank test. A p-value of less than 0.05 was deemed statistically significant.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Louis DN Perry A Reifenberger Gvon Deimling A Figarella-Branger D Cavenee WK. The 2016 World Health Organization classification of tumors of the central nervous system: a summary. Acta Neuropathol (Berl). (2016) 131:803–20. doi: 10.1007/s 00401-016-1545-1 27157931 · doi ↗ · pubmed ↗
- 2Verhaak RGW Hoadley KA Purdom E Wang V Qi Y Wilkerson MD. An integrated genomic analysis identifies clinically relevant subtypes of glioblastoma characterized by abnormalities in PDGFRA, IDH 1, EGFR and NF 1. Cancer Cell. (2010) 17:98. doi: 10.1016/j.ccr.2009.12.020 20129251 PMC 2818769 · doi ↗ · pubmed ↗
- 3Brennan CW Verhaak RGW Mc Kenna A Campos B Noushmehr H Salama SR. The somatic genomic landscape of glioblastoma. Cell. (2013) 155:462–77. doi: 10.1016/j.cell.2013.09.034 PMC 391050024120142 · doi ↗ · pubmed ↗
- 4Omuro A De Angelis LM. Glioblastoma and other Malignant gliomas: a clinical review. JAMA. (2013) 310(17):1842–50. doi: 10.1001/jama.2013.280319 24193082 · doi ↗ · pubmed ↗
- 5Stupp R Hegi ME Mason W Pvan den Bent MJ Taphoorn MJB Janzer RC. Effects of radiotherapy with concomitant and adjuvant temozolomide versus radiotherapy alone on survival in glioblastoma in a randomised phase III study: 5-year analysis of the EORTC-NCIC trial. Lancet Oncol. (2009) 10:459–66. doi: 10.1016/S 1470-2045(09)70025-7 19269895 · doi ↗ · pubmed ↗
- 6Weller M Wen PY Chang SM Dirven L Lim M Monje M. Glioma. Nat Rev Dis Primer. (2024) 10:33. doi: 10.1038/s 41572-024-00516-y 38724526 · doi ↗ · pubmed ↗
- 7Wolchok JD Chiarion-Sileni V Gonzalez R Rutkowski P Grob J-J Cowey CL. Overall survival with combined nivolumab and ipilimumab in advanced melanoma. N Engl J Med. (2017) 377:1345–56. doi: 10.1056/NEJ Moa 1709684 PMC 570677828889792 · doi ↗ · pubmed ↗
- 8Gandhi L Rodríguez-Abreu D Gadgeel S Esteban E Felip E De Angelis F. Pembrolizumab plus chemotherapy in metastatic non-small-cell lung cancer. N Engl J Med. (2018) 378(22):2078–92. doi: 10.1056/NEJ Moa 1801005 29658856 · doi ↗ · pubmed ↗