Therapeutic targeting of ependymoma as informed by oncogenic enhancer profiling

Abstract

Genomic sequencing has driven precision-based oncology therapy; however, the genetic drivers of many malignancies remain unknown or non-targetable, so alternative approaches to the identification of therapeutic leads are necessary. Ependymomas are chemotherapy-resistant brain tumours, which, despite genomic sequencing, lack effective molecular targets. Intracranial ependymomas are segregated on the basis of anatomical location (supratentorial region or posterior fossa) and further divided into distinct molecular subgroups that reflect differences in the age of onset, gender predominance and response to therapy. The most common and aggressive subgroup, posterior fossa ependymoma group A (PF-EPN-A), occurs in young children and appears to lack recurrent somatic mutations. Conversely, posterior fossa ependymoma group B (PF-EPN-B) tumours display frequent large-scale copy number gains and losses but have favourable clinical outcomes. More than 70% of supratentorial ependymomas are defined by highly recurrent gene fusions in the NF-κB subunit gene RELA (ST-EPN-RELA), and a smaller number involve fusion of the gene encoding the transcriptional activator YAP1 (ST-EPN-YAP1). Subependymomas, a distinct histologic variant, can also be found within the supratetorial and posterior fossa compartments, and account for the majority of tumours in the molecular subgroups ST-EPN-SE and PF-EPN-SE. Here we describe mapping of active chromatin landscapes in 42 primary ependymomas in two non-overlapping primary ependymoma cohorts, with the goal of identifying essential super-enhancer-associated genes on which tumour cells depend. Enhancer regions revealed putative oncogenes, molecular targets and pathways; inhibition of these targets with small molecule inhibitors or short hairpin RNA diminished the proliferation of patient-derived neurospheres and increased survival in mouse models of ependymomas. Through profiling of transcriptional enhancers, our study provides a framework for target and drug discovery in other cancers that lack known genetic drivers and are therefore difficult to treat.

Extended Data Figure 1. DNA fingerprint analysis of ependymoma sequence data

(a-b) Unsupervised clustering of ChIP-seq, RNA-seq, WES, WGS, RNA-seq, and Illumina DNA methylation profiles with genotypes that have an average heterozygosity score greater than 0.25 in the Heidelberg (n = 25 independent samples) (a) and Toronto Cohort (n = 18 independent samples) (b).

Extended Data Figure 2. Summary of genome sequencing and copy number data

(a) Number of somatic single nucleotide variants (SNVs) detected per ependymoma sample. (b) Frequency of somatic mutations detected across the Heidelberg ependymoma cohort (n = 24 independent samples). (c) Unsupervised hierarchical clustering of copy number alterations detected by WGS in primary ependymoma samples (n = 24 independent samples).

Extended Data Figure 3. Preprocessing and clustering of ependymoma H3K27ac profiles

(a-b) Box plots of H3K27ac enhancer profiles (n = 556,676 enhancer loci evaluated per sample) pre-quantile normalization for both Heidelberg (n = 24 independent samples) (a) and Toronto (n = 18 independent samples) (b) cohorts as compared to Roadmap Epigenomics and ENCODE cohort (n = 98 independent samples). Box plots are shown with the center (median), upper and lower quartile range, and dotted line indicating minima and maxima per sample. (c-d) Box plots of H3K27ac enhancers post-quantile normalization for both Heidelberg (n = 24 independent samples) (c) and Toronto (n = 18 independent samples) (d) cohorts as compared to Roadmap Epigenomics cohort (n = 98 independent samples). (e-f) Unsupervised hierarchical clustering of enhancer profiles as measured using the Top 10,000 variant enhancer loci identified in the Roadmap Epigenomics cohort with the Heidelberg (n = 122 independent samples) (e) and Toronto cohort (n =116 samples) (f) and compared in a pair-wise fashion using a Spearman Correlation.

Extended Data Figure 4. Ependymoma enhancer supporting data

