The DNA methylation landscape of glioblastoma disease progression shows extensive heterogeneity in time and space

Abstract

Glioblastoma is characterized by widespread genetic and transcriptional heterogeneity, yet little is known about the role of the epigenome in glioblastoma disease progression. Here, we present genome-scale maps of DNA methylation in matched primary and recurring glioblastoma tumors, using data from a highly annotated clinical cohort that was selected through a national patient registry. We demonstrate the feasibility of DNA methylation mapping in a large set of routinely collected FFPE samples, and we validate bisulfite sequencing as a multipurpose assay that allowed us to infer a range of different genetic, epigenetic, and transcriptional characteristics of the profiled tumor samples. On the basis of these data, we identified subtle differences between primary and recurring tumors, links between DNA methylation and the tumor microenvironment, and an association of epigenetic tumor heterogeneity with patient survival. In summary, this study establishes an open resource for dissecting DNA methylation heterogeneity in a genetically diverse and heterogeneous cancer, and it demonstrates the feasibility of integrating epigenomics, radiology, and digital pathology for a national cohort, thereby leveraging existing samples and data collected as part of routine clinical practice.

Figure 1. DNA methylation landscape of glioblastoma disease progression

(a) Integrative analysis of longitudinal DNA methylation data (RRBS) with matched magnetic resonance (MR) imaging data (morphology, segmentation), clinical annotation data (treatment, progression, IDH mutation status), and histopathological data (segmentation, morphology, immunohistochemistry) using statistical methods and machine learning. TMZ: Temozolomide; RTX: Radiation therapy; PC: Palliative care. (b) Overview of the progression cohort, summarizing the disease courses of 112 primary glioblastoma patients with IDH-wildtype status and at least two tumor resections (ordered by time of first surgery). (c) DNA methylation profiles for primary (t=1) and recurring (t=2) tumors at three glioblastoma-linked gene loci (SFRP2, PLXNB2, MGMT). Gene annotations and epigenome segmentations for astrocytes were obtained from the WashU EpiGenome Browser (session URL: http://epigenome-segmentation.computational-epigenetics.org). (d) DNA methylation levels at individual CpGs indicative of the glioma CpG island methylator phenotype (G-CIMP), shown for IDH-mutated reference samples (which are G-CIMP positive) and for the IDH-wildtype primary glioblastoma samples from the study cohort (which are G-CIMP-negative). The DNA methylation fold change between IDH-mutated and IDH-wildtype samples is based on published data for G-CIMP marker CpGs.

Figure 2. Glioblastoma transcriptional subtypes inferred from DNA methylation profiles

(a) Overview of the machine learning approach for predicting transcriptional subtypes from DNA methylation. (b) Transcriptional subtype heterogeneity in the progression cohort, as indicated by class probabilities of the RRBS-based subtype classifier. Samples are grouped and ordered by their dominant subtype. (c) Distribution of class probabilities across different regions of the same tumor (indicated by Roman numbers) and across different surgeries (indicated by Arabic numbers) for two patients with multisector samples. N denotes the total number of tumor samples profiled for each patient. (d) Riverplot depicting transitions in the predicted transcriptional subtype between primary and recurring tumors. The number of tumor samples (N) in each state is indicated. Only patients for which the sample-specific subtype classifier achieved high accuracy (ROC AUC > 0.8) were included in the analyses for panel d to f. (e) Kaplan-Meier plots showing progression-free and overall survival for patients of the progression cohort stratified by predicted transcriptional subtypes (left) and subtype switching (right). The number of patients (N) is provided in Supplementary Table 3. (f) Heatmap displaying the DNA methylation levels of the most differential CpGs between the three transcriptional subtypes, excluding CpGs with more than 100 missing values. Only tumor samples that were classified with high class probabilities (>0.8) for their dominant subtype were included in this analysis. (g) LOLA region set enrichment analysis for differentially methylated CpGs between transcriptional subtypes (binned into 1-kilobase tiling regions). Region sets from astrocytes or embryonic stem cells with an adjusted p-value below 0.001 are shown. (h) Schematic depicting the calculation of ‘DNA methylation inferred regulatory activity’ (MIRA) scores. High MIRA scores reflect local demethylation at the binding sites of a specific transcription factor across the genome, which indicates high regulatory activity of that factor. (i) DNA methylation profiles (top row) and corresponding MIRA scores (bottom row) for three sets of transcription factor binding sites that were hypomethylated in the mesenchymal subtype (CTCF, EZH2, KDM4A) and for the binding sites of three key regulators of pluripotency (POU5F1/OCT4, NANOG, SOX2). The number of tumor samples (N) in each group is provided in Supplementary Table 4.

Figure 3. DNA methylation and the tumor microenvironment

(a) Level of tumor-infiltrating immune cells across transcriptional subtypes (Cla: classical, Mes: mesenchymal, Pro: proneural), as measured by quantitative immunohistochemistry for the indicated marker proteins. (b) Illustrative immunohistochemical stainings for FOXP3 and CD45ro in three individual patients (one for each subtype). (c) Kaplan-Meier plots showing progression-free and overall survival for patients in the progression cohort stratified by the level of CD163 positive and CD68 positive immune cell infiltration. The number of patients (N) is provided in Supplementary Table 3. (d) Level of tumor-infiltrating immune cell for samples originating from the patients’ first surgery (primary tumor), second surgery (recurring tumor), or third surgery (where available) in the progression cohort, and for the first (and only) surgery in the validation cohort. (e) Illustrative immunohistochemical stainings for three marker proteins (CD68, CD8, CD163) in three individual patients undergoing with changing levels of tumor-infiltrating immune cells between primary and recurring tumors. (f) Level of tumor-infiltrating immune cells (CD3, CD8, CD68) and proliferating cells (MIB-1 positive) compared between patients that were assigned to different progression types based on magnetic resonance (MR) imaging: Classic T1 (claT1), cT1 relapse / flare-up (cT1), and T2 diffuse (T2). (g) ROC curves for DNA methylation based prediction of immune cell infiltration levels. N denotes the number of tumor samples in each group. The number of tumor samples (N) for the association tests (panels a, d, f) is provided in Supplementary Table 4.

Figure 4. DNA methylation and histopathological tumor characteristics

(a) Fraction of proliferating (MIB-1 positive) cells across transcriptional subtypes (Cla: classical, Mes: mesenchymal, Pro: proneural). (b) Kaplan-Meier plots showing progression-free and overall survival stratified by the fraction of proliferating (MIB-1 positive) cells. (c) ROC curves for DNA methylation based prediction of the fraction of proliferating (MIB-1 positive) cells. N denotes the number of tumor samples in each group. (d) Clustered heatmap for the column-scaled DNA methylation levels of the most predictive genomic regions (5-kilobase tiling regions) from the classifier predicting the fraction of proliferating (MIB-1 positive) cells. (e) Histogram showing the distribution of DNA methylation levels across the most predictive genomic regions from the classifier predicting the fraction of proliferating (MIB positive) cells. (f) Average nuclear eccentricity (AVG) and its coefficient of variation (COV) for tumors that shift to a sarcoma-like phenotype during disease progression and those that retain a stable histological phenotype. (g) Illustrative H&E stains of matched primary and recurring tumors from one patient that shifted to a sarcoma-like phenotype. (h) Comparison of additional tumor properties (as indicated by the y-axis labels) between tumors that shifted to a sarcoma-like phenotype during disease progression and those that retained a stable histological phenotype. (i) ROC curves for DNA methylation based prediction of average nuclear eccentricity (AVG) and its coefficient of variation (COV). N denotes the number of tumor samples in each group. (j) Kaplan-Meier plots showing progression-free and overall survival stratified by whether or not their tumors shifted to a sarcoma-like phenotype. The number of tumor samples (N) for the association tests (panels a, f, h) is provided in Supplementary Table 4. The number of patients (N) in the Kaplan-Meier analysis (panels b, j) is provided in Supplementary Table 3.

Figure 5. DNA methylation heterogeneity in glioblastoma disease progression

(a) Illustration of ‘proportion of discordant reads’ (PDR) and ‘epi-allele entropy’ (EPY) as complementary measures of epigenetic tumor heterogeneity. Genomic regions can have high values for one but not the other measure (Locus 1 and 2, top row) or for both measures simultaneously (Locus 3, top row). Regions that undergo extensive changes in their epi-allele composition (i.e., in the distribution of subclonal DNA methylation patterns) are referred to as eloci (Locus 2 and 3, bottom row). (b) Sample-wise PDR and EPY scores averaged across promoter regions for first surgery (t=1, primary tumor) and second surgery (t=2, recurring tumor) of the progression cohort and first (and only) surgery of the validation cohort. The samples with the 20% highest and lowest scores are shown in green and orange, respectively. N denotes the number of tumor samples. (c) Relative epi-allele frequencies across promoter regions for all samples highlighted in green and/or orange in panel b. Samples are ordered by epi-allele diversity, and samples with high epi-allele diversity (left) tend to show high PDR and/or EPY (green dots). For clearer visualization, the “0000” majority epi-allele with a frequency of ~70% to 90% is not displayed. 0: unmethylated, 1: methylated (d) Kaplan-Meier plots showing progression-free and overall survival stratified by PDR and EPY values, respectively. The number of patients (N) is provided in Supplementary Table 3. (e) Scatterplot depicting the (negative) association of the number of differentially methylated promoters between primary and recurring tumor samples (normalized to the number of assessed promoters) and the time between first and second surgery. r: Pearson correlation. The p-value was calculated with a one-sided Pearson’s test. N denotes the number of tumor samples.

Figure 6. DNA methylation differences between primary and recurring tumors

(a) Scatterplot depicting the association of DNA methylation at promoter regions between primary and recurring tumors. Promoters that were differentially methylated in at least five patients are highlighted. r: Pearson correlation. N denotes the number of data points (samples × promoters). (b) Barplots (top) depicting the number of patients that show significant gain or loss of DNA methylation at the differentially methylated promoters highlighted in panel a; dot and line plots (bottom) showing the change in DNA methylation associated with disease progression (percentage points, pp) for patients following (red) or not following (blue) the cohort trend. Trend lines were calculated using the loess() function in R. (c) Histogram showing the number of “trend” and “anti-trend” patients based on the Manhattan distance between the maximal trend at differentially methylated promoters (DNA methylation values of 0% or 100%, respectively) and the observed difference in DNA methylation for each patient. (d) Kaplan-Meier plots showing progression-free and overall survival for trend and anti-trend patients according to panel c. (e) Gene set enrichment analysis of those genes that recurrently gain (top) or lose (bottom) promoter DNA methylation during disease progression. N denotes the number of genes. (f) Kaplan-Meier plots showing progression-free and overall survival stratified by progression-linked changes in DNA methylation at the promoters of Wnt signaling genes (top and bottom 30% of patients included). (g) ROC curve and clustered heatmap for DNA methylation based prediction of overall survival, comparing patients with survival above the median in the progression cohort (“long”) to patients with survival below the median in the validation cohort (“short”). N denotes the number of tumor samples in each group. (h) LOLA region set enrichment analysis for the top-2,000 tiling regions (5 kilobases) most predictive of overall survival. Region sets from astrocytes, brain tissue, or embryonic stem cells with an adjusted p-value below 0.001 are shown. The number of patients (N) in the Kaplan-Meier analysis (panels d, f) is provided in Supplementary Table 3.