Defining a Correlative Transcriptional Signature Associated with Bulk Histone H3 Acetylation Levels in Adult Glioblastomas

Glioblastoma (GB) is the most prevalent primary brain cancer and the most aggressive form of glioma because of its poor prognosis and high recurrence. To confirm the importance of epigenetics in glioma, we explored The Cancer Gene Atlas (TCGA) database and we found that several histone/DNA modifications and chromatin remodeling factors were affected at transcriptional and genetic levels in GB compared to lower-grade gliomas. We associated these alterations in our own cohort of study with a significant reduction in the bulk levels of acetylated lysines 9 and 14 of histone H3 in high-grade compared to low-grade tumors. Within GB, we performed an RNA-seq analysis between samples exhibiting the lowest and highest levels of acetylated H3 in the cohort; these results are in general concordance with the transcriptional changes obtained after histone deacetylase (HDAC) inhibition of GB-derived cultures that affected relevant genes in glioma biology and treatment (e.g., A2ML1, CD83, SLC17A7, TNFSF18). Overall, we identified a transcriptional signature linked to histone acetylation that was potentially associated with good prognosis, i.e., high overall survival and low rate of somatic mutations in epigenetically related genes in GB. Our study identifies lysine acetylation as a key defective histone modification in adult high-grade glioma, and offers novel insights regarding the use of HDAC inhibitors in therapy.


Introduction
Gliomas constitute a heterogenous group of primary brain tumors that include the especially aggressive glioblastomas (GB) [1]. Current treatments are primarily based on the application of the Stupp protocol, which consists of surgical resection of the tumor mass, radiotherapy, and chemotherapy with temozolomide (TMZ) [2]; however, they are still insufficient to avoid recurrence in nearly all cases. Epigenetics dysregulation can be relevant in the formation and maintenance of GB because glioma stem cells, which are thought to be responsible of GB recurrence, are maintained in an undifferentiated and self-renewal state due to a permanent epigenetic block [3]. Thus, the epigenetic landscape of the tumoral mass can have a significant impact on the patient outcome. For example, the silencing of the O-6-methylguanine-DNA methyltransferase (MGMT) gene by hypermethylation of its promoter improves the efficacy of alkylant agents such as TMZ [4]. Moreover, the most important molecular alteration with diagnosis utility to date, the mutation of arginine 132 in the isocitrate dehydrogenase (IDH1) gene, is sufficient to remodel the DNA methylome [5], a phenomenon that may explain the better outcomes in mutated gliomas compared to wild type tumors. Samples were processed following the procedures described in [23] for Western blotting and retrotranscription-quantitative PCR (RT-qPCR) assays and RNA-seq analysis. For Western blotting normalization, we first normalized the signal intensities of the band with the signal intensity of the same external protein extract loaded in all gels; we next normalized the values of each histone modification with the corresponding value of total histone H3 per sample, as we were interested on the fraction of the total protein with such covalent modification independently on potential variations of total histone H3 across gliomas. Samples with low signal intensity for total histone H3 were discarded. Western blotting assays were performed using the same cohort previously reported in [23], adding 7 GBs not included at the time in these assays. For RT-qPCR, we used TBP as the housekeeping gene. For Western blotting, the following antibodies were used: H3K9,14ac (06-599), H3K4me3 (07-473), H3K27me3 (07-449) (1:1000, EMD Millipore, Darmstadt, Germany), total H3 (1:8000, ab1791, Abcam, Cambridge, UK), and HRP-conjugated secondary antibodies (1:7500, A0545 and A4416, Sigma-Aldrich, Darmstadt, Germany). The sequences of primer pairs used in the RT-qPCR assays are listed as follows: A2ML1, 5 -CATTGTTGGCCCAGCTTACC-3 and 5 -ATGTGCGCTGGAAATTCTCAG-3 , CD83, 5 -GAGGGTGGTGAAGAGAGGATG-3 and 5 -CTCTTCTTTACGCTGTGCAGG-3 , CDYL2, 5 -CTCAGATACA-GTGTCCGCCAG-3 and 5 -TCGCCGGACTTCTTTCATGAT-3 , GABRA4, 5 -GATGGTCATGCATGCCCTTTG-3 and 5 -ACTTGATACGGTTTGCCCAATC-3 , GFAP, 5 -AGAGGAAGATTGAGTCGCTGG-3 and 5 -TGTCAGGTCTGC-AAACTTGGA-3 , HPCAL4, 5 -GAAGCTCAACTGGGCCTTTG-3 and 5 -GGTCGTCCTTATCCTGGTCC-3 , PROM1 (CD133) 5 -ACTCAGCGTCTTCCTATTCAG-3 and 5 -AAAATCACGATGAGGGTCAGC-3 , SLC17A7, 5 -GCTTCGGGATCTTCTGGTACC-3 and 5 -CAGAAGTTGG-CCACGATGATG-3 , SLITRK4, 5 -AATCCTGACTGTGGCTCCAT-3 and 5 -ATCTTTCTCCTTGTCGGCCA-3 , TNFSF18, 5 -TCTTTGCTCCTTCAGTTGGC-3 and 5 -ATACAGCCGCACCTCAAAAG-3 . The sequences of TBP primers are described in [23].

