Transcriptomic and epigenetic profiling of ‘diffuse midline gliomas, H3 K27M-mutant’ discriminate two subgroups based on the type of histone H3 mutated and not supratentorial or infratentorial location

Central pathology review

High-grade glioma cases were reviewed centrally to confirm the diagnosis according to the 2007 WHO classification and its 2016 update as previously described [10, 20].

Specific immunostainings were performed to detect nuclear expression of the trimethylation mark at position K27 of the histone 3 tail (1:1000, polyclonal rabbit antibody, Diagenode, Belgium) as well as nuclear expression of the K27M form of histone H3 (1:1000, polyclonal rabbit antibody, Millipore, CA).

Derivation and culture of glioma stem-like cells (GSCs)

GSCs were derived from DIPG tumors at diagnosis as previously described [19]. Briefly, tumor cells were mechanically dissociated from biopsies within 24 h of surgery, and further cultured as an adherent monolayer in laminin-coated flask (Sigma) in neural stem cells medium consisting of NeuroCult NS-A proliferation medium (Stemcell technologies) supplemented with heparin (2 μg/mL, Stemcell technologies), human-basic FGF (20 ng/ml, Peprotech), human-EGF (20 ng/ml, Peprotech), PDGF-AA (10 ng/ml, Peprotech), and PDGF-BB (10 ng/ml, Perprotech). Medium was renewed every other day, and passaging performed when cells reached 80% confluence using Accutase (Thermo).

Case selection for overall survival analysis and gene expression profiling by microarray

Case and sample selections for methylation analysis

Gene expression profiling was also conducted by either microarray or RNA sequencing for 5 of these tumors. Eight glioma stem-like cell (GSC) cultures derived from patient biopsies at diagnosis and matching primary tumors were analyzed similarly [19].

Methylation profiling

DNA was extracted from tumors and genome-wide DNA methylation analysis was performed using either the Illumina HumanMethylation450 BeadChip (450 k) or EPIC arrays. DNA methylation analysis was performed with custom approaches as previously described [12, 23]. DNA methylation profiles from 50 K27M pHGG were compared to defined supratentorial tumor subgroups, i.e. G34R-H3.3 mutated (n = 10), MYCN (n = 10) and PDGFRA/pedRTK1 (n = 10) subgroup tumors. For t-SNE analysis (t-Distributed Stochastic Neighbor Embedding, Rtsne package version 0.11), 428,230 uniquely mapping autosomal probes in common between the 450 k and EPIC arrays were used. The input for the t-SNE calculation is 1-Pearson correlation, weighted by variance. Clustering analyses were performed using the beta values of the top 10,000 most variably methylated probes by standard deviation. Methylation probes in the heatmap representation were reordered by unsupervised hierarchical clustering using Pearson correlation distance and median linkage.

Microarray gene expression profiling

Gene expression analysis was conducted on an Agilent platform as previously described [2] but using RUV4 correction of batch effects [7] implemented in the R package ruv. GE data from DIPG were collected from one of our previous study [2] and microarray analysis was performed for 75 additional pHGG tumors located outside the brainstem. PCA, k-means and t-SNE analysis were performed using the same parameters as for RNA-seq data on the probes associated with the highest standard deviation. One hundred and twenty genes accounting for 0.79% of the entire probeset were selected.

RNA-seq gene expression profiling