(a) Number of unique H3K27ac peaks detected by MACS1.4 (p < 1e-9 cutoff) with increasing sample number in the Heidelberg cohort (n = 24 independent samples) (b) Boxplot of gene expression values comparing typical enhancer (n = 9826 genes) versus super enhancer (SE) (n = 1682 genes) associated genes. Statistical analysis was assessed using a two-sided Wilcoxon-Rank Sum Test. Box plots of gene expression values are shown with the center (median), upper and lower quartile range, and dotted line indicating minima and maxima (c) Frequency of enhancer and SE regions as a function of size in basepairs (d) Dot plots illustrating the numbers of SE detected in the Heidelberg (n = 24 independent samples), Toronto (n = 18 independent samples), and Normal brain (n = 7 independent samples) cohort. The center bar in the dot plot indicates the mean. (e) Heatmap illustrating significant gained and lost enhancer loci in both ependymoma cohorts as compared to normal brain samples. Comparisons were evaluated using a two-sided Wilcoxon-rank sum test with false discovery rate (FDR) correction and a cut-off of less than FDR < 0.05 (f) Example plots of normalized and scaled H3K27ac RPKM profiles at example ependymoma candidate genes in Heidelberg-ependymomas and normal brain (n = 32 independent samples). (g) Gene expression comparing ependymoma SE associated genes derived from Johnson et al., (2010) (n = 83 independent samples) with normal brain (n = 172 independent samples). Statistical analysis was assessed using a two-sided Wilcoxon-Rank Sum Test. (h) Table comparing the number and percent confirmation between the Heidelberg (n = 24 independent samples) and Toronto ependymoma cohorts (n = 18 independent samples) (i) G-profiler pathway enrichment analysis of ependymoma specific SE associated genes in the Toronto cohort (n = 18 independent samples), with statistical significance determined using a hypergeometric test (j) Overlap analysis as measured by a two-sided binomial test between tumour specific ependymoma SEs and cancer census genes from the Catalogue of Somatic Mutations in Cancer database (k) Classification of tumour specific ependymoma SE genes also found in the COSMIC database as tumour suppressor genes (n = 12), oncogenes (n = 26), or unknown (n = 21).

Extended Data Figure 5. Subgroup specific enhancers of ependymoma

(a-b) Heatmap of all subgroup-specific active enhancers detected in ependymoma in independent samples in the Heidelberg (n = 24 independent samples) and Toronto (n = 18 independent samples) cohort (c) Box plot of gene expression for ependymoma SE-SSEA associated genes in the Heidelberg cohort (n = 24 independent samples). Comparisons were made using a two-sided Wilcoxon-Rank Sum test. Box plots of gene expression values are shown with the center (median), upper and lower quartile range, and dotted line indicating minima and maxima. (d-f) Venn diagrams of the number and percent of subgroup-specific SE loci validated between the Heidelberg and Toronto cohort (g-h) Non-negative factorization of ependymoma SE profiles in the Heidelberg (n = 24 independent samples) and Toronto (n = 18 independent samples) cohort (i) Normalized H3K27ac profiles for subgroup specific genomic example loci in the Heidelberg cohort with at minimum 3 biological replicates per subgroup with the exception of ST-EPN-SE shown as a biological duplicate. (j) This has been corrected to read: G-profiler pathway enrichment analysis of ependymoma subgroup-specific SE associated genes in the Heidelberg cohort (n = 24 independent samples) with statistical significance determined using a hypergeometric test. (k-n) H3K27ac profiles surrounding the EPHB2 (k) and CCND1 (m) locus in the Heidelberg cohort with at minimum 3 biological replicates per subgroup with the exception of ST-EPN-SE shown as a biological duplicate. EPHB2 (l) and CCND1 (n) expression by RNA-seq across ependymoma subgroups in the Heidelberg cohort with horizontal bar indicating the median value and each dot representing an independent ependymoma sample (n = 24 independent samples).

Extended Data Figure 6. Workflow describing the functional validation of ependymoma SE genes

(a) Workflow of super enhancer target gene prioritization for functional evaluation. (b) Bar chart comparing the top ranked SE associated genes against top ranked genes detect by RNA-seq defined as significantly gained or over-expressed compared to normal brain controls across all ependymoma samples (n = 42 independent samples). Significant genes were identified by a two-sided Wilcoxon rank sum test with FDR correction and ranked by FDR corrected p-value cut-off of less than a minimum of 0.05.

Extended Data Figure 7. RNA interference of ependymoma SE genes

(a) Individual shRNA time-course knockdown experiments in EP1-NS (ST-EPN-RELA) cells, using two shRNA constructs (shRNA.1 and shRNA.2) as compared to two controls (shCONTROL.1 and shCONTROL.2). Shown are time-course experiments for 19 genes performed in 6 technical replicates. (b) Ependymoma cell viability (EP1-NS) following treatment with shRNAs targeting super enhancer associated genes over a 7-day time course (alphabetically ordered). Cell viability data for treatment with non-targeting controls: shCONTROL.1 (black), shCONTROL.2 (grey), and for two gene-specific shRNA constructs: shRNA.1 (red) and shRNA.2 (pink).