DNA Methylome Analysis
In total, 8 GBs were randomly selected for this analysis: 6 females and 2 males with a median of age of 62.5 years (IQR = 11.25). Bisulfite-treated DNA from these samples was labelled and hybridized into Infinium MethylationEPIC beadchips (Illumina, San Diego, CA, USA) according to the manufacturer's recommendations at the Unidad de Genómica of Genyo (Granada, Spain).
The "minfi" software was used for data extraction, quality assessment, background correction with dye bias normalization (using the method of normal-exponential out-ofband or Noob, implemented in the package), and the filtering of probe sets located in SNPs and sexual chromosomes [24]; β-values (methylated Cy5/unmethylated Cy3 signals) were normalized using the BMIQ method implemented in the "RnBeads" R-based package [25], and the CpG sites were mapped into the GRCh38/hg38 genome assembly at the Unidad de Bioinformática of Genyo. Differential methylation was determined using the "limma" package [26]. We only considered the CpG sites that were mapped into the genes interrogated in the RNA-seq analysis, resulting in 13,354 genes containing both transcriptomics and epigenomics data in our cohort of study.

Primary Glioblastoma Cultures
For culturing, a glioblastoma sample from a surgical resection was dissociated enzymatically and the cells were maintained as neurospheres in DMEM/F-12 medium supplemented with B27 and N2 (Gibco, Thermo Fisher, Madrid, Spain) plus epidermal growth factor (EGF) and basic fibroblast growth factor (bFGF) (10 ng/mL, PeproTech, London, UK) [27]. Seven days later, neurospheres were attached to a laminin-coated flask using serum-containing medium [28]. These cultures were established in conditions that favored the proliferation and maintenance of glioma stem cells [28], while minimizing the heterogeneity of resected tissues and enabling their pharmacological manipulation. Subsequently, cells were seeded in 12-well plates (5 × 10 4 cells/well) and treated with 2 µM TSA (Sigma-Aldrich, Darmstadt, Germany) for 24 h. After treatment, cell viability was assessed using the XTT cell proliferation assay according to manufacturer's instructions (Canvax Biotex, Valladolid, Spain) and cells were processed for downstream experiments using the same procedures as tissue resections [23]. The only modification concerned RNA-seq analysis in which DNA libraries were obtained using the TruSeq Stranded mRNA kit (Illumina, San Diego, CA, USA) and sequenced in a 75 bp paired-end configuration (Unidad de Genómica, Cabimer, Sevilla, Spain).