RNA-seq was performed on 21 primary tumor samples. Libraries were prepared using the TruSeq stranded mRNA sample preparation kit according to the supplier recommendations and paired-end sequencing was conducted on Illumina NextSeq500 to generate a mean of 150 million reads of 75 base pairs by sample. Trimmed reads were then mapped using tophat2 (v2.1.0) and bowtie2 (v2.2.5) first to the reference transcriptome, then to the reference genome for the remaining reads. Genes with a row sum of raw counts over the studied samples equal to or below 10 were filtered out to remove non-expressed genes. We handled outliers as default using minReplicatesForReplace = 7 in DESeq() function used to estimate size factors, dispersion and model coefficients. Distances between samples were computed by using ‘1-Pearson correlation coefficient’ as the distance measure. PCA and t-SNE analysis were performed on the 250 genes associated with the highest variance to keep the same proportion of genes selected with the microarray analysis. All samples were projected on the two first principal components computed with rlog transformation of the counts of the 120 genes with the highest standard deviation. Using Rtsne package (v 0.11), we applied t-SNE on the same data matrix with the Pearson correlation as a distance and the following parameters: theta = 0, perplexity = min(floor((ncol(rlog_VariableGenes)-1)/3), 30), check_duplicates = FALSE, pca = FALSE, max_iter = 10,000, verbose = TRUE, is_distance = TRUE.

RNA-seq was also performed on 6 distinct GSC models using TruSeq stranded total RNA sample preparation kit according to the supplier recommendations (Illumina) and then processed similarly as primary tumors.

Histone ChIP-sequencing and data processing

ChIP-seq of H3K27me3 epigenetic modification was performed in 6 GSC models at Active Motif according to proprietary methods. The 75-nt sequence reads were generated on a Illumina NextSeq 500 platform, mapped using BWA algorithm and peak calling was performed using SICER1.1 algorithm [27] with cutoff FDR 1e-10 and gap parameter of 600 bp. False positive ChIP-seq peaks were removed as defined within the ENCODE blacklist [5]. Overlapping intervals between the different samples were merged, and the average number of normalized reads in the different samples were calculated for these 16,977 genomic intervals defined as ‘bound regions’. These bound regions were separated for further analysis in overlapping or not overlapping gene loci using Genecode annotation (gencode.v19.chr_patch_hapl_scaff_annotation.gtf). PCA for all samples were generated after scaling to unit variance using the PCA function from the FactoMineR package (v1.41) and plotted using Factoextra (v1.0.5).