Extended Data Figure 8. Validation of ependymoma subgroup specific SE genes

(a) H3K27ac profiles at the ependymoma specific SE locus IGF2BP1 in the Heidelberg cohort (n = 24 independent samples) with at minimum 3 biological replicates per subgroup with the exception of ST-EPN-SE shown as a biological duplicate. (b) IGF2BP1 gene expression derived from RNA-seq data of the Heidelberg cohort (n = 24 independent samples) and horizontal bar for each subgroup that indicates the mean. (c-d) Normalized survival of PF-EPN-A (S15) primary cultures (c) and EP1-NS (d) cell cultures following shRNA knockdown of IGF2BP1 with 2 independent non-overlapping shRNA constructs as compared to shCONTROL.1. Experiments performed as 6 technical replicates and independently validated in 3 biological replicates. Horizontal black bar indicates mean value for each treatment. (e) H3K27ac profiles at the ependymoma specific SE locus CACNA1H in the Heidelberg cohort with at minimum 3 biological replicates per subgroup with the exception of ST-EPN-SE shown as a biological duplicate. (f) H3K27ac profiles of a ST-EPN-RELA model (EP1-NS), a PF-EPN-A model (S15), and a normal neural stem cell control surrounding the CACNA1H locus performed in biological duplicates. (g) CACNA1H gene expression derived from RNA-seq data of the Heidelberg cohort (n = 24 independent samples) and horizontal bar for each subgroup that indicates the mean. (h-i) Normalized survival of PF-EPN-A (S15) primary cultures (k) and EP1-NS (l) cell cultures following shRNA knockdown of CACNA1H with 2 shRNA constructs as compared to shCONTROL.1. Experiments performed as 4 technical replicates and independently validated in 3 biological replicates. Horizontal black bar indicates mean value for each treatment. (j) Normalized cell survival of EP1-NS, S15, and NSC194 cells treated with increasing concentrations of Mibefradil. Shown are technical triplicates and replicated in biological triplicates (k) Overlay of ATAC-seq and H3K27ac seq data centered upon ATAC-seq peak regions identified in the ST-EPN-RELA cell culture EP1-NS. (l) CRISPR dCAS9 targeting of CACNA1H active enhancers impairs CACNA1H expression. H3K27ac-seq (top panel) and ATAC-seq (bottom panel) surrounding the CACNA1H locus, indicating regions targeted by CRISPR-dCAS9 sgRNA complexes. Region 1 (R1) indicates a negative control region devoid of H3K27ac (green), while Region 2-4 (R2-4) indicate experimental regions under evaluation. Experiments replicated in biological duplicates. (m) Gene expression for various sgRNA constructs relative to a 'dummy' targeting control (D103), negative control (green), and uninfected control. All group comparisons were made using a two-sided Wilcoxon-Rank sum test, error bars indicated as standard deviation in all dot plots of this figure, and horizontal bar indicating the mean value. Experiments were replicated in biological triplicates.

Extended Data Figure 9. Validation of ependymoma transcription factors

(a-b) Gene expression of 'high activity' transcription factors (ranked < 50) (a) and 'low activity' transcription factors (ranked > 50) (b) in ependymoma (n = 83 independent samples) versus normal brain tissue (n = 172 independent samples). Box plot is described by the median value (horizontal bar), interquartile range, and dotted line representing the data range. Comparison between groups was assessed using a two-sided Wilcoxon rank sum test. Boxplots indicate the range of data points. (c) Constituent enhancer activity in the central nervous system of developing zebrafish embryos derived from subgroup specific SEs identified in ependymoma

Extended Data Figure 10. Putative cell lineage programs of origin uncovered by TF mapping

