Mapping chromatin accessibility and active regulatory elements reveals pathological mechanisms in human gliomas

Abstract

Chromatin structure and accessibility, and combinatorial binding of transcription factors to regulatory elements in genomic DNA control transcription. Genetic variations in genes encoding histones, epigenetics-related enzymes or modifiers affect chromatin structure/dynamics and result in alterations in gene expression contributing to cancer development or progression. Gliomas are brain tumors frequently associated with epigenetics-related gene deregulation. We perform whole-genome mapping of chromatin accessibility, histone modifications, DNA methylation patterns and transcriptome analysis simultaneously in multiple tumor samples to unravel epigenetic dysfunctions driving gliomagenesis. Based on the results of the integrative analysis of the acquired profiles, we create an atlas of active enhancers and promoters in benign and malignant gliomas. We explore these elements and intersect with Hi-C data to uncover molecular mechanisms instructing gene expression in gliomas.

Fig. 1. Experimental workflow.

Freshly resected glioma specimens were acquired from patients representing major glioma malignancy groups: pilocytic astrocytomas (WHO grade I, PA, n = 11), diffuse astrocytomas (WHO grade II and III, DA, n = 7) and glioblastomas (WHO grade IV, GBM, n = 15). Clinical data about each patient was collected. Glioma tissues were homogenized to a single cell suspension and experimentally defined numbers of cells were split into portions that were processed for RNA-seq, ATAC-seq, ChIP-seq, or DNA methylation experiments. Differential analysis of epigenomic and transciptomic profiles across glioma grades was performed, followed by identification of related regulatory mechanisms. Results from all genome-wide experiments were made available in a genome browser (http://regulomics.mimuw.edu.pl/GliomaAtlas/).

Fig. 2. Construction of an atlas of regulatory sites in gliomas.

A Numbers of samples from gliomas of different grades for which ATAC-seq, ChIP-seq, and RNA-seq data were collected. B Representative peaks of ATAC-seq, ChIP-seq, RNA-seq at the NRXN2 gene in pilocytic astocytomas (PA04) and glioblastoma (GB01) samples. Note high expression of the NRXN2 in PA in comparison to GB, which correlates with the lack of the ATAC-seq signal in GB (shadowed blocks). C Total numbers of peaks identified in all samples across performed experiments. D Numbers of identified common and variable active regulatory elements, genes and transcripts. E Example of promoter-enhancer interactions identified using Hi-C data (Won et al., 2016). Brown dots represent “anchors”, while brown horizontal lines depict predicted interaction sites.

Fig. 3. Global characterization of chromatin structure and its relationship with gene expression.

A Commonality of histone marks, open chromatin regions and gene expression patterns between samples. Color scale represents the percentage of patients sharing a particular genomic element status. Bar height represents a number of genomic elements (scaled to 100%). B Venn diagram shows intersection of the promoter regions marked by H3K4me3, H3K27ac and open chromatin sites detected with ATAC-seq. C Evolutionary conservation of the identified promoters (green) and enhancers (purple) measured by PhastCons 100-Way scores. Mean scores for the identified regulatory regions are compared to mean scores obtained for random genomic intervals. P values were calculated with one-sided Mann–Whitney U test. Data are represented as boxplots in which the box shows the quartiles of the dataset, the middle line is the median and the whiskers extends to the largest or smallest value no further than 1.5 × the inter-quartile range. D Correlation between H3K4me3 (green), H3K27ac (purple) and transcript expression. Normalized ChIP-seq coverage is plotted against an expression level of transcripts of protein coding genes, divided into quantiles. Data are represented as violin plots with a nested boxplot, where shape indicates the distribution of data. The box shows the quartiles of the dataset, the middle white dot is the median, whiskers extend to 1.5 × IQR past the low and high quartiles. E Hierarchical clustering of samples based on gene expression profiles. Shades of blue indicate tumor grade (from light = normal brain to dark = Grade IV and black PG). In the table on the right side the subtype of glioblastoma samples (GBM) is reported as well as assignment of glioma samples to one of Pan Glioma RNA Expression Clusters (LGr) as well as IDH1 gene mutation status (IDH). F Correlation between H3K27ac ChIP-seq coverage on enhancers and transcript expression. Numbers above bars indicate numbers of enhancers in each group. The p value has been estimated empirically as explained in “Methods”.

Fig. 4. Chromatin activity profiles indicate the presence of the normal brain signature in pilocytic astrocytomas.

A, B Scatter-plots representing abundance of the H3K4me3 (A) and H3K27ac (B) marks obtained for PAs (X-axes) and DA/GBs (Y-axes). A single point on the scatter-plots represents an average abundance (ChIP-input) of a corresponding histone mark around a transcription start site [TSS + /-2 kb]. P values shown above the plots were estimated with the Wilcoxon rank-sum test. C Comparison of abundance values shown on 4A (X-axis) and 4B (Y-axis). The dashed frame indicates the strongest difference in H3K27ac abundance (top 10%) values. D Heatmaps showing H3K27ac abundance around the TSS for the top 10% regions (each row corresponds to a TSS, sorted by signal abundance). E ATAC-seq profiles of 4 PA samples and 4 DA/GB samples (see “Methods” for details) around TSS overlap with the regions selected on 4C. The blue and red lines correspond to PA and DA/GB samples, respectively. F Average H3K27ac profiles around TSS. The blue and red lines correspond to the top 10% regions selected in 4C and were computed for the PA and DA/GB samples, respectively. The gray lines correspond to regions not selected on 4C. The thick lines correspond to the average profiles, while colored areas give reference of the confidence interval for mean (CI). G Differences of the gene expression levels computed either for genes corresponding to the top 10% selected regions (left boxplot) or to the rest of the regions (right boxplot). P value shown above the plot and indicating differential regulation was estimated with the two-sided Wilcoxon rank-sum test. Data are represented as boxplots in which the box shows the quartiles of the dataset, the middle line is the median and the whiskers extends to the largest or smallest value no further than 1.5 × the inter-quartile range. H Results of Gene Ontology over-representation analysis performed for the genes corresponding to the top 10% regions. The barplot shows scaled Bonferroni corrected p values from the one-sided Fisher’s exact test and the vertical, dashed red line stands for the significance threshold (P = 0.05). I Logos of top DNA-binding motifs enriched in the top 10% regions. The enrichment was computed with regard to the rest of the analyzed regions. P values indicated above the logos stand for the significance of the enrichment and were computed with a tool from MEME-Suite (see “Methods” for details).

Fig. 5. A group of promoters targeted by PRC2 exhibits H3K4me3 hypermethylation and DNA hypomethylation in PAs.

A Hierarchical clustering of samples based on the presence of H3K4me3 at the promoters. Color scale and numbers on the heatmap indicate a similarity between H3K4me3 patterns at the promoters (dark blue, 1 - identical; white, 0 - maximally dissimilar). Red box – the identified cluster of PA samples. B Expression of transcripts associated with promoters exhibiting H3K4me3 hypermethylation in PAs. The p value was calculated with one-sided Mann–Whitney U test. C Enrichment of binding sites from the ENCODE transcription factor ChIP-seq data in PAs specific promoters compared to promoters active in all samples. The number of peaks for each DNA-binding protein is given in parentheses; enrichment values per megabase are shown next to each bar and indicated with proportional coloring (red – highest, blue-lowest). D DNA methylation in promoters exhibiting H3K4me3 hypermethylation in PAs and having an EZH2 binding site. Methylation levels were determined with WGBS. Samples with IDH mutations were excluded from the plot. Colors indicate sample grades (light orange – PA, dark orange – DA, red – GB/PG). E DNA methylation in promoters exhibiting H3K4me3 hypermethylation in PAs and having an EZH2 binding site in an independent cohort, determined by HumanMethylation450 BeadChip array. Colors indicate sample grades (light orange – PA, red – GB). F GO terms enriched in the set of the PA-specific genes with EZH2 binding sites in promoters. The enrichment was calculated using the PANTHER Overrepresentation Test . More detailed results are provided in Supplementary Fig. 5E. In panels (B), (D), and (E) data are represented as boxplots in which the box shows the quartiles of the dataset, the middle line is the median and the whiskers extends to the largest or smallest value no further than 1.5× the inter-quartile range.

Fig. 6. Intersection of the Atlas data with fetal brain Hi-C data reveals FOXM1-ANXA2R axis that impacts survival of GB patients and regulates glioma invasion.

A Graphical representation of chromatin contacts from Hi-C data. B Heatmap showing normalized average expression of 17 differentially expressed genes in (FDR < 0.01) selected in the contact analysis. C Barplot showing significance of enhancer-promoter co-activation for five genes computed as the Pearson correlation of H3K27ac enrichments in putative enhancers and H3K4me3 enrichments in corresponding promoters. D FOXM1 expression in gliomas. The adjusted p values computed with edgeR Bioconductor correspond to expression differences. Logo represents FOXM1 DNA-binding motif. Boxplots show the quartiles of the dataset, lines indicate median and the whiskers the largest or smallest value. E, F Scatterplots showing co-expression of FOXM1 and ANXA2R as log-transformed FPKMs in our (E) and TCGA (F) datasets. Pearson correlation coefficients shown above the plots. G Kaplan–Meier survival curves plotted for patients stratified by ANXA2R or FOXM1 expression; differences assessed with the log-rank test. H ANXA2R and FOXM1 expression in normal human astrocytes (NHA), patient-derived (WG4 and IPIN) and established glioma cells (LN229, LN18, U87) determined by RT-qPCR; data normalized to GAPDH mRNA; p values calculated on the raw data with two-sided t test. I Capture-C profile representing interactions around the viewpoint in the ANXA2R promoter (green arrow). Region enrichments indicate the putative enhancer (purple arrow). Fragment-based raw data are visualized as grey dots, whereas blue dots represent smoothed data. J Significant enrichment of FOXM1 binding at the ANXA2R promoter and enhancer, and CCNB1 and PKL1 promoters (positive controls). No significant FOXM1 enrichment detected at H19, myoglobin and CCND1 promoters (negative controls) or with a neutral IgG antibody. Results calculated as % of input, mean ± SD (n = 5), significance with ratio paired t test, two-sided. K, L Knockdown of FOXM1 in WG4 glioma cells verified by qPCR and Western blotting; densitometry of immunoblots determined from three experiments, mean ± SD, (K) two-sided paired t test, (L) two-sided t test. M Reduced ANXA2R expression in siFOXM1-transfected cells determined with qPCR; mean ± SD, n = 4, two-sided paired t test. N Knockdown of FOXM1 reduces migration and invasion of glioma. Cell migration after scratch quantified in four fields before and 24 h after scratch. Results shown as mean ± SD, n = 4. Invasion determined with Matrigel and Transwell^® inserts. Images from five fields acquired using fluorescence microscopy, cell nuclei counted with ImageJ software. Results shown as mean ± SD, n = 4 (in duplicate). p values calculated using paired t test, two-sided.