Exploring the Landscape of Histone Modifications in Glioblastoma Multiforme
Based on the established temozolomide-resistant cell lines, U251-TR and U118-TR, we performed transcriptomic sequencing and differential expression analysis. The differential analysis volcano plots for these cell lines are displayed in Figure S1A and S1B. A total of 869 and 5823 differentially expressed genes were identified in the U118 and U251 cell lines, respectively. The Venn diagram (Figure S1C) highlights 73 genes commonly differentially expressed across both cell lines. Functional enrichment analysis of these common genes revealed that they are predominantly involved in chromatin regulation and post-translational histone modifications (Figure S1D).
Subsequently, using single-sample enrichment analysis, we quantitatively assessed histone post-translational modification-related pathways (sourced from the GSEA database) across samples from TCGA-GBM (IDH1-wildtype), GTEx normal brain tissue, and CGGA-GBM (IDH-wildtype, including datasets 325 and 695). We observed that compared to normal brain tissue, GBM tissue samples exhibit active histone post-translational modifications (Fig. 1A). Further analysis of primary and recurrent GBM samples from the CGGA dataset identified significant differences in pathways related to H3K4 monomethylation, H3K9 demethylation, H3K27 methylation, H3K36 demethylation and methylation, and H4K5 and H4K12 acetylation (Fig. 1B, p < 0.05). Moreover, varying levels of activity in pathways such as H3K4 dimethylation, H3K4 monomethylation, H3K9 acetylation and demethylation, H3K14 acetylation, H3K27 methylation, H3K36 methylation, H4K5 acetylation, H4K12 acetylation, and H4K20 methylation were observed across different GBM subtypes (Fig. 1C, D, p < 0.05). Spearman correlation analysis revealed a significant positive correlation between various histone modification pathway activities and GBM tumor stemness, especially in pathways like H4K5 acetylation, H4K20 methylation, H4K16 acetylation, H4K12 acetylation, H3K9 tri-methylation, dimethylation, and acetylation, H3K36 methylation, H3K27 methylation, and H3K14 acetylation when correlated with GBM stemness score (mRNAsi) (Fig. 1E). Lastly, we examined inter-correlations among various histone modification pathways. Our results indicate that different types of modifications generally show significant positive correlations, such as H3K9 acetylation with H3K14 acetylation, and H4K5 acetylation with H4K12 acetylation, H4K16 acetylation, H3K4 dimethylation, H3K4 monomethylation, H3K36 methylation, and H4K20 methylation (Fig. 1F). This gives a comprehensive view of the histone modification landscape in GBM tumor samples.
Correlation of Histone PTMs with Temozolomide Resistance in GBM Cell Lines
Subsequently, to investigate the association between various histone PTMs and TMZ resistance, we curated transcriptomic data and corresponding TMZ IC50 values from the GBM cell lines within the CCLE database. We recalculated the ssGSEA scores for various histone PTMs. Based on the ascending order of IC50 values, we chronologically arranged different GBM cell lines (Fig. 2A). Using the median IC50 value of 7.08µM as a threshold, the cell lines were classified into IC50-high (TMZ-resistant) and IC50-low (TMZ-sensitive) groups. Spearman correlation analysis identified several pathways showing a robust positive correlation with IC50, namely H3K36 demethylation, H3K4 trimethylation, H3K4 monomethylation, and H3K9 acetylation, as visualized in a heatmap (Fig. 2B).
To further understand the predictive power of these histone PTMs for IC50 values, we implemented machine learning algorithms. Comparing the residuals of three different machine learning models (Fig. 2C, D), the Generalized Linear Model (GLM) consistently demonstrated the smallest residuals, outperforming both the Random Forest and SVM algorithms. This superiority was further confirmed by the ROC curve, with the GLM model attaining an AUC value of 0.833, surpassing the other two models (Fig. 2E). Moreover, the top ten pathway variables for each algorithm are displayed, with the three foremost in GLM being h3k4me, H3K9ac, and h3k4me3, whereas in the Random Forest model, they were H3k9ac, h3k14ac, and h3k9 demethylation (Fig. 2F). Consolidating these metrics, GLM exhibited superior predictive capability. We thus presented the feature importance values computed by the GLM algorithm, where h3k4me and H3K9ac scored 0.7 and 0.68, respectively (Fig. 2G).
To explore potential downstream targets regulated by H3k4me1 and h3k9ac, differentially expressed genes between IC50 subgroups (Fig. 2H) were subjected to WGCNA. Optimal soft power value was determined to be 6 based on scale independence and mean connectivity (Figure S2A, B), leading to gene clustering (Figure S2C). Ultimately, genes were clustered into five modules: turquoise, yellow, blue, brown, and grey. A heatmap demonstrating the correlation between these module vectors and h3k4me1, h3k9ac, and IC50 is presented (Fig. 2I). Modules with stronger correlation to IC50 (MEturquoise, MEyellow, and MEblue) predominantly showed significant associations with h3k9ac (Figure S2D-I), while such a trend was absent for h3k4me1. Further functional annotation via KEGG enrichment analysis revealed genes in the turquoise module primarily align with "Glycerophospholipid metabolism", those in the yellow module align with the TGF-beta signaling pathway, and those in the blue module were implicated in the mTOR signaling pathway (Fig. 2J). To gain a deeper understanding of other potential pathways influenced by H3k9ac, we conducted GSEA functional enrichment analysis, which prominently identified "WNT beta-catenin signaling" and "DNA repair" in the H3K9ac high-activity subgroup (Fig. 2K, L). Additional significant KEGG pathways were depicted in Figure S2J-O. Our comprehensive analysis underscores the profound influence of histone PTMs, especially h3k9ac, on TMZ resistance in GBM, prompting us to further explore the regulatory mechanisms underlying the expression of MGMT, a pivotal enzyme in DNA repair.
Correlation Between Histone PTMs and MGMT Expression in GBM
Considering that elevated MGMT expression is a pivotal mechanism for TMZ resistance in GBM cells and the methylation status of the MGMT promoter plays a significant role in regulating its expression31, we analyzed the correlation between various histone post-translational modifications (PTMs) and MGMT expression regulation. Figure 3A displays the expression levels of MGMT alongside the histone PTMs corresponding to the methylation status of the MGMT promoter. We confirmed that methylation of the MGMT promoter notably downregulates MGMT mRNA expression (Fig. 3B). In subgroups with varying MGMT expression levels, we observed a significant upregulation of H3K4me2, h3k4me1, h3k9ac, h3k9me3, h3k36me, and h4k20me in the low MGMT expression group (Fig. 3C). Examining the correlation between the activity of histone PTMs and the methylation status of the MGMT promoter (Fig. 3D), we noticed a significant upregulation of h3k9ac in the subgroup with MGMT promoter methylation (p = 0.0008). We suspected that the notable downregulation of H3K9ac in the low MGMT expression group might be associated with the decreased MGMT expression caused by the promoter methylation. Hence, we conducted a further analysis in the TCGA-GBM database. We categorized all samples based on the methylation status of the MGMT promoter and integrated the data into scatter plots (Fig. 3E-F). Overall, Spearman's analysis indicated a significant negative correlation between MGMT expression and h3k9ac activity (R=-0.38, p < 0.05) (Fig. 3G). In the subgroup with MGMT promoter methylation (Fig. 3H), no significant correlation was observed (R=-0.24, p = 0.083). In contrast, in the non-methylated MGMT promoter subgroup (Fig. 3I), there was a significant positive correlation between MGMT expression and H3K9ac activity (R = 0.27, p < 0.05). We validated these observations in the CGGA-GBM dataset (Fig. 3J-N), and consistent findings were obtained, with a significant positive correlation between MGMT and H3K9ac (R = 0.29, p < 0.05) only observed in the non-methylated MGMT promoter subgroup.
The observed outcomes might be attributed to the intricate interplay between histone modifications and DNA methylation32. It is conceivable that the methylation of the MGMT promoter may impede the recruitment or binding of certain transcription factors or coactivators that interact with H3K9ac. Alternatively, the upregulation of H3K9ac in the absence of MGMT promoter methylation might facilitate a more open chromatin structure, promoting MGMT expression33.
Interplay of MGMT Promoter Methylation and H3K9ac Levels in Determining GBM Outcomes
Considering the potential intricate regulation of MGMT expression by both the methylation of the MGMT promoter and H3K9ac levels, we sought to further explore the prognostic relevance of MGMT and H3K9ac levels in GBM patients treated with TMZ. Initially, among these patients (n = 110), those with elevated MGMT expression demonstrated a notably worse overall survival and progression-free survival (Fig. 4A, B). A similar trend was observed with H3K9ac, where patients with higher H3K9ac levels had significantly poorer outcomes (Fig. 4C, D). Moreover, recognizing that a major cause of GBM recurrence is the insensitivity of patients to adjuvant treatments post-surgery7, including TMZ, we reanalyzed and compared H3K9ac levels in primary and recurrent GBM patients treated with TMZ across different datasets from CGGA (CGGA_325 and CGGA_695). We observed a significant upregulation of H3K9ac in the 325 dataset (p < 0.05) (Fig. 4E, F). Crucially, the potential influence of MGMT promoter methylation on the effects of H3K9ac on GBM prognosis cannot be ignored. Thus, a combined survival analysis curve accounting for both MGMT promoter methylation status and H3K9ac activity was generated (Fig. 4G). Results from TCGA-GBM indicated that under MGMT promoter methylation, patients with higher H3K9ac levels had significantly worse prognosis. In contrast, in the absence of MGMT promoter methylation, H3K9ac levels did not significantly impact overall survival. Viewing this from another angle, under low H3K9ac levels, the methylation status of the MGMT promoter significantly affected patient prognosis, whereas under high H3K9ac levels, the status of the MGMT promoter did not. This phenomenon was observed in the CGGA-GBM dataset (Fig. 4H). The observed patterns suggest a multifaceted regulatory mechanism wherein both MGMT promoter methylation and H3K9ac influence the prognosis of GBM patients.
Reciprocal Regulation of MGMT and H3K9ac in TMZ Resistance
To further elucidate the interplay between MGMT and H3K9ac in the development of TMZ resistance in GBM, we established two TMZ-resistant cell lines, U118-TR and U251-TR. Figure 5A outlines the workflow used to develop these resistant strains. Subsequently, we confirmed the resistance of these cell lines to TMZ treatment using MTT assays (Fig. 5B, C). The half-maximal inhibitory concentration (IC50) of TMZ for U118, U118-TR, U251, and U251-TR were found to be 174.3, 215.3, 142.1, and 210.0µg/ml respectively. Notably, the IC50 values of the resistant cell lines were significantly higher than the control, which is consistent with our transcriptome sequencing samples (Figure S1).
MGMT mRNA expression levels were evaluated using qPCR across the cell lines, revealing a notable upregulation in the resistant strains (Fig. 5F). At the protein level, MGMT expression was also heightened in the resistant lines (Fig. 5G). Additionally, nuclear protein extraction demonstrated an augmented presence of the H3K9ac modification in the TMZ-resistant strains (Fig. 5G). Immunofluorescence analysis further pinpointed the localization of H3K9ac in the cell nucleus, with a pronounced increase in the resistant strains (Fig. 5H).
Importantly, leveraging the chip-seq and corresponding RNA-seq data from the GSE11381629 dataset in the GEO database (Figure S3A), we visualized H3K9ac and MGMT data on IGV (Fig. 5I, Figure S3B). There was a significant elevation of H3K9ac at the MGMT DNA locus in the TMZ-resistant strains, correlating with an increase in MGMT RNA expression. Comparing the most enhanced H3K9ac regions with annotated regulatory elements, we observed a concentration of H3K9ac around the promoter-proximal transcription factor binding sites and three subsequent enhancer segments. It is postulated that H3K9ac may facilitate chromatin accessibility, thereby modulating transcription factor binding and subsequently regulating MGMT transcription. Moreover, we explored histone PTMs in the resistant strains on DNA segments of several other TMZ resistance-associated genes, such as ATRX, MPG, and PARP1 (Figure S3B), yet no conspicuous upregulation of H3K9ac was discerned.
Our findings emphasize the pivotal role of H3K9ac in augmenting MGMT expression in TMZ-resistant GBM cells, hinting at a complex regulatory mechanism that warrants deeper exploration.
SP1 and H3K9ac's Concerted Regulation of MGMT in TMZ-Resistant GBM Cells
To elucidate the transcription factors (TFs) involved in the regulation of MGMT transcription and ascertain if they operate in conjunction with H3K9ac, we initially leveraged the TRRUST online database (https://www.grnpedia.org/trrust/) to predict transcription factors targeting MGMT. Visualization was performed using Cytoscape (Fig. 6A), identifying several transcription factors including EP300, SP1, DNMT1, RELA, NFKB1, MECP2, MBD1, JUN, HDAC1, and TP53. Subsequent analysis explored potential pathways these TFs might collaboratively partake in (Fig. 6B, C). To understand the correlation between the RNA expression levels of these transcription factors, the H3K9ac pathway, MGMT expression levels, and IC50 values, heatmap results revealed only NFKB1 and SP1 demonstrated significant positive correlations with H3K9ac, MGMT, and IC50. Validation using the CGGA-GBM dataset further supported only SP1 showing a significant positive correlation with both H3K9ac and MGMT (Fig. 6E). Consequently, our attention shifted to the regulation of MGMT by SP1. Employing the online tool, JASPAR, we scanned the TF binding site on MGMT and its surrounding sequence (where H3K9ac notably upregulates in the drug-resistant group). Successful prediction highlighted potential SP1 binding sites in the region. The predicted motifs and associated sequences are displayed in Fig. 6F. The first potential binding site was located at 5‘-AATCCTCATGCCAACC-3' with a binding index of 5.965 and a relative score of 0.818, positioned within the TF binding site. The second potential binding site was at 5'-TTTCTTCCTCT-3', with a binding index of 2.696 and a relative score of 0.815, situated near the TF binding site. Subsequent qPCR and mass spectrometry analysis confirmed an upregulation of SP1 at both the RNA and protein levels in the constructed TMZ-resistant strains (Fig. 6G, H). Notably, primers designed based on the MGMT TF binding site sequence for chip-qPCR also indicated significant enrichment of both H3K9ac and SP1 in this region in the drug-resistant cell strains (Fig. 6I). These results underscore SP1's central role in regulating MGMT transcription, potentially modulated by H3K9ac.
Single-cell Transcriptomic Profiling Uncovers H3K9ac, SP1 and MGMT Dynamics in GBM Subtypes
To elucidate the expression patterns and correlations between H3K9ac and MGMT across distinct cellular subpopulations, we conducted single-cell transcriptomic analysis in glioblastoma multiforme (GBM). Utilizing the TISCH2 single-cell transcriptome database, the GSE131298_10X dataset comprising 9 GBM samples with 13,553 identified cells was visualized. These cells were categorized into 26 clusters (Fig. 7A), which were further annotated as malignant tumor cells, immune cells, and other cell types (Fig. 7B). Predominantly, malignant tumor cells were sub-typed into neural-progenitor-like (NPC-like), oligodendrocyte-progenitor-like (OPC-like), astrocyte-like (AC-like), and mesenchymal-like (MES-like) categories (Fig. 7C), with the proportional distribution of these subtypes in each sample depicted in Figure S4A. Mapping the expression of H3K9ac-related genes allowed us to represent H3K9ac pathway activity on the dimensionally reduced cellular landscape (Fig. 7D). Differential analysis with violin plots confirmed the upregulation of the H3K9ac pathway in malignant tumors (Figure S4B). Further, MES-like cells exhibited marginally higher average levels of H3K9ac compared to other subtypes (Figure S4C). Given their known aggressive nature and resistance to chemotherapeutic interventions34, the mesenchymal subtype's association with the H3K9ac pathway was intriguing. Through GSEA, we compared the activity of the DNA-repair pathway and the WNT-β catenin pathway (a previously identified potential downstream pathway, Fig. 2K, L) across different cell subgroups (Fig. 7E, F), indicating a stronger alignment of H3K9ac with the DNA-repair pathway.
To delve into the expression dynamics of H3K9ac, SP1, and MGMT during tumor recurrence, paired primary and recurrent GBM samples from the GSE174554 dataset were studied. These patients underwent standard TMZ therapy as recommended by Stupp35 but exhibited rapid relapse. Single-cell transcriptomic data underwent rigorous quality control and preprocessing (Figure S4D). After selecting optimal PC numbers (7, Figure S4E, F) and resolution (0.1, Figure S4G), eight clusters were identified (Figure S4H, I). Dimensional reduction visualized different primary and paired samples (Fig. 7G, H). With references from the singleR package (Figure S4J) and GBM subtype-associated markers (Figure S4K), cellular annotations were established (Fig. 7I). When comparing H3K9ac, SP1, and MGMT levels between primary and recurrent samples, all showed increased expression in the recurrent subgroups (Fig. 7J). Pseudotemporal analysis using the Monocle R package (Fig. 7K) depicted the trajectories of primary and recurrent cells over time, with primary samples primarily in early phases (Fig. 7L). Earlier cell populations were characterized as Glioblastoma Stem-like Cells (Figure S4L). As tumor recurrence progressed and using GAPDH as a reference, the stemness of the tumor decreased, while MGMT, SP1, and H3K9ac levels trended upwards (Fig. 7M). Some cellular subpopulations persisted as mesenchymal GBM, followed by differentiation into astrocytic and classic types, subsequently differentiating into neuronal GBM (Figure S4L, M, Fig. 7N).