(a-c) Immunohistochemical staining of foxj1 at day 13.5 of mouse embryonic development (E13.5). Staining in discrete regions encompassing the choroid plexus and ependymal layer are shown in the forebrain (b) and hindbrain (c). (d) Log2 normalized gene expression of FOXJ1 in ependymoma (n = 83 independent samples) as compared to independent sample cohorts of the following tissue types: Normal Brain (n = 172), Pediatric Glioma (n = 53), Glioblastoma (n = 84), Atypical Rhabdoid Teratoid Tumours (n = 18), Medulloblastoma (n = 62), and Pilocytic Astrocytoma (n = 41). Horizontal bar indicates the mean value. (e) Subgroup specific gene expression of FOXJ1 derived from Pajtler et al., Cancer Cell 2015 (n = 209 independent samples) Error bars indicate standard deviation, interquartile range, and horizontal bar indicating the median. (f) Gene set enrichment analysis demonstrating significant enrichment of the FOXJ1 transcriptional program derived from E14.5 mouse embryos specifically in PF-EPN-B tumours (n = 209 independent samples). FDR corrected significance evaluated by gene set enrichment analysis. (g) Significant FOXJ1 gene expression correlations with proteins known to regulate cilia assembly and function. P-values for significant positive/negative correlations have been corrected for multiple testing using the Bonferroni method. (h-m) FOXJ1 gene set enrichment plots of PF-EPN-A (h), PF-EPN-B (i), PF-EPN-SE (j), ST-EPN-RELA (k), ST-EPN-YAP1 (l), and ST-EPN-SE (m) ependymomas. FDR corrected significance evaluated by gene set enrichment analysis, n = 209 independent samples.

Figure 1. H3K27ac profiles define active regulatory elements of ependymoma

(a) Unsupervised hierarchical clustering of the top 10,000 variant enhancer loci detected in ependymoma as compared to the Roadmap Epigenomics Consortium samples, n = 143 independent samples. (b-c) Inflection plot indicating ependymoma super enhancers (SEs) identified. (d-e) Venn diagram depicting the number of shared enhancers (d) and SEs (e) between the Heidelberg (n = 24) and Toronto (n = 18) ependymoma independent sample cohorts. (f) Quantitative RT-PCR knockdown efficiency of 15 ependymoma SE associated genes (n = 3 technical replicates, error bars indicated as standard deviation. Results were reproduced in independent biological duplicates). (g) Percent of top ependymoma SE genes that demonstrate greater than 50% decrease in viability over 7 days. Cell survival from knockdown of each gene assayed and independently replicated as biological triplicates.

Figure 2. Active enhancers delineate subgroups of ependymoma

(a-b) Unsupervised hierarchical clustering of all H3K27ac enhancer loci in Heidelberg (n = 24) and Toronto (n = 18) independent sample cohorts (c) Combined t-SNE analysis of the top 10,000 variably methylated Illumina 450K CpG probes. (d) Combined t-SNE analysis of all enhancer loci. n = 43 independent samples. (e-f) t-SNE analysis of the H3K27ac marked ependymoma SE regions in ependymoma. n = 42 independent samples. (g-l) Inflection plot indicating the SEs with subgroup specific enhancer activity (SE-SSEA) in ependymoma. n = 24 independent samples. (m) G-profiler pathway analysis of ependymoma subgroup SE associated genes with significant enrichment indicated as the FDR corrected p-value. n = 24 independent samples.

Figure 3. Transcription factor circuitries of ependymoma

(a) DNA motifs enriched within shared ependymoma typical enhancers that overlay with ATAC-seq peaks derived from the EP1-NS cell culture model as determined by HOMER motif analysis (see Methods and Supplementary Table 18) (b) Heatmap of transcription factors (TFs) ranked by predicted activity using core circuitry analysis (left panel) and presence or absence of self-loop activity (right panel). n = 18 independent samples of Toronto cohort. (c-f) shRNA constructs targeting SE associated genes ordered by normalized cell survival. Highlighted in red are shRNAs targeting SE associated core TFs. Each gene assayed with 6 technical replicates and replicated in 3 independent biological experiments. (g-l) Connections between subgroup specific transcription factors (TF) integrated with gene expression in subgroups of ependymoma. n = 24 independent samples.

Figure 4. Active regulatory maps identify candidate drugs against ependymoma

(a) Pie graph of candidate drug compounds detected by integrating shared SEs with the Washington University Drug Gene Interaction Database. (b-d) Ependymoma cells and NSC (NSC1) controls treated with JQ1 (b), AZD1775 (c), and AZD4547 (d) for 72h and assessed using an Alamar Blue stain. Error bars as standard deviation. Experiment performed as six technical replicates and replicated in biological triplicates. (e) Kaplan - Meier curve for immunodeficient mice bearing H.612 ependymomas, treated with either vehicle or AZD4547 (25 mg/kg/d). Significance of endpoint difference was assessed using a log-rank test. Median survival ratio of Treatment (AZD4547) : Control (vehicle) is 44 days : 33 days, and reported as a ratio of 1.333 with a 95% confidence interval of 0.4677 to 3.801.