Additional Statistical Analysis and Bioinformatics
Mann-Whitney U tests, Student's t-tests, Pearson's product moment correlations, χ 2 test, and heatmap plotting were performed using the native R environment. Principal component analysis (PCA) was conducted with the "rgl" package (http://cran.r-project. org/package=rgl) (accessed on 6 June 2022), and survival analysis using the "survminer" and "survival" packages to calculate both the Kaplan-Meier curves and log-rank p-values (accessed on 20 October 2022) in the R environment (R version 4.0.5). In this analysis, we divided the samples into three parts after trial and error, and compared the days to death associated with the third of data containing the highest expression values and the third of data containing the lowest expression values for each dataset.
We used the genetics and transcriptomics information from The Cancer Gene Atlas (TCGA) consortium as deposited in the Genomic Data Commons (GDC) website (https: //portalgdc.cancer.gov, accessed on 6 June 2022), and from the REMBRANDT cohort as deposited in GeneBank with the accession number GSE108474. Differential expression between GB and lower-grade glioma (LGG) datasets was previously calculated [23,29].
Epigenetic genes were obtained from the EpiFactors database [30], which included enzymes and associated cofactors related to epigenetic regulation. Of the 815 proteins contained in the database, we only considered 590 related to histone modification (acetylation, methylation, phosphorylation, and ubiquitylation), DNA modification, chromatin remodeling, and readers, either containing enzymatic activity or acting as cofactors, after excluding histones, protamines, chaperones, epitranscriptomics and unclassified factors.

Histone Acetylation Levels Is Defective in Glioblastoma Compared to Lower-Grade Gliomas
In an attempt to explain glioma malignancy, we compared the gene expression profiles between GB and LGG using the transcriptomics datasets contained in the GDC portal (hereafter referred to as TCGA datasets), under the assumption that the differential expression between these two major types of glioma might identify the key factors that enhance the characteristic aggressivity of GB. This analysis revealed an extensive differential expression:~23% of the whole transcriptome (fold change > 2, adjusted p-value < 0.05), which affected~11% of genes related to DNA and histone modifiers and readers, chromatinremodeling genes, and cofactors according to the curation of EpiFactors [30] (hereafter referred to as epigenetic genes). Of these, genes modulating histone acetylation were more consistently affected because they exhibited the largest and significant differential expression, according to the absolute averages of fold change and adjusted p-values for the common genes also differentially expressed in the REMBRANDT cohort (adjusted p-value < 0.05) ( Figures 1A and S1). In addition to the gene expression variations, we also observed a slightly higher rate of somatic mutations in these histone acetylation regulators in the TCGA gliomas compared to other epigenetic categories; KAT6B was the only gene containing the same mutation (T1203Rfs*21) at a single position in more than one case (four cases), but no mutated genes were found in "histone phosphorylation", "reader", or "chromatin remodeling" categories ( Figure 1B). Overall, this analysis indicates that enzymes and cofactors regulating histone modifications, especially histone acetylation, were affected in gliomas and might be linked to relevant clinical outcomes. Figure 1. Histone H3 deacetylation in glioblastomas: (a) Genes with differential expression between 211 the datasets of "GBM" and "LGG" projects of TCGA that were also confirmed in the REMBRANDT 212 cohort. Individual values of fold change and adjusted p-values are shown in Figure S1. FDR, false 213 discovery rate. (b) For the number of somatic mutations, we only consider nonsynonymous substi-214 tutions in protein coding regions in the differentially expressed genes of (a). (c) Upper panel, repre-215 sentative Western blots for total and covalently modified histones in tumor resections where each 216 lane represents a patient's sample; the number denotes the histological grade. Lower panel, quanti-217 fication of the blots after normalization by the total H3 signal. n = 13 for grade II, n = 7 for grade III, 218 n = 57 for grade IV (GB). **, p-value < 0.005, Mann-Whitney U test related to grade II. Raw blots are 219 shown in Figure S2 with the location of the image crops. (d) Pearson correlation coefficients (r) for 220 H3K9,14ac levels and other histone modifications (this study) and histone H3 variants [23]. The plot 221 shows the correlation between H3K9,14ac and H3.3, the only one to be significant. 222 Figure 1. Histone H3 deacetylation in glioblastomas: (A) Genes with differential expression between the datasets of "GBM" and "LGG" projects of TCGA that were also confirmed in the REMBRANDT cohort. Individual values of fold change and adjusted p-values are shown in Figure S1. FDR, false discovery rate. (B) For the number of somatic mutations, we only consider nonsynonymous substitutions in protein coding regions in the differentially expressed genes of (A). (C) Upper panel, representative Western blots for total and covalently modified histones in tumor resections where each lane represents a patient's sample; the number denotes the histological grade. Lower panel, quantification of the blots after normalization by the total H3 signal. n = 13 for grade II, n = 7 for grade III, n = 57 for grade IV (GB). **, p-value < 0.005, Mann-Whitney U test related to grade II. Raw blots are shown in Figure S2 with the location of the image crops. (D) Pearson correlation coefficients (r) for H3K9,14ac levels and other histone modifications (this study) and histone H3 variants [23]. The plot shows the correlation between H3K9,14ac and H3.3, the only one to be significant.
To determine whether these alterations correlated with effective changes in histone modifications, we performed Western blotting assays in two collections of glioma surgical resections that are described in our previous work [23]: after scanning the bulk levels of well-known covalent modifications of histone H3 in tumor resections, only the acetylated lysines 9 and 14 (H3K9,14ac or H3ac) were significantly different between GB and grade II gliomas after normalizing by total histone H3 levels; grade III represented an intermediate state, although caution should be taken in this interpretation because of the low number of samples and the sex bias. For simplicity, Figure 1C shows the combined results for both cohorts. However, the bulk levels for trimethylation of either H3K4 or H3K27 were not significantly different across grades ( Figure 1C) and did not correlate with the levels of acetylated histone H3 ( Figure 1D). In contrast, the bulk levels of the replication-independent variant H3.3, also reduced in adult high-grade gliomas [23], were significantly correlated with H3K9,14ac signal intensities ( Figure 1D), suggesting that H3.3 was acetylated at these lysines.

A Transcriptional Signature Can Be Correlated with Differential Bulk Levels of Histone Acetylation in Glioblastoma
Our Western blotting results reveal general histone hypoacetylation in high-grade tumors. To gain further insights regarding the biological consequences of such alteration in GB, we selected the three GBs showing the lowest and highest levels of H3K9,14ac for which RNA was available in our cohorts, namely, H3ac low and H3ac high (Figure 2A). No differences were observed for bulk levels of trimethylation of K4 or K27 between the selected samples ( Figure 2A). In the PCA analysis, these two types of GB were similarly separated from low-grade gliomas despite the comparable bulk histone acetylation levels between H3 high GB and low-grade gliomas. This result suggests that different levels of bulk histone H3 acetylation did not have a profound impact on the GB transcriptomes ( Figure 2B). After performing differential expression analysis (adjusted p-value < 0.05), we identified 249 differentially expressed genes (DEGs) between H3ac high and H3ac low resections, in which 218 genes (~88%) and 31 genes were enriched in H3ac high and H3ac low tumors, respectively (Table S1). This 249-gene signature was specifically associated with H3K9,14ac levels, as this mark showed better correlation coefficients with the median of gene expression than the other histone modifications analyzed in this work ( Figure 2C). In the genes enriched in H3ac low GB, we found significant terms for regulation of neurogenesis and neuronal differentiation (e.g., LHX2, FZD8, GPC2, INSM1), whereas mature neuronal markers were observed in H3ac high GB that encoded for synaptic proteins (e.g., SYT1, SYN2) and neurotransmitter receptors (e.g., GABA receptor subunits, GRM1; Figure 2D), suggesting that bulk levels of H3K9,14ac might indicate divergent stages of neuronal differentiation and different abundance of healthy neurons within resections.
To understand why some particular genes showed an expression correlated with bulk tumoral levels of H3K9,14ac, we examined whether their gene expression differed from the general GB transcriptome. In addition, we also used information regarding DNA methylation to infer the general epigenetic status of these genes: we randomly selected some samples from our cohort to hybridize bisulfate-treated DNA into Infinium MethylationEPIC beadchips. We did not explore the binding of acetylated histone to DNA, as previous work reported that differential gene expression and differential histone acetylation occupancy across gliomas converged in a low number of common genes [31]. Furthermore, indirect effects on gene expression cannot be disregarded in our histone acetylation-dependent transcriptional signature; therefore, we reason that DNA methylation might be suitable to recapitulate the epigenetic status of the loci of interest. To this aim, we averaged the expression/methylation values of each gene that were common to the RNA-seq and DNA methylation arrays, and we compared the median values of the enriched genes in H3ac high GB with those of the whole population; the number of enriched genes in H3ac low GB were too low to extract meaningful conclusions. The genes associated with high levels of histone acetylation were mostly expressed below the median of gene expression in GB ( Figure 2E), and, according to this low expression, were more methylated at their CpG islands (CGIs) ( Figure 2F); this is not unsurprising because CGIs usually colocalize within promoters, and have been more tightly linked with the regulation of transcription [32]. Despite their increased methylation, these CGIs still exhibited the expected CpG hypomethylation compared to other genomic features such as shores, shelves, and open seas ( Figure 2F). However, this analysis was restricted to the significance cutoff used to define DEGs, probably missing general patterns that might comprise more genes, i.e., although not individually significant, they might be part of changes affecting gene subpopulations. Therefore, we investigated the median values of both gene expression and CGI methylation in windows of 500 genes across the whole GB transcriptomes once arranged according to significance and direction of change, as retrieved in the pairwise comparison between H3ac high and H3ac low GBs. In this manner, we revealed a general extended association between gene responsiveness to elevated histone acetylation at lysines 9 and 14 and low basal gene expression and CGI hypermethylation ( Figure 2G), suggesting that well-expressed genes did not need to increase their histone acetylation state to increase their rate of expression [33].
To confirm that the differential levels of H3ac between GBs were linked to the reported transcriptional variations in the RNA-seq analysis, we manipulated H3ac levels by pharmacological means. To this end, we inhibited HDAC activity in GB-derived cultures. Treatment with 2 µM of the HDACi TSA for 24 h resulted in a striking induction of histone hyperacetylation ( Figures 2H and S3A) without significant affectation in the viability of cells: untreated, n = 6, fold change = 1.00 ± 0.28; TSA-treated, n = 6, fold change = 1.23 ± 0.21; SDS-treated, n = 2, fold change = 0.00 ± 0.24 in XTT assays. This timing for RNA sampling enabled the identification of the transcriptional changes associated with the TSA treatment prior to any potential cell loss that might introduce bias in the results. In these conditions, the RNA-seq between treated and control cells revealed an extensive transcriptional rearrangement consisting of 1995 and 1261 upregulated and downregulated genes, respectively (Table S2); we validated the upregulation of selected genes in independent treated cultures ( Figure S3B). Among the DEGs, we noticed the enrichment of downregulated genes related to cell division and mitosis, and the enrichment of upregulated genes related to synaptic transmission ( Figure S3C); the latter result is reminiscent of the GO enrichment documented in the genes that were highly expressed in H3ac high tumors. We also observed an important fraction of epigenetically related genes among DEGs. In contrast, TSA upregulation only affected the expression of~2.9% epigenetic genes; this percentage was increased to~11.0% within downregulated genes with TSA, probably as part of a homeostatic response against the pharmacological treatment. These results suggest that HDAC inhibition might have relevant consequences in the reconfiguration of the GB epigenetic landscape. Importantly, these DEGs behaved similar to genes in H3ac high and H3ac low tumors: we plotted this subset of genes across the whole GB transcriptome, once ordered according to the significance and direction of change as a result of the pairwise comparison between H3ac high and H3ac low , and we observed that the genes upregulated with TSA tended to be predominantly enriched in the H3ac high GB, whereas the genes downregulated after the treatment tended to be enriched in H3ac low GB ( Figure 2I) (χ 2 = 754.94, d.f. = 19, p-value < 0.00001 for downregulation; χ 2 = 533.38, d.f. = 19, p-value < 0.00001 for upregulation). 9 and 14 and low basal gene expression and CGI hypermethylation ( Figure 2G), suggest-275 ing that well-expressed genes did not need to increase their histone acetylation state to 276 increase their rate of expression [33]. Transcriptional signatures associated with differential levels of histone H3 acetylation in glioblastomas: (A) H3K9,14ac levels in the selected samples for the differential expression analysis; GBs expressing lower and similar levels of acetylated protein compared to grade II gliomas (H3ac low and H3ac high , respectively). The data are expressed as the mean ± SEM. **, p-value < 0.005; Mann-Whitney U test between both types of GB. (B) Tridimensional PCA of the whole transcriptomes for selected gliomas of our cohort, in which RNA was optimal for deep sequencing. (C) Median of the Pearson correlation coefficients (r) between the expression of each gene in the 249-gene signature (divided upon direction of change) and the bulk level of each histone modification. *, p-value < 0.05. (D) Significantly enriched GO terms (FDR < 0.05, DAVID) in the DEG between H3ac low and H3ac high ; the plot shows −log-adjusted p-value (bars) and number of genes (next to each bar) for each GO term. (E) Box-and-whisker plot of the average expression across GB of the H3ac high -signature. Quartiles of expression of the whole GB transcriptome are shown. (F) Average CpG methylation expressed as β-values (0, fully unmethylated; 1, fully methylated) of the H3ac high -signature compared to the whole GB transcriptome, and divided in CGI, shores (<2 kb from CGIs), shelves (>2 kb and <4 kb from CGIs), and open seas (>4 kb from CGIs). The p-value < 0.05 was calculated using Student's t-test between the signature and the rest of genes. (G) The GB transcriptome was ranked according to the significance and direction of the gene expression change in the pairwise comparison between H3ac high and H3ac low and divided into bins of 500 genes; median values were calculated for each bin. Adjusted p-values were calculated using Student's t-test with Bonferroni correction between each bin of genes and whole population of genes; the plot indicates all significant p-values grouped in the genes showing relative higher expression in H3ac high compared to H3ac low GBs. (H) Left panel, Western blots for total H3 and H3K9,14ac from untreated and treated GB-derived cultures with TSA. The raw blot is shown in Figure S3. Right panel, quantification of the blots after normalization by total H3 signal. The data are expressed as the mean ± SEM. n = 6 for each condition. **, p-value < 0.005, Mann-Whitney U test related to the control (vehicle) condition. (I) Distribution of DEGs obtained in the pairwise comparisons between TSA-treated and untreated cultures across the whole GB transcriptomes, ordered as in (G) and divided into bins of 1000 genes. The number of down-and upregulated genes from our lists was counted in each bin and fold enrichment was calculated (number of observed genes/number of expected genes if equally distributed).

The Transcriptional Signatures Associated with Histone Acetylation Can Be Linked to Differential Overall Survival Only in the Samples Showing the Most Extreme Expression Levels
When considering adjusted p-value cutoffs (<0.05), we found a significant overlap of 33 genes enriched in H3ac high that were also upregulated with TSA treatment (χ 2 = 29.01, d.f. = 1, p-value < 0.00001, Tables S1 and S2). Among them, we found the following genes that have been reported in previous studies: alpha-2-macroglobulin-like 1 (A2ML1), in which protein levels are decreased after in vitro treatments of glioma cells lines [34][35][36]; CD83, a dendritic cell maturation marker in which expression is reduced within glioblastoma microenvironment [37,38]; dedicator of cytokinesis 8 (DOCK8), which is associated with progression from low-grade to higher-grade gliomas, and can be lost in chromosome imbalances [39]; gamma-aminobutyric acid type A receptor subunit alpha 4 (GABRA4), a potential inhibitor of cancer progression that is downregulated across glioma grades and negatively correlated with inflammation [40]; hippocalcin-like 4 (HPCAL4), an unfavorable risk factor defined as a core gene in gliomagenesis that becomes upregulated in GB compared to LGG [41,42]; myocyte enhancer factor 2C (MEF2C), which acts as an oncogene in glioma cells by activating the JAGGED-NOTCH pathway [43]; solute carrier family 17 member 7 (SLC17A7), a tumor suppressor gene with bivalent properties (with co-occupancy of H3K4me3 and H3K27me3 marks) that is overexpressed in patients with glioma-related seizures and downregulated in oligodendrogliomas [44][45][46]; SLIT and NTRK-like family member 4 (SLITRK4), which exhibits a complex expression patterns across brain tumors [47]; synaptotagmin 1 (SYT1), which is defined as a hub gene in GB, associated with poor prognosis and copy number variation [41,47,48]; and TNF superfamily member 18 (TNFSF18), which is moderately beneficial as part of immunotherapy approaches in murine gliomas [49]. In conclusion, high levels of H3K9,14ac might be simultaneously associated with genes with poor and good outcomes.
However, this subset of 33 might have a net prognostic value as a whole. First, the median expression of this subset enabled the classification of our cohort into two main branches as highly and weakly acetylated tumors ( Figure 3A), confirming that this signature was tightly dependent on histone H3 acetylation. However, both groups were not different in terms of overall survival (log-rank p-value = 0.5); only when comparing those samples exhibiting the most differential expression values, the median expression of the 33-gene signature had a trend towards a good outcome ( Figure 3B). The same approach was unsuccessful using the genes enriched in H3ac high (218-gene signature) and H3ac low (31-gene signature), as we obtained log-rank p-values of 0.22 and 0.35, respectively. Taking the TCGA samples for which the "days to death" parameter was available, we used the 33-gene signature to divide the GB into two groups: TCGA-subset low and TCGAsubset high ( Figure 3C). The former subset (TCGA-subset low ) accumulated more somatic mutations than the TCGA-subset high GBs, in which genes involved in the regulation of histone acetylation, followed by histone methylation and chromatin remodeling, were more frequently mutated (Figure 3C), suggesting a relationship between somatic mutation of the epigenetic enzymatic machinery and the altered gene expression patterns of the H3ac-linked signature. Independent of this classification, the median expression of the 33-gene signature in the TCGA database was significantly associated with a good prognosis after comparing samples containing the most extreme values ( Figure 3D).

Discussion
In this work, we describe the transcriptional and genetic alterations that affect epigenetic regulators in GB. Such alterations may influence the tumoral epigenetic landscape with potential consequences for survival, quality of life, and response to treatment. By comparing overall modifiers according to their type of epigenetic activity (DNA and histone modifications, readers of these modifications, and chromatin remodeling factors), we observed that the most consistent alterations occurred in genes controlling histone acetylation, independent of whether they participate in writing or erasing this post-translational modification. More attention has been paid for HDAC, reporting either decreased or increased gene expression levels across histological grades and normal tissue [50][51][52][53]. In agreement with references [50] and [53], we found the downregulation of HDAC4 and HDAC11 in GB to be related to LGG. The class II HDAC4 was found to inhibit the expression of the proapoptotic p21/CDKN1A gene [13], the high expression levels of which may have a poor prognostic value in the GB mesenchymal subtype [54], and may promote radioresistance [55]; therefore, the general downregulation in GB of TCGA and REMBRANDT cohorts compared to LGG may probably be part of a homeostatic response to refrain tumor malignancy. Apart from a recent study suggesting a role in glioma progression and therapeutical resistance [56], little is known about the role of KAT6B in glioma, a lysine acetyltransferase of the MYST family with oncogenic or tumor suppressor properties depending on the cancer [57]. However, the KAT6B activity seems to be more critical in regulating the acetylation of lysine 23 [58], and cannot explain, in principle, the decreased levels of bulk H3K9,14ac in our samples. Overall, complex interactions between the altered activities of KAT, HDAC, and cofactors may result in a net balance of H3K9,14ac reduction in GB.
Although gliomas showed increased histone H3 acetylation levels compared to normal brain tissue in a Western blotting analysis, no conclusive result was obtained when comparing gliomas of different grades in a previous report [50]. In our work, we were able to observe decreased levels of bulk histone H3 acetylation in GB resections compared to low-grade gliomas. Of note, such reduction paralleled the levels of H3.3 in our previous work using the same extract proteins [23]. Interestingly, the mutation of H3-3A (also known as H3F3A) in pediatric high-grade gliomas interferes with the activity of EZH2 methyltransferase/PRC2 complex, leading to an overall decrease in the levels of H3K27me3 concomitant to an increase in H3K27ac [59]; the observation for such hyperacetylation was extended to other histone lysines [60]. According to these results, it seems that the depletion of wild type H3.3 in adult gliomas (causing a loss of function [61]) associates with histone deacetylation, whereas the H3.3 mutation in pediatric gliomas (causing a gain-of-function) is linked to histone hyperacetylation.
However, differential bulk levels in H3K9,14ac were not sufficient for a profound differential gene expression among GB, suggesting that additional histone covalent modifications (or alternative histone acetylation marks, such as H3K27ac, a mark that is part of transcriptional activation domains in the glioma chromatin [31,33]) are required for understanding the GB transcriptome arrangement. However, we observed that genes with low (but not absent) expression were more likely to exhibit variation in the levels of H3K9,14ac, suggesting that an increase in H3ac does not necessarily influence the transcription of highly expressed genes, as they might already contain sufficient acetylated histones at their regulatory regions. This fact should be considered regarding the use of HDACis (or alternatively, KAT enhancers) to manipulate the expression of tumor suppressor genes. Nonetheless, most well-known HDACis (such as TSA) are relatively unspecific and in-hibit entire classes of HDAC that not only regulate histones, but also transcription factors that may also affect gene expression; additionally, histone-independent actions should not be disregarded [3]. Such wide actions may explain contradictory outputs regarding the induction of genes of good and poor prognosis simultaneously (see examples in the Results section), despite the downregulation of mitosis-and proliferation-related genes. This phenomenon may result in ineffective actions over the tumoral mass in failed clinical trials, in addition to the poor pharmacokinetics of most HDACis [15,16]. Moreover, HDACi treatment evokes a complex homeostatic response (e.g., resulting in the downregulation of regulatory genes with opposite actions such as KAT6A and HDAC10 in our cellular preparations; Table S2) that may be key to understanding the induced chemoresistance against these epigenetic drugs after long-term regimes [17,62]. Therefore, functional assays are required to assess the tumorigenic effects of candidates. Furthermore, chromatin immunoprecipitation will determine the direct targets that will be helpful for the design of more selective compounds for pharmacological modulation of histone acetylation within the signature. Overall, knowing the precise targets of these epigenetic inhibitors is crucial in predicting the clinical outcome after treatment and in improving therapeutics in patients.

Supplementary Materials:
The following supporting information can be downloaded at https: //www.mdpi.com/article/10.3390/cells12030374/s1. Figure S1: Differential expression changes in epigenetically related genes in TCGA and REMBRANDT gliomas. Figure S2: Raw Western blots of histone H3 modifications in glioma resections. Figure S3: Histone H3 acetylation in TSA-treated glioma culture. Table S1: Complete RNA-seq results from differential expression analysis between H3aclow/H3achigh glioblastomas. Table S2: Complete RNA-seq results from differential expression analysis between TSA/untreated cultures. Institutional Review Board Statement: The study was conducted according to the guidelines of the Declaration of Helsinki and approved in 28 October 2016 and 31 May 2019 by the local ethics committee (Comité de Ética de la Investigación de Cádiz) according to the national and regional laws and regulations concerning biomedical research on human samples, personal data protection, and the use of biobank services.
Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
In-house datasets generated and analyzed during the current study are available in the GEO repository under accession numbers GSE147391 and GSE185861.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.