Merging of the 3 biological replicates of H3.1- or H3.3-K27M subgroups was performed using bigWigMerge tool (UCSC kent utils, http://hgdownload.soe.ucsc.edu/admin/exe/macOSX.x86_64/). Heatmaps of H3K27me3 ChIPseq enrichment across genomic loci were calculated using deepTools version 1.5.11. ComputeMatrix was used with regions of either +/− 5 kb or +/− 10 kb around the center of the genomic intervals for ‘bound regions’ or TSS for differentially expressed genes, respectively. Heatmaps were plotted with or without k-means (k = 5) and their average profiles of ChIP-Seq enrichment in the same − 10/+ 10 kb genomic intervals were also generated for each k-means group. The bigWig files (all signal) were annotated with chipSeeker package using UCSC hg19 known gene annotation and visualizated by peakAnno and Vennpie.

Survival curve comparisons

The distribution of overall survival (OS) was calculated according to the Kaplan-Meier method and all survival function estimate comparisons were performed in PRISM software using a log-rank test. OS was calculated from the date of histo-radiological diagnosis until death of patient from disease or last contact for patients who were still alive.

Histone H3 K27M midline pHGG and K27M DIPG display similar gene expression profiles and survival but differ significantly from other high-grade gliomas

We conducted microarray gene expression profiling of the single center cohort from Necker Enfants Malades hospital of 119 pHGG with histone H3 genotype previously determined either by Whole Genome Sequencing [25] or targeted Sanger sequencing [2] (Table 1 and Additional file 1: Table S1). In total, 131 distinct gene expression microarrays were hybridized as 12 samples were analyzed twice to control for potential batch effects. The 12 duplicated samples were located in close proximity on PCA plot when projected on the 2 first principal components, confirming the appropriate removal of a batch effect in our dataset (data not shown). Then, we selected a set of genes associated with the highest standard deviation for subsequent tumor classification analysis (n = 120). Gene Ontology over-representation analysis showed an enrichment of genes involved in brain development (Bonferroni adjusted p-value 8.29e-09) and morphogenesis (adjusted p-value 2.67e-07); many homeobox genes belong to this later set of genes reflecting probably the differences in tumor location.

Principal component analysis (PCA) was performed to highlight the principal sources of variation among pHGG tumors (Additional file 2: Table S2: weight in principal components 1 and 2 associated to each gene). The DIPG group appeared rather homogenous as all samples clustered together in the PCA plot, separated from the tumors originating from thalamic or cortical areas of the brain which were more scattered (Fig. 1a). Considering the mutational status of histone H3 genes, the results clearly showed that H3-K27M tumors, whichever their pontine or thalamic location, could be separated from the wild-type and G34R/V tumors on the first principal component (Fig. 1b). Indeed, all H3-K27M mutated thalamic and spinal tumors were close to DIPG samples, whereas histone H3 wild-type thalamic tumors are distributed on the right side of the plot among histone H3 wild-type non-thalamic midline and cortical tumors. Interestingly, 5 DIPG without any mutation in H3F3A, HIST1H3B/C and HIST2H3A/C were located within the K27M DIPG subgroup. These tumors all showed H3K27-trimethylation loss by immunohistochemistry (Additional file 3: Figure S1A). Unsupervised K-means analysis was performed on the same dataset (using k = 2 which showed the best BIC value) and led to a similar conclusion as the k-mean group 1 corresponded to H3-K27M and H3-wild type samples presenting H3K27-trimethylation loss whereas k-mean group 2 contained all other pHGG tumors (Additional file 3: Figure S1B). Additionally, this GE dataset was analyzed in parallel with another dimension reduction and visualization technique for high-dimensional data, t-SNE, as it was shown to be more robust than PCA with respect to outliers. Moreover, t-SNE has also been frequently used for pediatric brain tumor classification based on DNA methylation profiles [23]. This analysis of GE profiles led to similar observations, re-iterating the similarity of H3-K27M thalamic midline and DIPG (Additional file 3: Figure S1C-D). Histone H3-G34R/V samples were tightly clustered together in particular in this analysis, thus reflecting a strong similarity in their gene expression profiles.

In agreement with these observations, a huge number of differentially expressed genes were identified between H3-K27M and H3-wild-type tumors as well as between H3-K27M and H3-G34R tumors (adjusted p-value < 0.01) in comparison with the other contrasts (Additional file 3: Figure S1E and Additional file 4: Table S3). In contrast, only 14 genes were significantly modulated between H3-G34R and H3-WT tumors which were not discriminated by PC1 (Fig. 1b). Gene ontology analysis showed an important enrichment of modulated genes associated with neurogenesis (4.37e-12 – 1.41e-11) and neuron differentiation (3e-07 - 4.87e-13) signaling pathways in both contrasts (H3-K27M vs. H3-WT and H3-K27M vs. H3-G34R) , as well as an upregulation of genes involved in ion transmembrane transport (1.6e-04) and apoptotic processes (5.9e-18) and an downregulation of genes linked to cell cycle (8.3e-11) and gliogenesis (6.87e-10) in the case of the comparison between H3-wild-type and H3-K27M tumors. Also, geneset enrichment analysis (GSEA) identified the ontology GO_oligodendrocyte_differentiation (Enrichment Score 0.70) enriched in upregulated genes in H3-K27M tumors and GO_Cerebral_cortex_neuron_differentiation (Enrichment Score − 0.79) enriched in genes upregulated in H3-G34R. Some of these biological processes were previously identified as significantly enriched in differentially expressed genes when comparing H3.1- and H3.3-K27M tumors [2].

We next conducted a survival analysis on this cohort (n = 119) using only location information and then both location and H3 mutation status for patient stratification (Fig. 1c and d). DIPG and thalamic tumors were associated with similar poor prognosis, i.e. 11.1 and 10.8 months median OS, respectively. Non-thalamic midline tumors exhibited the best prognosis (median OS not reached), whereas tumors arising in the cortex presented an intermediate outcome with a median survival around 30.5 months (p-value < 0.0001, Fig. 1c). Focusing on Kaplan-Meier estimates for midline tumors, our data clearly indicate that H3-WT non-thalamic midline have a significantly higher overall survival, whereas the other midline malignant gliomas (mostly thalamic), with or without alteration of histone H3 genes, display equivalent poor survival (p-value < 0.0001, Fig. 1d).

Methylation profiling separates HIST1H3B and H3F3A K27M tumors

Previous studies have shown that genome-wide DNA methylation data can provide a robust classification of pediatric brain tumors into clinically meaningful epigenetic subgroups mostly characterized by recurrent genetic alterations [14, 15, 22]. Consequently, we compared the methylation profiles of K27M-mutated diffuse midline gliomas (including DIPG) to G34R-mutated tumors and well-characterized supratentorial tumors without mutation in histone H3 genes, i.e. MYCN and PDGFRA tumor subgroups (Fig. 2a, Table 2) [15]. Eighty primary tumor samples were used in this analysis and t-SNE visualization of the DNA methylation data was conducted. We confirmed that H3-G34R, PDGFRA and MYCN subgroups constitute 3 distinct homogenous entities, as they defined three distinct clusters. All H3-K27M samples were located on the opposite side of the 2D representation, reflecting important differences in the methylome compared to these three well-defined pHGG subgroups. This observation was thus concordant with our results on GE profiling by microarray.

In addition, the same methylation profiling splits H3-K27M samples in two subgroups that corresponded to either H3.1 or H3.3 mutated tumors. The obvious separation of these tumors in an analysis containing other very distinct biological entities clearly indicated the significant difference between them. Also, the unique H3.2-K27M sample appeared closer to H3.1-K27M than H3.3-K27M samples (Fig. 2a).

The same classification by t-SNE was repeated for the subset of H3-K27M mutated midline gliomas. First, t-SNE analysis did not reveal a segregation of these samples according to their location, as all DIPG and thalamic midline were scattered in the 2D plot (Fig. 2b). Conversely, when considering samples based on the mutated histone H3 gene, the t-SNE analysis clearly highlighted two non-overlapping subgroups corresponding to H3.1/H3.2-K27M and H3.3-K27M classes (Fig. 2c). This observation indicates that H3.3-K27M DIPG are closer to other midline H3.3-K27M HGG than to H3.1-K27M DIPG. The histone H3 variant affected by the K27M substitution thus has a stronger correlation with the modulation of DNA methylation profile than the tumor location across the midline.

An additional analysis was performed on a subset of 21 primary DIPG tumors and 8 glioma stem-like cells (GSCs) deriving from these same biopsies. The sample classification by unsupervised hierarchical clustering confirmed the previous result, with two main clusters corresponding to H3.3-K27M and H3.1/2-K27M samples (Fig. 3a). Additionally, we observed that the majority of the GSCs clustered with their corresponding primary tumors indicating the close similarity of their methylome profile. This consequently underlined that GSC population remained very similar to their primary counterpart with respect to DNA methylation, reflecting the variation in DNA methylation observed between the two subgroups of H3.1 and H3.3 mutated tumors.

RNA-seq profiling also discriminates HIST1H3B/C and H3F3A K27M mutated gliomas

DNA methylation is a relatively stable component of the epigenome involved in the establishment and maintenance of distinct gene expression patterns. Consequently, we decided to evaluate if the different DNA methylation profiles were associated with distinct transcriptome profiles. PCA analysis of GE measurements by microarrays did not clearly discriminate H3.1 and H3.3 mutated tumor samples (Fig. 1b). Indeed, even if H3.1-K27M DIPG were closer to each other in the 2-dimensional PCA plot, they are surrounded by H3.3-K27M samples. However, as microarray data were generated in several batches, we could not exclude that this could obscure the dataset, despite the use of a batch correction method. Therefore, we took advantage of a RNA-seq study of 21 new H3-K27M DIPGs samples which appeared more suitable as it provides an exhaustive measurement of transcriptome in contrast to microarray analysis. Grouping of tumors based on their RNA-seq expression profiles in either t-SNE or PCA classifications confirmed the discrimination of H3.1-K27M from H3.3-K27M tumors observed in our DNA methylation study (Fig. 3b and Additional file 5: Figure S2).

H3.1- and H3.3-K27M mutations are associated with different genomic distribution of the H3K27me3 epigenetic mark

To complement the description of the epigenetic landscape of DIPGs, we decided to assess the direct epigenetic consequences of the K27M mutation in H3.1-K27M and H3.3-K27M tumors, i.e. the loss and redistribution of the trimethylation mark at position K27. We thus took advantage of 6 GSC models of different genotypes as an expandable source of material and interrogated the genome-wide distribution of this mark with ChIP-seq. We conducted a PCA analysis on all genomic regions showing enrichment for this mark and could again separate the samples based on the type of mutated histone (Fig. 3c). This implied that despite the similar biochemical impact of both H3.1-K27M and H3.3-K27M mutated histones on PRC2, interferences of the two different mutated histones H3 on PRC2 is not similarly distributed in the genome and this may induce meaningful changes in the epigenetic landscapes of these DIPGs.

The average signal of all H3K27me3 peak regions was equivalent in H3.1- and H3.3-K27M tumors (Additional file 6: Figure S3A, top). We plotted the ChIP-seq signal of the union of H3K27me3 peaks identified from all the GSC samples and did not observed a specificity of the genomic regions presenting a H3K27me3 deposition between H3.1 and H3.3 K27M subgroups but rather a difference in the level of the signal enrichment (Additional file 6: Figure S3A, bottom). However, the genome distribution of H3K27me3 showed an enrichment of peaks in distal intergenic regions in H3.3-K27M GSC and in contrast in intronic and at a lesser extent exonic regions in H3.1-K27M cells (Fig. 4a).

Subsequently, we divided the active genomic loci presenting H3K27me3 deposition according to their location. Overall, the global signal was significantly higher in regions overlapping genes versus intergenic regions (Additional file 6: Figure S3B-C). The heatmap displayed that the majority of regions overlapping genes presented a lower level of H3K27me3 in H3.1-K27M samples (Additional file 6: Figure S3C). This observation was concordant with a higher number of expressed genes in H3.1, which are moreover associated with higher expression levels (Additional file 6: Figure S3D).

The clustering of genomic intervals by k-means highlighted the existence of distinct subgroups of loci, smaller gene clusters with a significantly higher level of H3K27me3 in H3.1-mutated cells (C1, C2 and C3 for genic and C1’ and C3’ for intergenic loci, Fig. 4b and c respectively). Consequently, the vast majority of loci were enriched in K27me3 in H3.3-mutated cells.

Thereafter, we focused on the differentially expressed genes identified by RNA-seq between the 2 subgroups (adj p-value< 0.01). The results showed that upregulated genes in H3.1-K27M are associated with a lower average signal of H3K27me3 and vice et versa for downregulated genes (Fig. 5a). The ChIP-seq read coverage are shown for two representative genes, OLIG2 and SLFN11, even if distinct localization of H3K27me3 is observed between both, located mainly upstream and downstream of transcribed region for SLFN11 and overlapping gene body for OLIG2 (Fig. 5b and c). However, gene expression modulations between H3.1- and H3.3-mutated tumors was not restricted to changes in H3K27me3 deposition at the close vicinity of the gene body as shown for HOXD8 which is upregulated and associated with an increased H3K27me3 deposition in H3.1-K27M tumors (Fig. 5d).